IPSL-Perm-LandN: improving the IPSL Earth System Model to represent permafrost carbon-nitrogen interactions

Permafrost soils have the potential to release large amounts of soil carbon to the atmosphere under climate change. However, in the Sixth Coupled Model Intercomparison Project (CMIP6), only two Earth System Models (ESM) represented permafrost carbon, both sharing the same land surface model. This makes future permafrost carbon dynamics highly uncertain and underscores the urgent need to include permafrost carbon in ESMs to enable more reliable future projections of climate change and remaining carbon budget estimates. Here, we present IPSL-Perm-LandN, an improved version of the Institut Pierre-Simon Laplace (IPSL) ESM (used for CMIP6) aiming at better representing high-latitude land ecosystems. The main developments are the inclusion of an explicit nitrogen cycle and of key permafrost physical and biogeochemical processes. The latent heat associated with soil water freeze/thaw is taken into account in the energy budget, as well as soil thermal insulation by soil organic matter and a surface organic layer (e.g. litter or moss). Soil organic carbon and nitrogen are vertically resolved with depth-dependent decomposition dynamics, a key feature for representing the effect of gradual permafrost thaw on soil biogeochemistry. Cryoturbation is represented as a diffusion process that buries organic matter in the deeper soil layers. Compared to the previous version of the model used for CMIP6, we show that the extent of the permafrost region has improved significantly and that the simulated active layer thickness in the Arctic is in better agreement with observations. Permafrost soil carbon stocks have increased 20-fold to reach 1006 PgC in the top 3 m of soil, which is consistent with observation-based estimates. We simulate that the permafrost region has been a net carbon sink over the past 150 years (+0.32 ± 0.04 PgC yr⁻¹ on average between 2005 and 2014), primarily due to carbon uptake from boreal forests. This is comparable with recent pan-Arctic carbon balance estimates, when accounting for unrepresented processes in our model (fire and riverine carbon losses). Overall, the inclusion of permafrost processes has improved the response of the model to anthropogenic perturbations in high latitudes over the past century, marking a step forward in the representation of Arctic ecosystems.

1 Introduction

The permafrost region, located mainly in cold high-latitude areas, is home to complex interactions between physical and biogeochemical processes. It contains large amounts of thermally protected soil organic carbon that has accumulated over millennia (IPCC SROCCC, 2019; Hugelius et al., 2014). Anthropogenic greenhouse gas emissions and the resulting climate warming lead to permafrost thawing, which threatens these vulnerable carbon stocks (Smith et al., 2022; Burke et al., 2020; Biskaborn et al., 2019). Subsequent decomposition of the newly unfrozen permafrost carbon would lead to CO₂ and CH₄ emissions, further amplifying global warming in a positive feedback loop known as the permafrost carbon-climate feedback (Schuur et al., 2015; Schaefer et al., 2014). On the other hand, increased CO₂ fertilisation from rising atmospheric CO₂ concentrations and longer growing seasons caused by warming could increase vegetation productivity in negative feedback loops, partially offsetting the positive climate feedback from warming-induced soil carbon losses (Gier et al., 2024; Arora et al., 2020). Nitrogen also impacts carbon cycle feedbacks in both directions. It can reduce vegetation productivity through nitrogen limitation (positive feedback, Gier et al., 2024; Street and Caldararu, 2022; Davies-Barnard et al., 2020), but can also increase plant carbon uptake through increased soil nitrogen availability due to soil warming and permafrost thaw (negative feedback Burke et al., 2022; Salmon et al., 2016; Finger et al., 2016; Koven et al., 2011). However, the timing and magnitude of these feedbacks remain highly uncertain (Schuur et al., 2022; Schädel et al., 2018). Therefore, the resulting overall response of the carbon cycle to anthropogenic emissions in permafrost regions is a major unknown in future projections of the global carbon cycle (Kleinen and Brovkin, 2018; McGuire et al., 2018; MacDougall et al., 2012).

Earth system models (ESMs) are numerical representations of the Earth system that simulate the coupled dynamics and exchanges of energy, water and carbon between the atmosphere, the ocean and continental surfaces. Based on the representation of physical and biogeochemical mechanisms at a large range of scales, they are essential tools for studying the past, present and future dynamics of the Earth's climate and carbon cycle. In particular, their use for climate projections plays a key role in informing adaptation and mitigation policies and is at the basis for IPCC Assessment Reports. Compared to simpler models, they take into account the feedbacks between the processes that control the exchange of energy, water and carbon, and are the most comprehensive representation of the Earth system currently available. They can be driven by different socio-economic and greenhouse gas emission-related scenarios to explore possible futures, and can isolate individual feedbacks to quantify their contribution to the global response (e.g. Arora et al., 2020). ESMs are therefore particularly well suited to studying the future dynamics of the permafrost carbon cycle as they provide a mechanistic description of the complex interactions between climate and the carbon cycle.

However, despite the urgent need to accurately predict the future permafrost carbon dynamics, the physical and biogeochemical mechanisms of permafrost are still not well represented in ESMs (Matthes et al., 2025; Schädel et al., 2024; Natali et al., 2021; Burke et al., 2020; Slater and Lawrence, 2013). Reducing the uncertainties surrounding permafrost carbon cycle feedbacks is becoming especially important as ESMs move towards emission-driven simulations, in which the atmospheric CO₂ concentrations will be largely determined by the simulated carbon cycle dynamics (Steinert and Sanderson, 2025; Park et al., 2025; Sanderson et al., 2024). Such emission-driven simulations are particularly relevant for producing policy-oriented climate projections and for properly accounting for feedbacks between the carbon cycle and climate. Although efforts have been made to include physical permafrost processes in land surface models (LSMs, the land component of ESMs), including soil freeze/thaw cycles and the influence of hydrology on soil thermal properties (De Vrese et al., 2023; Steinert et al., 2021; Cuntz and Haverd, 2018; Hagemann et al., 2016; Ekici et al., 2014; Gouttevin et al., 2012), multilayer snow schemes including snow hydrological and thermal effects and snow compaction (Decharme et al., 2016; Wang et al., 2013; Ekici et al., 2014), or soil organic matter and moss insulation (Yokohata et al., 2020; Guimberteau et al., 2018; Chadburn et al., 2015; Lawrence et al., 2008), large discrepancies remain between models. Most of the CMIP6 ESMs perform poorly in simulating critical permafrost properties such as the active layer thickness (ALT, the maximum annual thaw depth) or snow insulation, partly due to shallow and poorly resolved soil profiles (Burke et al., 2020).

Furthermore, the representation of the permafrost carbon cycle in ESMs is still in its infancy. Among the CMIP6 models, only two ESMs (CESM2 and NorESM2-LR) included a vertically resolved representation of soil carbon – an essential feature for simulating permafrost carbon dynamics – and both shared the same land surface model (CLM5) (Schädel et al., 2024). The lack of such a vertical soil carbon discretisation prevents most models from representing the large soil carbon content of the permafrost region as well as the effect of gradual and abrupt (e.g. through fire or thaw slumps) permafrost thaw on soil carbon dynamics and the permafrost carbon-climate feedback (Schädel et al., 2024; Gier et al., 2024; Varney et al., 2022; Turetsky et al., 2020). Therefore, most models used for the calculation of remaining carbon budgets do not include permafrost carbon and the permafrost contribution must be added from external estimates (Rogelj et al., 2019). The inclusion of nitrogen processes in ESMs and their coupling to the carbon cycle has been a major advance in the last decade, although only half of the CMIP6 ESMs representing the carbon cycle had an explicit representation of the nitrogen cycle (six out of the eleven ESMs from Arora et al., 2020). An accurate representation of the nitrogen cycle is particularly important for high latitudes where vegetation is generally considered to be nitrogen-limited and where mineral nitrogen release from permafrost thaw could affect both vegetation productivity and soil organic carbon decomposition (Street and Caldararu, 2022; Wooliver et al., 2019; Beermann et al., 2017; Keuper et al., 2017). The complex interactions between carbon and nitrogen in permafrost regions could lead to very different model responses and their inclusion in ESMs is therefore key to evaluating and reducing uncertainties in future projections of permafrost carbon dynamics (Burke et al., 2022; Lacroix et al., 2022; Koven et al., 2015a).

This paper describes and evaluates a new version of the IPSL Earth system model – called IPSL-Perm-LandN – designed to better simulate high-latitude processes and permafrost carbon dynamics, based on the CMIP6 version IPSL-CM6A-LR (Boucher et al., 2020). New developments include vertically resolved coupled carbon and nitrogen cycles and key physical and biogeochemical permafrost processes in ORCHIDEE, the land surface component of the model (Vuichard et al., 2019; Guimberteau et al., 2018; Krinner et al., 2005). In particular, the model accounts for nitrogen limitation of vegetation photosynthetic activity and decomposition of soil carbon and litter (Vuichard et al., 2019). It represents permafrost freeze/thaw cycles (based on Gouttevin et al., 2012), soil insulation by snow, soil organic matter and surface organic layers (e.g. litter, moss,  Gaillard et al., 2025b), vertically resolved soil organic carbon and nitrogen with depth-dependent dynamics, thermal protection of soil organic matter when frozen and its mixing along the vertical profile (bio- and cryoturbation).

IPSL-Perm-LandN marks an important step in the representation of high-latitude ecosystems in the IPSL ESM by integrating first-order permafrost processes. These new developments significantly improve the simulation of permafrost physics and carbon cycle dynamics in the IPSL ESM. It is expected to be continuously improved by integrating new mechanisms (e.g. fire/permafrost interactions or abrupt thaw) and by better constraining the processes already included.

2 Model description

2.1 General presentation

IPSL-Perm-LandN is based on IPSL-CM6A-LR, the version of the Earth system model developed by the Institut Pierre-Simon Laplace (IPSL) modeling center for the 6th phase of the Coupled Model Intercomparison Project (CMIP6) (Boucher et al., 2020; Lurton et al., 2020; Eyring et al., 2016). It is composed of the atmospheric model LMDZ (version 6A-LR) (Hourdin et al., 2020), the oceanic model NEMO and the land surface model ORCHIDEE. The ocean model includes the ocean physics NEMO-OPA (Madec et al., 2016), the sea ice dynamics and thermodynamics NEMO-LIM3 (Rousset et al., 2015; Vancoppenolle et al., 2009) and the ocean biogeochemistry NEMO-PISCES (Aumont et al., 2015) models. The coupling between the atmosphere and the surface is done every 15 min while the other components of IPSL-Perm-LandN are coupled at a frequency of 90 min. The resolution of the atmospheric model is 144 × 143 points in longitude and latitude, corresponding to a resolution of 2.5° × 1.3° (average resolution of 157 km), and 79 vertical levels extending up to 80 km. The resolution of the ocean model is 1° and 75 vertical levels.

This new configuration of the IPSL Earth System Model aims to better represent high-latitude ecosystems and climate as well as permafrost physics and carbon cycle. The main modifications compared to IPSL-CM6A-LR concern the land surface model ORCHIDEE. While IPSL-CM6A-LR included ORCHIDEE-v2, a carbon-only version of the land component, IPSL-Perm-LandN uses ORCHIDEE-v3 which includes the implementation of a fully prognostic nitrogen cycle (Vuichard et al., 2019) and several key permafrost physical and biogeochemical processes (Gaillard et al., 2025b; Zhu et al., 2019; Guimberteau et al., 2018).

Section 2.2 and 2.3 briefly recall the main characteristics of the atmosphere and ocean components. A more complete description can be found in Boucher et al. (2020).

2.2 Atmospheric model LMDZ

The atmospheric general circulation model used in IPSL-Perm-LandN is LMDZ6A-LR (Hourdin et al., 2020). It solves the primitive equations using a finite-difference formulation (Sadourny and Laval, 1984), and advects water vapour, solid and liquid water and trace gases with a monotonic second-order finite volume scheme (Hourdin and Armengaud, 1999; Van Leer, 1977). LMDZ6A-LR physical parameterisations are based on LMDZ5B (Hourdin et al., 2013), the version of LMDZ included in IPSL-CM5B that participated in CMIP5. The turbulent scheme is based on the turbulent kinetic energy prognostic equation of Yamada (1983), a thermal plume model (Hourdin et al., 2002; Rio and Hourdin, 2008) and a parameterization of cold pools (Grandpeix and Lafore, 2010; Grandpeix et al., 2010). Convection has been improved since LMDZ5B with a better representation of the transition from stratocumulus to cumulus clouds (Hourdin et al., 2019) and the inclusion of a statistical triggering for deep convection (Rochetin et al., 2014a, b). The radiative transfer scheme includes the Rapid Radiative Transfer Model (RRTM) for thermal infrared radiation and a six-bands versions of Fouquart and Bonnel (1980) scheme for solar radiation. Gravity waves generated by mountains, convection (Lott and Guez, 2013) and fronts (de la Cámara et al., 2016; de la Cámara and Lott, 2015) are represented, as well as the quasi-biennal oscillation. Further details on the LMDZ6A model can be found in Hourdin et al. (2020).

2.3 Ocean model NEMO

The version 3.6 of NEMO (Nucleus for European Models of the Ocean) is the ocean component of IPSL-Perm-LandN and includes both physical and biogeochemical processes. The ocean physics are represented by NEMO-OPA (Madec et al., 2016) and are based on the Navier-Stokes equations and a nonlinear equation of state (Roquet et al., 2015). The vertical mixing of momentum and tracers uses a turbulent energy scheme (Blanke and Delecluse, 1993; Gaspar et al., 1990) and parameterisations of mixing caused by internal tides (de Lavergne et al., 2019; de Lavergne, 2016) and submesoscale processes (Fox-Kemper et al., 2011).

Sea ice is described by the NEMO-LIM (version 3.6) model (Rousset et al., 2015; Vancoppenolle et al., 2009). NEMO-LIM uses a distribution of ice thickness (Bitz et al., 2001; Lipscomb, 2001), allowing the representation of thin to thick ice. Sea ice can be transported horizontally and snow can accumulate above it. Vertically, two ice layers and one snow layer are represented. Within the ice layers, the ice is represented by an elastic-viscous plastic continuum (Bouillon et al., 2013; Hunke and Dukowicz, 1997). It can dynamically exchange energy and salinity with the ocean, allowing for a prognostic evolution of the coupled system. Notably, ice albedo parameters are used for model tuning as well as the snow thermal conductivity (Boucher et al., 2020).

The ocean biogeochemistry is based on PISCES-v2 (Aumont et al., 2015) and simulates the lower trophic levels of marine ecosystems, including phytoplankton and zooplankton, and the biogeochemical cycles of carbon and main nutrients (phosphorus, nitrogen, silicon and iron). The carbon cycle includes a representation of carbonate chemistry. Nutrients are supplied to the ocean by atmospheric deposition, river inputs and sediment mobilisation. Carbon compounds can be exchanged with the atmosphere through physical and biogeochemical processes, and buried at the bottom of the ocean. The parameterisation of nitrogen fixation has been modified compared to IPSL-CM6A-LR, which has an impact on the biological carbon pump at high temperatures.

2.4 Land surface model ORCHIDEE

2.4.1 General description

ORCHIDEE-v3 (ORganizing Carbon and Hydrology in Dynamic EcosystEms) is a state-of-the art process-based land surface model that calculates energy, water, carbon and nitrogen exchanges between the surface and the atmosphere, as well as terrestrial physical and biogeochemical processes. It is composed of two main sub-models : SECHIBA that describes exchanges of energy and water between the atmosphere, the biosphere and the soil, and STOMATE that simulates the phenology and carbon and nitrogen dynamics of the terrestrial biosphere (Vuichard et al., 2019; Zaehle and Friend, 2010; Krinner et al., 2005). Fast processes (e.g. latent and sensible heat fluxes, photosynthesis, ecosystem respiration) are computed every 15 min while slow processes (e.g. carbon and nitrogen allocation) are computed daily.

Vegetation is represented by plant functional types (PFTs), i.e. groups of species sharing similar characteristics (Prentice et al., 1992). These PFTs share the same equations for most processes, but with different parameters. ORCHIDEE-v3 represents 15 PFTs, classified into forests, grasses, crops and bare soil, describing a variety of ecosystems (Table A1). PFTs can coexist in every grid box and the fraction occupied by each PFT is read from a prescribed map (which can change on a yearly basis) (Lurton et al., 2020). For each PFT, carbon and nitrogen are contained in seven plant pools (leaves, below- and above-ground sapwood and heartwood, fruits and fine roots), five litter pools (above- and below-ground metabolic and structural, and woody litter) and three soil pools (active, slow, and passive).

ORCHIDEE-v3 represents energy exchanges between the surface and the atmosphere and takes into account shortwave and longwave radiative fluxes, turbulent latent and sensible heat fluxes, and a ground flux (Ducoudré et al., 1993). The turbulent fluxes are calculated separately for each PFT and then summed for each grid box. This coupling with the atmosphere is regulated by vegetation properties such as its albedo and its height (which impacts on surface roughness). Within the ground, heat transfers are represented by a heat diffusion equation and depend on the mineral and organic soil properties (thermal capacity, thermal conductivity, porosity) and soil hydrology. Mineral soil properties are extrapolated from the soil texture map of Zobler (1986). Soil thermal dynamics is based on an 18-layer vertical scheme, extending down to 90 m (Table A2). The thickness of each layer increases with depth, with thinner layers near the surface. A zero flux condition is imposed at the bottom boundary.

The model also represents exchanges of water between the surface and the atmosphere. Water reaches the land through rain or snowfall, and can be lost through evaporation of water stored in the soil but also intercepted by the canopy, transpiration by vegetation, snow sublimation, surface runoff and percolation and transfer to groundwater (i.e. drainage). Internal water exchanges between land components can also occur through various mechanisms, such as snow melt, or plant root uptake. Soil moisture is resolved on a 11-layer scheme (the same as for soil thermics) down to 2 m, where a free drainage bottom boundary condition is imposed (de Rosnay et al., 2002). Therefore, the bedrock differs between soil thermics (90 m) and hydrology (2 m), and a deeper hydraulic scheme is under development. Water is transferred from one layer to another according to a one-dimensional Fokker-Planck equation (Ducharne et al., 2018). Below 2 m, the calculation of soil thermal properties uses the water content of the deepest hydrological layer. Vegetation has a major influence on water exchanges by regulating evapotranspiration through stomatal closure and soil water uptake.

The representation of the carbon and nitrogen cycles have already been described in detail in Vuichard et al. (2019), Zaehle and Friend (2010) and Krinner et al. (2005). The following sections are limited to the description of relevant processes for high latitudes and new developments. A more detailed description of ORCHIDEE-v3 can be found in Appendix A.

2.4.2 Latent heat of soil water phase change

The improvements to permafrost physics (Sect. 2.4.2, 2.4.3 and 2.4.4) have been described in Gaillard et al. (2025b) and are summarised here for the sake of completeness. The ground temperature in ORCHIDEE-v3 is calculated using a one-dimensional Fourier equation with a boundary condition at the surface allowing heat exchanges with the atmosphere (Eq. 5 in Gouttevin et al., 2012):

$$\begin{array}{}\text{(1)}& c_{\mathrm{app}}{\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{th}}{\displaystyle \frac{\partial T}{\partial z}}\right)\end{array}$$

where T is the soil temperature (K) and K_th the soil thermal conductivity (W m⁻¹ K⁻¹). c_app is apparent volumetric soil thermal capacity (J K⁻¹ m⁻³). It incorporates volumetric soil thermal capacity and a term representing the latent heat of soil water phase changes during melting and freezing:

$$\begin{array}{}\text{(2)}& c_{\mathrm{app}}=c_p-{\mathit{\rho }}_{\mathrm{ice}}L{\displaystyle \frac{\mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{ice}}}{\mathrm{\Delta }T}}\end{array}$$

where c_p is the volumetric soil thermal capacity (J K⁻¹ m⁻³), ρ_ice the ice density (kg m⁻³), L the latent heat of fusion (J kg⁻¹) and Θ_ice the volumetric ice content (m³ m⁻³).

Taking into account the latent heat of water phase change is essential to correctly simulate the soil thermal dynamics in the permafrost region. It acts as a buffer – also called zero-curtain effect – absorbing energy from thawing ice in spring and summer and releasing energy when the water refreezes in autumn and winter, thus reducing the amplitude of the seasonal cycle of ground temperature.

2.4.3 Modifications of soil thermal properties by soil organic carbon

Soil organic carbon (SOC) has been shown to be an important driver of surface-atmosphere energy exchanges at high latitudes and of permafrost thermal dynamics (Zhu et al., 2019; Loranty et al., 2018). Its effect is taken into account in our model by weighting the soil thermal properties by the SOC volume fraction (f_SOC). f_SOC is calculated as:

$$\begin{array}{}\text{(3)}& f_{\mathrm{SOC}}={\displaystyle \frac{C_{\mathrm{SOC}}}{C_{\mathrm{SOC}\max }}}\end{array}$$

where C_SOC is the SOC density (kgC m⁻³) and C_SOC max= 500 kgC m⁻³ is a reference value. C_SOC max has been tuned to simulate a realistic high latitude climate (Gaillard et al., 2025b), ensuring that its value remains in the range of soil carbon densities from the SoilGrids database (Poggio et al., 2021; Batjes et al., 2019). The heat diffusion equation (Eq. 1) then uses the total soil thermal conductivity and capacity (mixing mineral and organic soil properties).

Solid and dry soil thermal conductivities and the dry thermal capacity are computed as weighted averages of those of mineral and organic soils (Guimberteau et al., 2018):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{solid}}=(\mathrm{1}-f_{\mathrm{SOC}}){\mathit{\lambda }}_{\mathrm{solid}\mathrm{mineral}}+f_{\mathrm{SOC}}{\mathit{\lambda }}_{\mathrm{solid}\mathrm{SOC}}\text{(5)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{dry}}=(\mathrm{1}-f_{\mathrm{SOC}}){\mathit{\lambda }}_{\mathrm{dry}\mathrm{mineral}}+f_{\mathrm{SOC}}{\mathit{\lambda }}_{\mathrm{dry}\mathrm{SOC}}\text{(6)}& {\displaystyle }{c}_{\mathrm{dry}}=(\mathrm{1}-f_{\mathrm{SOC}})c_{\mathrm{mineral}}+f_{\mathrm{SOC}}{c}_{\mathrm{SOC}}\end{array}$$

where c_mineral (J K⁻¹ m⁻³) and λ_mineral (W m⁻¹ K⁻¹) are the thermal capacities and conductivities of solid/dry mineral soils, which depend on the dominant soil texture of the grid box. Solid refers to the solid fraction of the soil (excluding pores) while the dry fraction also includes the pores filled with air (not those filled with water). The total thermal capacity is then calculated for each soil layer as:

$$\begin{array}{}\text{(7)}& c=c_{\mathrm{dry}}+{\mathrm{\Theta }}_{\mathrm{liq}}{c}_{\mathrm{liq}}+{\mathrm{\Theta }}_{\mathrm{ice}}{c}_{\mathrm{ice}}\end{array}$$

where Θ_liq (unitless) and Θ_ice (unitless) are the volumetric liquid water and ice contents computed by the model and c_liq and c_ice are the thermal capacities (J K⁻¹ m⁻³) of liquid water and ice, respectively equal to 4.18 × 10⁶ and 2.11 × 10⁶ J K⁻¹ m⁻³. The thermal conductivity of dry organic carbon (c_dry) is fixed at 2.5 × 10⁶ J K⁻¹ m⁻³. For each soil layer, the thermal conductivity is computed as:

$$\begin{array}{}\text{(8)}& \mathit{\lambda }=\mathit{\text{Ke}}{\mathit{\lambda }}_{\mathrm{sat}}+(\mathrm{1}-\mathit{\text{Ke}}){\mathit{\lambda }}_{\mathrm{dry}}\end{array}$$

where:

$$\begin{array}{}\text{(9)}& {\mathit{\lambda }}_{\mathrm{sat}}={\mathit{\lambda }}_{\mathrm{solid}}^{(\mathrm{1}-{\mathrm{\Theta }}_{\mathrm{sat}})}{\mathit{\lambda }}_{\mathrm{liq}}^{\left({\mathrm{\Theta }}_{\mathrm{sat}}\frac{{\mathrm{\Theta }}_{\mathrm{liq}}}{{\mathrm{\Theta }}_{\mathrm{liq}}+{\mathrm{\Theta }}_{\mathrm{ice}}}\right)}{\mathit{\lambda }}_{\mathrm{ice}}^{\left({\mathrm{\Theta }}_{\mathrm{sat}}\frac{{\mathrm{\Theta }}_{\mathrm{ice}}}{{\mathrm{\Theta }}_{\mathrm{liq}}+{\mathrm{\Theta }}_{\mathrm{ice}}}\right)}\end{array}$$

with λ_liq and λ_ice the thermal conductivities of liquid water and ice, respectively equal to 0.57 and 2.2 W m⁻¹ K⁻¹, and Θ_sat (unitless) the volumetric moisture content at saturation, which depends on the dominant mineral soil texture. The thermal conductivity of dry organic carbon (λ_dry) is fixed at 0.25 W m⁻¹ K⁻¹.

Ke is the Kersten number defined for unfrozen soil as:

$$\begin{array}{}\text{(10)}& \mathit{\text{Ke}}=\left\{\begin{array}{ll}{\log }_{\mathrm{10}}\left(S_{\mathrm{r}}\right)+\mathrm{1}& \mathrm{if}{S}_{\mathrm{r}}>\mathrm{0.1}\\ \mathrm{0.7}{\log }_{\mathrm{10}}\left(S_{\mathrm{r}}\right)+\mathrm{1}& \mathrm{if}\mathrm{0.05}<S_{\mathrm{r}}\le \mathrm{0.1}\\ \mathrm{0}& \mathrm{if}{S}_{\mathrm{r}}\le \mathrm{0.05}\end{array}\right.\end{array}$$

where S_r is the saturation ratio and is calculated as $S_{\mathrm{r}}=\frac{{\mathrm{\Theta }}_{\mathrm{liq}}}{{\mathrm{\Theta }}_{\mathrm{sat}}}$. For (fully or partially) frozen soils, Ke =S_r.

The modification of soil thermal parameters by soil organic carbon creates a coupling between the carbon cycle and soil thermodynamics, eventually impacting surface-atmosphere energy transfers. Importantly, the porosity calculated by the thermal module of ORCHIDEE-v3 differs from that used in the hydrological scheme (which is equal to that of a mineral soil), which prevents a direct feedback between soil moisture and soil temperature through soil porosity.

2.4.4 Modification of soil thermal properties by a surface organic layer

In the Arctic, the surface organic layer (SOL) formed by litter and groundcover vegetation (moss, lichens) may significantly reduce surface-atmosphere energy exchanges through their insulative properties and therefore thermally protect permafrost soils from warmer summer air temperatures (Loranty et al., 2018; Porada et al., 2016). In IPSL-Perm-LandN, we decided to modify the thermal capacity and conductivity of the upper soil layers to mimic the effect of such a surface organic layer on soil thermal dynamics. This assumption is made in some land surface models (Wu et al., 2016; Chadburn et al., 2015), whereas some other models explicitly represent an organic layer on top of the soil column (Park et al., 2018; Porada et al., 2016). We further assumed that the surface organic layer covers a fraction f_SOL of each grid box containing boreal PFTs, as bryophytes are widespread in these ecosystems (Lewis et al., 2017; Barry et al., 2013).

The calculation of the effect of the surface organic layer on soil thermal transfers is carried out in two steps. First, a virtual column (not explicitly represented in the model) is defined over a fraction f_SOL of the grid box, representing moss, lichen and/or decomposing litter (dashed red in Fig. A1). The thermal capacity of the virtual column is calculated as a weighted average of the surface organic layer and soil thermal capacities:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& c_{\mathrm{virtual}\mathrm{column}}\mathrm{SID}=c_{\mathrm{SOL}}\mathrm{SOLT}+c_{\mathrm{soil}}\mathrm{SID}\\ & \iff c_{\mathrm{virtual}\mathrm{column}}=c_{\mathrm{SOL}}{\displaystyle \frac{\mathrm{SOLT}}{\mathrm{SID}}}+c_{\mathrm{soil}}\end{array}\end{array}$$

where c_SOL is the volumetric thermal capacity of the surface organic layer (J K⁻¹ m⁻³), c_soil is the volumetric soil thermal capacity (as calculated in Sect. 2.4.3, J K⁻¹ m⁻³), f_SOL is the fraction of the grid box that contains the surface organic layer, SOLT is the surface organic layer thickness and SID is the soil integration depth, i.e. the depth down to which the properties of the soil organic layer are mixed with those of the soil.

Then, the total thermal capacity of the grid box (c_tot), which takes into account the fraction not covered by the surface organic layer, is calculated as the weighted average of c_{virtual column} and c_soil:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{c}_{\mathrm{tot}}& =f_{\mathrm{SOL}}{c}_{\mathrm{virtual}\mathrm{column}}+(\mathrm{1}-f_{\mathrm{SOL}})c_{\mathrm{soil}}\\ & =f_{\mathrm{SOL}}{c}_{\mathrm{SOL}}{\displaystyle \frac{\mathrm{SOLT}}{\mathrm{SID}}}+c_{\mathrm{soil}}\end{array}\end{array}$$

The approach for thermal conductivity is similar but takes into account that it is an intensive property (i.e. its value is independent of the size of the system). The thermal conductivity virtual column (λ_{virtual column}) is the equivalent thermal conductivity of the surface organic layer and soil layers in series:

$$\begin{array}{}\text{(13)}& {\mathit{\lambda }}_{\mathrm{virtual}\mathrm{column}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{SOL}}{\mathit{\lambda }}_{\mathrm{soil}}\mathrm{SID}}{\mathrm{SOLT}{\mathit{\lambda }}_{\mathrm{soil}}+\mathrm{SID}{\mathit{\lambda }}_{\mathrm{SOL}}}}\end{array}$$

where λ_SOL is the thermal conductivity of the soil organic layer and λ_soil is the thermal conductivity of the soil.

The total thermal conductivity of the grid box is the equivalent thermal conductivity of the surface organic layer column and the soil column in parallel:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\mathit{\lambda }}_{\mathrm{tot}}=f_{\mathrm{SOL}}{\mathit{\lambda }}_{\mathrm{virtual}\mathrm{column}}+(\mathrm{1}-f_{\mathrm{SOL}}){\mathit{\lambda }}_{\mathrm{soil}}\\ & =f_{\mathrm{SOL}}{\mathit{\lambda }}_{\mathrm{SOL}}{\mathit{\lambda }}_{\mathrm{soil}}{\displaystyle \frac{\mathrm{SID}}{\mathrm{SOLT}{\mathit{\lambda }}_{\mathrm{soil}}+\mathrm{SID}{\mathit{\lambda }}_{\mathrm{SOL}}}}+(\mathrm{1}-f_{\mathrm{SOL}}){\mathit{\lambda }}_{\mathrm{soil}}\end{array}\end{array}$$

Finally, the mineral soil capacity and conductivity are replaced by c_tot and λ_tot in all the soil layers between the surface and SID.

In this study, we chose f_SOL=1, SOLT =0.03 m and SID =0.03 m for evaluating the model. This value of SOLT is consistent with the moss thickness measured in Soudzilovskaia et al. (2013). SID was chosen small enough to allow the soil organic layer to influence surface-atmosphere energy exchanges, but to limit the modification of soil thermal properties to the very top layers.

In addition, the thermal properties of the surface organic layer depend on its water content (Soudzilovskaia et al., 2013; O'Donnell et al., 2009). They are parameterized using observations made on mosses, using the upper soil water content of each soil layer down to SID as a proxy for the water content of the surface organic layer. The thermal capacity of the soil organic layer is calculated as:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& c_{\mathrm{SOL}}=\\ & \left\{\begin{array}{l}{c}_{\mathrm{SOL}\mathrm{dry}}+\mathit{\theta }(c_{\mathrm{SOL}\mathrm{wet}}-c_{\mathrm{SOL}\mathrm{dry}})\\ \mathrm{if}T<\mathrm{2}\mathit{°}C\\ c_{\mathrm{SOL}\mathrm{dry}}+\mathit{\theta }\left[\right(\frac{T}{\mathrm{2}}+\mathrm{1})c_{\mathrm{SOL}\mathrm{wet}}-\frac{T}{\mathrm{2}}{c}_{\mathrm{SOL}\mathrm{frozen}}-c_{\mathrm{SOL}\mathrm{dry}}]\\ \mathrm{if}-\mathrm{2}\mathit{°}C\le T\le \mathrm{0}\mathit{°}C\\ c_{\mathrm{SOL}\mathrm{dry}}+\mathit{\theta }(c_{\mathrm{SOL}\mathrm{frozen}}-c_{\mathrm{SOL}\mathrm{dry}})\\ \mathrm{if}T>\mathrm{0}\mathit{°}C\end{array}\right.\end{array}\end{array}$$

where c_SOL dry, c_SOL wet and c_SOL frozen (J m⁻³ K⁻¹) are the thermal capacities of dry, wet and frozen surface organic layers, respectively, and θ is the volumetric moisture content (unitless).

The thermal conductivity of the soil organic layer is calculated as:

$$\begin{array}{}\text{(16)}& {\mathit{\lambda }}_{\mathrm{SOL}}={\mathit{\lambda }}_{\mathrm{SOL}\mathrm{dry}}+\mathit{\theta }({\mathit{\lambda }}_{\mathrm{SOL}\mathrm{sat}}-{\mathit{\lambda }}_{\mathrm{SOL}\mathrm{dry}})\end{array}$$

where λ_SOL dry is the thermal conductivity of a dry surface organic layer and λ_SOL sat is the thermal conductivity of a saturated surface organic layer, calculated as:

$$\begin{array}{}\text{(17)}& {\mathit{\lambda }}_{\mathrm{SOL}\mathrm{sat}}={\mathit{\lambda }}_{\mathrm{SOL}\mathrm{liq}}^{\left({\mathit{\theta }}_{\mathrm{sat}}\frac{{\mathit{\theta }}_{\mathrm{liq}}}{{\mathit{\theta }}_{\mathrm{liq}}+{\mathit{\theta }}_{\mathrm{ice}}}\right)}{\mathit{\lambda }}_{\mathrm{SOL}\mathrm{frozen}}^{\left({\mathit{\theta }}_{\mathrm{sat}}\frac{{\mathit{\theta }}_{\mathrm{ice}}}{{\mathit{\theta }}_{\mathrm{liq}}+{\mathit{\theta }}_{\mathrm{ice}}}\right)}\end{array}$$

The values of surface organic layer thermal properties are taken from in situ measurements and laboratory experiments (Soudzilovskaia et al., 2013; O'Donnell et al., 2009). Thermal capacities are set to c_SOL dry=0.29×10⁶, $c_{\mathrm{SOL}\mathrm{wet}}=\mathrm{4.29}\times {\mathrm{10}}^{\mathrm{6}}$ and $c_{\mathrm{SOL}\mathrm{frozen}}=\mathrm{3.26}\times {\mathrm{10}}^{\mathrm{6}}$ J m⁻³ K⁻¹ (Soudzilovskaia et al., 2013; Druel et al., 2017). Thermal conductivities are equal to λ_SOL dry= 0.05, λ_SOL wet= 0.56 and λ_SOL frozen= 1.40 W m⁻¹ K⁻¹ (O'Donnell et al., 2009; Porada et al., 2016).

2.4.5 Snow

ORCHIDEE-v3 uses a 3-layer snow scheme of intermediate complexity with dynamic layer thickness, which was already used in IPSL-CM6A-LR. Snow strongly influences the surface-atmosphere energy transfer at high latitudes due to its insulating properties. Heat diffusion within the snowpack is accounted for by a heat-transfer equation:

$$\begin{array}{}\text{(18)}& c_p{\displaystyle \frac{\partial T_j}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_j{\displaystyle \frac{\partial T_j}{\partial z}}\right)+{\displaystyle \frac{\partial R}{\partial z}}\end{array}$$

where T_j is the snow temperature of the layer j, c_p is the snow heat capacity (J K⁻¹ m⁻³), κ_j is the thermal conductivity of the snow (W m⁻¹ K⁻¹) and takes into account vapour transfer in the snow, z is the vertical coordinate and t is time. $\frac{\partial R}{\partial t}$ is the solar-radiative energy source and depends on the incoming solar radiative energy and the snow depth.

Water phase change can occur within the snowpack as snow melts or refreezes, further affecting soil hydrology and surface-atmosphere water exchange. In particular, snow can melt in the upper layer of the snowpack due to solar radiation, infiltrate down to the next layer and may refreeze, releasing latent heat and heating lower layers. Snow compaction is also represented and depends on the weight of the overlying snow. It modifies the density and thickness of snow layers over time. Finally, the snow albedo is included and depends on the snowfall rate and the liquid water content of the snowpack.

Further details on these processes and their implementation can be found in Wang et al. (2013).

2.4.6 Soil carbon and nitrogen dynamics

Soil organic carbon and nitrogen dynamics in ORCHIDEE follow a CENTURY-based scheme (Parton et al., 1993) which is schematised in Figs. A2 and A3. Plant residues are divided into structural and metabolic litter pools according to their lignin content. Litter decomposition follows a first-order kinetics with pool-dependent decomposition factors, and depends on temperature, moisture and lignin content. Part of the decomposed carbon is respired as CO₂ and the remaining flux is transferred to soil organic carbon (SOC) pools. Importantly, the model only represents CO₂ emissions and does not include CH₄ dynamics. Active, slow and passive SOC pools have different turnover times and can exchange carbon with each other, each time with an associated loss of CO₂ through microbial respiration. SOC decomposition also follows a first-order kinetics with a dependence on soil temperature, moisture and texture (i.e. soil sand, silt and clay content):

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\displaystyle \frac{\partial C_i}{\partial t}}{|}_{\mathrm{decomposition}}& =k_i\cdot f\left(T\right)\cdot f\left(\mathrm{moisture}\right)\\ & \cdot f\left(\mathrm{texture}\right)\cdot C_i\end{array}\end{array}$$

where C_i is the carbon content of the pool i (kgC m⁻², where i corresponds to active, slow or passive) and k_i is the decomposition factor (s⁻¹).

Nitrogen is decomposed at the same rate as carbon. Nitrogen fluxes are driven by carbon fluxes and the C:N ratios of the pools (Fig. A3). The nitrogen flux between a pool A and a pool B (kgN m⁻² s⁻¹) is expressed as the product of the corresponding carbon flux and of the N:C ratio of the receiving pool:

$$\begin{array}{}\text{(20)}& f_{\mathrm{nitrogen},\mathrm{A}\mathrm{\to }\mathrm{B}}=f_{\mathrm{carbon},\mathrm{A}\mathrm{\to }\mathrm{B}}\cdot \mathrm{N}:{\mathrm{C}}_{\mathrm{B}}\end{array}$$

The nitrogen associated with the carbon lost by respiration is assumed to be mineralised. If the decomposed organic nitrogen cannot meet the demand of the receiving pools, mineral nitrogen is immobilised to complete the nitrogen flux. If the amount of nitrogen in the mineral pool is not sufficient, nitrogen is taken from the atmosphere to complete the required immobilisation flux. Conversely, if there is an excess of decomposed nitrogen, it is mineralised and transferred to the mineral nitrogen pool. Furthermore, decomposition rates are independent of C:N ratios. These ratios are dynamic and depend on the concentration of soil mineral nitrogen (NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$), with a lower nitrogen demand (higher C:N ratios) when mineral nitrogen is scarce, and a higher demand (lower C:N ratios) when mineral nitrogen stocks are high.

Soil mineral nitrogen follows the DNDC model which accounts for ammonium (NH${}_{\mathrm{4}}^{+}$), nitrates (NO${}_{\mathrm{3}}^{-}$), nitrogen oxides (NO_x) and nitrous oxide (N₂O) (Li et al., 1992, 2000; Zhang et al., 2002). It represents nitrification, denitrification, mineralisation and immobilisation, ammonium adsorption and desorption, plant uptake (NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ only), gaseous emissions and leaching (Fig. A4). Plant uptake is expressed as:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{N}_{\mathrm{up}}& =v_{\max }\times N_{\min }\times \left(k_{N_{\min }}+{\displaystyle \frac{\mathrm{1}}{{\mathrm{K}}_{N_{\min }}+N_{\min }}}\right)\\ & \times f\left({\mathrm{NC}}_{\mathrm{plant}}\right)\times C_{\mathrm{root}}\end{array}\end{array}$$

where N_up is the plant nitrogen uptake (gN m⁻² d⁻¹), N_min is the amount of mineral nitrogen available (NH${}_{\mathrm{4}}^{+}+$ NO${}_{\mathrm{3}}^{-}$, gN m⁻²), v_max is the maximum rate of nitrogen uptake (gN gC⁻¹ d⁻¹), $k_{N_{\min }}$ (m² gN⁻¹) and $K_{N_{\min }}$ (gN m⁻²) are Michaelis–Mentens coefficients, C_root the root carbon mass per unit area (gC m⁻²) and f(NC_plant) the dependency of plant nitrogen uptake to NC_plant, expressed as:

$$\begin{array}{}\text{(22)}& f\left({\mathrm{NC}}_{\mathrm{plant}}\right)=\max \left({\displaystyle \frac{{\mathrm{NC}}_{\mathrm{plant}}-{\mathrm{nc}}_{\mathrm{leaf},\max }}{{\mathrm{nc}}_{\mathrm{leaf},\min }-{\mathrm{nc}}_{\mathrm{leaf},\max }}},\mathrm{0}\right)\end{array}$$

where nc_leaf,min and nc_leaf,max are the minimum and maximum leaf N:C ratios, respectively (PFT-dependent), and NC_plant is defined as the mean N:C ratio of leaves, roots and labile nitrogen pools:

$$\begin{array}{}\text{(23)}& {\mathrm{NC}}_{\mathrm{plant}}={\displaystyle \frac{N_{\mathrm{leaf}}+N_{\mathrm{root}}+N_{\mathrm{labile}}}{C_{\mathrm{leaf}}+C_{\mathrm{root}}+C_{\mathrm{labile}}}}\end{array}$$

Further details can be found in Zaehle and Friend (2010) and Vuichard et al. (2019).

A major improvement from IPSL-CM6A-LR to IPSL-Perm-LandN is the vertical discretisation of soil organic carbon and nitrogen on an 18-layer scheme (the same as for soil thermal dynamics), with depth-dependent decomposition rates depending on environmental conditions. This is particularly important in permafrost regions where the upper soil layers can thaw while deeper layers remain frozen, keeping organic matter thermally protected. Soil mineral nitrogen, however, is not vertically resolved and remains represented on a single soil layer in each grid box. It can exchange nitrogen with all the organic nitrogen layers through mineralisation or immobilisation.

Organic carbon and nitrogen can be exchanged between soil layers through bio- or cryoturbation. This process is described by a diffusion equation:

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\partial C_i}{\partial t}}{|}_{\mathrm{cryoturbation}}=D{\displaystyle \frac{\partial C_i^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\end{array}$$

where C_i is the carbon or nitrogen content of the pool i at a given depth and time, and D is the diffusive mixing rate. In the permafrost region (defined as ALT ≤ 3 m), D is set to 10⁻³ m² yr⁻¹ in the active layer and decreases linearly to zero between ALT and 3 × ALT. Thus, the permafrost region where cryoturbation occurs is dynamic. Elsewhere, D is set to 10⁻⁴ mm² yr⁻¹ in the top 2 m of soil to represent bioturbation.

The depth-dependent decomposition of soil organic matter depends on environmental conditions. In particular, it is modulated as a function of temperature (f(T) in Eq. 19):

$$\begin{array}{}\text{(25)}& f\left(T\right)=\left\{\begin{array}{ll}\exp \left(\log \left(Q_{\mathrm{10}}\right)\left(\frac{T-T_{\mathrm{ref}}}{\mathrm{10}}\right)\right)& \text{if}T>\mathrm{0}\mathit{°}\mathrm{C}\\ (T+\mathrm{1})\cdot \exp \left(-\log \left(Q_{\mathrm{10}}\right)\left(\frac{T_{\mathrm{ref}}}{\mathrm{10}}\right)\right)& \text{if}-\mathrm{1}\mathit{°}\mathrm{C}<T\le \mathrm{0}\mathit{°}\mathrm{C}\\ \mathrm{0}& \text{if}T\le -\mathrm{1}\mathit{°}\mathrm{C}\end{array}\right.\end{array}$$

where Q₁₀=2 and T_ref= 30 °C. Above 0 °C, decomposition follows a Q₁₀ function (Q₁₀= 2), then decreases linearly to zero between 0 and −1 °C. Below −1 °C no decomposition can take place.

Decomposition also increases monotonically with soil moisture (f(moisture) in Eq. 19):

$$\begin{array}{}\text{(26)}& \begin{array}{rl}f\left(\mathrm{moisture}\right)& =\max (\mathrm{0.25};-\mathrm{1.1}\cdot {\mathrm{moisture}}^{\mathrm{2}}+\mathrm{2.4}\\ & \cdot \mathrm{moisture}-\mathrm{0.29})\end{array}\end{array}$$

where moisture represents the humidity profile (unitless) and is between 0 and 1. Below 2 m (the depth to which hydrology is resolved), a constant soil moisture profile is used, taken from the lowest layer.

Overall, for each soil layer, the organic matter dynamics follows the equation below:

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{\displaystyle \frac{\partial C_i(z,t)}{\partial t}}& =I_i(z,t)-k_i\cdot f\left(T\right)(z,t)\cdot f\left(\mathrm{moisture}\right)(z,t)\\ & \cdot f\left(\mathrm{texture}\right)\cdot C_i(z,t)+D(z,t){\displaystyle \frac{\partial C_i^{\mathrm{2}}(z,t)}{\partial z^{\mathrm{2}}}}\end{array}\end{array}$$

where I_i are the carbon or nitrogen inputs to the pool i, the second term corresponds to decomposition and the third term to vertical mixing.

2.4.7 Initialisation of soil organic carbon and nitrogen

IPSL-Perm-LandN is unable to build up the observed large permafrost carbon stocks from scratch during spinup (even covering several thousands of years) due to the constant pre-industrial climate forcing of the spinup (i.e. no glacial/interglacial cycles), the long timescales required for carbon burial, missing processes (dust deposition, peat development) and the lack of deep permafrost deposits. Consistent representation of permafrost soil carbon is critical to avoid biases in its insulating effect or underestimation of future permafrost CO₂ emissions. Therefore, soil organic carbon and nitrogen pools are initialised with the contemporary observation-based product SoilGrids, which provides a global map of soil organic carbon and nitrogen with a detailed depth resolution (version 2.0, Poggio et al., 2021; Batjes et al., 2019). This allows the unfrozen soil layers to reach an equilibrium state driven by the carbon cycle and climate dynamics, while the organic matter in the frozen layers cannot be decomposed throughout the spinup. SoilGrids gathers observations from about 240 000 locations and uses more than 400 covariates. The original product has a horizontal resolution of 250 m and 6 vertical layers down to 2 m (0–5, 5–15, 15–30, 30–60, 60–100 and 100–200 cm). It has been conservatively regridded to the ORCHIDEE horizontal grid (2.5° × 1.25°) using the CDO “remapcon” command, and vertically interpolated to the 18-layer scheme. Below 2 m, initial organic carbon and nitrogen have been set to zero. Organic carbon and nitrogen stocks were divided into active, slow and passive fractions following the fractions given in Koven et al. (2015b) (2 % in active, 29 % in slow and 69 % in passive pools). As there is no global gridded map of soil mineral nitrogen, the mineral nitrogen pool is initialised to zero prior to the spinup.

3 Methods

3.1 Simulations and forcings

3.1.1 Spinup

Before running IPSL-Perm-LandN under varying forcings, it is necessary to bring the carbon and nitrogen pools into equilibrium, although such a target is questionable given that parts of the carbon cycle were not at equilibrium in 1850 (e.g. permafrost soils and peatlands were accumulating carbon, Schuur et al., 2022). This is done by performing a spinup in pre-industrial configuration. The spinup protocol starts with a spinup using ORCHIDEE offline (i.e. not coupled to the atmosphere and the ocean) under pre-industrial conditions for 2600 years. The model is forced by a 50-year cyclic climate from the spinup of IPSL-CM6A-LR (piControl simulation of CMIP6), which has identical atmosphere and ocean physics to that of IPSL-Perm-LandN. The PFT map (Lurton et al., 2020) and nitrogen deposition (National Center for Atmospheric Research-Chemistry-Climate Model Initiative) and fertilisation (Hurtt et al., 2020) remain at their 1850 values. Biological nitrogen fixation follows the approach of Cleveland et al. (1999) and is fixed in time (Vuichard et al., 2019). ORCHIDEE is then coupled to LMDZ (atmosphere) and NEMO (ocean) to form IPSL-Perm-LandN. The model is restarted from the offline ORCHIDEE spinup for land variables and from a spinup of IPSL-CM6A-LR for atmosphere and ocean variables. Importantly, the restart state of the ocean is from a 4000-year simulation, providing initial already equilibrated ocean physics and carbon pools. The spinup is run in concentration-driven configuration for 670 years. The land forcings remain the same and the atmospheric and oceanic forcings are fixed at their pre-industrial values. In particular, the atmospheric CO₂ concentration is set to 284 ppm. After the spinup, the coupled model is considered to be sufficiently close to equilibrium to avoid significant drifts in global climate variables and in the land and ocean net carbon fluxes in historical simulations (see Table A3).

3.1.2 Historical simulations

Three historical simulations (1850–2014) were performed with IPSL-Perm-LandN following the CMIP6 protocol in order to quantify the uncertainty in the simulated processes due to internal model variability. They differed only in their restart state, as the model was restarted from three distinct pre-industrial climate states (years 420, 450, and 480). These restart points were verified to be significantly different in terms of global temperature, thus providing three distinct restart states within the internal variability of IPSL-Perm-LandN. The forcings are provided by the CMIP6 input4MIP project (https://aims2.llnl.gov/search/input4MIPs/, last access: 15 May 2025), including greenhouse gas concentrations, which were taken as global averages from Meinshausen et al. (2017). Tropospheric and stratospheric ozone radiative forcings came from Checa-Garcia et al. (2018) and Hegglin et al. (2016). Tropospheric aerosols were not simulated interactively by IPSL-Perm-LandN and were prescribed from a historical LMDZOR-INCA simulation (i.e. a coupled surface-atmosphere simulation with tropospheric chemistry). In addition, stratospheric (volcanic) aerosols were prescribed from the version 3 of the dataset from Thomason et al. (2018) as a latitude-height time-varying climatology. Finally, the solar forcing is provided by Matthes et al. (2017).

Atmospheric nutrient deposition to the ocean (iron, phosphorus, and silicate) was provided by LMDZOR-INCA simulations. Wet and dry oceanic deposition of nitrogen (inorganic nitrate and ammonium) came from the National Center for Atmospheric Research-Chemistry-Climate Model Initiative nitrogen deposition rates. The river supply of biogeochemical elements to the ocean was sourced from Mayorga et al. (2010) for dissolved inorganic and organic nitrogen, dissolved inorganic and inorganic phosphorus, and silicate. Dissolved inorganic carbon and alkalinity were provided by the simulations using the Global Erosion Model of Ludwig et al. (1996). The river supply of iron was calculated from the river supply of inorganic carbon, assuming a constant Fe/dissolved inorganic carbon ratio.

Land cover (i.e. the PFT map), wood harvest and nitrogen fertilisation are provided by the land use harmonisation database Hurtt et al. (LUH2, 2020). Nitrogen deposition is provided by th National Center for Atmospheric Research-Chemistry-Climate Model Initiative and BNF follows the approach of Cleveland et al. (1999).

A complete description of the implementation of the forcings can be found in Lurton et al. (2020).

3.2 Evaluation data

Surface air temperature data is taken from ERA5 reanalysis (Copernicus Climate Change Service, 2019) for absolute values and NOAAGlobalTemp (Huang et al., 2023) and HadCRUT (Morice et al., 2021) for temperature anomalies compared to 1850–1900. Total precipitation data come from ERA5 and MSWEP (Beck et al., 2019) and snowfall data from ERA5 only. Snow cover data come from the ESA-CCI CryoClim product (Solberg et al., 2023). Sea surface temperature and salinity come from the World Ocean Atlas (Locarnini et al., 2024; Reagan et al., 2024). Sea ice concentration is taken from the National Snow and Ice Data Center (DiGirolamo et al., 2022). The extent of the permafrost region is taken from ESA-CCI (Westermann et al., 2024a) and active layer thickness data come from ESA-CCI and the CALM network (Westermann et al., 2024b; Brown et al., 2000). GPP comes from the FLUXCOM network (Jung et al., 2020), RH from Konings et al. (2019), Warner et al. (2019) and Hashimoto et al. (2015) (the latter two based on Bond-Lamberty and Thomson, 2010), and NBP from the 2023 Global Carbon Budget (GCB2023, Friedlingstein et al., 2023) and the CAMS inversion product (Chevallier et al., 2023). Ocean net air-sea carbon flux come from GCB2023. Gridded data of vegetation biomass is taken from the ESA-CCI product (Santoro and Cartus, 2021) and soil carbon comes from HWSD (Wieder et al., 2014), SoilGrids (Poggio et al., 2021; Batjes et al., 2019) and NCSCD (Hugelius et al., 2013). Anthropogenic fossil emissions are from GCB2023 (Friedlingstein et al., 2023).

IPSL-Perm-LandN is compared to ESMs from the Coupled Climate-Carbon Cycle Model Intercomparison Project (C4MIP, Jones et al., 2016), which are part of the broader CMIP6 ensemble (C4MIP models are listed in Arora et al., 2020). These models represent interactive land and ocean carbon cycle and can therefore represent carbon cycle feedbacks. Data for C4MIP models has been retrieved from the IPSL ESGF node (https://esgf-node.ipsl.upmc.fr/projects/esgf-ipsl/, last access: 15 May 2025) at the time of the study. For each model, the first 10 members are used, except for UKESM1-0-LL and NorESM2-LM where only 4 and 3 members were available, respectively. For IPSL-CM6A-LR, the 33 members are used.

3.3 Evaluation metrics

3.3.1 Permafrost region

A necessary but tricky step in the study of permafrost modeling is to clearly define permafrost in the model. A first clarification is needed to avoid the common confusion between the permafrost region and the permafrost area (Obu, 2021). The permafrost region is defined as the total area covered by permafrost zones (continuous, discontinuous, sporadic and isolated patches). However, each permafrost zone is not completely underlain by permafrost and the actual area underlain by permafrost is smaller than the permafrost region. This area actually underlain by permafrost is called the permafrost area, and takes into account, for example, that there is more permafrost in the continuous than in the sporadic zone. Many observation products provide both the permafrost region and the permafrost area (Obu, 2021; Obu et al., 2019; Gruber, 2012). In Earth System Models, however, each pixel of the grid either contains permafrost or does not. A finer description of permafrost would require the representation of sub-grid land surface heterogeneity and the estimation of a permafrost fraction for each pixel, which is not the case in current ESMs despite promising developments (Shirley et al., 2022; Cai et al., 2020; Beer, 2016; Fowler et al., 2024; Torres-Rojas et al., 2022). Thus ESMs can only represent the permafrost region as the total area where grid boxes contain permafrost. However this modeled permafrost region is slightly different from the one estimated from observations. As the ESMs represent the dominant environmental conditions over each grid box, areas with small amounts of permafrost are likely to be missing permafrost. On the contrary, in areas with observed permafrost fractions greater than 50 %, the majority of the area is underlain by permafrost and the models should consider them as pixels containing permafrost. Thus, continuous and discontinuous permafrost zones (>50 % of permafrost) should be similar between models and observations while disagreement is expected for sporadic permafrost and isolated patches (<50 % of permafrost).

Apart from this, a second source of uncertainty comes from the way in which is decided whether a model grid box contains permafrost or not. Comparing 10 different definitions of permafrost in ESMs, Steinert et al. (2024) found large differences within each model of the CMIP6 ensemble and showed that the spread due to permafrost definition could even be larger than the inter-model spread. Among the classical permafrost definitions, those based on ground-air temperature coupling show a better agreement between models but miss the complexity introduced by ground thermodynamics by implicitly assuming the same ground thermodynamics for all models. More relevant definitions are based on ground thermal properties and are closer to the original definition of permafrost. A direct application of this definition in models would be to define the zero annual amplitude depth (D_zaa) as the minimum soil depth at which the temperature variation within a year is less than 0.1 °C. If the temperature at D_zaa is less than or equal to 0 °C for at least two consecutive years, there is assumed to be permafrost in the grid box (Burke et al., 2020). However the D_zaa can be deep, especially in models with a deep soil column such as IPSL-Perm-LandN. With this definition, if deep permafrost is modeled, the grid box is marked as containing permafrost. This can be problematic if the lower soil layers are poorly represented. For instance, the lower ground boundary condition in IPSL-Perm-LandN does not represent the heat coming from the Earth's mantle, resulting in an incorrect geothermal gradient. This can cause deep ground to remain unrealistically frozen and to overestimate the area of permafrost using this definition. This is why in this study, we chose another commonly used permafrost definition, based on the active layer thickness (ALT) (McGuire et al., 2018; Koven et al., 2013). If the ALT is less than 3 m, i.e. if the annual maximum thaw depth is less than 3 m, the grid box is said to contain permafrost. This definition includes surface permafrost but excludes deep permafrost (i.e. below 3 m), which is fine for two reasons:

• IPSL-Perm-LandN poorly represents deep soil temperature profile and focusing on surface permafrost avoids overestimating the permafrost region.

• The vast majority of soil organic carbon is in the top 3 m of soil in IPSL-Perm-LandN and soil carbon decomposition following permafrost thaw would occur within the top 3 m of soil.

Thus we chose to define the permafrost region (ℛ_permafrost) as the total area where ALT <3 m, i.e.:

$${\mathcal{R}}_{\text{permafrost}}=\sum _{\text{ilon=1}}^{\mathrm{144}}\sum _{\text{ilat=1}}^{\mathrm{143}}\mathit{\delta }\text{(ilon,ilat)}\cdot \mathcal{A}\text{(ilon,ilat)}\cdot f_{\text{land}}\text{(ilon,ilat)}$$

with f_land(ilon,ilat) the fraction of land in the grid box, 𝒜(ilon,ilat) the grid box area and

$$\mathit{\delta }\text{(ilon,ilat)}=\left\{\begin{array}{ll}\text{1}& \text{if}\mathrm{ALT}<\mathrm{3}\mathrm{m}\\ \text{0}& \text{otherwise}\end{array}\right..$$

The permafrost region is calculated for IPSL-Perm-LandN, IPSL-CM6A-LR, C4MIP models and the ESA-CCI observation product (Westermann et al., 2024a). As some C4MIP models have a poorly resolved soil thermal profile, an exponential vertical interpolation at 3 m depth is performed instead of taking the temperature of the nearest soil layer. If the interpolated 3 m-temperature is less than or equal to 0 °C, the ALT is less than 3 m and the grid box contains permafrost. For IPSL-Perm-LandN, the yearly maximum ALT is directly available and is used to calculate the size of the permafrost region (altmax <3 m). The ESA-CCI observation product provides the permafrost fraction (f_perm) for each pixel, which allows the calculation of the permafrost region (area where f_perm>0), the permafrost area (area weighted by f_perm) and the region of continuous and discontinuous permafrost (area where f_perm>0.5).

3.3.2 Active layer thickness

The spatially-averaged time evolution of the active layer thickness is computed using a mask of the permafrost region. This mask is defined as the simulated 2005–2014 permafrost region, using the definition ALT < 3 m.

3.3.3 Compatible CO₂ emissions

Instead of prescribing anthropogenic CO₂ emissions to IPSL-Perm-LandN, the historical simulations are run with an imposed atmospheric CO₂ concentration. This prevents the simulated land and ocean carbon fluxes from feeding back onto climate, removing a source of uncertainty for the study of atmospheric processes, despite the use of a spatially homogeneous CO₂ concentration with no vertical gradient. However, these fluxes can be used in addition to atmospheric CO₂ changes to calculate the fossil fuel emissions that are compatible with the prescribed CO₂ concentration scenarios. The rate of compatible fossil fuel emissions is equal to the sum of the rate of atmospheric CO₂ change, the net atmosphere-land and atmosphere-ocean fluxes, i.e. :

$$\begin{array}{}\text{(28)}& E_{\mathrm{FF}}=G_{\mathrm{ATM}}+F_{\text{A–O}}+F_{\text{A–L}}\end{array}$$

with E_FF the rate of anthropogenic fossil fuel emissions (PgC yr⁻¹), G_ATM the rate of change of atmospheric CO₂ concentration (PgC yr⁻¹), F_A–O the net atmosphere-ocean flux (PgC yr⁻¹, positive for ocean uptake) and F_A–L the net atmosphere-land flux (PgC yr⁻¹, positive for land uptake). Land-use change emissions are included in the NBP, and therefore in F_A–L.

4 Results and discussion

4.1 Atmosphere physics

Over the period 1940–2014, the mean annual land surface air temperature (SAT) of IPSL-Perm-LandN is about 1.5 °C colder than the ERA5 reanalysis (Fig. 1a). During the last decade of the simulation (2005–2014), the mean land SAT of IPSL-Perm-LandN is 13.46 ± 0.14°C while ERA5 has a warmer land SAT of 14.84 °C. IPSL-Perm-LandN is consistently very close to IPSL-CM6A-LR as both share the same radiative scheme, and is at the lower bound of the C4MIP range, although the models generally tend to correctly simulate temperature changes (i.e. ΔSAT) rather than absolute temperatures. The cold land SAT bias in IPSL-Perm-LandN is mainly due to underestimated tropical and mid-latitude temperatures across all seasons while the Arctic land SAT is closer to ERA5 estimates, due to canceling cold and warm biases in spring and autumn, respectively (Figs. 1b and A5). These biases could impact permafrost freeze and thaw but are unevenly distributed across the region (Fig. A6). Although the absolute land temperature is too cold, the land SAT anomaly relative to 1850–1900 is close to observations. Over land (emerged land excluding Greenland and Antarctica), IPSL-Perm-LandN has warmed by +1.60 ± 0.14 °C (mean 2005–2014 warming compared to 1850–1900) while the observations show a warming of +1.40°C for NOAAGlobalTemp and +1.16 °C for HadCRUT (Fig. 1c). In contrast to the absolute temperature, the land SAT change compared to 1850–1900 is at the upper limit of the range of the C4MIP models. This relatively high warming mainly comes from the tropics and the Arctic where land SAT change (ref. 1850–1900) is overestimated compared to both NOAAGlobalTemp and HadCRUT (Figs. 1d and A7a). In particular, the Arctic amplification is overestimated in IPSL-Perm-LandN with a high latitude warming twice as large as in the observations. This Arctic warming bias was already present in IPSL-CM6A-LR and is amplified in IPSL-Perm-LandN. In addition, when including the oceans to compute the global surface air temperature (GSAT) anomaly, IPSL-Perm-LandN deviates from the observations and starts to warm faster from 1990 onwards, driven by a strong oceanic warming in the Arctic ocean (Fig. A8). The mean global warming for 2005–2014 relative to 1850–1900 is +1.27 ± 0.12 °C for IPSL-Perm-LandN and +0.84 °C (NOAAGlobalTemp) and +0.80 °C (HadCRUT) for observation-based datasets. This departure from observations in the recent period was already present in IPSL-CM6A-LR and depends on the reference period used to compute the anomaly (Boucher et al., 2020).

Figure 1 Historical surface temperature over land. (a) Mean land surface air temperature (SAT) over the historical period for IPSL-Perm-LandN, IPSL-CM6A-LR and ERA5 reanalysis. Colored dots represent the mean land SAT (2005–2014) for IPSL-Perm-LandN, IPSL-CM6A-LR, ERA5 and C4MIP models. Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (b) Mean land SAT (2005–2014) over the Arctic (>60 ° N), mid-latitudes (30–60° S and 30–60° N) and the tropics (30° S–30° N) for IPSL-Perm-LandN and C4MIP models. (c) Anomaly of mean land SAT relative to 1850–1900 for IPSL-Perm-LandN, IPSL-CM6A-LR, NOAAGlobalTemp and HadCRUT reanalyses. Colored dots represent the mean land SAT anomaly (2005–2014) for IPSL-Perm-LandN, IPSL-CM6A-LR, NOAAGlobalTemp, HadCRUT and C4MIP models. Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (d) Mean land land SAT anomaly over the Arctic (>60° N), mid-latitudes (30–60° S and 30–60° N) and the tropics (30° S–30° N) for IPSL-Perm-LandN, IPSL-CM6A-LR and C4MIP models, compared to 1850–1900.

The mean total precipitation (liquid + solid) in IPSL-Perm-LandN for the period 2005–2014 is shown in Fig. 2a. The latitudinal distribution of precipitation is very close to the observations in the Arctic and mid-latitudes (Fig. 2a and c). In the tropics, although the model correctly represents the ITCZ, it has a pronounced peak at 5° S, which is much lower in the observations. Such a double ITCZ is a known bias in many CMIP6 models and could be due to the representation of deep convection as well as model resolution (Ma et al., 2023). The mean total snowfall is represented in Fig. 2b. Its latitudinal distribution follows that of ERA5 and, in particular, the Arctic snowfall is well represented in the recent period (Fig. 2b and d). However, the good agreement between IPSL-Perm-LandN and ERA5 masks a slight overestimation of Arctic snowfall over land and a slight underestimation over the ocean. In addition, the mean seasonality of both total precipitation and snowfall is well captured by the model in the Arctic (Fig. A9). This slight overestimation of Arctic snowfall does not lead to significant snow cover biases (Fig. A10). However, snow cover is underestimated by 10 % to 20 % in the permafrost region in April–May and October-November, which could lead to reduced ground insulation and faster thawing and refreezing of permafrost in spring and autumn. In the mid-latitudes, the seasonal cycle of snowfall is well represented while total precipitation is overestimated by up to 0.16 mm d⁻¹ (∼ 6 %), except in late summer (Fig. A9). Although total precipitation has a double ITCZ in the tropics, the amplitude and phase of its seasonal cycle are in agreement with observations. In general, both total precipitation and snowfall are close to those of IPSL-CM6A-LR.

Figure 2 Historical global precipitation and snowfall. (a) Left: map of mean total precipitation (liq + sol) over 2005–2014 for IPSL-Perm-LandN. Right: zonal mean of total precipitation over 2005–2014 for IPSL-Perm-LandN, ERA5 reanalysis and MSWEP observation product. (b) Left: map of mean snowfall (2005–2014) for IPSL-Perm-LandN. Right: zonal mean of snowfall (2005–2014) for IPSL-Perm-LandN and ERA5. (c) Difference in mean total precipitation (liq + sol) between IPSL-Perm-LandN and ERA5 over 2005–2014. (d) Difference in mean snowfall between IPSL-Perm-LandN and ERA5 over 2005–2014.

4.2 Ocean physics

The sea surface temperature (SST) mean pattern computed over the historical period in IPSL-Perm-LandN is quite similar to that of IPSL-CM6A-LR, as the same version of the ocean model NEMOv3.6 was used. The main bias in IPSL-Perm-LandN is a negative SST anomaly in the North Atlantic ocean compared to observations from the World Ocean Atlas over the period 2005–2014, which is associated with the position of the North Atlantic drift and due to a weaker AMOC than IPSL-CM6A-LR (Fig. 3a and d). This bias was already present in IPSL-CM6A-LR but was less pronounced (Boucher et al., 2020) (Fig. A11a). The maximum temperature negative anomaly around 45° N (in the box 60–15° W, 40–55° N) for the period 2005–2014 is −7.2 °C in IPSL-Perm-LandN while it was −5.5°C for IPSL-CM6A-LR. Such a cold bias is a common feature of CMIP6 models and is stronger in winter (Zhang et al., 2023). Other classical SST biases of CMIP6 models are present in IPSL-Perm-LandN : warm biases in eastern ocean borders (although not very strong along South America), cold mid-latitudes and a warm bias near Antarctica (Zhang et al., 2023; Boucher et al., 2020). Sea surface salinity (SSS) also shows similar patterns as IPSL-CM6A-LR (Figs. 3b and A11b). A negative salinity anomaly is observed in the North Atlantic – in the same region as the cold SST bias – but has been reduced in IPSL-Perm-LandN, although exact reasons are yet unclear. As in IPSL-CM6A-LR, the eastern equatorial Pacific ocean is too salty compared to the World Ocean Atlas. This could be due to an underestimation of precipitation in the area, which would reduce the dilution effect (Fig. 2c). Similarly, positive and negative salinity biases are consistent with precipitation biases, suggesting that SSS biases could be driven by precipitation.

Figure 3 Historical ocean physics for IPSL-Perm-LandN. Difference in annual mean sea surface (a) temperature and (b) salinity between IPSL-Perm-LandN the World Ocean Atlas (2005–2014). (c) Mean annual maximum mixed layer depth (2005–2014). (d) Atlantic meridional overturning stream function, on average over 2005–2014.

The Atlantic Meridional Overturning Circulation (AMOC) cell has a very similar latitudinal extent and a maximum around 40° N but its strength is lower than for IPSL-CM6A-LR (Figs. 3d and A11d). The sign of the AMOC stream function changes around 2200 m depth while it changes around 2500 m for IPSL-CM6A-LR. In the short observational dataset available, this change is diagnosed to occur around 4500 m. This shallow AMOC cell is a known bias of the IPSL model (Boucher et al., 2020). The maximum mixed layer depth (MLD) is maximum in the Labrador and Nordic seas, indicating areas of dense water production (Figs. 3c and A11c). The location of the MLD maxima is consistent with observations in spite of a large variability among members (Boucher et al., 2020). The MLD of IPSL-Perm-LandN is shallower than that of IPSL-CM6A-LR, which is consistent with a weaker AMOC and suggests a reduced production of dense water in the northern North Atlantic.

The March sea ice extent – generally the annual sea ice maximum extent – is overestimated by the IPSL-Perm-LandN when compared to NSIDC observations (Fig. A12). Over the historical period, the March sea ice extent decreases from 20.3 Mkm² (1850–1900) to 16.7 Mkm² (2005–2014), while observations show a slower decrease over the last decades and yet a weaker total sea ice extent of 14.8 Mkm² (2005–2014). In the last years of the historical simulation, the model comes closer to the satellite observations. On the contrary, the March sea ice extent was very close to the observations in IPSL-CM6A-LR. Almost all of the difference is explained by the presence of sea ice at the Labrador sea-Atlantic junction in winter with fractions close to 1 in IPSL-Perm-LandN, while this area is almost ice-free in IPSL-CM6A-LR (Fig. A12a and b). This is consistent with the strong cold SST bias, the strong reduction of MLD in the Labrador sea and the weakening of the AMOC previously observed in IPSL-Perm-LandN. The annual minimum sea ice area (in September) is also slightly overestimated by IPSL-Perm-LandN, but less than for winter sea ice. The decreasing trend in the simulations is consistent with observed trends.

4.3 Permafrost physics

In IPSL-Perm-LandN, the permafrost region covers 16.5 Mkm² at the end of the historical simulation (2005–2014) (Fig. 4a). This is higher than the ESA-CCI mean permafrost area (regridded to the resolution of IPSL-Perm-LandN) (14.0 Mkm²), but just below the upper limit of uncertainty, and lower than the ESA-CCI permafrost region (mean 19.3 Mkm²). This was expected as the ESA-CCI permafrost area represents the area underlain by permafrost, that the model cannot represent and which is smaller than the permafrost region. In addition, as the ESA-CCI permafrost region is the region covered by all permafrost zones, it results in a larger estimate than the models that cannot capture sporadic permafrost and isolated patches. However, the simulated permafrost region is slightly higher than the ESA-CCI continuous and discontinuous permafrost region (permafrost fraction >50 %, mean 14.17 Mkm²) that the model is expected to simulate, mainly due to overestimated permafrost extent over the Tibetan Plateau. The modeled permafrost region is also within the range of C4MIP models estimates, although they have not been regridded and the permafrost representation of each model is superimposed to the effect of its spatial resolution. Higher resolution models should, in principle, be closer to observations as they capture finer permafrost patterns. Notably, there is a clear improvement in the representation of permafrost compared to IPSL-CM6A-LR which had an extremely small permafrost region. This is mainly due to the inclusion of the latent heat of soil water phase change in IPSL-Perm-LandN. Its absence in IPSL-CM6A-LR resulted in overestimated ALT and underestimated permafrost region (Steinert et al., 2024). In the recent period, the permafrost region is very close for all three simulation members, with only small differences at the southern permafrost edges (Fig. 4b). Overall, there is a very good agreement between IPSL-Perm-LandN and the ESA-CCI product (permafrost fraction >50 %). In Eurasia, the permafrost region compares well with the 50 % permafrost contour from ESA-CCI observations, with a slight overestimation over the southern boundary, which could be due to a legacy effect of the spring cold bias in this region (Fig. A6). As expected, the model also predicts too much permafrost over the Tibetan Plateau, which has a known cold bias in surface air temperature (Boucher et al., 2020). In North America, simulated permafrost in IPSL-Perm-LandN is present in the north, but is absent at the southern edge, in Canada, which is not clearly related to a warm temperature bias (Fig. A6) but is a known bias in many CMIP6 models (Burke et al., 2020). Overall, the permafrost region has decreased by 2.4 Mkm² (−15.0 %) over the historical period compared to 1850–1900.

Figure 4 Historical permafrost region. (a) Permafrost region in the northern hemisphere over the historical period for IPSL-Perm-LandN, IPSL-CM6A-LR and C4MIP models, and permafrost area and permafrost region for ESA-CCI observation product. Colored dots represent the mean permafrost region (2005–2014) for IPSL-Perm-LandN, IPSL-CM6A-LR and C4MIP models and the mean permafrost area, permafrost region and region of >50 % permafrost for ESA-CCI. Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (b) Map of the permafrost region in IPSL-Perm-LandN (2005–2014). Dark blue: all three members diagnose permafrost. Light Blue: two members diagnose permafrost. Orange: only one member diagnoses permafrost. Red contour: 50 % permafrost fraction (continuous and discontinuous) from ESA-CCI.

The simulated mean ALT of IPSL-Perm-LandN is in good agreement with CALM observations in eastern and northern Canada, and northern and eastern Siberia (Fig. 5). However, it is too deep in western Siberia, western Alaska and along the MacKenzie river in western Canada. It also compares well to the ESA-CCI product over most of the permafrost region. At the southern edge of Canadian permafrost, there is no permafrost in IPSL-Perm-LandN and the ALT is expectedly too deep. Within the modeled permafrost region, the simulated ALT is also too deep in Western Alaska and Western Siberia, the latter being partly due to the underestimation of ALT in this area by the ESA-CCI product (Fig. A13).

Figure 5 Active layer thickness (2005–2014). (a) Background: map of ALT for IPSL-Perm-LandN (2005–2014). Colored circles: CALM observations. (d) Background: difference of ALT between IPSL-Perm-LandN and ESA-CCI (2005–2014).

4.4 Global land carbon cycle dynamics

4.4.1 Growth Primary Production (GPP)

On a global scale, gross primary production (GPP) increases slowly until the 1960's and much faster thereafter for both IPSL-Perm-LandN and IPSL-CM6A-LR (Fig. 6a). As in other ESMs, this bent curve is mainly driven by the fertilisation effect caused by an increase in anthropogenic CO₂ emissions (Piao et al., 2009; Schimel et al., 2015) as well as increased nitrogen atmospheric deposition and fertilisation (Huntzinger et al., 2017; O'Sullivan et al., 2019). The change in the slope around the 1960s is more pronounced for IPSL-Perm-LandN than for IPSL-CM6A-LR, primarily driven by the explicit representation of the nitrogen cycle in IPSL-Perm-LandN and its effect in the tropics and mid-latitudes. In IPSL-Perm-LandN, the global GPP reaches 132 PgC yr⁻¹ in the last decade of the simulation, higher than estimates from Jung et al. (2020) but within the range of C4MIP ESMs, although there is a large variability across models. GPP is overestimated in the Arctic and mid-latitudes compared to data-driven products, and within the observational range in the tropics (Figs. 6b and A14a). This is likely due to IPSL-Perm-LandN simulating larger organic nitrogen stocks in the mid-latitudes and the Arctic than in the tropics, leading to higher mineralisation under warming, and therefore to greater sensitivity of nitrogen limitation to warming (Fig. A15). Compared to IPSL-CM6A-LR, GPP has largely increased in the Arctic (+3.3 PgC yr⁻¹) and mid-latitudes (+15.4 PgC yr⁻¹), and decreased in the tropics (−10.1 PgC yr⁻¹), resulting in an overall global increase of 8.6 PgC yr⁻¹. These differences are explained by the introduction of an explicit nitrogen cycle, which replaces an empirical GPP downregulation in IPSL-CM6A-LR (limitation of Vcmax under increasing atmospheric CO₂ to mimic nutrient limitation without explicitly representing it), and to vertically-resolved soil biogeochemistry in IPSL-Perm-LandN. In addition, IPSL-CM6A-LR has been largely tuned using different data sources (FLUXNET, atmospheric CO₂, NDVI, Peylin et al., 2016), while the new model including the nitrogen cycle has not been extensively calibrated (see Appendix A4). The seasonal cycle was improved in the tropics compared to IPSL-CM6A-LR, with a seasonality closer to data-driven estimates (Fig. A14a). In the northern mid-latitudes, the shape of the seasonal cycle is consistent with the observations but its amplitude is too large. IPSL-Perm-LandN captures the onset of vegetation growth well, but overestimates GPP during the summer peak and vegetation senescence in autumn. In contrast, IPSL-CM6A-LR was very close to data-driven products throughout the year. In the Arctic, IPSL-Perm-LandN overestimates the amplitude of the seasonal cycle, but also shows a delayed decrease in GPP in late summer and autumn. This was already the case for IPSL-CM6A-LR and is partly due to a warm autumn bias in the Arctic which allows vegetation to survive later in the season (Figs. A6 and A16).

Figure 6 GPP over the historical period. (a) Global GPP over the historical period for IPSL-Perm-LandN, IPSL-CM6A-LR, C4MIP models, Jung-RS and Jung-RSMETEO observation products (Jung et al., 2020). Colored dots represent the mean GPP (2005–2014) for IPSL-Perm-LandN, IPSL-CM6A-LR, C4MIP models and observations products. Plain (resp. empty) circles represent models with (resp. without) an explicit nitrogen cycle. Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (b) Total GPP (2005–2014) over the Arctic (>60° N), mid-latitudes (30–60° S and 30–60° N) and the tropics (30° S–30° N) for IPSL-Perm-LandN, IPSL-CM6A-LR and C4MIP models.

4.4.2 Soil heterotrophic respiration (RH)

For IPSL-Perm-LandN, the soil heterotrophic respiration (RH) follows the same bent shape as GPP over the historical period (Fig. 7). This was expected as enhanced GPP leads to increased litter and soil carbon, resulting in higher RH. In the last decade of the simulation, RH reaches 47.4 PgC yr⁻¹ for IPSL-Perm-LandN, close to IPSL-CM6A-LR (45.5 PgC yr⁻¹) and data-driven products (43.4 PgC yr⁻¹ for Konings et al., 2019, 48.8 PgC yr⁻¹ for Warner et al., 2019 and 51.9 PgC yr⁻¹ for Hashimoto et al., 2015). Similar to IPSL-CM6A-LR, IPSL-Perm-LandN is one of the ESMs with the globally simulated RH value that is closest to these data-driven products over the recent period (Guenet et al., 2024). However, even if the global RH is close to IPSL-CM6A-LR, the use of a discretised soil carbon profile, the inclusion of permafrost and of an explicit nitrogen cycle in IPSL-Perm-LandN leads to very different regional RH patterns. As with GPP, RH has increased in the Arctic and mid-latitudes, and decreased in the tropics compared to IPSL-CM6A-LR. The modeled RH for IPSL-Perm-LandN is in good agreement with the Warner et al. (2019) and Hashimoto et al. (2015) products globally, but is slightly overestimated over forests (Fig. A17). In tropical and mid-latitude grassland ecosystems, RH tends to be underestimated. As expected, given their correlation, GPP and RH show the same regional biases when confronted with independent observational products.

Figure 7 Soil heterotrophic respiration over the historical period. (a) Global soil heterotrophic respiration (RH) over the historical period for IPSL-Perm-LandN, IPSL-CM6A-LR, C4MIP models, Konings et al. (2019), Warner et al. (2019) and Hashimoto et al. (2015) observational products. Colored dots represent the mean RH (2005–2014) for IPSL-Perm-LandN, IPSL-CM6A-LR, C4MIP models and observations. Plain (resp. empty) circles represent models with (resp. without) an explicit nitrogen cycle. Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (b) Total RH (2005–2014) over the Arctic (>60° N), mid-latitudes (30–60° S and 30–60° N) and the tropics (30° S–30° N) for IPSL-Perm-LandN, IPSL-CM6A-LR and C4MIP models.

4.4.3 Net land-atmosphere carbon flux (NBP)

For both models, the net land-atmosphere carbon flux (NBP, positive for land uptake), including land-use change emissions (E_LUC), is negative until the 1970's, mainly because of the negative contribution of land-use change (Tharammal et al., 2019). Thereafter, NBP increases, driven by CO₂ fertilisation and nitrogen fertilisation to reach 1.83 ± 0.34 PgC yr⁻¹ for IPSL-Perm-LandN in the last decade (2005–2014), making the land a net carbon sink over the last 50 years (Fig. 8a). This value is very close to estimates from the 2023 Global Carbon Budget (1.86 ± 1.13 PgC yr⁻¹), which uses offline Dynamic Global Vegetation Models (DGVMs) for the land carbon sink and bookkeeping models for land-use change emissions. The NBP is slightly larger for IPSL-Perm-LandN than IPSL-CM6A-LR (1.52 ± 0.78 PgC yr⁻¹) due to increased net carbon uptake in the Arctic and mid-latitudes, while the tropical NBP remains similar (Fig. 8b). The inverse modeling approach used by the Copernicus Atmosphere Monitoring Service (CAMS) shows a higher global NBP (2.89 PgC yr⁻¹) and a different latitudinal distribution. This is due to the fact that atmospheric inversions account for lateral carbon fluxes (between the land and the ocean), whereas land surface models (and hence ESMs) typically do not model this flux and have a near-zero land-atmosphere carbon flux in the pre-industrial period. In contrast, the global pre-industrial river flux is estimated to be around 0.65 PgC yr⁻¹ (Regnier et al., 2022). Subtracting the contribution of lateral fluxes from the inversions generally helps to reconcile both approaches, leading to more comparable NBP values (Ciais et al., 2021). However, there is still significant uncertainty in these estimates and the 2023 CAMS estimate has a relatively large land sink (Friedlingstein et al., 2023). In the tropics, the CAMS product diagnoses a net carbon source while all C4MIP ESMs rather show a positive to near-neutral NBP. At mid- and high-latitudes, the NBP is positive and much larger in the inversion than in C4MIP models, indicating a large net carbon sink that more than compensates for the tropical net carbon source. Such discrepancies between models and inversions are a known knowledge gap and an area of active research (Friedlingstein et al., 2023; Bastos et al., 2020). Recently, the work of O'Sullivan et al. (2024) has shown the key role of forest disturbances at mid-to-high latitude to reconcile the estimates of the northern carbon sink between atmospheric inversions and DGVMs. The seasonal cycle of NBP for IPSL-Perm-LandN is consistent with that of CAMS despite differences in amplitude (Fig. A18a). In general, IPSL-Perm-LandN has a smaller amplitude than CAMS in the tropics and a larger amplitude in the extra-tropics. This difference is greater during periods of negative NBP, especially during autumn and winter of the northern hemisphere. Over the last decade, the mean NBP is positive over most of the globe, with the notable exception of regions of high deforestation (eastern and southern Brazil, equatorial African forest, Indonesia) (Fig. A18b). Large sinks are simulated over Europe, Amazonian forest, western African forest, eastern China and the boreal forests of Canada, Alaska and Siberia. By removing the contribution of land-use change emissions in the NBP, we can estimate the land carbon sink (S_LAND in GCB2023), which is positive almost everywhere with deforested areas close to neutrality (Fig. A18c). However, we can only approximate S_LAND as it is calculated using fixed pre-industrial vegetation in GCB2023, whereas the vegetation evolves over time in our simulations. Therefore, the large spread in E_LUC hinders a more precise assessment of the land carbon sink in our simulations (Bastos et al., 2021; Friedlingstein et al., 2023).

Figure 8 Net land-atmosphere carbon flux over the historical period. (a) Global net land-atmosphere carbon flux (NBP) over the historical period for IPSL-Perm-LandN, IPSL-CM6A-LR, C4MIP models, CAMS inversion product and the Global Carbon Budget 2023. Positive (resp. negative) values correspond to a land carbon sink (resp. a source). Colored dots represent the mean NBP (2005–2014) for IPSL-Perm-LandN, IPSL-CM6A-LR, C4MIP models, CAMS and GCB2023. Plain (resp. empty) circles represent models with (resp. without) an explicit nitrogen cycle. Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (b) Total NBP (2005–2014) over the Arctic (>60° N), mid-latitudes (30–60° S and 30–60° N) and the tropics (30° S–30° N) for IPSL-Perm-LandN, IPSL-CM6A-LR, C4MIP models and CAMS.

4.4.4 Land carbon stocks

The global vegetation biomass (above- and below-ground, averaged over 2005–2014) amounts to 479 PgC for IPSL-Perm-LandN, which is lower than the ESA-CCI observation-based product (607 PgC, estimated in 2010), mainly due to the lower tropical biomass, and close to the mean of C4MIP models (Fig. 9). Compared to IPSL-CM6A-LR, the tropical biomass remained almost unchanged while it doubled at mid- and high-latitudes. Vegetation in the permafrost region remains limited (23 PgC) and smaller than the ESA-CCI estimate (37 PgC). Comparison with other models is provided for information, but it should be noted that the permafrost mask used here is that of IPSL-Perm-LandN while the permafrost region may differ between models. The total amount of litter carbon has remained similar since CMIP6 (108 PgC for IPSL-Perm-LandN and 107 PgC for IPSL-CM6A-LR) but its distribution has changed, with less carbon in the tropics and more in mid- and high-latitudes. However, there are no global scale observations and the spread across models is large, making it difficult to assess the performance of IPSL-Perm-LandN. Finally, IPSL-Perm-LandN simulates a total amount of SOC of 1985 PgC in 0–1 m (3001 PgC in 0–3 m), distributed between the tropics (521 PgC for 0–1 m, 639 PgC for 0–3 m), mid-latitudes (934 PgC for 0–1 m, 1376 PgC for 0–3 m) and the Arctic (530 PgC for 0–1 m, 985 PgC for 0–3 m). Compared to IPSL-CM6A-LR (total SOC of 550 PgC), it has largely increased in all latitudes, with the highest increases in the mid-latitudes and the Arctic. These changes are due to the discretisation of SOC along a vertical profile, the initialisation of the soil organic carbon and nitrogen pools by observation-based products and the representation of an explicit nitrogen cycle. Observed total SOC is 1204 PgC for HWSD and 2498 PgC for SoilGrids in 0–1 m (3384 PgC in 0–3 m). The large spread across these products and the resulting uncertainty in soil organic carbon content hampers a constrained assessment of ESMs on a global scale. IPSL-Perm-LandN is naturally closer to SoilGrids, which was chosen to initialise the SOC and SON pools due to the large number of observations, the robustness of the machine learning algorithm and the availability of gridded SOC and SON on 6 soil layers. Furthermore, the choice of a product with a high amount of SOC seems justified as global SOC gridded datasets tend to underestimate SOC content when compared to field data (e.g. Tifafi et al., 2018). Compared to CMIP6 models contributing to C4MIP (Arora et al., 2020), IPSL-Perm-LandN is the model with the highest amount of SOC, mainly due to large pools in the mid-latitudes and the Arctic. Permafrost SOC amounts to 511 PgC in the first meter of soil (1006 PgC in 0–3 m) and has largely increased compared to IPSL-CM6A-LR (46 PgC), which is a significant improvement of IPSL-Perm-LandN. It is again similar to the SoilGrids product (760 PgC in 0–1 m, 1028 PgC in 0–3 m) and larger than both NCSCD (282 PgC in 0–1 m, 668 PgC in 0–3 m) and HWSD (186 PgC). It is also very close to specific estimates of permafrost SOC stocks from Mishra et al. (2021) (${\mathrm{1014}}_{-\mathrm{170}}^{+\mathrm{186}}$ PgC) and Hugelius et al. (2014) (1035 ± 150 PgC), both assessing the amount of SOC in the first 0–3 m, which is also what IPSL-Perm-LandN aims to represent. However, IPSL-Perm-LandN does not represent the carbon stored in Yedoma and Arctic river deltas, missing an additional 327–466 and 96 ± 55 PgC, respectively (Schuur et al., 2022). Deep deposits outside Yedoma or the carbon stored in subsea permafrost are also not represented by the model, but remain challenging to estimate (Schuur et al., 2022; Sayedi et al., 2020).

Figure 9 Carbon stocks at the end of the historical period (2005–2014). (a) Mean vegetation biomass, (b) mean SOC 0–1 m, (c) mean litter biomass and (d) mean SOC 0–3 m for IPSL-Perm-LandN, IPSL-CM6A-LR and other C4MIP models over the tropics (30° S–30° N), mid-latitudes (30–60° S and 30–60° N), the Arctic (>60° N) and the permafrost region. The vegetation biomass observation product is from ESA-CCI and SOC observation products are HWSD, SoilGrids and NCSCD. The other C4MIP models ensemble is composed of CNRM-ESM2-1, CESM2, UKESM1-0-LL, CanESM5, ACCESS-ESM1-5, MIROC-ES2L, NorESM2-LM and MPI-ESM1-2-LR. The error bar shows the full range of C4MIP models.

4.5 Permafrost carbon dynamics

For IPSL-Perm-LandN, the permafrost region is a carbon sink over the historical period, with a net land uptake of 0.32 ± 0.04 PgC yr⁻¹ over the last decade (2005–2014) (Fig. 10). The NBP is higher for IPSL-Perm-LandN than IPSL-CM6A-LR (0.24 ± 0.04 PgC yr⁻¹ for 2005–2014), which is mainly due to differences in the initial state (similar temporal evolution), with IPSL-Perm-LandN being a small carbon sink and IPSL-CM6A-LR a small carbon source. C4MIP models are divided into three groups and show a wide spread due to their differing representations of permafrost and soil carbon processes, as well as due to variations in the permafrost region between models. A first group shows a small land carbon sink, including CESM2 and NorESM2-LM, both of which share the land surface model CLM5, as well as CanESM5 which is known to have a small land NBP (Swart et al., 2019). On the other hand, a second group including UKESM1-0-LL, MIROC-ES2L and CNRM-ESM2-1 has a strong NBP over the permafrost region. IPSL-Perm-LandN belongs to the third group with a moderate permafrost sink, and which includes MPI-ESM1-2-LR, ACCESS-ESM1-5 and IPSL-CM6A-LR. The net land sink simulated in IPSL-Perm-LandN contradicts a recent study based on the upscaling of flux measurements, which concludes that the carbon cycle in the permafrost region is close to neutrality (Ramage et al., 2024). However, the main processes contributing to CO₂ emissions in this study are boreal fires (−0.10 PgC yr⁻¹) and carbon losses from rivers (−0.16 PgC yr⁻¹), two processes that are not represented in IPSL-Perm-LandN. In contrast, Ramage et al. (2024) find boreal forests to be the main contributor to carbon uptake with a net flux of 0.27 PgC yr⁻¹, which is close to the NBP of IPSL-Perm-LandN, although the region considered is slightly different. Therefore, the NBP of IPSL-Perm-LandN could be explained by the lack of important high-latitude CO₂-emitting processes in the model that cannot counterbalance the carbon uptake by boreal forests, which is of the correct order of magnitude. The persistent carbon uptake in IPSL-Perm-LandN leads to an accumulation of land carbon (+17.0 PgC over the historical period, Fig. A19a and b). Most of this carbon enters the soil - especially the slow SOC pool - and is partly buried by cryoturbation. Instead, in boreal areas outside the permafrost region, most of the carbon uptake is stored in vegetation (Fig. A19f). The difference with the permafrost region is particularly striking and is likely due to warmer temperatures, increased soil nitrogen uptake and an abrupt deepening of ALT outside the permafrost region (Fig. A20), and a change in dominant vegetation type (Fig. A21). The increase in land carbon in the permafrost region over the historical period is also found in the majority of C4MIP models, except CESM2 and NorESM2-LM which show a net carbon loss, and UKESM1-0-LL and CanESM5 which show almost no change (Fig. A19d).

Figure 10 Permafrost net land-atmosphere carbon flux (NBP) over the historical period for IPSL-Perm-LandN, IPSL-CM6A-LR and C4MIP models. Positive (resp. negative) values correspond to a land carbon sink (resp. a source). Colored dots represent the mean NBP over the last decade (2005–2014). Black dots correspond to estimates from Ramage et al. (2024). Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR.

The introduction of a vertical discretisation for SOC in IPSL-Perm-LandN allows a better representation of soil carbon dynamics, especially in permafrost soils. The global SOC profile is very close to observations from SoilGrids below 0.5m but shows lower soil carbon in the upper 0.5 m, probably due to overly high turnover rates of the active and slow carbon pools (Fig. A22a). The agreement between IPSL-Perm-LandN and SoilGrids at deeper levels is partly due to the model initialisation by this observation-based product. However, the proportions of active, slow and passive differ from their initial value and vary with depth, while the total SOC concentration remains close to SoilGrids. In general, surface SOC contains a higher proportion of active and slow soil carbon, which tends to decrease with depth and to switch to higher passive carbon fractions. This is also consistent with observations showing older carbon in deeper soil layers because of the time required for SOC to be buried by bioturbation or cryoturbation, leaving only the most stable fraction (Balesdent et al., 2018). In the permafrost region, the SOC vertical profile is flatter than SoilGrids, with less carbon in the upper layers and more at depth (Fig. A22b). As SoilGrids was used to initialise IPSL-Perm-LandN, this simulated SOC profile shows the efficiency of cryoturbation in burying soil carbon and increasing its concentration at depth. Although the absolute SOC concentration is larger, the shape of the SOC profile is closer to NCSCD, which is specifically designed for Arctic regions. Both the observation products and the simulated SOC profile show significant amounts of soil carbon that could lead to carbon emissions as permafrost thaws. In particular, the active and slow SOC fractions are larger in the permafrost region than globally and represent a reservoir of reactive carbon on timescales of days to centuries. Finally, the grid boxes of the model can be grouped by classes (bins) of active layer thickness and the mean SOC profile is calculated for each group (Fig. A22c). For regions of shallow ALT (in purple), the ground remains frozen for most of the year, with only surface layers thawing in summer. In this case, the SOC profile is very similar to its initial value as decomposition is almost non-existent. Conversely, in areas of deep ALT, the soil is mainly unfrozen and the profile is representative of the carbon cycle dynamics of the model. In particular, the profile is flattened compared to the initialisation and shows the effect of cryoturbation, with greater deeper soil SOC concentration. In between, for intermediate ALT, the deep SOC is still close to its initial value while the upper soil responds to the carbon cycle dynamics.

4.6 Ocean carbon cycle

The total net ocean-atmosphere carbon flux (fgco2) of IPSL-Perm-LandN increases slightly until the 1950's and more rapidly thereafter, to reach a mean value of 2.16 ± 0.05 PgC yr⁻¹ over the 2005–2014 period (Fig. A23). This is close to the lower bound of GCB2023 estimates (2.52 ± 0.4 PgC yr⁻¹). fgco2 is also lower than in IPSL-CM6A-LR (2.55 ± 0.16 PgC yr⁻¹ over 2005–2014), due to the effect of the initial state (IPSL-CM6A-LR slightly out of equilibrium with a pre-industrial fgco2 of 0.25 PgC yr⁻¹ compared to 0.045 PgC yr⁻¹ for IPSL-Perm-LandN). This difference is mainly due to the equatorial oceans with larger CO₂ degassing in IPSL-CM6A-LR. The pattern of CO₂ fluxes is consistent with observations (e.g. Fay et al., 2024) with degassing in equatorial ocean and carbon uptake in mid-to high latitudes. Compared to IPSL-CM6A-LR, there is an enhancement of fgco2 pattern in the southern mid- and high-latitudes (i.e. larger uptake in areas of CO₂ uptake and larger release in areas of CO₂ release). Large compensating differences are also evident in the North Atlantic with a reduced carbon sink in the Labrador sea and an increased carbon uptake in the Norwegian and Greenland seas. This is broadly consistent with the observed changes in ocean dynamics in the North Atlantic.

4.7 Compatible CO₂ emissions

After a slow but steady increase from 1850 to 1950, simulated fossil fuel compatible emissions rose much faster during the second half of the 20th century and beyond, reaching 8.3 PgC yr⁻¹ in the last decade (2005–2014) for IPSL-Perm-LandN (Fig. 11a). For both models, they are very close to the fossil fuel emissions diagnosed by the Global Carbon Budget 2023 from different emission datasets, suggesting a relatively accurate simulation of the historical total (land+ocean) carbon sink, except for the simulated plateau in the 1940s. This plateau is due to the stabilisation of the atmospheric CO₂ concentration during this period (Bastos et al., 2016, Fig. 1 and Rubino et al., 2013, Fig. 6b), which led to a decrease of G_ATM, causing a stagnation of E_FF in concentration-driven C4MIP models. However no such stagnation is observed in E_FF estimates from GCB, suggesting that a concomitant increase in carbon sinks occurred during this period (Liddicoat et al., 2021). No C4MIP model represents such an increase, and the dynamics of carbon sinks in this period is still not fully understood (Liddicoat et al., 2021; Bastos et al., 2016). Overall, hypotheses on the origin of this plateau are a decadal variability in the ocean carbon sink not accounted for in reconstructions, a terrestrial sink missing from land surface model estimates, or land-use change processes not included in current datasets (Bastos et al., 2016).

Figure 11 Historical compatible CO₂ emissions. (a) Compatible CO₂ emissions over the historical period for IPSL-Perm-LandN, IPSL-CM6A-LR and GCB2023. Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (b) Cumulative compatible CO₂ emissions over the historical period for IPSL-Perm-LandN, IPSL-CM6A-LR and GCB2023. (c) Total carbon sink for IPSL-Perm-LandN, IPSL-CM6A-LR, C4MIP models and GCB2023. The land carbon sink is plotted on the x axis and the ocean carbon sink is on the y axis. The black line corresponds to the total CO₂ sink estimated by the Global Carbon Budget 2023 (2005–2014). The envelope corresponds to one standard deviation. Plain (resp. empty) circles represent models with (resp. without) an explicit land nitrogen cycle.

The cumulative compatible fossil fuel emissions of IPSL-Perm-LandN from 1850 to 2014 are 406 PgC, which is very close to GCB estimates (404 PgC) (Fig. 11b). Cumulative compatible emissions are overestimated between 1850 and 1950 but the plateau in compatible emissions in the 1940s allows GCB estimate to catch up with IPSL-Perm-LandN. Over the second half of the 20th century and the 21th century (up to 2014), the model is comparable to GCB. This shape is typical of most of C4MIP models with a slowdown of the rate of increase of cumulative emissions in the 1940s and an acceleration from the 1960s onwards (Liddicoat et al., 2021). Cumulative emissions are lower for IPSL-Perm-LandN than for IPSL-CM6A-LR (446 PgC) due to a lower historical total (land+ocean) carbon sink. Most of this difference results from lower compatible emissions from 1850 to 1950 (mostly due to a lower ocean uptake) and from a stronger plateau in the 1940s (due to higher land losses). In principle, the difference in cumulative compatible emissions (EgC) between IPSL-Perm-LandN and IPSL-CM6A-LR could be multiplied by the transient climate response to cumulative emissions (TCRE, °C EgC⁻¹) of IPSL-Perm-LandN to infer the strength of the permafrost carbon-climate feedback (ΔT, °C) (assuming negligible change in the carbon-concentration feedback from the inclusion of permafrost). However, differences between IPSL-Perm-LandN and IPSL-CM6A-LR arise from both the inclusion of new permafrost processes and an explicit nitrogen cycle, leading to superimposed effects in the permafrost region, and to different carbon cycle dynamics in the tropics and mid-latitudes. Therefore, differences in cumulative emissions and TCRE between both versions are not solely due to the inclusion of permafrost in IPSL-Perm-LandN, which prevents a direct assessment of the historical permafrost carbon-climate feedback.

Inter-model differences in compatible emissions arise from the representation of the land and ocean sinks. The total sink of IPSL-Perm-LandN over the last historical decade (3.98 PgC yr⁻¹) is lower than the mean GCB estimate (4.57 PgC yr⁻¹) but within its range of uncertainty, with the land and ocean taking up carbon at a similar rate (Fig. 11c). Compared to IPSL-CM6A-LR, the ocean sink has been reduced while the land sink has increased, resulting in a comparable total sink in the last decade. Overall, ESMs generally underestimate the total carbon sink, either because of low land or ocean carbon sinks, or both. In particular, CanESM5 is known to have a low land carbon sink (Swart et al., 2019) while CNRM-ESM2-1 has a low ocean sink, due to a legacy drift in the net air-sea carbon flux from the spinup. A group of models, including CESM2, UKESM1-0-LL, ACCESS-ESM1-5 and MPI-ESM1-2-LR, has moderate land and ocean sinks, resulting in a slightly underestimated total carbon sink. NorESM2-LM, MIROC-ES2L and IPSL-CM6A-LR are within the range of uncertainty of the total carbon sink from GCB. The reasons for the general underestimation of the total sink by ESMs are very model dependent, but the lack of representation of forest dynamics and demography, the representation of land use change and of the nutrient cycles could explain part of this underestimation (O'Sullivan et al., 2022).

4.8 Limitations of IPSL-Perm-LandN in simulating permafrost ecosystems

Although IPSL-Perm-LandN includes several key permafrost processes, it lacks some important features of high-latitude ecosystems. First, soil hydrology in IPSL-Perm-LandN is limited to a depth of 2 m, and deep water freezing and thawing are based on the water content of the deepest hydrological layer. This can result in unrealistic changes of soil thermal properties associated with water content changes. IPSL-Perm-LandN would benefit from a deeper soil hydrology, particularly in warmer permafrost regions where the active layer can exceed 2 m in depth. In addition, IPSL-Perm-LandN only represents gradual thaw and misses abrupt thaw processes that could be a major source of permafrost carbon loss in the future (Turetsky et al., 2020). Incorporating such processes would require the inclusion of excess ice and soil subsidence in IPSL-Perm-LandN, and could draw on developments made in CLM (Cai et al., 2020; Lee et al., 2014). In IPSL-Perm-LandN, soil organic carbon and nitrogen are vertically resolved but mineral nitrogen is not. Therefore, vegetation can access mineral nitrogen released throughout the soil column, regardless of the depth at which the release occurs. This will impact the future response of the permafrost carbon cycle, as deep nitrogen released at the thaw front will be made directly available for vegetation, possibly leading to overestimated plant nitrogen uptake and productivity. In addition, althoug IPSL-Perm-LandN takes into account their thermal effect, it lacks a comprehensive representation of non-vascular vegetation (e.g. mosses, lichens). They play a critical role in boreal and Arctic ecosystems, regulating soil moisture and accounting for a significant proportion of net primary productivity (Turetsky et al., 2012, 2010). Some moss species can also fix atmospheric nitrogen, providing a nutrient source for other plants, especially in generally nitrogen-limited boreal forests (Markham, 2009). Furthermore, shrubs are also not included in IPSL-Perm-LandN, despite their important physical (lower albedo, shading effect) and biogeochemical (carbon uptake, competition for nutrients and water) impacts particularly in tundra ecosystems. Their interactions with snow are important drivers of the soil thermal dynamics (Domine et al., 2022; Loranty, 2022; Myers-Smith and Hik, 2013). The first attempt to include such high-latitude PFTs in ORCHIDEE was made by Druel et al. (2017) and their complete integration is currently under development. Finally, boreal fires are also a key missing process in IPSL-Perm-LandN that affects permafrost physical properties through immediate ground warming and the burning of insulating vegetation, as well as its biogeochemistry through the combustion of vegetation and soil organic matter. Their absence in IPSL-Perm-LandN is one of the reasons for the overestimated carbon sink during the last decade of the historical simulation compared to Ramage et al. (2024).

5 Conclusions

This work describes IPSL-Perm-LandN, an ESM aiming at better representing the physics and biogeochemistry of high latitudes, and its response to natural and anthropogenic forcings during the historical period. Compared to IPSL-CM6A-LR – the previous version of the model –, the permafrost region has greatly extended and is now close to observations. Soil thermal dynamics has also improved, as shown by the good agreement of the model's active layer thickness with field measurements. Permafrost now holds much larger amounts of soil organic carbon, with a vertical profile close to observations, which is a prerequisite for assessing future permafrost carbon emissions under climate change. In the historical period, the permafrost region is a net carbon sink in IPSL-Perm-LandN, whereas more recent estimates rather show a neutral net flux. However, this is consistent with the processes represented in our model, which does not yet include boreal fires and riverine carbon losses. Overall, the representation of physical and biogeochemical permafrost has greatly improved the response of the model in the Arctic during the historical period. The model developments presented in this study are essential for evaluating potential future permafrost physical and biogeochemical changes. In particular, the vertical discretisation of soil carbon and nitrogen and related soil biogeochemical processes enables the assessment of the permafrost carbon-climate feedback associated with gradual thaw in IPSL-Perm-LandN. Such a feedback analysis under future climate change will be conducted in a forthcoming article. Additionally, an emission-driven version of IPSL-Perm-LandN is under development and will enable the strength of the permafrost carbon-climate feedback to be properly assessed. Most of the new permafrost processes described in this study will be integrated into the IPSL ESM for CMIP7 Fast Track (CMIP7-FT), including the latent heat of soil water phase change, soil insulation by soil carbon and surface organic layers, and the explicit nitrogen cycle. The vertically-resolved soil biogeochemistry will likely only be included for the broader CMIP7 phase, due to the long spinup required and the time constraints of CMIP7-FT. A number of other processes are currently under development, including boreal fire disturbance, peatlands, lake and river biogeochemistry and Arctic vegetation. Medium to long-term developments include the representation of abrupt thaw and associated carbon emissions, excess ice and permafrost small-scale heterogeneity. Collectively, these processes would provide a more comprehensive and realistic picture of future permafrost changes and allow to capture the more complex dynamics of permafrost ecosystems beyond gradual thaw.

Appendix A: Additional description of ORCHIDEE

A1 Carbon assimilation

Carbon assimilation by photosynthesis is based on the scheme proposed by Yin and Struik (2009), which is an extension of the model of Farquhar et al. (1980), developed for C₃ plants. It calculates carbon assimilation as the minimum of the rubisco-limited rate of CO₂ assimilation and the electron-transport-limited rate of CO₂ assimilation. Both the maximum rate of rubisco-limited carboxylation (V_c,max, µmol CO₂ m${}_{\left[\mathrm{leaf}\right]}^{-\mathrm{2}}$ s⁻¹) (i.e. unstressed photosynthetic capacity at optimum temperature) and the maximum rate of electron-transport under saturated light (J_max, µmol e⁻ m${}_{\left[\mathrm{leaf}\right]}^{-\mathrm{2}}$ s⁻¹) follow the formulation of Kattge et al. (2009).

Table A1 PFTs and their dominant locations in ORCHIDEE.

Photosynthetic activity depends on the leaf nitrogen content and can be reduced under nitrogen starvation. Thus, the introduction of an explicit nitrogen cycle allows the model to represent nitrogen limitation of photosynthesis. However, the leaf C:N ratio is dynamic and varies within a limited range as a function of root nitrogen supply and biomass allocation requirements, preventing a strict nitrogen limitation. The C:N ratio of all other vegetation nitrogen pools is determined by the leaf C:N ratio multiplied by a pool-dependent factor.

Table A2 Soil layer structure. Layer node depth (z_n,i), thickness (Δz_i) and depth at layer interface (z_l,i). All in meter.

The dependence of GPP on leaf nitrogen content introduced in IPSL-Perm-LandN replaces the downregulation of maximum photosynthetic capacity as a function of CO₂ used in IPSL-CM6A-LR. In this earlier version of the model, GPP was artificially reduced at high CO₂ concentrations to mimic a nutrient limitation effect. This downregulation mechanism was modeled as a logarithmic function of the CO₂ concentration relative to 380 ppm, following Sellers et al. (1996).

A2 Carbon allocation

The allocation of carbon to the different tissues of the plant (leaves, roots, sapwood, heartwood, and fruits) follows the pipe model theory (Shinozaki et al., 1964), which states that a unit of leaf mass is associated with the downward continuation of non-photosynthetic tissue that has a constant cross-sectional area (Lehnebach et al., 2018). In other words, the production of one unit of leaf mass requires a proportional amount of sapwood to transport water and nutrients from the roots to the leaves, and a proportional amount of roots to take up the water and nutrients from the soil. The allocation scheme dynamically simulates the leaf area depending on the cost of maintaining a unit leaf area, which takes into account the effects of external stresses such as water and nitrogen availability. For instance, more carbon is allocated to roots compared to leaves in case of drought, or nitrogen limitation. The total nitrogen required to sustain the carbon assimilation is then allocated to the different tissues. If nitrogen uptake is insufficient to sustain the carbon uptake, the leaf C:N ratio of the newly growing tissues increases within a certain range. If nitrogen is still deficient, carbon uptake is reduced proportionally to match nitrogen availability. Only the leaf nitrogen concentration is explicitly simulated, while nitrogen is allocated to other tissues in proportion to the leaf nitrogen content.

Table A3 Global equilibrium values and trends of main model variables after spinup. Variables are divided into six sections: land carbon fluxes, land nitrogen fluxes, land carbon stocks, land nitrogen stocks, atmospheric physics and oceanic carbon fluxes. IPSL-CM6A-LR values and trends are presented for comparison. Means and standard deviations are calculated over a 150-year period centered on the start of the historical simulations. Trends are calculated by linear regression (least squares method) over the same period. The years chosen as starting points for the historical simulations are shown for completeness.

Table A4 Equilibrium values and trends of main model variables after spinup in permafrost region. Variables are divided into five sections : land carbon fluxes, land nitrogen fluxes, land carbon stocks, land nitrogen stocks and atmospheric physics. IPSL-CM6A-LR values and trends are presented for comparison. Means and standard deviations are calculated over a 150-year period centered on the start of the historical simulations. Trends are calculated by linear regression (least squares method) over the same period. The years chosen as starting points for the historical simulations are shown for completeness.

A3 Autotrophic respiration

Based on Ruimy et al. (1996), autotrophic respiration is divided into maintenance and growth respiration. Maintenance respiration represents the respiration of the biomass already present, and therefore depends on the amount of biomass of each PFT. It also varies linearly with temperature, with a PFT-dependent coefficient. Maintenance respiration is subtracted from photosynthetic carbon assimilation before allocation, up to a certain threshold (80 % of GPP). If maintenance respiration is higher than this threshold, carbon is taken directly from the tissues. Maintenance respiration also increases with the amount of leaf nitrogen, as in Sitch et al. (2003). In addition, a prescribed fraction of the resulting allocatable carbon (i.e. after maintenance respiration) is lost through growth respiration, which represents the respiration of newly assimilated carbon. The remaining carbon after maintenance and growth respiration (i.e. NPP) is allocated to plant tissues.

A4 Calibration of soil organic matter decomposition and mineral nitrogen losses

When running a spinup under pre-industrial conditions with ORCHIDEE in offline mode, more than 60 % of the initial global soil carbon content (initialised with the SoilGrids product) was lost in 2000 years, with significant losses from the high latitudes. Such low soil organic carbon stocks would lead to a low insulation effect and possible underestimation of soil carbon losses under warming. Therefore the decomposition constant of the passive pool – which contains more than $\mathrm{2}/\mathrm{3}$ of the initial carbon – was decreased by a factor of four.

This did indeed reduced soil carbon losses during spinup but also led to the immobilisation of large amounts of nitrogen, eventually resulting in a strong nitrogen limitation of photosynthesis. Global GPP decreased up to 50 PgC yr⁻¹ under pre-industrial conditions, while it was 95 PgC yr⁻¹ for IPSL-CM6A-LR, and no less than 85 PgC yr⁻¹ for the C4MIP models. To reduce the strength of the nitrogen limitation, we increased the soil mineral nitrogen content available for plant uptake by reducing NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ losses through nitrification and gaseous emissions. This was done by changing the values of the following parameters: N2O_NITRIF_P = 0.0004 gN-N₂O gN-NO${}_{\mathrm{3}}^{-\mathrm{1}}$, NO_NITRIF_P = 0.0016 gN-NO gN-NO${}_{\mathrm{3}}^{-\mathrm{1}}$, CHEMO_0 = 19 (unitless), EMM_FAC = 0.125 (unitless), CTE_BACT = 9 × 10⁻⁵ (unitless) and K_NITRIF = 1.25 d⁻¹. The long turnover times of organic matter prevented a comprehensive statistical optimisation and a manual optimisation of critical parameters had to be performed using a limited number of simulations (with an offline ORCHIDEE configuration). The reduction of mineral nitrogen losses was also motivated by a study showing the overestimation of losses by denitrification in CMIP6 models (Feng et al., 2023).

Overall, with a decreased decomposition constant of the passive carbon and nitrogen pools and reduced mineral nitrogen losses, soil organic carbon and nitrogen losses during the spinup are limited and these pools approach equilibrium faster, while GPP remains close to the value of IPSL-CM6A-LR. Calibration of the model is difficult due to the long turnover times of soil carbon and nitrogen dynamics, and the feedbacks between processes controlling them, and is therefore a source of uncertainty.

Figure A1 Scheme of the integration of surface organic layer thermal properties. Soil thermal properties are modified to include a surface organic layer over a fraction f_SOL of the soil surface.

Figure A2 Schematic of the soil organic carbon dynamics in ORCHIDEE. Red sandglasses correspond to organic carbon decomposition. The red text shows the associated drivers where the indices i refers to the associated pool. Black arrows show internal organic carbon transfers between pools. Blue arrows show CO₂ emissions. Grey text corresponds to the fractions of carbon fluxes that are transferred to another pool or lost as CO₂, and depends on a fixed factor f_j (j corresponds to the associated pool) and soil texture. $\mathrm{L}/\mathrm{N}$ corresponds to the lignin to nitrogen ratio of plant residues. Litter is decomposed into below- and above-ground pools in the model but for clarity, they have been grouped together on this schematic.

Figure A3 Schematic of the soil organic nitrogen in ORCHIDEE. Black arrows show internal organic nitrogen transfers between pools. Purple arrows show exchanges between soil organic and mineral nitrogen pools (immobilisation or mineralisation). f_carbon,i is the corresponding carbon flux between the associated soil organic carbon pools and N:C_X is the N:C ratio of the receiving pool, with X=A, S or P (active, slow or passive). Litter is decomposed into below- and above-ground pools in the model but for clarity, they have been grouped together on this schematic.

Figure A4 Schematic of the soil mineral nitrogen dynamics in ORCHIDEE. Red (resp. green) arrows represent soil mineral nitrogen losses (resp. inputs). Black arrows represent internal land nitrogen transfers. Purple arrows show exchanges between the mineral and organic soil nitrogen pools (immobilisation or mineralisation). The blue arrow shows plant nitrogen uptake. The dashed grey box represents the mineral nitrogen species available for plant uptake (NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$).

Figure A5 Seasonal cycle of historical land surface air temperature. Seasonal cycle (2005–2014) of mean land surface air temperature for IPSL-Perm-LandN, IPSL-CM6A-LR and ERA5 in (a) the Arctic, (b) mid latitudes, and (c) the tropics. Light orange lines represent the three historical members for IPSL-Perm-LandN. Light blue envelopes correspond to one standard deviation between members of IPSL-CM6A-LR. Greenland has been excluded from the Arctic land SAT to only account for non-glaciated land.

Figure A6 Arctic surface air temperature bias of IPSL-Perm-LandN (2005–2014). Mean monthly difference in 2 m air temperature between IPSL-Perm-LandN and ERA5 for the period 2005–2014.

Figure A7 Historical surface temperature over land. (a) Zonal mean of mean land SAT anomaly (2005–2014) relative to 1850–1900 for IPSL-Perm-LandN, IPSL-CM6A-LR, NOAAGlobalTemp and HadCRUT. Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (b) Map of mean land SAT anomaly (2005–2014) relative to 1850–1900 for IPSL-Perm-LandN. Greenland and Antarctica have been excluded for all panels to only account for non-glaciated land.

Figure A8 Historical global surface temperature. (a) Mean global surface air temperature (GSAT) over the historical period for IPSL-Perm-LandN, IPSL-CM6A-LR and ERA5 reanalysis. Colored dots represent the mean GSAT (2005–2014) for IPSL-Perm-LandN, IPSL-CM6A-LR, ERA5 and C4MIP models. Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (b) Mean global GSAT (2005–2014) over the Arctic (>60° N), mid-latitudes (30–60° S and 30–60° N) and the tropics (30° S–30° N) for IPSL-Perm-LandN, IPSL-CM6A-LR and C4MIP models. (c) Anomaly of mean GSAT relative to 1850–1900 for IPSL-Perm-LandN, IPSL-CM6A-LR, NOAAGlobalTemp and HadCRUT reanalyses. Colored dots represent the mean GSAT anomaly (2005–2014) for IPSL-Perm-LandN, IPSL-CM6A-LR, NOAAGlobalTemp, HadCRUT and C4MIP models. Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (d) Mean GSAT anomaly over the Arctic (>60° N), mid-latitudes (30–60° S and 30–60° N) and the tropics (30° S–30° N) for IPSL-Perm-LandN, IPSL-CM6A-LR and C4MIP models. (e) Zonal mean of mean GSAT anomaly (2005–2014) relative to 1850–1900 for IPSL-Perm-LandN, IPSL-CM6A-LR, NOAAGlobalTemp and HadCRUT. Light orange lines represent the three historical members for IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (f) Map of mean GSAT anomaly (2005–2014) relative to 1850–1900 for IPSL-Perm-LandN. The products NOAAGlobalTemp and HadCRUT provide global mean surface temperature (GMST) anomaly, defined as land surface air temperature anomaly over land and sea surface temperature anomaly over the ocean. GMST and GSAT differ by at most 10 % (IPCC AR6 WGI, 2021).

Figure A9 Seasonal cycle of historical precipitation and snowfall. Seasonal cycle (2005–2014) of mean total precipitation for IPSL-Perm-LandN, IPSL-CM6A-LR, ERA5 and MSWEP in (a) the Arctic, (b) mid latitudes, and (c) the tropics. Seasonal cycle (2005–2014) of mean snowfall for IPSL-Perm-LandN, IPSL-CM6A-LR, ERA5 and MSWEP in (d) the Arctic, (e) mid latitudes, and (f) the tropics. Light orange lines represent the three historical members for IPSL-Perm-LandN. Light blue envelopes correspond to one standard deviation between members of IPSL-CM6A-LR.

Figure A10 Snow cover bias of IPSL-Perm-LandN (2005–2014). Annual mean and monthly differences of fractional snow cover (fraction of the ground covered by snow) between IPSL-Perm-LandN and CryoClim for the period 2005–2014. Hatched areas show non significant differences.

Figure A11 Historical ocean physics for IPSL-CM6A-LR. Difference in annual mean sea surface (a) temperature and (b) salinity between IPSL-Perm-LandN the World Ocean Atlas (2005–2014). (c) Mean annual maximum mixed layer depth (2005–2014). (d) Atlantic meridional overturning stream function, on average over 2005–2014.

Figure A12 Historical sea ice. Mean march sea ice fraction (2005–2014) for (a) IPSL-Perm-LandN and (b) IPSL-CM6A-LR. (c) Time series of sea ice extent (total area enclosed within the 15 % sea ice fraction) over the Northern Hemisphere for IPSL-Perm-LandN, IPSL-CM6A-LR and NSIDC observations. The upper and lower curves represent March and September sea ice extents, respectively. Light orange lines represent the three historical members of IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR.

Figure A13 Active layer thickness for CALM and ESA-CCI (2005–2014). Background: map of ALT observation from ESA-CCI (2005–2014). Colored circles: CALM observations.

Figure A14 GPP over the historical period. (a) Mean seasonal cycle (2005–2014) of total GPP over the Arctic (>60° N), the northern hemisphere mid-latitudes (30–60° N), the tropics (30° S–30° N) and the southern hemisphere mid-latitudes (30–60° S). Orange: IPSL-Perm-LandN. Blue: IPSL-CM6A-LR. Plain dark: Jung RS-METEO. Dotted dark: Jung RS. (b) Mean GPP difference between IPSL-Perm-LandN and Jung-RSMETEO product (2005–2014).

Figure A15 Land carbon and nitrogen stocks after the spinup. (a) Total land carbon, including soil, litter and vegetation carbon pools. (b) Soil carbon. Hatching shows the three soil carbon pools (active, slow, passive) for each latitudinal band and permafrost area. (c) Litter carbon. (d) Vegetation C. (e) Total land nitrogen, including soil organic and mineral nitrogen, litter and vegetation. (f) Soil organic nitrogen. Hatching shows the three soil organic nitrogen pools (active, slow, passive) for each latitudinal band and permafrost area. (g) Soil mineral nitrogen. (h) Litter nitrogen. (i) Vegetation nitrogen. Stocks are averaged over 150 years of the piControl simulation surrounding the start years of the historical simulations for IPSL-Perm-LandN, and over the last 150 years of the piControl simulation for IPSL-CM6A-LR. Stocks are given by latitudinal band and over the permafrost area.

Figure A16 Arctic surface air temperature bias of IPSL-CM6A-LR (2005–2014). Mean monthly difference in 2 m air temperature between IPSL-CM6A-LR and ERA5 for the period 2005–2014.

Figure A17 Maps of soil heterotrophic respiration over the historical period. Mean RH difference (2005–2014) between IPSL-Perm-LandN and (a) Warner et al. (2019) and (c) Hashimoto et al. (2015). Mean RH difference (2005–2014) between IPSL-CM6A-LR and (b) Warner et al. (2019) and (d) Hashimoto et al. (2015).

Figure A18 NBP over the historical period. (a) Mean seasonal cycle (2005–2014) of total NBP over the tropics (30° S–30° N), the southern hemisphere mid-latitudes (30–60° S), the northern hemisphere mid-latitudes (30–60° N) and the Arctic (>60° N) for IPSL-Perm-LandN (orange), IPSL-CM6A-LR (blue) and CAMS (2005–2014, dotted dark). (b) Map of IPSL-Perm-LandN NBP (2005–2014). It corresponds to S_LAND−E_LUC in GCB2023. (c) Land carbon sink (sum of NBP and land use change emissions) for IPSL-Perm-LandN (2005–2014). It corresponds to S_LAND in GCB2023.

Figure A19 Changes in permafrost carbon stocks over the historical period. (a) Permafrost cumulative land carbon stocks since 1850 for IPSL-Perm-LandN over the historical period. A positive value corresponds to a land carbon gain. (b) Change in total land C in IPSL-Perm-LandN compared to 1850–1900. A positive value corresponds to a land carbon gain. The red contour shows the limits of the permafrost region in IPSL-Perm-LandN. (c) Change in total litter C in IPSL-Perm-LandN compared to 1850–1900. A positive value corresponds to a carbon gain by the litter. (d) Permafrost cumulative land C stocks (2005–2014) compared to 1850–1900 for IPSL-Perm-LandN and C4MIP models. (e) Change in total SOC in IPSL-Perm-LandN compared to 1850–1900. A positive value corresponds to a SOC gain. (f) Change in total vegetation C in IPSL-Perm-LandN compared to 1850–1900. A positive value corresponds to a carbon gain by the vegetation.

Figure A20 Physical and biogeochemical potential drivers of vegetation change in IPSL-Perm-LandN (2005–2014). (a) Annual maximum active layer thickness, (b) 2 m air temperature, (c) change in 2 m air temperature compared to 1850–1900, (d) soil moisture stress, (e) net soil nitrogen mineralisation and (f) plant nitrogen uptake, averaged over 2005–2014. The red contour shows the limits of the permafrost region in IPSL-Perm-LandN.

Figure A21 PFT coverage of the permafrost region in IPSL-Perm-LandN for the 2005–2014 period. (a) Dominant PFT for each grid cell. See Table A1 for a description of PFTs. (b) Fraction of the grid cell occupied by the dominant PFT. The red contour shows the limits of the permafrost region in IPSL-Perm-LandN. (c) Area and fraction of the permafrost region occupied by each PFT.

Figure A22 SOC historical profile. (a) Mean global SOC profile (2005–2014) for IPSL-Perm-LandN and SoilGrids. Horizontal bars represent the proportion of passive (blue), slow (purple) and active (red) SOC in each soil layer. (b) Mean permafrost SOC profile (2005–2014) for IPSL-Perm-LandN, SoilGrids and NCSCD. (c) Mean permafrost SOC profile (2005–2014) binned by ALT. For all profiles, the first seven soil layers have been averaged.

Figure A23 Net ocean-atmosphere carbon flux over the historical period. (a) Global net ocean-atmosphere carbon flux (fgco2) over the historical period for IPSL-Perm-LandN, IPSL-CM6A-LR, C4MIP models and the Global Carbon Budget 2023. Colored dots represent the mean fgco2 (2005–2014) for IPSL-Perm-LandN, IPSL-CM6A-LR, C4MIP models and GCB2023. Plain (resp. empty) circles represent models with (resp. without) an explicit land nitrogen cycle. Light orange lines represent the three historical members of IPSL-Perm-LandN. The light blue envelope corresponds to one standard deviation between members of IPSL-CM6A-LR. (b) Total fgco2 (2005–2014) over the Arctic (>60° N), mid-latitudes (30–60° S and 30–60° N) and the tropics (30° S–30° N) for IPSL-Perm-LandN, IPSL-CM6A-LR and C4MIP models. (c) Map of IPSL-Perm-LandN mean fgco2 (2005–2014). (d) Difference in mean fgco2 between IPSL-Perm-LandN and IPSL-CM6A-LR.

Figure A24 Ocean and land net carbon fluxes after the spinup. (a) Net sea-air carbon flux, globally and by latitudinal bands. (b) Net land-atmosphere carbon flux (Net Biome Production, NBP), globally, by latitudinal bands and over the permafrost region. Fluxes are averaged over 150 years of the piControl simulation surrounding the start years of the historical simulations for IPSL-Perm-LandN, and over the last 150 years of the piControl simulation for IPSL-CM6A-LR. Positive (resp. negative) fluxes correspond to a land or oceanic carbon sink (resp. source).

Appendix B: Equilibrium state after spinup

In this section, we analyse the equilibrium state reached at the end of the spinup, ensuring that climate and carbon cycle drifts are reasonably small. After about 400 years of coupled spinup, the model is considered to be close enough to equilibrium to start historical simulations. The three historical members were started in years 419, 449 and 479 of the coupled spinup, from different phases of the internal variability of the model. The metrics presented in this section are averaged over 150 years surrounding the start years of the historical simulations, to look for potential drifts in the pre-industrial state that could affect these simulations. A detailed description of the model state after spinup can be found in Table A3.

The global mean surface temperature (GMST) after spinup is 12.28 ± 0.12 °C (mean ± std over the three simulation members) and is at equilibrium (trend of +0.0003 °C yr⁻¹). IPSL-Perm-LandN is slightly colder than the IPSL-CM6A-LR piControl simulation used for CMIP6, which has a GMST of 12.54 ± 0.12 °C, and this is consistent at all latitudes. No significant regional trends are observed, indicating that equilibrium is reached everywhere (not shown). In the permafrost region, the mean temperature is −12.10 ± 0.33 °C, significantly colder than IPSL-CM6A-LR (−10.30 ± 0.51) °C. Global mean precipitation (2.95 ± 0.01 mm d⁻¹) and snowfall (0.260 ± 0.004 mm d⁻¹) also show negligible trends.

The global net air-sea carbon flux fgco2 is 0.045 ± 0.086 PgC yr⁻¹ (Fig. A24a). This positive value corresponds to a small remaining oceanic carbon sink after the spinup. Although full equilibrium is not reached, the net air-sea carbon flux is much lower than for IPSL-CM6A-LR (0.25 ± 0.09 PgC yr⁻¹). The resulting bias in the historical period simulations is consequently one order of magnitude lower than IPSL-CM6A-LR. However this well-balanced global net air-sea flux masks regional variability, with an oceanic carbon source in the tropics (−1.09 ± 0.04 PgC yr⁻¹) that is counterbalanced by carbon sinks in mid- and high-latitudes (1.13 ± 0.07 and 0.01 ± 0.02 PgC yr⁻¹ respectively). This is an expected behaviour as the large-scale oceanic circulation induces CO₂ outgassing in the tropics, and an oceanic CO₂ sink in cooling poleward flowing subtropical surface waters as well as in equatorward flowing subpolar surface waters.

The global net land-atmosphere carbon flux (NBP) after the spinup is 0.038 ± 0.510 PgC yr⁻¹, which is equal to 1 % of the present-day land carbon sink (Friedlingstein et al., 2023) (Fig. A24b). Over the historical period (1850–2014), this drift is responsible for a cumulative land carbon accumulation of 6.3 PgC, which is negligible compared to land carbon changes during this period. On the contrary, in the IPSL-CM6A-LR piControl simulation, the land is a carbon source with a negative NBP of −0.19 ± 0.64 PgC yr⁻¹. The absolute value of NBP is an order of magnitude smaller in IPSL-Perm-LandN, indicating that the model is closer to equilibrium. In addition, the net land-atmosphere carbon exchange is close to equilibrium at all latitudes with a small carbon sources in the tropics (−0.029 ± 0.440 PgC yr⁻¹) and small carbon sinks at mid- (0.044 ± 0.250 PgC yr⁻¹) and high latitudes (0.023 ± 0.070 PgC yr⁻¹). The permafrost region is also a small sink, with a net carbon flux of 0.034 ± 0.48 PgC yr⁻¹. Approximately three quarters of the remaining imbalance in the global net land carbon flux is due to a drift in soil carbon cSoil and one quarter to a drift in vegetation carbon cVeg, both slowly increasing over time. The positive drift in total soil carbon results from opposite trends in the model carbon pools, with the active and slow soil carbon pools gaining carbon (resp. +0.004 and +0.06 PgC yr⁻¹) and the passive pool losing carbon (−0.03 PgC yr⁻¹). Soil organic nitrogen trends are similar for individual pools with the active and slow pools gaining nitrogen (resp. +0.0002 and +0.001 PgN yr⁻¹) and the passive pool losing nitrogen (−0.02 PgN yr⁻¹), but resulting in a net soil organic nitrogen loss of −0.02 PgN yr⁻¹.

At the end of the spinup, the land stores 3567 PgC, distributed into 485 PgC in vegetation, 100 PgC in the litter and 2982 PgC in the soil (Fig. A15). Most of the vegetation biomass is found in the tropics (302 PgC), followed by mid-latitudes (147 PgC) and the Arctic (37 PgC). Soil organic carbon is divided into 641 PgC in the tropics, 1364 PgC in mid-latitudes and 845 PgC in high latitudes. Most of the soil carbon is found in the so-called passive pool (85 %), the rest being stored in the slow pool (14 %) and a tiny fraction in the active pool (1 %). Compared to IPSL-CM6A-LR, which has a downregulation of GPP with CO₂, litter and vegetation stocks are slightly higher but the main difference is the almost 6-fold increase in soil carbon, especially in mid- and high latitudes. In particular, permafrost soil carbon was almost non-existent in IPSL-CM6A-LR and now amounts to 1006 PgC. Total land nitrogen is 225 PgN, most of which is stored in the soil in organic form (215 PgN). In addition, the vegetation contains 8 PgN, the litter 1 PgN and the soil also stores 1 PgN of mineral nitrogen, which cannot be compared to IPSL-CM6A-LR which did not include a representation of the nitrogen cycle.