Evaluating the impact of peat soils and snow schemes on simulated active layer thickness at pan-Arctic permafrost sites

Permafrost stability is significantly influenced by the thermal buffering effects of snow and active-layer peat soils. In the warm season, peat soils act as a barrier to downward heat transfer mainly due to their low thermal conductivity. In the cold season, the snowpack serves as a thermal insulator, retarding the release of heat from the soil to the atmosphere. Currently, many global land models overestimate permafrost soil temperature and active layer thickness (ALT), partially due to inaccurate representations of soil organic matter (SOM) density profiles and snow thermal insulation. In this study, we evaluated the impacts of SOM and snow schemes on ALT simulations at pan-Arctic permafrost sites using the Energy Exascale Earth System Model (E3SM) land model (ELM). We conducted simulations at the Circumpolar Active Layer Monitoring (CALM) sites across the pan-Arctic domain. We improved ELM-simulated site-level ALT using a knowledge-based hierarchical optimization procedure and examined the effects of precipitation-phase partitioning methods (PPMs), snow compaction schemes, and snow thermal conductivity schemes on simulated snow depth, soil temperature, ALT, and CO2 fluxes. Results showed that the optimized ELM significantly improved agreement with observed ALT (e.g. RMSE decreased from 0.83 m to 0.15 m). Our sensitivity analysis revealed that snow-related schemes significantly impact simulated snow thermal insulation levels, soil temperature, and ALT. For example, one of the commonly used snow thermal conductivity schemes (quadratic Sturm or SturmQua) generally produced warmer soil temperatures and larger ALT compared to the other two tested schemes. The SturmQua scheme also amplified the model’s sensitivity to PPMs and predicted deeper ALTs than the other two snow schemes under both current and future climates. The study highlights the importance of accurately representing snow-related processes and peat soils in land models to enhance permafrost dynamics simulations.


Introduction
Peat soils and snow play crucial roles in the thermal insulation of permafrost and significantly impact the permafrost thermal dynamics (Decharme et al 2016).Peat soils, characterized by high soil organic matter (SOM) and soil water content, are commonly found in permafrost regions and are known to regulate permafrost thermal states and carbon dynamics (Schuur et al 2015).Hugelius et al (2014) reported that peat deposits extend below 3 m depth in many permafrost regions, highlighting the need to understand the intricate relationship between SOM and permafrost stability.Owing to the low thermal conductivity of high SOM content, peat soils act as a thermal buffer against downward heat transfer during the warm season, thus partially insulating permafrost soil from warming (Perreault andShur 2016, Du et al 2023).Many global land models, including DOE's Energy Exascale Earth System Model (E3SM) land model (ELM), prescribe SOM using observationally-based total soil carbon content datasets (e.g.NCSCD; Hugelius et al 2013b) at coarse spatial resolutions of ∼50 km.Also, these land models usually approximate the vertical SOM density profile along soil depth using one single distribution profile over the vast boreal region (Lawrence and Slater 2008) without considering the stratification of peatland soils that is significantly different from location to location (Hugelius et al 2013a, 2014, Mishra et al 2021).This oversimplification can lead to overestimated summer soil temperatures and active layer thickness (ALT, i.e. the maximum thaw depth overlying the permafrost) (Hugelius et al 2020).Given the distinct hydrological and thermal properties of organic matter, incorporating the stratified peatland soil structure and appropriately estimated vertical organic matter distribution is crucial for capturing the effects of peatland organic matter on permafrost warming (Jafarov and Schaefer 2016, Tao et al 2017, Tran et al 2017).Currently, many global land models overestimate summer soil temperatures and thus ALT in permafrost regions, partially attributed to the lack of representation of spatial variability of peatland organic matter density and stratification (Fisher et al 2016, Tao et al 2017, Loranty et al 2018, Hugelius et al 2020).At the site scale where measurements are available for validating model results, model accuracy is often poor, primarily due to the significant discrepancy between global surface datasets and in-situ conditions, including both the total SOM content and the highly uncertain vertical profile of soil carbon (Tao et al 2017).
Snow, acting as a thermal insulator to the ground during the cold season and impeding energy release to the atmosphere, plays a comparably important role as air temperature in affecting permafrost soil warming (Tao et al 2019, Mekonnen et al 2021).That is, a thicker snowpack generally provides better ground insulation, leading to warmer winter soil temperatures and more soil moisture once the snowpack thaws and snowmelt water infiltrates into soils.Inaccurately simulated snow thermal insulation can result in biased soil temperature simulations (Lawrence and Slater 2010, Jafarov et al 2014, Tao et al 2017, 2019), leading to large biases in simulated cold-season CO 2 and CH 4 emissions (Tao et al 2021a(Tao et al , 2021b) ) and soil carbon stocks (Koven et al 2009, Gouttevin et al 2012).Moreover, the biased winter soil temperature significantly influences microbial decomposition and nitrogen mineralization throughout the non-growing season and impacts microbial and plant community shifts (Schimel et al 2004, Riley et al 2018, 2021), thereby affecting the estimation of the pan-Arctic ecosystem carbon budget (Xu and Zhuang 2023).
Furthermore, snow dynamics also influence soil moisture and, thereby, plant communities (Liu et al 2018, Qi et al 2020).For instance, Liu et al (2018) found that snow depth could alter the fungi-tobacteria and the grass-to-forb ratios and shift the correlation relationship between microbial and plant biomass from positive to negative by alleviating soil water stress in temperate steppes.Moreover, soil moisture increases resulting from infiltrated snowmelt water impact soil thermal conductivity due to the larger thermal conductivity of soil water than air in the soil pore spaces, thus impacting subsurface heat transfer and the active layer (Tao et al 2017, Clayton et al 2021).Soil water content also affects latent heat of fusion during freeze-thaw cycles, thereby influencing soil temperature changes during shoulder seasons and subsequently impacting the active layer (Clayton et al 2021, Tao et al 2021a).On the other hand, the high water retention and resulting ice content of peat soils in the wintertime insulates permafrost due to the high heat capacity associated with high ice content, and thus slower soil temperature decreases during the cold season (Du et al 2023).Finally, high SOM-associated hydraulic conductivity may facilitate increased drainage fluxes and advective heat transport, particularly over hillslope landscapes, also impacting the ALT in regions where advective heat transport dominates permafrost thawing (Gao and Coon 2022).
To simulate snow thermal insulation correctly, a series of processes need to be well represented in the land models of Earth System Models (ESMs), including climate forcing, precipitation partitioning (i.e.rainfall and snowfall fraction before reaching the land surface), and snow physics (i.e.snow accumulation, compaction, snowmelt, and water and heat transfer throughout the snowpack) (Domine et al 2019).Current land models usually use a spatially and temporally static air temperature threshold to partition total precipitation into rain and snow, including ELM.Only a few hydrology and land surface models have tested coupled relationships between air temperature and relative humidity to partition the total precipitation (e.g.Behrangi et al 2018, Wang et al 2019).Uncertainties in these empirical static threshold-based precipitation-phase partitioning methods (PPMs) usually lead to large uncertainties in simulated snow depth and Snow Water Equivalent (SWE) (Harder and Pomeroy 2014, Wayand et al 2016, Jennings and Molotch 2019, Sun et al 2019).
In addition to climate forcing and the PPM employed, snow densification and compaction are critical processes that control snowpack evolution and snow density profile (Wang et al 2013).A recent study has shown that a more realistic representation of perennial snow densification leads to better simulations of glaciers and ice sheets (van Kampenhout et al 2017).Yet the impact of how snow compaction impacts permafrost remains unknown because snow insulation is not only controlled by accumulated snow mass but also by the small but highly variable thermal conductivity of the snowpack.Empirical equations between snow thermal conductivity (STC) and snow density are commonly used to parameterize snow thermal insulation (Sturm andJohnson 1992, Calonne et al 2011).The STC parameterizations currently used in ESMs (e.g. Yen 1965a, Jordan 1991, Sturm et al 1997) are usually derived from limited in-situ or laboratory observations, and thus need a close examination if used over permafrost regions.Different STC parameterizations produce substantial differences in simulated soil temperatures and, in regions with permafrost, active layer dynamics (Dutch et al 2022).
Here, we aim to evaluate and enhance the ELMsimulated permafrost ALT at observational sites distributed over the pan-Arctic region by (1) optimizing the peat SOM density profiles and (2) examining the impacts of PPMs, snow compaction mechanisms (including destructive metamorphism of new snow and snow densification due to overburden pressure), and STC schemes on ELM simulated snow depth, soil temperature, and ALT.We also discuss the feasibility of optimizing a global land model at a regional scale using spatially distributed observational sites, knowing existing model deficiencies.

CALM sites and datasets
We used ALT measurements from the Circumpolar Active Layer Monitoring (CALM) (www2.gwu.edu/∼calm/; Brown et al 2000) to evaluate model simulations.Figure 1 shows all the CALM sites across the pan-Arctic.We first consolidated CALM data for 219 individual locations out of 239 CALM sites by combining overlapping site locations.These observations were then averaged within ELM's 0.5-degree resolution grid cells by aggregating multiple sites within the same grid cells, resulting in a dataset of 160 CALM observed grid cells that could be directly compared with ELM simulation results.Subsequently, we refined our dataset by focusing solely on the Tundra and Taiga regions above 60 • N, thereby excluding sites in topographically complex regions (e.g. the Tibetan Plateau) and/or sites that have large measurement uncertainties (e.g. the sites in Mongolia), reducing the number of CALM observed grid cells to 108.Finally, to ensure data consistency, we retained only those grid cells that have observations available since 1990, resulting in a final dataset of 100 grid cells.We then conducted a series of ELM simulations at these 100 consolidated grid cells (section 2.2).
We used a snow climate classification map (Sturm et al 2009) and a permafrost map (Brown et al 2002) to categorize the types of snow-climate classes and permafrost regions at the CALM sites.In addition, we used the Snow Telemetry (SNOTEL) site measurements (https://wcc.sc.egov.usda.gov/)over Alaska to evaluate simulated snow depth, SWE, and soil temperature.Among the useful SNOTEL sites that have SWE or snow depth observations, we identified five sites that collocate with the assembled CALM sites, including three Tundra and two Taiga sites (figure 1).

Experiment design and optimization procedure
We used the E3SM land model (ELMv1-ECA) (Riley et al 2018, Golaz et al 2019, Zhu et al 2019, 2020, Tao et al 2021a) to run ensemble simulations at the assembled 100 grid cells that cover CALM sites.The ELM snow module simulates essential snow physical processes, including snow aging, accumulation, compaction, melting, sublimation, dew and frost on the top of the snow, and water and energy exchange at the snow surface (Oleson et al 2013, Golaz et al 2019).It also represents mass and heat transfer through the snowpack, which is discretized in up to five dynamic snow layers depending on snow depth.Although ELM does not explicitly represent factors such as wind slab and depth hoar layers that cause heterogeneity in the snow structure (Zhang et al 1996), ELM calculates the ice and water content and snow grain radius for each snow layer.Thus, the snow column has a vertically changing snow density profile.ELM estimates snow coverage fraction based on previous coverage and new snowfall, modified by the subsequent depletion processes (Swenson and Lawrence 2012).The ELM estimates heat transfer along the snow and soil column by solving a one-dimensional heat diffusion equation with the Crank-Nicholson method (Tao et al 2021a).The soil column comprises 15 layers with a bottom depth of ∼42 m, and the thickness for each soil layer increases exponentially with depth.The ELM calculates soil thaw depth at each time step, with the annual maximum thaw depth representing the simulated ALT.The simulation indicates permafrost-free if the simulated coldest soil temperature throughout the soil column remains consistently above 0 • C year-round or the predicted ALT extends deeper than or equal to the bottom depth of the ELM 10th soil layer (∼4 m).More model details can be found in Tao et al (2021a).
We used two reanalysis climate forcing datasets, including the Climatic Research Unit and Japanese Reanalysis (CRUJRA; Harris 2019) and Global Soil Wetness Project Phase 3 (GSWP3; Kim 2017), which ended in 2014 and then was extended to 2017 by the bias-corrected GSWP3 with the W5 × 10 5 global meteorological forcing data processed for Inter-Sectoral Impact Model Intercomparison Project, i.e.GSWP3-W5 × 10 5 (Lange et al 2019)).Simulations followed the spinup (i.e.cycling 20 years from 1901 to 1920 until the system reached equilibrium states) and transient run (i.e. from 1901 to 2017) protocols described The influence of SOM content on soil thermal and hydrologic properties (e.g.soil thermal conductivity, specific heat capacity, porosity, etc) is represented within current ELM (following Lawrence and Slater 2008) using the Northern Circumpolar Soil Carbon Database (NCSCD) (Hugelius et al 2013b) for high latitudes at a 0.5 • × 0.5 • spatial resolution, which significantly differ from the site-level SOM content at CALM sites.Also, the ELM distributes the NCSCD total soil carbon content across the soil layers according to a soil carbon vertical distribution profile for polar and boreal regions reported by Zinke et al (1986), which concentrates more carbon towards the surface than the tropical and temperate profile (Lawrence and Slater 2008).However, this vertical distribution profile of soil carbon content, although typical for polar and boreal regions, becomes highly uncertain regarding site-level simulations (Hugelius et al 2020).Therefore, to better simulate site-level permafrost dynamics using global land models such as ELM, it is critical to improve the representation of total soil carbon stocks and also the vertical distribution of peat soils.Similarly, to better simulate snow thermal insulation at the site level using global land models, processes that impact snow morphology, accumulation, ripening, melting, and the thermal properties of the snowpack need to be better examined.Here, we designed a knowledge-based hierarchical optimization procedure (figure 2), aiming to sequentially evaluate model components within ELM relevant to soil and snow dynamics.
First, to validate and examine the sensitivity of model simulated snow depth and SWE to climate forcing and partitioned precipitation (i.e.rain and snow fractions), we ran ELM with 18 PPMs (table 1 and  figure 3) based on a comprehensive review of rain-orsnow thresholds and ranges as examined by Jennings and Molotch (2019).Table 1 lists all the PPMs tested in this study, and a list of all the experiments is included in supplementary table S1.Details about these PPMs can be found in Jennings and Molotch (2019).With different PPMs, the snowfall fraction relationship with air temperature varies significantly, especially during the shoulder seasons (e.g.freezing and thawing seasons) when the atmospheric temperature is around the range (e.g.−2 • C to 5 • C) that the PPMs are most sensitive to (figure 3).
Then, we tested the impact of snow compaction parameterizations within ELM since the snow compaction mechanisms further impact snow depth and density.We tested a scheme with an increased snow compaction rate by doubling the fractional compaction rate and the upper limit on destructive metamorphism and decreasing the viscosity coefficient for the overburden pressure (see '1. Snow compaction schemes' in the supplementary file).Also, we tested three STC schemes, which usually increase with snow Further, soil thermal properties, largely influenced by SOM density, are crucial in regulating subsurface heat transfer through diffusion.To address the SOM representation at multiple sites distributed across the pan-Arctic, we tested ten peat stratification structures by making ELM soil layers purely peat soils from the top layer down to the 10th layer (illustrated in the supplementary table S1).Specifically, we made the topsoil purely organic matter layer in the 1st experiment by assigning its organic matter density as the bulk density of peat (i.e.130 kg m −3 ) (Farouki 1981), made the top two layers as peat soils in the 2nd experiment, and so on for all ten layers for the 10th experiment.Note this modification still does not fully address distinct hydro-thermal processes for peat soils, but it better accounts for the thermal conductivity and hydraulic properties than the baseline simulation that uses SOM density derived from the coarse global datasets.An additional ten simulations were conducted with a modified plant functional type (PFT) dataset that assigned the bare ground friction (again extracted from a 0.5 • spatial resolution global dataset) to Arctic grass to better representing site-level vegetation conditions.Among these simulations, we searched for the minimum ALT bias for each site and used the SOM density profile and modified PFT.
Lastly, variations in soil temperature regimes influence plant growth, which can, in turn, influence snow accumulation processes through various mechanisms (e.g.altering precipitation interception and the surface albedo).Specifically, Arctic grass, broadleaf deciduous boreal shrub, needleleaf evergreen boreal trees, and deciduous boreal trees are the dominant plant types over the tundra and taiga areas.The Arctic grass and broadleaf deciduous boreal shrub are found in tundra sites, whereas taiga sites usually have a combination of the four dominant plant types.These plant types, possessing different leaf areas and canopy sizes that cause differences in precipitation and light interception, significantly impact soil temperature and thereby nitrogen mineralization, affecting the feedback strength between plant growth and snow accumulation (as discussed in section 3.3).
All the simulations (table S1) were conducted following the spinup and transient simulations for all the modified soil and snow schemes (table S1).Finally, we identified the best simulation by comparing simulated ALT climatology to CALM observations.This approach not only provides an optimized modeling system and simulation results that agree well with observations, thereby enhancing the modeling capability to predict future permafrost dynamics confidently, but also enables a comprehensive assessment of snow impacts on ELM-simulated ALTs.

Evaluation of ELM simulated ALT at CALM sites
We identified optimal combinations of SOM density profiles and PFTs for each site that provided the best agreement with CALM observations, i.e. resulting in a minimum RMSE in ALT (table S2).The simulated ALTs demonstrated significant improvement compared to baseline results at the 100 CALM sites (figure 5(a)).Specifically, the optimized ELM achieved a significant improvement in simulated ALT with a RMSE of 0.15 m, which was reduced from 0.83 m in baseline ELM.We also found that the baseline simulation overestimated summer soil temperatures and thus predicted permafrost-free conditions for many CALM sites in discontinuous or sporadic permafrost regions.In contrast, the optimized simulation improved the number of permafrost sites by 21 (figure 5(b)).Note that this optimized SOM density might not be realistic because the current ELM lacks representation of a moss layer and a seasonally dynamic surface litter layer.(Although ELM represents litter as an important biogeochemical component, the model doesn't account for litter's hydrologic and thermal effects.Without this surface litter and moss layer, the optimization procedure mimicked the litter and moss thermal buffer effects by concentrating more organic matter in topsoils and thus might overestimate SOM in the top ∼50 cm.Results thus revealed that for about half the sites, the model needs to include peat soils down to ∼50 cm to facilitate sufficiently slow heat transfer during warm seasons to simulate better summer soil temperature and thus ALT (figures S1 and S2).However, incorporating deeper peat soils would increase porosity in the deep, thick soils and facilitate slower temperature changes due to the large heat capacity of ice (associated with the large porosity) during cold seasons.Therefore, further incorporation of peat soils into deep soil layers, in contrast to incorporating peat soils in the top layers, would slightly increase rather than decrease ALTs at some sites (figure S1).
Although the optimized SOM density resulted in different profiles from site to site, the identified simulations all assign bare ground to the default Arctic grass PFT fraction.Also, instead of searching for the best snow scheme for each site, we identified the one snow scheme set that provided a minimum RMSE in site-averaged ALT climatology.Results demonstrated that the default snow scheme, combined with the optimized surface datasets, still showed superior performance than the other tested snow schemes (table S2).Next, we examine the sensitivity of ELMsimulated snow depth, SWE, soil temperature at different layers, and the net ecosystem exchange (NEE, i.e.net exchange of CO 2 with the atmosphere) to the employed snow schemes at the five SNOTEL sites that collocate with CALM sites.We also validate the simulated snow depth, SWE, and soil temperatures at the sites where measurements are available.Next, we examine ALT sensitivity to PPMs and STC schemes at all the assembled CALM sites.To evaluate snow influence on ALT simulations solely, we conducted comparisons with simulations using different snow schemes but maintaining the optimized surface datasets (e.g.SOM density and PFTs).

Snow depth and SWE
Generally, simulations with CRUJRA forcing underestimate snow depth significantly at all five SNOTEL sites, regardless of PPMs and snow compaction schemes (figure 6).In contrast, GSWP3 forcing provided better-simulated multi-year averaged seasonal cycles of snow depth with improved magnitudes of snow depth peaks, although simulations show earlier peak timing at three sites.The shaded areas in figure 6, indicating uncertainty associated with PPMs in simulated snow depth, are generally large during the fall-to-winter transition season (e.g.September-October) when the snow starts to accumulate, followed by the thawing season (e.g.May-June).This uncertainty is also large around the snow depth peak months (e.g.February to April) at two sites (961 and 1177), probably owing to the climate forcing characteristics (e.g.Tair, RH, etc).When using the scheme with increased snow compaction, simulations greatly underestimated snow depths at all the sites, although they did not show significant differences in snow mass (i.e.SWE), as shown by the SWE comparison at site 1174 (figure 6).Consequently, the increased snow compaction scheme resulted in a much larger snow density, facilitating a large STC, given the empirical relationships between snow density and STC (figure 4).

Soil temperature and ALT
The increased snow compaction scheme reduced snow thermal insulation due to increased snow density and thus STC (figure 7), causing faster heat release from soil to the atmosphere during cold seasons, and leading to greatly underestimated winter soil temperature (ELM_GSWP3_Jordan&IncComp).As a result, even though the ELM_GSWP3_Jordan&IncComp simulated SWE agrees with observations very well (figure 6), its simulated soil temperatures are overly cold (second only to the coldest simulations with CRUJRA), owing to its underestimated snow depth and overestimated snow density.Although winter soil temperatures are less important than warm-season climate on warm-season soil temperatures and ALT, they affect the onset of the thawing season and plant leaf-out, thereby significantly impacting plant carbon uptake and ecosystem carbon dynamics (as discussed in section 3.2.3).In addition, with the same forcing and the default PPM0, simulations that used the lowest STC scheme (ELM_GSWP3_Yen) and the quadratic Sturm equation (ELM_GSWP3_SturmQua) generally show warmer soil temperatures than the simulation using the default Jordan scheme (ELM_GSWP3_Jordan) during cold seasons (figure 7).Constrained by the correctly simulated snow depth, here ELM_GSWP3_Jordan still shows better simulations of soil temperature, e.g. with smaller RMSE of 3.4 • C for the 5 cm soil temperature at site 1177, compared to 3.9 • C and 4.7 • C for ELM_GSWP3_Yen and ELM_GSWP3_SturmQua, respectively.
On average, ELM_GSWP3_SturmQua simulated multi-year averaged annual soil temperature is much warmer than the other two simulations, e.g. is 1.9 • C and 0.7 • C warmer than ELM_GSWP3_Jordan and ELM_GSWP3_Yen, respectively, at site 1177 (figure 7).As a result, ELM_GSWP3_SturmQua generally overestimated ALTs, whereas ELM_GSWP3_Jordan shows good agreement with CALM observed ALTs (figure 8).In addition, ELM_GSWP3_Yen, being similar to the default Jordan scheme for small snow density scenarios (figure 4), also demonstrates a similar performance in simulating ALT to ELM_GSWP3_Jordan at some sites but overestimates ALTs at many other sites (figure 8(b)).For the five SNOTEL sites where we examined the accuracy of simulated snow depth, the site-averaged mean bias in ALT is 0.10 m, 0.24 m, and 0.46 m for ELM_GSWP3_Jordan, ELM_GSWP3_Yen, and ELM_GSWP3_SturmQua, respectively, indicating better simulations with the default Jordan STC scheme.

Ecosystem carbon fluxes
Snow thermal insulation influences cold-season microbial activities (e.g.decomposition and nitrogen mineralization) and affects plant growth via accumulating nutrients for the plant to use in the following year (Riley et al 2018(Riley et al , 2021)).As figure 9 shows, ELMsimulated NPP with different STC schemes demonstrates large differences in warm seasons.These NPP differences are not only related to varying soil temperature owing to different STC schemes but are also associated with cold-season nitrogen mineralization resulting from slow but persistent microbial processes, which is clearly shown by cold-season (September to May) differences in NEE (i.e.microbial respiration during the cold season) (figure 9).As a result, simulations using the SturmQua STC predicted more CO 2 emissions (i.e. more positive NEE) than other simulations due to the larger microbial respiration persistently lasting throughout the cold season.

Predicted ALT changes under RCP8.5
With the default PPM0 and the Jordan STC scheme, ELM predicted that 54% of the currently simulated permafrost sites would become permafrostfree by the end of the 21st century.This number is 78% when using the SturmQua scheme.Figure 10 shows the time series of predicted ALT ensembles with PPMs and the Jordan and SturmQua schemes under RCP 8.5.The ensemble mean of predicted ALTs with the Jordan (SturmQua) scheme estimated an increasing rate of ALT of 1.0 ± 0.1 cm year −1 (0.4 ± 0.2 cm year −1 ) and 3.8 ± 0.3 cm year −1 (4.1 ± 1.5 cm year −1 ) for tundra and taiga sites, respectively.The SturmQua scheme generally predicted deeper ALTs than the Jordan scheme.Notably, the SturmQua scheme predicted that Taiga sites would transition to a permafrost-free state by 2041, which is 36 years earlier than the prediction made by the Jordan scheme (figure 10(b)).Nevertheless, due to the limited number of Taiga sites, especially across the vast extent of the Siberian region, the predicted ALT trends exhibit significant uncertainty and may not accurately represent certain Taiga areas that currently lack CALM sites.
Additionally, plant growth indirectly impacts snow physical processes by affecting surface albedo and precipitation interception by the canopy, thus facilitating a negative feedback loop between snow thermal insulation and plant growth, i.e. more snow thermal insulation leads to warmer soils and higher NPP, which, in turn, reduces snow depth and snow thermal insulation.This negative feedback between snow thermal insulation and plant growth gradually gained prominence after ∼2050 for the SturmQua simulations (figure 10(a)) at the Tundra sites, and thus, the SturmQua scheme projected ALT declined to a similar magnitude as predicted by the Jordan scheme at the end of the 21st century.

Summary and discussion
This study evaluated the impact of SOM and snow schemes on ALT simulations at permafrost sites across the pan-Arctic.Specifically, by testing site-level peat SOM vertical distribution structure and examining snow schemes (e.g.PPMs, snow compaction, and STC schemes), we enhanced ELM simulations of permafrost ALT at 100 grid cells CALM sites via a knowledge-based hierarchical optimization strategy (i.e.reduced ALT RMSE from 0.83 m for the baseline to 0.15 m).We found that, because the model does not adequately represent the thermal buffer layers, the optimized simulations prescribed peat soils up to ∼50 cm depth at 11 sites, around ∼80 cm at 18 sites, and ∼140 cm at 21 sites (figure S2) to mimic the thermal insulation from the surface litter and moss layer during warm seasons.These thick peat layers in permafrost regions are not uncommon, as highlighted by Hugelius et al (2014), who found peat deposits could extend to depths of up to 3 m.Nevertheless, the optimized SOM content profile and snow schemes may still have biases, even though ELM with the optimized schemes produced improved simulations of site-level ALT.The purpose of the developed method here is to accurately represent ALT, given the current model structure.We are actively working on incorporating the moss and surface litter layers into ELM to mechanistically account for their thermal insulation effects.However, before a more sophisticated model becomes available (which usually introduces new parameters that need to be determined at the global scale), the optimization scheme used in this study is a reasonable step toward better estimation and prediction of permafrost thermal dynamics.
We found that modeled soil temperatures are very sensitive to snow-related schemes, including the PPMs, snow compaction scheme, and the STC scheme.We also noticed a saturation of snow thermal insulation capacity (i.e.modeled snow thermal insulation does not increase with snow depth beyond an effective snow depth of about 50 cm) (Slater et al 2017).Thus, simulated ALT results demonstrated a larger sensitivity to STC than to simulated snow depth.For the five sites where we examined the accuracy of simulated snow depth, the Jordan scheme simulated site-averaged mean ALT showed the smallest bias.Through the sensitivity analysis, we demonstrated that the default ELM snow scheme provides reasonable results at pan-Arctic CALM sites with the modified SOM density profiles and PFTs.However, large uncertainties are found in the ensemble simulations with different snow schemes, and in general, the SturmQua amplified model sensitivity to PPMs and predicted deeper ALTs than the other two snow schemes under both current and future climates.
This study elucidates a roadmap to optimizing and advancing land modeling within ESMs across permafrost regions, with an emphasis on the critical importance of reasonably representing snow-related processes and peat soil representations.Our findings underscore the importance of a nuanced examination of coupled hydrology, biogeochemistry, and plant modeling when simulating and predicting permafrost dynamics.This examination includes close scrutiny of climate forcing, precipitation-phase partitioning, snow accumulation processes, STC schemes, energy and water exchange at the land-atmosphere interface, subsurface heat and water transfer, and the plant-snow feedback mechanism.Specifically, a calibration or optimization procedure solely for J Tao et al the subsurface heat transfer process is not reliable unless the snow physical properties and thermal processes are correctly simulated.This comprehensive scrutiny is essential for ensuring both accuracy and an understanding of uncertainty propagation within land models, thereby improving the understanding of the complex interactions in pan-Arctic permafrost ecosystems.Nevertheless, without mechanistic representations of all thermal buffer layers in the model and the availability of comprehensive measurements to support the knowledge-based hierarchical optimization strategy, potential biases and uncertainties may still persist in the optimized schemes.
In addition to the snow schemes we examined, the snow fractional cover, which impacts surface albedo, absorbed solar radiation, outgoing longwave radiation, and thus net radiation, also plays a role in affecting permafrost thermal states (Zhang et al 1996, Zhang 2005) and ecosystem carbon cycling.The future warming climate might cause precipitation to fall more as rainfall than snowfall (Berghuijs et al 2014, Bintanja and Andry 2017, Hyncica and Huth 2019), resulting in decreased snow depth and snow fractional cover in some areas.Given the projected faster warming of the cold season compared to the warm season (Tao et al 2021b), future studies are needed to examine the net impact of regionalscale snow conditions on permafrost thermal states and ecosystem biogeochemical cycling.

Figure 1 .
Figure 1.Study domain and sites.The background snow-climate classification map is based on (Sturm et al 2009), and the permafrost outline map is adapted from (Brown et al 2002) and Diamond symbols represent CALM sites across the pan-Arctic domain (above the northern latitude of 60 • .The circles represent the SNOTEL sites that have measurements relevant to this study and are paired with CALM sites.

Figure 2 .
Figure 2. Iterative knowledge-based hierarchical optimization strategy employed by this study.Black arrows indicate how the forcing data and the physical process parameterizations (blue boxes) directly impact the variables (yellow boxes).Green dashed lines represent the indirect influence of certain variables on other variables, i.e. the feedback loops.A more detailed description of this optimization procedure can be found in the supplementary file.

J
Schemes Snowfall fraction of total precipitation

Figure 3 .
Figure 3. Four types of PPMs were tested in this study (table 1).(a) Air temperature (Tair)-based interpolation method (ELM uses PPM0 by default).(b) Threshold-based deterministic method with Tair, dew point temperature (Tdew), or wet bulb temperature (Twet).(c) Binary Logistic Regression probability method with two variables, i.e.RegBi-a function of Tair and RH.(d) Binary Logistic Regression probability method with three variables, i.e.RegTri-a function of Tair, RH, and Ps.The RH values are fixed at 80% and 60% for solid and dashed lines, respectively.

Figure 4 .
Figure 4. We tested three snow thermal conductivity (STC) schemes commonly used by global land models.By default, ELM uses the STC scheme by Jordan (1991).

Figure 5 .
Figure 5. (a) Comparison between CALM measured and ELM simulated ALT with baseline and optimized surface datasets.Sites the ELM baseline simulations erroneously predict permafrost-free conditions, in contrast to the optimized simulations that better represent permafrost with valid ALT values, are excluded from the scatter comparison.(b) The number of sites with valid ALT values observed by the CALM network is compared to the number of permafrost sites simulated by both the baseline and optimized ELM models, showing a significant improvement in capturing site-level permafrost active layer.

Figure 6 .
Figure 6.Comparison of observed and simulated snow depth at five SNOTEL Alaska sites collocated with CALM sites.Among the five sites, only site 1174 has SWE measurements; thus, SWE comparison is only shown at this site.Shaded areas indicate uncertainty associated with precipitation-phase partitioning methods (PPMs), which influence the simulations of snow depth and SWE together with climate forcing (e.g.CRUJRA and GSWP3) and the snow compaction schemes used.The error bars indicate interannual variability for both simulations and observations.The simulations with increased snow compaction (SnowIncComp) reduced snow depth and increased density, leaving snow mass (i.e.SWE) not much changed.

Figure 7 .
Figure 7.Comparison between soil temperature time series observed at SNOTEL Alaska site 1177 and simulated by ELM with two forcings (CRUJRA and GSWP3) and different snow thermal schemes.Jordan, Yen, and SturmQua indicate the three thermal conductivity schemes (section 2.2), and IncComp means the snow scheme with an increased compaction rate.The shaded areas (for ELM_CRUJRA_Jordan, ELM_GSWP3_Jordan, and ELM_GSWP3_Jordan&IncComp) indicate uncertainty associated with PPMs, i.e. the minimum-to-maximum soil temperature envelopes simulated with the tested PPMs (table1).

Figure 8 .
Figure 8.(a) ALT climatology comparison between CALM measurements and ELM simulations with three STC schemes, i.e.Jordan (by default), Yen (1965b), and the Sturm quadratic equation (Sturm et al 1997), at the five SNOTEL sites that collocated with CALM sites.Error bars indicate standard deviations of simulated ALTs with different PPMs.(b) Scatter comparison between CALM observed ALT multi-year average and ELM simulation ensemble mean at the same CALM sites as in figure 5(a).Here, the vertical bars are standard deviations of simulation ensembles at each site, representing uncertainty associated with PPMs.

Figure 9 .
Figure 9. Sensitivity of ELM simulated NPP and NEE to climate forcing and STC schemes.NPP differences among the simulations are attributed to the influence of employed STC schemes on soil temperatures and microbial activities (e.g.decomposition and nitrogen mineralization), which can be reflected by cold-season NEE differences (September to May).

Figure 10 .
Figure 10.ELM predicted ALT ensemble mean under RCP8.5 climate for (a) Tundra sites and (b) Taiga sites (figure 1).Shaded areas indicate uncertainty associated with the PPMs.Time series end when half of the ensemble simulations predict permafrost-free.
et al 2015 Climate change and the permafrost carbon feedback Nature 520 171-9 Slater A G, Lawrence D M and Koven C D 2017 Process-level model evaluation: a snow and heat transfer metric Cryosphere 11 989-96 Sturm M, Holmgren J, König M and Morris K 1997 The thermal conductivity of seasonal snow J. Glaciol.43 26-41 Sturm M, Holmgren J and Liston G 2009.Global seasonal snow classification system.Version 1.0 (UCAR/NCAR-Earth Observing Laboratory) Sturm M and Johnson J B 1992 Thermal-conductivity measurements of depth hoar J. Geophys.Res.97 2129-39 Sun N, Yan H X, Wigmosta M S, Leung L R, Skaggs R and Hou Z S 2019 Regional snow parameters estimation for large-domain hydrological applications in the Western United States J. Geophys.Res.-Atmos.124 5296-313 Swenson S C and Lawrence D M 2012 A new fractional snow-covered area parameterization for the community land model and its effect on the surface energy balance J. Geophys.Res.-Atmos.117 D21107 Tao J, Koster R D, Reichle R H, Forman B A, Xue Y, Chen R H and Moghaddam M 2019 Permafrost variability over the Northern Hemisphere based on the MERRA-2 reanalysis Cryosphere 13 2087-110 Tao J, Reichle R H, Koster R D, Forman B A and Xue Y 2017 Evaluation and enhancement of permafrost modeling with the NASA catchment land surface model J. Adv.Model.Earth Syst.9 2771-95 Tao J, Zhu Q, Riley W J and Neumann R B 2021a Improved ELMv1-ECA simulations of zero-curtain periods and cold-season CH4 and CO2 emissions at Alaskan Arctic tundra sites Cryosphere 15 5281-307 Tao J, Zhu Q, Riley W J and Neumann R B 2021b Warm-season net CO2 uptake outweighs cold-season emissions over Alaskan North Slope tundra under current and RCP8.5 climate Environ.Res.Lett.16 055012 Tran A P, Dafflon B and Hubbard S S 2017 Coupled land surface-subsurface hydrogeophysical inverse modeling to estimate soil organic carbon content and explore associated hydrological and thermal dynamics in the Arctic tundra Cryosphere 11 2089-109 van Kampenhout L, Lenaerts J T M, Lipscomb W H, Sacks W J, Lawrence D M, Slater A G and van den Broeke M R 2017 Improving the representation of polar snow and firn in the community earth system model J. Adv.Model.Earth Syst.9 2583-600 Wang T, Ottlé C, Boone A, Ciais P, Brun E, Morin S, Krinner G, Piao S and Peng S 2013 Evaluation of an improved intermediate complexity snow scheme in the ORCHIDEE land surface model J. Geophys.Res.Atmos.118 6064-79 Wang Y H, Broxton P, Fang Y H, Behrangi A, Barlage M, Zeng X B and Niu G Y 2019 A wet-bulb temperature-based rain-snow partitioning scheme improves snowpack prediction over the drier western united states Geophys.Res.Lett.46 13825-35 Wayand N E, Stimberis J, Zagrodnik J P, Mass C F and Lundquist J D 2016 Improving simulations of precipitation phase and snowpack at a site subject to cold air intrusions: Snoqualmie Pass, WA J. Geophys.Res.-Atmos.121 9929-42 Xu Y M and Zhuang Q L 2023 The importance of interactions between snow, permafrost and vegetation dynamics in affecting terrestrial carbon balance in circumpolar regions Environ.Res.Lett.18 044007 Yen Y C 1965b Effective thermal conductivity and water vapor diffusivity of naturally compacted snow J. Geophys.Res.70 1821-5 Yen Y-C 1965a Heat transfer characteristics of naturally compacted snow Research Report 166 (U.S. Army Cold Regions Research and Engineering Laboratory) p 9 Zhang T J 2005 Influence of the seasonal snow cover on the ground thermal regime: an overview Rev. Geophys.43 RG4002 Zhang T, Osterkamp T E and Stamnes K 1996 Influence of the depth hoar layer of the seasonal snow cover on the ground thermal regime Water Resour.Res.32 2075-86 Zhu Q, Riley W J, Iversen C M and Kattge J 2020 Assessing impacts of plant stoichiometric traits on terrestrial ecosystem carbon accumulation using the E3SM land model J. Adv.Model.Earth Syst.12 e2019MS001841 Zhu Q, Riley W J, Tang J Y, Collier N, Hoffman F M, Yang X J and Bisht G 2019 Representing nitrogen, phosphorus, and carbon interactions in the E3SM land model: development and global benchmarking J. Adv.Model.Earth Syst.11 2238-58 Zinke P J, Stangenberger A G, Post W M, Emanuel W R and Olson J S 1986 Worldwide Organic Carbon and Nitrogen Data ONRL/CDIC-18 (Carbon Dioxide Information Centre)