The use of radiocarbon ¹⁴C to constrain carbon dynamics in the soil module of the land surface model ORCHIDEE (SVN r5165)

Despite the importance of soil as a large component of the terrestrial ecosystem, the soil compartments are not well represented in land surface models (LSMs). Indeed, soils in current LSMs are generally represented based on a very simplified schema that can induce a misrepresentation of the deep dynamics of soil carbon. Here, we present a new version of the Institut Pierre Simon Laplace (IPSL) LSM called ORCHIDEE-SOM (ORganizing Carbon and Hydrology in Dynamic EcosystEms-Soil Organic Matter), incorporating the ¹⁴C dynamics into the soil. ORCHIDEE-SOM first simulates soil carbon dynamics for different layers, down to 2 m depth. Second, concentration of dissolved organic carbon and its transport are modelled. Finally, soil organic carbon decomposition is considered taking into account the priming effect.

After implementing ¹⁴C in the soil module of the model, we evaluated model outputs against observations of soil organic carbon and modern ¹⁴C fraction (F¹⁴C) for different sites with different characteristics. The model managed to reproduce the soil organic carbon stocks and the F¹⁴C along the vertical profiles for the sites examined. However, an overestimation of the total carbon stock was noted, primarily on the surface layer. Due to ¹⁴C, it is possible to probe carbon age in the soil, which was found to be underestimated. Thereafter, two different tests on this new version have been established. The first was to increase carbon residence time of the passive pool and decrease the flux from the slow pool to the passive pool. The second was to establish an equation of diffusion, initially constant throughout the profile, making it vary exponentially as a function of depth. The first modifications did not improve the capacity of the model to reproduce observations, whereas the second test improved both estimation of surface soil carbon stock as well as soil carbon age. This demonstrates that we should focus more on vertical variation in soil parameters as a function of depth, in order to upgrade the representation of the global carbon cycle in LSMs, thereby helping to improve predictions of the of soil organic carbon to environmental changes.

1 Introduction

The complexity of the mechanisms involved in controlling soil activity (Jastrow et al., 2007) and therefore the carbon flux from the soil to the atmosphere makes predicting the response of these systems to climate change extremely complex. Thus our ability to predict future changes in carbon stocks in soils using global climate models is currently heavily criticized (Todd-Brown et al., 2013; Wieder et al., 2013). Indeed, Earth system models (ESMs) are increasingly used today in order to predict the future evolution of the climate. For instance, results of a set of ESMs are taken into account within the Intergovernmental Panel on Climate Change (IPCC; Taylor et al., 2012) for assessment of the impacts of climate change and design of mitigation strategies. Hence, their predictions need to be as accurate as possible. These models represent the physical, chemical and biological processes within and between the atmosphere, ocean and terrestrial biosphere. They allow us to follow and understand both the effect of the climate on carbon storage and vice versa. However, ESMs are continuously under development and some key processes in the global carbon cycle are still missing or not represented with the necessary details. One of the components of an ESM is the land surface model (LSM). This component primarily manages the carbon cycle, energy and water on land and simulates the carbon exchange between the land surface and the atmosphere, namely the gross primary production (GPP) and the autotrophic and heterotrophic respiration.

Despite the importance of soils as a large component of the global carbon storage, soil compartments are not well represented in LSMs (Todd-Brown et al., 2013). Indeed, carbon dynamics in soil described in LSMs are based on the “Century” (Parton et al., 1987) or RothC models (Coleman et al., 1997) where soil carbon is represented as several pools with different turnover rates for each pool. Carbon is decomposed in each pool, one part of which is then transferred from one pool to another and the other part is lost through heterotrophic respiration. In addition, soils are generally represented as a single-layer box in LSMs that do not take into account the evolution and variation in soil organic processes as a function of depth (Todd-Brown et al., 2013).

One way to reconcile this simplified representation of carbon dynamics of the models with the complexity of the data collected in the field is to integrate isotopic tracers into the models themselves and thus facilitate the comparison between model outputs and data (He et al., 2016). Moreover, thanks to an additive constraint on the model structure, this may improve the model performance. For instance, radiocarbon is an important tool for studying the dynamics of soil organic matter (Trumbore, 2000). Indeed, ¹⁴C data acquired from soil organic matter provide complementary information on the dynamics (temporal dimension) of soil organic matter. This tracer has the major advantage of being integrator of carbon dynamics on long timescales (a few decades to several centuries). It is therefore a very powerful tool to constrain conceptual schemes that may not be directly compared to variables measured in the field (Elliott et al., 1996). Different authors have already successfully implemented radiocarbon in soil models and were able to clearly show that the introduction of pools with turnover time of thousands of years were unnecessary to fit radiocarbon data (Ahrens et al., 2015), whereas Braakhekke et al. (2014) showed that after a re-parameterization of the models based on radiocarbon data the prediction of their model was quite different with more carbon in topsoil and less in deep soil compared to the model without radiocarbon.

Radiocarbon is produced naturally at a constant rate in the upper atmosphere through bombardment of cosmic rays. It thus provides information on the dynamics of organic matter that has been stabilized by interaction with mineral surfaces and stored long enough for significant radioactive decay (Trumbore, 2000), as the half-life of ¹⁴C is about 5730 years. We must also take into account radiocarbon produced during atmospheric tests of thermonuclear weapons in the early 1960s (Delibrias et al., 1964; Hua et al., 2013). Atmospheric bomb testing in the late 1950s and early 1960s led to an abrupt doubling of atmospheric ¹⁴C concentration in a span of 2–3 years. Through exchange with ocean and terrestrial reservoirs, it has decreased but still remains above the natural background. As with any other carbon isotope, this ¹⁴C was metabolized by the vegetation and transferred to soil. By measuring ¹⁴C activity of a soil sample, it is possible to evaluate the amount of carbon introduced into the soil since the 1960s (Balesdent and Guillet, 1982; Scharpenseel and Schiffmann, 1977).

In this study, we present a new version of the IPSL LSM called ORCHIDEE-SOM (ORganizing Carbon and Hydrology in Dynamic EcosystEms-Soil Organic Matter) incorporating ¹⁴C dynamics in the soil. Thanks to this tracer, we can evaluate the SOC dynamics, in particular by looking at the ¹⁴C peak produced by atmospheric weapons testing and observed in the soils at four different sites having different biomes.

2 Materials and methods

2.1 ORCHIDEE-SOM overview

ORCHIDEE is the LSM of the IPSL Earth system model (Krinner et al., 2005). It is composed of three different modules. First, SECHIBA (Ducoudré et al., 1993; de Rosnay and Polcher, 1998), the surface–vegetation–atmosphere transfer scheme, describes the soil water budget and energy and water exchanges. The time step of this module is 30 min. Second, the module of the vegetation dynamics has been taken from the dynamic global vegetation model LPJ (Sitch et al., 2003). The time step of this module is 1 year. Finally, the STOMATE (Saclay Toulouse Orsay Model for the Analysis of Terrestrial Ecosystems) module simulates vegetation phenology and carbon dynamics with a time step of 1 day.

ORCHIDEE can be run coupled to a global circulation model where the boundary conditions of the model are provided by the atmospheric modules (temperature, precipitation, atmospheric CO₂ concentration, etc.). In return, ORCHIDEE provides the land surface carbon, energy and water fluxes. However, since our study focuses on changes in the land surface rather than on the interaction with climate, we ran ORCHIDEE in the offline configuration. In this case, atmospheric conditions such as temperature, humidity and wind are read from a meteorological dataset. The climate data CRUNCEP used for our study (6-hourly climate data over several years) were obtained from the combination of two existing datasets: the Climate Research Unit (CRU; Mitchell et al., 2004) and the National Centers for Environmental Prediction (NCEP; Kalnay et al., 1996).

Figure 1 Overview of the different fluxes and processes in soil as presented in the version of ORCHIDEE-SOM adapted from Camino-Serrano et al. (2018).

Our starting point is a ORCHIDEE-SOM version based on the SVN r3340 (Krinner et al., 2005), which is presented in detail in Camino-Serrano et al. (2017). Figure 1 represents how the soil is described in this new version. Indeed, the major particularity of ORCHIDEE-SOM is that it simulates the dynamics of soil carbon for 11 layers from the surface to 2 m depth. First, litter is divided into four pools: metabolic or structural litter pools which can be found below or above ground. Only the belowground litter is modelled on 11 levels, from surface to 2 m depth, as the aboveground litter layer has a fixed thickness of 10 mm. Second, SOC is divided into three pools (active, passive and slow), following Parton et al. (1988), which differ in their turnover rates and which are discretized into 11 layers up to a depth of 2 m. Then, dissolved organic carbon (DOC) is represented as two pools and also discretized over 11 layers up to a depth of 2 m: labile DOC has a high decomposition rate and recalcitrant DOC has a low decomposition rate (Camino-Serrano et al., 2018). Finally, another particularity of this version of ORCHIDEE-SOM is that the SOC decomposition is modified to account for the priming effect following Guenet et al. (2016). Briefly, priming is described following Eq. (1).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial {\mathrm{SOC}}_{i,z}}{\partial t}}={\mathrm{DOC}}_{\mathrm{Recycled},i,j}\left(t\right)-k_{\mathrm{SOC},i}\\ \text{(1)}& {\displaystyle }& {\displaystyle }\times \left(\mathrm{1}-e^{-c\times {\mathrm{LOC}}_z\left(t\right)}\right)\times \mathrm{SOC}{\left(t\right)}_{i,z}\times \mathit{\theta }\left(t\right)\times \mathit{\tau }\left(t\right)\end{array}$$

with DOC_recycled being the unrespired DOC that is redistributed into the pool i considered for each soil layer z in g C m⁻² day⁻¹, k_SOC being a SOC decomposition rate constant (days⁻¹), and LOC being the stock of labile organic C defined as the sum of the C pools with a higher decomposition rate than the pool considered within each soil layer z. We therefore considered that for the active carbon pool LOC is the litter and DOC, but for the slow carbon pool LOC is the sum of the litter, DOC and so on. Finally, c is a parameter controlling the impact of the LOC pool on the SOC mineralization rate, i.e. the priming effect. The equation was parameterized based on soil incubations data and evaluated over litter manipulation experiments (Guenet et al., 2016).

Since the soil profile is divided into 11 layers, SOC and DOC transport following the diffusion must also be described. SOC diffusion is actually a representation of bioturbation processes (animal and plant activity), whereas DOC relies more on non-biological diffusion. Both diffuse through concentration gradients.

This is represented using Fick's law (Braakhekke et al., 2011; Elzein and Balesdent, 1995; O'Brien and Stout, 1978; Wynn et al., 2005):

$$\begin{array}{}\text{(2)}& F_{\mathrm{D}}=-D\times {\displaystyle \frac{{\partial }^{\mathrm{2}}C}{\partial z^{\mathrm{2}}}},\end{array}$$

where F_D is the flux of carbon transported by diffusion in g C m⁻³ day⁻¹, D is the diffusion coefficient (m² day⁻¹) and C is the amount of carbon in the pool (DOC or SOC) subject to transport (g C m⁻³). The diffusion coefficient is assumed to be constant across the soil profile in ORCHIDEE-SOM but the diffusion parameters (D) used in the equations for SOC and DOC can differ. All the transport processes goes up to 2 m, corresponding to the soil depth fixed in the model. For DOC, at 2 m the DOC can be exported through drainage.

2.2 ORCHIDEE-SOM-¹⁴C

In ORCHIDEE-SOM, the different compartments (soil carbon input, litter, SOC, DOC and heterotrophic respiration) are presented as a matrix with a single dimension referring to the total carbon. In order to introduce the ¹⁴C, a new dimension has been added to all the variables cited above. Thus, all processes that apply to the total soil carbon are now also represented for ¹⁴C. We label this new version including ¹⁴C as ORCHIDEE-SOM-¹⁴C.

Several ways of reporting ¹⁴C activity levels are available. We chose to use the fraction modern, with the F¹⁴C symbol as advocated by Reimer et al. (2004) rather than absolute concentration of ¹⁴C (reported as Bq).

$$\begin{array}{ll}{\displaystyle }{\mathrm{F}}^{\mathrm{14}}\mathrm{C}& {\displaystyle }=\left({\displaystyle \frac{A_{\mathrm{S}}}{\mathrm{0.95}{A}_{\mathrm{OX}\mathrm{1}}}}\right)\times \left({\displaystyle \frac{\mathrm{0.975}}{\mathrm{0.981}}}\right)\\ \text{(3)}& {\displaystyle }& {\displaystyle }\times \left[\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{OX}\mathrm{1}}}{\mathrm{1000}}}\right)/\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{S}}}{\mathrm{1000}}}\right)\right]\end{array}$$

with $A{=}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{12}}$C; S for sample and OX1 for “Oxalic Acid 1”, the ¹⁴C international standard. F¹⁴C is twice normalized: (i) it takes into account isotopic fractionation by being normalized to a δ¹³C = −25 ‰, and (ii) it corresponds to a deviation towards an international standard (i.e. 95 % of OX1 as measured in 1950; Stuiver and Polach, 1977). By propagating F¹⁴C from atmosphere at the origin of vegetal photosynthesis to soil respired CO₂, there is no need to focus on ¹³C isotopic fractionation all along the organic matter mineralization with F¹⁴C.

To ease the readability of the paper, we will further express F¹⁴C as F¹⁴C $=A_{\mathrm{sample}}/A_{\mathrm{ref}}$ with A_sample being the A of the measured (or modelled) data and A_ref an international reference. Normalizations are included in A_ref and F¹⁴C will be written as F¹⁴ to simplify notation involving superscripts and subscripts.

Since we focus on SOC dynamics, we did not include the ¹⁴C in plants but did include ¹⁴C in the litter. The ¹⁴C-litter is obtained by multiplying the atmospheric value by the total carbon in the litter:

$$\begin{array}{}\text{(4)}& \mathrm{Litter}\left({}^{\mathrm{14}}\mathrm{C}\right)=F_{\mathrm{atm}}^{\mathrm{14}}\times \mathrm{Litter}\left(C\right),\end{array}$$

where F${}_{\mathrm{atm}}^{\mathrm{14}}$ is the F¹⁴C of atmosphere at the time of leaf growth (Fig. 2).

Figure 2 Evolution of the F¹⁴C of atmospheric CO₂ in Argentina, Republic of the Congo and France (data from Hua et al., 2013).

Thus, from the litter, all processes defined in Sect. 2.1 that apply to total soil carbon are also represented for ¹⁴C.

We also take into account the radioactive decay of ¹⁴C. For that, we calculate the amount of ¹⁴C as follows:

$$\begin{array}{}\text{(5)}& {}^{\mathrm{14}}\mathrm{C}={}^{\mathrm{14}}\mathrm{C}-K_{\mathrm{decrease}}{\times }^{\mathrm{14}}\mathrm{C},\end{array}$$

where K_decrease is the radioactive decay constant (= Ln2/5730; Godwin, 1962).

The F¹⁴C of the soil is then calculated back for carbon, per pool:

$$\begin{array}{}\text{(6)}& {\mathrm{F}}_{\mathrm{Pool},z}^{\mathrm{14}}={\displaystyle \frac{{}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{Pool},z}}{{\mathrm{C}}_{\mathrm{Pool},z}}}\end{array}$$

with Pool representing the active, slow or passive pool.

Finally, we calculate a mean F¹⁴C value per soil layer, according to depth:

$$\begin{array}{}\text{(7)}& {\mathrm{F}}_{\mathrm{Mean},z}^{\mathrm{14}}={\displaystyle \frac{{\displaystyle \begin{array}{c}{\mathrm{F}}_{\mathrm{active},z}^{\mathrm{14}}{\times }^{\mathrm{14}}{\mathrm{C}}_{\mathrm{active},z}+{\mathrm{F}}_{\mathrm{slow},z}^{\mathrm{14}}{\times }^{\mathrm{14}}{\mathrm{C}}_{\mathrm{slow},z}\\ +{\mathrm{F}}_{\mathrm{passive},z}^{\mathrm{14}}{\times }^{\mathrm{14}}{\mathrm{C}}_{\mathrm{passive},z}\end{array}}}{{}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{active},z}{+}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{slow},z}{+}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{passive},z}}}.\end{array}$$

2.3 Site descriptions

2.3.1 French sites

Two Luvisol (WRB, 2006) profiles located in northern France were selected: the Feucherolles and Mons sites. In Mons (49.87^∘ N, 3.03^∘ E; Luvisol) the soils sit under grassland, and are developed from several metres of loess and are therefore well drained. The mean annual air temperature is 11 ^∘C and the annual precipitation is about 680 mm (Keyvanshokouhi et al., 2016). In Feucherolles (48.9^∘ N, 1.97^∘ E), the soil sits under oak forest, and clay and gritstone deposits are found at approximately 1.5 m depth. The mean annual air temperature is 11.2 ^∘C and the annual precipitation is about 660 mm (Keyvanshokouhi et al., 2016). Both soils are neutral to slightly acidic and are characterized by the presence of a clay accumulation Bt horizon with clay content reaching 30 % for Feucherolles and 27 % for Mons, while the upper horizons are poorer in clay (17 % for Feucherolles and 20 % for Mons).

The ¹⁴C data from the soils of both sites were obtained after chemical treatment done at Laboratoire des Sciences du Climat et de l'Environnement (LSCE) using a protocol adapted to achieve carbonate leaching without any loss of organic carbon; ¹⁴C activity was measured by AMS (aerosol mass spectrometry) at the French Laboratoire de mesure du ¹⁴C (LMC14) facility (Cottereau et al., 2007). Details on measurements and sampling can be found in Jagercikova et al. (2017).

2.3.2 Republic of the Congo site

The studied site is located in Kissoko (4.35^∘ S, 11.75^∘ E). It belongs to the SOERE F-ORE-T (Site de l'ObservatoirE de Recherche en Environnement sur le Fonctionnement des écosystèmes fOREsTiers) field observation sites of Pointe Noire, Republic of the Congo. The mean annual air temperature is about 25 ^∘C with low seasonal variation (±5 ^∘C), and average annual precipitation of 1400 mm, and a dry season between June and September. The deep acidic sandy soil is a ferralic Arenosol (WRB, 2006). The soil is characterized by a sand content greater than 90 % (Laclau et al., 2000). A soil profile was taken under native savanna vegetation dominated by C₄ plants (Epron et al., 2009). The soil was sampled in May 2014 at different depths: 0–5, 5–10, 10–15, 15–20, 20–30, 30–40, 40–50, 50–60, 60–80, 80–100 and 100–120 cm. All samples were crushed and air-dried. Once in the laboratory, they were homogenized, crushed, randomly subsampled and sieved at 200 µm. Then ¹⁴C measurements were made the same way as the two French sites, using the LSCE chemical treatment and the French LMC14 facility following recommendations by Cottereau et al. (2007).

2.3.3 Argentinian site

The Province of Misiones is located in northeastern Argentina. The climate is subtropical humid without a dry season, an annual mean temperature of 20 ^∘C and 1850 mm of mean annual rainfall (Morrás et al., 2009). The profile used in this study is located in the southern part of Misiones (27^∘ S, 55^∘ W). Native vegetation is a forest dominated by C₃ plants. The soil selected is an Acrisol (WRB, 2006). It is a red clay soil, strongly to very strongly acid with a clay content varying from 40 % at the surface to 60 % at 1 m depth. The ¹⁴C measurements were made using a new compact radiocarbon system called ECHoMICADAS (Environment, Climate, Human, Mini Carbon Dating System; Tisnérat-Laborde et al., 2015). Briefly, the soil was sampled in May 2015 at different depths: 0–5, 5–10, 10–15, 15–20, 20–30, 30–40, 40–50, 50–60, 60–80 and 80–100 cm. All samples were crushed and air-dried. Once in the laboratory, they were homogenized, crushed, randomly subsampled and sieved at 200 µm. Then ¹⁴C measurements were made using the ECHoMICADAS following the recommendations of Tisnérat-Laborde et al. (2015).

For the four sites, the SOC (kg m⁻³), for each depth z, was calculated using carbon content and bulk density data using the following equation:

$$\begin{array}{}\text{(8)}& {\mathrm{SOC}}_z={\mathrm{OCC}}_z\times {\mathrm{BD}}_z,\end{array}$$

where OCC (wt/wt) is the carbon content and BD (kg m⁻³) is the bulk density.

2.4 Different model tests

After the implementation of radiocarbon in the model, different tests were carried out (Table 2). Here we represent the outputs provided by three simulations:

First, simulation using the initial version ORCHIDEE-SOM-¹⁴C (labelled “Control” in figures and tables) in which no changes were made. The diffusion was kept constant throughout the profile ($D={\mathrm{1.10}}^{-\mathrm{4}}$ m² yr⁻¹) and the other parameters are those of the detailed version in Camino-Serrano et al. (2018).

Second, simulation using the initial version ORCHIDEE-SOM-¹⁴C in which we modified some parameters following He et al. (2016) (labelled “He et al., 2016, parameterization” in figures and tables). In brief, the authors used ¹⁴C data from 157 globally distributed soil profiles sampled to 1 m depth to evaluate CMIP5 models. Their results show that ESMs underestimated the mean age of soil carbon by a factor of more than 6 and overestimated the carbon sequestration potential of soils by a factor of nearly 2. So, the suggestion (that we apply in this simulation) for the IPSL model was to multiply the turnover time of the passive pool by 14 and the flux from slow pool to passive pool by 0.07 (Table 2). The diffusion was kept constant throughout the profile ($D={\mathrm{1.10}}^{-\mathrm{4}}$ m² yr⁻¹) but the turnover time of the passive pool increased from 462 to 6468 years and the flux from the slow pool to the passive pool decreased from 0.07 to 0.0049.

Third, simulation using the initial version ORCHIDEE-SOM-¹⁴C in which we assume that the diffusion varies as a function of depth (“Depth-varying diffusion constant” in figures and tables) according to the equation below:

$$\begin{array}{}\text{(9)}& D\left(z\right)={\mathrm{5.42.10}}^{-\mathrm{4}}{e}^{(-\mathrm{0.04}z)},\end{array}$$

where D is the diffusion (m² yr⁻¹) at a specific depth and z is the depth. This equation of diffusion varying as a function of depth is following Jagercikova et al. (2014) and assumes that bioturbation is higher in the topsoil than in deep soil.

Table 1 General description of the studied sites. The mean bulk density, pH and clay fraction values calculated from the different soil layer depths available from the data were used as input for each site. For the Mons and Feucherolles sites, min and max values of pH and clay fraction are provided between brackets.

Table 2 The main differences between the three simulations.

2.5 Model simulations

In order to reach a steady state of the soil module, we ran the model over 12 700 years (spin-up). The state at the last time step of this spin-up was used as the initial state for the simulations. For this, the CRUNCEP meteorological data for the period 1901–1910 were used. This has been applied for Misiones, Feucherolles and Mons. However, for Kissoko, a first spin-up similar to the other sites was carried out but a second one (over approximately 4200 years) was also done after the end of the first to take into account the change of the land cover from a tropical forest to a C₄ savanna at this site (Schwartz et al., 1992). The atmospheric CO₂ concentration has been set to 296 ppm (year 1901; Keeling and Whorf, 2006) for the spin-ups and the F¹⁴C has been set to one corresponding to pre-industrial values. For each site, specific pH, clay content and bulk density values were used (Table 1). It should be noted that for these last data, only one value (the mean value on the profile) is provided as input for the model.

The simulations were outputted at a yearly time step, from 1900 to 2011. A yearly atmospheric CO₂ concentration value (Keeling and Whorf, 2006) is read for the sites. The same specific pH, clay content and bulk density values were used (Table 1).

Figure 2 shows the evolution of the F¹⁴C values in the atmosphere used in our model for Argentina, Republic of the Congo and France (Fig. 5 from Hua et al., 2013). The values provided are classified into five zones, three in the Northern Hemisphere (NH) and two in the Southern Hemisphere (SH), corresponding to different levels of ¹⁴C. For France, the values correspond to the NH zone 2, for the Republic of the Congo to the SH zone 3 and finally for Argentina to the SH zone 1–2. Thus, for our simulations, a yearly value is read for each site.

An F¹⁴C value of 1.8 represents a doubling of the amount of ¹⁴C in atmospheric CO₂. In Fig. 2, it can be noted that the values recorded in France (NH) are higher than those in the Republic of the Congo and Argentina (SH). This is due to the preponderance of atmospheric tests in the NH and the time required to mix air across the Equator.

2.6 Statistical analysis

Simulating carbon processes in soil requires comparison between the model outputs and the measurements to test the model accuracy and possibly implement further improvement. Statistical analysis based on the statistics of deviation were done to evaluate the model–measurement discrepancy according to Kobayashi and Salam (2000) (where a detailed description of the method is provided). Here, we only reproduce the different equations used; x refers to the model outputs and y to the measurements, while i refers to soil depth. The intervals of soil depth of the model outputs and the measurements were homogenized by linearly interpolating the data to common depth intervals defined for each site. The simulations and data were then compared for each depth interval.

$$\begin{array}{}\text{(10)}& \mathrm{RMSD}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n(x_i-y_i{)}^{\mathrm{2}}}\end{array}$$

RMSD is the root mean squared deviation, which represents the mean distance between simulation and measurement.

$$\begin{array}{ll}{\displaystyle }\mathrm{MSD}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n(x_i-y_i{)}^{\mathrm{2}}=(\stackrel{\mathrm{‾}}{x}-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\\ \text{(11)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n{\left[\left(x_i-\stackrel{\mathrm{‾}}{x})-(y_i-\stackrel{\mathrm{‾}}{y}\right)\right]}^{\mathrm{2}}\end{array}$$

MSD, the mean squared deviation, is the square of RMSD. The lower the value of MSD, the closer the simulation results are to the measurements.

$$\begin{array}{}\text{(12)}& \mathrm{SB}=(\stackrel{\mathrm{‾}}{x}-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}},\end{array}$$

where $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{y}$ are the means of x_i (model outputs) and y_i (measurements), respectively.

SB is a part of the MSD (Eq. 14) and represents the bias of the simulation from the measurement.

$$\begin{array}{}\text{(13)}& {\mathrm{SD}}_{\mathrm{s}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n(x_i-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}}\end{array}$$

SD_s is the standard deviation of the simulation.

$$\begin{array}{}\text{(14)}& {\mathrm{SD}}_{\mathrm{m}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n(y_i-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}}\end{array}$$

SD_m is the standard deviation of the measurements.

$$\begin{array}{}\text{(15)}& r={\displaystyle \frac{\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n\left(x_i-\stackrel{\mathrm{‾}}{x})-(y_i-\stackrel{\mathrm{‾}}{y}\right)}{{\mathrm{SD}}_{\mathrm{m}}{\mathrm{SD}}_{\mathrm{s}}}}\end{array}$$

r is the correlation coefficient between the simulation and measurements.

$$\begin{array}{}\text{(16)}& \mathrm{SDSD}=({\mathrm{SD}}_{\mathrm{s}}-{\mathrm{SD}}_{\mathrm{m}}{)}^{\mathrm{2}}\end{array}$$

SDSD is the difference in the magnitude of fluctuation between the simulation and measurements.

$$\begin{array}{}\text{(17)}& \mathrm{LCS}=\mathrm{2}{\mathrm{SD}}_{\mathrm{s}}{\mathrm{SD}}_{\mathrm{m}}(\mathrm{1}-r)\end{array}$$

LSC represents the lack of positive correlation weighted by the standard deviations.

The MSD can therefore be rewritten as

$$\begin{array}{}\text{(18)}& \mathrm{MSD}=\mathrm{SB}+\mathrm{SDSD}+\mathrm{LCS}.\end{array}$$

For the different simulations, the MSD and its components were calculated according to the total soil carbon and to the F¹⁴C.

Table 3 The correlation coefficient (r) between model outputs and measurements for carbon stock (kg C m⁻²) over the soil profile, for the four sites. The results of the initial version of the model ORCHIDEE-SOM-¹⁴C (Control), those from the version including the modification according to He et al. (2016) (He et al., 2016, parameterization) and diffusion varying according to depth (Depth-varying diffusion constant) are provided.

Figure 3 Total soil carbon (kg C m⁻³) according to depth for the four sites. The results of the initial version of the model ORCHIDEE-SOM-¹⁴C (Control), those from the version including the modification according to He et al. (2016) (He et al., 2016, parameterization) and diffusion varying according to depth (Depth-varying diffusion constant) are shown.

3 Model results and evaluation

3.1 Outputs from simulation using the initial version of the model ORCHIDEE-SOM-¹⁴C (Control)

3.1.1 Simulated total soil carbon

Results from the initial version of ORCHIDEE-SOM-¹⁴C show that in all the studied sites, the model succeeds in reproducing the trend of the total carbon profiles, with more carbon at the surface which then decreases according to depth (Fig. 3). Moreover, total soil carbon stock simulated down to 2 m depth is in accordance with data in the case of Misiones and Feucherolles where the major difference mainly lies on the surface. This results in correlation coefficients of 0.44 and 0.2, respectively (Table 3). For the sites of Kissoko and Mons, an overestimation of the total soil carbon is found to a depth of 50 cm for Kissoko and up to a depth of 120 cm for Mons. Correlation coefficients are 0.14 and 0.49 for Kissoko and Mons, respectively (Table 3).

Figure 4 Mean squared deviation (MSD) and its components for total soil carbon (kg C m⁻⁶): lack of correlation weighted by the standard deviation (LCS), squared difference between standard deviations (SDSD) and the squared bias (SB). For the four sites, the results of the initial version of the model ORCHIDEE-SOM-¹⁴C (Control), those from the version including the modification according to He et al. (2016) (He et al., 2016, parameterization) and diffusion varying according to depth (Depth-varying diffusion constant) are shown.

Metrics presented in Fig. 4 showed that this version (ORCHIDEE-SOM-¹⁴C) represents relatively well the observation from Feucherolles (MSD =206 kg C m⁻⁶), whereas the others are highly overestimated (Kissoko, MSD =1343 kg C m⁻⁶; Misiones MSD =2180 kg C m⁻⁶; Mons MSD =3355 kg C m⁻⁶). By detailing the different components of the MSD (Fig. 4), we note that for Mons and Kissoko, standard bias (SB) is the major component of the MSD, contributing 70 % and 60 %, respectively. This reflects that the average of total soil carbon over the soil profile simulated by the model is primarily the origin of the deviation of the model outputs from data. The mean total soil carbon estimated by the model (Table 3) is almost 3 times higher than the mean total carbon measured for Mons (2.37 kg C m⁻² against 0.8 kg C m⁻², respectively) and it is more than 5 times that measured for Kissoko (2.44 kg C m⁻² against 0.42 kg C m⁻², respectively). For Mons a net primary production (NPP) of 6.7 t ha⁻¹ yr⁻¹ was estimated by the technical institute for pasture in this region of France based on the annual yields, whereas the model predicts a NPP of 7.5 t ha⁻¹ yr⁻¹. The large overestimation of the SOC stocks may therefore be due to an overestimation of the NPP. This significant gap recorded in the case of the Kissoko site, where the measured SOC is very low, is probably due to an overestimation of decay rates by ORCHIDEE in sandy soils. The correlation coefficient for Mons is relatively high compared to other sites (Table 3), whereas Fig. 3 shows that the model performance was not very good for this site. This is mainly due to a large SB, whereas other MSD components were rather low.

However, the main components of MSD for Feucherolles and Misiones are both SB (46 % and 56 % for Feucherolles and Misiones, respectively) and also LCS (53 % and 31 % for Feucherolles and Misiones, respectively). This means that for these two sites, the deviation between model outputs and measurements is mainly due to a variation in carbon stock estimation throughout the profile. The mean total soil carbon estimated in both these cases (Table 3) is only slightly higher than those measured (2.03 kg C m⁻² estimated against 2.14 kg C m⁻² measured for Misiones and 0.7  kg C m⁻² estimated against 0.68 kg C m⁻² measured for Feucherolles).

The vertical profiles of the SOC stock were fairly represented by the model. The overestimation, especially at the top, suggests that the distribution of the litter following the root profile and/or the vertical transport of SOC by diffusion are not correctly described in the model.

Figure 5 Modern fraction F¹⁴C according to depth for the four sites. The results of the initial version of the model ORCHIDEE-SOM-¹⁴C (Control), those from the version including the modification according to He et al. (2016) (He et al., 2016, parameterization) and diffusion varying according to depth (Depth-varying diffusion constant) are shown.

3.1.2 Simulated F¹⁴C

Regarding the ¹⁴C activity, bulk F¹⁴C profiles show a classical pattern with higher ¹⁴C activity on the top, slightly influenced by the peak bomb-enriched years. Subsequently, profiles show decreasing ¹⁴C activity with depth (Fig. 5).

The estimated profiles (Control) follow the same trend with a decrease from the surface to depth. However, there is a significant difference between the estimated values and those measured throughout the profile. The statistical analyses (Fig. 6) provide MSD values: 0.02 for Mons and Misiones, 0.03 for Kissoko and 0.09 for Feucherolles. The major component of the MSD in the four sites is the LCS, with a proportion reaching 90% for Mons, 80% for Misiones and 70% for Republic of the Congo, but only 55% for Feucherolles. The high proportions of LCS suggest that the model fails to reproduce the shape of the profile. The lower values estimated by the models reflect a more modern carbon age than in reality. This can be explained, first, by the fact that the root profile puts too much fresh organic carbon in deep soil. Afterwards, in ORCHIDEE, root profile is assumed to follow an exponential function without modulation due to environmental conditions.

Table 4 The correlation coefficient (r) between model outputs and measurements and the mean values (provided by the model and the measurements) over the profile according to F¹⁴C for the four sites. The results of the initial version of the model ORCHIDEE-SOM-¹⁴C (Control), those from the version including the modification according to (He et al., 2016) (He et al., 2016, parameterization) and diffusion varying according to depth (Depth-varying diffusion constant) are provided.

Figure 6 Mean squared deviation (MSD) and its components: lack of correlation weighted by the standard deviation (LCS), squared difference between standard deviations (SDSD) and the squared bias (SB) calculated for modern fraction F¹⁴C. For the four sites, the results of the initial version of the model ORCHIDEE-SOM-¹⁴C (Control), those from the version including the modification according to He et al. (2016) (He et al., 2016, parameterization) and diffusion varying according to depth (Depth-varying diffusion constant) are shown.

Table 5 F¹⁴C profile obtained for each site.

SB's contribution to the MSD does not exceed 7% for Misiones, Kissoko or Mons but reaches about 40% for Feucherolles. This reflects that the mean value of the F¹⁴C estimated by the model and that obtained after the measurements are not very different, except for the Feucherolles site (Table 4). Indeed, the average value estimated for Misiones is 0.920, very close to that measured at 0.930, 0.995 for Kissoko against 0.985 measured and 0.860 for Mons against 0.815 measured. Yet, the difference is greater for the Feucherolles site, the estimated value being 0.915 while the measurement is 0.725. This difference might be caused by the low F¹⁴C value measured at 150 cm (0.257), which the model is not able to capture. This suggests that modelled deep soil carbon is much younger than the observed total soil carbon, probably because ORCHIDEE-SOM simulates a relatively small proportion of passive pool in the lower soil horizons (Fig. 7), while an increasing proportion of passive carbon with soil depth could be expected.

In brief, SOC stocks are generally overestimated and soil carbon age in deep soils (as shown by the F¹⁴C) is underestimated, suggesting that the turnover rate of the passive pool is subject to improvements in ORCHIDEE-SOM.

3.2 Outputs from simulation using the initial version of the model ORCHIDEE-SOM-¹⁴C including the suggestion of He et al. (He et al., 2016, parameterization)

3.2.1 Simulated total soil carbon

Figure 3 shows profile outputs after the suggestion of He et al. (2016) was implemented into ORCHIDEE-SOM-¹⁴C (green dotted curves). Resulting profiles follow the same trend as observations but in this case (labelled “He et al., 2016, parameterization”), the overestimation is very high across the whole profile. This is further confirmed by the metrics analysis (Fig. 4). MSD values markedly increased, resulting in an even higher variance. Obviously, the major component of MSD in all cases is the SB (varying from 80% to 87%) reflecting an even more marked overestimation of the mean total carbon estimates: 7.38 kg C m⁻² against 2.14 kg C m⁻² for Misiones, 2.44 kg C m⁻² against 0.42 kg C m⁻² for Kissoko, 2.33 kg C m⁻² against 0.66 kg C m⁻² for Feucherolles and 9.99 kg C m⁻² against 0.8 kg C m⁻² for Mons.

3.2.2 Simulated F¹⁴C

He et al. (2016) parameterization outputs (Fig. 5, green dotted curves) for F¹⁴C are once again even further away from observations, and MSDs (Fig. 6) are much higher, except for Feucherolles. The MSD components for the Feucherolles site show that the LCS increases from 0.05 to 0.06, whereas the SB decreases from 0.04 to 0.03, again reflecting a variation in the profile more than a difference from the means.

Improvement of the model–measurement fit for the F¹⁴C at 150 cm in Feucherolles confirms that the deep soil carbon simulated by the control version of ORCHIDEE-SOM-¹⁴C was excessively young, since the longer residence time of the passive pool reported by He et al. (2016) resulted in a higher proportion of passive pool across the soil profile (Fig. 7), thus improving deep soil carbon age. Nevertheless, this test only improves the simulation of deep soil carbon in Feucherolles. On the contrary, this increase in carbon residence time increases model deviation from observations for all other cases (Figs. 5 and 6).

Indeed, taking the priming effect into account in this new version of ORCHIDEE has contributed to a 50 % decrease in carbon storage over the historical period. The correction of He et al. (2016) was also aimed at reducing this storage and is of the same order of magnitude as the priming effect. Thus, applying the correction of He et al. (2016) to this version of the model, which takes into account the priming effect, contributes to a double correction for the same target, which then generates this important difference between model outputs and measurements. Moreover, the work of He et al. (2016) is done under the standard parameterization of ORCHIDEE based on “Century”, while ORCHIDEE-SOM was re-parameterized after adding several different processes, the priming effect among them (Camino-Serrano et al., 2018), which makes it difficult to compare results between the two studies.

3.3 Outputs from simulation using the initial version of the model ORCHIDEE-SOM-¹⁴C with diffusion varying according to depth (Depth-varying diffusion constant)

3.3.1 Simulated total soil carbon

Fick's law of diffusion is classically used in models to represent bioturbation assuming that soil fauna activity may be represented following the Fick's law of diffusion (Elzein and Balesdent, 1995; Guenet et al., 2013; Koven et al., 2013; O'Brien and Stout, 1978; Wynn et al., 2005). Using a fixed diffusion constant (D in Eq. 2) implicitly suggests that soil fauna activity is uniform over the entire soil profile. This is generally the case for several models of diffusion, in particular at the level of an ecosystem (Bruun et al., 2007; Guimberteau et al., 2018; O'Brien and Stout, 1978). However, soil faunal activity vary naturally with depth and the diffusion constant should therefore be depth dependent (Jagercikova et al., 2014).

With depth-varying diffusion constant, the carbon profile (orange dashed curves) was improved compared to the initial outputs (Control). The overestimation at the surface decreases at the four sites (Fig. 3). In particular, the Misiones outputs fit very well the observed profiles. This is confirmed with lower MSDs for the four sites for this version compared to the control (Fig. 4).

The total SOC stocks simulated according to this third simulation are closer to the measured values and describing the vertical transport of SOC by varying diffusion according to depth significantly improves the model outputs.

Figure 7 Relative proportion of each of the soil carbon pools summing the total soil carbon at each soil layer. The results of the initial version of the model ORCHIDEE-SOM-¹⁴C (Control, left column), those from the version including the modification according to He et al. (2016) (He et al., 2016, parameterization; middle column) and diffusion varying according to depth (Depth-varying diffusion constant, right column) are shown.

3.3.2 Simulated F¹⁴C

Regarding the F¹⁴C outputs, the simulations using the initial version ORCHIDEE-SOM-¹⁴C in which we assume that the diffusion varies as a function of depth (Depth-varying diffusion constant) results in an improvement of the F¹⁴C profiles (orange dashed curves), in particular for the sites Misiones, Mons and Kissoko (Fig. 5). Statistical analyses prove it with significantly lower MSDs. In addition, the proportion of LCS is 98 %, 92 % and 88 % for Mons, Misiones and Kissoko, respectively, highlighting an estimated average very close to the measurements with a clear disparity, less marked than with the first two simulations, throughout the profile (Fig. 6). Overall, the simulated F¹⁴C to 2 m of depth according to this third simulation are in better agreement with the measured values, and thus incorporating diffusion that varies with depth significantly improves the model outputs.

Using a diffusion coefficient that varies as a function of depth seems to correct the overestimation of the surface total soil carbon by increasing the proportion of labile soil carbon pools in the first soil layers.

When we sum the total soil carbon at each soil layer and look at the relative proportion of each of the soil carbon pools (Fig. 7), we note that it is mainly the distribution of the litter according to depth which varies. In fact, the structural litter proportion is multiplied by about 2 in all four cases, and this proportion remains relatively constant across the profile. This increase in litter proportion has also resulted in a decrease in the passive pool, more pronounced at the surface but also important at depth (except for Feucherolles where the decrease is only marked at the bottom). It suggests that the vertical carbon distribution, which is largely modified by the diffusion coefficient, greatly impacts the SOC and ¹⁴C profiles, which is in line with Dwivedi et al. (2017) who found that the vertical carbon input profiles were important controls over the ¹⁴C depth distribution.

In this study, the vertical transport of SOC and litter through diffusion has been improved by varying diffusion according to depth. Further model development should explore the impact of the other processes defining the soil carbon pools vertical distribution and especially the distribution of the litter according to the root profile.

Overall, by using radiocarbon (¹⁴C) measurements we have been able to diagnose internal model biases (underestimation of deep soil carbon age) and to propose further model improvements (depth-dependent diffusion). Therefore, the use of radiocarbon (¹⁴C) tracers in global models emerges as a promising tool to constrain not only SOC turnover times in the long-term (He et al., 2016) but also internal SOC processes and fluxes that have no direct comparison with field measurements. Nevertheless, the model evaluation performed here on only four sites should be considered as proof of concept and more in-depth evaluations are needed, in particular using a large ¹⁴C database available at a global scale (Balesdent et al., 2018; Mathieu et al., 2015). Indeed, the F¹⁴C is largely controlled by pedo-climatic conditions such as clay content, climate and mineralogy (Mathieu et al., 2015) and the range of situations we covered here is relatively limited.

4 Conclusions

ORCHIDEE-SOM-¹⁴C, is one of the first land surface models (LSMs) that incorporates the ¹⁴C dynamics into the soil (Koven et al., 2013). Its starting point is ORCHIDEE-SOM, a recently developed soil model. We evaluated the new model ORCHIDEE-SOM-¹⁴C for four sites in different biomes. The model almost managed to reproduce the soil organic carbon stocks and the ¹⁴C content along the vertical profiles at all four sites. However, an overestimation of the total carbon stock throughout the profile was noted, with the greatest deviation at the surface. By using radiocarbon (¹⁴C) measurements, we have been able to diagnose internal model biases (underestimation of deep soil carbon age) and to propose further model improvements (depth-dependent diffusion). These results demonstrate the importance of depth-dependent diffusion for improving model outputs with regards to observations. This suggests that, from now on, model improvements should mainly focus on a depth-dependent parameterization. We limited our work here to depth-varying diffusion, but other parameters are also depth dependent and should be represented as such in the next version of the model. For instance, belowground litter production in the model is simply represented by an exponential law without any representation of the effect of resource distribution on root profile (e.g. water or nutrients). This is a complex task in a LSM running at large scale with a classical resolution of 0.5^∘, but the soil modules of LSMs are quite sensitive to the NPP (Camino-Serrano et al., 2018; Todd-Brown et al., 2013), and a better constraint on the profile of the belowground litter production would likely improve the model performance. Furthermore, here we used only one averaged value over the soil profile for soil boundary conditions (texture, pH, bulk density) but those variables are known to impact the F¹⁴C (Mathieu et al., 2015) and change with depth (Barré et al., 2009) and depth-varying boundary conditions may also help to improve the model. Finally, the next step will deal with the comparison of model outputs to data at larger scales to be able to run the new version ORCHIDEE-SOM-¹⁴C at both regional and global scales.