Estimation of Permafrost SOC Stock and Turnover Time Using a Land Surface Model With Vertical Heterogeneity of Permafrost Soils

We developed vertically resolved soil biogeochemistry (carbon and nitrogen) module and implemented it into a land surface model, ISAM. The model captures the vertical heterogeneity of the northern high latitudes permafrost soil organic carbon (SOC). We also implemented Δ14C to estimate SOC turnover time, a critical determinant of SOC stocks, sequestration potential, and the carbon cycle feedback under changing atmospheric CO2 concentration [CO2] and climate. ISAM accounted for the vertical movement of SOC caused by cryoturbation and its linkage to frost heaving process, oxygen availability, organo‐mineral interaction, and depth‐dependent environmental modifiers. After evaluating the model processes using the site and regional level heterotrophic respiration, SOC stocks, and soil Δ14C profiles, the vertically resolved soil biogeochemistry version of the model (ISAM‐1D) estimated permafrost SOC turnover time of 1,443 years, which is about 3 times more than the estimation based on the without vertically resolved version of ISAM (ISAM‐0D). ISAM‐1D‐simulated SOC stocks for permafrost regions was 319 Pg C in the top 1 m soil depth by the 2000s, about 80% higher than the estimates based on ISAM‐0D. ISAM‐1D SOC stock and turnover time were compared well with the observations. However, the longer SOC turnover time preserves less SOC stocks due to the lower carbon use efficiency (CUE) for SOC than ISAM‐0D and thus respires more SOC than being transferred downward by cryoturbation. ISAM‐1D simulated reduced SOC sequestration (3.7 Pg C) compared to ISAM‐0D (4.8 Pg C) and published Earth system models (ESMs) over the 1860s–2000s, due to weaker [CO2]‐carbon cycle and stronger climate‐carbon cycle feedbacks, highlighting the importance of the vertically heterogeneous soil for understanding the permafrost SOC sinks.


Introduction
The northern high latitudes permafrost soils store more than half of the global total soil organic carbon (SOC) (approximately 1,300 Pg C; Hugelius et al., 2014) and show extreme spatial and vertical heterogeneity (Beer, 2016;Mishra et al., 2017) due to various cryopedogenic processes (Ping et al., 2008(Ping et al., , 2015, including the vertical movement of SOC due to seasonal freeze-thaw cycles and associated movement of water and soil particles in all dimensions (Beer, 2016;Nicolsky et al., 2008;Woo, 2012). These vertically heterogeneous and waterlogging permafrost conditions slow SOC decomposition and humification of the plant litter input, resulting in a lower amount of loss of SOC stock from the upper soil layers and relatively longer SOC turnover time, which is the time it would take for the SOC pools to empty without further inputs as defined and described in supporting information Text S6. Under changing environmental conditions these processes affect two basic feedbacks and hence SOC stocks (Friedlingstein et al., 2006;He et al., 2016): (1) carbon-concentration (or CO 2 fertilization) feedback, which enhances plant productivity and litter input flux under the higher atmospheric CO 2 concertation ([CO 2 ]) conditions and hence increase SOC, and (2) carbon-climate feedback, which increases SOC losses due to higher decomposition under warmer climate and hence reduce SOC.
However, the uncertainties in model-based estimates of the effect of these two feedback processes and hence SOC storage amount in the permafrost regions are large (Burke et al., 2012;Harden et al., 2012;Koven et al., 2011;MacDougall et al., 2012;Schaefer et al., 2014). This is mainly caused by the single soil layer with multiple pools structure of the terrestrial ecosystem models (TEMs) (Clark et al., 2011;Jain et al., 2009;Jain & Yang, 2005;Melillo et al., 1993;Sitch et al., 2003;Thornton et al., 2007), which is The objective of this study is to estimate the spatial and the vertical distribution of SOC stock and turnover time in the permafrost regions using a land surface model, ISAM (Barman et al., 2014a(Barman et al., , 2014bEl-Masri et al., 2015;Yang et al., 2009). This study builds upon and extends the approaches of the modeling studies discussed above. While ISAM considers the similar soil and snow processes and their interactions with biogeophysics and biogeochemistry under various environmental conditions in the permafrost regions, as accounted for in the vertically unresolved soil biogeochemistry version of ISAM (Barman & Jain, 2016;El-Masri et al., 2015), in this study we implement in ISAM more advanced parameterizations for calculating the cryoturbation of SOC. These parameterizations explicitly estimate ice lens formation and link the cryoturbation rate to the freeze-thaw cycle of the soil above the permafrost table. In addition, we implement Δ 14 C to calculate SOC turnover time based on radioactive decay based on inputs of Δ 14 C from atmospheric "bomb carbon" (He et al., 2016). After calibration and evaluations of the newly implemented processes, as well as examine the uncertainty in the model estimated SOC stock and turnover time due to uncertainties in the model parameters, the improved model was applied to estimate SOC stock and turnover time of the permafrost regions due to [CO 2 ] and climate change and analyze their changes over the historical time period 1861-2010.
In the following section 2, first, we describe ISAM original version (refer to as ISAM-0D hereafter) and extended version (refer to as ISAM-1D hereafter). Then in section 3, we describe the data used to calibrate and evaluate model estimated SOC stock and turnover time. Section 4 describes the results for model calibration and evaluation, and model estimated SOC stocks and turnover time results for the permafrost regions over the period 1901-2010. In section 5 we discuss the importance of model improvements for calculating the vertical heterogeneity of SOC, as well as some model limitations and conclude our major findings.

Model Description
The ISAM model applied here to study permafrost SOC dynamics in the permafrost regions is a fully coupled biogeophysical and biogeochemical cycle model with 0.5°× 0.5°spatial and multitemporal resolution from half-hour to year (El-Masri et al., 2013;Jain et al., 2009Jain et al., , 2013Yang et al., 2009). Each grid cell is occupied by a combination of fractional vegetation, bare soil, and glacier (Meiyappan & Jain, 2012). The original version of the model, ISAM-0D, includes a comprehensive description of the biogeophysical (energy and water) processes, which are resolved with 15 discretized soil layers extending up to 50 m. The five bottom layers, which extend between 3.5 and 50 m, are the deep thermally active soil layers, which are added to avoid numerically biased estimation of the deep soil temperature and only active for the soil thermal conduction calculation. All other calculations, including soil hydrology, soil carbon, and nitrogen, are performed in the first 10 layers, which extend up to 3.5 m. ISAM biogeochemical (carbon and nitrogen) processes resolving across the atmosphere-vegetation-soil with a single 1 m soil layer (Barman et al., 2014a(Barman et al., , 2014bEl-Masri et al., 2015) ( Figure 1). The biogeochemistry in ISAM-0D are cycled among the vegetation (tree and aboveground foliage and tree-woody biomass), aboveground litter (dead foliage and dead wood), belowground litter (dead roots), aboveground and belowground SOM (microbial biomass and hummus), and inorganic N (NH 4 + and NO 3 − ) pools   (Figure 1). Carbon assimilation, heat, and water fluxes are calculated through coupled canopy photosynthesis, energy, and hydrological processes. Carbon assimilation is allocated into different vegetation pools and then transfer to litter and SOM pools. The C cycle is then coupled with a complete N cycle, which accounts for major N processes, including N deposition, N fixation, N mineralization, N immobilization, nitrification, denitrification, and leaching (Yang et al., 2009).
The model incorporates several structural extensions and processes in recently published papers (Barman & Jain, 2016;El-Masri et al., 2015), which are important for describing the distinct biogeophysical and biogeochemical processes in the permafrost region; wind compaction of snow and depth hoar formations (Barman & Jain, 2016) by updated thermal conduction and soil water movement processes. The biogeochemical cycle component accounts for C-gain based dynamic plant phenology, dynamic rooting distribution and depth parameterizations (El-Masri et al., 2015).

Extended Version of ISAM
Here we extend the single soil layer of the biogeochemistry component of ISAM-0D to the same 10 discretized soil layers extend up to 3.5 m depth. The carbon and nitrogen pools in each of the 10 soil layers are the same as has been considered in the ISAM-0D (Figure 1). The major permafrost region plant functional types (PFTs) considered in the extended model, ISAM-1D, are shown in Figure S1. The ISAM-1D includes three major processes essential to better represent the vertical distribution of SOC stocks and turnover time in the permafrost region. These processes are (1) bioturbation and cryoturbation, including the formation of ice lens-induced soil transport processes, as well as the response of permafrost systems to changing environmental conditions; (2) vertically resolved SOC decomposition and related environmental modifiers, including soil temperature, soil moisture, oxygen availability, and the strength of organo-mineral interaction; and (3) moss insulation due to changing soil temperature and soil state heterogeneity, including moisture and texture. Besides these processes, a 14 C tracer is incorporated to evaluate the model estimated SOC turnover time ( Figure 1). In the following, we provide a brief description of the individual process added to ISAM-1D; the detailed description can be found in Text S1. The description of the newly added variables and equations in the model are provided in Tables S1 and S2.

Global Biogeochemical Cycles
The vertical transport of SOC among the soil layers are accounted for by the bioturbation and cryoturbation processes. In the model bioturbation is represented by the diffusion term and cryoturbation by advection terms (Elzein & Balesdent, 1995, Equations S1 and S2) (Text S1). Bioturbation, which is the movement of soil particles through living soil organisms, matches the diffusion process that determines the movement direction through the SOC gradient from topsoils with high SOC concentration to the deeper soils. Cryoturbation, on the other hand, has several forms based on the particle moving direction, including horizontal, vertical, and mixed (Rempel, 2007). In this study, we consider cryoturbation in the vertical direction, because the vertical transport serves as an important mechanism to protect SOC (Ping et al., 2015). The cryoturbation process in this study is represented by advection (Equation S4, Rempel, 2007) as a function of the ice lens growth rate (Equation S5). This is because cryoturbation can bring carbon from a low carbon density layer to a high carbon density layer, depending on the direction of the ice lens nudging. This SOC gradient effect can only be accounted for well by the advection rate. In addition, the model accounts for the cryoturbation rate (ν) as a linear function of the movement rate of the ice lens with γ as the slop parameter of the liner function (Equation S4). This parameterization is useful for accounting for the effect of the cryoturbation on the SOC profile instead of implementing the detailed microscopic ice-soil particle interaction processes (Rempel, 2007). A depth limitation parameter (Z maxcryo ) is also added in the model, which prevents the penetration of SOC into deeper soils due to cryoturbation. The ice lens growth rate is calculated using a rigid-ice module (O'Neill & Miller, 1985) (Text S2), which simulates the uplift of water-saturated soil due to expansion caused by the ice lens growth. To calculate a more realistic rate of ice lens growth with time, the hydraulic conductivity in frozen soils accounts for the effects of ice within the soil matrix by using an ice impedance factor (Equation S7, Swenson et al., 2012). This cryoturbation scheme improves the modeled changes in the cryoturbation rate under different climate conditions. For example, a drier climate may slow the cryoturbation rates due to limited soil water availability for ice lens growth, leading to a lower frost heave rate and cryoturbation activity (Ping et al., 2015). The model accounts for a reduction in soil respiration due

10.1029/2020GB006585
Global Biogeochemical Cycles to frozen conditions using an empirical temperature modifier as a function of frozen Q 10 temperature coefficient (Equations S11 and S12) ; see Text S4 for detailed description).
The model grid cells can have frozen and unfrozen fractions of SOC, and the model treats decomposition and vertical oxygen diffusion separately for the frozen and unfrozen soils. The decomposition in the frozen soil is controlled by the frozen Q 10 temperature coefficient (see Text S4) and assumes no oxygen diffusion. The model calculates O 2 concentration in each soil layer for both saturated and unsaturated soil conditions (Equations S15-S21) by accounting the process of oxygen diffusing into the soil from the atmosphere and the soil oxygen consumption through heterotrophic respiration (Wania et al., 2010) (Text S5). To determine the interface between the soil surface and air-water (Swenson et al., 2012), we incorporated the perched water level (Equation S23) on the permafrost during the early growing season, that is, the beginning of the vegetation activity, due to the thawing of the ice.
We have added three different types of SOC processes into the vertically resolved model structure: (1) SOC decomposition and degradation, (2) vertically resolved litter input, that is, the flux entering the soil pools based on the above ground and below ground plant dead material, and (3) SOC accumulation in soil layers once the SOC exceed the threshold value described by the bulk density (Equations S8-S10 and Text S3). The model treats the SOC decomposition, degradation processes, and the corresponding temperature (Equations S11-S12), moisture (Equation S13), and the environmental modifiers in the model's C-N coupled scheme ) at each soil layer separately. We have added two new modifiers: (1) A depth modifier (Equation S14, see detailed description in Text S6) accounting for the resistance of SOC decomposition due to organo-mineral interaction (Masiello et al., 2004;Torn et al., 1997), and (2) an oxygen modifier (Equation S15, Riley et al., 2011;Koven et al., 2013) to slow down the SOC decomposition due to oxygen limitation. Text S4 provides a detailed technical description of all the environmental modifiers.
To improve modeled changes in soil microbial activities and active layer thickness, the revised model considers the top first 10 cm soil layer as a moss layer, containing the accumulated amount of poorly decomposed litter. The thermal and hydraulic conductivities of this soil layer are replaced by the moss conductivity values (O'Donnell et al., 2009). The added moss layer, which covers most of the boreal and tundra areas, serves as a thermal insulator in the permafrost region, restricting seasonal soil temperature changes; a cooler soil temperature during the summer inhibiting the soil microbial activity, while warmer soil temperature during the early winter prolonged the soil microbial activity (Text S7 for details). When the SOC is transferred downward through the processes discussed above, and the calculated stored amount of SOC (Equation S8) in the SOC layer exceeds the maximum limit (e.g., the Bemmelen conversion factor, which is 58% of C content for pure organic soil), the excess amount is transferred to the deeper soil layer (Equation S9) (Lorenz & Lal, 2005). Without accounting this transfer, ISAM-1D overestimated the vertical transfer of soil thermal energy, resulting in the underestimation of the active layer thickness, surface organic layer thicknesses, SOC storage, and Δ 14 C.
The model estimated Δ 14 C, described as the ratio of 14 C and 12 C (as SOC), soil profiles are used to calculate SOC turnover time (Torn et al., 2009; described in Text S6). To calculate Δ 14 C (Equation S26, Stuiver & Polach, 1977), we accounted for 14 C in each SOC pool of ISAM, which is explicitly calculated by including a set of equations (Equations S1b and S8b). The equations for the 14 C are the same as for SOC, except that 14 C equations account for radioactive decay with a half-life (λ) of 5,568 years. In general, fractionations are part of the 14 C equations, but we assume no fractionation, because the observations data of Δ 14 C we used compare ISAM modeled Δ 14 C to be already fractionation-corrected . The decomposition rate and environmental modifiers for 14 C are assumed to be the same as for the SOC. For all 14 C simulations, we forced the model with a linearly interpolated long term atmospheric 14 C time series for the period of 30 kyr before 0 BP (before present) (Reimer et al., 2013) to 2010 (Hua et al., 2013).
The more detailed light reaction scheme is incorporated in ISAM-1D for better estimation of the CO 2 fertilization effect in the permafrost region. The light reaction of plant photosynthesis in the ISAM-0D model accounts for electron transport rate (J, mol m −2 s −1 ) calculation using a linear function of the light response curve (Chen et al., 2011). In ISAM-1D, we considered the electron transport rate to be limited by the maximum capability (J max25 ) due to the limitation of enzyme activity, which is dependent on leaf temperature. The response curve of J max25 to leaf temperature is fitted to the measured response curve from Yamasaki et al. (2002) using a parabolic equation (Text S8 for detail).

Field Observation Data Set
Two different field data sets are compiled for the model calibrations and evaluations. First, the data for modern SOC and Δ 14 C profiles from five nonpeat sites in the permafrost region (He et al., 2016, See also Figure S2) are used to calibrate and evaluate the model parameters (described in section 3.2). These sites cover a variety of climate conditions for boreal forest and tundra biomes in the permafrost region. More basic site information can be found in Table 1. We selected three specific sites (Boreal Forest: Donnelly Flats and Hess Creek, and Shrub/Tundra: Gydansky Western Siberia) data for SOC and Δ 14 C profiles to calibrate and two boreal forest sites (Old Black Spruce and Frost Fire Black Spruce) data to validate eight model parameters listed in Table 2. All these parameters are used to calculate vertical SOC and soil Δ 14 C distributions in ISAM. (Table S3).
Second, we compiled 354 samples of Gelisol SOC profiles Vitharana et al., 2017) for evaluating the modeled SOC stock and turnover time at a regional scale. These profiles cover all major ecoregions (Omernik & Griffith, 2014) in permafrost regions ( Figure S2), with SOC stock ranging from 2-120 kg C m −2 for the top 1 m. Since the model does not account for peat soil, we first exclude all Histels soil profiles from the compiled data. Next, we exclude soil samples from the coastal region where climate forcing data is not available to drive the model. Finally, we used the remaining 197 Gelisol samples ( Figure S2).
We also compiled the following auxiliary surface variable data for each site to drive the model: the surface PFT types, soil bulk density, and soil texture profiles (Michaelson et al., 2013;Mishra et al., 2017). In case   (2003) Note. The listed low-and high-end range values are taken from the literature. The ISAM-1D calibrated values are the mean for five sites (see Table S3 for values of the individual site) and have been used for the rest of the study. a For estimating the range of γ, we examine the range of cryoturbation rate by varying γ and select the γ value that produces a 5 times bigger cryoturbation rate  compared to the diffusive rate. b We used the standard deviation of SOC observation from 197 soil samples  to estimate the bulk density (bd) range and applied it to the empirical relationship of Périé and Ouimet (2008) to constrain the uncertainty for bulk density of SOC. the data for one or more variables is missing at particular sites, we used the Global Soil Dataset for Earth system modeling (GSDE, Shangguan et al., 2014) data. We then used a mass-preserving interpolation algorithm (Malone et al., 2009) to calculate the SOC content and soil bulk density at each soil layer.

Experimental Design for Model Calibration, Evaluation, and Uncertainty Analyses
The model calibration is performed in the following two steps: First, we spin-up the model to estimate steady-state vertical depth profiles for SOC and prebomb Δ 14 C for three calibration sites. In this step the model is forced by repeating the 1901-1920 CRUNCEP 6-hourly climate data with a fixed preindustrial (1860) atmospheric CO 2 concentration (286.41 ppm) and N deposition levels (Kanakidou et al., 2016).
We calibrated and evaluated the values of eight specific parameters (Table 2), controls all major processes, which we have incorporated in ISAM-1D, including cryoturbation, rooting depth, SOC vertical transport, the temperature sensitivity of decomposition, and organo-mineral interaction. These eight parameters are (1) diffusive rate (D, Equation S3) determining the strength of bioturbation, (2) slop parameter (γ, Equation S4) controlling cryoturbation rate as a function of the movement rate of ice lens, (3) maximum depth (z maxcryo , Equation S3) to which SOC can be penetrated to, (4) depth modifier slope parameter (s, Equation S14) controls the strength of organo-mineral-interaction, (5) temperature coefficient (Q 10 , Equations S11 and S12) controls the decomposition rates of labile litter, and SOC for each 10°temperature change, (6) soil bulk density (bd, Equation S9) determines the accumulation of organic layer, (7) root growth direction parameter (aa, Equation S22), and (8)  In the first step, we prescribed the values of eight parameters based on ISAM-0D. In the second step, the model is run iteratively from 1901-2010 by repeating the CRUNCEP climate forcing, atmospheric CO 2 concentrations (Friedlingstein et al., 2019) and historical N deposition (Kanakidou et al., 2016) data and the parameter values are adjusted within the predefined range, which we obtain based on published research articles (Table 2). We optimized the cost function value of each parameter, which is the root-mean-square error (RMSE) between the observed and modeled SOC and Δ 14 C profiles. The final parameter values are achieved by using the automatic FFSQP (FORTRAN implementation of the Feasible Sequential Quadratic Programming) numerical optimization package (Zhou et al., 1997). Under certain cases, FFSQP is unable to minimize the cost function. In such cases, the parameters are manually adjusted using the trial and error method.
To evaluate the SOC stocks and Δ 14 C profiles at the regional scale, we used the same model settings as described above and calculated SOC stocks, Δ 14 C, and SOC turnover time for 197 sites using ISAM-0D and ISAM-1D. These site-specific simulations are carried out using the mean of calibrated site values for all eight parameters (Table S3).
We determined the uncertainty in the model estimated SOC stocks and turnover time due to uncertainty in the eight model parameters. A total of 16 simulations (8 for low range and 8 for high range of values) were performed for all five sites (Table S3). For each simulation, we spin-up the model by repeating the climate and other soil data over the period 1901-2010, the same method as we have used for the calibration and evaluation cases. By changing each parameter values one at a time and keeping other parameter values fixed at the calibrated values, we obtained the model estimated high and low range of values of SOC stocks and turnover time.

Evaluation of the Model Results for Carbon Fluxes
After incorporating processes described above, we reevaluated ISAM-1D estimation of major carbon fluxes using site-level data sets for gross primary productivity (GPP) (Table S4), net primary productivity (NPP) and heterotrophic respiration (Rh) ( Table S5). All sites are located above 45°N. Specifically, the microbial CUE, the efficiency with which microorganisms transfer the available organic matter into labile and resistant pools, is calibrated to match the observed Rh (Text S9). The microbial CUE in ISAM is determined by the microbial efficiency parameter, which is the fraction of litter carbon transferring into the foliar and 10.1029/2020GB006585

Global Biogeochemical Cycles
woody microbial pools (Figure 1c) (Gougoulias et al., 2014). In the model we calibrate the fraction value using the site level Rh data as described in Text S9. The calibrated microbial CUE parameter for the labile pool is 0.3, which is 50% less compared to ISAM-0D, leading to a decrease of carbon sequestration ability due to an increase of the SOC respiration, particularly from the labile pool rather than transfer it to the resistant pools.

Estimating Changes in SOC Stock and the SOC Turnover Time for Permafrost Regions Over the Period 1901-2100
We conducted two regional simulations, one using ISAM-0D and another using ISAM-1D, to estimate the current SOC stocks in the entire global permafrost region ( Figure S2), which contains the entire permafrost and the boreal forest area in North America, Asia, and Europe. The ISAM-0D and ISAM-1D do not account for peat dynamics and the spodosol mechanisms. Therefore, our study region excludes the SOC from peatlands and spodosol. Both model simulations are forced with the same spin-up settings as for the model calibration case (i.e., repeating the 1901-1920 CRUNCEP climate data with 1860 level of atmospheric CO 2 concentration, 286.41 ppm, and N deposition). Next, the historical simulations were performed using 1901-2010 data for climate, atmospheric CO 2 concentrations, and N deposition. The SOC turnover time is then calculated as the SOC mass-weighted average of all layers above 1-m soil (Trumbore & Harden, 1997

Model Calibration and Validation Results
All five sites have distinct soil vertical structure and composition. Therefore, the vertical distribution of SOC are quite different for all five sites as represented by the large variations in calibrated values for two specific parameters, z maxcryo and γ, which are varying by 5 and 18 times across different sites (Table S3). The z maxcryo of 0.27 m at Hess Creek has the shallowest SOC profile (Figure 2e) with the lowest accumulation of labile SOC. In contrast, z maxcryo of 2.98 m value at Donnelly Flats Creek Control site suggests that the cryoturbation process is active even in the deeper soils. The γ is the highest at Donnelly Flats Creek Control site due to a higher frost heave rate resulted from faster ice lens growth. This is common for the silt soil to form high SOC content in the subsoil (Kaiser et al., 2007;Ping et al., 2015). The ice lens growth at other calibrated sites is small due to the lower cryosuction rates. Thus, a smaller calibrated γ value is sufficient for model simulation to match the observed SOC stocks and Δ 14 C profiles.
After the implementation of the frost heave process in representing cryoturbation, ISAM is able to capture explicitly the increased cryoturbation rate in response to the higher variation of soil moisture ( Figure S3) at calibrated sites with higher frost heave rates (Donnelly Flats Creek Control and Gydansky Western Siberia). This results in a relatively higher SOC stocks (Figures 2a and 2c) and Δ 14 C (Figures 2b and 2d) in subsoils at both sites compared to calibrated sites.
Despite a large heterogeneity at different sites, the ISAM-1D model is able to capture the observed trend of a higher Δ 14 C at the topsoil due to the accumulation of the large amount of the fresh litter residue with higher 14 C and relatively negative Δ 14 C at the subsoil due to downward transport of relatively lower amount of fresh litter input and decomposed SOC from topsoil to subsoil, which is resistant to further decomposition (Trumbore & Harden, 1997), at calibrated sites ( Figure 2) and evaluated sites (Figure 3). The reduction of Δ 14 C along with the soil depth indicates longer SOC turnover time in the subsoil compared to the topsoil. Also, the newly added input to the topsoil decomposes faster and contributes to the Rh.
Further evaluation of the ISAM-1D eight parameters, which were calibrated and evaluated using five sites data, for the mean SOC stocks for the top 1 m soil over the 197 sites show that the ISAM-1D results of 30.3 kg C m −2 are consistent with the Northern Circumpolar Soil Carbon Database (NCSCD), 30 kg C m −2 , and Mishra et al. (2017), 33 kg C m −2 data sets than the results of ISAM-0D, 15 kg C m −2 (Figure 4a). At the same time, ISAM-1D estimated mean permafrost SOC turnover time, which is derived from the observed Δ 14 C profile, of the 1,443 years for the 197 permafrost soil sample sites, which is 10.1029/2020GB006585 Global Biogeochemical Cycles also in much better agreement with the value of the 1,408 years derived from the observed Δ 14 C profile than ISAM-0D estimated value of 439 years (Figure 4b). The ISAM-1D estimates lower decomposition rates than ISAM-0D, resulting in more SOC amount in the topsoil. Excess amount of SOC is transported to the subsoil layers over time where it is stabilized due to the organo-mineral interactions and the colder temperature

10.1029/2020GB006585
Global Biogeochemical Cycles below the permafrost table, thus increasing SOC stock and SOC turnover time in subsurface soil layers. These results suggest that ISAM-1D can also be applied to the entire permafrost region for the calculations of SOC stocks and turnover time.

ISAM-1D Estimated SOC Stocks and Turnover Time for the Permafrost Region
The estimated top 1 m SOC over the permafrost region for the 2000s are 319 and 157 Pg C based on ISAM-1D and ISAM-0D (Figures 5a and 5b). Over 1860-2000s, ISAM-1D and ISAM-0D estimated a total increase of 1.2% (3.7 Pg C) and 3.9% (4.8 Pg C) (Table 3), mostly in boreal forest region (Figures 6a and 6b). These results suggest that ISAM-1D results are in much better agreement with NCSCD estimates of 327.36 Pg C (after excluding the peatland and spodosol SOC) than ISAM-0D results (Figures 5a-5c).
The split of 1 m SOC between the topsoil and subsoil for the year 2000 is 37% and 63% (Figures 5b, 5e and 5h) for ISAM-1D, and 30% and 70% (Figures 5a, 5d, and 5g) for ISAM-0D compared to NCSCD values of 43% and 57% (Figures 5c, 5f, and 5i), suggesting that both versions of ISAM estimated values for topsoil is underestimated and for subsoil overestimated compared to NCSCD data. However, the split of observed total SOC stocks (33 Pg C) from the 197 collected soil samples is 38% and 62%. These results are aligned with ISAM-1D estimated results, but about 5% lower (higher) than the topsoil (subsoil) estimates based on NCSD data. The ISAM estimates and these two data sets are agreed within the uncertainty ranges (6.5-55.2 kg C m −2 ) for NCSCD (6.7-53.7 kg C m −2 ) and 197 collected soil data sets (8.9-54.2 kg C m −2 ).
Most importantly, ISAM-1D estimated spatial distribution of SOC stock shows a better alignment with NCSCD spatial pattern than estimated based on ISAM-0D, especially in the eastern Siberia tundra region with low SOC content due to shallower rooting depth. ISAM-1D estimation also captured the high SOC for a boreal forest than for grassland and tundra in the eastern and southern Siberia region due to accounting the frost heave, organo-mineral interaction processes and deeper rooting depth (El-Masri et al., 2015).

Global Biogeochemical Cycles
In regards to SOC turnover time, ISAM-1D estimated 1 m Δ 14 C-based SOC turnover time for the 2000s is 486 years (Figure 7b) with much higher SOC turnover time for the subsoil (815 years, Figure 7f) than for the topsoil (168 years, Figure 7d). However, the estimated amounts are one order of magnitude larger than estimated based on ISAM-0D (Figure 7a). Over 1860-2000s, ISAM-1D estimated changes of SOC turnover time is about 5 times higher than ISAM-0D (Figures 8a and 8b), suggesting that ISAM-1D estimates lesser sequestration of SOC (or higher decay rate) in comparison to ISAM-0D. While ISAM-0D estimated Δ 14 Cbased 1 m SOC turnover time increases over the permafrost region (Figure 8a) for topsoil and subsoil (Figures 8c and 8e) over the period 1860-2000, ISAM-1D estimated Δ 14 C-based 1 m SOC turnover time increases for the tundra region and decreases for boreal forests (Figures 8a and 8b), particularly for Siberia and northern Canada boreal regions. In addition, ISAM-1D SOC turnover time (Figure 8d) for topsoil decreases for most of the regions, implicating a faster decomposition for the near future, whereas SOC turnover time (Figure 8f) increases for subsoil in the tundra and most of the boreal region, suggesting the near future the subsoil would sequester more SOC.

Effects of Carbon-Concentration and Carbon-Climate Feedbacks on Carbon Socks and Fluxes
The carbon-concentration feedback leads to a higher SOC over the time period 1860-2000s. For example, ISAM estimates higher SOC stocks in the southern Siberia boreal forest region due to a strong CO 2 fertilization effect (Figures 9a and 9c). Some of this sequestration amount is respired to the atmosphere due to carbon-climate feedback as shown by the increase of decomposition caused by higher soil temperature

10.1029/2020GB006585
Global Biogeochemical Cycles ( Figure S4). However, the estimated effects of these two feedbacks based on both versions of ISAM are quite distinct. For example, ISAM-1D (ISAM-0D) estimates higher (lower) SOC sequestration amount for topsoil than subsoil (Figures 9c-9f), because ISAM-1D (ISAM-0D) estimated effect of carbon-concentration feedback leads to a larger increase of SOC sequestrations in topsoil layers (Figures 9e, 9g, 9i, and 9k). ISAM-1D based respired soil carbon due to carbon-climate feedback for topsoil and subsoil is about the same amount and have similar regional patterns (Figures 9f,  9h, 9j, and 9l), both the topsoil and subsoil shows more SOC losses from boreal forest regions in Canada and the western and southern Siberia and from the tundra region in the eastern Siberia. In contrast, ISAM-0D estimated losses in topsoil and subsoil layers are only in tundra regions.
ISAM-1D estimates a gain of 2.3% SOC (7.3 Pg C) due to carbon-concentration feedback, whereas ISAM-0D estimates 6.4 Pg C, almost 2 times more (4.1%) carbon-concentration feedback effect than ISAM-1D. At the same time, ISAM-1D estimates a loss of SOC of 1.1% (3.6 Pg C) and ISAM-0D 1.0% (1.5 Pg C) due to carbon-climate feedback. Specific for ISAM-1D, the higher litter input caused by a warmer climate is preserved through organo-mineral interaction. In addition, the warmer climate triggers a more frequent freeze-thaw cycle and accelerates the vertical transport of SOC via the cryoturbation process. However, the part of the accumulated amount of SOC in the subsoil is lost due to a higher decomposition rate, as a result of warmer temperatures (Figure 9l). The SOC in the topsoil is lost in one part due to SOC decomposition, but in a lesser amount than in subsoil, and in another part due to the transport of carbon to the subsoil (Figure 9h).

Uncertainties in ISAM-1D Estimated SOC Stock and SOC Turnover Time Owing to Calibrated Parameters
Observation and model estimation for SOC stock and turnover time show good agreement within the uncertainty range, suggesting that they are statistically indistinguishable. The overall estimated uncertainty in ISAM-1D SOC for the topsoil ranges 64.2% to 46.4% and for subsoil ranges −78.1% to 132.5% ( Figure S5a). The parameter s, describing the strength of organo-mineral interaction, and the root distribution parameter (bb) introduce the maximum positive (46%) and negative (−64%) uncertainties in estimated topsoil stock; and the slope parameter for cryoturbation rate (γ) and the maximum depth for cryoturbation (z maxcryo ) introduce the maximum positive (133%) and negative (−78%) uncertainties in estimated subsoil stock.
The uncertainty in ISAM-1D SOC turnover time for topsoil and subsoil range −44% to 69.5% and −44.2% to 93.2% ( Figure S5b), respectively. Since SOC turnover time is the indicator of the SOC sequestration ability under the changing environmental conditions, the parameter uncertainty analysis suggests that the major parameters that introduce maximum uncertainty in ISAM-1D estimated topsoil and subsoil SOC turnover time and SOC stock are the same, with the exception of γ and D. The γ is more important for determining SOC, especially subsoil SOC, while D is more important for determining SOC turnover time.
The uncertainties in the model parameter values can be reduced by improving the cryoturbation process and its rates, particularly in the deeper soil layers. First, Miller's frost heaving model for cryoturbation implemented in ISAM-1D is based on the "rigid-ice" (O'Neill & Miller, 1985) assumption; while this model is able to capture the SOC turnover at different soil depth and the frost heave rates at the annual scale, its estimation of seasonal to diurnal variation under extreme climate change conditions can be further improved by accounting for detailed parameterization of ice lens dynamics, which can improve the cryoturbation rates. Next, the implementation of vertical soil dissolved organic carbon (DOC) flux, which currently is not accounted for in ISAM-1D, could further improve the parameterization of organo-mineral interaction and hence cryoturbation rates, particularly in the deeper soil layers (Riley et al., 2014). Although DOC is low in permafrost regions in the current environment, the vertical movement of DOC through soil water flow might become substantial during the thawing of permafrost in the future. We plan to account for this effect in our future modeling efforts based on recent efforts to model DOC in permafrost regions (e.g., Bowring et al., 2019).

Modeled Steady-State SOC Sock and SOC Turnover Time
After incorporating a number of processes to resolve the vertical soil heterogeneity, including the spatially variable cryoturbation rates, the vertically resolved litter input, and environmental modifiers, ISAM-1D improved the estimation of spatial heterogeneity for both permafrost SOC stocks and turnover time for topsoil and subsoil pools. The spatial heterogeneity of the cryoturbation rate in ISAM-1D, which is calculated through linking the frost heave process with the freezing and thawing of the soil ice, is able to capture more accurately the redistribution of SOC stocks in the permafrost region soils affected by environmental factors. For example, ISAM-1D (Figure 10a) estimated the highest cryoturbation rate to be found in the regions caused by more active freezing and thawings, such as the northern tundra region and the Hudson Bay. ISAM-1D also estimated the change of cryoturbation rates over the 1860-2000s (Figure 10b), showing a larger region with an increase of cryoturbation rate than the region with a decreasing rate, resulting in a faster

10.1029/2020GB006585
Global Biogeochemical Cycles movement of topsoil SOC stocks into the subsoils. The increase of cryoturbation rate is linked to the increase of soil temperature ( Figure S4), which activates stronger freezing and thawing activity of soil ice and hence increases the frost heave strength and the estimated cryoturbation rate.
ISAM-1D estimated Δ 14 C-based SOC turnover time is also consistent with the published radiocarbon dates of buried SOC stocks (Bockheim, 2007), which is longer than ISAM-0D, resulting in an increase of steadystate SOC. While SOC turnover time is longer, ISAM-1D calibrated CUE for labile pools is lower than ISAM-0D, resulting in a reduced amount of SOM transferred to the resistant pools in topsoil. Therefore, more SOC is decomposed and then respired to the atmosphere from the labile pools of topsoil layers, resulting in a transfer of lesser amount to the deeper soil layers by the cryoturbation process. While moving downward along the vertical soil layers, decomposition rates of SOC are reduced due to organo-mineral interaction accounted for in ISAM-1D, resulting in an increase of SOC turnover time and thus stock in subsoils ( Figure 8). Since the ISAM-1D considers the vertically transport of litter input, it suffers the same fate and hence further increases the SOC stocks and turnover time.
In order to tease out the effect of calibrated CUE for labial pools from other vertically resolved processes incorporated in the ISAM-1D, we performed another model simulation using ISAM-1D with the original CUE value (i.e., the value used in ISAM-0D) over the 197 soil sample sites . The comparison of the two versions of the ISAM-1D, one with calibrated CUE and another with original CUE, shows that the model estimated SOC sink for the calibrated case is reduced by about 20% ( Figure S9).
In contrast to ISAM-0D, which assumed fixed values of the modifier for the SOC decomposition calculation, ISAM-1D calculations with vertically resolved biogeophysical modifiers for soil decompositions show the improvement in the model estimated decomposition and Rh compared to Soil Respiration Database (SRDB) observation data (Bond-Lamberty & Thomson, 2018). This is because ISAM-1D SOC decomposition rates were higher in the topsoil than in ISAM-0D due to the warmer daytime topsoil average teamperature than the surface air temperature, which was used in the ISAM-0D. In the case of soil moisture modifier, the ISAM-0D applied only the first 30 cm of soil moisture, thus lead to a dry bias for litter decomposition, while ISAM-1D applied the vertically resolved moisture at each layer, which was wetter. In addition, the biogeochemical modifiers, including oxygen and nitrogen availability, also enhance the ISAM-1D steady-state SOC stocks and turnover time. For example, the oxygen modifier in ISAM-1D under the steady-state conditions shows a mean oxygen modifier of 0.8, and most of the oxygen-limited region is found in the northern permafrost region ( Figure S6a) due to the shallower water table, which limits the oxygen diffusion into the soil column. In the case of nitrogen modifier, ISAM-1D shows a stronger N limitation in the northern region compared to the southern boreal region ( Figure S7a).

Modeled Carbon Cycle Feedback Effects on Soil Carbon Sequestration Potential Over Historical Time
ISAM-1D estimates longer SOC turnover time, suggesting an increase in SOC, a positive response. However, this response changes under carbon-concentration and carbon-climate feedbacks over the historical period (1860s-2000s). ISAM-1D estimates a weaker carbon-concentration and stronger carbon-climate feedbacks compared to ISAM-0D. Thus, overall lesser carbon sequestration was estimated based on ISAM-1D as reflected by the changes in SOC (Table 3). This is because SOC losses from labile pools of the topsoil layers are about 10 times more compared to the subsoil ( Figures S8a and S8b), whereas ISAM-0D sequestration rates were higher, because more SOC was transferred from labile pools to the resistant pools as discussed next.
A comparison of ISAM-0D and ISAM-1D results for carbon dynamics of the labile pool ( Figures S8a and S8b) shows shorter SOC turnover time in ISAM-1D than in ISAM-0D for the labile pool, suggesting that the labile pool in ISAM-1D case is receiving additional carbon input, but respiring back to the atmosphere at much faster rates as discussed above. This results in a less carbon remain in labile pools for further transform or transport to the subsoil pools with a longer SOC turnover time and higher SOC sequester potential (Figures S8c and S8d). In contrast, resistant pools show the opposite response. Overall, ISAM-1D shows a 1.2% gain of SOC, which compared to the estimation of 3.1% gain of SOCS from ISAM-0D (Table 3), confirming a large decrease of carbon sequestration ability in comparison to ISAM-0D. In addition, the changes in SOC stock in ISAM-1D were also impacted by oxygen availability. Compared to the 1860s, oxygen availability is increased by the 2000s (Figure S6b), resulting in an increase of SOC decomposition and the loss of more SOC stock. Thus, further weakens the carbon-concentration feedback effect.
In regards to carbon-climate feedback, warmer conditions accelerate the decomposition of SOC in the topsoil layer. In contrast, warmer condition enhances the soil freezing-thawing frequency, resulting in an enhancement of the frost heave rate and thus the cryoturbation rate (Equations S4-S6). A higher cryoturbation rate transports more SOC to the deeper soil layers where it is preserved due to stronger organo-mineral interaction. However, ISAM-1D shows almost similar warming for the topsoil and subsoil layers ( Figure S10), but the studies also show that higher temperatures are more sensitivity to the resistant pool (Equation S12) compared to the labile pool (Knorr et al., 2005;Schmidt et al., 2011), which we incorporated in ISAM by using a variant of Q 10 formulations for labial and resistant pools (Equation S11 for labial pool and Equation S12 for a resistant pool). Since the resistant pools in subsoils have a relatively higher amount of SOC, therefore, there is an increased amount of carbon loss from the subsoil layers. In addition, higher oxygen availability under warmer temperature further accelerates carbon decomposition, which means strengthening the carbon-climate feedback.

Future Perspective of the Work
Despite the advances made in ISAM-1D, the model performance perhaps can be further improved by improving the parameterization of modeling cryoturbation rates, and hence SOC stock and turnover time, particularly in the deeper soil layers. First, Miller's frost heaving model for cryoturbation implemented in ISAM-1D is based on the "rigid-ice" (O'Neill & Miller, 1985) assumption; while this model is able to capture the SOC turnover time at different soil depth and the frost heave rates at the annual scale, its estimation of seasonal to diurnal variation under extreme climate change conditions can be further improved by accounting for detailed parameterization of ice lens dynamics. Second, the implementation of vertical soil DOC flux, which currently is not accounted for in ISAM-1D, could further improve the parameterization of organo-mineral interaction and hence cryoturbation rates, particularly in the deeper soil layers (Riley et al., 2014). Although DOC is low in permafrost regions, the vertical movement of DOC through soil water flow might become important during the thawing of permafrost in the near future and would be accounted for in our future modeling efforts. A recommendation for further reducing this uncertainty for large scale TEMs is to explicitly account for the spatial heterogeneity of soil mineralogy (Torn et al., 1997).

Overall Conclusion
Overall, ISAM-1D model results for the permafrost region show substantial improvements in modeling both SOC stocks and turnover time compared to previous study results (He et al., 2016). Also, ISAM-1D estimated 10.1029/2020GB006585

Global Biogeochemical Cycles
SOC profiles are consistent with observations and thus can reproduce the different responses of topsoil and subsoil to the change of the environmental conditions. The revised ISAM-1D model structure, added processes, and effects of environmental factors improve the model performance for carbon-concentration and carbon-climate feedback effects. The model estimates lesser carbon sequestration in the permafrost region soils over the historical time period (1860-2000s). While the effect of these to feedbacks is small over the historical time, we expect ISAM-1D to calculate much less soil carbon sequestration rates under the future scenarios with higher CO 2 and higher temperature changes.
Overall, ISAM-1D framework with the improved representation of SOC dynamics and turnover time is a better tool to quantify the future permafrost SOC dynamics, greenhouse gas release, and the permafrost carbon cycle-climate feedback compared to the current ESMs.

Data Availability Statement
ISAM simulated decadal mean permafrost SOC stocks and turnover time for the 1860s and 2000s can be downloaded from the ISAM website (http://climate.atmos.uiuc.edu/ShuEtAl_Permafrost/). All the results are stored in NetCDF files. In addition, ISAM model code is available upon request from the corresponding author (email: jain1@illinois.edu). The ISAM-0D and ISAM-1D results discussed in this paper and listed in the Data Access section are available for download from the ISAM website (http://climate.atmos.uiuc.edu/ ShuEtAl_Permafrost/).