The role of snow cover affecting boreal-arctic soil freeze – thaw and carbon dynamics

Northern Hemisphere permafrost affected land ar­ eas contain abont twice as mnch carbon as the global at­ mosphere. This vast carbon pool is vnlnerable to acceler­ ated losses throngh mobilization and decomposition nnder projected global warming. Satellite data records sparming the past 3 decades indicate widespread redactions ( ~ 0.81.3 days decade “ ̂ ) in the mean aimnal snow cover extent and frozen-season dnration across the pan-Arctic domain, coinci­ dent with regional climate warming trends. How the soil car­ bon pool responds to these changes will have a large impact on regional and global climate. Here, we developed a conpled terrestrial carbon and hydrology model framework with a de­ tailed 1-D soil heat transfer representation to investigate the sensitivity of soil organic carbon stocks and soil decompo­ sition to climate warming and changes in snow cover condi­ tions in the pan-Arctic region over the past 3 decades (19822010). Om resnlts indicate widespread soil active layer deep­ ening across the pan-Arctic, with a mean decadal trend of 6.6 ± 12.0 (SD) cm, corresponding to widespread warming. Warming promotes vegetation growth and soil heterotrophic respiration particnlarly within snrface soil layers (< 0.2 m). The model simnlations also show that seasonal snow cover has a large impact on soil temperatmes, whereby increases in snow cover promote deeper (> 0.5 m) soil layer warm­ ing and soil respiration, while inhibiting soil decomposition from snrface (< 0.2 m) soil layers, especially in colder cli­ mate zones (mean armnal T < -10°C ). Onr resnlts demon­ strate the important control of snow cover on northem soil freeze-thaw and soil carbon decomposition processes and the necessity of considering both warming and a change in precipitation and snow cover regimes in characterizing per­ mafrost soil carbon dynamics.


Introduction
The northem high latitndes contain abont twice as mnch car bon as the global atmosphere, largely stored in permafrost and seasonally thawed soil active layers (Hngelins et al., 2014). This vast carbon pool is vnlnerable to accelerated losses throngh mobilization and decomposition nnder re gional warming, with potentially large global carbon and climate impacts Schaefer et al., 2011;Schnnr et al., 2015). The northem high latitndes have expe rienced a mnch stronger warming rate than the global aver age over recent decades (Serreze and Francis, 2006), and this warming trend is projected to continne, along with a gen eral increase in snrface precipitation (Solomon et al., 2007). A better nnderstanding of how the northem soil carbon pool responds to these changes is critical to predict climate feed backs and associated impacts to northem ecosystems.
The potential vnlnerability of soil carbon to mobilization and accelerated decomposition with climate warming, par ticnlarly in permafrost areas, will largely depend on changes in soil moistme and thermal conditions (Grosse et al., 2011;Schaefer et al., 2011;Schnnr et al., 2015). Widespread soil thawing and permafrost degradation in the boreal and Arctic have been reported (e.g., Jorgenson et al., 2006; Pnblished by Copernicns Publications on behalf of the European Geoseienees Union. 5812 Y. Yi et al.: The role of snow eover affeeting boreal-aretie soil freeze-thaw et al., 2010a, b). This has triggered a series of changes in bo real and Arctic ecosystems, inclnding changes in lake and wetland areas (Smith et al., 2005;Watts et al., 2012), tnndra shmb cover expansion (Tape et al., 2006;Stnrm et al., 2005), thermokarst and other distnrbances (Grosse et al., 2011)which are likely having a profonnd inllnence on both snrface and snbsnrface hydrology -and biogeochemical cycles. In particnlar, increases in soil temperatnre and associated soil thawing potentially expose vast soil organic carbon stocks, formerly stabilized in perennial frozen soils, to mobilization and decomposition, which may promote large positive cli mate feedbacks (Schaefer et al., 2011;Schnnr et al., 2015).
Previons stndies have highlighted the importance of both snrface air temperatme and snow cover conditions affecting the soil thermal regime among many other factors (Stieglitz et al., 2003;Zhang, 2005;Osterkamp, 2007;Lawrence and Slater, 2010;Romanovsky et al., 2010a). Changes in the rate of accnmnlation, timing, dnration, density and amonnt of snow cover dnring the winter season play an important role in determining how soil responds to snrface warming dne to strong insnlation effects of snow cover on gronnd tempera tnre and its role in the snrface energy bndget (Zhang, 2005). Both snrface warming and a changing precipitation regime can modify seasonal snow cover conditions, leading to a non linear soil response to warming (Lawrence and Slater, 2010). Increases in winter precipitation and a deepening of the snowpack may enhance soil warming, while a rednced snowpack, dne to precipitation decreases or warming-enhanced snow snblimation, may promote soil cooling. Changes in snow cover dmation and condition can also alter the amonnt of energy absorbed by the gronnd and modify the rate of soil warming (Enskirchen et al., 2007). The Arctic is expected to experience continned warming and precipitation increases nnder projected climate trends (Solomon et al., 2007); how these climate trends will affect soil moistme and thermal dy namics is a key qnestion affecting potential changes in north em soil carbon dynamics and associated climate feedbacks.
Satellite data records over the past 3 decades (1979-2011) indicate widespread rednctions (~ 0.8-1.3 days decade"^) in mean annnal snow cover extent and frozen-season dnra tion across the pan-Arctic domain, coincident with regional warming (Brown and Robinson, 2011;Kim et al., 2012). An earlier onset of spring snowmelt and soil thaw has been ob served from both in sitn gronnd and satellite measnrements, while the onset of snow cover and soil freezing in the fall show more variable trends (Brown and Robinson, 2011;Kim et al., 2012). More active snowmelt dnring the snow sea son, largely in the early snow season, has also been observed from satellite observations of regional snow cover extent and snrface freeze-thaw cycles (Kim et al., 2015). On the other hand, snow depth trends in the boreal-Arctic region show large spatial variability. For example, several stndies have shown a general snow depth increase in eastem Siberia (e.g.. Park et al., 2014) and a decrease in westem North America in recent decades (Dyer and Mote, 2006).
The objective of this stndy is to assess how northem soil thermal and carbon dynamics respond to snrface warming and changes in snow cover conditions dnring the satellite era (since 1979). To that end, we developed a conpled hydrology and carbon model framework with detailed soil heat transfer representation adapted for the pan-Arctic basin and Alaska domain. We nsed this model to investigate recent climaterelated impacts on soil thermal and carbon dynamics over the past 3 decades . We condncted a sensitivity analysis by mnning the model with different configmations of snrface meteorology inpnts to evalnate how soil thermal conditions and soil carbon dynamics respond to changes in air temperatme and precipitation dnring the same period.

Model deseription
A conpled hydrology and carbon model was nsed to inves tigate the sensitivity of the soil thermal regime and soil car bon decomposition to changes in snrface air temperatme and snow cover conditions. The hydrology model acconnts for the effects of soil organic layers, changes in snrface snow cover properties and soil water phase change on the soil freeze-thaw process in permafrost landscapes (Rawlins et al., 2013). These factors represent important controls on soil thermal dynamics within the active layer (Nicolsky et al., 2007;Lawrence andSlater, 2008, 2010), enabling an im proved estimation of snbsnrface soil temperatnre and moistnre profiles, particnlmly in permafrost meas, and a represen tation of essential environmental constraints on soil cmbon decomposition.
The hydrology model nsed for this investigation is an ex tension of previons efforts regarding Imge-scale pan-Arctic water balance modeling (PWBM; Rawlins et al., 2003Rawlins et al., , 2013. Recent npdates to the model inclnde an improved simnlation of snow or gronnd and snbsnrface temperatnre dynamics nsing a 1-D heat transfer eqnation (Rawlins et al., 2013) instead of the empirical thaw depth estimation based on the Stefan solntions nsed in Rawlins et al. (2003). The npdated PWBM model has 23 soil layers down to 60 m below snrface, with increasing layer thickness at depth. Up to five snow layers are nsed to acconnt for the effects of seasonal snow cover evolntion on the gronnd thermal regime, and changes in seasonal snow density and thermal condnctivities are also considered. Other model improve ments inclnde acconnting for the impact of soil organic car bon content on soil thermal and hydranlic properties (Ap pendix Sect. A l, Eq. A3); this impact is an important featnre of boreal and Arctic soils (Lawrence and Slater, 2008). Fmther details on the npdated hydrology model are provided in Appendix Sect. AL A satellite-based terrestrial carbon llnx (TCF) model  was conpled to the hydrology model for this Y. Yi et al.: The role of snow eover affeeting boreal-aretie soil freeze-thaw 5813 investigation. The TCF model uses a light use efficiency al gorithm driven by satellite estimates of FPAR (fraction of vegetahon canopy intercepted photosynthehcally active radi ation) to eshmate vegetation productivity and litterfall inputs to a soil decomposihon model. In the original TCF model, soil carbon stocks and respirahon fluxes were eshmated us ing a simplifled three-pool soil orgarhc carbon (SOC) decom posihon framework with enviromnental constraints on soil decomposihon rates derived from either satellite-estimated surface soil moistme and temperature flelds (Kimball et al., 2009) or reanalysis data . This approach as sumes that the major source of soil heterotrophic respiration (/?h) comes from surface (< 10 cm) liher and surface orgaihc layers. However, the contribution of deeper soils to total may be non-negligible, especially in high-latitude boreal and Arche tundra landscapes with characteristic carbon-rich soils Schnnr et al., 2015). Therefore, in this study, we incorporated a more detailed soil decomposition model represenhng SOC stocks, extending to 3 m below the surface and representing differences in litterfall and soil or ganic maher substrate quality within the soil proflle (Thorn ton et al., 2002). The resulting soil decomposition model used for this study includes three litterfall pools, three SOC pools with relatively fast turnover rates and a deep SOC pool with a slow turnover rate (Fig. SI in the Supplement). The three litterfall pools were distributed within the top 20 cm of the soil layers; the three fast SOC pools were distributed within the top 50 cm of the soil layers, and the deep SOC pool extended from 50 cm to 3 m below the surface. Substanhal SOC may be stored in permafrost soils below 3 m depth (Hngelins et al., 2014) and may potentially undergo mobilization with conhnued warming. However, this contri bution to total land-atmosphere carbon (CO2) exchange was assumed negligible for the recent historical period examined (Schaefer et al., 2011) and was not considered in this study. Further details on the carbon model used in this study are provided in Appendix Sect. A2.

Data sets
The modeling domain for this investigahon encompasses the pan-Arche drainage basin and Alaska, represenhng a total land area extent of approximately 24.95 millionkm^. The model was run at a 25 km Northem Hemisphere Equal-Area Scalable Earth Grid (EASE-Grid) spatial resolution and daily hme step from 1979 to 2010. Further details on the model validation data sets and inputs used for this study are pro vided below.

In situ data
In situ measurements from 20 eddy covariance (EC) tower sites across the pan-Arctic domain were obtained from the La Thuile FLUXNET data set (Baldocchi, 2008) and were used to evaluate the model-simulated daily carbon fluxes and soil temperature and moistme flelds (Supplement Ta  ble SI). These tower sites represent major vegetation commrmity types across the study domain and have at least 1 year of observations available. For validahon, the model was driven using tower-observed meteorology. The tower daily carbon flux observations are derived from half-homly EC CO2 flux measurements that have been processed and aggre gated using consistent gap fllling and quality control proce dures (Baldocchi, 2008). Limited surface (< 15 cm) soil tem perature and moisture measmements were also provided at a portion of the tower sites but with unknown soil sampling depths and very few measurements at the tundra sites. There fore, we selected one boreal forest and one tundra site with detailed in situ measurements (including carbon fluxes, soil temperatme and soil moistme) for addihonal model evalua tion (Table 1). The boreal forest site represents a single tower, whereas the tundra site includes three towers, represenhng three different tundra commmhty types. The tundra site is located within the Imnavait Creek wa tershed in the northem foothills of the Brooks Range, Alaska (68°37'N, 149° 18'W), and underlain with continuous per mafrost (Enskirchen et al., 2012). Mean aimual air tempera ture and precipitahon at the site is -7 .4 °C and 318 mm with about 40 and 60 % of armnal precipitation occurring as rain and snow, respechvely. There me three towers in three dif ferent tundra commmhty types, including dry heath, moist acidic tussock and wet sedge tundra. The surface soil organic layer thickness varies from 34.0 ± 2.4 cm in wet sedge tundra to 2.3 ± 0.3 cm for dry heath tundra. The maximum achve layer thaw depth varies from ~ 40 cm at the dry heath site to ~ 70 cm at the tussock tundra site (Enskirchen et al., 2012). Soil temperature and moistme at 5 cm depth were measured within each tundra tower footprint. All observations includ ing cmbon fluxes and soil temperature and moisture me avail able from 2008 to 2011.
The boreal forest site used in this study is part of a net work of tower EC sites spamhng a lire chronosequence in cenhal Marhtoba (55°54'N, 98°31'W ) at various stages of succession following large stand replacement flres (Goulden et al., 2011). We chose one of the two oldest chronosequence tower sites burned in 1930 for model validation because this site had more continuous measmements of carbon fluxes and surface meteorology and high-quality data (indicated by the tower metadata) during the observation period (2002)(2003)(2004)(2005). This site is dominated by mature closed-canopy black sprace stands. The mean armnal air temperature and precipitation at this site are -3 .2 °C and 520 mm, respectively. Soil temper atures were measured at the surface (0 cm) and at multiple (6,11,16,18,29,41 and 55 cm) soil depths, while soil mois ture was also measmed at mulhple (11,18,28,41 and 55 cm) depths.

5814
Y. Yi et al.: The role of snow eover affeeting boreal-aretie soil freeze-thaw Table 1. Characteristics o f two selected tundra and boreal forest tower sites used for model validation. Three tundra types are represented by the tower measurements at Imnavait Creek, A laska, including dry heath, moist acidic tussock and wet sedge tundra. The boreal forest site encompasses a set o f tower eddy covariance (EC) sites and measurements spanning a regional fire chronosequence at various succession stages in central M anitoba, Canada.

Model inputs
Primary model drivers include daily surface meteorology and satellite-based normalized difference vegetation index (NDVI) records. Daily average and minimum air tempera tnre, precipitation, wind speed, atmosphere vapor pressme deficit (VPD) and downward solar radiation were obtained from a new version of the WATCH Forcing Data (WFD) ap plied to the ERA-Interim reanalysis (WFDEI; Weedon et at., 2014). Tliis data set was created by extracting and interpo lating the ERA-Interim reanalysis to 0.5° x 0.5° spatial res olution with sequential elevation correction of snrface mete orological variables and monthly bias correction from gridded observations inclnding CRU TS (Climatic Research Unit Time Series; v3.1 and v3.2) and GPCC (Global Precipita tion Climatology Centre; v5 and v6) data sets (for precipi tation only). The daily WFDEI snrface meteorology data is available from 1979 to 2010 and allows more thorough com parisons of hydrological model outputs with other relevant satellite products than the previons WFD data set (Weedon et at., 2014). The third-generation Global Inventory Mon itoring and Modelling Stndies (GlMMS3g) NDVI data set (Xn et at., 2013) was nsed to estimate litterfall seasonality and FPAR, as critical inpnts to the TCF model (Yi et at., 2013). The GlMMS3g data set was assembled from different NOAA advanced very high-resolntion radiometer (AVHRR) sensor records, acconnting for various deleterions effects in clnding calibration loss, orbital drift and volcanic eraptions. The NDVI data have a 15-day temporal repeat and 8 km spa tial resolution, extending from 1982 to 2010. For the model simnlations, both WFDEI and GlMMS3g forcing data sets were regridded to a consistent 25 km EASE-Grid format and the bimonthly GlMMS3g data was interpolated to a daily time step. The NDVI data from 1982 were nsed as drivers for model spin-np and simnlations prior to the start of the GlMMS3g observation record (i.e., 1979-1981).
Other ancillary model inpnts included a merged 8 km land cover data set  combining the 500 m MODIS International Geosphere-Biosphere Programme (IGBP) land cover map (Friedl et al., 2010) and the Circnmpolar Arctic Vegetation Map (CAVM; Walker et al., 2005). The CAVM was nsed to identify tnndra vegetation within the circnmpo lar region as a supplement to the IGBP classification, which does not provide a specific category for tnndra and foresttnndra transition biome types . The dominant land cover type within each 25 km EASE-Grid cell was cho sen based on the merged 8 km land cover data set and reclas sified according to the original PWBM land cover classifi cation (Rawlins et al., 2013;Fig. S2). Tnndra, forest-tnndra and taiga-boreal biomes acconnt for approximately 70 % of the total pan-Arctic drainage basin area (Fig. S2).
Soil organic carbon inventory data (GSDT, 2000;Hngelins et al., 2014) were nsed to prescribe the SOC fraction in each model soil layer. The fraction of SOC has a large im pact on soil thermal and hydranlic properties and is there fore an important control on characterizing soil freeze-thaw and moisture processes (Lawrence and Slater, 2008;Nicol sky et al., 2007). The IGBP Global Soil Data Task (GSDT, 2000) and the Northem Circnmpolar Soil organic Carbon Database (NCSCD; Hngelins et al., 2014) SOC data were distributed throngh the top 11 model soil layers (< 1.4 m depth) across the stndy area following Rawlins et al. (2013) and Lawrence and Slater (2008). The NCSCD data, which provide an npdated estimate of SOC in permafrost affected areas, were nsed to prescribe the SOC fraction for permafrost areas, while the GSDT data were applied to non-permafrost areas. Generally, the organic carbon fraction within the top 5 soil layers (< 23 cm depth) is high, with mean values of 53.7 and 39.4% for the two deeper snrface soil layers (13-23 cm depth) averaged over the pan-Arctic domain.

Model parameterization
A dynamic litterfall allocation scheme based on satellite NDVI data (Appendix Sect. A2) was nsed to prescribe the daily litterfall fraction throngh each annnal cycle to acconnt for litterfall seasonality, particnlarly for deciduous vegeta tion types (Randerson et al., 1996;White et al., 2000). The GlMMS3g NDVI bimonthly data were first aggregated to a monthly time step and then nsed to characterize monthly leaf loss and tnmover rates of fine roots dnring the active growth period based on Eq. (A7). The monthly litterfall frac tion was then evenly distribnted at a daily time step within each month. This approach generally allocates more litterfall dnring the latter half of the growing season, while the model simnlations show generally more soil heterotrophic respira tion dnring the latter portion of the year (Fig. S3). A compari son of model simnlations against tower measnrements shows an overall improved net ecosystem exchange (NEE) season ality relative to a previons TCF model application where lit terfall was distribnted evenly over the armnal cycle .

Model sensitivity analysis
We condncted a model sensitivity analysis to examine how the estimated soil thermal regime and SOC decomposition respond to changes in snrface air temperatme and snow con ditions over the most recent 3 decades. Three sets of daily model simnlations were mn by (1) varying air temperatme (T) and precipitation (P) inpnts; (2) varying T inpnts alone (temperatnre sensitivity analysis), and (3) varying P inpnts alone (precipitation sensitivity analysis). Daily mean T (in clnding daily mean and miitimnm temperatnre) and P clima tology was first derived from the iititial 3-yem (1979)(1980)(1981) WFDEI meteorological record and nsed in the model sensi tivity mns. The daily climatology, based on 3-yem (1979)(1980)(1981) meteorological records rather than a single yem (i.e., 1979), was nsed to minimize effects from characteristically large climate llnctnations in the northem high latitndes. For precipitation, we first created a monthly climatology from the daily record (1979)(1980)(1981) and then scaled the daily WFDEI precipitation by maintaining the monthly climatology valne (Lawrence and Slater, 2010): where y, m and d represent a particnlar year, month and day; P{m) is the precipitation monthly climatology averaged from 1979 to 1981 and P(y,m) is the monthly total pre cipitation for a particnlm yem and month; P{y,m,d) and P'{y,m,d) me the original and scaled daily precipitation, respectively, for a particnlar year, month and day. Dne to a relatively short record (i.e., 1979-1981) and Imge variability in northem latitnde precipitation, the ratio of too large for a particnlar month with very low precipitation rates. In this case, the daily precipitation was not adjnsted to avoid nmeasonable estimates. We then ran the model with different configmations of the daily snrface meteorology data sets. Model simnlations derived nsing the dynamic WFDEI daily snrface meteorology from 1979 to 2010 (i.e., varying T and P) were nsed as the model baseline simnlation. For the temperatme sensitivity analysis, we ran the model nsing the dynamic daily WFDEI temperatnre records from 1979 to 2010 bnt holding P as the climatology valne from 1979 to 1981. For the precipitation sensitivity analysis, we ran the model nsing the dynamic daily WFDEI precipitation records bnt with the T daily climatology. Since VPD is dependent on air temperatnre, we also created a daily VPD climatol ogy (1979)(1980)(1981) as an additional inpnt to the carbon model, assnming negligible changes in relative hnmidity dnring the stndy period for the precipitation sensitivity analysis. There was no sigitificant trend in solar radiation dnring the stndy period; we therefore nsed the historical (i.e., 1979-2010) so-Im radiation data for the three sets of simnlations. The model was iititialized nsing a two-step process prior to the three sets of simnlations. The model was first spnn-np nsing the daily snrface climatology (1979)(1980)(1981) inclnding T, VPD, and P for 50 yems to bring the top 3 m soil tem peratnre into dynamic eqnilibrinm; the model was then mn nsing the same climatology and simnlated soil temperatnre and moistme fields over several thonsand yems to bring the SOC pools to eqnilibrinm.
We mainly nsed correlation analysis to evalnate the cli matic controls on simnlated soil temperatnre and cmbon llnxes. The ontpnts from the model baseline simnlations (i.e., varying T and P) from 1982 to 2010 were nsed for this anal ysis. The period from 1979 to 1981 was exclnded in order to rednce the impact of the spin-np process on model simnla tions. We first calcnlated the correlation coefficients between the time series of each climate variable and modeled soil tem peratnre or carbon llnxes at each grid cell from 1982 to 2010. The resnlting correlation coefficients were then averaged for each climate zone classified nsing the armnal mean air tem peratnre  and biimed into 2.5 °C intervals. The climate variables nsed in the correlation analysis inclnded air temperatme, snow water eqnivalent (SWE) and snow cover extent (SCE). The model did not simnlate SCE directly, and the SCE was estimated nsing the following eqnation: where SNOWD is the simnlated snow depth (m), and the smface ronghness was set as 0.1 m (Lawrence and Slater, 2010).

Model validation
The model simnlations were generally consistent with ob served daily cmbon llnxes from the 20 EC tower sites across the pan-Arctic domain (Table 2), with mean R valnes of 0.84 ±0.11 (SD) for gross primary prodnctivity (GPP) and 0.63 ± 0.17 for NEE, and mean RMSE differ ences of 1.44 ± 0.50 g C m "^ d"^ for GPP and 1.04 ± 0 .3 6 g C m "^d "^ for NEE. The model resnlts showed relatively Imge discrepancies with the tower-based cmbon llnxes for tnndra sites; however, Imge nncertainties me associated with the tower measmements in tnndra meas dne to the charac teristically harsh environment and extensive missing data. The simnlated temperatnre and moistnre fields also captnre Table 2. Coefficient o f determination (R^) and root mean square error (lyvISE) differences between model-simulated daily carbon fluxes and in situ tower EC measurement-based observations across the study area. The mean o f tower-observed daily GPP flux is also shown.
The uncertainty o f the estimates including mean, and RM SE values was indicated as a standard deviation when there were multiple sites represented for each plant functional type.

PFT
Tower sites 12 2.18 ± 1 .2 3 0 .7 0 ± 0 .1 7 1.46 ± 0 .5 9 0.34 ± 0 .1 5 1.06 ± 0 .4 0 DBF 2 2.11 ± 0 .9 6 0.82 ± 0 .0 2 1.31 ± 0 .6 0 0.59 ± 0 .0 4 1.29 ± 0 .3 9 MXF 3 1.99 ± 1 .0 2 0.77 ± 0 .0 3 1.46 ± 0 .4 5 0 .5 8 ± 0 .1 1 1.00 ± 0 . the seasonality of the in sitn snrface (< 15 cm) soil measnre ments representing variable soil depths (not shown), despite large nncertainties in the snrface meteorology inpnts (partic nlarly precipitation or snowfall) and soil parameters, inclnd ing the definition of textme and peat fraction within the soil profile. Additional assessment of the model simnlations was condncted nsing detailed in sitn measnrements at selected tnndra and boreal forest validation sites (Table 1) as snmmarized below. The model simnlations compared favorably with in sitn measnrements at the tnndra validation sites for snrface soil temperatnre (R =0.93, RMSE = 3.12 °C) and carbon llnxes, inclnding GPP (R = 0 .7 2 , RMSE = 0.76 g C m "^d "^) and NEE {R = 0.79, RMSE = 0.50 g C m" 2 d" i) bnt had a rela tively larger discrepancy dnring the winter when the model showed lower valnes of NEE (e.g., less CO2 emissions) than the measnrements (December to Febmary, DJF; Fig. 1). The simnlated maximnm soil thaw depth (~ 50 cm averaged from 2008 to 2011) was also consistent with site measnrements, ranging from 40 to 70 cm at three locations within the tnndra validation site (Enskirchen et al., 2012). An apparent cold bias ranging from -2 to -5 °C in the simnlated soil tem peratme dnring the fall and winter period of 2009 and 2010 ( Fig. la) reflects lower model-simnlated snow depth and as sociated rednctions in thermal bnffering between the atmo sphere and nnderlying soil layers. This cold bias in the simn lated soil temperatnres resnlts in early freezing of simnlated soil water content (Fig. S4). Compared with the tower obser vations, the simnlated daily snrface soil temperatnres gener ally show large temporal variations, particnlarly dnring the snmmer (Jnne to Angnst, JJA). There were also considerable differences among in sitn soil temperatmes at the different tnndra sites. Snmmer (JJA) soil temperatme at the wet sedge tnndra location was generally lower than for the other tnndra vegetation types, which may reflect higher soil water con tent and specific heat capacity and greater latent heat loss from evapotranspiration, leading to slower soil warming at this site. Overall, the model simnlations compare well with the tower-observed cmbon llnxes dnring the growing season bnt significantly nnderestimate NEE and soil respiration dming the dormant season. Model nnderestimation of soil respi ration dnring the dormant season may reflect less liqnid soil water represented by the model nnder frozen (<0°C ) tem peratnres than the tower measnrements ( Fig. S4) as well as a lack of model representation of wind-indnced CO2 exchange between the atmosphere and snrface snowpack (Luers et al., 2014). The model generally shows emlier seasonal onset and offset of photosynthesis relative to the in sitn measmements, while partitioning of the tower NEE measmements dnring the shonlder season may be snbject to large nncertainties nn der partial snow cover conditions (Enskirchen et al., 2012).
The model simnlations also compmed favorably against observations at the boreal forest validation sites (Fig. 2), captnring observed seasonality in soil temperatnres (R>0.95, RM SE<2.00 °C) at different soil depths and daily variations in tower-observed carbon llnxes for GPP {R = 0.89, RMSE = 1.24 g C d " i ) and NEE {R = 0.73, RMSE = 0.65 g C m "^d "^). Similm to the tnndra sites, snow depth also has a Imge impact on simnlated soil temperatmes at the boreal forest sites bnt is snbject to Imge nncertainties from both model snowfall inpnts and forest canopy snow interception processes. The timing of simnlated thaw and freeze of soil water at different depths is generally consistent with the tower measnrements, with later seasonal thawing and freezing occnrring in deeper soils (Fig. S5). The tower site soil moistnre measmements show Imger variability than the model simnlations dnring the growing season and likely reflect differences in the model pmameterization of snrface moss or peat and mineral soil hydranlic condnctivities relative to local site conditions. The model-simnlated NEE llnxes dnring the non-growing season stem mainly from soil heterotrophic respiration and were averaged aeross three tundra eommunity types, ineluding dry heath, moist aeidie tussoek and wet sedge tundra exeept for the NEE measurements during the winter. N EE measurements were not eolleeted at the tussoek tundra site during the winter; therefore, the winter N EE measurements were averaged for the dry heath and wet sedge tundra sites only.
are largely consistent with the in sitn tower observations, generally diminishing towards the end of the year and then gradnally recovering with soil warming toward the onset of the growing season. Both the model and in sitn tower NEE flnxes show large temporal variations dnring the growing season, largely dne to GPP rednctions cansed by high vapor pressnre deficits or water stress. The model-simnlated SCE was generally consistent with satellite-observation-based global climate data records docnmenting weekly SCE changes (Brown and Robinson, 2011;Fig. 3). The model simnlations show a similar mean seasonal cycle as the satellite observations, with spring snowmelt mostly occnrring from April to May and fall onset of sea sonal snow cover occnrring in October over the 1982 to 2010 record (Fig. 3a). The model-simnlated SCE shows consis tent changes with the satellite observations in spring, indicat ing realistic simnlation of the snow melting process. How ever, the model generally nnderestimates SCE in the fall and winter. The model did not directly simnlate SCE, which was calcnlated from simnlated snow depth nsing an empir- ical eqnation (Eq. 2). Based on Eq. (2), the modeled SCE will never approach 100 %, while the satellite data indicates nearly complete winter snow cover over the stndy domain. Larger model SCE differences from the satellite observations are expected when the snow cover is relatively shallow and patchy owing to the relatively coarse spatial resolntion of both model simnlations and satellite observations. Moreover, the satellite SCE data set is presented as a binary classifica tion at a weekly time step, which may not adeqnately depict transient SCE llnctnations nnder active snrface melting and freezing processes in the fall (Kim et al., 2015).

Climatic control on simnlated permafrost and soil temperatnres
The simnlated permafrost area is generally consistent with reported estimates from previons stndies. The simnlated mean permafrost area from 1982 to 2010 is approximately 11.3 million km^, which is within the range of observationbased estimates (11.2-13.5 millionkm^) of the combined area for continnons (90-100%) and discontinnons (50-90%) permafrost extent over the northem polar region (> 45° N) (Zhang et al., 2000). The simnlated active layer depth (ALD) shows an overall increasing trend across the pan-Arctic domain over the 1982 to 2010 record (Fig. 4a, b). No strong bias was observed for the model ALD simulations compared to in situ obser vations for 53 pan-Arctic sites from the Circnmpolar Active Layer Monitoring (CALM) program (Brown et a l, 2000); these results showed a mean model bias of -9.48 cm, repre senting approximately 16.5 % of the estimated ALD but with low model correspondence (R =0.31 , p<0. 1) relative to in situ observations (Fig. S6). The discrepancy between modelsimulated ALD results and in situ observations may be partly due to a spatial scale mismatch between the coarse-resolution model simulations and the local CALM site measurements, as well as imcertainties in the reanalysis surface meteorology data used as model forcings (Rawlins et al, 2013). Previ ous studies have shown large local spatial variations in ALD due to strong surface heterogeneity including microtopogra phy, vegetation and soil moisture conditions (Romanovsky et al., 2010a, b;Mishra and Riley, 2014). Simulated widespread ALD deepening is consistent with generally decreasing snow cover extent in the pan-Arctic region (Fig. 4c). Simulated ALD trends over the 1982-2010 record range from -4.32 to 8.05 cm yr" \ with a mean value of 0.66cm yr"^ A no table model ALD deepening trend occurs in discontinuous permafrost areas with relatively large mean ALD values. However, in portions of Alaska, the model simulations indi cate slightly decreasing ALD trends across the study period ( Fig. 4b), despite a strong reduction in the local snow cover extent (Fig. 4c). This mainly reflects a large decrease in the simulated snowpack (Fig. 4d) due to a decreasing trend in WFDEI precipitation or snowfall data, resulting in less ther mal insulation of tmderlying soil, which may offset warming effects from decreasing snow cover extent. The regional differences in snow cover effects on modelsimulated ALD can be explained by different climatic con trols on warm-season (May to October) soil temperatures. The correlation analysis between climate variables and warm-season soil temperatures (Fig. 5) indicates that sur face warming has a dominant control on upper (<0.5m ) soil temperatures in all climate zones, and also on deeper (> 0.5 m) soil temperatures in warmer climate zones (mean annual Tan > -4 °C). A deep snowpack has a strong warm ing effect on simulated deeper (> 0.5m ) soil temperatures in colder climate zones (mean annual Tan < -4 °C) but with limited warming effects on surface soil temperatures across all pan-Arctic climate zones. Correspondingly, the effects of seasonal snow cover duration on model soil temperatures vary across different climate zones and soil depths. In colder climate areas, a longer snow cover duration has a relatively Strong warming effect on deeper (> 0.5 m) soil temperatnres bnt with negligible warming effects on snrface soil layers, tn warmer areas, a shorter snow cover season promotes warmer soils, particnlarly witliin snrface soil layers, dne to stronger air and soil thermal conpling. Additional analysis also indi cates that earlier snow cover seasonal onset in the fall has a stronger warming effect on modeled soil temperatnres in colder climate areas, while earlier offset of seasonal snow cover in the spring has a stronger warming effect on modeled soil temperatnres in warmer climate areas.

Climatic control on simnlated carbon llnxes
The model simnlations indicated that air temperatnre has an overall dominant control on the two main components of the NEE tlnx (i.e., net primary prodnctivity, NPP, and R\^) across all pan-Arctic climate zones, while snow has a larger con trol on estimated annnal NEE flnxes in colder climate areas (Fig. 6). These resnlts indicate that warming generally pro motes vegetation photosynthesis and soil heterotrophic res piration in the pan-Arctic region. However, a rednced posi tive correlation between NPP and air temperatnre in warmer climate zones (mean Fair > 0 °C) also indicates that warmingindnced dronght may rednce vegetation prodnctivity to some extent (Kim et al., 2012;Yi et at., 2014). No signiflcant cor relation {p > 0. t) between NEE and air temperatme was ob served for most areas (mean Fair < 5 °C) dne to NEE being a residnal between two large flnxes (i.e., NPP and Fh) with similar temperatnre responses. A predominantly positive cor relation (mean R =0.32; p <0.1) between annnal NEE and SWE in colder regions (mean Fair < -4 °C) is mainly dne to a strong positive correlation (R > 0.60, p < O.Ot) between SWE and NEE flnxes dnring the cold season (November to April; Fig. S7). A deeper snowpack promotes warmer soil condi tions (Fig. 5b) and associated SOC decomposition and het erotrophic respiration, which contribntes signiflcantly to an nnal NEE, especially in colder climate areas (Zimov et at., 1996). No signiflcant correlation (p > 0 .t) between armnal SCE or SWE and warm-season (MJJASO) carbon flnxes was observed.
While snow cover has a negligible effect on total estimated carbon flnxes dnring the warm season, it has a strong control on the composition of soil Fh (Fig. 7). An overall, deeper snowpack promotes soil decomposition and respiration from deeper (> 0.5m ) soil layers while inliibiting contribntions from snrface (< 0.2m ) soil layers, especially in colder cli mate areas. This response is dne to a stronger warming effect of snow cover on deeper soil layers in colder areas (Fig. 5). Comparatively, even thongh air temperatnre has a strong con trol on total warm-season Fh flnxes, it has a limited effect on the contribntion of different soil depths to total soil decom position and respiration except in the warmer climate areas (mean armnal Fair > 0 °C). tn the cold season, a deeper snow pack also promotes soil decomposition in deeper (>0.2m) soil layers more than in snrface (0-0.2 m) soil layers.

Sensitivity of simulated soil thermal dynamies and soil earbon deeomposition to elimate variations
The model sensitivity analysis nsing different snrface me teorology inpnts indicated that warming and rednced snow cover extent promoted widespread ALD deepening across the pan-Arctic domain over the 1982 to 2010 record (Fig. 8). tn Enrasia, strong winter warming rednced model-simnlated SWE and SCE, while increasing winter precipitation gener ally increased SWE and SCE. tn North America, regional trends in winter snowpack and SCE were more variable dne to variable trends in winter air temperatnre and precipitation. Therefore, the resnlting model-simnlated trends in SWE and SCE based on varying temperatme and precipitation inpnts showed strong spatial heterogeneity across the pan-Arctic domain. The model sensitivity analysis based on varying temperatme inpnts alone indicated overall ALD deepening in permafrost areas, corresponding with widespread warm ing and rednced SCE. However, the sensitivity analysis based on varying precipitation alone showed more variable trends in the simnlated ALD resnlts. Areas with strong decreasing winter precipitation and snowpack trends, snch as interior Alaska and eastem Siberia, showed a decreasing ALD trend, attribnted to rednced snow insnlation effects. The resnlts also indicated that changing air temperatnre had an overall dom inant effect on the simnlated ALD trends, thongh changing precipitation also contribnted to ALD changes in some areas. The model sensitivity analysis indicated that varying pre cipitation acconnts for more of the change in the simnlated Rh contribntion from different soil depths (i.e., soil Rh frac tion; Figs. 9 and tO, and Fig. SB), which is consistent with the above resnlts indicating strong control of snow cover on the soil Rh fraction at different soil depths. The model sensitivity resnlts also indicated that changing air temperatnre has min imal impact on the simnlated soil Rh fraction, while increas ing (decreasing) winter snowpack in permafrost areas gener- ally corresponded to increasing (decreasing) soil Rh fraction from deeper (>0.5m ) soil layers and decreasing (increas ing) soil Rh contribntions from snrface (0-0.2 m) soil lay ers (Fig. 9). This is particnlarly trae in cold climate regions (mean annnal Fair < -1 0 °C; Fig. 10). The simnlated Rh frac tion from the deeper soil layers ( < -10 °C) showed stronger deereasing trends in the fraetion from surfaee soil layers and inereasing soil R\^ eontributions from deeper soil layers, likely due to inereasing winter preeipitation and snow eover (Figs. 8 and 9) and eonsistent with field studies involving snow eover manipulations and assoeiated impaets on soil respiration (e.g., Nowinski et al.,2010).

Impact of climate variations on soil active layer properties
Our results show that reeent strong surfaee warming trends in the pan-Aretie region have promoted widespread soil thaw ing and ALD deepening (Fig. 8), while ehanging preeipi tation and snow depth have had a relatively smaller impaet on ALD trends (Figs. 4 and 8). We find a mean inereasing ALD trend of 0.66 ± 1.20 em yr"^ aeross the pan-Aretie re- gion over the past 3 deeades, whieh is similar to values re ported in previous studies Romanovsky et al., 2010a), albeit representing different time periods. This overall ALD deepening trend aeross the pan-Aretie domain eorresponds with widespread warming and warming-indueed deereases in SCE (Fig. 4e) and inereasing non-frozen-season duration (Kim et al., 2012). Our analysis indieates that air temperature has a dominant eontrol on upper (<0.5m ) soil layer temperatures during the warm-season, with an inereas ing eontrol in warmer elimate zones (Fig. 5a). The model simulations also suggest that most pan-Aretie permafrost ar eas, espeeially eontinuous permafrost areas, have a relatively shallow (< 1 m) aetive layer (e.g.. Fig. 4a). Therefore, rapid warming of the upper soil layers eorresponds with general AFD deepening.
Previous studies have also shown that summer air tem perature is a primary eontrol on AFD trends, while the re lationship between snow eover and AFD is more variable Romanovsky et al., 2010a, b). Our re sults demonstrate that deeper snowpaek eonditions promote warming of deep (>0.5m ) soil layers, espeeially in eolder elimate areas (Fig. 5b), and this effeet exeeeds the impaet of surfaee warming on deeper soil layers (e.g., > 1 m). Pre vious studies indieate that ehanges in snow depth ean influenee borehole (10-20 m) permafrost temperatures as mueh as ehanges in air temperature (Stieglitz et al., 2003;Ro manovsky et al., 2010a, b). Regional simnlations from the improved Community Fand Model (CFM) also indieate that snow state ehanges ean explain 50 % or more of soil tem perature trends at Im depth over the reeent 50-year reeord (Fawrenee and Slater, 2010). On the other hand, the impaet of ehanging snow eover duration on soil temperatures may vary aeross different elimate zones (Fig. 5e) due to the influenee of both air temperature and preeipitation or snow fall on snow eover duration. A shorter snow eover season may eool the soil in eolder elimate zones due to less in sulation from eold temperatures but may warm the soil in warmer elimate zones by promoting greater atmospherie heat transfer into soils (Fawrenee and Slater, 2010;Euskirehen et al., 2007). Our results indieate that reeent regional trends to ward eontinued warming, earlier spring snowmelt onset and a shorter snow eover season will likely enhanee soil warm ing and permafrost degradation in relatively warmer (mean annual Fair > -5 °C) regions of the pan-Aretie domain.

Impact of climate variations on soil carbon dynamics
Snow eover is an important eontrol on the annual earbon bud get in eold regions (annual mean Fair<-4°C; Fig. 6b-e), even though air temperature has a dominant eontrol on both annual NPP and Rh fluxes aeross all elimate zones (Fig. 6a). Strong snow eover buffering of underlying soil temperatures sustains soil respiration even under very eold winter air tem- (b) m e a n Tair <-10°C -1 0°c < m e a n ra irr< 0°C

Biogeosciences
'ii m e a n ^i r > -f)°C ■ Runi BRun2 RunS 0.08 > 0.06 -■ Runi iR u n 2 RunS response is due to different controls of snrface air temper atnre and snow cover on soil temperatnres at different soil deptlis (Zliang, 2005;Romanovsky et at., 20t0a, b;Lawrence and Slater, 2010). Snrface warming dnring tlie snmmer lias a dominant control on upper soil layer temperatmes (< 0.5 m; Fig. 5a), while a deeper winter snowpack has a persistent warming effect on deeper soil temperatnres in colder cli mate areas (Fig. 5b;Gonttevin et at., 2012). Therefore, smface warming likely promotes more heterotrophic respiration from snrface litter and soil layers, while a deeper snowpack promotes soil respiration from deeper soil layers. This is par ticnlarly important for soil carbon dynamics in permafrost areas, where a large amonnt of soil carbon occurs in deep perennial frozen soils (Hngelins et al., 2014). Previons stnd ies inclnding field experiments have primarily focused on the effects of snrface wanning on permafrost soil carbon decom position (e.g., Schnmet al., 2007;Koven et al., 2011;Schae fer et al., 2011), while onr resnlts show that snow cover may play a larger role than air temperatnre in iirflnencing deeper soil temperatnres and permafrost stability. This is also sup ported by a recent snow addition experiment in Alaskan tnn dra areas (Nowinski et al., 2010), which showed that a deeper snow treatment resulted in a larger contribntion of deep and old soil carbon decomposition to total soil heterotrophic res piration.

Limitations and nncertainties
peratmes, and the resnlting winter soil respiration can be a large component of the armnal NEE bndget (Sullivan et al., 2010). Field experiments have shown that winter soil respira tion in tnndra areas can offset total net carbon uptake dnring the growing season and thus switch the ecosystem from a net carbon sink to a carbon source (Zimov et at., 1996;En skirchen et at., 2012;Luers et at., 2014). Om resnlts also indicate that cold-season (November-April) Rh acconnts for ~ 25 % of total armnal Rh over the entire pan-Arctic domain, while this estimate may be conservative since onr model may nnderestimate soil respiration in tnndra areas (Fig. lb). The model simnlations indicate very low (< 5 %) mrfrozen wa ter below ~ -3 °C at the tnndra sites, while previons stnd ies and the tower measnrements (Fig. S4) indicate that sub stantial unfrozen water may remain even nnder very low soil temperatnres (e.g., ~ -1 0 °C), snstaiiting soil microbial ac tivities (Romanovsky and Osterkamp, 2000). On the other hand, winter warming may change the depth and stractme of insulating snow cover, affecting nnderlying soil tempera tnres, which could alter soil N mineralization rates and soil microbial activities that iirflnence ecological processes dnr ing the growing season (Schimel et al., 2004;Stnrm et al., 2005;Monson et al., 2006).
Even thongh air temperatnre has a dominant control on Rh dnring the warm season (from May to October), snow cover strongly iirflnences the contribntion of different soil depths to total soil decomposition and Rh (Fig. 7). This nonlinear Although soil temperatnre and moistnre are the two major environmental controls on soil carbon decomposition, other factors may also iirflnence soil decomposition rates and per mafrost carbon feedback potential bnt are not represented by onr modeling stndy (Hobbie et al., 2000). A number of chemical and biological factors can affect the temperatnre sensitivity of soil carbon decomposition in northem soils, in clnding enzyme abundance, microbial population size and oxygen availability (Waldrop et al., 2010). Previons stnd ies also show that soil carbon decomposition rates may be depth-dependent. Acconnting for vertical changes in soil bio geochemical properties and processes (inclnding the size and substrate quality of the soil active layer and permafrost car bon pool, and the degree of N mineralization with decompos ing permafrost carbon) may have significant impacts on the sign and magnitude of the projected high-latitnde carbon re sponse to future warming . Finally, changing wintertime soil microclimate will alter the amonnt and timing of plant-available nutrients {N) in tnndra ecosys tems and may drive a positive feedback between snow, soil temperatme, microbial activity, and plant community com position (Schimel et al., 2004;Stnrm et al., 2005).
A number of processes, notably fire disturbance, shmb ex pansion and thermokarst, are not inclnded in this stndy bnt may be important factors affecting regional permafrost and soil carbon dynamics (Grosse et al., 2011;Schnnr et al., 2015). A warming climate has been linked with increasing boreal-arctic fire activity and severity (Grosse et al., 2011). Fire can change the snrface vegetation composition and consnme a large portion of the soil organic layer, which can dra matically alter the snrface energy balance and soil thermal properties, and canse rapid permafrost degradation (Harden et al. 2006;Jafarov et al., 2013). Both field experiments and satellite measnrements indicate a "greeiung" Arctic with in creasing shmb abnndance dne to climate warming (Tape et al., 2006). Shmb expansion in Arctic tnndra can change the snow distribntion and snrface albedo, affecting the snrface energy balance and nnderlying active layer and permafrost conditions (Stnrm et al., 2005). The development of snrface water ponding with thermokarst in ice-rich permafrost areas can alter the local snrface hydrology, affecting permafrost and soil carbon decomposition (Schnnr et al., 2007;Grosse etal., 2011).
Another important featnre of the Arctic is strong snrface heterogeneity, characterized by widespread lakes, ponds, wetlands and waterlogged soils as a resnlt of both topography and restricted snrface drainage dne to nnderlying permafrost. Changes in both snrface and snbsnrface hydrology are tightly conpled with local permafrost conditions and potential car bon and climate feedbacks (Smith et al., 2005;Watts et al., 2012;Yi et al., 2014;Schnnr et al., 2015). Cnrrent largescale model simnlations, inclnding this stndy, generally oper ate on the order of tens of kilometers or even larger, and may not adeqnately represent the effects of snrface heterogeneity on simnlated permafrost hydrologic processes and soil car bon decomposition processes Rawlins et al., 2013;Schnnr et al., 2015). For example, most mod els prescribe a dominant vegetation type or a single valne for the orgaruc layer thickness commensnrate with the model spatial resolntion, which likely introdnces large nncertain ties to the model-simnlated moistnre and heat llnxes and thns the permafrost properties. Next generation satellites, inclnd ing the NASA SMAP (Soil Moistnre Active Passive) mis sion provide for finer-scale (i.e., 3-9 km resolntion) monitor ing and enhanced (L-band) microwave sensitivity to snrface (~< 5 c m ) soil freeze-thaw and moistnre conditions (Entekhabi et al., 2010) and may enable improved regional hy drological and ecological model parameterizations and sim nlations that more accmately represent active layer condi tions. Finer-spatial-scale observations nsing lower-freqnency (snch as P-band) synthetic apertnre radar (SAR) measnre ments from airborne sensors snch as AirMOSS (Airbome Microwave Observatory of Snbcanopy and Snbsnrface instrament; Tabatabaeenejad et al., 2015) may also provide improved iirformation on snb-grid-scale processes and snb snrface soil thermal and moistme profiles, providing critical constraints on model predictions of soil active layer changes and soil carbon and permafrost vnlnerability.

Conclusions
We developed a conpled hydrology and terrestrial carbon llnx modeling framework to evalnate the sensitivity of soil ther mal and carbon dynamics to snow cover and recent climate variations across the pan-Arctic basin and Alaska dnring the past 3 decades . Onr resnlts indicate that smface warming promotes widespread soil thawing and active layer deepeiung dne to a strong control of snrface air temper atnre on npper (<0.5 m) soil temperatnres dnring the warm season (from May to October). Recent trends indicating ear lier spring snowmelt and shorter seasonal snow cover dnra tion with regional warming (Dyer and Mote, 2006;Brown and Robinson, 2011;Kim et al., 2012) will most likely en hance soil warming in relatively warmer climate zones (mean aimnal Fair > -5 °C) and promote permafrost degradation in these areas. Even thongh air temperatnre has a dominant con trol on soil decomposition dnring the warm season, snow cover has a strong control on the contribntion of different soil depths to the total soil heterotrophic respiration llnx. A deeper snowpack inhibits snrface (< 0.2 m) litter and soil or garuc carbon decompositionbnt enhances soil decomposition and respiration from the deeper (>0.5 m) soil carbon pool. This nonlinear relationship between snow cover and soil de composition is particnlarly important in permafrost areas, where a large amoimt of soil carbon is stored in deep peren nial frozen soils that are potentially vnlnerable to mobiliza tion and accelerated losses from near-term climate change. Onr resnlts demonstrate the important control of snow cover in affecting active layer properties and soil carbon decom position processes across the pan-Arctic and the necessity of considering both warming and a change in precipitation and snow cover regimes in characterizing permafrost soil car bon dynamics. In addition, further improvements in regional assessment and monitoring of precipitation and snow cover across the northem high latitndes are needed to improve the quantification and nnderstanding of linkages between snow and permafrost carbon dynamics. The PWBM model (Rawlins et al., 2013) simulates snow and gronnd thermal dynamics by solving a 1-D heat trans fer eqnation with phase change (Nicolsky et al., 2007):

Biogeoseienees
where T(z, t) is the temperatnre (°C) and C(T,z) and X(T, z) are the volnmetric heat capacity (Jm "^K "^) and thermal condnctivity (W m"^ K "^) of soil, respectively; L is the vol nmetric latent heat of the fusion of water (Jm "^); i; is the volnmetric water content, and 6 is the mrfrozen liqnid water fraction. The Dirichlet boundary conditions at the snow or gronnd snrface Zs, i.e., T (zs, t) = Tair(f), and a heat boundary condition at the lower boundary zb, i.e., X-^TQ, t) = g, were nsed to solve the heat eqnation, where Tair is the observed air temperatnre and g is the geothermal heat llnx (K m "^). The volnmetric water content (C) can be obtained by solving the Richard's eqnation. The mrfrozen liqnid water fraction (0) was estimated empirically as where the constant T* is the freezing point depression, and b is a dimensionless parameter obtained from unfrozen water curve fitting (Romanovsky and Osterkamp, 2000). The bulk thermal properties of soil (i.e., C and X) are a combination of the thermal properties of soil solids, air, and thawed and frozen states of soil water (Rawlins et al., 2013). Particnlarly, for the soil solids, the volnmetric heat capacity (Cs) and thermal condnctivities (Xs) vary with the fraction of organic carbon of the soil, defined as where / is the fraction of orgaruc carbon in the soil. Cm and Co are the volnmetric heat capacities of the mineral and or ganic soils, respectively, and and Xg are the thermal con dnctivities of the mineral and organic soils, respectively.
Up to five snow layers were nsed to characterize the snowpack dynamics and solve the snow temperatme profile, with varying depth for each layer depending on the snow depth. A two-layer snow density model similar to Schaefer et al. (2009) was nsed to characterize the impact of the bottomdepth hoar layer on the snow thermal condnctivity for tnn dra and taiga, with fixed snow thermal condnctivity for this layer. For the npper snow layer, both the snow heat capacity and thermal condnctivity vary with snow density. Following Liston et al. (2007), the temporal evolntion of the snow den sity is mainly affected by new snowfall and compaction dne to winds: dP; dt (A4) where ps is the snow density (kgm "^), U represents the wind-speed contribntion to the snow density changes with negligible iirflnence for wind speed below 5 m s "^; Tf and Ts are the freezing and snow temperatmes, respectively; a\, fl2 and h me empirical dimensionless pmameters. The snow thermal condnctivity (Xsnow) is an empirical estimate of snow density based on Stnrm et al. (1997): = 0.138 -l.Olps + 3.233p," More details on the numerical solution of the heat transfer eqnation and the pmameterization of the snow model can be found in Rawlins et al. (2013) and Nicolsky et al. (2007).

A2 Carbon model description
A satellite-based light use efficiency (LUE) approach was nsed to estimate vegetation prodnctivity: where GPP is the gross primary prodnctivity (g C m "^ d"^); e (gC M J"^) is the LUE coefficient converting absorbed photosynthetically active solm radiation (APAR) to vege tation biomass, and FPAR defines the fraction of incident PAR (M Jm "^ d"^) absorbed by the vegetation canopy (i.e., APAR). A maximnm LUE coefficient (smax, g C M J"^) was prescribed for each land cover type and was rednced for snboptimal environmental conditions (inclnding low air temper atnre, soil moistme and frozen conditions) to estimate e (Yi et al., 2013). Vegetation net primary prodnctivity (NPP) was es timated as a fixed portion of GPP for each biome type based on an assumption of conservatism in vegetation cmbon use efficiency within similar plant functional types. A dynamic carbon allocation of litterfall estimated from NPP, based on Randerson et al. (1996) and White et al. (2000), was nsed to characterize litterfall seasonality. The total litterfall was partitioned into three components, in clnding leaves, fine roots, and woody components with pre scribed ratios for each plant functional type based on field experiments (White et al., 2000; Table S2). Daily constant tnmover rates were prescribed for the woody components of litterfall inclnding stems and coarse roots (White et al., 2000), while the NDVI time series were nsed to character ize tnmover rates of the other two variable components of litterfall dnring leaf senescence and active growth periods (Randerson et al., 1996). Approximately half of the fine root tnmover was assumed to occur dnring the active growing sea son, and the monthly variable fraction of litterfall was calcn- -[NDVI(f + 1) + 0.5 ■ NDVI(f + 2)], where LTvari(?) and LTvar2 ( 0 represent the litterfall frac tion associated with leaf loss (i.e., LL(t)) and vegetation active growth, respectively; LTieaf and LTfroot are the pre scribed fractions of leaf and fine-root components for each plant functional type, respectively (Table S2). The estimated monthly litterfall fraction was then distribnted evenly over the month.
To acconnt for the contribntion of deep soil organic carbon pools to the total heterotrophic respiration {R\x), we extended the original terrestrial carbon llnx (TCF) soil decomposition model to incorporate soil organic carbon down to 3 m below the snrface, and multiple litter and soil orgaruc carbon (SOC) pools were nsed to characterize the progressive decompo sition of fresh litter to more recalcitrant materials. Follow ing Biome-BGC (BioGeochemical Cycles; Thornton et al., 2002), the new soil decomposition model includes three lit terfall pools, 3 SOC pools with relatively fast tnmover rates and a deep SOC pool with slow tnmover rates (Fig. SI). The litterfall carbon inpnts were first allocated to the three litter fall pools according to the substrate quality of each litterfall component, i.e., labile, cellulose and lignin fractions of esti mated leaf, fine root, and woody litterfall (Table S3; White et al., 2000), and then transferred to the SOC pools throngh progressive decomposition.
For each carbon pool (C;), the carbon balance of the de composition process was defined as (AS) where Ri is the carbon inpnt from litterfall allocated to pool i (only nonzero for the 3 litterfall pools), Tji is the fraction of carbon directed from pool j to pool i with fraction rj lost as respiration, and ki {k/) is the decomposition rate of carbon pool i{j). The heterotrophic respiration (R^) is then com puted as the sum of respiration flnxes from the decomposi tion process: The soil decomposition rate (ki) for each pool is derived as the product of a theoretical maximnm rate constant (kmaxj. Fig. SI) and dimensionless mnltipliers for soil temperatnre (Tmuit) and moistnre (W muit) constraints to decomposition nnder prevailing climate conditions: where Tmuit and Wmuit vary between 0 (fully constrained) and 1 (no constraint), as defined in Yi et al. (2013).