Improved ELMv1-ECA simulations of zero-curtain periods and cold-season CH₄ and CO₂ emissions at Alaskan Arctic tundra sites

Field measurements have shown that cold-season methane (CH₄) and carbon dioxide (CO₂) 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 CH₄ and CO₂ 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 and fill the knowledge gap of how cold-season CH₄ and CO₂ 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 the methane 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 (MAE) in zero-curtain durations at 12 cm depth was reduced by 62 % on average. Furthermore, the MAEs of simulated cold-season carbon emissions at three tundra sites were improved by 72 % and 70 % on average for CH₄ and CO₂, respectively. Overall, CH₄ emitted during the early cold season (September and October), 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 (September to May) in the model, compared with around 49.4 % (43 %–58 %) in observations. From 1950 to 2017, both CO₂ emissions during the zero-curtain period and during the entire cold season showed increasing trends, for example, of 0.17 and 0.36 gC m⁻² yr⁻¹ at Atqasuk. This study highlights the importance of zero-curtain periods in facilitating cold-season CH₄ and CO₂ emissions from tundra ecosystems.

1 Introduction

Cold-season carbon emissions from the Arctic tundra could potentially offset warm-season net carbon uptake under 21st century warming climate (Commane et al., 2017; Oechel et al., 2014, 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 CO₂ losses over Arctic tundra ecosystems (Oechel et al., 2014; Natali et al., 2019). Also, CH₄ emitted from September to May was found to contribute more than 50 % of the annual total CH₄ 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 changing climate, prevailing earth system land models do not accurately reproduce cold-season CH₄ and CO₂ 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; Wang et al., 2019; 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).

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., 2015, 2017; Yan et al., 2018; Moyano et al., 2013; Tang and Riley, 2019; Rafique et al., 2016; Bhanja and Wang, 2020; Kim et al., 2019). Specifically, the temperature sensitivities of soil carbon decomposition is often represented with a Q₁₀ 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 Q₁₀ are controversial (Davidson and Janssens, 2006). Some studies found a uniform Q₁₀ across biomes and climate zones, e.g., as 1.4 (Mahecha et al., 2010). Other studies demonstrated that Q₁₀ varies with environmental conditions, ecosystem types, and soil texture (Meyer et al., 2018; Graf et al., 2011; Kim et al., 2019), showing 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 Q₁₀ 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, 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 modeling studies have closely investigated the HR–moisture relationships in frozen conditions.

The strong dependency of CO₂ and CH₄ emission on soil temperature and moisture in ESM land models (Riley et al., 2011; Koven et al., 2017; Lawrence et al., 2015) requires accurate estimates of these two closely related soil variables, especially in cold regions where both increases and decreases in soil temperature could lead to soil “drying” due to drainage or freezing processes. However, current land models tend to significantly underestimate soil temperature during the cold season over permafrost regions (Dankers et al., 2011; Tao et al., 2017; Nicolsky et al., 2007; Yang et al., 2018b). One possible reason is that while many land models account for latent heat released during soil water freezing, they do not treat and distribute this heat appropriately and/or do not simulate soil moisture correctly (Yang et al., 2018a; Nicolsky et al., 2007). Latent heat released during freezing might be sufficient to offset heat conduction towards the surface, thus maintaining the subsurface soil temperature around the freezing point (i.e., 0 ^∘C) for weeks or even months during the fall (i.e., the so-called zero-curtain period; ZCP) (Outcalt et al., 1990). The ZCP conditions allow for continued soil heterotrophic respiration at notable rates and thus CO₂ and CH₄ production and emissions from subsurface soils (Kittler et al., 2017; Arndt et al., 2019; Commane et al., 2017). For instance, Zona et al. (2016) reported that a substantial portion of cold-season CH₄ emissions occurred during the ZCP from Alaskan upland tundra sites. Nevertheless, many land models cannot accurately capture the ZCP length due to inaccurately simulated soil moisture and/or inadequate representation of latent heat, thus underestimating soil temperature and cold-season emissions of CO₂ (Commane et al., 2017) and CH₄ (Zona et al., 2016). We note that snow representation can also play a major role in correctly simulating winter soil temperatures (Slater et al., 2017; Lawrence and Slater, 2010), although we do not focus on this process here.

We hypothesize that the underestimation of modeled cold-season CO₂ and CH₄ emissions in ESMs land models primarily results from underestimated soil temperatures during the cold season, the poor representations of environmental controls on heterotrophic respiration in subfreezing soils, and the lack of appropriate representation of cold-season methane transport processes. Here we apply the Energy Exascale Earth System Model (E3SM) land model version 1 equilibrium chemistry approximation configuration (ELMv1-ECA) (Golaz et al., 2019; Zhu et al., 2019; Burrows et al., 2020) to explore these hypotheses. We apply ELMv1-ECA to (i) improve simulations of subsurface soil temperatures, ZCPs, and CO₂ and CH₄ emissions over the permafrost tundra ecosystem; (ii) investigate the underlying processes that influence cold-season carbon emissions from freezing and subfreezing soils, including source characterization and transport pathways; and (iii) estimate historical trends (from 1950 to present) of cold-season CO₂ and CH₄ emissions at multiple Alaskan tundra sites.

The paper is organized as follows. (1) We describe the study sites and the data used in the study. (2) We present the theoretical background of essential modules of ELMv1-ECA relevant to this study and our modifications to the model's representations of phase-change, SOC decomposition, and methane dynamics. (3) We then describe the model configuration and experimental design. (4) We assess the modified phase-change scheme by comparing simulated soil temperatures and ZCPs against observations. (5) With the revised phase-change scheme, we analyze how the parameterization of decomposition schemes and methane module impact simulated CO₂ and CH₄ emissions at the site scale. (6) Finally, we summarize the main findings and discuss needed observations and model development to further improve predictability.

2 Study sites and data

We assembled daily observations of CO₂ and CH₄ fluxes from 2013 to 2017 at five eddy-covariance flux tower sites in Alaska's North Slope tundra (Fig. 1) from the Arctic-Boreal Vulnerability Experiment (ABoVE) project (2015–2017) (Oechel and Kalhori, 2018) and CH₄ fluxes from the Carbon in Arctic Reservoirs Vulnerability Experiment (CARVE) flight campaign (2013–2014) (Zona et al., 2015, 2016). The CARVE CO₂ measurements were not available at the data archive we used here; therefore, monthly wintertime CO₂ flux data from 2013 to 2014 at the same towers assembled by Natali et al. (2019) are included to complete CO₂ observations. The five sites include three eddy covariance (EC) towers at Barrow (i.e., the Barrow Environmental Observatory (BEO) tower, the Biocomplexity Experiment South (BES) tower, and the Climate Monitoring and Diagnostics Laboratory (CMDL) tower), one tower at Atqasuk (ATQ), and another tower 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 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. 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 included in the Supplement.)

Figure 1 Red dots indicate the five ABoVE flux tower sites used in this study. Cyan circles are GIPL-UAF permafrost sites.

ABoVE and CARVE provide soil temperature and moisture measurements at various depths from 5 to 40 cm. The 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, last access: 19 November 2021) (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 scale and 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.

The observed soil moisture is only available at two or three depths that are quite different from model layer node depths and also show discontinuities in time. Thus, evaluation of ELMv1-ECA simulated liquid water content 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 observed maximum soil moisture to infer site-scale soil porosity and then organic carbon content at IVO (see Sect. 3.2), which is used to prescribe thermal and hydraulic soil properties. Note that carbon substrate for respiration is simulated dynamically in the model (see Appendix B).

3 Methodology

3.1 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). In the appendix, we describe in detail its subsurface soil 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 CH₄ and CO₂ emissions and summarize the modifications made to ELMv1-ECA, emphasizing the model enhancements.

3.1.1 Phase-change scheme

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; T_f), adjusts the soil liquid and ice content by mass and energy conservation, and then readjusts temperatures after accounting for the heat deduction or compensation resulting from melting or freezing (see the detailed description in the Appendix A). The underlying assumption here is, taking the freezing process as an example, that the available liquid water at the initially solved temperature ($T_{\mathrm{i}}^{n+\mathrm{1}}$) will be completely frozen, releasing latent heat (H_i) to bring up $T_{\mathrm{i}}^{n+\mathrm{1}}$ back to T_f. Then, the estimated phase-change rate will be tuned down and the current temperature (i.e., T_f) will be readjusted if the to-be-increased ice mass is larger than the required mass change (−H_m) (see Eq. A4 in the Appendix A), which, however, only occasionally occurs. When the liquid water available to be frozen decreases enough, the released latent heat is not sufficient to compensate for the required energy deficit ($T_{\mathrm{f}}-T_{\mathrm{i}}^{n+\mathrm{1}}$), 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, as well as cold-biased winter temperatures (Nicolsky et al., 2007; Yang et al., 2018a).

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 modeling 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), introduces the dependency of available liquid water on the phase-change rate (Le Moigne et al., 2012). The phase-change efficiency for freezing, ${\mathit{\epsilon }}_{\mathrm{liq},\mathrm{i}}^n$ (see Eq. A9), is identical to the degree of moisture saturation or the volumetric fraction of soil liquid water content (i.e., ${\mathrm{Sf}}_{\mathrm{liq},\mathrm{i}}^n={\mathit{\theta }}_{\mathrm{liq},\mathrm{i}}^n/{\mathit{\theta }}_{\mathrm{sat},\mathrm{i}}$, where ${\mathit{\theta }}_{\mathrm{liq},\mathrm{i}}^n$ is soil liquid water content and θ_sat, i is porosity). We applied the phase-change efficiency to the initially estimated energy and mass change involved, i.e., H_i and H_m (see Eq. A4 in the Appendix) when freezing or thawing processes occur.

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. A10) instead of T_f. 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 T_f, the virtual temperature describes the lowest temperature that can 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 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.

3.1.2 Environmental modifiers to the decomposition rate

We revisited ELMv1-ECA's representation for soil heterotrophic respiration dynamics in subfreezing soils and then 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 (f_T), soil moisture (f_W), oxygen stress (f_O), and a depth scalar (f_D) (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 f_W, therefore, predicts zero respiration rates for subfreezing soils given a specific lower limit of soil water potential ψ_min (−10 MPa; Eq. B2) (Oleson et al., 2013), as shown by Fig. S1a in the Supplement. We thus decreased the ψ_min further to prevent zero respiration within the active layer when soil becomes subfreezing during cold-season months (Fig. S1b) as long as the soil water potential ψ_i exceeds the prescribed ψ_min.

For wet soils, the factor that primarily limits the decomposition rates is oxygen availability (Sierra et al., 2017; Yan et al., 2018), since increases in soil moisture result in decreased dissolved oxygen. ELMv1-ECA approximates oxygen stress (f_O) as a ratio of available oxygen to the demand by decomposers, which, however, is highly uncertain and unstable (Oleson et al., 2013). Some existing moisture scalars incorporate the oxygen stress together to account for the inhibition of decomposition in wet anoxic conditions, e.g., a moisture function proposed by Yan et al. (2018) and several functions tested in Sierra et al. (2015), including STANDCARB (Harmon and Domingo, 2001), Daycent (Kelly et al., 2000), Skopp (Skopp et al., 1990), and Moyano (Moyano et al., 2013). We thus also tested these existing moisture functions by replacing the original moisture scalar with them in the ELMv1-ECA. Particularly for the moisture function of Yan et al. (2018), we implemented it for each soil layer using the soil properties (i.e., porosity and clay content) of each layer and also tested it with different shape parameters and optimal wetness thresholds. When using the moisture scalars with built-in oxygen stress within ELMv1-ECA, the total environmental impacts on decomposition, $i.e.$, f_total=f_Tf_Wf_Of_D, will be modified accordingly as f_total=f_Tf_Wf_D to avoid double counting of the oxygen stress.

ELMv1-ECA uses a Q₁₀-based standard exponential function to account for the temperature effect on SOC decomposition (Eq. B1), with Q₁₀ as 1.5 and T_ref as 25 ^∘C. Here, rather than striving for a single value of Q₁₀ or a spatial map of Q₁₀ as discussed in the introduction, or a particular individual temperature function, we seek a group of environmental modifiers (f_total) that can correctly represent moisture and temperature sensitivity on heterotrophic respiration. Specifically, we assembled and tested 814 cases of f_total using the newly modified moisture scalars and a variety of other empirical moisture and temperature functions, as documented by Sierra et al. (2015) and Yan et al. (2018). A full list of the specific moisture and temperature scalars tested is provided in Table S2.

3.1.3 Cold-season methane process

The ELMv1-ECA methane model solves the reaction and diffusion equation for CH₄ and O₂ fluxes with the Crank–Nicolson method. It includes the representations of CH₄ production, oxidation, and three pathways of transport, including aerenchyma tissues, ebullition, and 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 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 scalar affect substrate availability (3.1.2). Both (3.1.1) and (3.1.2) influence soil heterotrophic respiration and could potentially lead to improvements in simulated CO₂ and CH₄ production but not necessarily CH₄ emissions, which are also controlled by oxidation and transport mechanisms.

Here, we first modified the ELMV1-ECA CH₄ transport mechanism in cold seasons by mimicking possible pathways for CH₄ emissions from freezing and subfreezing soils. Specifically, we mimicked 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 CH₄ 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 Eq. (C3). Also, during the cold season over the tundra ecosystem, snow on the land surface provides strong resistance to CH₄ transport to the atmosphere in ELMv1-ECA. 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 (Appendix C).

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 the default formation and values listed in Table 1. We also conducted sensitivity tests on three variables three key parameters related to CH₄ oxidation and transport processes and tested seven parameterizations (Table S3). The three CH₄ process-related parameters include two key variables in the original CH₄ model that have been reported to have a large uncertainty (Riley et al., 2011), i.e., $f_{{\mathrm{CH}}_{\mathrm{4}}}$ (a fraction of anaerobically mineralized carbon atoms becoming CH₄; Eq. C1) and R_o, max (the maximum oxidation rate constant; Eq. C4), as well as the newly introduced variable ε_aere (a factor representing remnants of aerenchyma tissues during cold seasons and possible pathways via ice cracks; Eq. C3). The sensitivity tests on CH₄ process-related parameters were applied to the model with identified carbon decomposition schemes that predicted good simulations of CO₂ flux (see Sect. 3.3).

Table 1 Specific modifications made to ELMv1-ECA.

3.2 Climate forcing, model configuration, and experiment design

We conducted transient simulations at 30 min temporal resolution driven by climate forcing from $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ 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, subsurface temperatures were initialized as 274 K for the 1st to 5th soil layers, 273 K for the 6th to 10th layer, and 272 K for the 11th to the 15th 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_OriDecomOriCH4 using original phase-change scheme, original decomposition scheme, and methane module (Table 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 modeling 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 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 notably overestimated (see Fig. 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 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., Eq. 3 in Lawrence and Slater, 2008), assuming the observed maximum volumetric water content was porosity (see Fig. 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 Fig. S2b). The derived SOC content is also consistent with the organic layer thickness 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.

Table 2 List of designed experiments.

^a “NewCH4” uses the newly modified methane module with default parameterization (see Table S3 in the Supplement).
^b “NewDecom” means replacing the original temperature- and moisture-dependency functions on decomposition rates with (814) new functions of environmental modifiers as listed in Table S2 in the Supplement. “NewCH4” here means the newly modified methane module with seven parameterizations (Table S3). The seven CH₄ parameterizations were iteratively applied to the model together with 160 common carbon decomposition schemes that provide good performance for CO₂ flux (i.e., NSE_CO₂ >0.5) at all the sites. This resulted in 1934 “NewPC_NewDecomNewCH4” simulations in total (i.e., $\mathrm{814}+\mathrm{160}\times \mathrm{7}=\mathrm{1934}$).
^c “NewPC_Optimized” is the optimal simulation among the 1934 NewPC_NewDecomNewCH4 cases at each site.
^d “NewPC_Generic” means the simulation with the identified generic parameterization 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.

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 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). The ELMv1-ECA's default plant type function (PFT) dataset was derived from satellite-based data by Lawrence et al. (2007). We analyzed the vegetation composition from the closet survey plot to the flux tower and examined the rationality of ELMv1-ECA's percentage of PFT for the site-scale simulation through testing different PFT datasets derived from this vegetation survey (Davidson and Zona, 2018). We found that these PFT datasets generally are not superior to the original PFT dataset, which generally reproduced satellite-based gross primary production (GPP; Fig. S4). We thus 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 the model. The surface CH₄ emission is a weighted average of simulated saturated and unsaturated components using predicted inundation and non-inundation fractions. To compare simulated CH₄ emissions with ABoVE measurements at the site scale, we use the estimated inundation fractions at the footprint of ABoVE eddy-covariance flux towers (see details in Xu et al., 2016).

Table 2 lists the experiments conducted in this study. We modified each model component (i.e., the heat transfer model, carbon decomposition model, and methane model) serially. All the experiments ran through 1901 to 2017 with spin-up as described earlier, although the evaluation and optimization were conducted only using results from 2013 to 2017. We first ran simulations with the 814 environmental modifiers together with the modified methane model with default parameterization (Table S3). Then, we selected the environmental modifiers that provided satisfactory performance in simulating CO₂ flux, and we repeated simulations with the seven CH₄ parameterizations (Table S3). Among all the simulations results, we identified an optimal simulation for each site (see details in Sect. 3.3).

3.3 Evaluation metrics, optimization method, and trend analysis

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, 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 duration and time. Similarly, we estimated the trends of cold-season CH₄ and CO₂ emissions through linear regression analysis. A p value of 0.05 is used to determine if the computed trend is statistically significant. Results for ZCP duration and trend vary with soil depths; thus, we choose a common modeling depth at which the ZCPs show significant trends for all the sites, to give an example.

To evaluate ELMv1-ECA simulated soil temperature and moisture, we calculated the RMSE for each soil layer, i.e., $\sqrt{{\sum }_{t=\mathrm{1}}^N{\left({\widehat{E}}_{\mathrm{t}}-O_{\mathrm{i}}\right)}^{\mathrm{2}}/N}$, where ${\widehat{E}}_{\mathrm{t}}$ and O_t are simulated and observed soil temperature or moisture, respectively, and t is a daily time step. We used the mean absolute error (MAE), i.e., $\frac{\mathrm{1}}{N}{\sum }_{t=\mathrm{1}}^N\left|{\widehat{E}}_{\mathrm{t}}-O_{\mathrm{i}}\right|$ 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.

Here the modeled 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, 2017). 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. Then, the soil layers containing or below the permafrost table have a zero-day ZCP. We computed the MAE of ALT simulated with both the 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.

We used Nash–Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970) to examine the performance of the ELMv1-ECA simulated time series of CH₄ and net CO₂ fluxes in comparison with assembled observations (Sect. 2) at the monthly timescale. The NSE ranges from negative infinity to one, calculated as Eq. (1):

$$\begin{array}{}\text{(1)}& \mathrm{NSE}=\mathrm{1}-\left({\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{t=\mathrm{1}}^N{\left({\widehat{E}}_{\mathrm{t}}-O_{\mathrm{t}}\right)}^{\mathrm{2}}\right)/{\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}},\end{array}$$

where t means monthly time step, N is the total number of time steps, ${\widehat{E}}_{\mathrm{t}}$ and O_t are simulated and observed flux at time step t, respectively, and σ_o is the standard deviation of observations. Note we only used observed monthly averages when the number of daily observations was more than 20 d. 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 CH₄ and CO₂ fluxes, we combined both ${\mathrm{NSE}}_{{\mathrm{CH}}_{\mathrm{4}}}$ and ${\mathrm{NSE}}_{{\mathrm{CO}}_{\mathrm{2}}}$ in the form of $\mathrm{dist}=\sqrt{(\mathrm{1}-{\mathrm{NSE}}_{{\mathrm{CH}}_{\mathrm{4}}}{)}^{\mathrm{2}}+(\mathrm{1}-{\mathrm{NSE}}_{{\mathrm{CO}}_{\mathrm{2}}}{)}^{\mathrm{2}}}$, representing the distance from (${\mathrm{NSE}}_{{\mathrm{CH}}_{\mathrm{4}}}$, ${\mathrm{NSE}}_{{\mathrm{CO}}_{\mathrm{2}}}$) to (1, 1) in a coordinate plane with x axis as ${\mathrm{NSE}}_{{\mathrm{CH}}_{\mathrm{4}}}$ and y axis as ${\mathrm{NSE}}_{{\mathrm{CO}}_{\mathrm{2}}}$. Then, the optimal simulation thereby is the one having the shortest distance to the ideal scenario (1, 1). We also define a satisfactory model performance in terms of simulating CH₄ and CO₂ fluxes as the case with both ${\mathrm{NSE}}_{{\mathrm{CH}}_{\mathrm{4}}}$ and ${\mathrm{NSE}}_{{\mathrm{CO}}_{\mathrm{2}}}$ larger than 0.5.

We optimized the model simulations through two steps. Specifically, we first evaluated the simulations using (814) environmental modifiers to the base decomposition rate that assembled commonly used empirical soil temperature- and moisture-dependency functions (Table S2). These simulations used the newly modified methane model with the default parameters (Table S3). We selected the common decomposition schemes that provided satisfactory results of CO₂ flux for all the sites (i.e., ${\mathrm{NSE}}_{{\mathrm{CO}}_{\mathrm{2}}}$ >0.5). Then, we iteratively repeated simulations with the common carbon decomposition schemes along with the seven CH₄ parameterizations (Table S3). Among all these simulations (NewPC_NewDecomNewCH4; Table 2), we identified an optimal simulation for each site that has the smallest distance from (${\mathrm{NSE}}_{{\mathrm{CH}}_{\mathrm{4}}}$, ${\mathrm{NSE}}_{{\mathrm{CO}}_{\mathrm{2}}})$ to (1, 1) (i.e., dist); the environmental modifier to the base decomposition rate and the methane parameterization used in the optimal simulation is the optimized parameterization for this site.

Further, among the common parameterizations of environmental modifiers and CH₄ parameterizations that show satisfactory performance both in CH₄ and CO₂ fluxes for all the sites, we identified a generic scheme as the one providing the minimum Euclidean distance in a site-performance space, calculated as $\sqrt{{\sum }_{i=\mathrm{1}}^n{\mathrm{dist}}_{\mathrm{i}}^{\mathrm{2}}}$, where n is number of sites. The generic scheme then is the common satisfactory scheme that provides the best overall performance for all the sites and can be applied for the regional simulation over Alaskan North Slope tundra.

In general, we use NSE to evaluate the model's performance in capturing seasonality (i.e., time series) of CH₄ and net CO₂ fluxes and optimize CH₄ and CO₂ simulations. We use RMSE and MAE to assess the model's capability in simulating the magnitudes of soil temperature, moisture saturation, ZCP durations, and cumulative CH₄ and CO₂ emissions.

4 Results and discussion

4.1 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_OriDecomOriCH4; Table 2).

Results for the BES/CMDL and IVO site are shown in Fig. 3; results for other sites are shown in Fig. S4 in the Supplement. At BES/CMDL, the baseline (i.e., OriPC_OriDecomOriCH4; Table 2) simulated soil temperatures (Ts) with the default phase-change scheme (Ts_OriPC; blue lines; Fig. 2a) 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 (bluish shaded area). Consequently, liquid water saturation (S_f_OriPC, green lines; Fig. 2a) quickly drops to a lower bound (i.e., the supercooled liquid water content divided by porosity), and the freezing process generally completes within a short period (days for top layers to 1 month at the most for deeper layers). The baseline model soil temperature drops (Ts_OriPC) sharply after the freezing process ends (i.e., S_f_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 (S_f_NewPC; magenta lines; Fig. 2a) within the ZCPs than the baseline simulation (S_f_OriPC; green lines) especially in the fourth and fifth layer. Hence, the gradually released latent heat maintains soil temperatures around the freezing point for a longer period (Ts_NewPC; red lines; Fig. 2a), effectively extending the ZCPs (reddish shaded area), which agree better with observations (gray shaded area) 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 (Fig. S4 in the Supplement). At IVO, however, while the new phase-change scheme greatly improved simulated results relative to the baseline simulation (Fig. 2b), the model still slightly underestimated ZCP durations and also underestimated winter (December to April) soil temperature (red vs. black). 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 (Fig. S5). In general, the improvements in ZCP are larger in deeper layers than topsoils, with the top layer showing only marginal improvement.

Figure 2 Comparison of multi-year (2013–2017) averaged daily soil temperatures observed (Ts_Obs, black) and simulated with the original (Ts_OriPC, blue) and improved (Ts_NewPC, red) phase-change schemes at BES/CMDL (a) and IVO (b). Simulated moisture saturation with the original (Sf_OriPC; green) and improved (Sf_NewPC; magenta) schemes is shown on the right y axes. Here, the moisture saturation means soil unfrozen (liquid) water content. The horizontal axes indicate days from July to June, with ticks representing the first day of each month. Shaded areas represent durations of zero-curtain periods observed (ZCP_Obs, gray) and simulated (ZCP_OriPC, blue; ZCP_NewPC, red). No baseline ZCP is shown in the sixth layer for BES/CMDL and the seventh layer for IVO because the maximum annual temperature is below 0 ^∘C.

Figure 3 Comparison between observed and ELMv1-ECA simulated durations of ZCP for the original (OriPC; open diamonds) and improved (NewPC; solid circles) phase-change schemes over four annual cycles (July to June) from 2013 to 2017. “ly” means model layer. Simulated ZCP durations with NewPC demonstrate significant improvements compared to OriPC (solid dots vs. open diamonds), especially for the fourth to the deepest layer above permafrost. Note that a zero-day ZCP means that the maximum daily temperature during an annual cycle is below 0 ^∘C. The pairs of zero vs. non-zero-day ZCP (e.g., OriPC_ly 7 at IVO and OriPC_ly 6 at other sites) indicate that baseline results underestimated ALT. The bias (simulation − observation) of multi-year averaged ALT simulated by the two experiments is provided in each panel.

Simulated ZCP durations with the revised phase-change scheme (NewPC) demonstrated notable improvements over the baseline (original) phase-change scheme (OriPC) (solid circles vs. open diamonds) (Fig. 3), showing greatly reduced mean absolute errors (MAEs) (Table 3). For example, at 12 cm depth (fourth 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 d for the sixth 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 sixth-layer temperature remained below freezing, leading to a zero-day ZCP (diamonds on the x axis in Fig. 3). 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 Fig. 3 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 sixth layer (and also the seventh layer for IVO), while the baseline results were zero days.

Table 3 Mean absolute error (MAE) of simulated ZCP (days) with the original phase-change scheme (Ori_PC) and newly resized phase-change scheme (NewPC), and the relative improvement (%) of using the new phase-change scheme compared to the baseline results, calculated as 100 % × (MAE_ZCP_OriPC − MAE_ZCP_NewPC) $/$ MAE_ZCP_OriPC.

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 (S_f_NewPC; magenta lines; Fig. 2) 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 S5). For example, at ATQ (Fig. S7), the RMSEs of the liquid water content were reduced by 5.4 %, 35.3 %, 42.6 %, and 25.4 % for the third through sixth layers, respectively (Table S5).

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 September and October (i.e., the 2 months that the ZCPs usually cover) generally increased with depth for surface layers (within about 20 cm of the surface, i.e., first to fourth layer) and were above 80 % for the intermediate layers (fifth to eight) at all the sites (Fig. 4). At the two Barrow sites where observed soil temperatures were available, the relative improvements for the deepest (13th) 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-Azuri et al., 2018). Improvements between September 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 (September to May) showed, on average, 44 % to 53 % reduction in RMSEs from the first to sixth layer at relatively warmer sites (i.e., ATQ and IVO) and a reduction from 19 % to 69 % for the top 13 layers for the two Barrow sites. 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 (Fig. S9) possibly resulting from inaccurate forcing (particularly snowfall), land cover, microtopography, and/or windblown snow redistribution.

Figure 4 Relative improvement in the RMSE of simulated soil temperature with the new phase-change scheme (RMSE_Ts_NewPC) compared to that with the original scheme (RMSE_Ts_OriPC), calculated as 100 % × (RMSE_Ts_OriPC − RMSE_Ts_NewPC) $/$ RMSE_Ts_OriPC.

Simulations with the new phase-change scheme also show improved agreements between simulated and observed soil temperatures during the spring thawing season compared to the baseline results (red vs. blue in Fig. 2). Compared to observations, the newly simulated 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 (Fig. 2). 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 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.

The improved simulations of soil temperature, liquid water content, and ZCP duration greatly impacted soil HR and methane production but did not necessarily guarantee improvements in CO₂ and CH₄ emissions. In the next section, we closely evaluate simulated CH₄ and CO₂ fluxes with different parameterizations of environmental modifiers and the modified CH₄ parameterizations as described in Sect. 3.1.

4.2 Evaluation of CO₂ and CH₄ fluxes

Here we evaluate the simulated monthly CO₂ and CH₄ fluxes at the site scale against EC tower observations. Figure 6 displays the NSEs of ELMv1-ECA simulations using different carbon decomposition schemes and CH₄ process-related parameters, i.e., NewPC_NewDecomNewCH4 (gray dots) (see configurations in Table 2). (Time series of all the simulations are provided in Fig. S7.) The failure of simulated CH₄ emissions to capture the methane seasonality at IVO (as indicated by Figs. S8 and S10) 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, (3) a representation of microbial dynamics, and (4) the geological micro-seepage emission of CH₄, as reported in previous studies (Anthony et al., 2012; Etiope and Klusman, 2010; Russell et al., 2020). For instance, Lyman et al. (2020) showed large temporal variability of CH₄ at natural gas well pad soils, similar to the observations at IVO (Anthony et al., 2012). The advective heat transport not only impacts soil temperature, but also affects soil moisture redistribution, substrate availability, and microbial activity. Also, methanogen seasonal dynamics would cause hysteretic effects on CH₄ emission response to soil temperatures (Chang et al., 2020; 2021; Chadburn et al., 2020). In the future, we will incorporate a representation of methanogen seasonal dynamics and simulate microbial population and activity levels to address the hysteresis of CH₄ emissions with temperature. We will also explore more on the contribution of geological micro-seepage emission. The four mechanisms discussed above (i.e., wetland dynamics, advective heat transport, microbial dynamics, and geological micro-seepage CH₄ emission) currently missing in our model are likely necessary to simulate CH₄ emissions at this site, and we therefore do not include CH₄ analysis at IVO in the following sections.

Figure 5 Scatter plot between the Nash–Sutcliffe efficiency (NSE) of simulated monthly CH₄ and CO₂ emissions. An ideal simulation has both NSEs of CH₄ and CO₂ as one (i.e., the upper right corner). The boxes encompass simulations with satisfactory performance (NSE >0.5). Optimized (red) – the best simulation for each site; generic (cyan) – the simulation with a common parameterization of carbon decomposition scheme and CH₄ parameter scheme that provides best overall performance for all the sites. See Table 2 for the configuration for each experiment. The gray dots represent all the tested (1934) simulations indicated in the annotation of Table 2. Symbols outside the plotting ranges indicate poor performance, e.g., (−34.9, −0.3) for baseline at IVO, and thus are not shown in the figure.

Figure 6 Observed and simulated monthly CO₂ (a.1, b.1, c.1, d) and net CH₄ (a.2, b.2, c.2) flux with the baseline model (ELM_Baseline) and the experiments with updated models (see Table 2 for the configuration for each experiment). Shaded gray areas indicate the minimum-to-maximum bound of simulations within the good performance zone (as shown in Fig. 5). Black open circles are observed monthly averages with the number of daily observations less than 10 d, which are not used for the computation in Fig. 5. “ELM_NewPC_Optimized” means the best simulation for each site (red diamonds in Fig. 5). Intermediate results including “ELM_NewPC_OriDecomNewCH4” is included in Fig. S8 in the Supplement.

The improved phase-change scheme, and thus improved simulations of ZCP durations and soil temperature and moisture, resulted in greatly improved performance for CO₂ emissions at BES/CMDL and BEO and in slightly better performance for CH₄ emissions at ATQ, compared to the baseline (blue for NewPC_OriDecomOriCH4 vs. green for baseline; Fig. 5), even though the carbon decomposition and methane model remained the same as the baseline. Incorporating the revised CH₄ model with the default parameter (discussed in Sect. 3.1.3) improved simulated CH₄ emissions at BES/CMDL, BEO, and ATQ (magenta for NewPC_OriDecomNewCH4 vs. blue for NewPC_OriDecomOriCH4), especially during the cold season (Fig. S8). The improved NSEs for CH₄ emissions mainly resulted from increased emissions over early winter (September and October) and slight but persistent enhancements throughout the rest of the cold season (magenta in Fig. S8), which were related to our modifications to CH₄ transport mechanisms. Further, with the identified optimal parameterization of environmental modifiers to the base decomposition rate and methane parameters, results demonstrate substantial improvements to the simulation of net CO₂ flux and CH₄ emissions compared to baseline results (red vs. others; i.e., shortest distance from (${\mathrm{NSE}}_{{\mathrm{CH}}_{\mathrm{4}}}$, ${\mathrm{NSE}}_{{\mathrm{CO}}_{\mathrm{2}}})$ to (1, 1)). Among the 121 common schemes providing good performance for both CO₂ and CH₄ emissions (i.e., both ${\mathrm{NSE}}_{{\mathrm{CH}}_{\mathrm{4}}}$ and ${\mathrm{NSE}}_{{\mathrm{CO}}_{\mathrm{2}}}$ larger than 0.5, indicated by the gray dots within the boxes in Fig. 5), we identified a generic scheme by selecting the common parameterization that provided the best overall performance for all the sites (except IVO) (cyan; Fig. 5). The specific environmental modifier functions and methane parameters for the optimal and generic scheme are provided in Table S6.

Figure S8 illustrates the uncertainty associated with the model representations of environmental influences on heterotrophic respiration and methane parameters. The optimal simulations at the study sites used either the modified ELMv1-ECA moisture scalar or Yanetal (see Table S6), i.e., two groups of moisture-dependency functions implemented for each soil layer. For the Sierra et al. (2015) empirical moisture functions, the influence of liquid moisture content on heterotrophic respiration is uniformly applied to all active soil layers, even though the soil properties (e.g., porosity and saturated soil water potential) are quite different vertically. ELMv1-ECA's moisture scalars (including the original scheme) that use soil water potential, in contrast, reasonably explained the varying influence along with the vertical soil profile (i.e., relationships between soil liquid water content and soil temperature varies with soil clay fraction as demonstrated by Fig. 1 in Niu and Yang, 2006). The Yanetal moisture functions also used soil-layer-dependent porosity and clay content to calculate relevant parameters (Yan et al., 2018). The simulations with moisture functions documented in Sierra et al. (2015) generally overestimated CO₂ and CH₄ emissions, especially during the warm season when the thaw depth is deep and soil wetness is high, thus permitting large moisture modifier scalar to be applied to the base decomposition rates for all the soil layers regardless of soil properties.

Reducing the minimum soil water potential ψ_min for moisture scalar effectively prevents the possibility of zero respiration in subfreezing soils during wintertime (Fig. 6). 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 CO₂ and CH₄ emissions at BES/CMDL and BEO were larger than those at ATQ (Fig. 5). Since the temperature at ATQ was not cold enough to make the supercooled liquid water content small enough to give a zero 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_OriDecomNewCH4 in Fig. S8; discussed in Sect. 3.1.2) already released much more CO₂ and CH₄ throughout the cold season than in the baseline simulations, owning to the improved simulations of soil temperature and moisture, as well as the modifications for CH₄ transport. At IVO, although generally showing low NSEs for CH₄, some new simulations have greatly improved ${\mathrm{NSE}}_{{\mathrm{CO}}_{\mathrm{2}}}$ values that are larger than 0.5 (Fig. 5), compared with −0.3 for baseline. Indeed, the best result at IVO (with a ${\mathrm{NSE}}_{{\mathrm{CO}}_{\mathrm{2}}}=$0.78) significantly improved the simulation of summer CO₂ sink compared to baseline result (Fig. 6).

The optimal simulations used Daycent2 temperature-dependency function at ATQ and Q₁₀-based temperature functions at BES/CMDL and BEO with high Q₁₀ values (e.g., 2.0 and 2.5, respectively) (Table S6), mutually mediating the response of microbial respiration with moisture functions discussed above. At all three sites, the optimized parameterizations used a higher ε_aere (i.e., 0.05; Table S6), representing possible cold-season CH₄ emissions through ice cracks and remnants of aerenchyma tissues. This newly introduced variable is highly uncertain, although it can be calibrated at any other site against cold-season measurements. At BES/CMDL, the optimized parameterization used a decreased maximum CH₄ oxidation rate constant which has been reported to be highly uncertain, especially over high latitudes (Riley et al., 2011). The generic scheme (i.e., carbon decomposition + CH₄ parameters) overlaps with the optimized parameterization at BES/CMDL, which provided the best overall performance at all three sites (Table S6). Despite the small site number and the limited spatial representativeness of each site, the identified generic scheme might be applied to the Alaska North Slope tundra. Nevertheless, the generic scheme might induce uncertainty in simulations and might not be the optimal regional scheme over other ecosystems or given different climate forcing and soil conditions. Still, we conclude that the generic scheme can serve as a reasonable initial scheme for estimating CO₂ and CH₄ emissions over other high-latitude areas (e.g., Fig. S11). In the future, we will explore more sites from newly published CO₂ and CH₄ datasets from pan-Arctic ecosystems, e.g., BAWLD-CH4 (Kuhn et al., 2021) and FLUXNET-CH4 (Delwiche et al., 2021; Knox et al., 2019).

The extended ZCPs, the revised environmental modifier to decomposition, and the modified CH₄ transport mechanism and oxidation parameter together resulted in large improvements for both CO₂ and CH₄ emissions, especially over the cold season. Nevertheless, the optimal simulations still overestimated the contribution of the early-cold-season (September and October) CO₂ emissions at BEO and ATQ (top panel; Fig. 7) and underestimated CH₄ emissions during post-ZCP months (e.g., October to December) (bottom panel; Fig. 7). There are many reasons for the early-cold-season CO₂ overestimations, including model deficiencies, prescribed land parameters, and possibly inaccurate forcing. As for the underestimations of post-ZCP carbon emissions, one critical reason is the lack of sudden bursts of CO₂ and CH₄ within the model, i.e., the gases are pushed out of freezing soils during the freeze-up period (Mastepanov et al., 2008; Pirk et al., 2017). Currently, the ELMv1-ECA mimics this sudden burst mechanism by preventing CO₂ and CH₄ from dissolving in the soil ice fraction (Riley et al., 2011), which could capture some burst emissions (e.g., CH₄ emissions in October and September of 2013 at ATQ; Fig. 6); but it still shows an overall underestimation for sudden-burst emissions, especially at colder sites (e.g., BES/CMDL and BEO; Figs. 6 and 7). We will improve this mechanism in the future by explicitly simulating ice-encroaching soil pores and pushing out gases and liquid water out of the soil matrix. In the next section, we quantify the cold-season contribution of CO₂ and CH₄ emissions and then estimate the historical trends of seasonal CO₂ and CH₄ emissions from 1950 to 2017.

Figure 7 Comparison of multi-year (2013–2017) averaged monthly mean net CO₂ flux (a1, b1, c1, d) and CH₄ emissions (a2, b2, c2) from simulations and measurements at the study sites. The error bars represent standard deviation of monthly mean.

4.3 Cold-season contribution of CH₄ and CO₂ emissions and historical trends

Throughout this section, we only retain and discuss the identified optimal simulation results (i.e., ELM_NewPC_Optimized) for each site. To better verify the cold-season contribution of CH₄ and CO₂ emissions to the annual budget, a multi-year average approach was taken because of discontinuity in the observed time series.

The new simulation results with the optimal parameterization showed greatly enhanced performance in terms of capturing the averaged seasonal cycle (red; Fig. 8), especially for the cold-season months (September to May; Fig. 8), reducing site-averaged MAEs in cold-season total CH₄ and CO₂ emissions by 72 % and 70 % (Fig. 8, Table 4), respectively. Specifically, compared to baseline results which significantly underestimated the cold-season carbon emissions, the optimized simulation results showed 0.79 and 44.0 gC m⁻² increases in site-averaged cold-season CH₄ and CO₂ emissions, respectively. The optimized simulations reduced biases in early-cold-season (cold-season) CH₄ emissions by 80 % (74 %), 86 % (76 %), and 77 % (61 %) for BES/CMDL, BEO, and ATQ (Table 4), respectively. The observed cold-season CH₄ emissions contributed at least ∼40 % to the annual total at three of the study sites, of which about half occurred in early-cold-season months (September and October) (Fig. 8; Table 4), i.e., the 2 months hosting the major part of ZCPs for the top to intermediate soil layers. The simulated contributions of early-cold-season (September and October) CH₄ emissions to the cold-season total were 51 %, 65 %, and 55 % for the three sites, in comparison with the observed 47 %, 58 %, and 43 %, showing slight overestimations. Compared to the baseline-simulated percentage of cold-season contributions to the annual total CH₄ emissions (i.e., only 5 %, 6 %, and 15 %), the optimized simulation showed greatly improved agreements with observed contributions, i.e., 35 %, 35 %, 33 % vs. 45 %, 42 %, 45 % for BES/CMDL, BEO, and ATQ, respectively.

Figure 8 Multi-year (2013–2017) averaged cumulative CH₄ emissions (a2, b2, c2) and net CO₂ flux (a1, b1, c1, d) during the early cold season (September and October), cold-season period (September to May), warm-season period (June to August), and the annual cycle (September to August) at our study sites. Due to the large discontinuity in CO₂ observations, especially over the warm season (shown in Fig. 6), the observed annual CO₂ budget is highly uncertain. Still, the cold-season contributions of both CH₄ and CO₂ emissions are greatly improved by the optimized ELMv1-ECA (i.e., ELM_NewPC_Optimized).

Table 4 Total CH₄ emissions and net CO₂ flux over three seasonal periods, including the early cold season, cold season, and the warm season. “ELM_Optimized” here means ELM_NewPC_Optimized (Table 2). The percentage of each seasonal total CH₄ emissions to the annual total is included in the brackets.

The optimized simulations showed larger improvements in cold-season CO₂ emissions (Fig. 8; Table 4) for cold sites (i.e., BES/CMDL and BEO) than for the warmer site (i.e., ATQ and IVO). Specifically, compared to baseline results, the updated ELMv1-ECA reduced the biases in simulated cold-season CO₂ emissions from −56.1 gC m⁻² (64 % of the observation) to −12.1 gC m⁻² (14 % of the observation) for BES/CMDL and from −65.0 gC m⁻² (68 % of the observation) to −12.4 gC m⁻² (13 % of the observation) for BEO. In contrast, the optimized simulation showed slight overestimations for cold-season CO₂ emissions at ATQ and IVO (Table 4). Nevertheless, the optimized ELMv1-ECA provided greatly improved warm-season net CO₂ flux for all the four sites, reducing biases by 110 %, 78 %, 37 %, and 102 % compared to baseline results at BES/CMDL, BEO, ATQ, and IVO, respectively. Indeed, the updated model switched warm-season net CO₂ flux from baseline-simulated net emissions (positive net CO₂ flux) to net uptake (negative net CO₂ flux) at BES/CMDL, BEO, and IVO, correctly matching with observed warm-season net CO₂ flux (Fig. 8).

The observed multi-year averaged annual net CO₂ flux was 19.9 gC m⁻² (source), 31.8 gC m⁻² (source), −3.8 gC m⁻² (sink), and −16.7 gC m⁻² (sink) at BES/CMDL, BEO, ATQ, and IVO, respectively. However, due to the large discontinuity in CO₂ observations (Fig. 6), the calculated annual CO₂ budget is uncertain. Still, we can characterize the CO₂ budget with simulated results using the updated ELMv1-ECA. We find that the simulated cold-season CO₂ emissions were larger than the warm-season CO₂ net uptake during the analyzing period (2013–2017) at all four sites (Fig. 8), showing annual net CO₂ flux as 1.1 gC m⁻² (source), 36.6 gC m⁻² (source), 36.5 gC m⁻² (source), and 18.2 gC m⁻² (source) at BES/CMDL, BEO, ATQ, and IVO, respectively. The simulated CO₂ emissions over the early cold season (September and October) accounted for 50 %, 56 %, 66 %, and 35 % of the total emissions throughout the cold season for BES/CMDL, BEO, ATQ, and IVO, respectively.

Through trend analysis between 1950 and 2017, we found that the ZCP durations showed increasing trends at all three sites, 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.12 to 0.49 d yr⁻¹ along with the vertical soil profile. At BES/CMDL and BEO, only soil layers at 3 cm and 6 cm show statistically significant increasing trends, ranging from 0.10 to 0.13 d yr⁻¹. The CO₂ emissions during the 6 cm ZCP and during cold-season months (September to May) both showed increasing trends at all three sites (Table 6), ranging from 0.12 to 0.17 gC m⁻² yr⁻¹ for the 6 cm ZCP and from 0.30 to 0.40 gC m⁻² yr⁻¹ for the entire cold-season period. Annual CH₄ emissions showed a nonsignificant increasing trend at ATQ with a rate of 0.52 mgC m⁻² yr⁻¹, but neither annual nor cold-season CH₄ emissions show increasing trends at other sites. In the future, we will examine the generic model parameterization at more sites over the pan-Arctic; we will also optimize regional simulations against spatial datasets of CO₂ and CH₄ upscaled from in situ measurements over pan-Arctic permafrost domain (Natali et al., 2019; Virkkala et al., 2021; Zeng et al., 2020; Peltola et al., 2019) and discuss the uncertainty of estimated trends of the spatially averaged CO₂ and CH₄ emissions associated with snow impact and model parameterizations.

Table 5 Historical trend of ZCP durations (days year⁻¹) for each soil layer from 1950 to 2017. (Trends with p>0.05 are not statistically significant.)

Table 6 Historical trend (1950–2017) in site-scale heterotrophic respiration, CH₄ emission, and CO₂ flux during the ZCP duration at 6 cm (third layer), cold-season months (September–May), and the whole annual cycle (September–August). (Trends with p>0.05 are not statistically significant.)

5 Summary and discussion

In this study, we improved ELMv1-ECA simulated subsurface soil temperature, zero-curtain period durations, and cold-season CH₄ and CO₂ 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 mimicking emission pathways through ice cracks and remnants of aerenchyma tissues, reducing the maximum oxidation rate constant, and reducing upper boundary (snow) resistance that allows CH₄ 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 CH₄ and CO₂ net 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 the Arctic tundra by ESMs. The strategy of improving ELMv1-ECA phase-change scheme, environmental controls on microbial activity, and methane parameterizations can be easily applied to other global land models.

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 (September to October) 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 the availability of observations, showed a relative reduction in RMSE as high as 43 % for the fifth layer at ATQ and site-averaged improvements of 15 % and 21 % for the fourth and fifth layer, respectively. Simulated ZCP durations were also greatly improved, with, e.g., relative reductions in MAEs of 65 %, 65 %, 66 %, and 50 % for the fourth layer (about 12 cm) at BES/CMDL, BEO, ATQ, and IVO, respectively.

Based upon the improved simulations of soil temperature and moisture with the new phase-change scheme, the optimized parameterization for SOC decomposition scheme and the revised methane module, the site-averaged mean annual errors of cold-season  emissions were reduced by 72 % and 70 % for CH₄ and CO₂, respectively. We also found that CH₄ and CO₂ emissions over the early cold season, i.e., September and October, which usually account 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 durations showed increasing trends from 1950 to 2017, with larger trends in deeper soil layers. Also, both CO₂ emissions during the 6 cm depth zero-curtain period and the entire cold-season period (September to May) showed increasing trends. Note that the optimized parameterizations would be biased if there is a bias in simulated soil carbon and therefore should not be taken directly to other models without further analysis. Instead, the optimization procedure described in this study provides a roadmap that can be directly adopted to calibrate other models at different sites.

Although showing improvements compared to baseline results, the new simulations generally overestimated the contribution of the early-cold-season (September and October) CO₂ emissions at BEO and ATQ. 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 CO₂ and CH₄ emissions across seven landscape types in Barrow, Alaska, were recently reported (Wang et al., 2019; Grant et al., 2017a, b). Sensitivity analysis demonstrates large impacts of snow depth on simulated winter soil temperature, summer soil moisture, heterotrophic respiration, and CO₂ fluxes (Fig. S9); therefore, the simulation of snow should be the subject of future investigations.

The underestimated emissions during post-ZCP months (October to November) may be caused by the lack of sudden bursts of CO₂ and CH₄ during the freeze-up period. In addition, the single static multiplicative function used to parameterize the impact of environmental conditions on respiration might not be appropriate because the environmental impact also depends on maximum respiration rate, soil texture, soil carbon content and quality, and microbial biomass (Tang and Riley, 2019). Moreover, due to lacking representations of wetland hydro-ecological dynamics, the model uses simulated upland heterotrophic respiration to estimate CH₄ production (Riley et al., 2011), which might cause underestimations of CH₄ emissions, especially under wet conditions. 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 CO₂ and CH₄ emission simulations. Customizing the complex local ecosystem vegetation community might be 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 continuing increases in cold-season CO₂ and CH₄ emissions from the permafrost tundra ecosystem. The increasing rate of cold-season heterotrophic respiration (releasing CO₂) may become larger than the trend of warm-season vegetation CO₂ 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 CO₂ and CH₄ emissions proven by this study, can serve as a starting point to better predict permafrost carbon responses to future climate (Tao et al., 2021a). Finally, the complex water–carbon interactions require modeling systems with fully coupled hydrological–thermal–biogeochemical processes to better predict the carbon budget in permafrost regions under future climate.

Appendix A: Description of ELMv1 subsurface heat transfer and modifications

ELMv1-ECA approximates the subsurface heat transfer process with a one-dimensional heat diffusion equation:

$$\begin{array}{}\text{(A1)}& c{\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\lambda }{\displaystyle \frac{\partial T}{\partial z}}\right),\end{array}$$

where T is the soil temperature (K), c is the volumetric soil heat capacity (J m⁻³ K⁻¹), λ is soil thermal conductivity (W m⁻¹ K⁻¹), and z is the soil depth (m) of the ELMv1-ECA soil layers. The ELMv1-ECA soil column consists of 15 layers, with soil thickness increasing exponentially with depth. The bottom of soil column is down to 42 m, and the top 10 layers are hydrologically active with layer node depth as 0.0071, 0.0279, 0.0623, 0.1189, 0.2122, 0.3661, 0.6198, 1.0380, 1.7276, and 2.8646 m, respectively. The soil heat capacity and thermal conductivity is updated at each time step based on the fractions of soil matrix components, i.e., liquid water content, ice content, and soil solids. The impact of organic carbon on soil thermal and hydraulic properties was incorporated as a linear combination of the counterpart properties of mineral soil and organic matter (Lawrence and Slater, 2008). To solve the (Eq. A1), ELMv1-ECA employs the Crank–Nicolson method, resulting in a tridiagonal system equation (Oleson et al., 2013). We assume a zero-flux bottom boundary condition. The top boundary condition is estimated by solving the energy balance equation at the air and ground interface, with additionally an overlying five-layer snow model and a one-layer surface water model in between. When snow and surface water are present, ELMv1-ECA incorporates the snow layers and surface water layer into the tridiagonal system to solve the heat transfer along the entire column.

ELMv1-ECA incorporates freeze–thaw processes of soil water in a decoupled way. Specifically, the model determines the onset of thawing or freezing by soil temperature initially solved at time step n+1 without consideration of the phase-change process, denoted as $T_{\mathrm{i}}^{n+\mathrm{1}}$, i.e.,

$$\begin{array}{}\text{(A2)}& \begin{array}{ll}{T}_{\mathrm{i}}^{n+\mathrm{1}}>T_{\mathrm{f}}\mathrm{and}{w}_{\mathrm{ice},\mathrm{i}}^n>\mathrm{0},& \mathrm{thawing},\\ T_{\mathrm{i}}^{n+\mathrm{1}}<T_{\mathrm{f}}\mathrm{and}{w}_{\mathrm{liq},\mathrm{i}}^n>w_{\mathrm{liq},\max ,\mathrm{i}}^{n+\mathrm{1}},& \mathrm{freezing},\end{array}\end{array}$$

where T_f is the freezing temperature of water (0 ^∘C in kelvin, i.e., 273.15 K), $w_{\mathrm{ice},\mathrm{i}}^n$ and $w_{\mathrm{liq},\mathrm{i}}^n$ are the mass of ice and liquid water (kg m⁻²) of layer i, and $w_{\mathrm{liq},\max ,\mathrm{i}}^{n+\mathrm{1}}$ (kg m⁻²) is the supercooled liquid water that is allowed to coexist with ice given the subfreezing soil temperature $T_{\mathrm{i}}^{n+\mathrm{1}}$. This $w_{\mathrm{liq},\max ,\mathrm{i}}^{n+\mathrm{1}}$ varies with soil texture and temperature and is calculated by the freezing-point depression equation (Niu and Yang, 2006),

$$\begin{array}{}\text{(A3)}& {\displaystyle }{w}_{\mathrm{liq},\max ,\mathrm{i}}^{n+\mathrm{1}}=\mathrm{\Delta }{z}_{\mathrm{i}}{\mathit{\theta }}_{\mathrm{sat},\mathrm{i}}{\left[{\displaystyle \frac{{\mathrm{10}}^{\mathrm{3}}{L}_{\mathrm{f}}\left(T_{\mathrm{f}}-T_{\mathrm{i}}^{n+\mathrm{1}}\right)}{g{T}_{\mathrm{i}}^{n+\mathrm{1}}{\mathit{\psi }}_{\mathrm{sat},\mathrm{i}}}}\right]}^{-\mathrm{1}/B_{\mathrm{i}}},\end{array}$$

where Δz_i is the soil thickness of the ith layer (in mm), θ_sat, i represents the soil porosity (i.e., the saturated volumetric water content), L_f is the latent heat of fusion (J kg⁻¹), B_i is the Clapp and Hornberger exponent (Clapp and Hornberger, 1978), g is the gravitational acceleration (m s⁻²), and ψ_sat, i is the soil-texture-dependent saturated matric potential (mm).

The rate of phase change is initially assessed from the heat excess (or deficit) needed to change the estimated temperature to the freezing point. Specifically, the model first computes the energy (H_i) needed for adjusting current soil temperature ($T_{\mathrm{i}}^{n+\mathrm{1}}$) to T_f:

$$\begin{array}{}\text{(A4)}& \begin{array}{ll}{H}_{\mathrm{i}}=-c_{\mathrm{i}}\frac{\mathrm{\Delta }{z}_{\mathrm{i}}}{\mathrm{\Delta }t}{T}_{\mathrm{inc}}+\left(\mathrm{1}-f_{\mathrm{sno}}-f_{\mathrm{h}\mathrm{2}\mathrm{osfc}}\right)\frac{\partial h}{\partial T}{T}_{\mathrm{inc}}& i=\mathrm{1}\\ H_{\mathrm{i}}=-c_{\mathrm{i}}\frac{\mathrm{\Delta }{z}_{\mathrm{i}}}{\mathrm{\Delta }t}{T}_{\mathrm{inc}}& i>\mathrm{1},\end{array}\end{array}$$

where $T_{\mathrm{inc}}=T_{\mathrm{f}}-T_{\mathrm{i}}^{n+\mathrm{1}}$, h is ground heat flux, and f_sno and f_h2osfc are the snow and surface water fraction within the grid cell, respectively. The mass change involved is then computed as $H_{\mathrm{m}}=\frac{H_{\mathrm{i}}\mathrm{\Delta }t}{L_{\mathrm{f}}}$ (i.e., $-c_{\mathrm{i}}\frac{\mathrm{\Delta }{z}_{\mathrm{i}}}{L_{\mathrm{f}}}{T}_{\mathrm{inc}}$ for soils below the top interface layer). That is, the mass of ice increased/decreased by freezing/melting is −H_m, releasing/absorbing energy H_i to bring up/down the current soil temperature to T_f. Accordingly, the ice and liquid mass are adjusted as

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& w_{\mathrm{ice},\mathrm{i}}^{n+\mathrm{1}}\\ & =\left\{\begin{array}{l}min\left(w_{\mathrm{ice},\mathrm{i}}^n+w_{\mathrm{liq},\mathrm{i}}^n-w_{\mathrm{liq},\max ,\mathrm{i}}^{n+\mathrm{1}*},w_{\mathrm{ice},\mathrm{i}}^n-H_{\mathrm{m}}\right)\\ \mathrm{when}{w}_{\mathrm{liq},\mathrm{i}}^n+w_{\mathrm{ice},\mathrm{i}}^n\ge w_{\mathrm{liq},\max ,\mathrm{i}}^{n+\mathrm{1}*}\\ \mathrm{0}\mathrm{when}{w}_{\mathrm{liq},\mathrm{i}}^n+w_{\mathrm{ice},\mathrm{i}}^n<w_{\mathrm{liq},\max ,\mathrm{i}}^{n+\mathrm{1}*}\end{array}\right.\\ & w_{\mathrm{liq},\mathrm{i}}^{n+\mathrm{1}}=\max (w_{\mathrm{liq},\mathrm{i}}^n+w_{\mathrm{ice},\mathrm{i}}^n-w_{\mathrm{ice},\mathrm{i}}^{n+\mathrm{1}},\mathrm{0}).\end{array}\end{array}$$

The H_i then is adjusted to $H_{i^{*}}$, calculated as $H_{i^{*}}=H_{\mathrm{i}}-\frac{L_{\mathrm{f}}\left(w_{\mathrm{ice},\mathrm{i}}^n-w_{\mathrm{ice},\mathrm{i}}^{n+\mathrm{1}}\right)}{\mathrm{\Delta }t}$. The $H_{i^{*}}$ then is the ultimately determined latent heat and is used to further readjust soil temperature as in equation (Eq. A6),

$$\begin{array}{}\text{(A6)}& T_{\mathrm{i}}^{n+\mathrm{1}*}=T_{\mathrm{f}}+{\displaystyle \frac{\mathrm{\Delta }t}{c_{\mathrm{i}}\mathrm{\Delta }{z}_{\mathrm{i}}}}{H}_{i^{*}}=T_{\mathrm{f}}-{\displaystyle \frac{L_{\mathrm{f}}\left(w_{\mathrm{ice},\mathrm{i}}^n-w_{\mathrm{ice},\mathrm{i}}^{n+\mathrm{1}}\right)}{c_{\mathrm{i}}\mathrm{\Delta }{z}_{\mathrm{i}}}},\end{array}$$

in which the temperature adjusted to T_f is further increased by $-\frac{L_{\mathrm{f}}\left(w_{\mathrm{ice},\mathrm{i}}^n-w_{\mathrm{ice},\mathrm{i}}^{n+\mathrm{1}}\right)}{c_{\mathrm{i}}\mathrm{\Delta }{z}_{\mathrm{i}}}$ due to soil freezing since $w_{\mathrm{ice},\mathrm{i}}^{n+\mathrm{1}}\ge w_{\mathrm{ice},\mathrm{i}}^n$ or decreased due to melting when $w_{\mathrm{ice},\mathrm{i}}^{n+\mathrm{1}}<w_{\mathrm{ice},\mathrm{i}}^n$.

To improve this scheme, we can incorporate soil water freezing phase change into equation (Eq. A1) and rewrite the heat transfer equation as (Eq. A7) or (Eq. A8),

$$\begin{array}{}\text{(A7)}& {\displaystyle }c{\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\lambda }{\displaystyle \frac{\partial T}{\partial z}}\right)-L_{\mathrm{f}}{\mathit{\rho }}_{\mathrm{liq}}{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{liq}}}{\partial t}},\text{(A8)}& {\displaystyle }\left(c+L_{\mathrm{f}}{\mathit{\rho }}_{\mathrm{liq}}{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{liq}}}{\partial T}}\right){\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\lambda }{\displaystyle \frac{\partial T}{\partial z}}\right),\end{array}$$

where L_f is the latent heat of fusion (J kg⁻¹), θ_liq is soil liquid water content (m³ m⁻³), and ρ_liq is the density of liquid water (kg m⁻³). To solve (Eq. A8), we need to compute the derivative of the soil freezing characteristic curve (θ_liq(T)) with respect to temperature ($\frac{\partial {\mathit{\theta }}_{\mathrm{liq}}}{\partial T}$). As discussed above, we approximate the θ_liq(T) curve by combining the freezing-point temperature-depression equation (Fuchs et al., 1978) and the soil water retention curve (Clapp and Hornberger, 1978). This leads to the supercooled water formulation (Eq. A3) (Niu and Yang, 2006). Computing $\frac{\partial {\mathit{\theta }}_{\mathrm{liq}}}{\partial T}$ requires the soil freezing curves θ_liq(T) to be continuous and differentiable for a range of temperatures during the freezing process (Kurylyk and Watanabe, 2013; Hansson et al., 2004). Here, we follow the existing ELM framework discussed above to solve Eq. (A8). The original numerical representation for readjusting soil temperature (Eq. A6) obtained by the uncoupled two-step implementation (Eqs. A2 to A5) significantly overestimates soil water freezing rates. There are two reasons for the overestimation. First, the freezing point (T_f=0 ^∘C) is used to determine the occurrence of soil water phase change under all conditions. To further freeze supercooled soil liquid water, however, the soil temperature has to be colder than a virtual soil temperature (as we described below). Second, due to the steep slope of $\frac{\partial {\mathit{\theta }}_{\mathrm{liq}}}{\partial T}$ (especially close to T_f=0 ^∘C), the estimated ice mass increase (i.e., $w_{\mathrm{ice}}^{n+\mathrm{1}}-w_{\mathrm{ice}}^n$ or $w_{\mathrm{liq}}^n-w_{\mathrm{liq},\max }^{n+\mathrm{1}}$; see, Eq. A5) most often exceeds the required mass change, i.e., $H_{\mathrm{m}}=-c_{\mathrm{i}}\frac{\mathrm{\Delta }{z}_{\mathrm{i}}}{L_{\mathrm{f}}}\left(T_{\mathrm{f}}-T_{\mathrm{i}}^{n+\mathrm{1}}\right)$, and thus soil liquid water freezes quickly in a large chunk. Soon, the liquid water available to be frozen becomes too small to release sufficient latent heat to compensate for the required energy deficit ($T_{\mathrm{f}}-T_{\mathrm{i}}^{n+\mathrm{1}}$).

Thus, we revised the phase-change scheme mainly through incorporating a phase-change efficiency (ε) and replacing the constant freezing point T_f with the temperature of the freezing-point depression in Eq. (A2). The phase-change efficiency, introduced by Le Moigne et al. (2012) and adopted by Masson et al. (2013) and Yang et al. (2018a), is calculated as

$$\begin{array}{}\text{(A9)}& {\displaystyle }\mathit{\epsilon }=\left\{\begin{array}{ll}{\mathit{\epsilon }}_{\mathrm{liq},\mathrm{i}}^n=\frac{{\mathit{\theta }}_{\mathrm{liq},\mathrm{i}}^n}{{\mathit{\theta }}_{\mathrm{sat},\mathrm{i}}}& \mathrm{for}\mathrm{freezing}\\ {\mathit{\epsilon }}_{\mathrm{ice},\mathrm{i}}^n=\frac{{\mathit{\theta }}_{\mathrm{ice},\mathrm{i}}^n}{{\mathit{\theta }}_{\mathrm{sat},\mathrm{i}}}& \mathrm{for}\mathrm{melting}\end{array}\right.,\end{array}$$

where ${\mathit{\theta }}_{\mathrm{liq},\mathrm{i}}^n$ and ${\mathit{\theta }}_{\mathrm{ice},\mathrm{i}}^n$ is the soil liquid and ice volumetric water content of layer i at previous time step n, respectively, and θ_sat, i represents the soil porosity (i.e., the saturated volumetric water content). The temperature of the freezing-point depression, as a virtual temperature (Tv) reversely derived from the freezing-point temperature-depression equation, i.e., $\mathit{\psi }\left(T\right)=\frac{{\mathrm{10}}^{\mathrm{3}}{L}_{\mathrm{f}}\left(T_{\mathrm{f}}-T_{\mathrm{i}}\right)}{g{T}_{\mathrm{i}}}$ (Fuchs et al., 1978; Cary and Mayland, 1972), is calculated as

$$\begin{array}{}\text{(A10)}& T{v}_{\mathrm{i}}^{n+\mathrm{1}}={\displaystyle \frac{{\mathrm{10}}^{\mathrm{3}}{L}_{\mathrm{f}}{T}_{\mathrm{f}}}{{\mathrm{10}}^{\mathrm{3}}{L}_{\mathrm{f}}+g{\mathit{\psi }}_{\mathrm{i}}^n}},\end{array}$$

where L_f is the latent heat of fusion (J kg⁻¹) and g is the gravitational acceleration (m s⁻²). ${\mathit{\psi }}_{\mathrm{i}}^n$ is the soil water potential (mm), calculated as the soil water retention curve of Clapp and Hornberger (1978), i.e., ${\mathit{\psi }}_{\mathrm{i}}^n={\mathit{\psi }}_{\mathrm{sat},\mathrm{i}}{\left(\frac{{\mathit{\theta }}_{\mathrm{liq},\mathrm{i}}^n}{{\mathit{\theta }}_{\mathrm{sat},\mathrm{i}}}\right)}^{-B_{\mathrm{i}}}$, where ${\mathit{\theta }}_{\mathrm{liq},\mathrm{i}}^n=w_{\mathrm{liq},\mathrm{i}}^n/\mathrm{\Delta }{z}_{\mathrm{i}}$ as in Eq. (A3), B_i is the Clapp and Hornberger exponent, and ψ_sat, i is the soil-texture-dependent saturated matric potential (mm).

Then, through multiplying the initially estimated mass change (H_m) by the phase-change efficiency (ε), we replace the freezing point with an efficiency-weighted average of the initially solved soil temperature ($T_{\mathrm{i}}^{n+\mathrm{1}})$ and the freezing point. Then, the updated soil temperature is expressed as

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}& T_{\mathrm{i}}^{n+\mathrm{1}*}=T_{\mathrm{f}}+{\displaystyle \frac{\mathrm{\Delta }t}{c_{\mathrm{i}}\mathrm{\Delta }{z}_{\mathrm{i}}}}\left(-c_{\mathrm{i}}{\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{i}}}{\mathrm{\Delta }t}}{T}_{\mathrm{inc}}{\mathit{\epsilon }}_{\mathrm{i}}+{\displaystyle \frac{L_{\mathrm{f}}\left(w_{\mathrm{ice},\mathrm{i}}^{n+\mathrm{1}}-w_{\mathrm{ice},\mathrm{i}}^n\right)}{\mathrm{\Delta }t}}\right)\\ & =T_{\mathrm{f}}-\left(T_{\mathrm{f}}-T_{\mathrm{i}}^{n+\mathrm{1}}\right){\mathit{\epsilon }}_{\mathrm{i}}+{\displaystyle \frac{L_{\mathrm{f}}\left(w_{\mathrm{ice},\mathrm{i}}^{n+\mathrm{1}}-w_{\mathrm{ice},\mathrm{i}}^n\right)}{c_{\mathrm{i}}\mathrm{\Delta }{z}_{\mathrm{i}}}}\\ & =\left(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{i}}\right)T_{\mathrm{f}}+{\mathit{\epsilon }}_{\mathrm{i}}{T}_{\mathrm{i}}^{n+\mathrm{1}}+{\displaystyle \frac{L_{\mathrm{f}}\left(w_{\mathrm{ice},\mathrm{i}}^{n+\mathrm{1}}-w_{\mathrm{ice},\mathrm{i}}^n\right)}{c_{\mathrm{i}}\mathrm{\Delta }{z}_{\mathrm{i}}}}.\end{array}\end{array}$$

Here, $w_{\mathrm{ice},\mathrm{i}}^{n+\mathrm{1}}$ is calculated by (Eq. A5) as well, but with updated H_m (i.e., $-{\mathit{\epsilon }}_{\mathrm{i}}{c}_{\mathrm{i}}\frac{\mathrm{\Delta }{z}_{\mathrm{i}}}{L_{\mathrm{f}}}(T{v}_{\mathrm{i}}^{n+\mathrm{1}}-T_{\mathrm{i}}^{n+\mathrm{1}})$). The two changes can effectively improve the soil water freezing process simulations and prevent soil becoming irreversibly too cold quickly as simulated by the baseline phase-change scheme.

Appendix B: Description of ELMv1 decomposition cascade model and modifications

ELMv1-ECA explicitly simulates carbon cycle dynamics (both plant and soil) and accounts for the limitation of nutrient (i.e., nitrogen and phosphorus) availability for plant growth and the nutrient competition between plants and microbes (Burrows et al., 2020; Zhu et al., 2019; Golaz et al., 2019; Zhu et al., 2020). The ELMv1-ECA uses a century-like soil carbon decomposition cascade model with vertically resolved soil biogeochemistry (Koven et al., 2013b) and explicitly accounts for the influence of substrate and nutrient availability on soil respiration (both root and microbes) (Zhu et al., 2019).

Within the ELMv1-ECA century decomposition cascade model, the respiration fractions are parameterized as the fraction of 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 (f_T) and moisture (f_W), an oxygen scalar (f_O), and a depth scalar (f_D), in a multiplicative way, i.e., f_total=f_Tf_Wf_Of_D.

We use a Q₁₀-based standard exponential function to account for the temperature effect on decomposition,

$$\begin{array}{}\text{(B1)}& f_T=Q_{\mathrm{10}}^{\left(\frac{T-T_{\mathrm{ref}}}{\mathrm{10}}\right)},\end{array}$$

where Q₁₀=1.5 by default, which is consistent with ecosystem-level observations (Mahecha et al., 2010), and T_ref 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 (f_W) within ELMV1-ECA is given in the formulation, calculated as

$$\begin{array}{}\text{(B2)}& f_{\mathrm{W}}=\left\{\begin{array}{ll}\mathrm{0}& \mathrm{For}{\mathit{\psi }}_{\mathrm{i}}<{\mathit{\psi }}_{\min }\\ \frac{\log \left({\mathit{\psi }}_{\min }/{\mathit{\psi }}_{\mathrm{i}}\right)}{\log \left({\mathit{\psi }}_{\min }/{\mathit{\psi }}_{\max }\right)}& \mathrm{For}{\mathit{\psi }}_{\min }\le {\mathit{\psi }}_{\mathrm{i}}\le {\mathit{\psi }}_{\max }\\ \mathrm{1}& \mathrm{For}{\mathit{\psi }}_{\mathrm{i}}>{\mathit{\psi }}_{\max }\end{array}\right.,\end{array}$$

where ${\mathit{\psi }}_{\mathrm{i}}={\mathit{\psi }}_{\max }{\left(\frac{{\mathit{\theta }}_{\mathrm{liq},\mathrm{i}}}{{\mathit{\theta }}_{\mathrm{sat},\mathrm{i}}}\right)}^{-B_{\mathrm{i}}}$ is the soil water potential, and B_i 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 equation, i.e., ${\mathit{\psi }}_{\mathrm{i}}=\frac{L_{\mathrm{f}}\left(T_{\mathrm{f}}-T_{\mathrm{i}}\right)}{{\mathrm{10}}^{\mathrm{3}}T}$ (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. By default, ψ_min is −10 MPa, which predicts zero f_W under frozen conditions since ψ_i under a subfreezing soil temperature easily gets smaller than −10 MPa (Fig. S1a). We thus reduced the ψ_min to −10³ and −10⁶ MPa to alleviate the zero respiration problem in the frozen soils (see Fig. S1a). We also tested scalars documented by Sierra et al. (2015) and Yan et al. (2018) and included the full list in Table S2. ELMv1-ECA approximates oxygen stress (f_O) as the ratio of available oxygen to that demanded by decomposers and has a minimum value of 0.2 (Oleson et al., 2013).

The depth scalar $\left(f_{\mathrm{D}}=\exp \left(-\frac{z_{\mathrm{i}}}{Z_{\mathit{\tau }}}\right)\right)$ represents other unresolved depth-dependent processes (i.e., soil microbial dynamics, priming effects, etc.) (Koven et al., 2013b; Lawrence et al., 2015; Koven et al., 2015). Applying the depth scalar to decomposition rates would exponentially decrease the respiration fluxes along with the vertical soil layers. The Z_τ is the e-folding depth for decomposition, and by default Z_τ is 0.5 m (Oleson et al., 2013).

Appendix C: Description of ELMv1 methane model and modifications

The ELMv1-ECA methane model includes the representations of CH₄ production, oxidation, and three pathways of transport (i.e., aerenchyma tissues, ebullition, and aqueous and gaseous diffusion) and solves the transient reaction diffusion equation for CH₄. ELMv1-ECA estimates CH₄ production (P; mol m⁻³ s⁻¹) in the anaerobic portion of the soil column based on the upland heterotrophic respiration (HR; mol C m⁻² s⁻¹) from soil and litter, further adjusted by factors representing influence from soil temperature (f_T), pH (f_pH), redox potential (f_pE), and seasonal inundation condition (S) (Riley et al., 2011), expressed as

$$\begin{array}{}\text{(C1)}& P=\mathrm{HR}\times f_{{\mathrm{CH}}_{\mathrm{4}}}\times f_T\times f_{\mathrm{pH}}\times f_{\mathrm{pE}}\times S.\end{array}$$

The $f_{{\mathrm{CH}}_{\mathrm{4}}}$ is a fraction of anaerobically mineralized carbon atoms becoming CH₄ and is 0.2 by default. Detailed explanation on these factors can be found in Riley et al. (2011). The methane production P is directly related to the estimated HR and impacted by soil temperature, and thus the changes in the carbon decomposition model (Appendix B) and water and heat transfer model (Appendix A) directly influence methane production simulations. Besides, ELMv1-ECA considers the availability of carbon substrate as an important driver of methanogenesis activity and methane production (Riley et al., 2011; Xu et al., 2016).

The ultimately estimated CH₄ emissions are controlled by oxidation, transport mechanisms (i.e., aerenchyma transport, ebullition, and diffusion), and the upper boundary resistance. Detailed descriptions on CH₄ oxidation and transport mechanisms are provided in Riley et al. (2011). Here we modified CH₄ transport mechanisms for facilitating reasonable cold-season CH₄ emissions.

Vascular plant aerenchyma tissues serve as diffusive pathways to transport CH₄ from the soil to the atmosphere. The CH₄ transport via aerenchyma from soil layer z (A(z), mol m⁻² s⁻¹) is calculated as

$$\begin{array}{}\text{(C2)}& A\left(z\right)=\left(C\left(z\right)-C_{\mathrm{a}}\right)/\left({\displaystyle \frac{r_{\mathrm{L}}z}{{\mathrm{DpT}}_{\mathrm{aere}}{\mathit{\rho }}_{\mathrm{r}}\left(z\right)}}+r_{\mathrm{a}}\right),\end{array}$$

where C(z) and C_a are the gaseous CH₄ concentrations (mol m⁻³) in soil depth z and in the atmosphere, respectively; r_a is the aerodynamic resistance (s m⁻¹); D is the gas diffusion coefficient (m² s⁻¹); p is aerenchyma porosity (–); r_L is the ratio of root length to vertical depth (i.e., root obliquity); and ρ_r(z) is the root fraction in soil depth z (–). T_aere is the specific aerenchyma area (m² m⁻²) and is expressed as

$$\begin{array}{}\text{(C3)}& T_{\mathrm{aere}}={\displaystyle \frac{f_{\mathrm{N}}{N}_{\mathrm{a}}\mathrm{LAI}}{\mathrm{0.22}}}\mathit{\pi }{R}^{\mathrm{2}},\end{array}$$

where R represents the aerenchyma radius ($=\mathrm{2.9}\times {\mathrm{10}}^{-\mathrm{3}}$ m); N_a is the annual net primary production (NPP), and f_N is the belowground fraction of current NPP; and the factor 0.22 represents average observed tiller biomass (gC per tiller) (Wania et al., 2010; Schimel, 1995). Here, tillers mean segmented stems of plants in the order of Poales, including grasses (Poaceae) and sedges (Cyperaceae) (Wania et al., 2010).

In ELMv1-ECA, methane emissions through aerenchyma were turned off when the soil temperature is below 0 ^∘C. We first removed this temperature limitation; then, we integrated the emissions from ice cracks and remnants of aerenchyma tissues with Eq. (C3) by applying a small ε_aere during winter time. That is, $T_{\mathrm{aere}}=\frac{f_{\mathrm{N}}{N}_{\mathrm{a}}{\mathit{\epsilon }}_{\mathrm{aere}}}{\mathrm{0.22}}\mathit{\pi }{R}^{\mathrm{2}}$ when soil temperature is below the freezing point, where ε_aere represents possible ice crack fractions and remnants of aerenchyma tissues, and it is 0.01 by default. We also tested a larger value for this parameter together with $f_{{\mathrm{CH}}_{\mathrm{4}}}$ and the maximum CH₄ oxidation rate constant R_o, max (Table S3).

The R_o, max is a key variable controlling CH₄ oxidation rate (R_oxic), calculated as

$$\begin{array}{}\text{(C4)}& R_{\mathrm{oxic}}=R_{\mathrm{o},\max }\left[{\displaystyle \frac{C_{{\mathrm{CH}}_{\mathrm{4}}}}{K_{{\mathrm{CH}}_{\mathrm{4}}}+C_{{\mathrm{CH}}_{\mathrm{4}}}}}\right]\left[{\displaystyle \frac{C_{{\mathrm{O}}_{\mathrm{2}}}}{K_{{\mathrm{O}}_{\mathrm{2}}}+C_{{\mathrm{O}}_{\mathrm{2}}}}}\right]Q_{\mathrm{10}}{F}_v,\end{array}$$

where $C_{{\mathrm{CH}}_{\mathrm{4}}}$ and $C_{{\mathrm{O}}_{\mathrm{2}}}$ are CH₄ and O₂ concentrations, respectively; $K_{{\mathrm{CH}}_{\mathrm{4}}}$ and $K_{{\mathrm{O}}_{\mathrm{2}}}$ are the half-saturation coefficients (mol m⁻³) for $C_{{\mathrm{CH}}_{\mathrm{4}}}$ and $C_{{\mathrm{O}}_{\mathrm{2}}}$, respectively. Details about the CH₄ oxidation are provided in Riley et al. (2011). The maximum oxidation rate R_o,max (mol m⁻³ s⁻¹) by default is $\mathrm{1.25}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{1.25}\times {\mathrm{10}}^{-\mathrm{6}}$ for saturated and unsaturated conditions, respectively. We tested the R_o, max with the smaller value ($\mathrm{1.25}\times {\mathrm{10}}^{-\mathrm{6}}$) for saturated condition as well (Table S3).

Another key variable that is highly uncertain is snow resistance to gas emissions. When snow is present, the upper boundary layer resistance to gas emissions is added by a snow resistance accounting for diffusion through the snow based on the Millington–Quirk expression (Riley et al. 2011). Specifically, the gaseous and aqueous diffusivity in snow is calculated by (Eq. C5)

$$\begin{array}{}\text{(C5)}& D_{\mathrm{snow}}=\left\{\begin{array}{l}{\mathit{\epsilon }}_{\mathrm{diff}}{D}_{\mathrm{gas}}\frac{{\left({\mathit{\theta }}_{\mathrm{air}}/\left({\mathit{\theta }}_{\mathrm{air}}+{\mathit{\theta }}_{\mathrm{water}}\right)\right)}^{\mathrm{10}/\mathrm{3}}}{{\left({\mathit{\theta }}_{\mathrm{air}}+{\mathit{\theta }}_{\mathrm{water}}\right)}^{\mathrm{2}}}\mathrm{for}\mathrm{gaseous}\\ {\mathit{\epsilon }}_{\mathrm{diff}}{D}_{\mathrm{liq}}{\left({\mathit{\theta }}_{\mathrm{water}}\right)}^{\mathrm{2}}\mathrm{for}\mathrm{aqueous}\end{array}\right.,\end{array}$$

where D_gas and D_liq are gaseous and aqueous diffusion coefficients (m² s⁻¹) for CH₄, respectively (see Table 2 in Riley et al. 2011), ε_diff is a scale factor (1 by default) for the diffusion coefficients, θ_water is water-filled fraction of snow volume (–), and θ_air is air fraction in snow (–). The D_snow is calculated for each snow layer; then, the total snow resistance is estimated as $R_{\mathrm{snow}}={\sum }_{i=\mathrm{1}}^n\frac{\mathrm{d}z\left(i\right)}{D_{\mathrm{snow}}\left(i\right)}$, where dz(i) is the thickness of ith snow layer. The top boundary layer conductance for CH₄ then is calculated as $\frac{\mathrm{1}}{R_{\mathrm{snow}}+R_{\mathrm{pond}}+\mathrm{1}/C_{\mathrm{g}}}$, where R_pond is ponding resistance and C_g is ground conductance (details can be found in Oleson et al., 2013). We found the computed snow resistance generally was too large. We thus decreased snow resistance by introducing a new scale factor, i.e., ε_snowdiff (10² by default), and replaced ε_diff in Eq. (C5) with ε_snowdiff to increase the conductance at the upper boundary when snow is present.