<figure><img src="image_1.png" /></figure>

HAL
open science

Representing high-latitude deep carbon in the
pre-industrial state of the ORCHIDEE-MICT land
surface model (r8704)

Yi Xi, Philippe Ciais, Dan Zhu, Chunjing Qiu, Yuan Zhang, Shushi Peng, Gustaf Hugelius, Simon P K Bowring, Daniel S Goll, Ying-Ping Wang

<figure><img src="image_2.png" /></figure>

To cite this version:

Yi Xi, Philippe Ciais, Dan Zhu, Chunjing Qiu, Yuan Zhang, et al.. Representing high-latitude deep carbon in the pre-industrial state of the ORCHIDEE-MICT land surface model (r8704). Geoscientific Model Development, 2025, 18 (18), pp.6043 - 6062. 10.5194/gmd-18-6043-2025. hal-05282584

HAL Id: hal-05282584

https://hal.science/hal-05282584v1

Submitted on 25 Sep 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

# Representing high-latitude deep carbon in the pre-industrial state of the ORCHIDEE-MICT land surface model (r8704)

Yi Xi<fnref n="1" />, Philippe Ciais<fnref n="1" />, Dan Zhu<fnref n="2" />, Chunjing Qiu<fnref n="3" />,<fnref n="4" />, Yuan Zhang<fnref n="5" />, Shushi Peng<fnref n="2" />, Gustaf Hugelius<fnref n="6" />,<fnref n="7" />,
Simon P. K. Bowring<fnref n="1" />, Daniel S. Goll<fnref n="1" />, and Ying-Ping Wang<fnref n="8" />

**Correspondence:** Yi Xi (yi.xi@lsce.ipsl.fr)

Received: 5 November 2024 – Discussion started: 23 January 2025

Revised: 23 June 2025 – Accepted: 30 June 2025 – Published: 16 September 2025

**Abstract.** Field measurements, after extrapolation, suggest that deep Yedoma deposits (ice-rich, organic-rich permafrost, formed during the late Pleistocene) and peatlands (formed mostly during the Holocene) account for about 600 Pg C of soil carbon storage. Incorporating this old, deep, cold carbon into land surface models (LSMs) is crucial for accurately quantifying soil carbon responses to future warming. However, it remains underrepresented or absent in current LSMs, which typically include a passive soil carbon pool (a conceptual soil carbon pool with the longest turnover time) to represent all “old carbon” and lack the vertical accumulation processes that deposited deep carbon in the layers of peatlands and Yedoma deposits. In this study, we propose a new, more realistic protocol for simulating deep and cold carbon accumulation in the northern high latitudes (30–90° N), using the ORCHIDEE-MICT (Organizing Carbon and Hydrology in Dynamic Ecosystems–aMeliorated Interactions between Carbon and Temperature) model. This is achieved by (1) integrating deep carbon from Yedoma deposits whose formation is calculated using Last Glacial Maximum climate conditions and (2) prescribing the inception time and location of northern peatlands during the Holocene using spatially explicit

data on peat age. Our results show an additional 157 Pg C in present-day Yedoma deposits, as well as a shallower peat carbon depth (by 1–5 m) and a smaller passive soil carbon pool (by 35 Pg C, 43 %) in northern peatlands, compared to the old protocol that ignored Yedoma deposits and applied a uniform, long-duration (13 500 years) peat carbon accumulation across all peatlands. As a result, the total organic carbon stock across the Northern Hemisphere (30–90° N) simulated by the new protocol is 2028 Pg C, which is 226 Pg C higher than the previous estimate. Despite the significant challenge of simulating deep carbon with ORCHIDEE-MICT, the improvements in the representation of carbon accumulation from this study provide a model version to predict deep carbon evolution during the last glacial–deglacial transition and its response to future warming. The methodology implemented for deep carbon initialization in permafrost and cold regions in ORCHIDEE-MICT is also applicable to other LSMs.

<footnote n="1">¹Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ,
Université Paris-Saclay, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="2">²Sino-French Institute for Earth System Science, Institute of Carbon Neutrality, College of Urban and Environmental Sciences, Peking University, Beijing, China</footnote>

<footnote n="3">³Research Center for Global Change and Ecological Forecasting, School of Ecological and Environmental Sciences, East China Normal University, Shanghai, China</footnote>

<footnote n="4">⁴Institute of Eco-Chongming, East China Normal University, Shanghai, China</footnote>

<footnote n="5">⁵Key Laboratory of Alpine Ecology, Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing, China</footnote>

<footnote n="6">⁶Department of Physical Geography, Stockholm University, 10691 Stockholm, Sweden</footnote>

<footnote n="7">⁷Bolin Centre for Climate Research, Stockholm University, 10691 Stockholm, Sweden</footnote>

<footnote n="8">⁸CSIRO Environment, Private Bag 10, Clayton South, VIC 3169, Australia</footnote>
# 1 Introduction

Soil organic carbon (SOC) is a major component of the global carbon cycle, as the size and composition of this reservoir directly control the amount of carbon dioxide ($CO_2$) emitted from terrestrial ecosystems via heterotrophic respiration (Rh) (Jackson et al., 2017). In the absence of other disturbances, the imbalance between the responses of Rh and net primary production (NPP) to climate warming determines whether terrestrial ecosystems will become a carbon source or a carbon sink in the future (Crowther et al., 2016; Lu et al., 2013; Walker et al., 2018). Earth system models (ESMs) employed in the recent IPCC Sixth Assessment Report (AR6) predict that terrestrial ecosystems will continue to act as a net carbon sink until 2100 due to the relatively dominant response of NPP to warming but may shift to a net source after 2100 as SOC respiration exceeds vegetative carbon uptake (IPCC, 2021). Variations in SOC density also indirectly influence global carbon cycling processes by altering soil hydrology (e.g., soil porosity) and soil heat transfer (e.g., thermal diffusivity, especially in cold regions) (Zhu et al., 2019). Collectively, these findings suggest that accurately representing and estimating global SOC stocks is paramount to better understanding the overall response of terrestrial ecosystems to future climate warming.

Widely used global SOC stock gridded datasets such as HWSD (the Harmonized World Soil Database; Hiederer and Köchy, 2012), WISE30sec (the World Inventory of Soil property Estimates project; Batjes, 2016), and SoilGrids250m (global gridded soil information; Hengl et al., 2017) are currently based on site-level soil profile measurements, which are extrapolated to regional scales through their combination with digitized data (SoilGrids250m). Due to the challenges posed by deep soil coring and sampling, these products usually only produce global SOC estimates to relatively shallow depths, e.g., 1417 Pg C for 0–1 m soil depth from HWSD (Hiederer and Köchy, 2012) and 2060 ± 215 Pg C for 0–2 m soil depth from WISE30sec (Batjes, 2016). However, high-latitude soils underlain by permafrost and peat typically exhibit vertical profiles significantly deeper than 2 m. Accounting for SOC stocks at the 2–3 m depth range would add roughly ~470 Pg C (+23 %) to the global soil stocks, including 175 (142–220) Pg C in permafrost deep soils (Mishra et al., 2021), 199 Pg C in non-permafrost deep mineral soils (Jobbagy and Jackson, 2000), and 92 Pg C in northern peatlands outside permafrost regions (Jackson et al., 2017).

Going deeper than 3 m, a substantial storage of SOC can still be found. In particular, in Alaska, northern Siberia, and northwest Canada, there is an area of $0.48 \times 10^6$ km² covered by icy, organic-rich, silt-dominated deposits, known as Yedoma (Strauss et al., 2017). These deposits were deposited to an average depth of 19–25 m by windblown particles over ice-free regions during the Pleistocene (Anthony et al., 2014; Strauss et al., 2013). Yedoma likely covered a vast area during the late Pleistocene, though only ~30%

has survived thawing and warming during the transition to the Holocene (Strauss et al., 2017). Empirical extrapolations based on field measurements suggest an SOC storage of 115 (83–269) Pg C (Strauss et al., 2017) in present-day Yedoma deposits. Peatlands are another reservoir of deep organic carbon in the northern high latitudes. Most formed during the Holocene, following the drainage of thermokarst lakes fed by the thawed ice from Yedoma deposits (Anthony et al., 2014). Peat carbon can accumulate to depths exceeding 10 m, due to anaerobic constraints on decomposition (Hugelius et al., 2020). At present, northern peatlands cover an area of $3.7 (±0.5) \times 10^6$ to $4.0 \times 10^6$ km² and account for 415 (±150) to 547 (±74) Pg C (Hugelius et al., 2020; Yu et al., 2010). Despite large uncertainties in these estimates, which result from inadequate field measurements and rough estimation methods, the substantial SOC storage (~600 Pg C) in Yedoma deposits and northern peatland soils highlights the importance of incorporating deep carbon in the accounting of global soil carbon stocks.

Land surface models (LSMs) are valuable tools for assessing the response of terrestrial ecosystems to climate warming. However, the simulated response by most LSMs is subject to biases in simulated SOC due to limitations in existing model structures and missing processes, e.g., deep carbon from Yedoma deposits and peatland soils (Bradford et al., 2016; Koven et al., 2015; Niu et al., 2024; Shi et al., 2024; Wieder et al., 2013). Recent efforts to improve model structures and develop new processes have emerged in some model groups, such as in ORCHIDEE-MICT (Organizing Carbon and Hydrology in Dynamic Ecosystems–aMeliorated Interactions between Carbon and Temperature) – a version of the land surface component of the IPSL global climate model focused on resolving high-latitude processes (Guimberteau et al., 2018). Key developments for SOC in this model include the representation of vertically resolved organic carbon inputs from litter and depth-dependent decomposition processes to simulate the vertical discretization of soil carbon (Koven et al., 2009), the incorporation of effects of SOC on soil thermal properties to improve the simulated SOC in the top meter (Koven et al., 2009), the representation of vertical mixing of SOC via cryoturbation and bioturbation to vertically distribute SOC to a depth of ~3 m (Koven et al., 2009), the representation of sedimentation processes to simulate SOC below 3 m over Yedoma deposits (Zhu et al., 2016), and the incorporation of peat-specific hydrology and carbon accumulation processes to simulate the dynamic peat carbon accumulation (Qiu et al., 2018, 2019). However, these advancements, developed in separate branches of ORCHIDEE-MICT, have yet to be integrated together for a realistic, process-based simulation of vertical and horizontal SOC distributions across the northern high latitudes.

Taking ORCHIDEE-MICT (hereafter MICT) as a starting point, in this study we present a new protocol for the accumulation of both shallow and deep soil organic carbon in the
northern high latitudes (30–90°N) by (1) merging the developed processes for SOC from various branches of MICT, including the carbon accumulation for non-peat vegetation types (Guimberteau et al., 2018), peat carbon accumulation (Qiu et al., 2018, 2019), sedimentation processes for Yedoma deposits (Zhu et al., 2016), and tile-specific energy budgets (Xi et al., 2024), and (2) representing the observed history of peat carbon evolution by integrating peat inception time according to radiocarbon dates measured from peat cores (Hugelius et al., 2020; Loisel et al., 2017). The simulated results with the new protocol are compared with those simulated with the old carbon accumulation protocol to evaluate if the new protocol better reproduces observed SOC distributions.

## 2 Data and methods

To implement the model development, we used a branch version of ORCHIDEE-MICT, namely ORCHIDEE-MICT-teb (r8205), from Xi et al. (2024) as our base model. This version has PFT-specific (plant-functional-type-specific) energy budgets, soil thermic properties, and their interactions with carbon and hydrological processes and has been comprehensively evaluated for its performance in simulating carbon, energy, and hydrological processes (Xi et al., 2024). To simulate Yedoma carbon accumulation, we incorporated sedimentation processes from Zhu et al. (2016), adapting them to be compatible with the PFT-specific framework of the new model. The model version used by Zhu et al. (2016) relied on grid-cell-averaged energy budgets, which meant that Yedoma carbon was mixed with coexisting PFTs (bare soil, tree, grass, crop) to regulate the average energy budget. Given that Yedoma deposits typically have much higher carbon densities than general soils, accounting for their distinct carbon-energy-water feedbacks such as the thermal insulation effect of SOC (Koven et al., 2009; Zhu et al., 2019) is essential for realistically simulating both their historical development and their future evolution. Regarding peatland processes, the base version already included peatland carbon accumulation. We thus did not develop new peatland processes, but we revised the way peatland initiation is simulated during the Holocene. Instead of assuming a homogeneous peatland age across northern peatlands starting in the early Holocene, we employed spatially explicit peatland age maps based on observational data to define the onset of peat development at each grid cell. This was achieved by prescribing peatland cover maps at different periods during the Holocene, which the model uses as a forcing to simulate realistic peatland expansion over time. As a result, the improved model allows for the independent simulation of Yedoma carbon, peatland carbon, and conventional soil carbon dynamics within a single grid cell.

In Sect. 2.1, we describe the soil carbon models in detail, including the general soil carbon processes and the peat-

land carbon accumulation scheme in the base version (Xi et al., 2024), as well as the Yedoma carbon accumulation processes from Zhu et al. (2016). Section 2.2 then presents the setup of the carbon accumulation simulations, using both the base version (without Yedoma carbon accumulation and with uniform peatland inception timing) and the improved version (with Yedoma carbon accumulation and spatially varying peatland inception).

### 2.1 Soil carbon model description

The soil carbon model in the original MICT (v8.4.1) (Guimberteau et al., 2018) is based on the CENTURY model (Paton et al., 1988), with two litter carbon pools (metabolic and structural) and three soil carbon pools (active, slow, and passive) defined by their turnover rates. The turnover time at 5 °C without moisture limitation for the five carbon pools is 0.37, 1.4, 0.84, 31, and 1363 years, respectively. The vertical transfer and accumulation of soil carbon are driven by three main processes: root-density-dependent organic carbon inputs, depth-dependent decomposition of soil carbon regulated by vertically layered soil temperature and soil moisture (considering pore ice), and carbon diffusion via bioturbation by animal and plant activity and via cryoturbation in permafrost soils. The assumed soil depth for soil hydrology is set to 2 m (11 layers), and the total depth for soil thermal/-carbon processes is 38 m (32 layers). The hydrological and thermal/carbon processes share the first 11 layers, while the thermal layers below use the soil moisture from the 11th hydrology layer. Limited by the root depth and active layer thickness, however, the deepest soil carbon simulated by the model is ~3 m. Cryoturbation can bury carbon below the active layer in the model, but this buried carbon merely exceeds a depth of 3 m. More details on soil carbon processes in MICT can be found in Krinner et al. (2005), Koven et al. (2009), Zhu et al. (2016), Guimberteau et al. (2018), and Huang et al. (2018).

Qiu et al. (2018, 2019) introduced a new soil tile to host peatland into MICT. Due to the unique highly saturated, oxygen-deficient, and nutrient-limited soil environments in peatlands, this soil tile is covered by specific vegetation – a $C_3$ grass plant functional type (PFT) with a shallower rooting depth – and is prescribed with specific soil hydrology properties, including hydraulic conductivity and diffusivity. Further, the peat soil tile has no drainage flux at the base of the soil and accepts lateral surface water input from other biomes within the same grid cell (Qiu et al., 2018). The waterlogged soil greatly reduces decomposition rates and leads to the accumulation of deep organic matter in peatlands. Instead of modeling the upward accumulation of peat, Qiu et al. (2019) simulated a downward “sedimentation” of peat carbon in the model, inspired by Jafarov and Schaefer (2016). When the simulated SOC at soil layer $l$ is larger than a prescribed fraction (0.95) of observation-based carbon density ($C_{obs,l}$, kg m$^{-2}$), a prescribed fraction (0.05) of SOC will be
moved down to the layer below (l + 1), resulting in accumulation of carbon in peat soils over time. The Cobs,l threshold is calculated as

$$C_{\text{obs},l} = BD_l \times \alpha_{c,l} \times \Delta Z_l, \quad (1)$$

where BDl (kg m⁻³) is the soil bulk density at model layer l (estimated with observation data), αc,l (wt%) is the mass fraction of carbon in the soil (derived from observation data), and ΔZl (m) is the thickness of soil layer l. Compared with northern peatland sites and gridded simulations over the Northern Hemisphere (> 30°N), the model generally reproduces the measured SOC density, peat depth from peat cores, and total northern peatland area (3.9 × 10⁶ km² against 3.4 × 10⁶ to 4.0 × 10⁶ km²) and total northern peat carbon stocks (463 Pg C against 270–540 Pg C) from previous studies (Qiu et al., 2019).

The simulation of the sedimentation processes forming Yedoma deposits was introduced by Zhu et al. (2016). Similar to peat carbon accumulation, the basic idea for modeling the upward accumulation of carbon when Yedoma deposits form is to simulate an equivalent downward transport of organic carbon (OC) in the model by employing a vertical advection equation:

$$\frac{\partial C_i(z, t)}{\partial t} + u(t) \frac{\partial C_i(z, t)}{\partial z} = f_i(z, t) - g_i(z, t) \times C_i(z, t), \quad (2)$$

where Ci(z, t) (g m⁻³) is the OC concentration of the litter or soil carbon pool i at depth z and time t, u(t) (m d⁻¹) is the downward movement rate of OC, fi(z, t) (g m⁻³ d⁻¹) is the OC input to pool i, and gi(z, t) is the decomposition rate of pool i. The downward rate of OC, i.e., the sedimentation rate, can be inferred from the height–age slope of Yedoma deposits based on site measurements. Since the height of Yedoma deposits is at landscape scale, the inferred sedimentation rate includes ice wedges (~ 50% in volume). In this way, although the representation of ice wedges is absent in MICT, the volume of ice wedges is implicitly included in the model. To reach the mean deposit thickness of 19.4 m suggested by Strauss et al. (2013), four typical sedimentation rates (0.6, 0.8, 1.0, and 1.2 mm yr⁻¹) observed in Yedoma sites with corresponding accumulation duration (32 000, 24 000, 19 000, and 16 000 years; hereafter kyr is used to denote thousand years) were assessed to simulate the present-day Yedoma deposits. The soil carbon stock (125–145 Pg C) in present-day Yedoma region (0.42 × 10⁶ km²) simulated by the model is comparable to the estimate by Strauss et al. (2013) (83 +61/-57 Pg C).

## 2.2 Simulation setup

Conventional LSMs employ a “spin-up” simulation to equilibrate the biogeochemical cycles before calculating perturbations. A perfect spin-up with a perfect model would accumulate SOC stocks to exactly match the observations. For

instance, in the TRENDY (Trends in the land carbon cycle) project, the spin-up simulation is run with recycled climate from the first two decades of the 20th century (1901–1920), atmospheric CO₂ concentration from the year 1700 (276.59 ppm), and a land cover map fixed to the year 1700 (Sitch et al., 2015). This classical spin-up procedure is a strong simplification, since in reality, soils were not in a steady state under pre-industrial climatic conditions. To address this problem, we designed a new spin-up protocol, aiming at reproducing the formation of deep carbon in Yedoma deposits and peatland soils since the late Pleistocene (see a schematic representation in Fig. 1). Below, we describe both the classical (or old) spin-up simulation with the base version (Xi et al., 2024), as used in the TRENDY project, and the new spin-up protocol with the improved MICT version, which includes the representation of deep carbon accumulation in Yedoma deposits and observed history of peat carbon evolution.

### 2.2.1 Old spin-up simulation

Table 1 shows the protocol of the old spin-up simulation for carbon accumulation with the starting MICT version. It includes three stages: (1) 100 years with the full ORCHIDEE land surface model (FullIO) run to initiate the energy cycle, hydrological processes, and vegetation growth in the model (the upper grid cell of the old spin-up in Fig. 1); (2) 13 500 years with a standalone soil carbon module (SubC) run to accelerate carbon accumulation (using fixed monthly litter input, soil moisture, and soil thermal conditions from the preceding FullIO simulation); and (3) 100 years with the FullIO run to bring energy, hydrology, and carbon back into equilibrium (the lower grid cell of the old spin-up in Fig. 1). All three stages of the old spin-up simulation use the same forcing datasets, including climate data recycling 1901–1920 from CRU-JRA v2.3 (the version used in the Global Carbon Budget 2021) (Friedlingstein et al., 2022), the fixed CO₂ level in 1700 (276.59 ppm), and the fixed land cover map generated by combining the 1700 land cover map from TRENDY for 15 conventional PFTs (bare soil, 8 tree PFTs, 4 grass PFTs, and 2 crop PFTs; see more details in Table S1) (Lurton et al., 2020) and the peat map from Hugelius et al. (2020) for the peat C₃ grass PFT. All the simulations were run for the Northern Hemisphere (30–90°N) at a 2° × 2° spatial resolution. Thus, the old spin-up simulation simply runs the model for a very long period (13 500 years) to accumulate enough soil carbon. However, it does not include Yedoma deposits, and although the starting MICT version used in the old spin-up has merged the peat carbon accumulation, it uses a uniform duration (13 500 years) of peat carbon accumulation for all peatlands and therefore misses the evolution of peatland extent as in reality (Fig. 1).
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Schematic representation of soil carbon accumulation within one grid cell across different simulation stages in the old and new spin-up simulations. The old spin-up uses the ORCHIDEE-MICT-teb (MICT-teb) model version and includes 16 plant functional types: 15 conventional PFTs (bare soil, 8 tree PFTs, 4 grass PFTs, and 2 crop PFTs) and 1 peatland PFT (C₃ grass) within a grid cell (Table S1). In the old spin-up, the cover fraction of peatland is fixed throughout the simulation, and all peatlands are assumed to have formed simultaneously, 13 500 years ago. By contrast, the new spin-up merges ORCHIDEE-MICT-Yedoma (MICT-Yedoma) with MICT-teb, enabling the model to simulate Yedoma carbon accumulation during the late Pleistocene and subsequent carbon accumulation in peatland and conventional soils during the Holocene. The new spin-up includes the same 15 conventional PFTs, 1 peatland PFT, and an additional Yedoma PFT (C₃ grass) within a grid cell (Table S1). Peatland cover varies dynamically across nine Holocene periods (shown in the peat PFT domain), based on prescribed maps constructed from peat basal age observations and present-day peatland distribution.</figcaption></figure>

<table><caption>**Table 1.** Protocol for the old spin-up simulation. FullO indicates the full ORCHIDEE-MICT and SubC indicates the standalone soil carbon module.</caption><thead><tr><td>Simulation</td><td>Duration</td><td>Climate</td><td>CO<sub>2</sub></td><td>Land cover map</td></tr></thead><tbody><tr><td>FullO</td><td>100 years</td><td></td><td></td><td></td></tr><tr><td>SubC</td><td>13 500 years</td><td>CRU-JRA 1901–1920 cycling</td><td>276.59 ppm</td><td>TRENDY 1700 + peat incepted before 1–0 ka</td></tr><tr><td>FullO</td><td>100 years</td><td></td><td></td><td></td></tr></tbody></table>

## 2.2.2 New spin-up simulation

Since radiocarbon dating of Yedoma deposits and peat cores provides their age of deposition, we are no longer entirely blind to the history of soil carbon evolution. In the new spin-up simulation, we explicitly represented the timing of inception and the duration of accumulation for Yedoma deposits and peat carbon (Fig. 1). Table 2 summarizes the new spin-up simulation.

### Yedoma deposits

In the real world, Yedoma deposits accumulated over tens of thousands of years due to the cold climatic and environmental conditions in the late Pleistocene glacial period, with about 70% of them undergoing thermokarst degradation or fluvial erosion during the warmer Holocene (Grosse et al., 2013; Strauss et al., 2017). In the new spin-up, after merging the Yedoma sedimentation processes developed by
<table><caption>**Table 2.** Protocol for the new spin-up simulation. FullO indicates the full ORCHIDEE-MICT land surface model, and SubC indicates the standalone soil carbon module.</caption><thead><tr><th>Period</th><th>Simulation</th><th>Duration</th><th>Climate</th><th>CO<sub>2</sub></th><th>Land cover map</th></tr></thead><tbody><tr><td colspan="6">Yedoma simulation</td></tr><tr><td rowspan="3">21–18 ka</td><td>FullO</td><td>100 years</td><td rowspan="3">100-year LGM climate cycling</td><td rowspan="3">185 ppm</td><td rowspan="3">Simulated by dynamic vegetation model in MICT</td></tr><tr><td>SubC</td><td>19 000 years</td></tr><tr><td>FullO</td><td>100 years</td></tr><tr><td colspan="6">Peat simulation</td></tr><tr><td rowspan="2">13.5–11.5 ka</td><td>FullO</td><td>100 years</td><td rowspan="2">CRU-JRA 1901–1920 cycling</td><td rowspan="2">238 ppm</td><td rowspan="2">TRENDY 1700 + peat incepted during 16–12 ka + Yedoma PFT</td></tr><tr><td>SubC</td><td>2000 years</td></tr><tr><td rowspan="2">11.5–10.5 ka</td><td>FullO</td><td>100 years</td><td rowspan="2">CRU-JRA 1901–1920 cycling</td><td rowspan="2">250 ppm</td><td rowspan="2">TRENDY 1700 + peat incepted during and before 12–11 ka + Yedoma PFT</td></tr><tr><td>SubC</td><td>1000 years</td></tr><tr><td rowspan="2">10.5–9.5 ka</td><td>FullO</td><td>100 years</td><td rowspan="2">CRU-JRA 1901–1920 cycling</td><td rowspan="2">265 ppm</td><td rowspan="2">TRENDY 1700 + peat incepted during and before 11–10 ka + Yedoma PFT</td></tr><tr><td>SubC</td><td>1000 years</td></tr><tr><td rowspan="2">9.5–8.5 ka</td><td>FullO</td><td>100 years</td><td rowspan="2">CRU-JRA 1901–1920 cycling</td><td rowspan="2">265 ppm</td><td rowspan="2">TRENDY 1700 + peat incepted during and before 10–9 ka + Yedoma PFT</td></tr><tr><td>SubC</td><td>1000 years</td></tr><tr><td rowspan="2">8.5–7.5 ka</td><td>FullO</td><td>100 years</td><td rowspan="2">CRU-JRA 1901–1920 cycling</td><td rowspan="2">260 ppm</td><td rowspan="2">TRENDY 1700 + peat incepted during and before 9–8 ka + Yedoma PFT</td></tr><tr><td>SubC</td><td>1000 years</td></tr><tr><td rowspan="2">7.5–6.5 ka</td><td>FullO</td><td>100 years</td><td rowspan="2">CRU-JRA 1901–1920 cycling</td><td rowspan="2">260 ppm</td><td rowspan="2">TRENDY 1700 + peat incepted during and before 8–7 ka + Yedoma PFT</td></tr><tr><td>SubC</td><td>1000 years</td></tr><tr><td rowspan="2">6.5–5.5 ka</td><td>FullO</td><td>100 years</td><td rowspan="2">CRU-JRA 1901–1920 cycling</td><td rowspan="2">262 ppm</td><td rowspan="2">TRENDY 1700 + peat incepted during and before 7–6 ka + Yedoma PFT</td></tr><tr><td>SubC</td><td>1000 years</td></tr><tr><td rowspan="2">5.5–4 ka</td><td>FullO</td><td>100 years</td><td rowspan="2">CRU-JRA 1901–1920 cycling</td><td rowspan="2">265 ppm</td><td rowspan="2">TRENDY 1700 + peat incepted during and before 6–5 ka + Yedoma PFT</td></tr><tr><td>SubC</td><td>1500 years</td></tr><tr><td rowspan="2">4–0 ka</td><td>FullO</td><td>100 years</td><td rowspan="2">CRU-JRA 1901–1920 cycling</td><td rowspan="2">272 ppm</td><td rowspan="2">TRENDY 1700 + peat incepted during and before 5–0 ka + Yedoma PFT</td></tr><tr><td>SubC</td><td>4000 years</td></tr><tr><td></td><td>FullO</td><td>100 years</td><td>CRU-JRA 1901–1920 cycling</td><td>276.59 ppm</td><td>TRENDY 1700 + peat incepted during and before 5–0 ka + Yedoma PFT</td></tr></tbody></table>

Zhu et al. (2016) into the starting MICT version, we used a 1 mm yr⁻¹ sedimentation rate and a 19 kyr accumulation period as in their study to simulate carbon stocks in Yedoma deposits under the Last Glacial Maximum (LGM) climate (21 000–18 000 years before present; hereafter ka is used to denote thousand years before present). A map of present-day Yedoma from Strauss et al. (2021), covering 0.48 × 10⁶ km² of Yedoma deposits (Fig. 2), was used to limit the model to predict Yedoma deposits which have survived the deglaciation. This was done because we aim to represent Yedoma carbon during the pre-industrial period in this study. We did not test other sedimentation rates and accumulation periods here because Zhu et al. (2016) found that doing so has a small impact on simulated total SOC stock. Other deposits, such as those from thermokarst lake basins, were not included in this study because they were deposited via different processes. To accelerate the carbon accumulation, we used the “FullO–SubC–FullO” strategy (Table 1) in the 19 kyr simulation. We used the LGM climate forcing data output from the Institut Pierre-Simon Laplace (IPSL) Earth system model, IPSL-CM5A-LR (Kageyama et al., 2013); the LGM CO₂ level at 185 ppm; and the LGM vegetation growth simulated by the

dynamic vegetation model in MICT due to the lack of land cover maps during that period (Zhu et al., 2015). Consistently with the old spin-up simulation, the Yedoma simulation was also run at a 2° × 2° grid resolution but only for 60–90° N to save computation cost, given that no Yedoma deposits are found south of 60° N (Fig. 2; Strauss et al., 2021). The schematic representation of Yedoma carbon accumulation is shown with the two upper grid cells for the new spin-up in Fig. 1. Further details on how the Yedoma simulation is connected with the peatland simulation are provided in the following subsection.

## Peatlands

The radiocarbon dates obtained from peat cores have revealed the heterogeneous inception times of global peatlands (MacDonald et al., 2006; Loisel et al., 2017), suggesting that it is unreasonable to treat all peatlands with a uniform duration (13 500 years) of carbon accumulation as in the old spin-up (Treat et al., 2019, 2021). To address this, we used peat inception maps created according to radiocarbon dates from 2860 geolocated peat cores (Fig. S1 in the Supplement)
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Spatial distribution of the present-day Yedoma region from Strauss et al. (2021), shown at the initial spatial resolution of 0.05° × 0.05°. The spatial extent of Yedoma in the map covers only areas north of 60°N.</figcaption></figure>

presented in Loisel et al. (2017). These maps summarize the mean peat inception ages grouped by 1° × 1° grid cells. The data from these maps are prescribed in the model to represent the observed history of peat carbon evolution. The earliest timing of peatland inception according to these peat cores was around 16 ka. During that period, deglacial warming halted Yedoma accumulation, causing the thaw of massive ice wedges in many Yedoma deposits, which resulted in ground subsidence and the formation of thermokarst lakes (Anthony et al., 2014; Brosius et al., 2021; Strauss et al., 2017). Where thaw created drainage channels, subsequent lake drainage lowered water levels, enabling the growth of mosses and other plants whose accumulation would ultimately lead to peatland formation (Anthony et al., 2014). Accordingly, in the new spin-up simulation, we disabled Yedoma sedimentation over the present-day Yedoma regions after the LGM and initiated peat carbon accumulation in suitable areas where peat inception was observed. The suitable areas and observed inception time were determined by combining the point-level peat age data (Fig. S1) and our peatland fraction map (Hugelius et al., 2020) and were transformed into gridded peatland cohort maps covering 30–90° N for each 1000-year interval since 16 ka (Fig. S2).

Based on the gridded peatland cohort maps, we computed the number of peatland grid cells and total peatland area in different periods of the pre-Holocene deglaciation and Holocene period from 16 ka to the present (Fig. 3). Since few peatlands formed during the earliest period (16–12 ka) and the most recent one (5–0 ka), compared to each kyr from 12 to 5 ka, we used the area-weighted peatland inception time as the inception time of peatlands during these two periods and reduced the simulation time from 16 to 9 kyr to save computational cost. The peatland inception times of the 9 periods or cohorts are 13.5, 11.5, 10.5, 9.5, 8.5, 7.5, 6.5, 5.5,

and 4 ka, respectively (Table 2). We also used the “FullO–SubC” strategy to accelerate the carbon accumulation in the peat profiles. For each kyr, we used the 20-year cycled climate data from CRU-JRA v2.3 (1901–1920) as in the old spin-up simulation, the CO₂ levels for each period derived from the ice-core record (obtained from https://www1.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/domec/, last access: 3 November 2023), and the land cover maps combining 15 conventional PFTs using the 1700 land cover map from TRENDY, 1 peat PFT (C₃ grass) defined by the gridded peatland cohort maps during each kyr (Fig. S2), and 1 Yedoma PFT (C₃ grass) defined by the Yedoma map from Strauss et al. (2021). We reintroduced Yedoma during this phase because our study aims to simulate Yedoma deposits that still exist today, that is, deposits that survived through the Holocene. To do this, we transferred the litter and soil carbon from the final year of the Yedoma simulation to initialize the carbon pools of the Yedoma PFT in this phase. However, we turned off sedimentation, letting the Yedoma PFT follow the same processes as conventional PFTs and accumulate carbon slowly. The fraction of the Yedoma PFT during each kyr of peat inception was fixed since this Yedoma was assumed to have survived until today. When creating the land cover maps, we prioritized Yedoma and peat fractions in each grid cell and adjusted the fractions of 15 conventional PFTs from the TRENDY land cover map proportionally. As a result, compared to the old spin-up simulation, our new spin-up simulation included carbon accumulation for Yedoma deposits and the “realistic inception history” for peatland evolution from 16 ka to the present (Fig. 1). All the simulations for this phase were run for 30–90° N at a 2° × 2° spatial resolution, consistent with the old spin-up simulation.

## Conventional soils

There is no new setup for conventional soils, where the 15 conventional PFTs (Table S1 in the Supplement) grow in the new spin-up simulation compared to the old spin-up. However, it is noteworthy that in the starting MICT version from Xi et al. (2024), the SOC simulated by the model can feed back to soil temperature – an important mechanism for simulating the formation of high-latitude SOC stocks (Koven et al., 2009) – while in the original MICT version, Guiberteau et al. (2018) used the observation-based SOC maps to represent the insulation effect of SOC on soil thermal diffusivity.

## Changes in plant and soil carbon pools following conversion from conventional soils into peatlands

In the starting MICT version from Xi et al. (2024), when conventional soils are converted into peatlands within a grid cell, the newly established peatland fraction initially inherits the plant and soil carbon pools from the displaced conventional soils to ensure mass conservation. After this transition, the
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Evolution of northern peatlands (> 30° N) from 16 ka to the present according to gridded peatland cohort maps as shown in Fig. S2. (a) Number and area of new peatland grid cells (with the initial spatial resolution of 1° × 1°) incepted during every 1000 years from 16 ka to the present. (b) Number and area of total peatland grid cells by each 1000 years from 16 ka to the present. Peatland grid cells are defined as grid cells with a > 0 peatland fraction. In panel (b), the numbers and area of total peatland grid cells used in the old spin-up are also shown, which misses the evolution of peatland extent from 13.5 ka to the present.</figcaption></figure>

peatland fraction begins to grow the peatland-specific PFT and accumulate soil carbon according to peatland soil characteristics. As described in Sect. 2.1, compared to conventional PFTs (e.g., tree or grass), the peatland PFT features a shallower rooting depth and is assigned distinct hydrological properties. Moreover, the peat soil tile does not allow drainage at the base of the soils and also receives lateral surface water input from other non-peatland PFTs within the same grid cell. These waterlogged soil conditions substantially suppress decomposition, promoting the accumulation of soil organic carbon. Peatland tiles also differ regarding vertical soil carbon transfer. Conventional PFTs can accumulate carbon down to a depth of 3 m through root-distributed inputs, depth-dependent decomposition, and vertical diffusion via bioturbation and cryoturbation (see the first paragraph of Sect. 2.1). In contrast, the peatland PFT employs a more efficient scheme for vertical carbon transfer (detailed in the second paragraph of Sect. 2.1), enabling substantial carbon accumulation down to depths of up to 10 m.

## 3 Results and interpretation

In this section, we present the simulated Yedoma carbon accumulation from the improved model version in Sect. 3.1. In Sect. 3.2, we examine how peatland soil carbon accumulation responds to the implementation of spatially varying peatland initiation timing. Since the vertical carbon transfer process in Yedoma is driven by deposition rates derived from site-level measurements and peatland development is informed by prescribed maps extrapolated from point-based peat age data, we evaluate model performance by comparing simulated soil carbon from both the old and the new spin-up schemes, as well as against site-level observations. In Sect. 3.3, we provide a spatial evaluation using the deepest available gridded

soil carbon map (to our knowledge), although it only covers the top 0–3 m of soil.

### 3.1 Yedoma carbon accumulation

Figure 3 shows the Yedoma carbon accumulation from the new spin-up simulation. With a 1 mm yr⁻¹ sedimentation rate and a 19 kyr accumulation duration under LGM climate, the model accumulates ~20 m of organic carbon on average over the present-day Yedoma region (Fig. 4b). This is very close to the mean deposit thickness of 19.4 m suggested by Strauss et al. (2013). The mean simulated carbon density over the soil column was 271 kg C m⁻², with a large spatial variation given an interquartile range of 366 kg C m⁻² (Fig. 4a). Vertically, the Yedoma SOC tended to be near-evenly distributed along soil layers, with a profile-mean volumetric carbon density of ~14 kg C m⁻³ (Fig. 4c). The spatial patterning and vertical distribution of Yedoma SOC are consistent with the outputs from Zhu et al. (2016), and the total SOC stock from the new spin-up (141 Pg C) is within the range of their estimate (125–145 Pg C). When excluding the volume of ice wedges (~50% of soil column volume) used to compute the sedimentation rate in Eq. (2), the simulated SOC density from our model (~28 kg C m⁻³) compares well with measurements, e.g., 25.98 ± 1.5 kg C m⁻³ (excluding ice) in Anthony et al. (2014) and 15.4–26.8 kg C m⁻³ (excluding ice) in Palmtag et al. (2022). Given the significant uncertainties in the Yedoma area and the scarce measurements of deep Yedoma SOC profiles, these results suggest that the new spin-up reasonably reproduces the Yedoma carbon stock.

### 3.2 Carbon accumulation in peatlands and other soils after the LGM

Regarding carbon accumulation in peat and other soils, in the old spin-up simulation, the total organic carbon (TOC, including litter and soil organic carbon) for peat and 15 con-
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Yedoma carbon accumulation from the new spin-up simulation. Spatial distribution of the simulated (a) total organic carbon (kg m⁻² Yedoma deposits) and (b) maximum depth of Yedoma deposits. (c) The simulated mean vertical distribution of the carbon density of three soil carbon pools (active, slow, and passive) over the present-day Yedoma region by the end of the LGM. The spatial extent of Yedoma coverage data in maps (a) and (b) is limited to areas north of 60°N, given that no Yedoma deposits are found south of 60°N (Fig. 2; Strauss et al., 2021).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Time series of net ecosystem productivity (NEP) and total organic carbon over the Northern Hemisphere (> 30°N) from the old (a, c) and the new (b, d) spin-up simulations. Conventional PFTs indicate the 15 non-Yedoma and non-peat PFTs including bare soil, 8 tree PFTs, 4 grass PFTs, and 2 crop PFTs. The simulation periods using the standalone soil carbon module (SubC) are not plotted in all figures as the model only outputs the simulation results of the last year of each period to save storage space.</figcaption></figure>

Conventional PFTs over the Northern Hemisphere (30–90°N) gradually accumulates from 0 to ~287 Pg C during the first simulation stage (Fig. 5c), driven by a large positive net ecosystem productivity (NEP) resulting from the excess of net primary productivity (NPP) over heterotrophic respiration (Rh) (Figs. 5a and S3). Subsequently, the long (13 500 years) integration of the standalone soil carbon module resulted in a TOC increase of ~1500 Pg C, followed by a near-stabilization period (NEP = 0.6 Pg C yr⁻¹, of which 0.4 Pg C yr⁻¹ is from conventional PFTs and 0.09 Pg C yr⁻¹ is from peat PFT) during the third simulation stage.

With its explicit representation of peatland inception times, the new spin-up simulation shows a similar temporal trend in NEP from 16 ka to the present to the old spin-up (Fig. 5b). By the end of the simulation, annual NEP is a carbon sink of 0.4 Pg C yr⁻¹ (0.3 Pg C yr⁻¹ from conventional PFTs, 0.08 Pg C yr⁻¹ from peat PFT, and 0.007 Pg C yr⁻¹ from Yedoma PFT). Accordingly, the conventional PFTs and peat PFT accumulate a total soil carbon stock of 1239 and 632 Pg C, respectively, while the Yedoma PFT holds a mean carbon stock of 144 Pg C throughout the entire simulation period (Table 3). Since we only simulated sedimentation pro-
<table><caption>**Table 3.** Summary of area and organic carbon (OC) stocks for different PFT categories by the end of the old and new spin-up simulations. The OC stocks for four PFT groups (all PFTs, conventional PFTs, peat PFT, and Yedoma PFT) are shown by four depth intervals (0–1, 1–2, 2–3, and 3–38 m), two regions (60–90 and 30–60° N), and four carbon pools (litter, active soil, slow soil, and passive soil). The combined observation-based product (WISE30sec for global, non-permafrost, 0–2 m soil layers and NCSCD (the Northern Circumpolar Soil Carbon Database) for permafrost, 0–3 m soil layers) is also shown for comparison (see details in Sect. 3.3).</caption><thead><tr><th rowspan="2">PFT category</th><th rowspan="2"></th><th colspan="3">All PFTs</th><th colspan="2">Conventional PFTs</th><th colspan="2">Peat PFT</th><th colspan="2">Yedoma PFT</th></tr><tr><th>WISE30sec<br>+ NCSCD</th><th>Old</th><th>New</th><th>Old</th><th>New</th><th>Old</th><th>New</th><th>Old</th><th>New</th></tr></thead><tbody><tr><td colspan="11">Area (×10<sup>6</sup> km<sup>2</sup>)</td></tr><tr><td></td><td></td><td>68.6</td><td>68.6</td><td>68.6</td><td>65.0</td><td>64.5</td><td>3.6</td><td>3.6</td><td>0</td><td>0.5</td></tr><tr><td colspan="11">Total OC stock (Pg C)</td></tr><tr><td></td><td></td><td>1663.3</td><td>1801.8</td><td>2027.5</td><td>1170.4</td><td>1238.6</td><td>631.4</td><td>632.0</td><td>0</td><td>157.0</td></tr><tr><td colspan="11">Litter OC stock (Pg C)</td></tr><tr><td></td><td></td><td>—</td><td>86.5</td><td>91.1</td><td>67.3</td><td>66.4</td><td>19.2</td><td>23.4</td><td>0</td><td>1.4</td></tr><tr><td colspan="11">Total</td></tr><tr><td rowspan="5">Soil OC<br>(Pg C)</td><td rowspan="5">By depth<br>interval<br>(m)</td><td>0–1</td><td>846.4</td><td>882.8</td><td>959.5</td><td>746.5</td><td>790.0</td><td>136.3</td><td>145.3</td><td>0</td><td>24.1</td></tr><tr><td>1–2</td><td>524.8</td><td>427.2</td><td>490.9</td><td>303.1</td><td>341.5</td><td>124.1</td><td>131.9</td><td>0</td><td>17.4</td></tr><tr><td>2–3</td><td>292.1</td><td>135.5</td><td>128.3</td><td>53.5</td><td>40.6</td><td>82.0</td><td>81.0</td><td>0</td><td>6.8</td></tr><tr><td>3–38</td><td>—</td><td>269.8</td><td>357.6</td><td>0</td><td>0.03</td><td>269.8</td><td>250.3</td><td>0</td><td>107.3</td></tr><tr><td>60–90° N</td><td>793.3</td><td>795.4</td><td>1084.9</td><td>473.3</td><td>557.8</td><td>322.1</td><td>371.5</td><td>0</td><td>155.6</td></tr><tr><td rowspan="2">By region</td><td>30–60° N</td><td>870.0</td><td>919.9</td><td>851.5</td><td>629.8</td><td>614.4</td><td>290.1</td><td>237.2</td><td>0</td><td>0</td></tr><tr><td>Active</td><td>—</td><td>68.9</td><td>115.5</td><td>24.0</td><td>26.3</td><td>44.8</td><td>78.6</td><td>0</td><td>10.6</td></tr><tr><td rowspan="2">By pools</td><td>Slow</td><td>—</td><td>1016.6</td><td>1242.7</td><td>529.9</td><td>630.5</td><td>486.7</td><td>484.1</td><td>0</td><td>128.1</td></tr><tr><td>Passive</td><td>—</td><td>629.9</td><td>578.2</td><td>549.2</td><td>515.4</td><td>80.7</td><td>46.0</td><td>0</td><td>16.8</td></tr></tbody></table>

cesses over the present-day Yedoma extent, it is reasonable that the simulated Yedoma carbon stocks deposited during the LGM climate could persist during the warm deglacial period. Due to increased vegetation productivity under the warmer Holocene climate compared to the LGM, the simulated Yedoma carbon stock increases by 16 Pg C through the Holocene, which matches the estimates of 13 ± 1 Pg C from Strauss et al. (2017) well. For the peat PFT, we found that the peatland carbon accumulation rate at the end of the simulation from the new spin-up (0.08 Pg C yr⁻¹) was quite close to that from the old spin-up (0.09 Pg C yr⁻¹), both of which are comparable with the observed rate (0.10 (0.02–0.19) Pg C yr⁻¹) (Li et al., 2023). The TOC for peat PFT from the two spin-ups is almost identical (Table 3). For the 15 conventional PFTs, the inclusion of Yedoma PFT in the land cover map slightly reduced their area (0.7%), while TOC is 68 Pg C (9%) higher in the new spin-up simulation (Table 3). As a result, the TOC over the Northern Hemisphere (> 30° N) for all PFTs is 2028 Pg C in the new spin-up, 226 Pg C more than the old spin-up.

To better understand the difference in TOC between the two spin-up simulations, we computed carbon stocks across four depth intervals (0–1, 1–2, 2–3, and 3–38 m), Arctic and boreal regions (60–90 and 30–60° N), and all the carbon pools (litter, active soil, slow soil, and passive soil) (Table 3). The litter OC for conventional PFTs is almost identical in the new spin-up to that in the old spin-up, but it is 4.2 Pg C (+22%) higher for peat PFT in the new spin-up. For soil carbon, both conventional PFTs and peat PFT show higher SOC

for 0–2 m soil (82 Pg C (+8%) and 17 Pg C (+16%), respectively) but lower SOC for 2–38 m soils (13 Pg C (−24%) and 21 Pg C (−6%), respectively) in the new spin-up. The higher SOC for conventional PFTs is concentrated in the Arctic region (85 Pg C, +15%), predominantly in the slow SOC pool (101 Pg C, +19%). By contrast, the peat PFT also has higher SOC (49 Pg C, +15%) in the Arctic region but lower SOC (53 Pg C, −18%) in the boreal region, which offsets and even reverses the SOC gains in the Arctic region in the new spin-up. Breaking down the PFT classes into their CENTURY-derived soil carbon quality pools, the peat PFT has a higher active SOC (34 Pg C, +75%) but a lower passive SOC (35 Pg C, −43%) in the new spin-up outputs. The lower simulated SOC for deep soils and in the boreal region, as well as less passive SOC for the peat PFT, could be attributable to the shorter duration of carbon accumulation after constraining peat inception time in the new spin-up. By comparing the peatland cells by age cohort, we found that those young peatland cells, e.g., incepted after 8–7 ka, in the new spin-up show a 6%–31% lower total SOC (Fig. 6a) and a 50%–70% lower passive SOC (Fig. 6b–d), relative to the 13 500-year-old peatlands in the old spin-up. The higher SOC for 0–2 m soils in the high latitudes for both conventional PFTs and the peat PFT in the new spin-up could be associated with the setup of the new spin-up simulation and some existing model processes, which are discussed in the Discussion section.

Figure 7 shows the spatial pattern of total organic carbon and the vertical profile of SOC density for three soil types after peat carbon accumulation simulation from the old and
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Comparison of total (a), active (b), slow (c), and passive (d) SOC by age cohorts for peat PFT between the new and old spin-up simulations. The age cohorts are identified as the peatland inception time as shown in Fig. 3 and Table 2. The numbers show the relative difference between the new and old spin-ups.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Carbon accumulation for three PFT classes after the peat simulation from the old and new spin-ups. (a–c) Spatial distribution of total organic carbon over the Northern Hemisphere (> 30°N) from the old (a) and new (b) spin-up simulations, as well as the difference between new and old spin-ups (c). (d–f) Spatial distribution of the maximum depth of soil organic carbon (SOC) from the old (d) and new (e) spin-up simulations, as well as the difference between new and old spin-ups (f). (g–i) Vertical distribution of the carbon density of three soil carbon pools (active, A; slow, S; and passive, P) for the Yedoma PFT (g), peat PFT (h), and conventional PFTs (i) from the old and new spin-ups.</figcaption></figure>

new spin-ups. Consistently with the results shown in Table 3, the higher total OC in the new spin-up is mainly distributed in the Arctic region (Fig. 7), where the inclusion of Yedoma carbon and the larger carbon stocks by conventional PFTs (Fig. S4) and peat PFT (Fig. S5) results in a > 50 kg C m⁻² higher carbon density (Fig. 7a–c) and a > 5 m deeper maximum depth of SOC than the old spin-up (Fig. 7d–f). Analyzing the cumulative percentage of SOC at depths of 1, 2, and 3 m, we found that the higher SOC in very high latitudes (> 70°N) is mainly contributed by soil layers deeper than 3 m from the Yedoma PFT (Fig. 7g), while the proportion of 0–3 m SOC in the new spin-up is 10%–20% lower than the old one (Fig. S6). In the boreal region, the shorter duration of peat carbon accumulation in the new spin-up leads to a shallower SOC depth (1–5 m) and a lower carbon density (~10 kg C m⁻²) (Figs. 7h and S5). This shallower peat carbon results in a 5%–25% increase in its share of SOC in

the 0–3 m soil layers relative to the old spin-up in this region (Fig. S6).

## 3.3 Spatial evaluation against observation-based soil organic carbon maps

To further evaluate the soil carbon simulation with the new spin-up, we compared the spatial distribution of simulated total soil carbon with two observation-based soil datasets: WISE30sec (Batjes, 2016) and NCSCD (the Northern Circumpolar Soil Carbon Database) (Hugelius et al., 2013) (Fig. 8). WISE30sec provides global coverage of soil carbon stocks but only for 0–2 m soil layers. In contrast, NCSCD is a regional dataset for permafrost regions but includes estimated soil carbon down to 3 m depth. We therefore combined the two datasets using NCSCD for permafrost regions and WISE for non-permafrost regions to evaluate simulated soil carbon stock. While the combined 0–3 m dataset is still
not sufficient to evaluate the deep carbon stored in Yedoma (~20 m) and peatlands (up to 10 m), to our knowledge it remains the deepest available gridded SOC product.

For 0–2 m soil layers, the new spin-up generally reproduces the observation-based spatial distribution of total OC. However, it overestimates OC by 16–40 kg C m⁻² across northern North America and Europe, the Yedoma region, and eastern Asia, while underestimating OC by 8–24 kg C m⁻² in East Siberia (Fig. 8c). As a result, the simulated total OC for 0–2 m soil layers over the Northern Hemisphere (> 30° N) from the new spin-up is 1450 Pg C, 79 Pg C higher than the combined observation-based soil datasets (Table 3). Despite the spatial discrepancies, the difference in total OC between our simulation and the combined observation-based datasets is well within the uncertainty range of WISE30sec (±215 Pg C) and NCSCD (±200 Pg C).

For 2–3 m soil layers in permafrost regions, the new spin-up underestimates total OC from NCSCD by 8–16 kg C m⁻² (~80%) across nearly all permafrost regions (Fig. 8f). Although the inclusion of Yedoma carbon and the higher OC for conventional and peat PFTs at high latitudes in the new spin-up (Fig. 7c) reduces the underestimation of OC from the old spin-up (Fig. S7), the total OC stock for 2–3 m soil layers in permafrost regions from the new spin-up is still only about one-third of that in the NCSCD (92 Pg C vs. 292 Pg C). This suggests that significant uncertainties remain in simulating the vertical distribution of soil carbon with MICT or that some other processes such as thermokarst deposits (Anthony et al., 2014) or deltaic deposits (Tarnocai et al., 2009) must still be incorporated into the model. However, it is also important to note that the number of high-latitude pedon samples located at sites in northern North America and Siberia and extending down to 3 m is quite limited in WISE and NCSCD (Batjes, 2016; Hugelius et al., 2014), which could also contribute to the discrepancies between the observation-based and simulated soil carbon.

# 4 Discussion

## 4.1 Effects of SOC insulation on soil carbon accumulation

Compared to the old spin-up simulation, the new spin-up protocol presented here incorporates Yedoma carbon and a “realistic inception history” of peatland evolution from 16 ka to the present. As a result, the simulated results add 157 Pg C of Yedoma carbon to the total SOC stocks and reproduce the observation-derived spatial patterns and vertical distributions of Yedoma carbon (Fig. 4). For the peat PFT, due to the shorter duration of peat carbon accumulation in the new spin-up, the simulated deep SOC stocks and the passive SOC pool are reduced by 20 Pg C (7%) and 35 Pg C (43%), respectively (Table 3), with the reduction being more obvious

in younger peatlands (Fig. 6). However, the differences between the two spin-up simulations go beyond this.

In the old spin-up, the SubC is run once for 13 500 years, whereas the new spin-up divides this time span into nine simulation periods. As described in the “Data and methods” section, the SubC run uses constant boundary conditions including monthly litter input, soil moisture, and soil thermal conditions from the last year of the preceding FullO simulation. In other words, the 13 500-year SubC run in the old spin-up uses soil temperature ($T_{soil}$) in a low-SOC-stock context of 287 Pg C (Fig. 5c). In contrast, in the new spin-up, the accumulation of OC during each kyr feeds back into $T_{soil}$, leading to a lower $T_{soil}$ in the following kyr due to the insulation effect of SOC, especially in high latitudes (Fig. 7). This lower $T_{soil}$ reduces soil carbon decomposition, promoting carbon accumulation in the new spin-up simulation. The SOC-$T_{soil}$ positive feedback (more soil carbon leads to cooler soils and increased SOC accumulation) is an important driver of higher SOC stocks derived from conventional PFTs in the new spin-up (Table 3 and Fig. S4).

The $T_{soil}$ for 0–20 cm over peatlands north of 60° N in the new spin-up is also lower than in the old spin-up in high latitudes due to the same SOC-$T_{soil}$ positive feedback (Fig. 9b), which partly explains the higher SOC in 0–2 m peat soils and in the Arctic region despite a shorter duration of carbon accumulation in the new spin-up (Table 3 and Fig. S5). Another process contributing to higher SOC stocks and deeper SOC depth (Fig. S5) in the new spin-up is the mechanism of carbon distribution due to grid-specific changes in the PFT fraction. As explained in Sect. 2.2, in MICT, carbon from shrinking PFTs is redistributed to expanding PFTs to ensure carbon conservation. Thus, when the peat PFT of a given grid cell expands in the new spin-up, it inherits SOC from conventional PFTs. Since conventional PFTs in high latitudes have high SOC density and deep SOC depth up to 3 m (Figs. 7 and S4), the expanding peatlands in the new spin-up could exhibit higher SOC and deeper SOC depths than those in the old spin-up (Fig. S5). In contrast, in middle latitudes, where the SOC inherited from conventional PFTs is modest, the peat PFT shows a 53 Pg C (18%) lower SOC, highlighting the important impacts of representing the peatland inception time in the model. Moreover, since all peatlands use the constant observation-derived constraint for the maximum carbon density in each soil layer ($C_{obs,1}$ in Eq. 1), high-carbon-input peatland areas in the model could accumulate very deep SOC. For example, in some areas of East Siberia and Alaska, the maximum SOC depth for peat PFT in the old spin-up reaches 38 m (Fig. S5d), implying limitations to applying the current parameterization of peat carbon accumulation to all peatlands. By constraining the inception age of peatlands in the new spin-up, the shortened duration of carbon accumulation significantly reduces SOC depth for some young peatlands within these regions (Fig. S5e and f).
<figure><img src="image_12.png" /></figure>

**Figure 8.** Comparison of spatial distributions of total organic carbon for 0–2 and 2–3 m soil layers over the Northern Hemisphere (> 30° N) between two observation-based soil datasets (NCSCD for 0–3 m soil layers in permafrost regions and WISE for 0–2 m soil layers in non-permafrost regions) and the new spin-up simulation. The first two columns show total organic carbon from the two observation-based soil datasets (**a**, **d**) and new spin-up simulations (**b**, **e**), and the third column (**c**, **f**) shows the difference between simulations and observations.

<figure><img src="image_13.png" /></figure>

**Figure 9.** Comparison of annual soil temperature for 0–20 cm for conventional PFTs (**a**) and peat PFT (**b**) in the Arctic region (60–90° N) between the new and old spin-up simulations. We only show the soil temperature used in the SubC run for the old spin-up, i.e., the fixed value from the last year of the preceding FullO simulation as the background (the dotted red line) for comparison.

Regardless of whether peatlands are initiated simultaneously or spatially variably, the simulated total northern peatland carbon stock (~630 Pg C) is at the high end compared to previous estimates, e.g., 547 ± 74 Pg C from Yu et al. (2010), 415 ± 150 Pg C from Hugelius et al. (2020), 410 (315–590) Pg C from Treat et al. (2019), and 463 Pg C from the earlier ORCHIDEE-MICT-PEAT version in Qiu et al. (2019). This overestimation may also be attributed to the SOC–$T_{\text{soil}}$ feedback (or coupling). As evaluated in Xi et al. (2024), wetter and higher-SOC peat soils tend to maintain lower temperatures than the grid-cell average. When SOC–

$T_{\text{soil}}$ coupling is implemented at the PFT level (i.e., tiled coupling) in the base model (ORCHIDEE-MICT-teb), the simulated peat SOC stock increases to 620 Pg C, compared to 534 Pg C with grid-cell-level coupling, i.e., an 86 Pg C difference. This suggests that further evaluation and calibration of peat SOC in the improved model version are still needed in the future.

## 4.2 Implications and limitations

Incorporating high-latitude deep carbon into LSMs or ESMs is needed to quantify its response to climatic change, given
the vast carbon storage and the accelerated warming in these regions. However, it remains inadequately addressed, or even entirely absent, in current ESMs (Schädel et al., 2024). Although far from complete, the improvements to soil carbon accumulation from this study, compared to the original MICT, offer an acceptable model tool capable of simulating both shallow and deep soil carbon for conventional soils, peatlands, and Yedoma deposits consistently. The next step in the application of this model will be to investigate if the inclusion or modification of deep carbon in the new spin-up will impact projected soil carbon emissions and whether this will alter the shift timing of the terrestrial transition from a net carbon sink to carbon source in the years to come (Koven et al., 2011, 2015; Nitzbon et al., 2024). This is a critical scientific question in permafrost carbon feedback research, as the highly confident projections of active layer thickening under future climate change (IPCC, 2021) are expected to increase the exposure of deep carbon to higher temperatures, thereby accelerating permafrost carbon loss. Importantly, none of the ESMs used in the current IPCC AR6 report explicitly simulate deep carbon dynamics (IPCC, 2021).

Additionally, an interesting finding from this study is that the improvements of deep carbon representation alter the distribution of simulated SOC quality, i.e., the fractions of active, slow, and passive SOC (Fig. 10). In the CENTURY model, decomposed active SOC is partially released as CO2 and partially transferred to slow and passive carbon pools. Similarly, decomposed slow SOC partly moves to the passive carbon pool. Since the additional SOC in the new spin-up is concentrated in very high latitudes (Fig. 7c), where the cold environments inhibit the SOC flow from active and slow pools to the passive pool, this results in 2% and 4% higher active and slow SOC, respectively, but 6% lower passive SOC over the Northern Hemisphere (> 30°N) in the new spin-up (Fig. 10). The higher fractions of active and slow SOC are more evident in higher latitudes. In areas north of 70°N, the new spin-up shows 3% and 7% higher active and slow SOC, respectively, compared to the old spin-up (Fig. 10). As deep carbon becomes increasingly exposed to a warming-induced thickening of the active layer, the higher fractions of active and slow SOC in the high latitudes, as well as in northern peatlands as mentioned before (Fig. 6), simulated by the new spin-up could amplify the projected increase in soil carbon emissions.

This study contains several potential limitations. First, we only simulated the sedimentation process over present-day Yedoma deposits, while the formation and evolution of Yedoma deposits have taken place throughout the last glacial–deglacial transition. Coupling with paleoclimate models and applying the sedimentation-enabled model to specific periods could help reveal the role of Yedoma deposits in the variation in atmospheric CO2 concentration during the last glacial–deglacial transition (Ciais et al., 2012). Second, in simulating peat carbon accumulation since 16 ka, we used recycled 1901–1920 climate data to force the model

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Difference in SOC fractions between the new and old spin-up simulations for areas > 30, > 50, > 60, and > 70°N. The four areas include 100%, 79%, 55%, and 11% of the total SOC of the study area, respectively, in the new spin-up simulation.</figcaption></figure>

in the new spin-up, in order to make a comparison with the old spin-up simulation. Using the actual climate data could increase carbon accumulation during the cold pre-Holocene deglacial period, reduce accumulation during the warm mid-Holocene, and increase accumulation again during the cold late Holocene (Kaufman et al., 2020; Osman et al., 2021). Coupling with paleoclimate models (Kleinen et al., 2012, 2016; Treat et al., 2019) or using available datasets such as TraCE-21ka (He et al., 2013), along with incorporating the dynamic evolution of conventional PFTs (e.g., trees, grasses, and croplands) during the last 20 kyr, would provide a more realistic carbon accumulation history for peatlands and conventional soils. Third, as we used present-day peatland distributions to reconstruct Holocene peatland maps, our simulations do not account for peatlands that were historically drained for human land use. Within our study region this is particularly relevant to Northern Europe, where substantial peatland drainage for agriculture and forestry took place over the past two centuries (Fluet-Chouinard et al., 2023). A previous study has shown that drainage or peatland-to-cropland conversion during the years 850–2010 could have increased total historical carbon stocks in northern peatlands by 72 Pg C (Qiu et al., 2021). However, because the development presented in this study is primarily aimed at projecting future changes in deep carbon stored in peatlands, omitting historical carbon losses from drainage is likely to have a relatively limited influence on future projections, as those losses mostly contributed to past emissions.

## 5 Conclusion

This study presents a new spin-up protocol for high-latitude carbon formation over the late Pleistocene and Holocene by representing deep carbon accumulation in Yedoma deposits and peatlands more realistically in the ORCHIDEE-MICT land surface model. Compared to the existing spin-up approach, the simulation outputs of the new scheme add 157 Pg C to the present-day Yedoma region, while deep carbon below 3 m is reduced by 20 Pg C and passive soil carbon
is reduced by 35 Pg C in northern peatlands due to the shorter duration of peat carbon accumulation after an observation-constrained peat initiation time. Since the inclusion of SOC-T<sub>soil</sub> feedback in carbon accumulation simulation slightly increases the carbon stock in conventional soils by 68 Pg C, the total organic carbon stock across the Northern Hemisphere (> 30° N) simulated by the new spin-up is 2028 Pg C, which is 226 Pg C larger than the old spin-up. Although evaluation against an observation-based soil dataset reveals significant challenges in simulating soil carbon for 2–3 m with ORCHIDEE-MICT, the improvements in carbon accumulation modeling from this study provide a model version for the simulation of deep land carbon evolution during the last glacial–deglacial transition and for the projection of the land carbon response to future warming.

*Code and data availability.* The ORCHIDEE-MICT-teb-deepc model (r8704) code used in this study is open-source and distributed under the CeCILL (CEA CNRS INRIA Logiciel Libre) license. It is deposited at https://forge.ipsl.jussieu.fr/orchidee/wiki/GroupActivities/CodeAvalabilityPublication/ ORCHIDEE-MICT-teb-deepc (last access: 1 November 2024) and archived at https://doi.org/10.14768/704f7a4e-3a1d-4a51-8c70- 1920c32684ce (Xi, 2025a), with guidance to install and run the model at https://forge.ipsl.jussieu.fr/orchidee/wiki/Documentation/ UserGuide (last access: 4 October 2023). Due to limitations imposed by their size, the model input files have not been uploaded to the public repository; however, they can be accessed from the corresponding author upon reasonable request. The Yedoma map can be obtained from https://doi.org/10.1594/PANGAEA.940078 (Strauss et al., 2022). The radiocarbon dates of peat cores can be obtained from https://doi.org/10.1594/PANGAEA.902549 (Loisel et al., 2019). The initial WISE30sec dataset can be obtained from https://data.isric.org/geonetwork/srv/api/ records/dc7b283a-8f19-45e1-aaed-e9bd515119bc (Batjes, 2016). The initial NCSCD dataset can be obtained from https://bolin.su.se/data/ncscd/ (Hugelius et al., 2013). The resampled data of WISE30sec and NCSCD used in this study and codes to process data, generate all results, and produce all figures are archived at https://doi.org/10.5281/zenodo.15371113 (Xi, 2025b).

*Supplement.* The supplement related to this article is available online at https://doi.org/10.5194/gmd-18-6043-2025-supplement.

*Author contributions.* YX designed the new spin-up protocol and implemented the model development. PC, DZ, CQ, YZ, SP, and GF provided general scientific guidance to improve the research and interpret results. YX performed the simulations, did the analysis, and wrote the paper. All authors contributed to commenting on and writing the manuscript.

*Competing interests.* The contact author has declared that none of the authors has any competing interests.

*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

*Acknowledgements.* The authors are grateful to the ORCHIDEE group for their kind help on coding, preparing input files, and publishing the model version on SVN.

*Financial support.* This research was supported by the CALIPSO (Carbon Loss in Plant Soils and Oceans) project, funded through the generosity of Eric and Wendy Schmidt by recommendation of the Schmidt Futures program.

*Review statement.* This paper was edited by Marko Scholze and reviewed by Zicheng Yu and three anonymous referees.

## References

Anthony, K. M. W., Zimov, S. A., Grosse, G., Jones, M. C., Anthony, P. M., Iii, F. S. C., Finlay, J. C., Mack, M. C., Davydov, S., Frenzel, P., and Frolking, S.: A shift of thermokarst lakes from carbon sources to sinks during the Holocene epoch, Nature, 511, 452–456, https://doi.org/10.1038/nature13560, 2014.

Batjes, N. H.: Harmonized soil property values for broad-scale modelling (WISE30sec) with estimates of global soil carbon stocks, Geoderma, 269, 61–68, https://doi.org/10.1016/j.geoderma.2016.01.034, 2016 (data available at: https://data.isric.org/geonetwork/srv/api/records/ dc7b283a-8f19-45e1-aaed-e9bd515119bc, last access: 3 March 2024).

Bradford, M. A., Wieder, W. R., Bonan, G. B., Fierer, N., Raymond, P. A., and Crowther, T. W.: Managing uncertainty in soil carbon feedbacks to climate change, Nat. Clim. Change, 6, 751–758, https://doi.org/10.1038/nclimate3071, 2016.

Brosius, L. S., Anthony, K. M. W., Treat, C. C., Lenz, J., Jones, M. C., Bret-Harte, M. S., and Grosse, G.: Spatiotemporal patterns of northern lake formation since the Last Glacial Maximum, Quaternary Sci. Rev., 253, 106773, https://doi.org/10.1016/j.quascirev.2020.106773, 2021.

Ciais, P., Tagliabue, A., Cuntz, M., Bopp, L., Scholze, M., Hoffmann, G., Lourantou, A., Harrison, S. P., Prentice, I. C., Kelley, D. I., Koven, C., and Piao, S. L.: Large inert carbon pool in the terrestrial biosphere during the Last Glacial Maximum, Nat. Geosci., 5, 74–79, https://doi.org/10.1038/ngeo1324, 2012.

Crowther, T. W., Todd-Brown, K. E. O., Rowe, C. W., Wieder, W. R., Carey, J. C., Machmuller, M. B., Snoek, B. L., Fang, S., Zhou, G., Allison, S. D., Blair, J. M., Bridgham, S. D., Burton, A. J.,
Carrillo, Y., Reich, P. B., Clark, J. S., Classen, A. T., Dijkstra, F. A., Elberling, B., Emmett, B. A., Estiarte, M., Frey, S. D., Guo, J., Harte, J., Jiang, L., Johnson, B. R., Kröel-Dulay, G., Larsen, K. S., Laudon, H., Lavallee, J. M., Luo, Y., Lupascu, M., Ma, L. N., Marhan, S., Michelsen, A., Mohan, J., Niu, S., Pendall, E., Peñuelas, J., Pfeifer-Meister, L., Poll, C., Reinsch, S., Reynolds, L. L., Schmidt, I. K., Sistla, S., Sokol, N. W., Templer, P. H., Treseder, K. K., Welker, J. M., and Bradford, M. A.: Quantifying global soil carbon losses in response to warming, Nature, 540, 104–108, https://doi.org/10.1038/nature20150, 2016.

Fluet-Chouinard, E., Stocker, B. D., Zhang, Z., Malhotra, A., Melton, J. R., Poulter, B., Kaplan, J. O., Goldewijk, K. K., Siebert, S., Minayeva, T., Hugelius, G., Joosten, H., Barthelmes, A., Prigent, C., Aires, F., Hoyt, A. M., Davidson, N., Finlayson, C. M., Lehner, B., Jackson, R. B., and McIntyre, P. B.: Extensive global wetland loss over the past three centuries, Nature, 614, 281–286, https://doi.org/10.1038/s41586-022-05572-6, 2023.

Friedlingstein, P., Jones, M. W., O'Sullivan, M., Andrew, R. M., Bakker, D. C. E., Hauck, J., Le Quéré, C., Peters, G. P., Peters, W., Pongratz, J., Sitch, S., Canadell, J. G., Ciais, P., Jackson, R. B., Alin, S. R., Anthoni, P., Bates, N. R., Becker, M., Bel-louin, N., Bopp, L., Chau, T. T. T., Chevallier, F., Chini, L. P., Cronin, M., Currie, K. I., Decharme, B., Djeutchouang, L. M., Dou, X., Evans, W., Feely, R. A., Feng, L., Gasser, T., Gilfil-lan, D., Gkritzalis, T., Grassi, G., Gregor, L., Gruber, N., Gürses, Ö., Harris, I., Houghton, R. A., Hurtt, G. C., Iida, Y., Ilyina, T., Luijkx, I. T., Jain, A., Jones, S. D., Kato, E., Kennedy, D., Klein Goldewijk, K., Knauer, J., Korsbakken, J. I., Körtzinger, A., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lienert, S., Liu, J., Marland, G., McGuire, P. C., Melton, J. R., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S. I., Niwa, Y., Ono, T., Pier-rot, D., Poulter, B., Rehder, G., Resplandy, L., Robertson, E., Rödenbeck, C., Rosan, T. M., Schwinger, J., Schwingshackl, C., Séférrian, R., Sutton, A. J., Sweeney, C., Tanhua, T., Tans, P. P., Tian, H., Tilbrook, B., Tubiello, F., van der Werf, G. R., Vuichard, N., Wada, C., Wanninkhof, R., Watson, A. J., Willis, D., Wiltshire, A. J., Yuan, W., Yue, C., Yue, X., Zaehle, S., and Zeng, J.: Global Carbon Budget 2021, Earth Syst. Sci. Data, 14, 1917–2005, https://doi.org/10.5194/essd-14-1917-2022, 2022.

Grosse, G., Robinson, J. E., Bryant, R., Taylor, M. D., Harper, W., DeMasi, A., Kyker-Snowman, E., Veremeeva, A., Schirrmeister, L., and Harden, J.: Distribution of late Pleistocene ice-rich syngenetic permafrost of the Yedoma Suite in east and central Siberia, Russia, US Geol. Surv. Open File Rep. 2013, US Geological Survey, 1–37, https://pubs.usgs.gov/of/2013/1078/ (last access: 3 October 2023), 2013.

Guimberteau, M., Zhu, D., Maignan, F., Huang, Y., Yue, C., Dantec-Nédélec, S., Ottlé, C., Jornet-Puig, A., Bastos, A., Laurent, P., Goll, D., Bowring, S., Chang, J., Guenet, B., Tifafi, M., Peng, S., Krinner, G., Ducharne, A., Wang, F., Wang, T., Wang, X., Wang, Y., Yin, Z., Lauerwald, R., Joetzjer, E., Qiu, C., Kim, H., and Ciais, P.: ORCHIDEE-MICT (v8.4.1), a land surface model for the high latitudes: model description and validation, Geosci. Model Dev., 11, 121–163, https://doi.org/10.5194/gmd-11-121-2018, 2018.

He, F., Shakun, J. D., Clark, P. U., Carlson, A. E., Liu, Z., Otto-Bliesner, B. L., and Kutzbach, J. E.: Northern Hemisphere forcing of Southern Hemisphere climate during the last deglaciation, Nature, 494, 81–85, https://doi.org/10.1038/nature11822, 2013.

Hengl, T., Mendes de Jesus, J., Heuvelink, G. B. M., Ruiperez González, M., Kilibarda, M., Blagotić, A., Shangguan, W., Wright, M. N., Geng, X., Bauer-Marschallinger, B., Guévara, M. A., Vargas, R., MacMillan, R. A., Batjes, N. H., Leenaars, J. G. B., Ribeiro, E., Wheeler, I., Mantel, S., and Kempen, B.: SoilGrids250m: Global gridded soil information based on machine learning, PLOS ONE, 12, e0169748, https://doi.org/10.1371/journal.pone.0169748, 2017.

Hiederer, R. and Köchy, M.: Global Soil Organic Carbon Estimates and the Harmonized World Soil Database, JRC Publ. Repos., https://doi.org/10.2788/13267, 2012.

Huang, Y., Zhu, D., Ciais, P., Guenet, B., Huang, Y., Goll, D. S., Guimberteau, M., Jornet-Puig, A., Lu, X., and Luo, Y.: Matrix-Based Sensitivity Assessment of Soil Organic Carbon Storage: A Case Study from the ORCHIDEE-MICT Model, J. Adv. Model. Earth Syst., 10, 1790–1808, https://doi.org/10.1029/2017MS001237, 2018.

Hugelius, G., Tarnocai, C., Broll, G., Canadell, J. G., Kuhry, P., and Swanson, D. K.: The Northern Circumpolar Soil Carbon Database: spatially distributed datasets of soil coverage and soil carbon storage in the northern permafrost regions, Earth Syst. Sci. Data, 5, 3–13, https://doi.org/10.5194/essd-5-3-2013, 2013 (data available at: https://bolin.su.se/data/ncscd/, last access: 3 March 2024).

Hugelius, G., Strauss, J., Zubrzycki, S., Harden, J. W., Schuur, E. A. G., Ping, C.-L., Schirrmeister, L., Grosse, G., Michaelson, G. J., Koven, C. D., O'Donnell, J. A., Elberling, B., Mishra, U., Camill, P., Yu, Z., Palmtag, J., and Kuhry, P.: Estimated stocks of circumpolar permafrost carbon with quantified uncertainty ranges and identified data gaps, Biogeosciences, 11, 6573–6593, https://doi.org/10.5194/bg-11-6573-2014, 2014.

Hugelius, G., Loisel, J., Chadburn, S., Jackson, R. B., Jones, M., MacDonald, G., Marushchak, M., Olefeldt, D., Packalen, M., Siewert, M. B., Treat, C., Turetsky, M., Voigt, C., and Yu, Z.: Large stocks of peatland carbon and nitrogen are vulnerable to permafrost thaw, P. Natl. Acad. Sci. USA, 117, 20438–20446, https://doi.org/10.1073/pnas.1916387117, 2020.

IPCC: Climate Change 2021: The Physical Science Basis, in: Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J. B. R., Maycock, T. K., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, https://doi.org/10.1017/9781009157896, 2021.

Jackson, R. B., Lajtha, K., Crow, S. E., Hugelius, G., Kramer, M. G., and Piñeiro, G.: The Ecology of Soil Carbon: Pools, Vulnerabilities, and Biotic and Abiotic Controls, Annu. Rev. Ecol. Evol. Syst., 48, 419–445, https://doi.org/10.1146/annurev-ecolsys-112414-054234, 2017.

Jafarov, E. and Schaefer, K.: The importance of a surface organic layer in simulating permafrost thermal and carbon dynamics, The Cryosphere, 10, 465–475, https://doi.org/10.5194/tc-10-465-2016, 2016.

Jobbagy, E. G. and Jackson, R. B.: The Vertical Distribution of Soil Organic Carbon and Its Relation to Climate and Vegetation, Ecol. Appl., 10, 423–436, https://doi.org/10.2307/2641104, 2000.
Kageyama, M., Braconnot, P., Bopp, L., Caubel, A., Foujols, M.-A., Guilyardi, E., Khodri, M., Lloyd, J., Lombard, F., Mariotti, V., Marti, O., Roy, T., and Woillez, M.-N.: Mid-Holocene and Last Glacial Maximum climate simulations with the IPSL model – part I: comparing IPSL_CM5A to IPSL_CM4, Clim. Dynam., 40, 2447–2468, https://doi.org/10.1007/s00382-012-1488-8, 2013.

Kaufman, D., McKay, N., Routson, C., Erb, M., Dätwyler, C., Sommer, P. S., Heiri, O., and Davis, B.: Holocene global mean surface temperature, a multi-method reconstruction approach, Sci. Data, 7, 201, https://doi.org/10.1038/s41597-020-0530-7, 2020.

Kleinen, T., Brovkin, V., and Schuldt, R. J.: A dynamic model of wetland extent and peat accumulation: results for the Holocene, Biogeosciences, 9, 235–248, https://doi.org/10.5194/bg-9-235-2012, 2012.

Kleinen, T., Brovkin, V., and Munhoven, G.: Modelled interglacial carbon cycle dynamics during the Holocene, the Eemian and Marine Isotope Stage (MIS) 11, Clim Past, 12, 2145–2160, https://doi.org/10.5194/cp-12-2145-2016, 2016.

Koven, C., Friedlingstein, P., Ciais, P., Khvorostyanov, D., Krinner, G., and Tarnocai, C.: On the formation of high-latitude soil carbon stocks: Effects of cryoturbation and insulation by organic matter in a land surface model, Geophys. Res. Lett., 36, L21501, https://doi.org/10.1029/2009GL040150, 2009.

Koven, C. D., Ringeval, B., Friedlingstein, P., Ciais, P., Cadule, P., Khvorostyanov, D., Krinner, G., and Tarnocai, C.: Permafrost carbon-climate feedbacks accelerate global warming, P. Natl. Acad. Sci. USA, 108, 14769–14774, https://doi.org/10.1073/pnas.1103910108, 2011.

Koven, C. D., Schuur, E. a. G., Schädel, C., Bohn, T. J., Burke, E. J., Chen, G., Chen, X., Ciais, P., Grosse, G., Harden, J. W., Hayes, D. J., Hugelius, G., Jafarov, E. E., Krinner, G., Kuhry, P., Lawrence, D. M., MacDougall, A. H., Marchenko, S. S., McGuire, A. D., Natali, S. M., Nicolsky, D. J., Olefeldt, D., Peng, S., Romanovsky, V. E., Schaefer, K. M., Strauss, J., Treat, C. C., and Turetsky, M.: A simplified, data-constrained approach to estimate the permafrost carbon-climate feedback, Philos. T. Roy. Soc. A, 373, 20140423, https://doi.org/10.1098/rsta.2014.0423, 2015.

Krinner, G., Viovy, N., de Noblet-Ducoudré, N., Ogée, J., Polcher, J., Friedlingstein, P., Ciais, P., Sitch, S., and Prentice, I. C.: A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system, Glob. Biogeochem. Cy., 19, GB1015, https://doi.org/10.1029/2003gb002199, 2005.

Li, J., Jiang, M., Pei, J., Fang, C., Li, B., and Nie, M.: Convergence of carbon sink magnitude and water table depth in global wetlands, Ecol. Lett., 26, 797–804, https://doi.org/10.1111/ele.14199, 2023.

Loisel, J., van Bellen, S., Pelletier, L., Talbot, J., Hugelius, G., Karran, D., Yu, Z., Nichols, J., and Holmquist, J.: Insights and issues with estimating northern peatland carbon stocks and fluxes since the Last Glacial Maximum, Earth-Sci. Rev., 165, 59–80, https://doi.org/10.1016/j.earscirev.2016.12.001, 2017.

Loisel, J.: “MGK” database of peat inception (basal $^{14}$C ages) across the northern peatland domain, PANGAEA [data set], https://doi.org/10.1594/PANGAEA.902549, 2019.

Lu, M., Zhou, X., Yang, Q., Li, H., Luo, Y., Fang, C., Chen, J., Yang, X., and Li, B.: Responses of ecosystem carbon cycle to

experimental warming: a meta-analysis, Ecology, 94, 726–738, https://doi.org/10.1890/12-0279.1, 2013.

Lurton, T., Balkanski, Y., Bastrikov, V., Bekki, S., Bopp, L., Braconnot, P., Brockmann, P., Cadule, P., Contoux, C., Cozic, A., Cugnet, D., Dufresne, J.-L., Éthé, C., Foujols, M.-A., Ghattas, J., Hauglustaine, D., Hu, R.-M., Kageyama, M., Khodri, M., Lebas, N., Levavasseur, G., Marchand, M., Ottlé, C., Peylin, P., Sima, A., Szopa, S., Thiéblemont, R., Vuichard, N., and Boucher, O.: Implementation of the CMIP6 Forcing Data in the IPSL-CM6A-LR Model, J. Adv. Model. Earth Syst., 12, e2019MS001940, https://doi.org/10.1029/2019MS001940, 2020.

MacDonald, G. M., Beilman, D. W., Kremenetski, K. V., Sheng, Y., Smith, L. C., and Velichko, A. A.: Rapid Early Development of Circumarctic Peatlands and Atmospheric CH$_4$ and CO$_2$ Variations, Science, 314, 285–288, https://doi.org/10.1126/science.1131722, 2006.

Mishra, U., Hugelius, G., Shelef, E., Yang, Y., Strauss, J., Lupachev, A., Harden, J. W., Jastrow, J. D., Ping, C.-L., Riley, W. J., Schuur, E. A. G., Matamala, R., Siewert, M., Nave, L. E., Koven, C. D., Fuchs, M., Palmtag, J., Kuhry, P., Treat, C. C., Zubrzycki, S., Hoffman, F. M., Elberling, B., Camill, P., Veremeeva, A., and Orr, A.: Spatial heterogeneity and environmental predictors of permafrost region soil organic carbon stocks, Sci. Adv., 7, eaaz5236, https://doi.org/10.1126/sciadv.aaz5236, 2021.

Nitzbon, J., Schneider von Deimling, T., Aliyeva, M., Chadburn, S. E., Grosse, G., Laboor, S., Lee, H., Lohmann, G., Steinert, N. J., Stuenzi, S. M., Werner, M., Westermann, S., and Langer, M.: No respite from permafrost-thaw impacts in the absence of a global tipping point, Nat. Clim. Change, 14, 573–585, https://doi.org/10.1038/s41558-024-02011-4, 2024.

Niu, S., Chen, W., Liáng, L. L., Sierra, C. A., Xia, J., Wang, S., Heskel, M., Patel, K. F., Bond-Lamberty, B., Wang, J., Yvon-Durocher, G., Kirschbaum, M. U. F., Atkin, O. K., Huang, Y., Yu, G., and Luo, Y.: Temperature responses of ecosystem respiration, Nat. Rev. Earth Environ., 5, 559–571, https://doi.org/10.1038/s43017-024-00569-3, 2024.

Osman, M. B., Tierney, J. E., Zhu, J., Tardif, R., Hakim, G. J., King, J., and Poulsen, C. J.: Globally resolved surface temperatures since the Last Glacial Maximum, Nature, 599, 239–244, https://doi.org/10.1038/s41586-021-03984-4, 2021.

Palmtag, J., Obu, J., Kuhry, P., Richter, A., Siewert, M. B., Weiss, N., Westermann, S., and Hugelius, G.: A high spatial resolution soil carbon and nitrogen dataset for the northern permafrost region based on circumpolar land cover upscaling, Earth Syst. Sci. Data, 14, 4095–4110, https://doi.org/10.5194/essd-14-4095-2022, 2022.

Parton, W. J., Stewart, J. W. B., and Cole, C. V.: Dynamics of C, N, P and S in grassland soils: a model, Biogeochemistry, 5, 109–131, https://doi.org/10.1007/bf02180320, 1988.

Qiu, C., Zhu, D., Ciais, P., Guenet, B., Krinner, G., Peng, S., Au-rela, M., Bernhofer, C., Brümmer, C., Bret-Harte, S., Chu, H., Chen, J., Desai, A. R., Dušek, J., Euskirchen, E. S., Fortuniak, K., Flanagan, L. B., Friborg, T., Grygoruk, M., Gogo, S., Grünwald, T., Hansen, B. U., Holl, D., Humphreys, E., Hurkuck, M., Kiely, G., Klatt, J., Kutzbach, L., Largeron, C., Laggoun-Défarge, F., Lund, M., Lafleur, P. M., Li, X., Mammarella, I., Merbold, L., Nilsson, M. B., Olejnik, J., Ottosson-Löfvenius, M., Oechel, W., Parmentier, F. J. W., Peichl, M., Pirk, N., Pel-tola, O., Pawlak, W., Rasse, D., Rinne, J., Shaver, G., Schmid,
H. P., Sottocornola, M., Steinbrecher, R., Sachs, T., Urbaniak, M., Zona, D., and Ziemblinska, K.: ORCHIDEE-PEAT (revision 4596), a model for northern peatland CO₂ water, and energy fluxes on daily to annual scales, Geosci. Model Dev., 11, 497–519, https://doi.org/10.5194/gmd-11-497-2018, 2018.

Qiu, C. J., Zhu, D., Ciais, P., Guenet, B., Peng, S. S., Krinner, G., Tootchi, A., Ducharne, A., and Hastie, A.: Modelling northern peatland area and carbon dynamics since the Holocene with the ORCHIDEE-PEAT land surface model (SVN r5488), Geosci. Model Dev., 12, 2961–2982, https://doi.org/10.5194/gmd-12-2961-2019, 2019.

Qiu, C., Ciais, P., Zhu, D., Guenet, B., Peng, S., Petrescu, A. M. R., Lauerwald, R., Makowski, D., Gallego-Sala, A. V., Charman, D. J., and Brewer, S. C.: Large historical carbon emissions from cultivated northern peatlands, Sci. Adv., 7, eabf1332, https://doi.org/10.1126/sciadv.abf1332, 2021.

Schädel, C., Rogers, B. M., Lawrence, D. M., Koven, C. D., Brovkin, V., Burke, E. J., Genet, H., Huntzinger, D. N., Jafarov, E., McGuire, A. D., Riley, W. J., and Natali, S. M.: Earth system models must include permafrost carbon processes, Nat. Clim. Change, 14, 114–116, https://doi.org/10.1038/s41558-023-01909-9, 2024.

Shi, Z., Hoffman, F. M., Xu, M., Mishra, U., Allison, S. D., Zhou, J., and Randerson, J. T.: Global-Scale Convergence Obscures Inconsistencies in Soil Carbon Change Predicted by Earth System Models, AGU Adv., 5, e2023AV001068, https://doi.org/10.1029/2023AV001068, 2024.

Sitch, S., Friedlingstein, P., Gruber, N., Jones, S. D., Murray-Tortarolo, G., Ahlström, A., Doney, S. C., Graven, H., Heinze, C., Huntingford, C., Levis, S., Levy, P. E., Lomas, M., Poulter, B., Viovy, N., Zaehle, S., Zeng, N., Arneth, A., Bonan, G., Bopp, L., Canadell, J. G., Chevallier, F., Ciais, P., Ellis, R., Gloor, M., Peylin, P., Piao, S. L., Le Quéré, C., Smith, B., Zhu, Z., and Myneni, R.: Recent trends and drivers of regional sources and sinks of carbon dioxide, Biogeosciences, 12, 653–679, https://doi.org/10.5194/bg-12-653-2015, 2015.

Strauss, J., Schirrmeister, L., Grosse, G., Wetterich, S., Ulrich, M., Herzschuh, U., and Hubberten, H.-W.: The deep permafrost carbon pool of the Yedoma region in Siberia and Alaska, Geophys. Res. Lett., 40, 6165–6170, https://doi.org/10.1002/2013GL058088, 2013.

Strauss, J., Schirrmeister, L., Grosse, G., Fortier, D., Hugelius, G., Knoblauch, C., Romanovsky, V., Schädel, C., Schneider von Deimling, T., Schuur, E. A. G., Shmelev, D., Ulrich, M., and Veremeeva, A.: Deep Yedoma permafrost: A synthesis of depositional characteristics and carbon vulnerability, Earth-Sci. Rev., 172, 75–86, https://doi.org/10.1016/j.earscirev.2017.07.007, 2017.

Strauss, J., Laboor, S., Schirrmeister, L., Fedorov, A. N., Fortier, D., Froese, D., Fuchs, M., Günther, F., Grigoriev, M., Harden, J., Hugelius, G., Jongejans, L. L., Kanevskiy, M., Kholodov, A., Kunitsky, V., Kraev, G., Lozhkin, A., Rivkina, E., Shur, Y., Siegert, C., Spektor, V., Streletskaya, I., Ulrich, M., Vartanyan, S., Veremeeva, A., Anthony, K. W., Wetterich, S., Zimov, N., and Grosse, G.: Circum-Arctic Map of the Yedoma Permafrost Domain, Front. Earth Sci., 9, 758360, https://doi.org/10.3389/feart.2021.758360, 2021.

Strauss, J., Laboor, S., Schirrmeister, L., Fedorov, A. N., Fortier, D., Froese, D. G., Fuchs, M., Günther, F., Grigoriev, M. N.,

Harden, J. W., Hugelius, G., Jongejans, L. L., Kanevskiy, M. Z., Kholodov, A. L., Kunitsky, V., Kraev, G., Lozhkin, A. V., Rivkina, E., Shur, Y., Siegert, C., Spektor, V., Streletskaya, I., Ulrich, M., Vartanyan, S. L., Veremeeva, A., Walter Anthony, K. M., Wetterich, S., Zimov, N. S., and Grosse, G.: Database of Ice-Rich Yedoma Permafrost Version 2 (IRYP v2), PANGAEA [data set], https://doi.org/10.1594/PANGAEA.940078, 2022.

Tarnocai, C., Canadell, J. G., Schuur, E. A. G., Kuhry, P., Mazhitova, G., and Zimov, S.: Soil organic carbon pools in the northern circumpolar permafrost region, Global Biogeochem. Cy., 23, GB2023, https://doi.org/10.1029/2008gb003327, 2009.

Treat, C. C., Kleinen, T., Broothaerts, N., Dalton, A. S., Dommain, R., Douglas, T. A., Drexler, J. Z., Finkelstein, S. A., Grosse, G., Hope, G., Hutchings, J., Jones, M. C., Kuhry, P., Lacourse, T., Lähteenoja, O., Loisel, J., Notebaert, B., Payne, R. J., Peeteet, D. M., Sannel, A. B. K., Stelling, J. M., Strauss, J., Swindles, G. T., Talbot, J., Tarnocai, C., Verstraeten, G., Williams, C. J., Xia, Z., Yu, Z., Väliranta, M., Hättestrand, M., Alexanderson, H., and Brovkin, V.: Widespread global peatland establishment and persistence over the last 130,000 y, 116, 4822–4827, https://doi.org/10.1073/pnas.1813305116, 2019.

Treat, C. C., Jones, M. C., Brosius, L., Grosse, G., Walter Anthony, K., and Frolking, S.: The role of wetland expansion and successional processes in methane emissions from northern wetlands during the Holocene, Quaternary Sci. Rev., 257, 106864, https://doi.org/10.1016/j.quascirev.2021.106864, 2021.

Walker, T. W. N., Kaiser, C., Strasser, F., Herbold, C. W., Leblans, N. I. W., Woebken, D., Janssens, I. A., Sigurdsson, B. D., and Richter, A.: Microbial temperature sensitivity and biomass change explain soil carbon loss with warming, Nat. Clim. Change, 8, 885–889, https://doi.org/10.1038/s41558-018-0259-x, 2018.

Wieder, W. R., Bonan, G. B., and Allison, S. D.: Global soil carbon projections are improved by modelling microbial processes, Nat. Clim. Change, 3, 909–912, https://doi.org/10.1038/nclimate1951, 2013.

Xi, Y.: Model code for ORCHIDEE-MICT-teb-deepc (r8704), IPSL Data Catalog [code], https://doi.org/10.14768/704f7a4e-3a1d-4a51-8c70-1920c32684ce, 2025a.

Xi, Y.: Data and code for GMDD paper: Representing high-latitude deep carbon in the pre-industrial state of the ORCHIDEE-MICT land surface model (r8704), Zenodo [code and data set], https://doi.org/10.5281/zenodo.15371113, 2025b.

Xi, Y., Qiu, C., Zhang, Y., Zhu, D., Peng, S., Hugelius, G., Chang, J., Salmon, E., and Ciais, P.: Assessment of a tiling energy budget approach in a land surface model, ORCHIDEE-MICT (r8205), Geosci. Model Dev., 17, 4727–4754, https://doi.org/10.5194/gmd-17-4727-2024, 2024.

Yu, Z., Loisel, J., Brosseau, D. P., Beilman, D. W., and Hunt, S. J.: Global peatland dynamics since the Last Glacial Maximum, Geophys. Res. Lett., 37, L13402, https://doi.org/10.1029/2010gl043584, 2010.

Zhu, D., Peng, S. S., Ciais, P., Viovy, N., Druel, A., Kageyama, M., Krinner, G., Peylin, P., Oettle, C., Piao, S. L., Poulter, B., Schepaschenko, D., and Shvidenko, A.: Improving the dynamics of Northern Hemisphere high-latitude vegetation in the ORCHIDEE ecosystem model, Geosci. Model Dev., 8, 2263–2283, https://doi.org/10.5194/gmd-8-2263-2015, 2015.
Zhu, D., Peng, S., Ciais, P., Zech, R., Krinner, G., Zimov, S., and Grosse, G.: Simulating soil organic carbon in yedoma deposits during the Last Glacial Maximum in a land surface model, Geophys. Res. Lett., 43, 5133–5142, https://doi.org/10.1002/2016GL068874, 2016.

Zhu, D., Ciais, P., Krinner, G., Maignan, F., Jornet Puig, A., and Hugelius, G.: Controls of soil organic matter on soil thermal dynamics in the northern high latitudes, Nat. Commun., 10, 3172, https://doi.org/10.1038/s41467-019-11103-1, 2019.