Improved ELMv1-ECA Simulations of Zero-Curtain Periods and Cold-season CH4 and CO2 Emissions at Alaskan Arctic Tundra Sites

Field measurements have shown that cold-season methane (CH4) and carbon dioxide (CO2) emissions contribute a substantial portion to the annual net carbon emissions in permafrost regions. However, most earth system land models do not accurately reproduce cold-season CH4 and CO2 emissions, especially over the shoulder (i.e., thawing and freezing) seasons. Here we use the Energy Exascale Earth System Model (E3SM) land model version 1 (ELMv1-ECA) to tackle this challenge 10 and fill the knowledge gap of how cold-season CH4 and CO2 emissions contribute to the annual totals at Alaska Arctic tundra sites. Specifically, we improved the ELMv1-ECA soil water phase-change scheme, environmental controls on microbial activity, and cold-season methane transport module. Results demonstrate that both soil temperature and the duration of zero-curtain periods (i.e., the fall period when soil temperatures linger around 0°C) simulated by the updated ELMv1-ECA were greatly improved, e.g., the Mean Absolute Error in zero-curtain durations at 12 cm depth was reduced by 15 62% on average. Furthermore, the simulated cold-season emissions at three tundra sites were improved by 84% and 81% on average for CH4 and CO2, respectively. Overall, CH4 and CO2 emitted during the early cold season (Sep. and Oct.), which often includes most of the zero-curtain period in Arctic tundra, accounted for more than 50% of the total emissions throughout the entire cold season (Sep. to May). From 1950 to 2017, both CO2 emissions during the 12 cm depth zerocurtain period and during the entire cold season showed increasing trends, for example, of 0.26 gC m-2 year-1 and 0.38 gC m20 2 year-1 at Atqasuk. This study highlights the importance of zero-curtain periods in facilitating CH4 and CO2 emissions from tundra ecosystems. https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Cold-season carbon emissions from the Arctic tundra could potentially offset warm-season net carbon uptake under 21 st 25 century warming climate (Commane et al., 2017;Oechel et al., 2014;Oechel et al., 2000;Koven et al., 2011;Piao et al., 2008;Natali et al., 2019;Belshe et al., 2013;Fahnestock et al., 1998;Jones et al., 1999). Field measurements have indicated large cold-season CO2 losses over Arctic tundra ecosystems (Oechel et al., 2014;Natali et al., 2019). Also, CH4 emitted from September to May were found to contribute more than 50% of the annual total CH4 emissions from Alaska upland tundra sites (Zona et al., 2016;Taylor et al., 2018). Despite the importance of cold-season carbon emissions and their sensitivity to 30 changing climate, prevailing earth system land models do not accurately reproduce cold-season CH4 and CO2 emissions and their contributions to the annual budgets, largely because of the poorly understood mechanisms of cold-season soil heterotrophic respiration and therefore uncertain numerical representations (Natali et al., 2019;Zona et al., 2016;Commane et al., 2017). Thus, it remains challenging to assess the response of permafrost carbon dynamics to Arctic warming and to predict future annual carbon budgets with current Earth System Models (ESMs). 35 In ESM land models, soil environment influences soil microbial heterotrophic respiration (HR) and decomposition of soil organic carbon (SOC) mainly through applying prescribed temperature and moisture functions to modify base decomposition rates. These functions, however, rely heavily on empirical or semi-empirical relationships which are highly uncertain (Sierra et al., 2017;Sierra et al., 2015;Yan et al., 2018;Moyano et al., 2013;Tang and Riley, 2019;Rafique et al., 2016;Bhanja and 40 Wang, 2020;Kim et al., 2019). Specifically, the temperature sensitivities of soil carbon decomposition is often represented with a value (i.e., the increase in respiration rate from a 10°C increase in temperature) that is fixed at 1.5 or 2.0 (Meyer et al., 2018). However, the values of Q10 are controversial (Davidson and Janssens, 2006). Some studies found a uniform Q10 across biomes and climate zones, e.g., as 1.4 (Mahecha et al., 2010). Other studies demonstrated that Q10 varies with environmental conditions, ecosystem types, and soil texture (Meyer et al., 2018;Graf et al., 2011;Kim et al., 2019), showing 45 a large spatial heterogeneity with generally higher values in the high-latitudinal regions (Zhou et al., 2009). In addition, Wilkman et al. (2018) reported a temporal heterogeneity in Q10 over the Alaskan Arctic Tundra and suggested a higher value (e.g., 2.45) for early summer (e.g., June) but lower value (e.g., 1.58 to 1.67) for the peak growing season (e.g., July).
Dynamic decomposition temperature sensitivities are also consistent with theory of microbial dynamics (Tang and Riley, 2015). Also, the response of HR to changes in soil moisture is commonly expressed by empirical relationships in ESMs, 50 which vary substantially (Sierra et al., 2015;Yan et al., 2018;Moyano et al., 2013). Although in-situ measurements reveal that microbial respiration occurs under very cold conditions (e.g., even when soil temperature is lower than -15 °C) (Natali et al., 2019;Zona et al., 2016), most process-based models completely shut down microbial activity due to limited liquid water in freezing and subfreezing soils, and few modelling studies have closely investigated the HR-moisture relationships in frozen conditions. 55 https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License.

Study Sites and Data
We assembled daily observations of CO2 and CH4 fluxes from 2013 to 2017 at five eddy-covariance flux tower sites in 90 Alaska's North Slope tundra (Figure 1) from the Arctic-Boreal Vulnerability Experiment (ABoVE) project (2015 -2017) (Oechel and Kalhori, 2018) and Carbon in Arctic Reservoirs Vulnerability Experiment (CARVE) flight campaign (2013 -2014) (Zona et al., 2016). CARVE CO2 measurements were not available; therefore, monthly winter-time CO2 flux data at the ABoVE towers assembled by Natali et al. (2019) are included to complement CO2 observations from 2013 to 2014. The five sites include three eddy covariance (EC) towers at Barrow (i.e., the Barrow Environmental Observatory (BEO) tower, 95 the Biocomplexity Experiment South (BES) tower, and the Climate Monitoring and Diagnostics Laboratory (CMDL) tower), one tower at Atqasuk (ATQ) and another at Ivotuk (IVO) which is located at the foothills of the Brooks Range. BES and CMDL are collocated with each other with sensors installed at different heights (i.e., 2 m for BES and 5 m for CMDL).
Vegetation at Barrow is mainly moist acidic tundra. Instrument height at ATQ and IVO is 2 m and 4 m, respectively. ATQ is a well-drained upland site, and the vegetation consists of moist-wet coastal sedge tundra and moist-tussock tundra surfaces. 100 Vegetation at IVO is polar tundra. Table S1 provides basic information including geolocations, vegetation mosaic, and climatologic air temperature at the sites. (Tables numbered with a prefix "S" are include in the supplementary file, which will not be repeated in the following context throughout the manuscript.) ABoVE and CARVE provide soil temperature and moisture measurements at various depths from 5 cm to 40 cm. The 105 Permafrost Laboratory, Geophysical Institute of University of Alaska Fairbanks (GIPL-UAF), provides daily subsurface soil temperature observations down to various depths at permafrost sites across Alaska(http://permafrost.gi.alaska.edu/sites_map) (Romanovsky et al., 2009). We used the GIPL-UAF permafrost sites that are collocated with the ABoVE sites to complement the ABoVE observations at deeper depths, including BR2 (down to 15 m) and IV4 (down to 1 m). We first filled missing gaps vertically by fitting a polynomial to the soil temperature profile (Kurylyk and Hayashi, 2016) on a daily 110 scale, then screened out outliers by examining the daily time series. Further, we aggregated both the ABoVE and the GIPL-UAF soil temperature measurements to ELMv1-ECA soil layer node depths using the inverse distance weighting method (Tao et al., 2017), and then averaged the two sets of aggregated observations. We used the assembled subsurface temperature observation datasets to evaluate the ELMv1-ECA simulated soil temperature profiles and the zero-curtain periods.

115
Due to the discontinuity of observed soil moisture over time and along with the vertical depth, evaluating ELMv1-ECA simulated liquid water content at layer node-depth was limited. We matched soil-moisture observations to the vertically closest model layer, and then evaluated the simulated volumetric fraction of soil liquid water content at layers for time periods during which observations were available. In addition, we used ABoVE soil moisture measurements to derive sitescale soil porosity and organic carbon content (see Section 3.2). 120 https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License.

Modifications to E3SM Land Model (ELM)
The E3SM land model version 1 (ELMv1-ECA) couples essential biogeophysical and biogeochemical processes that solve terrestrial ecosystem energy, water, carbon, and nutrient dynamics (Golaz et al., 2019;Zhu et al., 2019). Figure 2 illustrates the coupling and interactions among the three components. In the appendix, we describe in detail its subsurface soil 125 thermodynamics, the carbon decomposition module, and the methane module that are of particular relevance to our study.
Here we identify the potential problems of ELMv1-ECA that are responsible for the underestimation of cold-season CH4 and CO2 emissions and summarize the modifications made to ELMv1-ECA, emphasizing the model enhancements, shown by the ellipses with red boundaries in Figure 2.

Phase Change Scheme 130
We first improved ELMv1-ECA's numerical representation of coupled water and heat transport with freeze-thaw processes via improving the phase-change scheme. The freeze-thaw processes of soil water within ELMv1-ECA is simulated in a decoupled way, i.e., it solves soil temperatures ignoring the latent heat associated with phase change, determines the mass change of soil water required to adjust the initially solved soil temperature to the freezing point (i.e., 0°C; ), adjusts the soil liquid and ice content by mass and energy conservation, and then readjusts temperatures after accounting for the heat 135 deduction or compensation resulted from melting or freezing (see the detailed description in the Appendix A). The underlying assumption here is, taking the freezing process as an example, the available liquid water at the initially solved temperature ( ) will be completely frozen, releasing latent heat ( ) to bring up back to . Then, the estimated phase-change rate will be tuned down and the current temperature (i.e., ) will be readjusted if the to-be-increased ice mass is larger than the required mass change ) (see (Eq. A4) in the Appendix A), which, however, only occasionally occurs. 140 When the liquid water available to be frozen becomes small enough, the released latent heat is not sufficient to compensate for the required energy deficit ( ), and then the freezing process stops. Consequently, the model freezes soil water quickly, resulting in an underestimated duration of the soil water phase-change processes and the zero-curtain periods, and also cold-biased winter temperatures (Nicolsky et al., 2007;Yang et al., 2018a).

145
Here, we employed a phase-change efficiency and the temperature of the freezing-point depression to effectively solve the problem of overestimating phase-change rates within the current ELMv1-ECA modelling structure. These modification factors are explained below. The phase-change efficiency, introduced by Le Moigne et al. (2012) and adopted by Masson et al. (2013) and Yang et al. (2018a) , is porosity). The underlying assumption here is that the liquid water of soil resists freezing as the freezing process proceeds and , decreases, analogous to how dry soils resist getting drier due to capillary force. We applied the phasechange efficiency to the initially estimated energy and mass change involved, i.e., and (see (Eq. A4) in the Appendix) when freezing or thawing process occur. 155 As in Nicolsky et al. (2007) and Yang et al. (2018a), the occurrence of a phase-change process is then determined by the temperature of the freezing point depression (i.e., an virtual temperature, see (Eq. A8)) instead of . The virtual freezing point depression temperature is reversely derived from the freezing point temperature-depression equation (Fuchs et al., 1978;Cary and Mayland, 1972). With an upper limit as , the virtual temperature describes the lowest temperature that can 160 hold current liquid water content in the freezing soils. That is, the soil temperature has to be lower than the current virtual temperature to allow the freezing process to occur further.
We describe in detail the revised phase-change scheme in the Appendix A. In short, we improved the phase-change scheme of ELMv1-ECA by incorporating two modifications: 1) applying a phase-change efficiency to implicitly account for the heat 165 compensation/deduction to the system from latent heat released/absorbed by soil water freezing/melting, and 2) replacing the constant freezing point with the temperature of the freezing point depression, as a virtual temperature, to determine the occurrence of phase change in subfreezing soils.

Environmental Modifiers to the Decomposition Rate
We revisited ELMv1-ECA's representation for soil heterotrophic respiration dynamics in subfreezing soils and then 170 scrutinized the environmental scalars of soil temperature and moisture. Within ELMv1-ECA's decomposition cascade model, the environmental factors that impact the decomposition rates of soil organic matter include soil temperature ( ), soil moisture ( ), oxygen stress ( ) and a depth scalar ( ) (See Appendix B). Within freezing and subfreezing soils, the soil water potential is related to temperature through the freezing point depression equation (Niu and Yang, 2006). The current moisture factor , therefore, can predict zero respiration rates for subfreezing soils given a minimum soil water 175 potential , as shown by Figure S1a in the supplementary file. (Figures numbered with a prefix "S" are include in the supplementary file, which will not be repeated in the following context throughout the manuscript.) We thus imposed a minimum threshold ( _ ) to prevent zero respiration within the active layer when soil becomes subfreezing during coldseason months ( Figure S1b).

180
For wet soils, the factor that primarily limits the decomposition rates is oxygen availability (Sierra et al., 2017), since increases in soil moisture result in decreased dissolved oxygen. ELMv1-ECA approximates oxygen stress ( ) as a ratio of available oxygen to the demand by decomposers, which, however, is highly uncertain and unstable (Oleson et al., 2013).
Adapting the concept and formulation of Yan et al. (2018), we incorporated oxygen stress into the moisture scalar to account for the inhibition of decomposition in wet anoxic conditions. The revised form of the moisture scalar * (Eq. B11) 185 gradually decreases when the degree of saturation exceeds an optimal wetness threshold ( ) that represents the most favorable soil moisture condition for decomposition, as shown by Figure S1b. We also tested several modified moisture schemes with different shape parameters ( in Eq. B11) and optimal wetness thresholds and minimum thresholds ( and _ in Eq. B11). When using the modified moisture scalar with the built-in oxygen stress, the total environmental impacts on decomposition, i.e., will be modified accordingly as to avoid double-counting of the 190 oxygen stress.
ELMv1-ECA uses a Q10-based standard exponential function to account for the temperature effect on SOC decomposition (Eq. B9), with as 1.5 and as 25°C. Here, rather than striving for a single value of Q10, or a spatial map of Q10 as discussed in the introduction, we seek an optimized scheme at the site scale and a generic scheme at the regional scale for the 195 total environmental modifier ( ) that combines moisture ( ) and temperature ( ) sensitivity. Specifically, we assembled and tested 200 cases of using the newly modified moisture scalars with different parameters , , and _ , temperature scalars with different values of Q10 and , and a variety of other empirical moisture and temperature functions, as documented by Sierra et al. (2015). A full list of the specific moisture and temperature scalars used is provided in Table S2. 200

Cold-season Methane Process
The ELMv1-ECA methane model solves the reaction and diffusion equation for CH4 and O2 fluxes with the Crank-Nicholson method. It includes the representations of CH4 production, oxidation, and three pathways of transport, including aerenchyma tissues, ebullition, aqueous and gaseous diffusion (Riley et al. (2011)). A short description of the ELMV1-ECA methane module is provided in Appendix C. The ELMv1-ECA methane model has been found to underestimate cold-season 205 methane emissions over northern wetlands (Xu et al., 2016). The modifications to the phase-change scheme impact simulations of soil water and heat transfer (3.1.1); the changes in environmental scaler affect substrate availability (3.1.2).
Both (3.1.1) and (3.1.2) influence carbon decomposition and soil heterotrophic respiration (Figure 2), and could potentially lead to improvements in simulated CO2 and CH4 production, but not necessarily CH4 emissions which are also controlled by transport mechanisms (Figure 2). Thus, we further refined the cold-season methane transport processes. 210 Here, we first modified the ELMV1-ECA CH4 transport mechanism in cold seasons by mimicking possible pathways for CH4 emissions from freezing and subfreezing soils. Specifically, we mimic the emissions from ice cracks by plant aerenchyma transport (Zona et al., 2016), approximating the gas diffusion through ice cracks to the similar mechanism of diffusion through the aerenchyma tissues. Although in-situ experiments demonstrated that during winter, produced CH4 in frozen soils is predominately emitted to the atmosphere through vascular plants aerenchyma tissues (e.g., Kim et al., 2007), here we integrate the possible transport pathways including ice cracks and remnants of aerenchyma tissues together through equation (Eq. C14).
During the cold season over the tundra ecosystem, snow on the land surface provides strong resistance to CH4 transport to 220 the atmosphere in ELMv1-ECA, as shown in Figure 2. But in reality, studies have shown methane can diffuse through snowpack at varying rates (Kim et al., 2007). We thus decreased snow resistance at the upper boundary by introducing a new scale factor when snow is present. Also, in ELMv1-ECA, the aqueous diffusion coefficients in freezing and subfreezing soils below the water table are based on the volumetric fraction of the liquid water content, which is quite small (i.e., the supercooled liquid water) and thus limits diffusion. We revised the formulation (Eq. C15), assuming a higher scaling factor 225 for frozen soils ( ) upon sensitivity experiments (not shown). Table 1 summarized all the specific modifications made to ELMv1-ECA. These modifications involve new parameters that are all tuneable and can be systematically optimized via calibration. Here, we seek to reproduce the first-order cold-season process relevant to this study with these default formation and values listed in Table 1.

Climate Forcing, Model Configuration, and Experiment Design 230
We conducted transient simulations at 30-minute temporal resolution driven by climate forcing from 0.5°×0.5° CRU JRA (Harris, 2019) from 1901 to 2017 at the four Alaska tundra site locations. Before the transient simulation, we conducted a 200-year Accelerated Decomposition (AD) spin-up period followed by a 200-year regular spin-up period (Koven et al., 2013b;Zhu et al., 2019) to initialize land carbon pools. Spin-up simulations start from a wet and cold condition. Specifically, sub-surface temperatures were initialized as 274 K for the 1 st to 5 th soil layers, 273 K for the 6 th to 10 th layer, and 272 K for 235 the 11 th to the 15 th layer, and volumetric soil water content was initialized fully saturated for all layers. In this manner, consistent vertical soil water content profiles were built in over the permafrost regions.
Baseline simulations were conducted with ELMv1-ECA default physics, parameters, and surface datasets, i.e., OriPC_OriDecom_OriCH4 using original phase-change scheme, original decomposition scheme and methane module (Table  240 2). To improve the model representation of the site-level soil environment, we first examined the global soil organic matter data at the ABoVE sites by evaluating ELMv1-ECA simulated subsurface soil temperature with the topsoil temperature prescribed to observations (as did in Tao et al., 2017). Using the top soil layer as the upper boundary, the modelling system excluded potential errors induced by inaccurate meteorological forcing and vegetation cover that impact the simulation of heat transfer from the atmosphere to the shallow soil (Tao et al., 2017). Then, the accuracy of simulated soil subsurface 245 temperature is directly determined by the factors impacting heat transfer along the "shallow-to-deep soil" gradient (Koven et al., 2013a), e.g., soil thermal properties which are mostly determined by SOC content (Tao et al., 2017;Lawrence and Slater, 2008). Results well reproduced the subsurface soil temperatures except at IVO, where summer soil temperatures were https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License.
notably overestimated (see Figure S2a). This result indicates that the SOC content at IVO was too small, leading to a large thermal conductivity, small soil porosity, and small heat capacity, altogether resulting in fast penetration of heat into the 250 subsurface soil during summer (Tao et al., 2017;Lawrence and Slater, 2008). Thus, we derived the organic matter density at IVO based on ABoVE soil moisture data through a linear relationship between SOC content and soil porosity (i.e., Equation 3 in Lawrence and Slater (2008)), assuming the observed maximum volumetric water content was porosity (see Figure S3 for details). With the newly derived profile of soil organic matter density at IVO, the simulation showed large improvements in summer soil temperatures compared to that using the original global carbon dataset (see Figure S2b). The derived SOC 255 content is also consistent with the soil survey data reported in Davidson and Zona (2018). Hereafter, the simulations at IVO presented in this paper use the newly derived organic carbon data without repeated clarification.
The representative spatial scale of the eddy flux tower is small compared to the grid cell of global surface datasets and the climate forcing data used by ELMv1-ECA, although the forcing dataset was interpolated to the site scale with a bilinear or 260 nearest-neighbor method. The site-scale vegetation cover also shows a large diversity of vegetation types according to the detailed vegetation survey at ABoVE flux tower footprints obtained in 2014 (Davidson and Zona, 2018). We analyzed the vegetation composition from the closet survey plot to the flux tower and examined the rationality of ELMv1-ECA's percentage of plant type function (PFT) for the site-scale simulation. We confirmed that ELMv1-ECA's PFT dataset was a good compromise between representing the site-scale ecosystem and other global parameters and surface datasets within 265 ELM. The simulated saturated and unsaturated CH4 emissions were weighted with the estimated inundation fractions at the footprint of ABoVE eddy-covariance flux towers (see details in (Xu et al., 2016)) in order to compare simulated CH4 emissions with ABoVE measurements at the site scale. Table 2 lists the experiments conducted in this study. We modified each model component (i.e., the heat transfer model, 270 carbon decomposition model, and methane model) serially. For the temperature-and moisture-dependency functions, we analyzed 200 environmental modifiers within the carbon decomposition module and identified an optimal scheme for each site and a generic scheme that can be applied for the regional simulation over Alaskan North Slope tundra (see next section).

Evaluation Method and Trend Analysis 275
We define the early cold season as September and October, the cold-season period as September to May which includes the two shoulder seasons (both thawing and freezing) as consistent with Zona et al. (2016), and the warm season from June to August. We define the zero-curtain period (ZCP) as the set of successive days when the soil temperature is within the range of [-0.75°C, 0.75°C] starting in fall (i.e., the freezing season) based on Zona et al. (2016). We computed the ZCP duration for each soil layer every year from 1950 to 2017 and estimated the historical trend as the regression slope between ZCP 280 duration and time. Similarly, we estimated the trends of cold-season CH4 and CO2 emissions through linear regression https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License. analysis. A p-value of 0.05 is used to determine if the computed trend is statistically significant. Results vary with soil depths; thus, we choose a common modelling depth, i.e., 12 cm, which locates within the active layer for all the sites, to give an example.

285
We used Nash-Sutcliffe Efficiency (NSE) (Nash and Sutcliffe, 1970) to examine the performance of the ELMv1-ECA simulated time series of CH4 and CO2 net fluxes in comparison with assembled observations (Section 2) at the monthly time scale. The NSE ranges from negative infinity to one, calculated as Eq. (1): where t means monthly time step, N is the total number of time steps, and is simulated and observed flux at time step t, respectively; and is the standard deviation of observations. Note we only used observed monthly averages when the 290 number of daily observations was more than 20 days. The model performance is generally considered satisfactory with an NSE > 0.50 (Moriasi et al., 2007), and perfect with an NSE as one. To simultaneously evaluate CH4 and CO2 fluxes, we combined both and in the form of 1 1 , representing the distance from ( , ) to (1, 1) in a coordinate plane with x-axis as and y-axis as . The optimal simulation thereby is the one having the shortest distance to the ideal scenario (1, 1). We also define a satisfactory model performance 295 in terms of simulating CH4 and CO2 fluxes as the case with both and larger than 0.5. The generic scheme then is the common satisfactory scheme that provides the best overall performance for all the sites.
To evaluate ELMv1-ECA simulated soil temperature and moisture, we calculated the RMSE for each soil layer, i.e., ∑ ⁄ where the and is simulated and observed soil temperature or moisture, respectively, and t is a 300 daily time step. We used the Mean Absolute Error (MAE, . . , ∑ to assess the simulated duration of ZCP of each soil layer. Note that, depending on the amount of soil liquid water content, the whole course of the freezing process may or may not entirely fall into the ZCP, i.e., the ending time of ZCP does not necessarily align with the end of the freezing process. The onset of freezing, though, is always later than the starting day of the ZCP, and the main course of the freezing process is still within the ZCP. 305 Here the modelled active layer thickness (ALT), i.e., maximum thaw depth during an annual cycle, is computed as the bottom depth of the deepest thawed soil layer (i.e., with a maximum annual temperature above 0°C) further extended down to the possible non-frozen fraction of the layer below, as in Tao et al. (2019;. We only derived the length of ZCP for soil layers with a maximum annual temperature above 0°C since limited phase-change processes occur in deeper layers. 310 Then, the soil layers containing or below the permafrost table have a zero-day ZCP. We computed the MAE of ALT https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License. simulated with both original (OriPC) and the new phase-change (NewPC) scheme. Also, we computed the relative improvement in simulated soil temperature (Ts) and ZCP compared to the baseline results. Specifically, we calculated 100% × (RMSE_Ts_OriPC -RMSE_Ts_NewPC) / RMSE_Ts_OriPC and 100% × (MAE_ZCP_OriPC -MAE_ZCP_NewPC) / MAE_ZCP_OriPC to quantify the enhancement by employing the new phase-change scheme. 315 In general, we use NSE to evaluate the model's performance in capturing seasonality (i.e., time series of CH4 and CO2 net fluxes) and use RMSE and MAE to assess the model's capability in simulating the magnitudes of soil temperature, moisture saturation, and ZCP durations.

Evaluation of Soil Temperature and Zero-curtain Period
We first evaluated the simulated daily soil temperature profiles against the observations from ABoVE and GIPL-UAF at the four site locations. Then, we examined improvements in simulations of soil temperature, soil moisture, and the durations of ZCPs by employing the newly revised phase-change scheme (i.e., "NewPC_OriDecom_OriCH4"; Table 2).

325
Results for the BES/CMDL and IVO site are shown in Figure 3; results for other sites are shown in supplementary Figure   S4. At BES/CMDL, the baseline (i.e., "OriPC_OriDecom_OriCH4"; Table 2) simulated soil temperatures (Ts) with the default phase-change scheme (Ts_OriPC; cyan lines; Figure 3a) decrease rapidly in fall due to the overestimated freezing rate (i.e., the slope of decreasing liquid water fraction), notably underestimating the duration of the ZCP (greenish shaded area). Consequently, liquid water saturation (Sf_OriPC, green lines; Figure 3a) quickly drops to a lower bound (i.e., the 330 supercooled liquid water content divided by porosity), and the freezing process generally completes within a short period (days for top layers to one month at the most for deeper layers). The baseline model soil temperature drops (Ts_OriPC) sharply after the freezing process ends (i.e., Sf_OriPC decreases to the lower bound). In contrast, the new phase-change scheme effectively slows freezing rates, showing relatively smaller slopes of decreasing liquid water saturation (Sf_NewPC; magenta lines; Figure 3a) within the ZCPs than the baseline simulation (Sf_OriPC; green lines) especially in the 4 th and 5 th 335 layer. Hence, the gradually released latent heat maintains soil temperatures around the freezing point for a longer period (Ts_NewPC; blue lines; Figure 3a), effectively extending the ZCPs (blue shaded area) which agree better with observations than the baseline results. The ZCP duration increases with depth and can extend into December for deep soil layers.
Similarly improved performance was found at the BEO and ATQ sites (supplementary Figure S4). At IVO, however, while the new phase-change scheme greatly improved simulated results relative to the baseline simulation (Figure 3b), the model 340 still slightly underestimated ZCP durations and also underestimated winter (December to April) soil temperature (blue vs. red). This result at IVO is consistent with the underestimation of late-season soil liquid water available to be frozen, and thereby to release sufficient latent heat ( Figure S5 Table 3). For example, at 12 cm depth (4 th layer), the relative improvements in MAE of the ZCP durations were 65%, 65%, 66%, and 50% for the four site locations (Table 3). The largest improvement in MAE was as large as 65 days for the 6 th layer at BES/CMDL, with a relative improvement of 84% (Table 3). This large improvement stems from the better-estimated ALT at this site; the OriPC simulated 6 th layer temperature remained below freezing, leading 350 to a zero-day ZCP (diamonds on the x-axis in Figure 4). The new phase-change scheme not only improved simulation of the ZCP and cold-season soil temperatures, but also affected the warm season dynamics and thus ALT estimates. As Figure 4 indicates, the NewPC improved simulated ALTs at all four site locations with reduced bias in multi-year averaged ALT, resulting in more reasonable ZCP durations for the 6 th layer (and also 7 th layer for IVO), while the baseline results were zero days. 355 The deeper active layer simulated by NewPC implies more soil water storage capacity, resulting in lower soil moisture in shallow soil layers and higher soil water in deep layers (Sf_NewPC; magenta lines; Figure 3) compared to baseline results.
The changes in soil liquid water content, in turn, impact phase-change and soil temperature simulations. Comparison with the observed soil liquid water content reveals a better agreement with observations (Table S3). For example, at ATQ ( Figure  360 S6), the RMSEs of the liquid water content were reduced by 5.4%, 35.3%, 42.6%, and 25.4% for the 3 rd through 6 th layers, respectively (Table S3).
The changes to model representations of phase change led to large reductions in soil temperature bias. The relative improvements in RMSE of simulated soil temperatures during Sep. and Oct. (i.e., the two months that the ZCPs usually 365 cover), generally increased with depth for surface layers (within about 20 cm of the surface, i.e., 1 st to 4 th layer), and were above 80% for the intermediate layers (5 th to 8 th ) at all the sites ( Figure 5). At the two Barrow sites where observed soil temperatures were available, the relative improvements for the deepest (13 th ) layer were 72.6% and 71.1%, on average, for the early winter and annual cycle, respectively. Therefore, incorporating the new phase-change scheme also resulted in improved bottom temperature boundary conditions, which is critical for accurately simulating permafrost dynamics (Sapriza-370 Azuri et al., 2018). Improvements between Septemper and December and the whole annual cycle also increased with soil depth, showing site-averaged reductions in RMSEs ranging from 47% to 63% and from 36% to 46% for the two periods, respectively. The whole cold-season period (Sep. to May) showed, on average, 44% to 53% reduction in RMSEs from the 1 st to 6 th layer at relatively warmer sites (i.e., ATQ and IVO), and from 19% to 69% for the top 13 layers for the two Barrow sites. 375 https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License.
Soil temperatures were still slightly underestimated during the thawing season (i.e., May) at all four sites, showing later onset of thawing indicated by the timing when warming soil temperatures cross 0°C and soil moisture starts to rise ( Figure   3). One possible reason for this bias is the lack of representation of advective heat transport. That is, the model does not represent the heat of spring rain that is advectively transported into soils (Neumann et al., 2019;Mekonnen et al., 2020); nor 380 does it account for advective heat transport associated with water fluxes in subsurface soils after the spring-rainwater mix with existing cold liquid water in soils. Also, after the freezing process ends, simulated deeper soil layer temperatures were underestimated (e.g., December through April). This bias might be caused by underestimated snow depth (not shown) resulting from inaccurate forcing (particularly snowfall), land cover, microtopography, and/or wind-blown snow redistribution. 385 The improved simulations of soil temperature, liquid water content, and ZCP duration greatly impacted soil HR and methane production ( Figure 6). Increases in the baseline modeled HR and CH4 production resulted from changes in soil temperature and moisture (Figure 6b1 and b2 vs. Figure 6c1 and c2) and mainly occured within the two-dimensional "zero-curtain zone" across the vertical soil profile spanning multiple months, i.e., from September to November (Figure 6c1). However, still very 390 small HR and CO2 and CH4 production were predicted during the following cold season months (Figure 6c3 and c4) due to the moisture scalar for subfreezing soils estimated by ELMv1-ECA's original moisture-dependency function on decomposition (Eq. B10), as discussed in Section 3.1.2. In addition, the sharp decreases of HR and CH4 production around the end of September were caused by the dramatically increased oxygen stress (i.e., the decreased oxygen scalar) to decomposition when freezing began (Figure 6c3 and c4). By replacing the original moisture scalar with the modified soil 395 moisture-dependency function scheme-2 with oxygen stress ((Eq. B11), also see Figure S1) along with the modified total environmental modifier, both the near-zero respiration and sharp drawdown trends in HR and CO2 and CH4 production were greatly alleviated (Figure 6c3 and c4 vs. Figure 6d3 and d4). In the next section, we analyze 200 environmental modifier schemes to the base decomposition rate (Table S2) that assembled commonly used empirical soil temperature-and moisturedependency functions as documented by Sierra et al. (2015) and the modified functions proposed in this study. 400

Evaluation of CO2 and CH4 Emissions
Here we evaluate the simulated monthly CO2 and CH4 fluxes at the site scale against EC tower observations. Figure 7 displays the NSEs of 200 ELMV1-ECA ensemble simulations with different combination of temperature and moisture scalers on soil decomposition, i.e., "NewPC_NewDecom_NewCH4" (grey dots) (see configurations in Table 2). (Daily time  405 series of all the simulations are provided in Figure S7). The failure of simulated CH4 emissions to capture the methane seasonality at IVO (as indicated by Figure S7) might occur because of the lack of 1) a reasonable wetland module that can adequately account for inundated hydro-ecological dynamics, 2) advective heat transport at the air-ground interface through rainfall infiltration and within subsurface soils through water transfer, and 3) the geological micro-seepage emission of CH4, https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License. as reported in previous studies (Anthony et al., 2012;Etiope and Klusman, 2010;Russell et al., 2020). For instance, Lyman 410 et al. (2020) showed large temporal variability of CH4 at natural gas well pad soils, similar to the observations at IVO (Anthony et al., 2012). Controlled experiments (not shown) that imposed observed soil temperature and moisture into the modelling system at all the layers with observations available do not demonstrate improvement for the simulation of CH4, although showing better performance for CO2. These results confirm that impacts from the soil environment (e.g., soil temperature and moisture) within the current water and heat transfer framework cannot explain the seasonal variability of 415 CH4 emissions. Thus, the three mechanisms discussed above (i.e., wetland dynamics, advective heat transport, and geological micro-seepage CH4 emission) currently missing in our model are likely necessary to simulate CH4 emissions at this site and we therefore do not include analysis at IVO in the following sections.
The improved phase-change scheme, and thus improved simulations of ZCP durations and soil temperature and moisture, 420 resulted in greatly improved performance for CO2 emissions at BES/CMDL and BEO, and slightly better performance for CH4 emissions at ATQ, compared to the baseline (cyan for "NewPC_OriDecom_OriCH4" vs. green for baseline; Figure 7), even though the carbon decomposition and methane modules remained the same. Incorporating the revised CH4 model (discussed in section 3.1.3) improved simulated CH4 emissions at BES/CMDL, BEO, and ATQ (blue for "NewPC_OriDecom_NewCH4" vs. cyan for "NewPC_OriDecom_OriCH4"), especially during the cold season (Figure 8). 425 The improved NSEs for CH4 emissions mainly resulted from increased emissions over early winter (Sep. and Oct.) and slight but persistent enhancements throughout the rest of the cold season (blue in Figure 8), which were related to our modifications to CH4 transport mechanisms. Further, with the identified optimal scheme of environmental modifiers to the base decomposition rate, results demonstrate substantial improvements to the simulation of cold-season CO2 and CH4 emissions compared to baseline results (yellow vs. others; i.e., shortest distance from ( , ) to (1, 1)). Among 430 the common schemes providing good performance for both CO2 and CH4 emissions (i.e., both and larger than 0.5, indicated by the gray dots within the boxes in Figure 7), we identified a generic scheme of environmental modifier to the decomposition rate by selecting the common scheme that provided best overall performance for all the sites (except IVO). The specific functions for the optimal and generic scheme of environmental modifiers are provided in Table S4. CO2 and CH4 emissions, especially during the warm season when the thaw depth is deep and soil wetness is high, thus permitting large moisture modifier scalar applied to the base decomposition rates.
Both the optimal and the generic decomposition scheme used the modified ELMv1-ECA moisture scalar (see Table S4), which assigns small thresholds for the moisture scalar and also incorporates oxygen stress when soil wetness exceeds a 450 favourable threshold (0.65 here) for decomposition. Imposing small thresholds for moisture scalar effectively prevents the possibility of zero respiration in subfreezing soils during wintertime. This change exerts more impact on cold sites, such as the two Barrow sites, due to the smaller supercooled liquid water under the colder temperature. Thus, the improved NSEs for CO2 and CH4 emissions at BES/CMDL and BEO were larger than those at ATQ (Figure 7; yellow or magenta vs. blue).
Since the temperature at ATQ was not cold enough to make the supercooled liquid water content small enough to give a zero 455 moisture scalar, the microbial respiration was not completely shut down with the original decomposition modifier at this site.
Indeed, at ATQ, where cold-season temperatures are relatively warmer than at BES/CMDL and BEO, simulations with the original ELMv1-ECA environmental modifier (i.e., "NewPC_OriDecom_NewCH4"; discussed in Section 3.1.2), already released much more CO2 and CH4 throughout the cold season than in the baseline simulations, owning to the improved simulations of soil temperature and moisture, and the modifications for CH4 transport. 460 The Q10-based temperature functions mediate the response of microbial respiration more over the warm season than the cold season due to the larger sensitivity of heterotrophic respiration to warm temperatures than to subfreezing temperatures (see Figure S1d). The different SOC decomposition Q10 values employed directly impact soil HR and thus CH4 and CO2 emissions, and also indirectly impact vegetation nutrient assimilation and thus primary productivity (Figure 2). Vegetation 465 growth, on the other hand, impacts CH4 emissions because the CH4 component transported to the surface via vegetation aerenchyma tissue generally dominates the total emissions and thus determines the seasonal peak and general seasonality of CH4 emissions. When temperature is below the reference temperature (i.e., , here is 25°C), a smaller Q10 permits larger HR and produces more CH4 and CO2, increases warm-season CO2 uptake via photosynthesis; and increases belowground biomass and aerenchyma tissue and thereby CH4 transport to the atmosphere. Thus, the seasonality of CH4 and CO2 net 470 emissions are closely linked through vegetation primary productivity, which vary from site to site. For cold sites (i.e., BES/CMDL and BEO), the sensitivity of simulated CH4 to Q10 values is larger than the sensitivity of CO2 net flux to Q10 because cold temperature suppresses vegetation growth (i.e., CO2 uptake); while for the warm site (i.e., ATQ), both CH4 and CO2 net flux are very sensitive to the Q10 values. Summarizing, the cold sites (i.e., BES/CMDL and BEO) better match CO2 and CH4 emissions observations with smaller Q10 values (1.7 or 1.8) than for the warmer site (i.e., 2.1 for ATQ; Table S4). 475 The generic decomposition scheme used a Q10 value of 2.0, which provided the best overall performance at all three sites (Table S4).
The extended ZCPs, the revised environmental modifier to decomposition, and the modified cold-season CH4 transport mechanism, together resulted in the largest improvements for both CO2 and CH4 emissions, especially over the cold season. 480 In the next section, we quantify the cold season contribution of CO2 and CH4 emissions and then estimate the historical trends of seasonal CO2 and CH4 emissions from 1950 to 2017.

Cold-season Contribution of CH4 and CO2 net emissions and Historical Trends
To better verify the cold-season contribution of CH4 and CO2 emissions to the annual budget, a multi-year average approach 485 was taken because of discontinuity in the observed time series. The new simulation results with the optimal decomposition scheme (yellow; Figure 9) showed greatly enhanced performance at three of the study sites in terms of capturing the averaged seasonal cycle, especially for the cold-season months (Sep. to May; Figure 9), reducing site-averaged MAEs in cold-season total CH4 and CO2 emissions by 84% and 81% , respectively. Specifically, compared to baseline results, the new simulation results showed 0.94 gC m -2 and 55.6 gC m -2 increases in site-averaged cold-season CH4 and CO2 emissions, 490 respectively. The observed cold-season CH4 emissions contributed at least ~40% to the annual total at three of the study sites, of which about half occurred in September and October ( Figure 10, Table 5), i.e., the two months hosting the major part of ZCPs for the top to intermediate soil layers. The simulated percentage of cold-season contributions to the annual totals were close to observed values, i.e., 38%, 41%, 28% vs. 45%, 42%, 45% for BES/CMDL, BEO, and ATQ, respectively. The simulated contribution of early cold season (Sep. and Oct.) CH4 emissions to the cold-season total was 495 62%, 52%, and 60% for the three sites, in comparison with the observed 47%, 58%, and 43%, showing slightly overestimations.
The new simulations accurately captured the observed cold-season contributions of both CH4 and CO2 emissions (Table 5), and the model improvements were larger for cold sites (i.e., BES/CMDL and BEO) than for the warmer site (i.e., ATQ), as 500 discussed above. Specifically, at ATQ, despite the small biases in the annual total CH4 emission (i.e., -0.16 gC m -2 ) and the early cold season component (i.e., -0.05 gC m -2 ), the new simulation underestimated the cold-season proportion of annual emissions, i.e., simulated 28% vs. observed 45%. In contrast, biases in contribution percentages were only 2% and 7% at BES/CMDL, and 3% and 1% at BEO for the early cold season and cold-season period, respectively. The updated ELMv1-ECA also provided improved cold-season CO2 emissions, showing small biases of -2.44 gC m -2 (3% of the observation) and 505 -1.5 gC m -2 (2% of the observation) for BES/CMDL and BEO, respectively. Compared to BES/CMDL and BEO, results for ATQ showed relatively larger bias of -23.9 gC m -2 (41% of the observation). The observed multi-year averaged annual CO2 net flux was 19.9 gC m -2 (source), 31.8 gC m -2 (source), and -3.8 gC m -2 (sink) at BES/CMDL, BEO, and ATQ, respectively.
However, due to the large discontinuity in CO2 observations, especially over the warm season (Figure 8), the calculated annual CO2 budget is uncertain. Still, we can characterize the CO2 budget with simulated results using the updated ELMv1-510 https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License.
ECA. We find that the simulated cold-season CO2 emissions were larger than the warm-season CO2 net uptake at all three sites ( Figure 10, Table 5). The released CO2 over the early cold season (September and October) accounted for 54%, 50%, and 72% of the total emissions throughout the cold season for BES/CMDL, BEO, and ATQ, respectively.
Through trend analysis between 1950 and 2017, we found that the ZCP durations showed increasing trends at all three sites, 515 with ZCP trends increasing with depth (Table 6). At ATQ, a warmer site than BES/CMDL and BEO, the trends of ZCP durations increase from 0.14 to 0.49 days yr -1 along the vertical soil profile. The CO2 emissions during the 12 cm ZCP and during cold-season months (September to May) both showed increasing trends at all three sites (Table 7), ranging from 0.19 to 0.26 gC m -2 yr -1 for the 12 cm ZCP, and from 0.33 to 0.38 gC m -2 yr -1 for the entire cold season period. The annual CO2 net flux showed positive trends, but they were not statistically significant. Annual CH4 emissions showed an increasing trend 520 at ATQ with a rate of 10.6 mgC m -2 yr -1 , but not at the two Barrow sites; cold-season CH4 emissions did not show significant trends at all the sites. In a companion paper, we discuss the regional trends of the spatially averaged CO2 emissions simulated by the updated ELMv1-ECA with the identified generic decomposition scheme.

Summary and Discussion
In this study, we improved ELMv1-ECA simulated subsurface soil temperature, zero-curtain period durations, and cold-525 season CH4 and CO2 net emissions at Alaskan North Slope tundra sites. We first improved the numerical representation of coupled water and heat transport with freeze-thaw processes via modifying ELMv1-ECA's phase-change scheme. Then, we revised the dependency of soil decomposition rates on soil temperature and moisture. We further refined the cold-season methane processes by updating upper boundary resistance that allows CH4 to be emitted from frozen soils through snow to the atmosphere. We also used the updated ELMV1-ECA to estimate historical trends of cold-season CH4 and CO2 net 530 emissions at the Alaskan tundra sites from 1950 to 2017. This study is among the first efforts toward improving simulations of zero-curtain periods and cold-season carbon emissions over Arctic tundra by ESMs. The strategy of improving ELMV1-ECA phase-change scheme and environmental controls on microbial activity can be easily applied to other global land models.

535
With the revised phase-change scheme, the updated ELMv1-ECA greatly improved site-scale simulations of soil temperature, soil moisture, and zero-curtain period. Specifically, the RMSE of daily subsurface soil temperature was substantially reduced compared to the baseline simulation, showing site-averaged improvements ranging from 58% to 87% over the early cold season (Sep. to Oct.) and from 36% to 46% over the annual cycle for soil layers within the active layer.
The evaluation of simulated liquid water content with the new phase-change scheme, although limited by availability of 540 observations, showed a relative reduction in RMSE as high as 43% for the 5 th layer at ATQ, and site-averaged improvements of 15% and 21% for the 4 th and 5 th layer, respectively. Simulated ZCP durations were also greatly improved, with, e.g., https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License. relative reductions in MAEs of 65%, 65%, 66%, and 50% for the 4 th layer (about 12 cm) at BES/CMDL, BEO, ATQ, and IVO, respectively.

545
Based upon the improved simulations of soil temperature and moisture with the new phase-change scheme, the identified an optimal SOC decomposition scheme, and the revised methane module, the site-averaged mean annual errors of coldseason emissions were reduced by 84% and 81% for CH4 and CO2, respectively. We also found that CH4 and CO2 emissions over the early cold season, i.e., September and October, which usually accounts for most of the zero-curtain period, contributed more than 50% of the total emissions throughout the cold season (September to May). Zero-curtain period 550 durations showed increasing trends from 1950 to 2017, with larger trends in deeper soil layers. Although the annual CO2 emissions did not show statistically significant trends, both CO2 emissions during the 12 cm depth zero-curtain period and the entire cold-season period (Sep. to May) showed increasing trends.
Although showing improvements compared to baseline results, the new simulations generally overestimated the contribution 555 of early cold season CO2 emissions. Many reasons could contribute to the overestimations, including poor representation of coupled biogeochemical and hydrological processes in the localized permafrost soil environment, the lack of accurate representation of inundated hydro-ecological dynamics, underestimation of snow accumulation due to micro-topographic effects and thus the snow insulation to the ground (e.g., Bisht et al., 2018), among others. Strong microtopographic impacts on CO2 and CH4 emissions across seven landscape types in Barrow, Alaska, were recently reported ; 560 Grant et al., 2017a;Grant et al., 2017b). In addition, the single static multiplicative function ( ) used to parameterize the total impact of environmental conditions on respiration might not be appropriate .
Also, inappropriately prescribed land cover at the site scale or inaccurate climate forcing (particularly air temperature and precipitation; Chang et al. (2019)) could all impact snow accumulation processes (Tao et al., 2017), which can significantly impact CO2 and CH4 emission simulations. Customizing the complex local ecosystem vegetation community might be 565 feasible at the site scale, however, it is less possible for regional or global land model simulations. This issue calls for the importance of upscaling methods to model (e.g., Pau et al., 2016;Liu et al., 2016) and measure (e.g., Natali et al., 2019;Virkkala et al., 2019) carbon and water cycle dynamics at the regional and global scales.
Given the persistent warming and the continued more severe warming in the cold season (Box et al., 2019), we envision 570 continuing increases in cold-season CO2 and CH4 emissions from the permafrost tundra ecosystem. The increasing rate of cold-season heterotrophic respiration (releasing CO2) may become larger than the trend of warm-season vegetation CO2 uptake under future climate. To accurately characterize cold-season emissions and warm-season net uptake, models have to correctly simulate both components, which, however, few models can do. The updated ELMv1-ECA, with the enhanced capacity to reproduce cold-season CO2 and CH4 emissions proven by this study, can serve as a starting point to better predict 575 permafrost carbon responses to future climate. Finally, the complex water-carbon interactions require modelling systems https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License. with fully coupled hydrological-thermal-biogeochemical processes to better predict the carbon budget in permafrost regions under future climate.

Author contributions
JT assembled observations, developed the methodology, conducted model simulations, analysed results, and wrote and revised the paper. QZ contributed to experiment design, editing the original and revised manuscript. WJR and RBN edited 585 the original and revised manuscript, and provided valuable discussion and guidance.

Competing interests
The authors declare that they have no conflict of interest.

Acknowledgments
We are grateful for valuable discussions with Jinyun Tang, Roisin Commane, Xiyan Xu, Kai Yang, and Chenghai Wang. We 590 thank the anonymous reviewers for their helpful comments. This material is based upon work supported by the U.S.

Appendices: Description of Relevant Modules within ELMv1-ECA
Here we describe the heat transfer in subsurface soils, the carbon decomposition, and the methane module within the 595 ELMv1-ECA that are of particular relevance to our study.

Appendix A Subsurface Heat Transfer
ELMv1-ECA approximates the subsurface heat transfer process with a one-dimensional heat diffusion equation: https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License. , (Eq. A1) where is the soil temperature (K), is the volumetric soil heat capacity (J m -3 K -1 ), is soil thermal conductivity (W m -1 K -1 ), and is the soil depth ( where is the freezing temperature of water (0°C in Kelvin, i.e., 273.15 K), , and , is the mass of ice and liquid water (kg m -2 ) of layer , and , , (kg m -2 ) is the supercooled liquid water that is allowed to coexist with ice given the subfreezing soil temperature . This , , varies with soil texture and temperature and is calculated by the freezing point depression equation (Niu and Yang, 2006), where ∆ is the soil thickness of the th layer (in mm), , represents the soil porosity (i.e., the saturated volumetric water 620 content), is the latent heat of fusion (J kg -1 ), is the Clapp and Hornberger exponent (Clapp and Hornberger, 1978), g is the gravitational acceleration (m s -2 ), and , is the soil texture-dependent saturated matric potential (mm).
https://doi.org/10.5194/tc-2020-262 Preprint.   Masson et al. (2013) and Yang et al. (2018a), is calculated as, where , and , is the soil liquid and ice volumetric water content of layer at previous time step , respectively, and The temperature of the freezing point depression, as a virtual temperature ( ) reversely derived from the freezing point temperature-depression equation, i.e., (Fuchs et al., 1978;Cary and Mayland, 1972), is calculated as, * 10 10 , (Eq. A8) where is the latent heat of fusion (J kg -1 ) and g is the gravitational acceleration (m s -2 ). is the soil water potential (mm), calculated as the soil water retention curve of Clapp and Hornberger (1978), i.e., , , , , where 645 , = , ∆ ⁄ as in (Eq. A3), is the Clapp and Hornberger exponent, and , is the soil texture-dependent saturated matric potential (mm).

Appendix B Decomposition Cascade Model
Within the ELMv1-ECA Century decomposition cascade model, the respiration fractions are parameterized as the fraction of 650 the decomposition carbon flux out of each carbon pool, including litter and soil organic matter. The base decomposition rate is modified by a function representing environmental controls on soil decomposition which accounts for the impacts of individual factors including temperature ( ) and moisture ( ), an oxygen scalar ( ), and a depth scalar ( ), in a multiplicative way, i.e., .

655
We use a Q10-based standard exponential function to account for the temperature effect on decomposition, , (Eq. B9) where Q10 = 1.5 on default, which is consistent with ecosystem-level observations (Mahecha et al., 2010), and is the reference temperature (25°C). During cold seasons when soil temperature becomes subfreezing, respiration continues but with more controls from liquid water stress. The original moisture scalar ( ) within ELMV1-ECA is given in the formulation, calculated as, 660 is the soil water potential, where is the Clapp and Hornberger exponent (Clapp and Hornberger, 1978). In frozen soil, the soil water potential is related to soil temperature through the freezing point depression https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License. equation, i.e., (Fuchs et al., 1978;Cary and Mayland, 1972) in the supercooled water formulation (Niu and Yang, 2006). Thus, the liquid water stress on decomposition is translated into dependency on temperature when soil temperature is below the freezing point. 665 ELMv1-ECA approximates oxygen stress ( ) as the ratio of available oxygen to that demanded by decomposers, and has a minimum value of 0.2 (Oleson et al., 2013). As described by section 3.1.2, we now incorporate the oxygen stress into the moisture scalar. The revised moisture scalar * is calculated as below, where is the degree of saturation, calculated as the ratio of soil volumetric liquid water content to porosity . . , is an optimal threshold beyond which the decomposition will be suppressed by oxygen stress, and b is a parameter controlling the shape of the decreasing limb, and _ is a minimum threshold for * . The depth scalar ( ) represents unresolved other depth-dependent processes (e.g., soil microbial dynamics, priming effects, etc.) (Koven et al., 2013b;Lawrence et al., 2015;Koven et al., 2015). Applying the depth scalar to 675 decomposition rates would exponentially decrease the respiration fluxes along with the vertical soil layers. The is the efolding depth for decomposition, and by default is 0.5 m (Oleson et al., 2013).

Appendix C Methane Model
The ELMv1-ECA methane model includes the representations of CH4 production, oxidation, and three pathways of transport 680 (i.e., aerenchyma tissues, ebullition, aqueous and gaseous diffusion), and solves the transient reaction diffusion equation for CH4. ELMv1-ECA estimates CH4 production ( ; mol m −3 s −1 ) in the anaerobic portion of the soil column based on the upland heterotrophic respiration (HR; mol C m −2 s −1 ) from soil and litter, further adjusted by factors representing influence from soil temperature ( ), pH ( ), redox potential ( ), and seasonal inundation condition ( ) (Riley et al., 2011), expressed as, 685 .
The is a fraction of anaerobically mineralized carbon atoms becoming CH4. Detailed explanation on other factors can be found in (Riley et al., 2011). The methane production will be directly impacted by the changes to water and heat transfer (Appendix A) and HR (Appendix B). The ultimately estimated CH4 emissions are also controlled by oxidation, transport mechanisms (i.e., aerenchyma transport, ebullition, and diffusion), and the upper boundary resistance. Detailed descriptions on CH4 oxidation and transport mechanisms are provided in (Riley et al., 2011). Here we modified CH4 transport 690 mechanisms for facilitating reasonable cold-season CH4 emissions.
Vascular plants aerenchyma tissues serve as diffusive pathways for CH4 to transport from soil layer ( , mol m −2 s −1 ) to the atmosphere, calculated as: where and is the gaseous CH4 concentration (mol m −3 ) in soil depth and in the atmosphere, respectively; is the 695 aerodynamic resistance (s m −1 ); is the gas diffusion coefficient (m 2 s −1 ); is aerenchyma porosity (-); is the ratio of root length to vertical depth (i.e., root obliquity); and is the root fraction in soil depth (-). is the specific aerenchyma area (m 2 m -2 ), and is expressed as, where represents the aerenchyma radius (=2.9×10 −3 m); is the annual net primary production (NPP), and is the belowground fraction of current NPP; and the factor 0.22 is the amount of carbon per tiller. We integrate the emissions from 700 ice cracks and remnants of aerenchyma tissues with (Eq. C14) by removing temperature limitation and applying a small during winter time, where now represents areas adding up ice crack fractions and remnants of aerenchyma tissues.
ELMv1-ECA estimates aqueous diffusion below water table as, 705 where is the diffusion coefficient (m 2 s -1 ), is the soil porosity, is a scaling factor for frozen soils, defined as where and is the volume (m 3 m -3 ) of liquid water and ice, respectively. In subfreezing soils when T < 0°, is largely limited by liquid water content. Upon sensitivity experiments, we alleviated this limitation by assuming a half deduction for the diffusion coefficient in saturated, frozen soils, i.e., = 0.5. We also decreased snow resistance by introducing new scale factors (Table 2) which intend to increase the conductance at the upper boundary when snow 710 presents.
# "NewPC_OptimalDecom_NewCH4" is the optimal simulation among the 200 "NewPC_NewDecom_NewCH4" cases at each site.
$"NewPC_GenericDecom_NewCH4" means the simulation with the identified generic scheme that can be applied to regional simulation.
The generic scheme is the common satisfactory scheme that provides the best overall performance for all the sites.

980
No baseline ZCP is shown in the 6 th layer for BES/CMDL and the 7 th layer for IVO because the maximum annual temperature is below 0°C.
https://doi.org/10.5194/tc-2020-262 Preprint. Discussion started: 22 October 2020 c Author(s) 2020. CC BY 4.0 License.     Table 2 for the configuration for each experiment). Gray line represents the ensemble mean of simulations within the good performance zone (as shown in Figure 7) with error bars as the standard deviation and the shaded area indicating the minimum-to-maximum bound. Red open circles are observed monthly averages with the number of daily 1020 observations less than 10 days, which are not used for the computation in Figure 7. Optimal -the best simulation for each site; Generic -the simulation with a common decomposition scheme that provides best overall performance for all the sites and will be applied to regional simulation.