# The relative importance of decomposition and transport mechanisms in accounting for soil organic carbon profiles

B. Guenet, T. Eglin, N. Vasilyeva, P. Peylin, Philippe Ciais, C. Chenu

► To cite this version:

B. Guenet, T. Eglin, N. Vasilyeva, P. Peylin, Philippe Ciais, et al.. The relative importance of decomposition and transport mechanisms in accounting for soil organic carbon profiles. Biogeosciences, 2013, 10 (4), pp.2379-2392. 10.5194/bg-10-2379-2013. hal-01567150

HAL Id: hal-01567150

https://agroparistech.hal.science/hal-01567150

Submitted on 10 Apr 2018

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

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

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

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

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

# The relative importance of decomposition and transport mechanisms in accounting for soil organic carbon profiles

B. Guenet<fnref n="1" />,<fnref n="2" />, T. Eglin<fnref n="1" />, N. Vasilyeva<fnref n="3" />, P. Peylin<fnref n="1" />, P. Ciais<fnref n="1" />, and C. Chenu<fnref n="3" />

*Correspondence to:* B. Guenet (bertrand.guenet@ua.ac.be)

Received: 5 September 2012 – Published in Biogeosciences Discuss.: 15 October 2012

Revised: 8 March 2013 – Accepted: 14 March 2013 – Published: 10 April 2013

**Abstract.** Soil is the major terrestrial reservoir of carbon and a substantial part of this carbon is stored in deep layers, typically deeper than 50 cm below the surface. Several studies underlined the quantitative importance of this deep soil organic carbon (SOC) pool and models are needed to better understand this stock and its evolution under climate and land-uses changes. In this study, we tested and compared three simple theoretical models of vertical transport for SOC against SOC profiles measurements from a long-term bare fallow experiment carried out by the Central-Chernozem State Natural Biosphere Reserve in the Kursk Region of Russia. The transport schemes tested are diffusion, advection and both diffusion and advection. They are coupled to three different formulations of soil carbon decomposition kinetics. The first formulation is a first order kinetics widely used in global SOC decomposition models; the second one, so-called “priming” model, links SOC decomposition rate to the amount of fresh organic matter, representing the substrate interactions. The last one is also a first order kinetics, but SOC is split into two pools. Field data are from a set of three bare fallow plots where soil received no input during the past 20, 26 and 58 yr, respectively. Parameters of the models were optimised using a Bayesian method. The best results are obtained when SOC decomposition is assumed to be controlled by fresh organic matter (i.e., the priming model). In comparison to the first-order kinetic model, the priming model reduces the overestimation in the deep layers. We also observed that the transport scheme that improved the fit with the data depended on the soil carbon mineralisation formula tion chosen. When soil carbon decomposition was modelled to depend on the fresh organic matter amount, the transport

mechanism which improved best the fit to the SOC profile data was the model representing both advection and diffusion. Interestingly, the older the bare fallow is, the lesser the need for diffusion is, suggesting that stabilised carbon may not be transported within the profile by the same mechanisms than more labile carbon.

## 1 Introduction

Soils are the major reservoir of terrestrial organic carbon (C) representing more than twice the amount of C stored in the atmosphere and three times the amount of C stored in terrestrial vegetation (Schimel, 1995; Schlesinger, 1990; MEA, 2005). In spite of the importance of the stock, the dynamics of soil C is not deeply understood (Sugden et al., 2004). Soil scientists have mainly focused on the surface horizons (Lueken et al., 1962; Sparling et al., 1982; Wu et al., 1993), which were considered to be the only depth of the soil which can emit CO₂ to the atmosphere in significant amounts. However, recent studies have shown that the amount of C stored in the deep layers (below 30 cm) could represent between 30 and 63 % of the total amount of soil C (Batjes, 1996; Jobbagy and Jackson, 2000; Tarnocai et al., 2009). Consequently, increasing attention has been paid to deep soil C and in particular to its dynamics (Fontaine et al., 2007; Salomé et al., 2010; Rumpel et al., 2010; Sanaullah et al., 2010).

Transport mechanisms of soil C into deep layers is still not well understood. The models applied at site-level generally represent both vertical advection and diffusion (Elzein and Balesdent, 1995; Bruun et al., 2007; Braakhekke et al.,

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

<footnote n="2">²Department of Biology, University of Antwerpen, Universiteitsplein 1, 2610, Wilrijk, Belgium</footnote>

<footnote n="3">³AgroParisTech, UPMC-CNRS-AgroParisTech UMR Bioemco7618, 78850 Thiverval-Grignon, France</footnote>
2011), but models also exist with only advection (Feng et al., 1999; Dörr and Münnich, 1989; Jenkinson and Coleman, 2008) or only diffusion (O'Brien and Stout, 1978; Wynn et al., 2005). To our knowledge, no formal and comprehensive comparison of the three transport mechanisms combined to different representations of decomposition has been performed, even if Bruun et al. (2007) suggested that the representation of both advection and diffusion mechanisms improved their model for a sandy soil. However, they compared a model with both advection and diffusion to an advection-only model, but they do not compare these models with a diffusion-only model.

The Soil Organic Matter (SOM) decomposition mechanisms proposed as equations that can be encapsulated in models are also very diverse (for review see, Manzoni and Porporato, 2009; Wutzler and Reichstein, 2008). Within all these approaches, the most used formulation is the first order kinetics as in CENTURY (Parton et al., 1988) or in RothC (Coleman and Jenkinson, 1996). In this formulation, the decay of each SOM pool is proportional to the pool's size, thereby considering that there are no interactions between two decomposing pools.

In particular, within the fourth Assessment Report (AR4) of the Intergovernmental Panel on Climate Change (IPCC), the climate-carbon models used during the Coupled Carbon Cycle Climate Model Intercomparison Project (C4MIP) represented the SOM decomposition with first order kinetics (Friedlingstein et al., 2006). This approach is now criticised (Fontaine and Barot, 2005; Wutzler and Reichstein, 2008), in particular, for its incapacity to represent the relationship existing between Fresh Organic Matter (FOM) inputs (e.g., roots exudates, litter, etc.) and the mineralisation of the SOM. This interaction seems to be a major mechanism of SOM stabilisation in the deep soil layers (Fontaine et al., 2007) even if it could be soil dependent (Salomé et al., 2010). However, here again, to our knowledge no clear comparison between first order kinetics and any of the alternative decomposition formulations linking FOM input to SOM mineralisation has been done.

To study how FOM may possibly interact with SOM mineralisation, experimental sites such as long-term bare fallow soils are interesting experiments. Instead of the complexity of real soils where FOM is permanently added and depends on ecosystem properties, in a bare fallow, the input of FOM has been stopped. Consequently, the relationship between FOM input and SOM mineralisation is switched-off in the bare fallow, whereas it remains switched on in the control plots. In our case, we used bare fallow soils sampled at different years including periods where FOM was assumed to be still available, because it had not yet decomposed, and periods where FOM stock was assumed to be close to zero (Barré et al., 2010).

In this study, we developed a suite of conceptual models to compare the three main transport schemes (advection only $T_A$; diffusion only $T_D$, both together $T_{AD}$) proposed in the

literature using measurement of soil C profiles obtained in a long-term bare fallow and a control plot near Kursk in Russia. We also aimed to cross the three transport schemes with the different formulations used to describe SOM mineralisation: a first order kinetics without relationship between FOM input and SOM mineralisation and a formulation inspired from Wutzler and Reichstein (2008) where the relationship between FOM input and SOM mineralisation is represented. Moreover, we used a three pools model with decomposition based on first order kinetics with both advection and diffusion to better understand how increasing the number of parameters may improve the fit with the data. We first optimised the parameters of each of the seven possible models, using the observed soil carbon profiles and a statistical optimisation method (least square minimization). We then compared the model outputs with measurements from all soil carbon profiles from the bare fallow.

# 2 Materials and methods

## 2.1 The Kursk long-term field experiment data

### 2.1.1 Site and soil plots description

Soils were sampled at the long-term field experiment the Central-Chernozem State Natural Biosphere Reserve named after V. V. Alekhin in the Kursk Region of Russia. The climatic zone is a forest steppe temperate, moderately cold with a mean annual precipitation of 587 mm and a mean annual air temperature of 5.4 °C (Central-Chernozem State Natural Biosphere Reserve, 1947–1997). The soil is a Haplic Chernozem defined as a silty loam Haplic Luvisol following the FAO classification. Two long-term plots were sampled within the site located in the Streletskyi section of the reserve at 51° N, 36° E, about 18 km south of the city of Kursk (Vinogradov, 1984).

The first plot is a long-term bare fallow soil where no fresh input entered the soils since 1947. The soil was weeded manually and tilled every year by horse traction at a depth corresponding to 17–18 cm until the middle of the 1970s and then using a machine at a depth of 22–24 cm. Before the start of the experiment, the soil was under a natural steppe that had been under hay-harvest and pasture for at least the last four centuries. The second plot is geographically close to the first one (about 50 m). It is the same natural steppe that has been absolutely reserved since the establishment of the Reserve in 1935 (Afanasyeva, 1966). It is a natural steppe since 1935 and was used previously for pasture. Dominant plant species are meadow bromegrass (Bromus riparius Rehm.), wild oats (Stipa pennata L.), narrow-leaved meadow grass (Poa angustifolia L.), intermediate wheatgrass (Elytrigia intermedia (Host) Nevski), meadowsweet rose (Filipendula vulgaris Moench), and green strawberry (Fragaria viridis
Duch.). This soil is considered as a control (undisturbed) plot for the bare fallow experiment.

## 2.1.2 Soil sampling and carbon measurements

Soils were sampled at the bare fallow plot in 1967, 1973 and 2005 at depths of every 10 cm down to 150 cm. The soil of the steppe was sampled in 2006 at depths of every 10 cm down to 150 cm. Soils were sampled five times in the steppe, in the bare fallow in 1967 and 2005 and three times in 1973. The corresponding profiles are hereafter called 20YBF, 26YBF and 58YBF for the bare fallow soil sampled in 1967, 1973 and 2005, respectively (i.e., after 20, 26 and 58 yr of bare fallow) and S for the steppe. C contents obtained by the Tyurin method in 1967 and 1973 are corrected by a multiplicative factor 1.13, determined particularly for this soil to match dry combustion method and, thus, avoid any underestimation of the C content (Vasilyeva et al., 2013). In years 2005 and 2006, soil C was measured by dry combustion (Vario Elementar, Analysensysteme, Hanau, Germany).

The bulk density of the soils was first measured for each layer until 120 cm depth in 1959, and this measurement was repeated in 2002 for the bare fallow and for the steppe plots. To take into account the difference in bulk density between the samples, the C stocks are expressed in kg C m⁻² as follows:

$$C_{\text{stock}} = \text{TOC} \times \text{BD} \times h \quad (1)$$

with TOC being the total organic carbon content of the layer considered, expressed in kg C kg⁻¹ soil, BD the bulk density expressed in kg soil m⁻³ and h the layer height expressed in m. Since bulk densities are available only until 120 cm, the soil profiles in kg C m⁻² are calculated only until 120 cm depth. The bulk densities measured in 1959 are used for the soil layer sampled in 1967 and in 1973, whereas the bulk densities measured in 2002 are used for the soil layers sampled in 2005 and in 2006. The chronosequence used in this study is interesting because it follows the decrease of the FOM stock from an equilibrium to another. Indeed, Barré et al. (2010) showed that FOM stock are generally close to zero after 40 yr of bare fallow. Thus, FOM stock is assumed to be high and at the equilibrium for the steppe and close to zero, but also at the equilibrium for the 58YBF. The 20YBF and the 26YBF profiles being intermediate situations.

## 2.2 Soil carbon decomposition models

The models used in this study split the total OM in two pools, the FOM and the SOM for each soil layer (Fig. 1). Input to the FOM pool comes from plant litter and the distribution of this input within the profile is assumed to depend upon depth from the surface (z) according to a negative exponential function, in order to represent the decreased injection of FOM from root mortality with increasing z. This is given by

the equation:

$$I(z) = \frac{I \times \exp(-\mu \times z)}{\sum \exp(-\mu \times z)} \quad (2)$$

where I(z) defines the input of FOM at depth z, I is a scalar corresponding to the total input of carbon. We used I = 3.4 kg C m⁻² yr⁻¹ which is the mean of the C input over the period 1954–1961 for the steppe considered in this study (Afanasyeva, 1966). μ is a parameter describing the exponential of FOM input from root mortality (or exudates) in the vertical profile.

Plant litter production is considered as constant and does not depend on climatic forcing. FOM mineralisation is assumed to be governed by first order kinetics, being proportional to the FOM pool, as given by:

$$\frac{\partial \text{FOM}_{\text{dec}}}{\partial t} = -k_{\text{FOM}} \times \text{FOM} \quad (3)$$

where FOM is the FOM carbon pool and kFOM is a scalar defining its decomposition rate. A fraction of decomposed FOM is humified (e) and another is respired as CO₂ (1 - e) (Fig 1).

$$\text{Humified FOM} = e \times \frac{\partial \text{FOM}_{\text{dec}}}{\partial t} \quad (4)$$

$$\text{Respired FOM} = (1 - e) \times \frac{\partial \text{FOM}_{\text{dec}}}{\partial t} \quad (5)$$

For other fluxes, we test 7 different model formulations, resulting from the combination of three formulations to describe SOM mineralisation (MIN), with three formulations to describe the vertical transport of carbon (T).

### 2.2.1 Formulation MIN1

The mineralisation of SOM is assumed to follow first order kinetics, depending only on the SOM amount, as given by:

$$\frac{\partial \text{SOM}_{\text{dec}}}{\partial t} = -k_{\text{SOM}} \times \text{SOM} \times tf \quad (6)$$

where SOM is the SOM carbon pool, kSOM its decomposition rate, and tf a parameter accelerating the mineralisation due to tillage in the bare fallow soils.

A fraction of decomposed SOM returns to the FOM assuming that a part of the SOM decomposition products is as labile as FOM (e) and another is respired as CO₂ (1 - e) (Fig. 1).

$$\text{Humified SOM} = e \times \frac{\partial \text{SOM}_{\text{dec}}}{\partial t} \quad (7)$$

$$\text{Respired SOM} = (1 - e) \times \frac{\partial \text{SOM}_{\text{dec}}}{\partial t} \quad (8)$$
### 2.2.2 Formulation MIN2

The mineralisation of SOM depends on decomposer activity and, therefore, is sensitive to FOM availability, according to:

$$ \frac{\partial \text{SOM}_{\text{dec}}}{\partial t} = -k'_{\text{SOM}} \times \text{SOM} \times (1 - \exp(-c \times \text{FOM})) \times \text{tf} \quad (9) $$

where $k'_{SOM}$ is a SOM decomposition rate, c a parameter controlling the importance of the FOM carbon pool in the SOM mineralisation, and tf a parameter accelerating the mineralisation due to tillage when necessary. Equation (9) is inspired from Wutzler and Reichstein (2008) who used $\frac{\partial \text{SOM}_{\text{dec}}}{\partial t} = -k'_{SOM} \times \text{SOM} \times (1 - \exp(-c \times \text{MB}))$, but unlike in their study, we do not explicitly represent microbial biomass, the MB term in Wutzler and Reichstein (2008) and rather we use a direct relationship between SOM mineralisation and FOM stock. Thus, Eq. (9) assumes that microbial biomass is always in equilibrium with FOM. A fraction of decomposed SOM returns to the FOM assuming that a part of the SOM decomposition products is as labile as FOM (e) and another is respired as CO₂ (1 - e) following Eqs. (7) and (8).

### 2.2.3 Formulation MIN3

Here, the mineralisation of SOM follows Eq. (6), but SOM is split into two pools with different mineralisation rates: the slow SOM and the passive SOM. The C-FOM mineralised is distributed into these pools according to:

$$ \text{Slow SOM} = f \times \text{Humified FOM} \quad (10) $$

$$ \text{Passive SOM} = (1 - f) \times \text{Humified FOM} \quad (11) $$

where f is the fraction of humified FOM which goes to the slow SOM pool.

A fraction of decomposed passive SOM returns to the FOM assuming that a part of the passive SOM decomposition products is as labile as FOM and another is respired as CO₂ following Eqs. (7) and (8). For the slow pool, a part is respired as CO₂ and the rest is distributed into the FOM assuming that a part of the slow SOM decomposition products is as labile as FOM C mineralised and the rest is going to the passive pools.

The three different formulations of C transport are based on an advection equation ($T_A$) or on Fick's law of diffusion ($T_D$) or on both transport mechanisms as defined by Elzein and Balesdent (1995) ($T_{AD}$).

### 2.2.4 Formulation $T_A$

Advection is defined by:

$$ F_A = A \times C \quad (12) $$

where $F_A$ is the flux of C transported downwards by advection, A the advection rate (mm yr⁻¹).

<figure><img src="image_4.png" /><figcaption>Fig. 1. Schematic representation of the fluxes and the pool in the model.</figcaption></figure>

### 2.2.5 Formulation $T_D$

The Fick's law is defined by:

$$ F_D = -D \times \frac{\partial^2 C}{\partial^2 z} \quad (13) $$

where $F_D$ is the flux of C transported downwards by diffusion, $-D$ the diffusion coefficient (cm² yr⁻¹) and C the amount of carbon in the pool subject to transport (FOM or SOM).

### 2.2.6 Formulation $T_{AD}$

Finally, advection and diffusion are represented following the scheme of Elzein and Balesdent (1995)

$$ F_{AD} = F_A + F_D \quad (14) $$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Presentation of the seven formulations used in this study.</figcaption></figure>

$$F_{AD} = A \times C - D \times \frac{\partial^2 C}{\partial z^2} \qquad (15)$$

where $F_{AD}$ is the flux of C transported downwards by advection and diffusion, A the advection rate, $-D$ the diffusion coefficient and C the amount of carbon in the pool subject to transport. We finally build the seven different models by forming pair of (MIN, T) formulations as illustrated by Fig. 2. Thus, the FOM and SOM pools dynamics correspond to:

$$\frac{\partial \text{FOM}}{\partial t \partial z} = I + \frac{\partial F_X}{\partial z} + e \times \frac{\partial \text{SOM}_{\text{dec}}}{\partial t} - \frac{\partial \text{FOM}_{\text{dec}}}{\partial t} \qquad (16)$$

$$\frac{\partial \text{SOM}}{\partial t \partial z} = \frac{\partial F_X}{\partial z} + e \times \frac{\partial \text{FOM}_{\text{dec}}}{\partial t} - \frac{\partial \text{SOM}_{\text{dec}}}{\partial t} \qquad (17)$$

All the models were developed using R 2.11.1 and run at a yearly time step. The models run from the ground (0 cm) until 200 cm and the vertical resolution is 5 mm for each layer. To compare with the data, C stocks are then summed each 20 layers to obtain 10 cm layers. Moreover only the layers until 120 cm are used since the data could be converted in kg C m⁻² only until 120 cm. The equations were solved using the deSolve library (Soetaert et al., 2010). This library solves a system of ordinary differential equations resulting from one dimensional partial differential equations that have been converted to ordinary differential equations by numerical differencing. We run the models with the steppe condition (i.e., with input of FOM) during 2000 yr to reach the equilibrium. The steppe was assumed to be at the equilibrium and

we used the C stocks obtained after the 2000 yr run for the steppe. Then we run the model for 58 yr without FOM input to reproduce the 58YBF plots condition. From this run, we extracted the data after 20 and 26 yr to reproduce the 20YBF and 26YBF plots condition, respectively.

## 2.3 Parameter optimisation

The 10 parameters used for each simulation are listed in Table 1. Eight of them are optimised for each model using a Bayesian inversion method with priors (see Tarantola, 1987) against the entire dataset (48 data points, i.e., 12 for each profile), with a statistical approach based on a Bayesian framework (Tarantola, 1987). Our approach is based on Santaren et al. (2007). The optimal parameter set corresponds to the minimum of the cost function

$$J(x) = \frac{1}{2} \left[ (y - H(x))^T R^{-1} (y - H(x)) + (x - x_b)^T P_b^{-1} (x - x_b) \right] \quad (18)$$

and contains both the mismatch between modelled and observed fluxes and the mismatch between prior and optimised parameters. $x$ is the vector of unknown parameters, $x_b$ the prior parameters, $H(x)$ the model outputs and $y$ the vector of observations. $P_b$ describes the prior parameter error variances/covariances, while $R$ contains the prior data error variances/covariances.

An efficient gradient-based iterative algorithm, called L-BFGS-B (Zhu et al., 1995) was used to minimise the cost function. This algorithm prescribes a range of values for each parameter. At each iteration, the gradient of the cost function $J(x)$ is computed, with respect to all the parameters. The L-BFGS-B algorithm does not provide uncertainties or
<table><caption>Table 1. Model parameters summary for the model MIN1-T<sub>A</sub>, MIN1-T<sub>D</sub>, MIN1-T<sub>AD</sub>, MIN2-T<sub>A</sub>, MIN2-T<sub>D</sub>, MIN2-T<sub>AD</sub>.</caption><thead><tr><th rowspan="2">Model parameter</th><th rowspan="2">Meaning</th><th rowspan="2">Optimized parameter</th><th rowspan="2">Prior range</th><th colspan="6">Posterior modes ± variance (prior modes)</th></tr><tr><th>MIN1-T<sub>A</sub></th><th>MIN1-T<sub>D</sub></th><th>MIN1-T<sub>AD</sub></th><th>MIN2-T<sub>A</sub></th><th>MIN2-T<sub>AD</sub></th><th>MIN2-T<sub>AD</sub></th></tr></thead><tbody><tr><td>I</td><td>Input of FOM (kg C m<sup>-2</sup> yr<sup>-1</sup>)</td><td>NO</td><td>N.A</td><td></td><td></td><td></td><td>3.4</td><td></td></tr><tr><td>μ</td><td>Exponential rate parameter of FOM input from root mortality</td><td>YES</td><td>0.2–2</td><td>1.382 ± 0.015<br>(0.357)</td><td>1.0 ± 0.714<br>(1.0)</td><td>1.466 ± 0.004<br>(0.228)</td><td>1.521 ± 0.004<br>(1.742)</td><td>0.954 ± 0.11<br>(0.357)</td><td>1.459 ± 0.002<br>(0.369)</td></tr><tr><td>k<sub>FOM</sub></td><td>Decomposition rate of FOM</td><td>YES</td><td>0.2–10</td><td>0.467 ± 0.163<br>(0.467)</td><td>0.833 ± 0.192<br>(0.833)</td><td>0.200 ± 0.059<br>(0.221)</td><td>0.313 ± 0.139<br>(0.313)</td><td>0.467 ± 0.163<br>(0.467)</td><td>0.375 ± 0.15<br>(0.375)</td></tr><tr><td>r</td><td>Respiration rate of FOM</td><td>NO</td><td>N.A</td><td></td><td></td><td></td><td>0.4</td><td></td><td></td></tr><tr><td>k<sub>SOM</sub></td><td>Decomposition rate of SOM</td><td>YES</td><td>2e<sup>-3</sup>–0.2</td><td>1.11e<sup>-2</sup> ± 0.1<br>e<sup>-4</sup>(2.50e<sup>-2</sup>)</td><td>5.46e<sup>-3</sup> ± 0.17<br>e<sup>-3</sup>(3.33e<sup>-2</sup>)</td><td>3.36e<sup>-3</sup> ± 0.45<br>e<sup>-3</sup>(3.52e<sup>-3</sup>)</td><td>6.74e<sup>-3</sup> ± 0.12<br>e<sup>-3</sup>(4.15<sup>-3</sup>)</td><td>1.04e<sup>-2</sup> ± 0.04<br>(2.50e<sup>-2</sup>)</td><td>1.02e<sup>-2</sup> ± 0.01<br>(3.65e<sup>-2</sup>)</td></tr><tr><td>e</td><td>Humification rate</td><td>YES</td><td>0.1–0.9</td><td>0.726 ± 0.132<br>(0.726)</td><td>0.5 ± 0.06 (0.5)</td><td>0.323 ± 0.027<br>(0.324)</td><td>0.579 ± 0.084<br>(0.579)</td><td>0.726 ± 0.132<br>(0.726)</td><td>0.669 ± 0.110<br>(0.669)</td></tr><tr><td>A</td><td>Advection coefficient (mm yr<sup>-1</sup>)</td><td>YES</td><td>0.01–1</td><td>0.282 ± 0.001<br>(0.298)</td><td>N.A</td><td>0.365 ± 0.335<br>(0.366)</td><td>0.361 ± 0.326<br>(0.361)</td><td>N.A</td><td>0.578 ± 0.31<br>(0.388)</td></tr><tr><td>D</td><td>Diffusion coefficient (cm<sup>2</sup> yr<sup>-1</sup>)</td><td>YES</td><td>1–100</td><td>N.A</td><td>66.61 ± 62.3<br>(50.0)</td><td>13.63 ± 4.65<br>(13.63)</td><td>N.A</td><td>36.32 ± 22.07<br>(55.74)</td><td>14.86 ± 5.53<br>(14.86)</td></tr><tr><td>C</td><td>Influence of the FOM carbon pool in the SOM mineralisation (priming parameter)</td><td>YES</td><td>0.1–160</td><td>N.A</td><td>N.A</td><td>N.A</td><td>64.91 ± 1046.31<br>(64.69)</td><td>47.57 ± 20.31<br>(12.93)</td><td>90.0 ± 1959.5<br>(88.5)</td></tr><tr><td>tf</td><td>Acceleration of SOM mineralisation due to tillage</td><td>YES</td><td>1–2.5</td><td>1.97 ± 0.97<br>(1.97)</td><td>1.5 ± 0.56 (1.5)</td><td>1.02 ± 0.26<br>(1.02)</td><td>2.23 ± 1.25<br>(2.23)</td><td>1.97 ± 0.96<br>(1.97)</td><td>1.01 ± 0.27<br>(1.01)</td></tr></tbody></table>

error correlations between optimised parameters, but, when the $J(x)$ is minimised, it calculates the posterior error covariance matrix on the parameters $P_a$ from the prior error covariance matrices and the Jacobian of the model at the minimum of the cost function, using the linearity assumption (Tartola, 1987). Absolute values of the error correlations close to 1 imply that the observations do not provide independent information to differentiate a couple of parameters.

The optimisation assumes that the errors associated to the model parameters and the observations can be described with Gaussian Probability Distribution Function (PDF). Assuming Gaussian PDF, this approach is sensitive to potential local minima. We, therefore, performed 20 optimisations starting from different prior parameter values randomly distributed in their allowed range of variation. We then selected the case that provides the lowest cost function. With this approach, we are much less sensitive to potential local minima. It makes use of prior information on the parameters, minimising an objective function that measures the distance between modelled and observed carbon vertical profiles and between prior and optimised parameter values (using a least squares approach).

The optimised parameters for the models using MIN1 and MIN2 (Table 1) are the SOM decomposition rate ($k_{som}$), the FOM decomposition rate ($k_{fom}$), the exponential rate parameter of FOM input from root mortality in the vertical profile ($\mu$), the fraction of SOM mineralised recycled in FOM ($e$), Fick's coefficient in models using $T_A$ or $T_{AD}$ ($A$), the advection rate in models using $T_D$ of $T_{AD}$ ($D$), the parameter controlling the FOM dependency of the SOM mineralisation

in models using MIN2 ($c$), and the tillage factor to increase mineralisation ($tf$). Prior estimates for each parameter are given on Table 1. We used such ranges to define priors because they are close to parameters already published (Baisden et al., 2002; Bruun et al., 2007; Braakhekke et al., 2011). However, since the $c$ parameter has never been estimated before, we considered the prior as non-informative and we set a very large prior error (50%). For the model using MIN3 the optimised parameters are presented in Table 2. They are the slow SOM decomposition rate ($k_{slow}$), the passive SOM decomposition rate ($k_{passive}$), the FOM decomposition rate ($k_{fom}$), the exponential rate parameter of FOM input from root mortality in the vertical profile ($\mu$), the fraction of SOM mineralised recycled in the other pools ($e$), the fraction of FOM mineralised distributed into the slow SOM pool ($f_{FS}$), the fraction of the slow SOM pool mineralised not respired and recycled in the FOM pool ($f_{SF}$) the Fick's coefficient ($D$), the advection rate ($A$).

As for the variance of the model-data residuals, note that with our formalism it should include both the measurement error and the model error. Given that the error on the measurements could be estimated from the existence of several replicates for each profile, we choose the measured standard deviation as error on the observations. At its minimum, $J(x)$ should be close to the half of the number of observations (reduced chi-square of one). Note finally that all errors (on the parameters and the observations) are assumed to be uncorrelated.

A compaction of the bare fallow soil is observed on the site, leading to 10 cm difference between the deepest studied
<table><caption>**Table 2.** Model parameters summary for the model MIN3-TAD.</caption><thead><tr><th>Model parameter</th><th>Meaning</th><th>Optimised parameter</th><th>Prior range</th><th>Posterior modes ± variance (prior modes)</th></tr></thead><tbody><tr><td>I</td><td>Input of FOM (kg C m<sup>-2</sup> yr<sup>-1</sup>)</td><td>NO</td><td>N.A</td><td>3.4</td></tr><tr><td>μ</td><td>Exponential rate parameter of FOM input from root mortality</td><td>YES</td><td>0.2–2</td><td>1.68 ± 0.02 (0.86)</td></tr><tr><td>k<sub>FOM</sub></td><td>Decomposition rate of FOM</td><td>YES</td><td>0.2–10</td><td>0.246 ± 0.137 (0.102)</td></tr><tr><td>r</td><td>Respiration rate of FOM</td><td>NO</td><td>N.A</td><td>0.4</td></tr><tr><td>k<sub>Slow SOM</sub></td><td>Decomposition rate of slow SOM pool</td><td>YES</td><td>2e<sup>-3</sup>–0.2</td><td>8.7e<sup>-3</sup> ± 8.6e<sup>-3</sup> (8.34e<sup>-5</sup>)</td></tr><tr><td>k<sub>Passive SOM</sub></td><td>Decomposition rate of passive SOM pool</td><td>YES</td><td>3.3e<sup>-4</sup>–1</td><td>1.7e<sup>-3</sup> ± 0.6e<sup>-3</sup> (2.53e<sup>-4</sup>)</td></tr><tr><td>r</td><td>Fraction of SOM C respired</td><td>YES</td><td>0.1–0.9</td><td>0.742 ± 0.017 (0.725)</td></tr><tr><td>A</td><td>Advection coefficient (mm yr<sup>-1</sup>)</td><td>YES</td><td>0.01–1</td><td>0.686 ± 0.402 (0.818)</td></tr><tr><td>D</td><td>Diffusion coefficient (cm<sup>2</sup> yr<sup>-1</sup>)</td><td>YES</td><td>1–100</td><td>18.11 ± 1.22 (57.83)</td></tr><tr><td>f<sub>FOM→Slow</sub></td><td>Fraction of humified FOM that goes to the slow SOM pool</td><td>YES</td><td>0.1–0.9</td><td>0.347 ± 0.030 (0.377)</td></tr><tr><td>F<sub>Slow→FOM</sub></td><td>Fraction of decomposed SOM from the slow pool that goes to the FOM pool*</td><td>YES</td><td>0.1–0.9</td><td>0.415 ± 0.043 (0.459)</td></tr><tr><td>tf</td><td>Acceleration of SOM mineralisation due to tillage</td><td>YES</td><td>1–2.5</td><td>2.5 ± 0.16 (2.1)</td></tr></tbody></table>

horizon of the control and the one of the bare fallow. To take into account graphically the compaction effect on soil depth, we defined the point at 0 m depth as the floor of the steppe and then the soil layers were assumed to be linearly compacted through time since 1947 to reproduce the observed final difference of 10 cm depth between the two bottom horizons. The compaction observed on site and its effects on transport are taken into account through the use of the bulk density in the stocks calculation in the dataset. The compaction effects are implicitly represented in the model as the optimisation was performed with the stocks expressed in kg C m⁻².

## 2.4 Comparison of model results with data

C stock measured and modelled for each soil layer are compared using statistical indicators developed first by Kobayashi and Salam (2000) and then improved by Gauch et al. (2003). These statistical indicators are the Mean Squared Deviation (MSD), the Squared Bias (SB), the Non-Unity slope (NU) and the Lack of Correlation (LC).

$$ \text{MSD} = \frac{\Sigma(m-o)^2}{n} \qquad (19) $$

with *o* the observed values, *m* the C stock calculated by the model and *n* the number of observations.

The MSD is decomposed into three additive components following Gauch et al. (2003): the Squared Bias (SB), the Non-Unity slope (NU) and the Lack of Correlation (LC).

$$ \text{SB} = (\bar{m} - \bar{o})^2 \qquad (20) $$

$$ \text{NU} = \left( 1 - \frac{\Sigma (m - \bar{m}) \times (o - \bar{o})}{\Sigma (m - \bar{m})^2} \right)^2 \times \frac{\Sigma (m - \bar{m})^2}{n} \qquad (21) $$

$$ \text{LC} = \left( 1 - \frac{\Sigma ((m - \bar{m}) \times (o - \bar{o}))^2}{\Sigma (o - \bar{o})^2 \times \Sigma (m - \bar{m})^2} \right) \times \frac{\Sigma (o - \bar{o})^2}{n} \qquad (22) $$

SB provides information about the mean bias of the simulation from the measurement. NU indicates the capacities of the model to reproduce the magnitude of fluctuation among the measurements. It can be considered as an indicator of the model's capacities to reproduce the scattering of the data. LC is an indication of the dispersion of the point over a scatterplot, i.e., the capacities of the model to reproduce the shape of the data.

Even though the previous statistical indicators are useful to compare one model with the data, they must be not used to compare the models together because they do not consider the differences in the parameters number. In order to do so, we used the Bayesian Information Criterion (BIC).

$$ \text{BIC} = \log(\text{MSD}) \times n + \log(n) \times p \qquad (23) $$

where MSD is the mean squared deviation calculated with Eq. (16), *n* the number of samples (i.e., 12 in our case, one per profile), and *p* the number of parameters.

## 3 Results

The optimised parameters and their associated variances are presented in the Tables 1 and 2. For some models, important correlation factors were observed (supplemental material Fig. 1). Considering the method used to optimise the parameters, these important correlation factors make complicated the presentation of the model output within an envelope. Therefore, we presented the model results using the optimised parameter without any envelope (Figs. 4, 5 and 6). The most important correlations were generally observed between μ and kSOM which control the input and the output of

<footnote>* Assuming that a part of the decomposition products from Slow SOM pools are as labile as FOM.</footnote>
<table><caption>**Table 3.** Bayesian Information Criterion (BIC).</caption><thead><tr><th></th><th>MIN1-T<sub>A</sub></th><th>MIN1-T<sub>D</sub></th><th>MIN1-T<sub>AD</sub></th><th>MIN2-T<sub>A</sub></th><th>MIN2-T<sub>D</sub></th><th>MIN2-T<sub>AD</sub></th><th>MIN3-T<sub>AD</sub></th></tr></thead><tbody><tr><td>Entire dataset</td><td>181.4</td><td>179.3</td><td>173.6</td><td>196.2</td><td>182.1</td><td>177.7</td><td>200.5</td></tr><tr><td>Steppe</td><td>63.3</td><td>39.2</td><td>57.0</td><td>67.2</td><td>59.5</td><td>61.6</td><td>70.7</td></tr><tr><td>20 yr bare fallow</td><td>41.2</td><td>51.7</td><td>46.3</td><td>54.1</td><td>45.8</td><td>40.8</td><td>56.5</td></tr><tr><td>26 yr bare fallow</td><td>48.4</td><td>56.5</td><td>50.6</td><td>57.6</td><td>50.0</td><td>47.4</td><td>53.3</td></tr><tr><td>58 yr bare fallow</td><td>54.1</td><td>59.8</td><td>58.4</td><td>54.6</td><td>62.3</td><td>63.3</td><td>70.8</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Components of mean squared deviation (MSD) for the seven formulations. The lowest the MSD value is, the best the fit is. The three components are lack of correlation (LC), non-unity slope (NU), and squared bias (SB).</figcaption></figure>

the SOM pools, but also between *c* and *k*<sub>SOM</sub> for MIN2-T<sub>D</sub> which both control the SOM mineralisation.

## 3.1 The representation of SOM decomposition

Figure 3 describes the MSD, SB, NU and LC statistical indicator values obtained by each different model for the entire dataset, or for each site. When MSD was calculated for the entire dataset, we observed generally very closed values of MSD between MIN1 and MIN2 except for T<sub>A</sub> formulation where the use of a first order kinetic to describe SOM mineralisation improve the fit with the data. With BIC that takes into account the effect of the parameter increase on the fit, the same pattern was observed (Table 3). The MIN3-T<sub>AD</sub> model presents a high MSD value. Even though, MSD values

are close for T<sub>D</sub> and T<sub>AD</sub> formulation, the values of the MSD components are different. MIN2 (priming model) decrease the SB values indicating that this formulation better represents the mean C stock over the profile, whereas MIN1 decrease the NU values, suggesting that first order kinetics better represent the data scattering.

For the youngest bare fallow plots, we always observed the same patterns. MSD values are reduced using priming model (MIN2) when soil C transport is represented using advection (T<sub>A</sub>) and both advection and diffusion (T<sub>AD</sub>). The BIC is following the same pattern (Table 3). In these cases, the improvement is mainly due to an important reduction of the SB values indicating the priming model (MIN2) better represent the mean C stock over the profile. Indeed, the Figs. 4 and 5 show that the MIN1 formulation overestimates the decomposition. However, when transport was represented using only advection (T<sub>A</sub>) the priming model (MIN2) over estimated the decomposition in the surface layers. The MIN3 formulation (three C pools) never presented the lowest MSD values when both advection and diffusion were represented. However, the SB value was reduced when SOM decomposition was described using MIN3.

For the oldest bare fallow, the priming model (MIN2) reduced MSD values when C transport was described using advection only (T<sub>A</sub>) or diffusion only (T<sub>D</sub>). When C transport mechanism was only advection, priming model (MIN2) reduced the LC value but increased the SB value indicating that MIN2 better reproduced the shape of the data, but poorly represented the mean C stock over the profile. For T<sub>D</sub> formulation, the priming model (MIN2) better represented the mean C stock over the profile (reduced SB values) and the data scattering (reduced NU values). However, the lowest BIC values were always obtained with MIN1 whatever the C transport mechanisms used (Table 3). The highest MSD value was obtained with the MIN3 formulation for the oldest bare fallow, but LC value was largely reduced indicating that such formulation better represented the shape of the data.

For the steppe profile, the priming model (MIN2) increased the MSD and the BIC values for the three transport mechanisms. In particular, when only diffusion was represented (T<sub>D</sub>) the MIN2 formulation overestimated the C stock (Figs. 4 and 5). The MIN3 formulation never improved the description of the steppe profile.
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Total Organic Carbon for the MIN1 transport scheme (first order kinetics) and the three transport schemes ($T_A$, $T_D$, $T_{AD}$). Measured C stocks are represented by the circles and modelled C stocks by the lines.</figcaption></figure>

## 3.2 The transport formulation

Considering the entire dataset, for each SOM mineralisation formulation, the best fit was always obtained with the same transport formulation ($T_{AD}$). This was also the case for the 58YBF, but for this profile, the best fit was always obtained with advection only. It is interesting to note that for the steppe and the youngest bare fallows (20YBF and 26YBF), the transport mechanisms that reduce the MSD and the BIC values depended on the SOM mineralisation formulation used. When using the MIN2 substrates interactions representation, the best fits were obtained with the formalisms including advection and diffusion (Fig. 3). In this case, for the steppe, the SB values were particularly reduced. However, the lowest BIC values were obtained when only diffusion is represented. For the 20YBF, it was the LC value that was reduced, suggesting that the shape of the curves is better represented with advection and diffusion. When using MIN1, the best fit for the steppe was obtained using only diffusion. The scattering of the data was better represented in this case (reduced NU values). For the 20YBF, the best fit was obtained using advection only ($T_A$). Here, the SB values were largely reduced. It is interesting to note that the mechanisms transport that produced the best fit with the data changed depending

on the FOM availability. Indeed, when SOM mineralisation was represented by MIN1, the best fit was obtained with diffusion ($T_D$) for the steppe and then with advection ($T_A$) for the bare fallows. With the priming model (MIN2), the lowest BIC values were always obtained with only diffusion for the steppe, then with both advection and diffusion ($T_{AD}$) for 20YBF and 26YBF. Finally, for the 58YBF only advection was needed to fit the data.

# 4 Discussion

## 4.1 The representation of SOM decomposition

Our goal was to better separate the role of vertical transport mechanisms such as diffusion/advection given different formulation of SOM decomposition using a simple conceptual model of SOM decomposition. For the data we used, we first showed that the substrates interactions representation proposed by Wutlzer and Reichstein (2008) was an interesting formulation to represent cases where FOM is not at the equilibrium such as in the young bare fallow soil (20YBF and 26YBF). For the entire set of cases, the priming model reduced the standard bias and, therefore, better reproduced the
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Total Organic Carbon for the MIN2 transport scheme (priming model) and the three transport schemes ($T_A$, $T_D$, $T_{AD}$). Measured C stocks are represented by the circles and modelled C stocks by the lines.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Total Organic Carbon for the MIN3 transport scheme (three pools model) and the $T_{AD}$ transport schemes. Measured C stocks are represented by the circles and modelled C stocks by the lines.</figcaption></figure>

mean C stock over the profile. For the youngest bare fallow, the priming formulation better performed when $T_D$ or $T_{AD}$ transport formulation are used. For the steppe, where the FOM stock is at equilibrium, the MIN1 better performed for the three transport mechanisms. However, the priming model improved the model output for the deep layers (> 40 cm) where FOM amount was low (Fig. 6). For the youngest bare fallow (20YBF), the profile representation was improved for

the priming model when advection and both advection and diffusion were represented. Indeed, when SOM decomposition was described by first order kinetics the model overestimated the SOM decomposition in deep soil layers. This result is in agreement with the study of Dörr and Münnich (1989). It suggests that the SOM decomposer activity is largely controlled by the availability of FOM at the Kursk site, and that using a single decomposition rate for all the soil layers
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Fresh organic matter for the four profiles calculated by the model. The steppe, the 20YBF, the 26YBF and the 58YBF are represented by the green line, the blue line, the black line and the red line, respectively. All the models share the same input scheme. We assume that in the models a fraction of the SOM decomposed may be as labile as the FOM and is, therefore, incorporated to the FOM pool.</figcaption></figure>

parameterised on the top layers might lead to largely over-
estimate SOM decomposition at depth (and underestimate it
near the surface). This conclusion is in agreement with the
results of Fontaine et al. (2007) who observed an important
increase of SOM mineralisation at depth when FOM was
added in an in vitro experiment. Finally, the MIN3 formula-
tion never obtained the lowest BIC, suggesting that the best
fit is not always obtained with the most parameterised model.
It shows that when the priming model improves the fit with
the data, it is not just an effect of increasing the number of
parameters.

4.2 Transport mechanisms

We found that for the MIN1 formulation the worst fit to the
data was always observed when only diffusion was repre-
sented, except for the steppe. In particular, the LC values
were higher indicating that the shape of the data curves was
not well represented when using only diffusion. When the
MIN2 (priming model) was used, the worst fit was always
obtained with advection except for 58YBF. In this case, the
LC values are high, indicating that using advection only was
not sufficient to reproduce the shape of the profile. The ad-
vection rates obtained after optimisation in this study were
ten times higher than those presented in Bruun et al. (2007),
but ten times lower than those presented in Braakhekke et
al. (2011). For the diffusion coefficient, the values obtained

here after optimisation are also higher than those of Bruun
et al. (2007) (one or two range of magnitude), but this diffu-
sion coefficient is a function of the bulk density (Braakhekke
et al., 2011). Thus, differences in the bulk density between
the soil used here and the one used by Bruun et al. (2007)
might explain the different diffusion coefficients. Baisden et
al. (2002) obtained good agreement between data and model
output with a model which only uses advection as transport
mechanism. The advection rates in our study are in agree-
ment with those observed by Baisden et al. (2002).

When substrates interactions are included in any of the
conceptual models, we observed that only diffusion must
be represented to properly fit the data for the steppe where
FOM is available. Then, for 20YBF and 26YBF profiles, a
representation of diffusion and advection was needed to fit
well with the data. Finally, for the 58YBF profile, only ad-
vection was needed. A shift from diffusion to advection, as
the most important mechanism to fit the data, was also ob-
served when SOM mineralisation was described by first or-
der kinetics. However, in this case advection was already the
most important mechanism for 20YBF and 26YBF soils. In
long-term bare fallows older than 40 yr, most of the labile C
has been mineralised (Barré et al., 2010). Consequently, the
SOM in this soil is quite different in decomposability from
that in the youngest bare fallow plots and from the steppe. For
example, particulate organic matter, i.e., decomposing plant
residues, which are labile components of SOM are depleted
from a temperate bare fallow in a few decades (Vasilyeva et
al., 2013). Diffusion is often used to account for transport
of plant debris and particulate organic matter by soil fauna,
whereas advection is used to represent C transport with the
liquid phase (O'Brien and Stout, 1978; Wynn et al., 2005;
Braakhekke et al., 2011). Soil fauna activity is closely re-
lated to SOM availability (Decaëns, 2010). Therefore, the
importance of soil fauna in the transport C in our sites de-
creased when FOM input were stopped and, therefore, when
SOM availability was reduced. This suggests that different
pools of SOM could be transported through different mech-
anisms. The more labile OM may be transported mainly by
bioturbation, whereas the more stabilised may be transported
with the liquid phase. To our knowledge, this is the first time
that different transport mechanisms are suggested for differ-
ent pools of C. This assumption must be tested against other
soil profiles from bare fallow experiments, but if confirmed it
suggests that soil models using different pools of C and aim-
ing to represent the C distribution within a profile must use
different transport mechanisms for labile and stable pools.

4.3 Transport mechanisms depending on the SOM
decomposition formulation

We observed that the transport mechanism inducing the
best fit for youngest bare fallow (20YBF and 26YBF)
where FOM stock are out from equilibrium was not always
the same for each decomposition formulation and might,
therefore, depend on the formulation used to described SOM
decomposition. Indeed, when SOM decomposition was de-
scribed with a first order kinetic, the best fit could only be ob-
tained with advection as the sole transport mechanism. How-
ever, a better fit was obtained with both advection and diffu-
sion as transport mechanisms when SOM decomposition de-
pended on the FOM stock. This suggests that the FOM input
regime may determine the most important transport mecha-
nism. Moreover, the structural importance of transport mech-
anisms in a model depends on how the SOM decomposition
is formulated. As a consequence, the use of such models to
understand and separate mechanisms not directly observable
may be a highly complex task, whose results could depend
on the underlying assumption in the formulation of the SOM
decomposition. For example, a first order kinetic model such
as the one used in this study assumes that the microbial com-
munity responsible for SOM decomposition is stable in terms
of biomass, but also in terms of structure and physiology dur-
ing the period considered. Several observations showed that
microbial community structure, biomass and physiology are
controlled by environmental conditions such as soil moisture
(e.g., Williams, 2007; Guenet et al., 2011) or temperature
(e.g., Pettersson and Bååth, 2003; Wu et al., 2009). More-
over, Hirsch et al. (2009) showed that the microbial commu-
nity structure differ between grassland and bare fallow soils.
The absence of an explicit representation of the microbial
community or biomass might explain also why our models
do not fit so well with the data when first order kinetics are
used. The second formulation obtained generally better fits,
but was not able to reproduce all the profiles perfectly. The
latter assumes that there is a constant nutrient limitation on
the microbial activity, which is implicitly represented in the
parameters of the model. Mikhailova et al. (2000) showed
that the N profiles in the Kursk site also differ between bare
fallow and the control. Thus, an explicit representation of the
N cycle in the profile might decrease the MSD values. Fi-
nally, the effects of temperature and soil moisture are not
represented in the models because not enough climate data
was available. The absence of such effects may explained at
least partially why the models and, in particular, the most
complex do not perfectly fit with the data after optimisation.

To our knowledge, this is the first time that the impor-
tance of SOM decomposition formulation on the transport
formulation is analysed and showed and it could have im-
portant consequences on the representation of C transport in
the model. Indeed, developing models inclusive of more and
more mechanisms should likely improve our capacities to re-
produce large scales datasets and improve our understand-
ing of the C cycle in the soil. Nevertheless, one must keep
in mind that these mechanisms will interact with each other
within the model structure and the choice of a certain repre-
sentation of different mechanisms will depend on how others
mechanisms are represented in the model.

4.4 Does a model perform better?

Regarding the MSD values for the entire dataset, we may
conclude that the better model over the six models tested
would be the MIN1-TAD and the MIN2-TAD which presented
very close values of MSD. However, the MIN2-TAD had one
more parameter and BIC must be used instead of MSD to
take into account the differences in the parameter number.
Using BIC MIN1-TAD would be the best model. Neverthe-
less, when each profile is analysed independently, MIN2-TAD
obtained the lowest BIC values for half of the profile. The an-
swer may also depend on the objectives fixed for the study.
Using the Gauch et al. (2003) evaluation methods, we evalu-
ated different characteristics of the models such as their ca-
capacities to reproduce the mean C stock value over the profile,
the standard deviation around this mean value and the shape
of the profile. For example, if the objective is to represent
the data scattering of the steppe profile very well, we should
use the model MIN1-TD. But if the objective is to evaluate
the mean C stock of the steppe, we should use the model
MIN2-TAD. Furthermore, the LC values were generally the
highest contributor to MSD. It suggests that all models were
not very good to precisely reproduce the shape of the data.
Finally, the correlation matrix (Supplement Fig. 1) showed
some cases with important correlations between parameters,
making it difficult for the calculation of an envelope. How-
ever, generally, correlation factors were low indicating that
our model was not over-parameterised in spite of the limited
amount of data compared to the number of optimised param-
eters.

Supplementary material related to this article is
available online at: http://www.biogeosciences.net/10/
2379/2013 bg-10-2379-2013-supplement.pdf.

*Acknowledgements.* The authors acknowledge the reserve staff for maintaining such an experimental site for so long and for the access to the data. The authors thank the GIS-Climat Environnement Société program CARBOSOIL and R2DS for financial support.

Edited by: E. Veldkamp

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

INSU
Institut national des sciences de l'Univers

The publication of this article is financed by CNRS-INSU.
## References

Afanasyeva, E. A.: Chernozemy sredne-russkoi vozvishennosti, Nauka, Moscow, 224 pp., 1966 (in Russian).

Baisden, W. T.: A multiisotope C and N modeling analysis of soil organic matter turnover and transport as a function of soil depth in a California annual grassland soil chronosequence, Global Biogeochem. Cy., 16, 1135, doi:10.1029/2001GB001823, 2002.

Barré, P., Eglin, T., Christensen, B. T., Ciais, P., Houot, S., Kätterer, T., van Oort, F., Peylin, P., Poulton, P. R., Romanenkov, V., and Chenu, C.: Quantifying and isolating stable soil organic carbon using long-term bare fallow experiments, Biogeosciences, 7, 3839–3850, doi:10.5194/bg-7-3839-2010, 2010.

Batjes, N. H.: Total carbon and nitrogen in the soils of the world, Eur. J. Soil Sci., 47, 151–163, 1996.

Braakhekke, M. C., Beer, C., Hoosbeek, M. R., Reichstein, M., Kruijst, B., Schrumpf, M., and Kabat, P.: SOMPROF: A vertically explicit soil organic matter model, Ecol. Model., 222, 1712–1730, 2011.

Bruun, S., Christensen, B. T., Thomsen, I. K., Jensen, E. S., and Jensen, L. S.: Modeling vertical movement of organic matter in a soil incubated for 41 years with ¹⁴C labeled straw, Soil Biol. Biochem., 39, 368–371, 2007.

Coleman, K., Jenkinson, D. S., Crocker, G. J., Grace, P. R., Klir, J., Körschens, M., Poulton, P. R., and Richter, D. D.: Simulating trends in soil organic carbon in long-term experiments using RothC-26.3., Geoderma, 81, 29–44, 1997.

Decaëns, T.: Macroecological patterns in soil communities, Glob. Ecol. Biogeogr., 19, 287–302, 2010.

Dörr, H. and Münnich, K.: Downward movement of soil organic matter and its influence on trace element transport, Radiocarbon, 31, 655–663, 1989.

Elzein, A. and Balesdent, J.: Mechanistic Simulation of Vertical-Distribution of Carbon Concentrations and Residence Times in Soils, Soil Sci. Soc. Am. J., 59, 1328–1335, 1995.

Feng, X., Peterson, J. C., Quideau, S. A., Virginia, R. A., Graham, R. C., Sonder, L. J., and Chadwick, O. A.: Distribution, accumulation, and fluxes of soil carbon in four monoculture lysimeters at San Dimas Experimental Forest, California, Geochim. Cosmochim. Ac., 63, 1319–1333, 1999.

Fontaine, S. and Barot, S.: Size and functional diversity of microbe populations control plant persistence and long-term soil carbon accumulation, Ecol. Lett., 8, 1075–1087, 2005.

Fontaine, S., Barot, S., Barré, P., Bdioui, N., Mary, B., and Rumpel, C.: Stability of organic carbon in deep soil layers controlled by fresh carbon supply, Nature, 450, 277–280, 2007.

Friedlingstein, P., Cox, P., Betts, R., Bopp, L., Von Bloh, W., Brovkin, V., Cadule, P., Doney, S., Eby, M., Fung, I., Bala, G., John, J., Jones, C., Joos, F., Kato, T., Kawamiya, M., Knorr, W., Lindsay, K., Matthews, H. D., Raddatz, T., Rayner, P., Reick, C., Roeckner, E., Schnitzler, K. G., Schnur, R., Strassmann, K., Weaver, A. J., Yoshikawa, C., and Zeng, N.: Climate-carbon cycle feedback analysis: Results from the C4MIP model intercomparison, J. Climate 19, 3337–3353, 2006.

Gauch Jr., H. G., Hwang, J. T., and Fick, G. W.: Model Evaluation by Comparison of Model-Based Predictions and Measured Values, Agron. J., 95, 1442–1446, 2003.

Guenet, B., Lenhart, K., Leloup, J., Giusti-Miller, S., Pouteau, V., Mora, P., Nunan, N., and Abbadie, L.: The impact of long-term CO₂ enrichment and moisture levels on soil microbial commu-

nity structure and enzyme activities, Geoderma, 170, 331–336, 2011.

Hirsch, P. R., Gilliam, L. M., Sohi, S. P., Williams, J. K., Clark, I. M., and Murray, P. J.: Starving the soil of plant input for 50 years reduces abundance but not diversity of soil bacterial communities, Soil Biol. Biochem., 41, 2021–2024, 2009.

Jenkinson, D. S. and Coleman, K.: The turnover of organic carbon in subsoils. Part 2. Modelling carbon turnover, Eur. J. Soil Sci., 59, 400–413, 2008.

Jobbagy, E. G. and Jackson, R. B.: The vertical distribution of soil organic carbon and its relation to climate and vegetation, Ecol. Appl., 10, 423–436, 2000.

Kobayashi, K. S. and Salam, M. U.: Comparing simulated and measured values using mean squared deviation and its components. Agron. J., 92, 345–352, 2000.

Lueken H., Hutcheon, W. L., and Paul, E. A.: The influence of nitrogen on the decomposition of crop residues in soil, Can. J. of Soil Sci., 42, 276–288, 1962.

Manzoni, S. and Porporato, A.: Soil carbon and nitrogen mineralisation: Theory and models across scales, Soil Biol. Biochem., 41, 1355–1379, 2009.

MEA: Millennium Ecosystem Assessment-Nutrient Cycling, World Resource Institute, Washington DC, 2005.

Mikhailova, E. A., Vassenev, R. B., Schwager, I. I., and Post, S. J.: Cultivation effects on soil carbon and nitrogen contents at depth in the Russian Chernozem, Soil Sci. Soc. Am. J., 64, 738–745, 2000.

O'Brien, B. J. and Stout, J. D.: Movement and turnover of soil organic matter as indicated by carbon isotope measurements, Soil Biol. Biochem., 10, 309–317, 1977.

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

Pettersson, M. and Bååth, E.: Temperature-dependent changes in the soil bacterial community in limed and unlimed soil, FEMS Microbiol. Ecol., 45, 13–21, 2003.

R Development Core Team: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria, ISBN 3-900051-07-0, http://www.R-project.org., 2010.

Rumpel, C. and Kögel-Knabner, I.: Deep soil organic matter – a key but poorly understood component of terrestrial C cycle, Plant Soil, 338, 143–158, 2010.

Salomé, C., Nunan, N., Pouteau, V., Lerch, T. Z., and Chenu, C.: Carbon dynamics in topsoil and in subsoil may be controlled by different regulatory mechanisms, Glob. Change Biol., 16, 416–426, 2010.

Sanaullah, M., Chabbi, A., Leifeld, J., Bardoux, G., Billou, D., and Rumpel, C.: Decomposition and stabilization of root litter in top- and subsoil horizons: what is the difference?, Plant Soil, 338, 127–141, 2010.

Santaren, D., Peylin, P., Viovy, N., and Ciais, P.: Optimizing a process-based ecosystem model with eddy-covariance flux measurements: A pine forest in southern France, Global Biogeochem. Cy., 21, GB2013, doi:10.1029/2006GB002834, 2007.

Schimel, D. S.: Terrestrial ecosystems and the carbon cycle, Glob. Change Biol., 1, 77–91, 1995.

Schlesinger, H. W.: Evidence from chronosequence studies for a low carbon-storage potential of soils, Nature, 348, 232–234, 1990.
Soetaert, K., Petzoldt, T., and Setzer, R. W.: Solving Differential Equations in R: Package deSolve, J. Stat. Softw., 33, 1–25, 2010.

Sparling, G. S., Cheschire, M. V., and Mundie, C. M.: Effect of barley plants on the decomposition of 14C-labelled soil organic matter, J. Soil Sci., 33, 89–100, 1982.

Sugden, A., Stone, R., and Ash, C.: Ecology in the underworld – Introduction, Science, 304, 1613–1613, 2004.

Tarantola, A.: Inverse Problem Theory: Methods of Data Fitting and Model Parameter Estimation, Elsevier Science Ltd., 630 pp., 1987.

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

Vasilyeva, N. A., Chenu, C., Tyugai, Z. N., and Milanovskiy, E. Y.: Century scale C metastability in full Chernozem profiles under changed organic matter input and tillage, in press, 2013.

Vinogradov, B. V.: Aerospace studies of protected natural areas in the USSR, in: Conservation, science and society, Nat. Resour. Res., XXI, Vol. 2, 435–448, UNESCO-UNEP, Paris, 1984.

Williams, M.: Response of microbial communities to water stress in irrigated and drought-prone tallgrass prairie soils, Soil Biol. Biochem., 39, 2750–2757, 2007.

Wu, J., Brookes, P. C., and Jenkinson, D. S.: Formation and destruction of microbial biomass during the decomposition of glucose and Ryegrass in Soil, Soil Biol. Biochem., 25, 1435–1441, 1993.

Wu, Y., Xiongsheng, Y., Haizhen, W., Na, D., and Jianming, X.: Does history matter? Temperature effects on soil microbial biomass and community structure based on the phospholipid fatty acid (PLFA) analysis, J. Soil Sediment, 10, 223–230, 2009.

Wutzler, T. and Reichstein, M.: Colimitation of decomposition by substrate and decomposers – a comparison of model formulations, Biogeosciences, 5, 749–759, doi:10.5194/bg-5-749-2008, 2008.

Wynn, J. G., Bird, M. I., and Wong, V. N. L.: Rayleigh distillation and the depth profile of 13C/12C ratios of soil organic carbon from soils of disparate texture in Iron Range National Park, Far North Queensland, Australia, Geochim. Cosmochim. Ac., 69, 1961–1973, 2005.

Zhu, C., Byrd, R. H., Lu, P., and Nocedal, J.: A limited memory algorithm for bound constrained optimisation, SIAM J. Sci. Stat. Comput., 16, 1190–1208, 1995.