# Enhanced prescription of soil organic and mineral content in the ORCHIDEE LSM to better simulate soil temperatures: application at nine high-latitude GEM and FLUXNET sites

Amelie Cuynet, Elodie Salmon, Efrén Lopéz-Blanco, Mathias Goeckede, Hiroki Ikawa, Hideki Kobayashi, Annalea Lohila, Catherine Ottlé

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

## ► To cite this version:

Amelie Cuynet, Elodie Salmon, Efrén Lopéz-Blanco, Mathias Goeckede, Hiroki Ikawa, et al. Enhanced prescription of soil organic and mineral content in the ORCHIDEE LSM to better simulate soil temperatures: application at nine high-latitude GEM and FLUXNET sites. Journal of Geophysical Research: Biogeosciences, 2025, 130, pp.e2025JG008776. hal-05059168

HAL Id: hal-05059168

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

Submitted on 17 Dec 2025

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

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

<figure><img src="image_2.png" /></figure>
# JGR Biogeosciences

## RESEARCH ARTICLE

10.1029/2025JG008776

**Key Points:**

* A new approach in the ORCHIDEE LSM incorporates pixelwise and depth-specific soil data from SoilGrids to calculate soil thermal properties

* Improving the distribution of organic and mineral matter within soils enhances the simulation of thermal dynamics

* Modeled soil temperatures show improved alignment with observed data, especially at depths between 0.2 and 1 m

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

A. Cuynet,
amelie.cuynet@lsce.ipsl.fr

**Citation:**

Cuynet, A., Salmon, E., López-Blanco, E., Goeckede, M., Ikawa, H., Kobayashi, H., et al. (2025). Enhanced prescription of soil organic and mineral content in the ORCHIDEE LSM to better simulate soil temperatures: Application at nine high-latitude GEM and FLUXNET sites. *Journal of Geophysical Research: Biogeosciences*, 130, e2025JG008776. https://doi.org/10.1029/2025JG008776

Received 16 JAN 2025
Accepted 22 NOV 2025

**Author Contributions:**

**Conceptualization:** Amélie Cuynet,
Catherine Ottlé

**Data curation:** Amélie Cuynet

**Formal analysis:** Amélie Cuynet,
Elodie Salmon, Catherine Ottlé

**Funding acquisition:** Catherine Ottlé

**Investigation:** Amélie Cuynet

**Methodology:** Amélie Cuynet

**Project administration:** Catherine Ottlé

**Resources:** Amélie Cuynet, Efrén López-Blanco, Mathias Goeckede, Hiroki Ikawa, Hideki Kobayashi, Annalea Lohila

**Software:** Amélie Cuynet

© 2025. The Author(s).
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

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

## Enhanced Prescription of Soil Organic and Mineral Content in the ORCHIDEE LSM to Better Simulate Soil Temperatures: Application at Nine High-Latitude GEM and FLUXNET Sites

Amélie Cuynet¹ ID, Elodie Salmon¹ ID, Efrén López-Blanco²,³, Mathias Goeckede⁴ ID,
Hiroki Ikawa⁵ ID, Hideki Kobayashi⁶, Annalea Lohila⁷ ID, and Catherine Ottlé¹ ID

¹LSCE-IPSL (CNRS-CEA-UVSQ), Université Paris-Saclay, Gif-sur-Yvette, France, ²Department of Ecoscience, Arctic Research Center, Aarhus University, Roskilde, Denmark, ³Department of Environment and Minerals, Greenland Institute of Natural Resources, Nuuk, Greenland, ⁴Max Planck Institute for Biogeochemistry, Jena, Germany, ⁵Hokkaido Agricultural Research Center, National Agriculture and Food Research Organization, Sapporo, Japan, ⁶Research Institute for Global Change, Japan Agency for Marine-Earth Science and Technology, Yokohama, Japan, ⁷Finnish Meteorological Institute, Climate System Research, Helsinki, Finland

**Abstract** Over the past two decades, numerous studies have emphasized the importance of including organic matter (OM) in land surface models (LSMs) to accurately represent soil thermal and hydrological properties. This is particularly relevant in Arctic regions, where organic-rich soils are widespread. Consequently, most LSMs incorporate parameterizations that account for OM effects, although these implementations are often simplified. Recent advancements in global soil data sets now enable more precise modeling of soil properties by providing detailed inputs for soil composition and physical characteristics. This study focuses on the refinement of the representation of soil organic and mineral content and the revision of the parameterizations of heat capacity, thermal conductivity, and porosity in the ORCHIDEE LSM, using data from the SoilGrids 250m v2.0 database. The updated model is evaluated across multiple Arctic and boreal sites and compared against two earlier versions: (a) a Bulk version that neglects OM effects on the thermal processes and (b) a simplified version with a basic OM prescription. Results show that incorporating OM into thermal processes modeling significantly improves soil temperature simulations, particularly under the soil surface in the critical zone. For some sites, root mean square errors (RMSE) are reduced by up to 50% compared to the Bulk version, especially during the snow-free summer months. These findings highlight the value of high-resolution soil data sets, such as SoilGrids, for improving simulations of thermal dynamics in carbon-rich Arctic soils.

**Plain Language Summary** Soil plays a critical role in shaping the Earth's climate by exchanging energy, water, and carbon with the atmosphere. Its behavior depends heavily on its chemical composition, particularly when it is rich in organic matter, which influences heat and water retention. Land surface models (LSMs) simulate soil's energy, water, and carbon budgets, but many don't fully capture the impact of organic matter on soil properties, leading to less accurate predictions. In this study, we improved the ORCHIDEE LSM by using detailed global soil data from a high-resolution database called SoilGrids. This allowed us to better prescribe soil organic and mineral content and to refine the modeling of heat transfers and soil temperatures. We tested this updated model at various instrumented Arctic and boreal sites and compared it to earlier versions. Results showed that including detailed organic matter data improved soil temperature predictions, especially in deeper layers and during snow-free periods, increasing accuracy by up to 50% in some cases. This underscores the value of precise soil data for better climate predictions, particularly in Arctic and boreal regions where organic-rich soils play a key role in the ecosystem.

## 1. Introduction

The importance of accounting for soil Organic Matter (OM) to correctly represent soil thermal processes in land surface models (LSMs) has been shown by previous studies (Koven et al., 2013; Lawrence & Slater, 2008; Rinke et al., 2008). The thermal properties of OM largely differ from those of the other soil components, that is mineral matter (MM), air or water. OM has a lower thermal conductivity associated with a relatively higher heat capacity
**Supervision:** Elodie Salmon,
Catherine Ottlé

**Visualization:** Amélie Cuynet

**Writing – original draft:** Amélie Cuynet,
Elodie Salmon, Catherine Ottlé

**Writing – review & editing:**
Amélie Cuynet, Elodie Salmon,
Efrén López-Blanco, Mathias Goeckede,
Hiroki Ikawa, Hideki Kobayashi,
Annalea Lohila, Catherine Ottlé

compared to mineral matter (MM). It means that in organic rich soils composed of a large proportion of OM, the
thermal insulation is greater than in mineral soils (Boelter, 1969; Letts et al., 2000).

OM is composed of Soil Organic Carbon (SOC) and includes other elements such as hydrogen, oxygen and nitrogen. It is present in the soils of most regions of the world in various amounts, mostly depending on the land cover type and use and climate history. In particular, soils in the Arctic and boreal regions can be particularly rich in OM, since low soil temperatures diminish the rate of decomposition of the biomass, thus enabling the aggregation of peat deposits (Decharme et al., 2016). That leads to the accumulation of large quantities of OM in the soils (Lawrence & Slater, 2008). Moreover, permafrost covers around 14 million km² of land in the northern hemisphere (Hugelius et al., 2024; Obu, 2021) and contains vast quantities of carbon that were estimated to be approximately 1,350 Pg (Hugelius et al., 2014), while northern latitude peatlands alone contain more than 400 Pg of carbon (Hugelius et al., 2020). The Arctic is warming at a rate up to four times larger than the rest of the globe since 1979 (Rantanen et al., 2022), which may lead to a large release of the stocked carbon as greenhouse gases into the atmosphere, due to permafrost thawing and peatlands degradation among other processes (Chadburn et al., 2022; Miner et al., 2022). To be able to correctly predict the future evolution of surface temperature and energy transfers to the atmosphere in the present context of global warming and permafrost thawing, it is essential to consider OM, especially in the Arctic.

The inclusion of OM content in the calculation of soil thermal properties in LSMs tends to have a more pronounced impact on soil temperatures during summer, when the soil is directly exposed to solar and atmospheric radiations. During winter, snow acts as a dominant insulating layer at the land surface by separating the soil from the atmosphere (Paquin & Sushama, 2015). As a result, the contribution of OM to insulation is less noticeable in winter, with snow primarily controlling the soil's thermal dynamics. However, in summer, adding OM effects to models that already account for snow processes prevents excessive soil warming, and generally leads to cooler simulated soil temperatures year-round (Bonan & Shugart, 1989; Lawrence & Slater, 2008). In addition, the porosity of organic rich soils can be more than twice as large as the porosity of mineral matter, enabling more water (liquid or solid) to be stored, which impacts the thermal properties as well (Lawrence & Slater, 2008).

Before the 2000s, in the soil science community, very few studies have focused on the combination of OM and mineral matter in soils to determine their thermal and hydraulic properties such as thermal and hydraulic conductivity, heat capacity or porosity (Campbell, 1985; Farouki, 1981), the other studies only assessing the properties of mineral matter or organic-rich soil separately. From 2000, the effects of OM on thermal and water dynamics have begun to be implemented in some LSMs, such as the Canadian Land Surface Scheme (CLASS) (Letts et al., 2000) or the National Center for Atmospheric Research Land Surface Model (NCAR LSM) (Beringer et al., 2001), where the authors considered layers of rich organic soil with different physical properties above mineral soil layers. Finally, the impact of OM on thermal processes started to be more commonly included in LSMs with other permafrost-related processes, while the importance of dealing with OM in LSMs was increasingly discussed in the literature (Koven et al., 2013; Lawrence & Slater, 2008; Rinke et al., 2008). OM was notably implemented in the Community Land Model (CLM) (Lawrence & Slater, 2008; Tran et al., 2017), ORCHIDEE (Gaillard et al., 2025; Guimberteau et al., 2018; Koven et al., 2009), the Joint UK Land Environment Simulator (JULES) (Chadburn et al., 2015), ISBA-SURFEX (Decharme et al., 2016) or the NOAH-MP LSMs (Sun et al., 2023). In most of the aforementioned models, improvements in the simulation results were observed when implementing the OM impacts on soil thermics.

Although thermal parameterizations vary between the different models, properties that depend on the OM content are all defined using the SOC volume fractions in the different model soil layers. The proportions of OM in each layer, which are used as weights to calculate the physical properties of the soil, are estimated with two different methods in the previous list of models: (a) a simple ratio between the SOC density and the bulk density of organic rich soil samples; (b) a more complex ratio integrating the volume occupied by coarse matter (Chen et al., 2012; Yuan et al., 2021).

Although this SOC fraction can be a good proxy of OM and therefore its physical properties, this sole information is not sufficient to completely describe the soil's composition and its properties, but this method is convenient and efficient when soil data is not easily available and reliable. The calculation of thermal parameters is influenced by the determination of other soil physical properties, such as bulk density and the amount of coarse and organic
matter. This influence has been demonstrated in several studies focusing, for example, on the estimation of soil thermal conductivities (Balland & Arp, 2005; He et al., 2021).

To calculate SOC fractions, most models rely on external databases that provide SOC densities and other soil properties at local or global scales, and usually with an in-depth discretization, such as the Harmonized World Soil Database (HWSD) (Decharme et al., 2016; FAO et al., 2012), the Northern Circumpolar Soil Carbon Database (NCSCD) (Chadburn et al., 2015; Hugelius et al., 2013) or more recently SoilGrids (Gaillard et al., 2025; Hengl et al., 2014). The availability of more site-level data, combined with the development of machine learning methods and other state-of-the-art approaches, led to a rapid evolution of soil information and it is now possible to have access to many more soil properties at very fine horizontal resolutions, multiple depths and global scale, like in SoilGrids250m v2.0 (Poggio et al., 2021), including the ones previously listed.

In the ORCHIDEE standard version, the impact of OM on the soil thermal processes is accounted for using in-depth SOC fractions prescribed from a fixed map derived from the previous version of SoilGrids, SoilGrids1km (Hengl et al., 2014). In this work, we want to refine the calculation of the thermal parameters by using an updated version of SoilGrids allowing us to model the impacts of OM on the thermal processes in a more rigorous way.

Based on these elements, this study has three main goals:

1. Improve the way to prescribe soil composition (OM, mineral matter, water, air...) and exploit the global soil database SoilGrids250m v2.0 in ORCHIDEE;

2. Improve the calculation of the soil thermal properties—heat capacity and thermal conductivities—by using the prescribed soil composition;

3. Assess the impact of these changes at boreal and arctic carbon-rich sites by evaluating ORCHIDEE simulated soil temperatures to local observations.

A brief presentation of the ORCHIDEE LSM and the model developments are presented in Section 2; the data and methods used for the study are respectively described in Sections 3 and 4; model evaluation on a selection of instrumented sites and comparison to the standard version is shown in Section 5, with a focus on the Zackenberg site (Greenland). Finally, discussion and conclusions are provided in Sections 6 and 7.

## 2. New Specification of the Soil Organic and Mineral Matter in ORCHIDEE

### 2.1. The ORCHIDEE Model

ORCHIDEE is the land surface model, part of the Earth System Model (ESM) of the Institute Pierre-Simon Laplace (IPSL). In this study, the ORCHIDEE-V3.0 version was chosen since it integrates recent developments around the modeling of the carbon cycle compared to the V2.0 version used for the Coupled Model Intercomparison Project Phase 6 (CMIP6) (Cheruy et al., 2020). The model simulates the energy and water transfers in the soil-atmosphere continuum and at the surface-atmosphere interface, as well as the carbon and nitrogen cycles and their interactions (Vuichard et al., 2019). Vegetation processes are parameterized for 15 different Plant Functional Types (PFTs) related to their phenology, leaf type, physiological activity (e.g., C3/C4 photosynthetic pathways) and climate (Box, 1996). Soil mineral composition is defined at the grid-cell level, based on the 12-classes USDA soil classification (Forbes et al., 1987). While energy and snow processes are computed at the grid-cell scale (one energy budget for the whole cell), water processes are resolved by accounting for the PFTs present in the grid cell, with a maximum of three soil columns separating bare soil from high (e.g., forests) and low (e.g., grasslands) vegetations. Carbon stocks and fluxes are, for their part, resolved for each PFT present within the grid cell. Different vertical grids are considered depending on the biogeophysical cycles considered. Soil hydrology and thermics share the same grid consisting of 11 layers from the surface and down to 2 m (default) with a geometrically increasing internode spacing. Thermal processes are resolved deeper, up to a maximum depth of 90 m, with 7 additional layers. Finally, for carbon processes, two options are offered: the first one based on a bucket scheme, where a bulk SOC content is modeled but does not impact the thermal processes which depend only on the soil mineral content and water amount, and a second option where the SOC content derived by SoilGrids (Hengl et al., 2014) is discretized at the same levels as the thermal processes and is used to calculate soil properties and thermal fluxes. In all versions, soil freezing and snow thermo-hydric processes are modeled in the same way, and SOC is not accounted for in the calculation of soil hydric transfers.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Soil composition in organic matter (OM), coarse matter (CM) and fine mineral matter (MM): definition of the volume fractions $\psi$, solid volume fractions $F$ and soil porosity $\theta_s$.</figcaption></figure>

ORCHIDEE can be run at various scales ranging from local to global. For the present study, the model requires the prescription of the atmospheric conditions (air pressure, temperature and humidity, wind speed and incoming radiation), a description of the land surface (fraction of the different PFTs, dominant soil texture) and the initialization of the energy, water, carbon and nitrogen stocks, as well as soil temperature. These conditions are derived from global data sets (i.e., atmospheric reanalysis, land cover and soil texture classifications) or from in situ measurements. In this work, the model was used at the site level and forced with in situ and atmospheric reanalysis data (described in Sections 3.1.1 and 3.1.3).

## 2.2. Soil Thermal Processes

The physical and chemical properties of the soil, such as soil moisture, porosity, SOC and mineral contents, influence its thermal properties. Indeed, soil thermal conductivities and heat capacities are sensitive to soil moisture and to the organic content since the air thermal conductivity ($0.025 \text{ W K}^{-1} \text{ m}^{-1}$) is more than 20 times lower than that of liquid water and similarly the solid thermal conductivity of mineral matter is usually more than 10 times higher than that of OM (see Table S3 in Supporting Information S1). Moreover, ice is nearly four times more conductive than liquid water. In a soil that combines all these components, resulting thermal conductivities may thus drastically change over space and time. Therefore, considering these physico-chemical properties in land ecosystem models is critical to simulate spatial and temporal soil thermal variability.

### 2.2.1. Soil Composition Definitions

To facilitate the understanding of the following section, Table S1 in Supporting Information S1 is provided in Supplementary, referencing all the symbols and parameters described here.

Soil is composed of a mixture of OM, fine mineral matter (MM) and coarse matter (CM). The soils being porous, they are also filled with air and water (solid or liquid). Figure 1 illustrates the composition of a soil volume unit. Volume fractions, corresponding to each type of matter and its porosity, are defined by $\psi_{OM}$, $\psi_{MM}$ and $\psi_{CM}$ for OM, MM and CM respectively. Other fractions, hereafter referred to as “solid volume fractions,” represent the effective volume occupied by the matter, excluding pores, and are defined by $F_{OM}$, $F_{CM}$ and $F_{MM}$. The remaining volume corresponds to the porosity ($\theta_s$) which represents the addition of the volumes occupied by water (liquid $F_{water}$ and solid $F_{ice}$) and air $F_{air}$.

### 2.2.2. Soil Representation in ORCHIDEE

In the ORCHIDEE default version, named Bulk hereafter, the bucket carbon module is used and the thermal and hydrological processes are not constrained with SOC. For each pixel, the non porous part of the soil has the thermal properties of its corresponding USDA soil texture class, that is prescribed by a fixed texture map (Forbes et al., 1987). The properties are provided in Tables S2 and S3 in Supporting Information S1.

In a second option, named DisC hereafter, SOC is used to constrain soil thermal properties, including dry and solid thermal conductivities, heat capacity, and porosity. This version stems from developments initially introduced in ORCHIDEE MICT (Guimberteau et al., 2018), which aimed to better represent the influence of SOC on vertically discretized soil thermal dynamics, particularly in high-latitude regions. This configuration is a simple
adaptation of the Bulk version by using the fraction of soil organic carbon to compute soil thermal properties. For each of the 18 soil layers i, this SOC fraction is estimated using Equation 1:

$$ \psi_{SOC,i} = \frac{\rho_{SOC,i}}{\rho_{SOC,max}} \qquad (1) $$

where $\rho_{SOC,i}$ is the Soil Organic Carbon content in an elementary volume of soil in layer i (kg C.m⁻³) and $\rho_{SOC,max}$ = 130 kg C.m⁻³ is assumed to be the maximum SOC content possible and is equal to the value of a typical soil carbon bulk density of peat (Guimberteau et al., 2018; Lawrence & Slater, 2008). This fraction can be interpreted as the total volume occupied by the soil organic carbon including its porosity within the matter in an elementary volume of soil. The values of $\rho_{SOC,i}$ can be those calculated in the model or can be derived from a fixed map, provided as an input to the model.

This method has several advantages: its simplicity, and the possibility to directly use the SOC quantities calculated by the model or derived from an external map. However, it relies on spurious hypotheses. Indeed, this method overlooks the fact that OM and SOC, although related, are not exactly the same, since OM is also constituted of other organic components such as hydrogen and nitrogen. The value $\rho_{SOC,max}$ from Equation 1 is derived from the OM density approximate value defined by Farouki (1981) to determine a theoretical maximum SOC density of organic rich soil by assuming a porosity of 0.9 and the SOC representing 100% of the solid soil's mass. The model then calculates some physical properties (porosity, heat capacity, dry and solid thermal conductivities) by performing weighted means of the properties of organic and mineral matter, the weights being respectively, the volume fraction of SOC $\psi_{SOC,i}$ and the assumed remaining soil volume $1 - \psi_{SOC,i}$. However, it would make more sense in terms of physics to calculate a fraction of OM $\psi_{OM,i}$ instead of a fraction of SOC $\psi_{SOC,i}$ since the physical properties used are those of OM and not SOC. To do so, we can replace $\rho_{SOC,i}$ by $\rho_{OM,i}$, which would be the effective mass of soil organic matter in the elementary volume. However, this formulation does not account for the soil compaction, and assumes that soil porosity is constant, while the physical properties of OM depend on its degree of decomposition, its structure and particle size, affecting the porosity (Huat et al., 2011; Robinson et al., 2022).

To correct this issue, we developed a new formulation, named NewDisC in the following, where we calculate the thermal properties of the soil according to the volume fractions of all components, taking advantage of the availability of high-resolution maps of soil properties—in particular soil bulk densities, carbon contents and coarse matter repartition. The SoilGrids250m v2.0 database (Poggio et al., 2021) provides such information globally.

Solid volume fractions are not directly available in the SoilGrids database and are therefore derived from other variables, namely the bulk density of fine matter (FM) $BD_{FM}$ [kg/m³] (FM here includes OM and MM while excluding CM), the concentration (or mass fraction) of SOC in the FM $f_{SOC}$ [kg/kg], and the volume fraction of CM $\psi_{CM}$ [m³/m³]. More details regarding the preparation of these variables are provided in Section 4.1.1.

It is possible to estimate the mass fraction of OM as a function of that of SOC. To determine the amount of OM from SOC, it is common to use a conversion factor $\alpha_{VB}$ called Van Bemmelen's factor (Heaton et al., 2016). Although its value can change depending on the soil types, a value of 0.58 has been chosen, as suggested in the literature (Decharme et al., 2019; Heaton et al., 2016; Post et al., 1982; Sun et al., 2016). The mass fraction of organic matter $f_{OM}$ can then easily be calculated from the mass fraction of SOC $f_{SOC}$ Equation 2:

$$ f_{OM} = \frac{f_{SOC}}{\alpha_{VB}} \qquad (2) $$

With that information, constant values such as the particle densities of OM and mineral matter, respectively $\rho_{OM}$ and $\rho_{MM}$ (see Table S2 in Supporting Information S1) and the porosity of mineral matter that depends on the soil texture (see Table S1 in Supporting Information S1) are also needed to compute the solid volume fractions.

$$ F_{OM} = \frac{V_{OM}}{V_{tot}} = f_{OM} \frac{BD_{FM}}{\rho_{OM}} (1 - \psi_{CM}) \qquad (3) $$
<figure><img src="image_5.png" /><figcaption>Figure 2. Location of the chosen sites.</figcaption></figure>

$$F_{MM} = \frac{V_{MM}}{V_{tot}} = (1 - f_{OM}) \frac{BD_{FM}}{\rho_{MM}} (1 - \psi_{CM}) \quad (4)$$

Assuming that the porosity of the coarse matter $\theta_{s,MM}$ in a given volume is similar to that of the fine mineral matter and is constant in a pixel, it is possible to calculate the solid volume fraction of coarse matter.

$$F_{CM} = \psi_{CM} (1 - \theta_{s,MM}) \quad (5)$$

These fractions will be used in the following sections, especially in order to calculate the thermal parameters such as capacity and conductivity, as well as the porosity in the thermal module.

The calculation of the fractions related to water, ice and air are obtained from the resolution of the soil water budget which was not modified in this study. More details can be found in De Rosnay et al. (2002).

### 2.2.3. Thermal Properties Modeling

The evolution of soil temperatures in ORCHIDEE and the calculation of the thermal properties in the current version are reported in Appendix A.

In NewDisC, the representation of soil thermal properties in ORCHIDEE is revised to improve physical consistency and incorporate detailed soil composition. Table A2 summarizes the modified equations compared to their previous version. Specifically, we modify the calculation of:

* Porosity $\theta_s$—now computed directly from solid volume fractions, rather than using an assumed average Equation A7.

* Heat capacity $C_{tot,i}$—includes contributions from air and adjusts the weighting of solid components using explicit volume fractions Equation A10.

* Solid thermal conductivity $K_{solid,i}$—revised to consider only the solid volume fractions instead of the porous volume Equation A12.

* Dry thermal conductivity $K_{dry,i}$—where the volume fraction of OM replaces the SOC fraction Equation A15.

## 3. Data

### 3.1. Databases

The evaluation of the new parameterization on ORCHIDEE simulations is performed by doing site simulations and comparing the results to in situ measurements. For that purpose, multiple databases are used and presented in this section.

#### 3.1.1. In Situ Data

This study evaluates the new developments at local scale across high-latitude sites equipped with thermal sensors to determine the impact of improved soil representation in cold environments, where soils often contain high levels of organic carbon (McGuire et al., 2009). The study sites are chosen accordingly: over 60°N, and well distributed geographically, as shown in Figure 2, with two sites in Greenland, two in Russia, four in Finland and one in the USA. Table 1 summarizes the characteristics of each site.

The Greenland sites (Zackenberg and Kobbefjord) data is issued from two sub-programs (GeoBasis and ClimateBasis) freely accessible from the GEM database (https://data.g-e-m.dk/ (Christensen et al., 2017)). The meteorological data sets include air temperature (°C), pressure (hPa), longwave and shortwave incoming radiation (W m⁻²), relative humidity (%), wind speed (m s⁻¹), accumulated precipitation (mm) and snow depth (m). Moreover, soil temperatures at different depths (°C) are also provided.

The FLUXNET2015 data set provides numerous measurements at the ecosystem scale, from more than 200 sites around the globe (Pastorello et al., 2020). For the purpose of this study, only a small sample of these sites has been
<table><caption>**Table 1**
Site Description</caption><thead><tr><th>Site</th><th>GL-Zaf</th><th>GL-Nuf</th><th>FI-Hyy</th><th>FI-Jok</th><th>FI-Lom</th><th>FI-Sod</th><th>RU-Cok</th><th>USA</th></tr></thead><tbody><tr><td>Location</td><td>Greenland</td><td>Greenland</td><td>Finland</td><td>Finland</td><td>Finland</td><td>Finland</td><td>Russia</td><td>Poker Flat</td></tr><tr><td>Name</td><td>Zackenberg</td><td>Nuuk</td><td>Hyytiälä</td><td>Jokioinen</td><td>Lompolojankka</td><td>Sodankyla</td><td>Chokurdakh</td><td>Research Range</td></tr><tr><td>Longitude</td><td>20°30′W</td><td>51°21′W</td><td>24.29°</td><td>23.51°</td><td>24.21°</td><td>26.64°</td><td>161.34°</td><td>147.49°</td></tr><tr><td>Latitude</td><td>74°30′N</td><td>64°07′N</td><td>61.85°</td><td>60.90°</td><td>68.00°</td><td>67.36°</td><td>68.61°</td><td>70.83°</td></tr><tr><td>Years</td><td>2009–2017</td><td>2007–2018</td><td>1996–2014</td><td>2000–2003</td><td>2007–2009</td><td>2001–2014</td><td>2002–2004</td><td>2003–2009</td></tr><tr><td>Vegetation/ecosystem</td><td>Sedge dominated fen (Eriophorum scheuchzeri and Dupontia psilosantha and abundant moss layer (Sanionia uncinata and Sphagnum lindbergii) (López-Blanco et al., 2020). Continuous permafrost, active layer depth: 0.5–1 m (Lund et al., 2014; López-Blanco et al., 2020)</td><td>Sedge dominated fen (Eriophorum angustifolium and Scirpus caespitosus and abundant moss layer (Sanionia uncinata and Sphagnum lindbergii) (López-Blanco et al., 2020)</td><td>Evergreen needleleaf Forests: Lands dominated by woody vegetation with a percent cover >60% and height exceeding 2 m. Almost all trees remain green all year. Canopy is never without green foliage</td><td>Croplands: Lands covered with temporary crops followed by harvest and a bare soil period (e.g., single and multiple cropping systems). Note that perennial woody crops will be classified as the appropriate forest or shrub land cover type</td><td>Permanent Wetlands: Lands with a permanent mixture of water and herbaceous or woody vegetation that cover extensive areas. The vegetation can be present in either salt, brackish, or fresh water. Sedge fen, no permafrost</td><td>Evergreen Needleleaf Forests: Lands dominated by woody vegetation with a percent cover >60% and height exceeding 2 m. Almost all trees remain green all year. Canopy is never without green foliage</td><td>Permanent Wetlands: Lands with a permanent mixture of water and herbaceous or woody vegetation that cover extensive areas. The vegetation can be present in fresh water</td><td>Open Shrublands: Lands with woody vegetation less than 2 m tall and with shrub canopy cover between 10% and 60%. The shrub foliage can be either evergreen or deciduous</td></tr><tr><td>Soil texture (USDA)</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>Sand</td><td>–</td><td>–</td></tr><tr><td>Soil temperature depths [cm]</td><td>2, 8, 12, 16, 20, 30, 40, 50</td><td>1, 5, 10, 30, 50, 75</td><td>2, 4, 12, 25, 50</td><td>5, 15, 30</td><td>7, 15, 30</td><td>5, 30</td><td>5, 10, 36</td><td>4, 8</td></tr><tr><td>Database</td><td>GEM</td><td>GEM</td><td>FLUXNET</td><td>FLUXNET</td><td>FLUXNET</td><td>FLUXNET</td><td>FLUXNET</td><td>FLUXNET</td></tr><tr><td>Link/doi</td><td>https://data.g-e-m.dk</td><td>https://data.g-e-m.dk</td><td>https://doi.org/10.18140/FLX/1440158</td><td>https://doi.org/10.18140/FLX/1440159</td><td>https://doi.org/10.18140/FLX/1440228</td><td>https://doi.org/10.18140/FLX/1440160</td><td>https://doi.org/10.18140/FLX/1440181</td><td>https://doi.org/10.18140/FLX/1440182</td></tr><tr><td>Reference</td><td>López-Blanco et al. (2020)</td><td>López-Blanco et al. (2020)</td><td>Suni et al. (2003)</td><td>Lohila et al. (2004)</td><td>Aurela et al. (2015)</td><td>Thum et al. (2007)</td><td>Merbold et al. (2009)</td><td>Van der Molen et al. (2007)</td></tr></tbody></table>
selected as seen in Table 1, with the following criteria: the site being above 60°N of latitude, and the availability of soil temperatures. The data sets contain time series of meteorological data and, for the selected sites, soil temperatures at various depths.

For both databases, the meteorological data are used for the preparation of the forcing files, while the soil temperatures serve as validation data.

### 3.1.2. SoilGrids250m Version 2.0

SoilGrids250m 2.0 is a worldwide database of soil characteristics, from the soil standard properties (organic carbon, bulk density) to its chemical composition for six standard soil depths (between 0–5 cm, 5–15 cm, 15–30 cm, 30–60 cm, 60–100 cm and 100–200 cm) (Hengl et al., 2014, 2017; Poggio et al., 2021). It is developed by the International Soil Reference Information Center (ISRIC, http://www.isric.org/) and is available under an Open Data Base License. It has been regularly updated during the last years, with multiple evolutions from SoilGrids 1 km to SoilGrids250m. Thanks to machine learning and the use of more than 150,000 soil profiles for the training, this database can now provide global predictions at a resolution of 250 m.

When evaluated with cross-validation measures on the variables provided by the data set, 61% of the variation could be explained with that model. Among the multiple variables provided by SoilGrids 2.0, three of them are used in this study: the soil organic carbon concentration (soc), the bulk density of fine matter (bdod) and the coarse matter fraction (cfvo).

### 3.1.3. WFDE5 CRU GPCC

In order to prepare the forcing files for each site, the WFDE5 data set has also been used (Cucchi et al., 2020). This data set corresponds to the ERA5 reanalysis bias-adjusted using the WATCH Forcing Data methodology (Weedon et al., 2014), combined to climatological data from the Climate Research Unit (CRU) (Harris et al., 2020) to correct air temperature and downward radiation. Rainfall and snowfall rates are also readjusted by using the number of wet days from the CRU data and the precipitation totals of the Global Precipitation Climatology Center (GPCC) data (Schneider et al., 2018). It has a spatial resolution of 0.5° and an hourly temporal resolution that spans from 1979 to 2019. The data set can be downloaded from the Climate Data Store of the Copernicus Climate Change Service (C3S): https://doi.org/10.24381/cds.20d54e34. For each studied site, the meteorological variables are extracted for the nearest grid cell.

## 4. Methods

### 4.1. Site Level Soil Composition and Atmospheric Inputs

To perform the site simulations with ORCHIDEE, multiple inputs are needed. For this particular work, the specification of the soil properties is paramount, since the soil's bulk density, its SOC content and its coarse matter fraction are necessary. Moreover, the soil texture as well as the vegetation have to be prescribed. Finally, as the model is run offline, it is necessary to generate the atmospheric forcings for each site based on the meteorological data.

#### 4.1.1. Specification of Soil Properties

The soil composition used in the thermal module (SOC, OM, MM and CM) for each pixel is derived from the SoilGrids database, presented in Section 3.1.2. The variable $\rho_{SOC,i}$, needed for DisC, and the variables $\psi_{CM}$, $BD_{FM}$ and $f_{SOC}$ used in NewDisC, presented in Section 2.2, are directly provided by the database, but with different units and a different soil discretization. The details regarding the full processing of the SoilGrids data are provided in Appendix A1.

Multiple differences arise from the use of the database between DisC and NewDisC. For DisC, only the SOC density provided by the database is used, instead of the combination of SOC concentration, bulk density and coarse matter fraction needed by NewDisC. Also, the data in DisC comes from the previous version SoilGrids1km (Hengl et al., 2014), while the most up-to-date version SoilGrids250m 2.0 (Poggio et al., 2021) is used with NewDisC. Finally, the vertical interpolation used in both configurations differs. In soils, the vertical distribution of SOC, and therefore OM, generally decreases with depth (Von Fromm et al., 2024). Accordingly, in DisC, SOC
is distributed vertically via a decreasing exponential law, inserting most of the SOC near the soil's surface. However, this approach does not conserve the total SOC mass per layer and may lead to an inaccurate distribution of SOC. In NewDisC, a different assumption is made: the soil properties are assumed to be homogeneous in each SoilGrids depth layer, as described in Section A1, leading to a piece-wise linear interpolation. This approach allows SOC to decrease with depth across the successive layers without introducing additional parameters to represent the intra-layer variation, which could introduce more uncertainties since the observed rate of SOC decrease is neither linear nor consistent across soil types and climatic regions (Von Fromm et al., 2024). This implies a different allocation of OM along depth in NewDisC compared to DisC, NewDisC containing less OM near the surface than DisC, but more in depth, leading to distinct thermal properties.

### 4.1.2. Prescription of Soil Texture and Vegetation for Each Site

The soil texture and vegetation type (given as a PFT) applied for each site during the simulation have been determined based on the information available for each site (see Table 1). It should be noted that wetlands, fens and shrublands are not described by the ORCHIDEE PFTs and are therefore represented as boreal C3 grass in the present study. Table S4 in Supporting Information S1 summarizes the prescribed soil textures and PFTs for each site.

### 4.1.3. Atmospheric Forcing

To correctly run offline simulations with ORCHIDEE, whose temporal resolution is 30 min, it is necessary to dispose of continuous and long-term time series of atmospheric forcings, while local forcings usually lack these properties.

To benefit from both long-term and continuous time series, and local in situ measurements, we generated atmospheric forcing files from the in situ meteorological station gap-filled with the WFDE5 reanalysis product. The resulting time series covers the period from 1980 to 2018. The missing values were inferred by a set of linear regressions between the reanalyzes and the observations calculated at different time scales to account for the temporal variability of the biases. The regression parameters were defined at hourly and monthly scales to capture both the diurnal and seasonal biases, resulting in the calculation of $24 \times 12 = 288$ sets of constant (a, b) parameters so that the missing variable (Missing_variable) can be estimated from the reanalysis (Reanalysis_variable):

$$ \text{Missing\_variable} = a(\text{month}, \text{hour}) \times \text{Reanalysis\_variable} + b(\text{month}, \text{hour}) \quad (6) $$

Finally, each forcing variable (Gapfilled_variable) is calculated as:

$$ \text{Gapfilled\_variable} = \begin{cases} \text{Obs\_variable}, & \text{if observations available and with acceptable quality} \\ \text{Missing\_variable}, & \text{or if } \left| \frac{\text{Obs\_variable} - \text{Missing\_variable}}{\text{Missing\_variable}} \right| > 50\% \\ \text{Missing\_variable}, & \text{else.} \end{cases} \quad (7) $$

A particular treatment was applied to precipitation, since rainfall and snowfall are not separated in the GEM database. To differentiate rainfall from snowfall, we used the air temperature and the snow depth local measurement. When the temperature is below 0°C, or when the temperature is above 0°C and the snow depth increases, the precipitation is considered as snow. Otherwise, it is considered as rain.

Regarding the gap-filling procedure, contrary to the other atmospheric variables bias-corrected at the hourly scale, we chose to correct the precipitation at the monthly scale, and thus, we derived only 12 sets of (a, b) parameters and applied the same methodology as for the other meteorological variables.

The procedure is a bit simpler with the FLUXNET2015 data sets, since they are continuous and the data are already gap-filled following Vuichard and Papale (2015). The procedure for extending the duration of the time series of the forcing files is exactly the same as the one described for the GEM data sets, without the steps consisting of recalculating the variables of interest from those available (in particular, distinguishing snowfall from rainfall).
## 4.2. Simulation Protocol

In order to check the ability of the new version NewDisC to correctly assess the soil thermal behavior in the presence of OM, it has been compared to the two other possible configurations previously presented: Bulk and DisC. The simulations were done with ORCHIDEE v3.0.

The soil temperature profile (up to 90 m) was initialized using the mean local ERA reanalysis air temperatures from 1960 to 1990 (Gouttevin et al., 2012). This approach provides an initial thermal state representative of near-present steady-state conditions and reduces the spin-up time required for the thermal module. To equilibrate the water and energy budgets, a 100-year spin-up was then performed by looping the atmospheric forcing 10 times over the years 1980–1989. This period was selected because the WFDE5 data set begins in 1980 and does not provide earlier forcing data (Cucchi et al., 2020). The historical simulations were then run from 1990 to 2018 for each configuration, and only the years where observation data were available were used for the analyses.

For each site, the atmospheric forcing is based on the in situ measurements gap-filled following the method presented in Section 4.1.3. The forcings cover the years 1980–2018. In the particular case of Zackenberg (GL-Zaf), two further simulations were carried out with the DisC and NewDisC configurations, but with the forcing provided directly by the WFDE5 reanalysis, in order to evaluate the sensitivity of the model to the atmospheric forcing and especially to the snowfall uncertainties.

For both DisC and NewDisC simulations, the SOC properties used in the thermal module are fixed and derived from different SoilGrids maps containing the variables described in Table A1. The details regarding the preparation of the maps for the model are provided in Section A1. For DisC, a map providing the organic carbon density on each ORCHIDEE level is used, while for NewDisC, three maps are used: soc, bdod and cfvo. In both cases, the selected soil profiles correspond to those of the site's nearest pixel on the maps.

All the parameters and the different configurations of the simulations are synthetized in Table S5 in Supporting Information S1.

# 5. Results

## 5.1. Case Study: Zackenberg (GL-Zaf)

### 5.1.1. Atmospheric Conditions and Local Soil Temperatures

Zackenberg is a long-term monitored site for meteorological conditions, soil properties, biome characterization and greenhouse gas emissions by the GEM program that will serve to validate simulation performances. The site is dominated by fens and is considered as a peatland in formation with an accumulation of less than 30 cm peat carbon that lies on top of early-Holocene delta sediments (Christensen et al., 2023; Christiansen, 2004).

Figure 3 presents the daily mean seasonal cycle of the air and soil temperatures (upper panel) as well as monthly rainfall and snowfall (lower panel), all averaged over the period 12 August 2009–29 May 2017. Air temperature and precipitation have been gap-filled using the procedure described in Section 4.1.3 with the partitioning of rain and snow inferred from the in situ air temperature and snow depth data. In the lower panel of Figure 3, we observe that early snow falls in September and lasts until May. Soil temperature amplitudes decrease with depth, due to the insulating effect of OM and snow (Bonan & Shugart, 1989; Lawrence & Slater, 2008). This effect is especially visible in winter between September and May, mostly because of the snow insulation (Sokratov & Barry, 2002). Moreover, there is a time lag between the air temperature and soil temperatures that is positive in April (soil warming is delayed) and negative in early September (soil freezing is delayed) where the soil temperatures remain close to 0°C while the air temperature increases/drops. This reflects the effect of the latent heat released or accumulated when the soil water freezes after the summer or thaws in spring (Boike et al., 1998).

Additionally, the site forcing has been compared to the WFDE5 CRU GPCC forcing. The two atmospheric forcings (Figure S1 in Supporting Information S1) differ essentially in terms of snowfall. The air temperatures differ, with a RMSE between the two annual mean cycles of 2.4°C, and larger diurnal and seasonal temperature amplitude for the local forcing than for the WFDE5 one. The precipitation amounts are quite different with WFDE5 showing larger mean annual precipitation of 345 mm compared to 229 mm with the site's forcing, and especially larger snowfall during spring, about three times higher than from site observations. Overall, the snow represents 90% of the annual precipitation with WFDE5, while it is only 60% with the site's data. For site
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Annual cycle of the air (dotted gray line) and soil daily mean temperatures (solid lines) measured at 8 depths, 0.02, 0.08, 0.12, 0.16, 0.2, 0.3, 0.4 and 0.5 m and of the monthly mean rainfall and snowfall (in mm/month) in Zackenberg, Greenland. The data is averaged over the period 2009–2017 that corresponds to the temporal range of the soil temperature measurements.</figcaption></figure>

measurements, underestimation of precipitation and especially the amount of snowfall is a common and well-known issue with rain gauges (Goodison et al., 1998; Grossi et al., 2017; Liljedahl et al., 2017). The underestimation can be attributed to multiple factors but is mostly due to the wind that affects the trajectory of the snowflakes, leading in some cases to a 100% underestimation (Pan et al., 2016). Such discrepancies in the forcing data can lead to different simulation results, by impacting the amount of snow on the soil, its thermal properties, and as a result, the soil surface and in-depth temperature. To evaluate the impact of snow amount uncertainties on the results, simulations have been performed using both the gap-filled and the WFDE5 forcing.

### 5.1.2. ORCHIDEE Simulations

Figure 4 presents the inter-annual mean (2009–2017) of the seasonal cycle of observed and simulated daily soil temperatures at various depths. For each depth, the soil temperatures simulated on the 18 soil layers of ORCHIDEE are linearly interpolated to the observed depths. The daily mean temperatures are then evaluated on the time period where the observations are available. The temperatures estimated with the three model versions Bulk, DisC and NewDisC with the site forcing are presented in solid lines, and the results obtained with the WFDE5 forcing appear in blue and green dotted lines. With the site forcing, independently of the used model version, the simulated temperatures respectively increase and decrease too early at the end of winter and summer. They also tend to be too low in winter, and too high in summer with errors of up to 5°C although the discrepancies vary depending on the simulations.

These errors are slightly reduced in the two versions accounting for SOC, even though they both still present higher temperatures in summer and lower temperatures in winter compared to the observations, at all levels up to 0.5 m. During the intermediate seasons, all versions tend to heat up in spring and to cool down in fall too quickly, highlighting possible underestimation of the latent heat exchanges during melting/freezing phases. When using the WFDE5 forcing, the transition periods when snow appears/melts are better simulated by the model, and NewDisC fits well the observations during the snow-free period, while presenting a slightly larger warm bias during winter.

To analyze the results quantitatively, a global RMSE comparing the observed and simulated soil temperatures is calculated, as well as the seasonal RMSEs (calculated for DJF, MAM, JJA and SON) for each simulation and each depth, as explained in Appendix B.

Figure 5 shows the total RMSE values obtained for each depth and simulation as a heatmap on the left panel, the best results being displayed in green and the worst in red, while the seasonal RMSEs are presented on the right panel as bar plots. As shown in Figure 5 (left panel), NewDisC shows better results at each depth except for the top
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Observed and simulated interannual averaged soil temperatures at daily time step in Zackenberg with the two model versions (DisC and NewDisC) forced with either the WFDE5 forcing (dotted lines) or the site's forcing (solid lines). The black solid line is the mean of the observed temperatures, the gray shading is the range between the maximum and minimum of the observed temperatures for each day of year. All data are processed and compared over the same time period, that is, 12 August 2009–29 May 2017.</figcaption></figure>

layer (0.02 m) with lower RMSE values by 0.14–1.39°C compared to the two other simulations with the site forcing. On the top layer, the RMSE value is 0.39°C higher for the NewDisC version compared to the Bulk version, but remains lower than that of the DisC version. At the seasonal scale (Figure 5, right panels), NewDisC performs better than the other parameterizations for MAM and SON. However, the DisC simulation tends to have lower RMSE values during summer (JJA), while having a cold bias during the rest of the year. Bulk simulation shows slightly smaller RMSEs in winter (DJF) compared to the SOC simulations, but errors are larger in summer.

The RMSE values for both site simulations that include SOC, that is DisC and NewDisC, are larger in winter (DJF) compared to those of the Bulk simulation with the site forcing. This can be explained by an underestimation of the snow cover during winter, directly linked to the atmospheric forcing used. The results obtained with the WFDE5 forcing show that the simulated soil temperatures are significantly improved compared to those obtained with the site forcing, particularly for the NewDisC model version, where the RMSE values are the lowest at most depths, except for the top soil layers (0.02 and 0.08 m). As expected, the cold winter bias that was observed with the site forcing is not observed anymore but becomes a warm bias with an overestimation of winter temperatures by 4–5°C with NewDisC. However, the MAM and SON temperatures are in better agreement with the observations and the temporal dynamics is also better phased compared to the original simulations. Furthermore, the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Annual (left panel) and seasonal (right panels) evaluation of simulated soil temperature using RMSE. Results are presented for the three simulations performed with the site forcing: Bulk_site, DisC_site and NewDisC_site and compared to the results obtained with the WFDE5 forcing: DisC_WFDE5 and NewDisC_WFDE5 at Zackenberg (GL) Trimester seasons are indicated by DJF, December, January, February; MAM, March, April, May; JJA, June, July, August; SON, September, October, November.</figcaption></figure>

summer temperatures using NewDisC fit the observed ones very well. It means that, with a more accurate estimation of the snowfall amount, the new parameterization enables a better prediction of the soil temperatures.

The use of the WFDE5 forcing improves the DisC results as well, but with different effects. Indeed, the winter temperatures are much more accurate than for the other simulations since the winter RMSE values drop from 5°C to 3°C, but the insulative effect during the other seasons is overestimated, in particular below 0.12 m, leading to larger total RMSE values.

This demonstrates the sensitivity of the model to the atmospheric forcing and especially the snowfall, as well as the ability of the new discretization to accurately determine soil temperature annual variability.

## 5.2. Multi-Site Results

In order to have a more complete overview of the abilities of the new parameterization to correctly model the thermal behavior of the soils in high latitude environments, the model evaluation was extended to the other studied sites presented in Section 3.1.1 (Figure 2 and Table 1). To have an overview of the results, two complementary metrics are used: the KGE and the skill score, whose properties are detailed in Appendix B. The KGE provides an evaluation of the model behavior by accounting for correlation, bias and variability of the simulation results compared to the observations, while the skill score pinpoints the absolute prediction error difference between two model versions. In this study, these metrics are applied to the soil temperature time series for each of the nine studied sites.

Figure 6 presents the KGE values obtained on simulated versus observed soil temperatures, for each site (represented by the different letters) and each depth measurement. In most cases, NewDisC presents higher KGE values, showing that the new model performs better than the previous versions. It is the case for FI-Jok, RU-Cok and for the lower depths of GL-Zaf. However, some other sites are better represented by DisC, such as FI-Hyy.

A more complete interpretation of the results can be done by estimating for each site and depth the seasonal skill scores, defined in Appendix B. These results are shown in Figure 7. The first thing that can be noted is that the skill score is most of the time positive for NewDisC, which proves its good performances and the interest of better accounting for SOC in ORCHIDEE. Compared to the Bulk simulation, the RMSE values generally tend to decrease down to 50%, with better results in summer (JJA), when the soil is not covered by snow and when the temperatures are mostly influenced by the soil's thermal insulative properties. The DisC model version can
<figure><img src="image_9.png" /><figcaption>**Figure 6.** KGE of the soil temperatures for each site (letters) and simulation (orange for Bulk, blue for DisC and green for NewDisC), and for each depth where observations are available.</figcaption></figure>

sometimes provide better results. It is notably the case in JJA with skill scores that can go over 75%. However, this version has some disadvantages, as negative values are observed in the other seasons, with a decrease in RMSE values up to -80% for GL-Nuf, for example. This means that this version of the model can sometimes produce large errors, which is less the case with NewDisC.

The analysis of the thermal conductivities simulated by the three model versions reveal that they fully explained the discrepancies shown. The soils are generally more conductive with NewDisC compared to DisC, especially in the top soil layers, because of the different vertical distribution of the OM which is generally lower at the surface in NewDisC compared to DisC. Indeed, the vertical interpolation to distribute OM in depth is different between the two versions, as seen in Section A1. An example of the estimated OM volumetric concentrations from both versions is provided in Figure S2 in Supporting Information S1 for the site of Zackenberg. In NewDisC, OM is assumed to be homogeneous within the SoilGrids layers, while DisC is using a decreasing exponential repartition, meaning that OM is more concentrated close to the surface in DisC compared to NewDisC, which is quite visible in Figure S2 of Supporting Information S1. The use of the new map from SoilGrids, combined to the choice to use the combination of bulk density, soil organic carbon concentration and coarse matter fraction also leads to different values of SOC densities before the vertical interpolation. It may be linked to the way the algorithm of SoilGrids works, that is a calculate first then interpolate later approach (Poggio et al., 2021).

In consequence, we observe larger soil temperatures in summer and deeper melting fronts with NewDisC compared to DisC, which translates into an increase of the simulated active layer thickness (ALT) in permafrost regions. As expected, both discretized versions calculate lower ALTs than the Bulk version (less insulated), which tended to largely overestimate the active layer thickness. For example, in Zackenberg, ALT is much larger in NewDisC compared to DisC, moving from 0.4 to 0.7 m, which is in better agreement with the observations (which are around 0.6 m according to Strand et al. (2021)), while the simulated ALT from Bulk is close to 1.4 m.

In general, we can see that the improvement of simulated temperatures is larger in depth, with increasing values of skill score with depth. This result was expected, because of the cumulative effect of the soil OM on the attenuation of the thermal wave propagating from the surface. Results in Greenland (GL-Zaf and GL-Nuf) are not as good as for the other sites. For GL-Zaf, results appearing in Figures 6 and 7 correspond to the site's forcing, and it has previously been shown that changing the forcing could lead to much better results. A similar behavior can be expected for GL-Nuf. Moreover, both sites are sedge dominated fens with an abundant moss layer (López-Blanco
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Seasonal skill scores of the soil temperatures for DisC (blue) and NewDisC (green) compared to Bulk for each site, and at each depth where observations are available.</figcaption></figure>

et al., 2020), which is not explicitly represented by this version of the model and can affect the results (Gaillard et al., 2025).

# 6. Discussion and Limitations of the Study

In this work, we aimed to improve the soil thermal parameterization of the ORCHIDEE model by considering the soil composition and the repartition of its different components, that is fine mineral matter, coarse matter, OM, liquid and solid water and air. Our new model (NewDisC) has been compared to two previous versions, Bulk and DisC, at site level. Our results helped providing answers to the following scientific questions.

## 6.1. Why Accounting for OM Is Critical in LSMs, and How Should OM Be Parameterized to Infer Soil Thermal Properties?

The objective of our study was to enhance soil thermal behavior predictions in northern high-latitude regions, by incorporating the influence of OM in ORCHIDEE. Soils in these regions, particularly in the Arctic, contain high levels of OM and therefore SOC, leading to soils that have significantly different thermal and hydrological properties compared to mineral soils. Organic carbon-rich soils affect subsurface thermal profiles, a point demonstrated in previous studies (Boelter, 1969; Letts et al., 2000), requiring models to consider this aspect. Our results indicate that OM inclusion leads to more accurate soil temperature predictions, facilitated by the use of high-resolution SOC and mineral content data from SoilGrids within ORCHIDEE.
We show that the accuracy of model predictions depends strongly on how OM is prescribed in the model. The new approach developed in this work, to determine soil thermal properties by accounting for the soil's bulk density, its carbon content and also its coarse matter volume, was shown to be essential in enhancing model performance. Our OM parameterization, compared to previous models that only considered SOC or OM content, generally yields more realistic soil temperatures on the studied sites.

## 6.2. What Are the Effects of OM on Soil Hydrology, and What Feedback Mechanisms Are Involved?

While OM parameterization has been integrated into the thermal module of ORCHIDEE, it has not been included into the hydrological module, which would necessitate larger developments. In this study, the calculation of water transfers uses parameters such as the porosity which are calculated with the only mineral soil properties, and thus, the hydrological effects observed are primarily due to thermal changes associated with OM presence. Figure S3 in Supporting Information S1 presents the annual cycles of surface runoff, deep drainage, total soil moisture and upper soil moisture for the three configurations in Zackenberg. The hydrological cycles tend to follow the same trends for the three configurations, but with different orders of magnitude. Soil moisture tends to be larger in the discretized versions, as well as the deep drainage. Introducing OM into the hydrological module could further impact soil moisture dynamics, as SOC-rich soils generally exhibit higher water-holding capacity, which also means that they would also require more water to reach saturation (Rawls et al., 2003).

A wetter simulated soil in ORCHIDEE should result in less conductive soils, while reducing the summer surface temperatures. In addition, wetter soils should increase evapotranspiration, with an additional negative feedback on the summer surface temperatures.

Future model developments should therefore explore the inclusion of SOC in the calculation of the hydrological processes. Given that SOC-rich soils present a larger porosity, this adjustment could also alter the soil's hydraulic properties, affecting water retention, runoff, and drainage, as well as freeze-thaw cycles. This integration could lead to a more comprehensive representation of soil behavior in organic-rich regions, accounting for interactions between thermal and hydrological modules. Other implementations have already included the effects of SOC on some hydrological parameters, for instance in ORCHIDEE MICT, where the porosity, the wilting point and the field capacity depend on the SOC content. It has been shown that higher carbon contents lead to larger plant available water capacities, closer to the observations, in addition to the improved soil temperatures (Guimberteau et al., 2018; Hudson, 1994). In the study of Lawrence and Slater (2008), the modification of the hydraulic parameters to account for SOC have shown an increase in column soil water content and saturation levels. However, at the surface, the saturation levels were lower because of a more efficient water transport in the soil, reducing soil evaporation and latent heat flux. The account of SOC in the calculation of the hydrological processes and their feedbacks with the soil thermics is the next step of our work.

## 6.3. How Do Initial Conditions and Atmospheric Forcing Influence Model Outcomes?

In this study, ORCHIDEE was run offline, constrained by atmospheric forcing data. It has been shown that precipitation tends to be the most uncertain variable among the atmospheric forcing data (Gelati et al., 2018). Uncertainties in the precipitation prescription can lead to underestimations of the snow cover fraction, snow water equivalent and snow depth, with a more pronounced impact in high-latitude regions, where the sensitivity to these atmospheric inputs peaks during winter months (Wang et al., 2020). This can in turn greatly affect the soil's temperatures, since snow has a strong impact on the ground thermal regime because of its low thermal conductivity and high albedo (Gouttevin et al., 2012). Our results align with these findings, demonstrating that atmospheric forcing data—particularly precipitation and air temperature—strongly influence soil temperature by altering snow cover duration and depth.

Initial conditions, such as soil temperature before spin-up, soil texture, and vegetation classification via PFT maps, also introduce uncertainties into the model simulations. Initial soil temperatures based on past large-scale climatological values, which may slightly differ from local conditions, are one source of model errors. In addition, the spin-up procedure, based on the repetition of the 1980–1989 forcing, may lead to slightly warmer soil conditions, as this period represents relatively warm recent climate conditions. Soil texture, a critical factor in defining thermal properties, must be accurately represented to avoid erroneous outputs. Likewise, SOC maps from SoilGrids, while improving in resolution and accuracy, still contain uncertainties (Poggio et al., 2021), particularly in remote locations where they are based on extrapolations from sparse measurements. Moreover, the
PFT classification used in this study does not account for northern peatlands, and sites such as Zackenberg or Kobbefjord, which are covered by fens and bogs are simply represented by boreal grasslands. Shrublands are also poorly simulated since shrubs are not explicitly represented and are simply considered as highly sparse forests.

New developments to represent such ecosystems in the model are in progress, aiming at solving these issues.

Furthermore, discrepancies between the simulation results and the observations can be linked to the disparity of the origin of the observation data on site. For example, all the observations and measurements for one site might not be taken exactly at the same place. It is the case in Zackenberg where the measurements are made in an area with a radius of more than 1 km (Christensen et al., 2021). Some differences can thus arise, depending on the vegetation heterogeneity of the site, its microclimate or its altitude.

## 6.4. What Are the Study's Limitations, and What Future Model Improvements Are Anticipated?

This study is based on ORCHIDEE-V3.0, which lacks some features included in the latest model version (V4.0), such as a more detailed snow scheme that discretizes the snowpack into 12 layers (Charbit et al., 2024) and has been shown to better simulate snow depth and temperature. Additionally, updated PFT maps (Harper et al., 2023) in ORCHIDEE V4.0 offer potentially improved vegetation representation of peatlands, shrubs and forests. Future work aims to integrate our SOC-driven thermal parameterization into the updated model version to evaluate its effectiveness with these added features.

Further model developments should focus on incorporating SOC into the hydrological processes to more accurately represent organic-rich soils. Including SOC effects in hydrology, along with enhanced snow modeling and vegetation mapping, could provide a more global understanding of soil-atmosphere interactions in boreal and Arctic regions.

Finally, only a few sites have been tested here. Further work will include a more general analysis of the impact of the inclusion of this new SOC parameterization in the model by using the SoilGrids data and analyzing the results at the global scale.

## 7. Conclusion

In this study, we introduced a new formulation for the soil organic and mineral content within the thermal module of ORCHIDEE, referred to as the NewDisC version. We evaluated its performance across multiple sites in the Arctic region and compared it to two other versions: the Bulk version, which does not account for OM in the thermal and hydrological modules, and the DisC version, which considers a very simple prescription of SOC in the thermal module. Our results demonstrate that both SOC-inclusive models (DisC and NewDisC) generally improved simulated soil temperatures compared to the Bulk version, with a more pronounced effect at greater soil depths. Moreover, the NewDisC version usually outperformed the DisC model, while providing a more accurate representation of the underlying physics. This improvement is attributed to the NewDisC version's ability to calculate pixel- and depth-wise porosities and thermal properties using high resolution global maps of soil properties such as SoilGrids250m 2.0.

Future work will involve extending the evaluation of this new version to the global scale and enhancing the hydrological module of ORCHIDEE to account for the presence of OM, which could further influence the soil thermal and hydrological dynamics and predictions of the evolution of boreal and arctic permafrosts.

# Appendix A: Physical Properties in the Current and Modified Versions of ORCHIDEE

## A1. SOC Input for the Thermal Module of ORCHIDEE

Table A1 lists the soil properties extracted from the SoilGrids database used as an input in ORCHIDEE, depending on the versions (DisC or NewDisC).

In DisC, a map containing in-depth discretized and pixel-wise SOC densities can be used in the thermal module. This map is derived from the "Soil Organic Carbon Stock" (ocs) product of SoilGrids1km (Hengl et al., 2014). It provides ocs on 6 layers: 0–5 cm, 5–15 cm, 15–30 cm, 30–60 cm, 60–100 cm and 100–200 cm. For each depth, this variable is converted into an organic carbon density by multiplying it by the depth of each layer. Then a pixel-
<table><caption>**Table A1**
Soil Properties Description (Hengl et al., 2017; Poggio et al., 2021)</caption><thead><tr><th>Soil property</th><th>Acronym</th><th>Units</th><th>Description</th><th>Version</th></tr></thead><tbody><tr><td>Organic carbon stock</td><td>ocs</td><td>ton/ha</td><td>Soil organic carbon stock per soil layer (from SoilGrids1km)</td><td>DisC</td></tr><tr><td>Bulk density</td><td>bdod</td><td>cg/cm³</td><td>Bulk density of the fine earth fraction oven dry (from SoilGrids250 m 2.0)</td><td>NewDisC</td></tr><tr><td>Organic carbon concentration</td><td>soc</td><td>dg/kg</td><td>Gravimetric content of organic carbon in the fine earth fraction of the soil (from SoilGrids250 m 2.0)</td><td>NewDisC</td></tr><tr><td>Coarse fragments</td><td>cfvo</td><td>cm³/dm³</td><td>Volumetric content of fragments larger than 2 mm in the whole soil (from SoilGrids250 m 2.0)</td><td>NewDisC</td></tr></tbody></table>

wise vertical interpolation is performed on the depth (z-axis) to fit a decreasing exponential law: $a \times \exp(bz) + c$ (with $0 < a < 10^7$, $-15 < b < 0$ and $0 < c < 10^6$) and a value is then calculated for each ORCHIDEE depth. The total amount of SOC on these new levels is then renormalized for mass conservation. Finally, it is assumed that there is no SOC below 2 m, and the organic carbon density is null below this depth.

In NewDisC, the soil properties are derived from the more recent SoilGrids250m 2.0 database 3.1.2. The maps are downloaded directly from the SoilGrids portal in the Virtual Raster Tile (VRT) format, then regridded to a lat/lon grid and saved as netCDF files with a coarser resolution of 0.1°. The variables $\psi_{CM}$, $BD_{FM}$ and $f_{SOC}$ presented in Section 2.2 are respectively obtained by converting the variables from Table A1 cfvo to m³/m³, bdod to kg/m³ and soc to kg/kg, provided on the 6 same layers as ocs (0–5 cm, 5–15 cm, 15–30 cm, 30–60 cm, 60–100 cm and 100–200 cm):

$$ \psi_{CM} = \frac{cfvo}{1000} \qquad (A1) $$

$$ BD_{FM} = bdod \times 10 \qquad (A2) $$

$$ f_{SOC} = \frac{soc}{10000} \qquad (A3) $$

Then, a vertical interpolation needs to be performed to fit the ORCHIDEE depths. When performing the interpolation, it is crucial to respect mass conservation. As the quantities of each element may vary greatly from one profile to another, we assume that each SoilGrids layer is completely homogeneous and that the properties remain independent of depth and constant on each SoilGrids layer.

For each layer, the masses of FM and SOC ($BD_{FM,tot}$ and $BD_{SOC}$) per unit of total volume (including the pores and the volume occupied by the coarse matter $\psi_{CM}$) are:

$$ BD_{FM,tot} = BD_{FM} \times (1 - \psi_{CM}) \qquad (A4) $$

$$ BD_{SOC} = f_{SOC} \times BD_{FM,tot} \qquad (A5) $$

The properties $BD_{FM,tot}$, $BD_{SOC}$ and $\psi_{CM}$ are then calculated for each ORCHIDEE soil layer by performing weighted arithmetic averages on the SoilGrids corresponding layers.

Similarly to DisC, we assume that the soil properties remain constant below 2 m and the values are the same globally. We assume that there is no SOC and no coarse matter below this depth so $f_{SOC}$ and $\psi_{CM}$ are set to 0, while $BD_{FM}$ is set to $1,620 \text{ kg.m}^{-3}$, which corresponds to the bulk density of a mineral soil with a porosity of 0.4 and a particle density of $2,700 \text{ kg.m}^{-3}$. Finally, three netCDF files containing the calculated properties on the ORCHIDEE levels are then generated and can be read by the model.

## A2. Porosity

In DisC, the soil porosity is determined by using the arithmetic weighted mean that uses the volume fraction of soil organic carbon $\psi_{SOC,i}$:
$$ \theta_{s,i} = (1 - \psi_{SOC,i})\theta_{s,MM} + \psi_{SOC,i}\theta_{s,OM} \quad (A6) $$

with the porosity of organic matter $\theta_{s,OM}$ set to 0.92, and the porosity of mineral matter $\theta_{s,MM}$ depending on the soil texture (see Table S2 in Supporting Information S1). This equation also applies to Bulk, with $\psi_{SOC,i} = 0$.

In NewDisC, this parameterization is revised. According to Figure 1 and assuming that there are no other materials in the soil aside from mineral matter, organic matter, water and air, the solid volume fractions of each component j on layer i $F_{j,i}$ and the total porosity of the soil unit $\theta_{s,i}$ add up to one. The exact soil porosity can therefore be determined using Equation A7 instead of Equation A6.

$$ \theta_{s,i} = 1 - (F_{MM,i} + F_{CM,i} + F_{OM,i}) \quad (A7) $$

<caption>**Table A2**</caption>

Modified Equations Between the Current Version of ORCHIDEE and NewDisC

<table><thead><tr><th colspan="2">Equations</th></tr></thead><tbody><tr><td>Heat capacity</td><td><strong>Bulk/DisC</strong><br/>C<sub>tot,i</sub> = C<sub>dry,i</sub> × (1 − θ<sub>s,MM</sub>) + F<sub>water,i</sub>C<sub>water</sub> + F<sub>ice,i</sub>C<sub>ice</sub> (A8)</td></tr><tr><td>Where</td><td>C<sub>dry,i</sub> = ψ<sub>SOC,i</sub>C<sub>dry,OM</sub> + (1 − ψ<sub>SOC,i</sub>) C<sub>dry,MM</sub> (A9)</td></tr><tr><td></td><td><strong>NewDisC</strong><br/>C<sub>tot,i</sub> = F<sub>OM,i</sub>C<sub>solid,OM</sub> + (F<sub>MM,i</sub> + F<sub>CM,i</sub>) C<sub>solid,MM</sub><br/>+F<sub>water,i</sub>C<sub>water</sub> + F<sub>ice,i</sub>C<sub>ice</sub> + F<sub>air,i</sub>C<sub>air</sub> (A10)</td></tr><tr><td>K<sub>solid,i</sub></td><td><strong>Bulk/DisC</strong><br/>K<sub>solid,i</sub> = (K<sup>qz</sup><sub>quartz</sub> × K<sup>1-qz</sup><sub>other minerals</sub>)<sup>1-ψ<sub>SOC,i</sub></sup> × K<sup>ψ<sub>SOC,i</sub></sup><sub>solid,OM</sub> (A11)</td></tr><tr><td></td><td><strong>NewDisC</strong><br/>K<sub>solid,i</sub> = (K<sup>F<sub>OM,i</sub></sup><sub>OM</sub> × K<sup>F<sub>quartz,i</sub></sup><sub>quartz</sub> × K<sup>(F<sub>MM,i</sub>+F<sub>CM,i</sub>-F<sub>quartz,i</sub>)</sup><sub>solid,other minerals</sub>)<sup>1/(F<sub>OM,i</sub>+F<sub>MM,i</sub>+F<sub>CM,i</sub>)</sup> (A12)</td></tr><tr><td>With</td><td>F<sub>quartz,i</sub> = qz × (F<sub>MM,i</sub> + F<sub>CM,i</sub>) (A13)</td></tr><tr><td>K<sub>dry,i</sub></td><td><strong>Bulk/DisC</strong><br/>K<sub>dry,i</sub> = K<sup>1-ψ<sub>SOC,i</sub></sup><sub>dry,MM</sub> × K<sup>ψ<sub>SOC,i</sub></sup><sub>dry,OM</sub> (A14)</td></tr><tr><td></td><td><strong>NewDisC</strong><br/>K<sub>dry,i</sub> = K<sup>ψ<sub>OM,i</sub></sup><sub>dry,OM</sub> K<sup>ψ<sub>CM,i</sub> + ψ<sub>MM,i</sub></sup><sub>dry,MM</sub> (A15)</td></tr></tbody></table>

Note. The Equations for Bulk and DisC are the Same, With $\psi_{SOC,i} = 0$ for Bulk. The Full Description of the Parameters Used in the Different Equations is Provided in Table S1 of Supporting Information S1.
### A3. Thermal Properties

This appendix presents the full set of equations used to compute soil thermal properties in the current version of ORCHIDEE (Bulk/DisC) and details the modifications introduced in NewDisC. Table A2 summarizes the modified equations, and more details are provided in the corresponding sections.

The evolution of the soil temperature (at each depth) is determined from a one-dimensional Fourier Equation A16, where T is the temperature, z is the depth, t is the time, $C_{tot,i}$ is the heat capacity of layer i, $K_{soil,i}$ is its thermal conductivity and LH is a term related to the latent heat exchanges related to the freezing/melting of soil water.

$$C_{tot,i} \frac{\partial T}{\partial t} = \frac{\partial}{\partial z} \left( K_{soil,i} \frac{\partial T}{\partial z} \right) + LH \quad (A16)$$

The latent heat term can be calculated thanks to Equation A17. More details can be found in Gouttevin et al. (2012).

$$LH = \rho_{ice} L \frac{\Delta\theta_{ice,i}}{\Delta T} \frac{\partial T}{\partial t} \quad (A17)$$

where $\theta_{ice,i}$ is the volumetric ice content of layer i, L the latent heat of fusion of ice and $\rho_{ice}$ the ice density.

### A4. Heat Capacity $C_{tot,i}$

In Bulk and DisC, the heat capacity is calculated following the works of Lawrence and Slater (2008) with Equation A8 and is a weighted sum of the dry volumetric heat capacities of the solid soil constituents $C_{dry}$, and of liquid water $C_{water}$ and ice $C_{ice}$. The air is not taken into account since its volumetric heat capacity is supposed to be accounted for in $C_{dry}$. The applied weights are the volumetric liquid and solid water contents $F_{water,i}$ and $F_{ice,i}$ for the water, and $1 - \theta_{s,MM}$ for the solid soil part, to only account for the volume occupied by the solid matter and not the pores. However, this is erroneous to multiply $C_{dry,i}$ by $1 - \theta_{s,MM}$ since $C_{dry,i}$ is a dry heat capacity and accounts for the porous volume. The use of this equation leads to an underestimation of the soil heat capacity.

For DisC, the capacity of the solid soil constituents $C_{dry,i}$ is also a weighted sum of the OM and MM heat capacities, while using $\psi_{SOC}$ as a weight, as seen in Equation A9. A first issue in this equation lies with the use of the SOC volume fraction $\psi_{SOC}$, which has been discussed previously, but also the use of the wrong parameter value for $C_{dry,OM}$, as the value used in the model is $2.5 \times 10^6 \text{ J K}^{-1} \text{ m}^{-3}$, which corresponds to the solid volumetric heat capacity of organic matter (Farouki, 1981).

In NewDisC, the calculation of the heat capacity is refined as an arithmetic weighted average of the solid volumetric heat capacities $C_{solid,i}$ of all the soil constituents Equation A10, therefore including the air's heat capacity to be more rigorous, although its heat capacity is negligible compared to the other components. Moreover, the contribution of the soil solid constituents is adjusted by using the solid volume fractions. The thermal properties of MM depend on the pixel's soil texture USDA class, and are provided in Table S2 of Supporting Information S1, while the soil thermal properties of the other constituents can be found in Table S3 of Supporting Information S1.

### A5. Thermal Conductivity $K_{soil,i}$

The thermal conductivity of layer i of a soil is calculated with the method proposed by Johansen (1975) via Equation A18.

$$K_{soil,i} = (K_{sat,i} - K_{dry,i}) Ke + K_{dry,i} \quad (A18)$$

To obtain $K_{soil}$, one thus needs to determine $K_{sat}$ which is the saturated conductivity and $K_{dry}$ is the dry conductivity of the soil (no water in the pores). Ke is the Kersten number and depends on the size of the mineral matter (coarse or fine) and its saturation level (see Section A8).
$K_{dry}$ and $K_{sat}$ directly depend on the soil composition and that is where the soil volume fractions intervene. There are several methods to calculate these two variables when the soil is not made of pure mineral matter and the most common are the use of arithmetic and geometric averages (Farouki, 1981; McGaw, 1969), whether it is to calculate $K_{dry}$ or $K_{sat}$. By default, it is the arithmetic method that is used in ORCHIDEE, but the geometric method, also used in other models such as JULES (Chadburn et al., 2015) or ISBA (Decharme et al., 2016) was also implemented and tested in some versions (Zhu et al., 2019), and is assumed to better represent the conductivity.

Indeed, according to Farouki (1981), the conductivity of the soil depends on the arrangement of the matter. In the parallel-flow model where the soil is divided into multiple columns, each being homogeneous and only containing one soil constituent, the conductivity $K$ is maximized and can be calculated directly by doing an arithmetic average using volume fractions $f_j$ of each component $j$ and their conductivity: $K = \sum_{j=1}^{n} f_j K_j$. In the series-flow model (each constituent occupies a different soil layer), the conductivity is minimized, and calculated as follows: $\frac{1}{K} = \sum_{j=1}^{n} \frac{f_j}{K_j}$. In practice, none of these cases are observed in the soils, where the components are more mixed. A good compromise between these two extreme cases is the geometric average, as suggested by Johansen (1975). The expression is then: $K = \prod_{j=1}^{n} K_j^f$. This formulation is chosen to calculate the dry and solid conductivities.

## A6. Dry Conductivity $K_{dry,i}$

In DisC, with the geometric average, the dry conductivity $K_{dry,i}$ of the soil is determined using Equation A14 (Johansen, 1975; Peters-Lidard et al., 1998). The mineral dry conductivity obtained with Equation A19 only depends on the mineral porosity $\theta_{s,MM}$ that is defined by the soil texture. Its values appear in Table S2 of Supporting Information S1. The dry conductivity of organic matter $K_{dry,OM}$ is set to 0.05 W K$^{-1}$ m$^{-1}$ (Lawrence & Slater, 2008) in the model.

$$K_{dry,MM} = \frac{0.135 \times (1 - \theta_{s,MM}) \times 2700 + 64.7}{2700(1 - 0.947 \times (1 - \theta_{s,MM}))} \quad (A19)$$

In NewDisC, very small adjustments are made. The dry conductivity takes into account the pores; it is therefore necessary to use the volume fractions $\psi$ as weights to calculate the geometric average, as was done in Equation A14. However, Equation A15 slightly differs from Equation A14 by considering the fractions of OM $\psi_{OM}$ instead of that of SOC $\psi_{SOC}$. $K_{dry,MM}$ stays the same and is calculated with Equation A19 as we assume that its value applies to both coarse matter and fine mineral matter.

## A7. Saturated Conductivity $K_{sat,i}$

The saturated conductivity $K_{sat,i}$, calculated with Equation A20, is a weighted geometric average of the conductivities of ice $K_{ice}$ and water $K_{water}$ on one hand, and of the solid soil conductivity $K_{solid,i}$ on the other hand (Lawrence & Slater, 2008). The pores are filled with either water, ice or both. The conductivity values used in the model are provided in Table S3 of Supporting Information S1.

$$K_{sat,i} = K_{solid,i}^{1-\theta_{s,i}} K_{water}^{F_{water}} K_{ice}^{\theta_{s,i}-F_{water}} \quad (A20)$$

In DisC, $K_{solid,i}$ is determined using Equation A11, with $qz$ the quartz content in the mineral matter, which depends on the soil texture, according to the USDA classification. Indeed, quartz (which is found in the fine and coarse mineral matter) has a conductivity 3 to 4 times as important as that of other minerals (Johansen, 1975). Its volume is therefore considered separately.

However, Equation A11 needs to be corrected. Indeed, the pores are not included in the solid soil and to calculate its conductivity, the solid volume fractions $F$ should be used instead of the soil volume fractions $\psi$ and normalized by the sum of all solid fractions (excluding the volume occupied by the pores), as suggested by Balland and Arp (2005), resulting in Equation A12 for the NewDisC version. The solid volume fraction of quartz is also recalculated from $qz$ to account for its high conductivity.
## A8. Kersten Number

Finally, the Kersten number that appears in Equation A18 depends on whether the water in the soil interstices is liquid or solid and on the degree of saturation $S_r$, which corresponds to the fraction of the porosity occupied by water and ice ($S_r = \frac{F_{water} + F_{ice}}{\theta_s}$). For unfrozen soils (Johansen, 1975; Peters-Lidard et al., 1998):

$$ Ke = \begin{cases} \log_{10}(S_r) + 1, & \text{if } S_r > 0.1 \text{ and fine matter} \\ 0.7\log_{10}(S_r) + 1, & \text{if } S_r > 0.05 \text{ and coarse matter} \\ 0, & \text{otherwise} \end{cases} \qquad (A21) $$

while for frozen soils, $Ke = S_r$. It is important to notice here that an incorrect estimation of the soil porosity $\theta_s$ may lead to wrong estimations of the Kersten number and affect the calculation of the soil conductivity.

## Appendix B: Metrics for Evaluation

### B1. Root-Mean-Square-Error (RMSE) Metrics

To evaluate the performances of each simulation, the RMSEs are estimated between the simulated soil temperatures $T_{simu,i}$ and the corresponding observed temperatures $T_{obs,i}$ for each day $i$ where observations are available.

$$ RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (T_{simu,i} - T_{obs,i})^2} \qquad (B1) $$

The seasonal RMSEs Equation B1 are defined using the days belonging to the 3 month period seasons. The seasons are: December to February (DJF), March to May (MAM), June to August (JJA) and September to November (SON).

### B2. Skill Score

The skill score Equation B2 is a metric that enables a direct comparison between two simulations. For instance, the NewDisC and DisC simulations are both compared to the Bulk simulation.

Values above 0 mean that the RMSE of the analyzed simulation is reduced compared to that of the Bulk simulation, and the comparative performance of the simulation increases when the skill score increases toward a maximum value of 1.

$$ \text{Skill score} = 1 - \frac{RMSE_{simu}}{RMSE_{Bulk}} \qquad (B2) $$

### B3. Kling-Gupta Efficiency (KGE)

The last metric to assess the performances of the different parameterizations is the KGE metric (Chaney et al., 2016; Gupta et al., 2009), described by Equation B3. It combines three metrics that compare the simulation results (simu) and the observations (obs): linear correlation ($r_{Pearson}$), variability ($\sigma_{simu}/\sigma_{obs}$) and mean bias ($\mu_{simu}/\mu_{obs}$), $\sigma_{simu}$ and $\sigma_{obs}$ being the standard deviations of the simulated and the observed temperatures, and $\mu_{simu}$ and $\mu_{obs}$ their mean values.

$$ KGE = 1 - \sqrt{(r_{Pearson} - 1)^2 + \left(\frac{\sigma_{simu}}{\sigma_{obs}} - 1\right)^2 + \left(\frac{\mu_{simu}}{\mu_{obs}} - 1\right)^2} \qquad (B3) $$

The ability of a model to fit the observations increases with the KGE, on the interval $]-\infty, 1]$.
## Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

## Data Availability Statement

For the Greenland sites, meteorological and soil data were extracted from the two sub-programs of the Greenland Ecosystem Monitoring (GEM) Programme (ClimateBasis and GeoBasis) (Christensen et al., 2017). The ClimateBasis data were provided by Asiaq—Greenland Survey, Nuuk, Greenland and are available via the GEM Data Portal (Greenland Ecosystem Monitoring, 2025a). The GeoBasis data were provided by the Department of Ecoscience, Aarhus University, Denmark, in collaboration with the Department of Geosciences and Natural Resource Management, Copenhagen University, and Asiaq—Greenland Survey, Nuuk, Greenland (Greenland Ecosystem Monitoring, 2025b). Meteorological and soil temperature data from FLUXNET2015 (Pastorello et al., 2020) are available on the FLUXNET data archive (Mammarella et al., 1996–2014; Lohila et al., 2000–2003; Aurela et al., 2007–2009, 2001–2014; Merbold et al., 2002–2005; Dolman et al., 2003–2014; FLUXNET2015-US-Prr, 2010–2014). The source code of ORCHIDEE and the guidelines to run the model are available online via the IPSL data catalogue (Cuynet, 2025).

## Acknowledgments

We would like to thank the ORCHIDEE Project Team for developing and maintaining the ORCHIDEE code and for providing the ORCHIDEE versions used in this study. This study benefited from the IPSL Data and Computing Center ESPRI (https://mesocentre.ipsl.fr) which is supported by CNRS, Sorbonne Université, CNES and Ecole Polytechnique, as well as national and international grants. We also thank all the data providers that have made this work possible, including the GEM and FLUXNET Site PIs and data management team, developers of the SoilGrids and WFDE5 CRU GPCC products, all of which are cited in papers referenced in this study. Most of the studies and results presented here were supported by the GREENFEEDBACK (Greenhouse gas fluxes and Earth system feedbacks) project funded by the European Union's HORIZON Research and Innovation programme under grant agreement number 101056921. This work was also partially funded by the ESA Climate Change Initiative, Contract 4000125156/18/I-NB. Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or CINEA. Neither the European Union nor the granting authority can be held responsible for them. Finally, the authors acknowledge the three anonymous reviewers who provided insightful comments that helped to improve the manuscript.

## References

* Aurela, M., Lohila, A., Tuovinen, J.-P., Hatakka, J., Penttilä, T., & Laurila, T. (2015). Carbon dioxide and energy flux measurements in four northern-boreal ecosystems at Pallas. *Boreal Environment Research*, 20, 455–473.

* Aurela, M., Lohila, A., Tuovinen, J.-P., Hatakka, J., Rainne, J., Mäkelä, T., & Lauria, T. (2007–2009). FLUXNET2015 FI-Lom polpojankka [Dataset]. https://doi.org/10.18140/FLX/1440228

* Aurela, M., Tuovinen, J.-P., Hatakka, J., Lohila, A., Mäkelä, T., Rainne, J., & Lauria, T. (2001–2014). FLUXNET2015 FI-Sod Sodankyla [Dataset]. https://doi.org/10.18140/FLX/1440160

* Balland, V., & Arp, P. A. (2005). Modeling soil thermal conductivities over a wide range of conditions. *Journal of Environmental Engineering and Science*, 4(6), 549–558. https://doi.org/10.1139/s05-007

* Beringer, J., Lynch, A. H., Chapin, F. S., Mack, M., & Bonan, G. B. (2001). The representation of arctic soils in the land surface model: The importance of mosses. *Journal of Climate*, 14(15), 3324–3335. https://doi.org/10.1175/1520-0442(2001)014(3324:TROASI)2.0:CO;2

* Boelter, D. H. (1969). Physical properties of peats as related to degree of decomposition. *Soil Science Society of America Journal*, 33(4), 606–609. https://doi.org/10.2136/sssaj1969.03615995003300040033x

* Boike, J., Roth, K., & Overduin, P. P. (1998). Thermal and hydrologic dynamics of the active layer at a continuous permafrost site (Taymyr peninsula, Siberia). *Water Resources Research*, 34(3), 355–363. https://doi.org/10.1029/97WR03498

* Bonan, G. B., & Shugart, H. H. (1989). Environmental factors and ecological processes in boreal forests. *Annual Review of Ecology and Systematics*, 20(1), 1–28. https://doi.org/10.1146/annurev.es.20.110189.000245

* Box, E. O. (1996). Plant functional types and climate at the global scale. *Journal of Vegetation Science*, 7(3), 309–320. https://doi.org/10.2307/3236274

* Campbell, G. S. (1985). *Soil physics with BASIC: Transport models for soil-plant systems*. Elsevier.

* Chadburn, S. E., Burke, E., Essery, R., Boike, J., Langer, M., Heikenfeld, M., et al. (2015). An improved representation of physical permafrost dynamics in the JULES land-surface model. *Geoscientific Model Development*, 8(5), 1493–1508. https://doi.org/10.5194/gmd-8-1493-2015

* Chadburn, S. E., Burke, E. J., Gallego-Sala, A. V., Smith, N. D., Bret-Harte, M. S., Charman, D. J., et al. (2022). A new approach to simulate peat accumulation, degradation and stability in a global land surface scheme (JULES vn5.8_accumulate_soil) for northern and temperate peatlands. *Geoscientific Model Development*, 15(4), 1633–1657. https://doi.org/10.5194/gmd-15-1633-2022

* Chaney, N. W., Herman, J. D., Ek, M. B., & Wood, E. F. (2016). Deriving global parameter estimates for the Noah land surface model using FLUXNET and machine learning. *Journal of Geophysical Research: Atmospheres*, 121(22), 13–218. https://doi.org/10.1002/2016JD024821

* Charbit, S., Dumas, C., Maignan, F., Ottlé, C., Raoult, N., Fettweis, X., & Conesa, P. (2024). Modelling snowpack on ice surfaces with the ORCHIDEE land surface model: Application to the Greenland ice sheet. *EGUsphere*, 2024(11), 1–38. https://doi.org/10.5194/tc-18-5067-2024

* Chen, Y., Yang, K., Tang, W., Qin, J., & Zhao, L. (2012). Parameterizing soil organic carbon's impacts on soil porosity and thermal parameters for Eastern Tibet grasslands. *Science China Earth Sciences*, 55(6), 1001–1011. https://doi.org/10.1007/s11430-012-4433-0

* Cheruy, F., Ducharne, A., Hourdin, F., Musat, I., Vignon, E., Gastineau, G., et al. (2020). Improved near-surface continental climate in IPSL-CM6A-LR by combined evolutions of atmospheric and land surface physics. *Journal of Advances in Modeling Earth Systems*, 12(10), e2019MS002005. https://doi.org/10.1029/2019MS002005

* Christensen, T. R., Lund, M., Skov, K., Abermann, J., López-Blanco, E., Scheller, J., et al. (2021). Multiple ecosystem effects of extreme weather events in the arctic. *Ecosystems*, 24(1), 122–136. https://doi.org/10.1007/s10021-020-00507-6

* Christensen, T. R., Scheller, J., Scheel, M., Rudd, D. A., Jackowicz-Korczynski, M., Mastepanov, M., & López-Blanco, E. (2023). Rapid shift in greenhouse forcing of emerging arctic peatlands. *Scientific Reports*, 13(1), 2828. https://doi.org/10.1038/s41598-023-29859-4

* Christensen, T. R., Topp-Jørgensen, E., Sejr, M. K., & Schmidt, N. M. (2017). Foreword: Synthesis of the Greenland ecosystem monitoring program. *Ambio*, 46(S1), 1–2. https://doi.org/10.1007/s13280-016-0860-z

* Christiansen, H. H. (2004). Meteorological control on interannual spatial and temporal variations in snow cover and ground thawing in two northeast Greenlandic Circumpolar-Active-Layer-Monitoring (CALM) sites. *Permafrost and Periglacial Processes*, 15(2), 155–169. https://doi.org/10.1002/ppg.489

* Cucchi, M., Weedon, G. P., Amici, A., Bellouin, N., Lange, S., Müller Schmied, H., et al. (2020). WFDE5: Bias-adjusted ERA5 reanalysis data for impact studies. *Earth System Science Data*, 12(3), 2097–2120. https://doi.org/10.5194/essd-12-2097-2020

* Cuynet, A. (2025). ORCHIDEE_3_NewDisC_r8617 [Software]. IPSL. https://doi.org/10.14768/3b1dc41a-0c9e-4292-8ea2-d2dac6a10c29
Decharme, B., Brun, E., Boone, A., Delire, C., Le Moigne, P., & Morin, S. (2016). Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model. *The Cryosphere*, 10(2), 853–877. https://doi.org/10.5194/tc-10-853-2016

Decharme, B., Delire, C., Minvielle, M., Colin, J., Vergnes, J.-P., Alias, A., et al. (2019). Recent changes in the ISBA-CTRIP land surface system for use in the CNRM-CM6 climate model and in global off-line hydrological applications. *Journal of Advances in Modeling Earth Systems*, 11(5), 1207–1252. https://doi.org/10.1029/2018MS001545

De Rosnay, P., Polcher, J. d., Bruen, M., & Laval, K. (2002). Impact of a physically based soil water flow and soil-plant interaction representation for modeling large-scale land surface processes. *Journal of Geophysical Research*, 107(D11), ACL–3. https://doi.org/10.1029/2001JD000634

Dolman, H., van der Molen, M., Parmentier, F.-J., Marchesini, L., Dean, J., van Huissteden, K., & Maximov, T. (2003–2014). FLUXNET2015 RU-Cok chokurdakh [Dataset]. https://doi.org/10.18140/FLX/1440182

FAO, IIASA, ISRIC, and ISSCAS. (2012). *JRC: Harmonized world soil database* (version 1.2). FAO, Rome, Italy and IIASA.

Farouki, O. (1981). *Thermal properties of soils* (Vol. 81(1)). US Army Corps of Engineers, Cold Regions Research and Engineering Laboratory.

FLUXNET2015-US-Prr. (2010–2014). FLUXNET2015 US-Prr poker flat research range Black spruce Forest [Dataset]. https://doi.org/10.18140/FLX/1440113

Forbes, T., Rossiter, D., & Van Wambeke, A. (1987). *Guidelines for evaluating the adequacy of soil resource inventories* (Vol. 4). Cornell University, Department of Agronomy.

Gaillard, R., Peylin, P., Cadule, P., Bastrikov, V., Chéruy, F., Cuynet, A., et al. (2025). Arctic soil carbon insulation averts large spring cooling from surface–atmosphere feedbacks. *Proceedings of the National Academy of Sciences*, 122(3), e2410226122. https://doi.org/10.1073/pnas.2410226122

Gelati, E., Decharme, B., Calvet, J.-C., Minvielle, M., Polcher, J., Fairbairn, D., & Weedon, G. P. (2018). Hydrological assessment of atmospheric forcing uncertainty in the Euro-Mediterranean area using a land surface model. *Hydrology and Earth System Sciences*, 22(4), 2091–2115. https://doi.org/10.5194/hess-22-2091-2018

Goodison, B., Louie, P., & Yang, D. (1998). *WMO solid precipitation measurement intercomparison* (Vol. 67). World Meteorological Organization-Publications-WMO TD.

Gouttevin, I., Krinner, G., Ciais, P., Polcher, J., & Legout, C. (2012). Multi-scale validation of a new soil freezing scheme for a land-surface model with physically-based hydrology. *The Cryosphere*, 6(2), 407–430. https://doi.org/10.5194/tc-6-407-2012

Greenland Ecosystem Monitoring. (2025a). *ClimateBasis: Meteorological standard observations, Zackenberg and Nuuk* [Dataset]. [CC-BY-SA-4.0]. *Greenland Ecosystem Monitoring*. Retrieved from https://data.g-e.m.dk

Greenland Ecosystem Monitoring. (2025b). *GeoBasis: Soil properties, Zackenberg and Nuuk* [Dataset]. [CC-BY-SA-4.0]. *Greenland Ecosystem Monitoring*. Retrieved from https://data.g-e.m.dk

Grossi, G., Lendvai, A., Peretti, G., & Ranzi, R. (2017). Snow precipitation measured by gauges: Systematic error estimation and data series correction in the central Italian Alps. *Water*, 9(7), 461. https://doi.org/10.3390/w9070461

Guimberteau, M., Zhu, D., Maignan, F., Huang, Y., Yue, C., Dantec-Nédélec, S., et al. (2018). ORCHIDEE-MICT (v8.4.1), a land surface model for the high latitudes: Model description and validation. *Geoscientific Model Development*, 11(1), 121–163. https://doi.org/10.5194/gmd-11-121-2018

Gupta, H. V., Kling, H., Yilmaz, K. K., & Martinez, G. F. (2009). Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling. *Journal of Hydrology*, 377(1–2), 80–91. https://doi.org/10.1016/j.jhydrol.2009.08.003

Harper, K. L., Lamarche, C., Hartley, A., Peylin, P., Ottlé, C., Bastrikov, V., et al. (2023). A 29-year time series of annual 300 m resolution plant-functional-type maps for climate models. *Earth System Science Data*, 15(3), 1465–1499. https://doi.org/10.5194/essd-15-1465-2023

Harris, I., Osborn, T. J., Jones, P., & Lister, D. (2020). Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. *Scientific Data*, 7(1), 109. https://doi.org/10.1038/s41597-020-0453-3

He, H., Liu, L., Dyck, M., Si, B., & Lv, J. (2021). Modelling dry soil thermal conductivity. *Soil and Tillage Research*, 213, 105093. https://doi.org/10.1016/j.still.2021.105093

Heaton, L., Fullen, M. A., & Bhattacharyya, R. (2016). Critical analysis of the Van Bemmelen conversion factor used to convert soil organic matter data to soil organic carbon data: Comparative analyses in a UK loamy sand soil. *Espaço Aberto*, 6(1), 35–44. https://doi.org/10.36403/e spacoaberto.2016.5244

Hengl, T., De Jesus, J. M., MacMillan, R. A., Batjes, N. H., Heuvelink, G. B., Ribeiro, E., et al. (2014). SoilGrids1km—Global soil information based on automated mapping. *PLoS One*, 9(8), e105992. https://doi.org/10.1371/journal.pone.0105992

Hengl, T., Mendes de Jesus, J., Heuvelink, G. B., Ruiperez Gonzalez, M., Kilibarda, M., Blagotić, A., et al. (2017). SoilGrids250m: Global gridded soil information based on machine learning. *PLoS One*, 12(2), e0169748. https://doi.org/10.1371/journal.pone.0169748

Huat, B. B., Kazemian, S., Prasad, A., & Barghchi, M. (2011). State of an art review of peat: General perspective. *International Journal of the Physical Sciences*, 6(8), 1988–1996. https://doi.org/10.5897/IJPS11.192

Hudson, B. D. (1994). Soil organic matter and available water capacity. *Journal of Soil and Water Conservation*, 49(2), 189–194. https://doi.org/10.1080/00224561.1994.12456850

Hugelius, G., Loisel, J., Chadburn, S., Jackson, R. B., Jones, M., MacDonald, G., et al. (2020). Large stocks of peatland carbon and nitrogen are vulnerable to permafrost thaw. *Proceedings of the National Academy of Sciences*, 117(34), 20438–20446. https://doi.org/10.1073/pnas.1916387117

Hugelius, G., Ramage, J., Burke, E., Chatterjee, A., Smallman, T., Aalto, T., et al. (2024). Permafrost region greenhouse gas budgets suggest a weak CO2 sink and CH4 and N2O sources, but magnitudes differ between top-down and bottom-up methods. *Global Biogeochemical Cycles*, 38(10), e2023GB007969. https://doi.org/10.1029/2023GB007969

Hugelius, G., Strauss, J., Zubrzycki, S., Harden, J. W., Schuur, E. A., Ping, C.-L., et al. (2014). Estimated stocks of circumpolar permafrost carbon with quantified uncertainty ranges and identified data gaps. *Biogeosciences*, 11(23), 6573–6593. https://doi.org/10.5194/bg-11-6573-2014

Hugelius, G., Tarnocai, C., Broll, G., Canadell, J., Kuhry, P., & Swanson, D. (2013). The Northern Circumpolar Soil Carbon Database: Spatially distributed datasets of soil coverage and soil carbon storage in the northern permafrost regions. *Earth System Science Data*, 5(1), 3–13. https://doi.org/10.5194/essd-5-3-2013

Ikawa, H., Nakai, T., Busey, R. C., Kim, Y., Kobayashi, H., Nagai, S., et al. (2015). Understory CO2, sensible heat, and latent heat fluxes in a black spruce forest in interior Alaska. *Agricultural and Forest Meteorology*, 214, 80–90. https://doi.org/10.1016/j.agrformet.2015.08.247

Johansen, O. (1975). *Thermal conductivity of soils* (PhD thesis). Norwegian Institute of Technology.

Koven, C. D., Friedlingstein, P., Ciais, P., Khvorostyanov, D., Krinner, G., & Tarnocai, C. (2009). On the formation of high-latitude soil carbon stocks: Effects of cryoturbation and insulation by organic matter in a land surface model. *Geophysical Research Letters*, 36(21), https://doi.org/10.1029/2009GL040150
Koven, C. D., Riley, W. J., & Stern, A. (2013). Analysis of permafrost thermal dynamics and response to climate change in the CMIP5 earth system models. *Journal of Climate*, 26(6), 1877–1900. https://doi.org/10.1175/JCLI-D-12-00228.1

Lawrence, D. M., & Slater, A. G. (2008). Incorporating organic soil into a global climate model. *Climate Dynamics*, 30(2–3), 145–160. https://doi.org/10.1007/s00382-007-0278-1

Letts, M. G., Roulet, N. T., Comer, N. T., Skarupa, M. R., & Verseghy, D. L. (2000). Parametrization of peatland hydraulic properties for the Canadian land surface scheme. *Atmosphere-Ocean*, 38(1), 141–160. https://doi.org/10.4324/9781315170206-8

Liljedahl, A. K., Hinzman, L. D., Kane, D. L., Oechel, W. C., Tweedie, C. E., & Zona, D. (2017). Tundra water budget and implications of precipitation underestimation. *Water Resources Research*, 53(8), 6472–6486. https://doi.org/10.1002/2016WR020001

Lohila, A., Aurela, M., Tuovinen, J.-P., Hatakka, J., & Laurila, T. (2000–2003). FLUXNET2015 FI-Jok jokioinen [Dataset]. https://doi.org/10.18140/FLX/1440159

Lohila, A., Aurela, M., Tuovinen, J.-P., & Laurila, T. (2004). Annual CO2 exchange of a peat field growing spring barley or perennial forage grass. *Journal of Geophysical Research*, 109(D18). https://doi.org/10.1029/2004JD004715

López-Blanco, E., Jackowicz-Korczyński, M., Mastepanov, M., Skov, K., Westergaard-Nielsen, A., Williams, M., & Christensen, T. R. (2020). Multi-year data-model evaluation reveals the importance of nutrient availability over climate in arctic ecosystem C dynamics. *Environmental Research Letters*, 15(9), 094007. https://doi.org/10.1088/1748-9326/ab865b

Lund, M., Hansen, B. U., Pedersen, S. H., Stiegler, C., & Tamstorf, M. P. (2014). Characteristics of summer-time energy exchange in a high Arctic tundra heath 2000–2010. *Tellus B: Chemical and Physical Meteorology*, 66(1), 21631. https://doi.org/10.3402/tellusb.v66.21631

Mammarella, I., Keronen, P., Kolari, P., Launiainen, S., Pumpanen, J., Rannik, Ü., et al. (1996–2014). FLUXNET2015 FI-Hyy hyytiala [Dataset]. https://doi.org/10.18140/FLX/1440158

McGaw, R. (1969). Heat conduction in saturated granular materials. Highway Research Board Special Report (p. 103).

McGuire, A. D., Anderson, L. G., Christensen, T. R., Dallimore, S., Guo, L., Hayes, D. J., et al. (2009). Sensitivity of the carbon cycle in the Arctic to climate change. *Ecological Monographs*, 79(4), 523–555. https://doi.org/10.1890/08-2025.1

Merbold, L., Kutsch, W. L., Corradi, C., Kolle, O., Rebmann, C., Stoy, P. C., et al. (2009). Artificial drainage and associated carbon fluxes (CO₂/CH₄) in a tundra ecosystem. *Global Change Biology*, 15(11), 2599–2614. https://doi.org/10.1111/j.1365-2486.2009.01962.x

Merbold, L., Rebmann, C., & Corradi, C. (2002–2005). FLUXNET2015 RU-Che cherski [Dataset]. https://doi.org/10.18140/FLX/1440181

Miner, K. R., Turetsky, M. R., Malina, E., Bartsch, A., Tamminen, J., McGuire, A. D., et al. (2022). Permafrost carbon emissions in a changing Arctic. *Nature Reviews Earth and Environment*, 3(1), 55–67. https://doi.org/10.1038/s43017-021-00230-3

Obu, J. (2021). How much of the earth's surface is underlain by permafrost? *Journal of Geophysical Research: Earth Surface*, 126(5), e2021JF006123. https://doi.org/10.1029/2021JF006123

Pan, X., Yang, D., Li, Y., Barr, A., Helgason, W., Hayashi, M., et al. (2016). Bias corrections of precipitation measurements across experimental sites in different ecoclimatic regions of western Canada. *The Cryosphere*, 10(5), 2347–2360. https://doi.org/10.5194/tc-10-2347-2016

Paquin, J.-P., & Sushama, L. (2015). On the Arctic near-surface permafrost and climate sensitivities to soil and snow model formulations in climate models. *Climate Dynamics*, 44(1–2), 203–228. https://doi.org/10.1007/s00382-014-2185-6

Pastorello, G., Trotta, C., Canfora, E., Chu, H., Christianson, D., Cheah, Y.-W., et al. (2020). The FLUXNET2015 dataset and the ONEflux processing pipeline for eddy covariance data. *Scientific Data*, 7(1), 225. https://doi.org/10.1038/s41597-020-0534-3

Peters-Lidard, C., Blackburn, E., Liang, X., & Wood, E. F. (1998). The effect of soil thermal conductivity parameterization on surface energy fluxes and temperatures. *Journal of the Atmospheric Sciences*, 55(7), 1209–1224. https://doi.org/10.1175/1520-0469(1998)055<1209:teo stc>2.0.co;2

Poggio, L., de Sousa, L. M., Batjes, N. H., Heuvelink, G. B. M., Kempen, B., Ribeiro, E., & Rossiter, D. (2021). SoilGrids 2.0: Producing soil information for the globe with quantified spatial uncertainty. *SOIL*, 7(1), 217–240. https://doi.org/10.5194/soil-7-217-2021

Post, W. M., Emanuel, W. R., Zinke, P. J., & Stangenberger, A. G. (1982). Soil carbon pools and world life zones. *Nature*, 298(5870), 156–159. https://doi.org/10.1038/298156a0

Rantanen, M., Karpechko, A. Y., Lipponen, A., Nordling, K., Hyvärinen, O., Ruosteenoja, K., et al. (2022). The Arctic has warmed nearly four times faster than the globe since 1979. *Communications Earth and Environment*, 3(1), 168. https://doi.org/10.1038/s43247-022-00498-3

Rawls, W., Pachepsky, Y. A., Ritchie, J., Sobecki, T., & Bloodworth, H. (2003). Effect of soil organic carbon on soil water retention. *Geoderma*, 116(1–2), 61–76. https://doi.org/10.1016/S0016-7061(03)00094-6

Rinke, A., Kuhry, P., & Dethloff, K. (2008). Importance of a soil organic layer for Arctic climate: A sensitivity study with an Arctic RCM. *Geophysical Research Letters*, 35(13). https://doi.org/10.1029/2008GL034052

Robinson, D., Thomas, A., Reinsch, S., Lebron, I., Feeney, C., Maskell, L., et al. (2022). Analytical modelling of soil porosity and bulk density across the soil organic matter and land-use continuum. *Scientific Reports*, 12(1), 7085. https://doi.org/10.1038/s41598-022-11099-7

Schneider, U., Becker, A., Finger, P., Meyer-Christoffer, A., & Ziese, M. (2018). GPCC full data monthly product version 2018 at 0.5°: Monthly land-surface precipitation from rain-gauges built on GTS-based and historical data. Global Precipitation Climatology Centre. https://doi.org/10.5676/DWD_GPCC/FD_M_V2018_050

Sokratov, S. A., & Barry, R. G. (2002). Intraseasonal variation in the thermoinsulation effect of snow cover on soil temperatures and energy balance. *Journal of Geophysical Research*, 107(D10), ACL-13. https://doi.org/10.1029/2001JFD000489

Strand, S. M., Christiansen, H. H., Johansson, M., Åkerman, J., & Humlum, O. (2021). Active layer thickening and controls on interannual variability in the nordic arctic compared to the circum-arctic. *Permafrost and Periglacial Processes*, 32(1), 47–58. https://doi.org/10.1002/ppp.2088

Sun, J., Yang, K., Lu, H., Zhou, X., Li, X., Chen, Y., et al. (2023). Land–atmosphere feedbacks weaken the cooling effect of soil organic matter property toward deep soil on the eastern Tibetan Plateau. *Journal of Hydrometeorology*, 24(1), 105–117. https://doi.org/10.1175/JHM-D-22-0074.1

Sun, S., Chen, B., Ge, M., Qu, J., Che, T., Zhang, H., et al. (2016). Improving soil organic carbon parameterization of land surface model for cold regions in the Northeastern Tibetan Plateau, China. *Ecological Modelling*, 330, 1–15. https://doi.org/10.1016/j.ecolmodel.2016.03.014

Suni, T., Rinne, J., Reissell, A., Altimir, N., Keronen, P., Rannik, U., et al. (2003). Long-term measurements of surface fluxes above a scots pine forest in Hyytiala, southern Finland, 1996–2001. *Boreal Environment Research*, 8(4), 287–302.

Thum, T., Aalto, T., Laurila, T., Aurela, M., Kolari, P., & Hari, P. (2007). Parametrization of two photosynthesis models at the canopy scale in a northern boreal Scots pine forest. *Tellus B: Chemical and Physical Meteorology*, 59(5), 874–890. https://doi.org/10.1111/j.1600-0889.2007.00305.x

Tran, A. P., Dafflon, B., & Hubbard, S. S. (2017). Coupled land surface–subsurface hydrogeophysical inverse modeling to estimate soil organic carbon content and explore associated hydrological and thermal dynamics in the Arctic tundra. *The Cryosphere*, 11(5), 2089–2109. https://doi.org/10.5194/tc-11-2089-2017
Van der Molen, M., Van Huissteden, J., Parmentier, F.-J. W., Petrescu, A., Dolman, A., Maximov, T. C., et al. (2007). The growing season greenhouse gas balance of a continental tundra site in the Indigirka lowlands, NE Siberia. *Biogeosciences*, *4*(6), 985–1003. https://doi.org/10.5194/bg-4-985-2007

Von Fromm, S. F., Hoyt, A. M., Sierra, C. A., Georgiou, K., Doetterl, S., & Trumbore, S. E. (2024). Controls and relationships of soil organic carbon abundance and persistence vary across pedo-climatic regions. *Global Change Biology*, *30*(5), e17320. https://doi.org/10.1111/gcb.17320

Vuichard, N., Messina, P., Luyssaert, S., Guenet, B., Zaehle, S., Ghattas, J., et al. (2019). Accounting for carbon and nitrogen interactions in the global terrestrial ecosystem model ORCHIDEE (trunk version, rev 4999): Multi-scale evaluation of gross primary production. *Geoscientific Model Development*, *12*(11), 4751–4779. https://doi.org/10.5194/gmd-12-4751-2019

Vuichard, N., & Papale, D. (2015). Filling the gaps in meteorological continuous data measured at FLUXNET sites with ERA-Interim reanalysis. *Earth System Science Data*, *7*(2), 157–171. https://doi.org/10.5194/essd-7-157-2015

Wang, Y., Xie, Z., Jia, B., Wang, L., Li, R., Liu, B., et al. (2020). Sensitivity of snow simulations to different atmospheric forcing data sets in the land surface model CAS-LSM. *Journal of Geophysical Research: Atmospheres*, *125*(16), e2019JD032001. https://doi.org/10.1029/2019JD032001

Weedon, G. P., Balsamo, G., Bellouin, N., Gomes, S., Best, M. J., & Viterbo, P. (2014). The WFDEI meteorological forcing data set: WATCH Forcing Data methodology applied to ERA-Interim reanalysis data. *Water Resources Research*, *50*(9), 7505–7514. https://doi.org/10.1002/2014WR015638

Yuan, Y., Ma, Y., Zuo, H., Yang, C., Yuan, L., & Chen, J. (2021). Modification and comparison of thermal and hydrological parameterization schemes for different underlying surfaces on the Tibetan Plateau in the warm season. *Journal of Geophysical Research: Atmospheres*, *126*(22), e2021JD035177. https://doi.org/10.1029/2021JD035177

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

## References From the Supporting Information

Carsel, R. F., & Parrish, R. S. (1988). Developing joint probability distributions of soil water retention characteristics. *Water Resources Research*, *24*(5), 755–769. https://doi.org/10.1029/WR024i005p00755

Pielke Sr, R. A. (2002). *Mesoscale meteorological modeling* (2nd ed., Vol. 78). Academic Press.