A multi-isotope model for simulating soil organic carbon cycling in eroding landscapes (WATEM_C v1.0)

There is increasing recognition that lateral soil organic carbon (SOC) fluxes due to erosion have imposed an important impact on the global C cycling. Field and experimental studies have been conducted to investigate this topic. It is useful to have a modeling tool that takes into account various soil properties and has flexible resolution and scale options so that it can be widely used to study relevant processes and evaluate the effect of soil erosion on SOC cycling. This study presents a model that is capable of simulating SOC cycling in landscapes that are subjected to erosion. It considers all three C isotopes (¹²C, ¹³C and ¹⁴C) with flexible time steps and a detailed vertical solution of the soil profile. The model also represents radionuclide cycling in soils that can assist in constraining the lateral and vertical redistribution of soil particles within landscapes. The model gives a three-dimensional representation of soil properties including ¹³⁷Cs activity, SOC stock, and δ¹³C and Δ¹⁴C values. Using the same C cycling processes in stable, eroding and depositional areas, our model is able to reproduce the observed spatial and vertical patterns of C contents, δ¹³C values, and Δ¹⁴C values. This indicates that at the field scale with a similar C decomposition rate, physical soil redistribution is the main cause of the spatial variability of these C metrics.

1 Introduction

Soil organic carbon (SOC) is the largest organic C pool on land, with approximately 1550 Pg C in the upper meter of soil (Lal, 2008). This is about 2 times the C in the atmosphere (ca. 760 Pg C). The annual C flux between soil and the atmosphere is ca. 60 Pg C, which makes atmospheric CO₂ sensitive to SOC cycling. SOC stock is a balance between input fluxes and output fluxes, which is controlled by various factors such as soil structure, soil parent material, soil pH, climate, and land use and management. Climate is an important controlling factor on SOC cycling as it is closely related to the rate of both C input and decomposition (e.g., Cox et al., 2000; Davidson and Janssens, 2006). The terrestrial net primary production and SOC decomposition rate generally decrease with increasing temperature (Koven et al., 2017). Globally, SOC stock decreases with increasing temperature (Jobbágy and Jackson, 2000). Land use is another important factor because different vegetation supplies SOC to the soil at different rates (Mahowald et al., 2017; Maia et al., 2010). The stable isotopic composition of SOC is affected by factors such as vegetation type (C₃ vegetation versus C₄ vegetation) and the Suess effect (Tans et al., 1979). Also, SOC can be enriched in ¹³C during the processes of SOC degradation due to preferential mineralization of ¹²C (Natelhoffer and Fry, 1988).

Recent studies have shown that lateral soil redistribution by erosion could also impose an important impact on SOC stock and soil–atmosphere C exchange (Chappell et al., 2016; Doetterl et al., 2016). During erosion events, soil aggregates are broken by raindrops and overland flow, which can enhance SOC decomposition (Van Hemelryck et al., 2011). In the eroding region, SOC in the topsoil is removed by erosion, resulting in depletion of SOC. SOC is added to soil minerals that are moved upwards from below due to soil erosion by inputs from plants (Harden et al., 1999). SOC deposited in depositional settings is buried to depth, and the buried SOC is well preserved (Van Oost et al., 2012). This lateral redistribution of SOC and the consequent disturbance of SOC cycling of both eroding and depositional regions result in spatial variability in SOC stocks and properties. It has been found that eroding sites are depleted of SOC compared to stable sites, while depositional sites are enriched in SOC compared to stable sites in agricultural fields (Li et al., 2007; Van Oost et al., 2005; Yoo et al., 2005). Soil redistribution could lead to differences in SOC stability between eroding and depositional areas. Berhe et al. (2008) found that SOC decomposes at faster rates in eroding areas compared to depositional areas using signatures of ¹⁴C. Wang et al. (2014) reported that SOC mineralization rates in eroding soil profiles are higher than those of depositional soil profiles from results of laboratory soil incubation. Radioactive C isotopes such as ¹⁴C give information on the SOC turnover time, and it is a useful tool to investigate long-term SOC cycling (Trumbore, 2009). SOC redistribution has been found to have an effect on SOC radioactive C isotope composition, with eroding areas more negative compared to depositional areas (Berhe et al., 2008).

Apart from the empirical studies mentioned above, various models have been developed to simulate soil erosion and SOC cycling. At the event scale, there are models simulating processes such as rainfall detachment, sediment entrainment and sediment transport (e.g., Hairsine and Rose, 1992a, b). Some models separate sediments into different grain sizes, and these models are suitable for simulating size selectivity in erosion and deposition (Nearing, 1989; Van Oost et al., 2004). These models are further modified to simulate the selectivity of SOC in erosion and deposition (Wilken et al., 2017). Models based on USLE (Universal Soil Loss Equation) utilize annual mean precipitation as model input to simulate long-term soil erosion (Renard et al., 1997). Given that atmospheric fallouts of ¹³⁷Cs are closely adsorbed on soil particles, ¹³⁷Cs inventories in soils are widely used to evaluate erosion rates (Gaspar et al., 2013; Quine et al., 1997). Soil erosion models have further added processes of ¹³⁷Cs deposition, decay and redistribution associated with soil particles so that they can be calibrated using observed ¹³⁷Cs data (Van Oost et al., 2003).

C turnover models have been developed under the condition of stable landscapes (i.e., free of erosion and deposition) to explore the effects of climate, land use and soil environment on SOC cycling. The decomposition of C is often represented by a first-order kinetic rate. Because SOC is a complex of different components, it is often represented by various pools with respect to C input and decomposition rates in models such as CENTURY (Parton et al., 1987), ICBM (Andren and Katterer, 1997) and RothC (Coleman and Jenkinson, 1995). C fractions obtained in laboratories have been related to C pools in models and have been used to calibrate model parameters to investigate the turnover of various C pools (Skjemstad et al., 2004; Wang et al., 2015a; Zimmermann et al., 2007).

These multiple-pool C models were further integrated with soil erosion models to make them applicable in eroding landscapes. For instance, a study adding erosion processes to the CENTURY model has been used to investigate the balance between the lateral SOC loss by erosion and in situ replacement of lost SOC by photosynthesis in eroding areas, and it has been found that proper management is important to maintain the dynamic replacement of lost C (Harden et al., 1999). At the depositional areas, Wang et al. (2015b) calibrated a profile-scale model integrating erosion and SOC cycling processes using observed SOC content and long-term depositional rate, and it was found that sedimentation rate plays an important role in determining burial efficiency of SOC in colluvial settings. Models have also been developed to investigate the relationships between erosion, crop productivity and SOC cycling (Bouchoms et al., 2019). At the field scale, models that combine SOC redistribution by erosion and SOC dynamics are now able to reproduce the spatial heterogeneity of SOC stock in fields under land uses with eroding areas depleted of SOC and depositional areas enriched in SOC (Liu et al., 2003; Rosenbloom et al., 2001, 2006; Van Oost et al., 2005; Yoo et al., 2005).

Carbon isotopes have also been included in SOC cycling models to constrain model parameters or explore controlling factors. Baisden et al. (2002) used C and N isotopes to simulate the turnover and transport of SOC along soil depth and showed that hydrological conditions had an important role in controlling the vertical transport of SOC. Also, an SOC cycling model integrating C isotope discrimination was utilized to explore the effects of SOC decomposition and physical mixing on the formation of the vertical increase in δ¹³C values with the soil depth (Acton et al., 2013). Ahrens et al. (2014) used ¹⁴C signatures to constrain model parameters of a multi-pool SOC model using the Bayesian method, and the model was further applied to quantify the contribution of sorption, dissolved organic carbon transport and microbial interactions in determining the Δ¹⁴C values of soil profiles (Ahrens et al., 2015). Although field studies have identified the effects of soil redistribution on the profile of SOC isotopes (Berhe et al., 2008), relevant models are not developed yet. SOC models including C isotopes applicable in eroding landscapes would be helpful to fully understand the C isotope profiles as well as the spatial variability of SOC isotopic composition at the landscape scale.

Here, we integrate SOC and soil erosion models and present a model tool that is capable of simulating SOC dynamics in an eroding landscape. The objectives of this modeling tool are the following: (i) it should be a multiple C pool model that is able to represent the complexity of SOC and to be related to the measurable SOC fractions; (ii) it should include various C isotopes so that it could not only represent these C metrics but also use them to constrain the model; and (iii) it should be flexible in terms of spatial and temporal scales so that it would be applicable in various cases regarding spatial and temporal settings.

2 Materials and methods

2.1 Study sites and field data

The first study site is located in the Belgian loam belt. The study area has a temperate climate with average annual precipitation of 750–800 mm and a mean annual temperature of approximately 9.5 ^∘C. Soils in the study area are mainly loess-derived Luvisols with a high silt content (>70 %) and relatively low clay (<15 %) and sand (<20 %) contents (Beuselinck et al., 2000). Arable lands with wheat, maize, sugar beet and potato are the main land use type of the study area. Soil samples were collected from cropland fields with rolling topography, from which soil cores were taken on plateau (with no erosion), convex slope (with erosion) and concave slope (with deposition) areas. SOC contents were measured with a vario MAX CN Macro elemental analyzer (Elementar Analysensysteme GmbH, Germany), while δ¹³C was measured with an ANCA 20-20 GSL mass spectrometer (Sercon Ltd, UK). The inorganic C was removed from soil samples using the HCl fumigation method proposed by Harris et al. (2001).

This study also used published SOC and Δ¹⁴C data in Berhe et al. (2008) to evaluate the developed model. The data were collected at Tennessee Valley in Marin County, northern California (37.9^∘ N latitude, 122.6^∘ W longitude). The climate at Tennessee Valley is Mediterranean, with a mean annual precipitation of 1050 mm and a mean annual temperature of 14 ^∘C. The dominant vegetation cover in the study area is represented by Mediterranean grasses and a coastal shrub (Baccharis pilularis, coyote brush). Soils at Tennessee Valley are derived from chert, greenstone and greywacke sandstone bedrock of the Franciscan assemblage. Soil profiles at the position of plateau (no erosion), convex slope (erosion), concave slope (deposition) and valley (deposition) areas were sampled along a slope. C content was measured with a Carlo Erba elemental analyzer, while the radiocarbon content was measured using accelerator mass spectrometry (AMS) following the methods of Trumbore et al. (1989).

2.2 WATEM_C

Here we present WATEM_C that simulates the redistribution of eroded soil and associated C within the catchment and its effects on the dynamics of SOC. The soil redistribution by water erosion is based on WATEM (Water And Tillage Erosion Model; Van Oost et al., 2000), while the simulation of C dynamics is based on a three-pool C model (Wang et al., 2015a). All three C isotopes (¹²C, ¹³C and ¹⁴C) are included in our model. Soil advection and diffusion through the soil profile are also included in the model. The model uses a flexible time step and vertical resolution of the soil profile so that it can be applied in various settings.

2.2.1 Water erosion

The RUSLE (Revised Universal Soil Loss Equation) (Renard et al., 1997) is used to simulate the long-term potential water erosion (E_pot; kg m⁻² yr⁻¹):

$$\begin{array}{}\text{(1)}& E_{\mathrm{pot}}=R\cdot K\cdot L\cdot S\cdot C\cdot P,\end{array}$$

where R is the rainfall erosivity (MJ mm ha⁻¹ h⁻¹ yr⁻¹), K is soil erodibility (kg h MJ⁻¹ mm⁻¹), L and S are slope length and steep factors, and C and P are factors for cover management and support practices.

The local erosion rate is considered to be equal to the potential erosion rate if the potential erosion rate does not exceed the local transport capacity. The local transport capacity (T_c; kg m⁻¹ yr⁻¹) is calculated as

$$\begin{array}{}\text{(2)}& T_{\mathrm{c}}=k_{\mathrm{tc}}\cdot E_{\mathrm{pot}},\end{array}$$

where k_tc is the transport capacity coefficient (m). In a grid cell, if its sediment inflow exceeds its local transport capacity, the amount of material transported through the grid equals the local transport capacity, while the remainder is deposited in the grid.

A routing algorithm was applied to transfer the mobilized sediments towards the catchment outlet. First, the grids of the study area were sorted in descending order based on the digital elevation model (DEM). Then, after comparing the local transport capacity of a grid cell with the incoming sediment and the locally produced sediment (Van Oost et al., 2000), sediments were routed downslope. Prediction of the flow direction was based on Takken et al. (2001).

The mobilization of SOC or ¹³⁷Cs by erosion (C_ero, kg m⁻² yr⁻¹) is estimated as

$$\begin{array}{}\text{(3)}& C_{\mathrm{ero}}=C_{\mathrm{top}}\cdot R_{\mathrm{ero}}\cdot {\mathrm{ER}}_{\mathrm{ero}},\end{array}$$

where C_top is the content of a carbon isotope or ¹³⁷Cs in the topsoil layer (%), R_ero is the local erosion rate (kg m⁻² yr⁻¹) and ER_ero is the C enrichment ratio in the eroded sediments that equals the ratio of C content in the eroded sediments to that in the source soils.

The deposited SOC or ¹³⁷Cs (C_depo, kg m⁻² yr⁻¹) can be calculated as

$$\begin{array}{}\text{(4)}& C_{\mathrm{depo}}=C_{\mathrm{sed}}\cdot R_{\mathrm{depo}}\cdot {\mathrm{ER}}_{\mathrm{depo}},\end{array}$$

where C_sed is the content of a C isotope or ¹³⁷Cs in the transported sediments (%), R_depo is the local deposition rate (kg m⁻² yr⁻¹) and ER_depo is the C enrichment ratio in the deposited sediments that equals the ratio of the C content in the deposited sediments to that in the bulk transported sediments reaching the depositional sites.

The enrichment ratios of SOC or ¹³⁷Cs in the mobilized sediments at the erosion sites or in the deposited sediments at the deposited sites are found to be closely related to the local erosion or deposition rates (Schiettecatte et al., 2008; Wang et al., 2010). Thus, the enrichment ratios of SOC or ¹³⁷Cs in the mobilized and deposited sediments are calculated as

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathrm{ER}}_{\mathrm{ero}}=a\cdot e^{b\cdot R_{\mathrm{ero}}}+\mathrm{1},\text{(6)}& {\displaystyle }{\mathrm{ER}}_{\mathrm{depo}}=-\mathrm{0.5}{e}^{d\cdot R_{\mathrm{depo}}}+\mathrm{1},\end{array}$$

where a, b and d are coefficients.

2.2.2 Soil C turnover

In our model, the three C isotopes (¹²C, ¹³C and ¹⁴C) are distinguished. As called in the CENTURY model (Parton et al., 1987), each C isotope is divided into three pools that are referred to as active, slow and passive pools. The decomposition of these C pools is described using the following differential equations.

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}{}^{n}A(z,t)}{\mathrm{d}t}}={}^{n}i\left(z\right)-{}^n{k}_{\mathrm{1}}r\left(z\right){}^{n}A(z,t)\text{(8)}& {\displaystyle \frac{\mathrm{d}{}^{n}S(z,t)}{\mathrm{d}t}}=h_{\mathrm{AS}}{}^n{k}_{\mathrm{1}}r\left(z\right){}^{n}A(z,t)-{}^n{k}_{\mathrm{2}}r\left(z\right){}^{n}S(z,t)\text{(9)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{}^{n}P(z,t)}{\mathrm{d}t}}& =h_{\mathrm{AP}}{}^n{k}_{\mathrm{1}}r\left(z\right){}^{n}A(z,t)\\ & +h_{\mathrm{SP}}{}^n{k}_{\mathrm{2}}r\left(z\right){}^{n}S(z,t)-{}^n{k}_{\mathrm{3}}r\left(z\right){}^{n}P(z,t)\end{array}\end{array}$$

Here, ⁿA(z,t), ⁿS(z,t) and ⁿP(z,t) (%) are the contents of the active, slow and passive pools of C isotope n at depth z (m) and time t (year), respectively; ⁿi(z) (Mg C yr⁻¹) is the input of C isotope n at depth z (m); h_AS is the humification coefficient from the active pool to the slow pool, h_AP from the active pool to the passive pool, and h_SP from the slow pool to the passive pool, respectively; r(z) is a coefficient modifying the variation of the C mineralization rate, which denotes the effect of local environmental factors (temperature, humidity, aeration, etc.) at depth z (m); and ⁿk₁, ⁿk₂ and ⁿk₃ (yr⁻¹) are the turnover rates at the reference condition (i.e., r(z)=1) of the active, slow and passive pools of C isotope n, respectively.

¹²C is preferentially lost through microbial respiration compared to ¹³C and ¹⁴C due to its lower atomic weight (Natelhoffer and Fry, 1988). We used a discrimination ratio to denote the difference in mineralization between isotopes, and thus the decomposition rate of a ¹³C pool (¹³k_m, yr⁻¹) can be calculated as

$$\begin{array}{}\text{(10)}& {}^{\mathrm{13}}{k}_m=R_{\mathrm{disc}\mathit{\_}\mathrm{13}}{\cdot }^{\mathrm{12}}{k}_m,\end{array}$$

where R_{disc_13} is the discrimination ratio between ¹³C and ¹²C, and ¹²k_m is the decomposition rate of the corresponding ¹²C pool.

Similarly, the decomposition rate of a ¹⁴C pool (¹⁴k_m, yr⁻¹) can be calculated as

$$\begin{array}{}\text{(11)}& {}^{\mathrm{14}}{k}_m=R_{\mathrm{disc}\mathit{\_}\mathrm{14}}\cdot {}^{\mathrm{12}}{k}_m,\end{array}$$

where R_{disc_14} is the discrimination ratio between ¹⁴C and ¹²C.

The r parameter represents the effect of environmental factors affecting C respiration at a given depth, and it is calculated as

$$\begin{array}{}\text{(12)}& r\left(z\right)=r_{\mathrm{0}}{e}^{-r_{\mathrm{e}}z},\end{array}$$

where r₀ is the value of the r parameter at the topsoil layer, and r_e (m⁻¹) is an exponential decreasing coefficient.

The input of the C isotopes from plant roots decreases exponentially with depth (Gerwitz and Page, 1974; Van Oost et al., 2005):

$$\begin{array}{}\text{(13)}& {}^{n}i\left(z\right)={}^n{i}_{\mathrm{0}}{e}^{-i_{\mathrm{e}}z},\end{array}$$

where ⁿi₀ is the input of C isotope n at the topsoil layer (Mg C yr⁻¹), ⁿi(z) (Mg C yr⁻¹) is the input of C isotope n at depth z (m) and i_e (m⁻¹) is an exponential decreasing coefficient.

The δ¹³C values are expressed in terms of per mill (‰) deviation:

$$\begin{array}{}\text{(14)}& {\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=\left({\displaystyle \frac{({}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}{)}_{\mathrm{Sample}}}{({}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}{)}_{\mathrm{PDB}}}}-\mathrm{1}\right)\cdot \mathrm{1000},\end{array}$$

where (¹³C∕¹²C)_Sample is the abundance ratio of ¹³C to ¹²C of the soil sample, and (¹³C∕¹²C)_PDB is the ratio of the Pee Dee Belemnite (PDB) as the original standard.

Thus, the ¹³C input can be calculated as

$$\begin{array}{}\text{(15)}& {}^{\mathrm{13}}i=({}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}{)}_{\mathrm{PDB}}\cdot \left(\mathrm{1}+{\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\left(t\right)}{\mathrm{1000}}}\right){\cdot }^{\mathrm{12}}i,\end{array}$$

where ¹³i is the ¹³C input (Mg C ha⁻¹ yr⁻¹) from plants, ¹²i is the ¹²C input from plants and δ¹³C(t) is the δ¹³C value of C input at time t (year).

We use the atmospheric Δ¹⁴C record as a proxy for the isotopic ratio of C input via root and leaf litter input (Hua and Barbetti, 2004). In this paper, the following definition of Δ¹⁴C (‰) is used (Stuiver and Polach, 1977):

$$\begin{array}{}\text{(16)}& {\mathrm{\Delta }}^{\mathrm{14}}\mathrm{C}=\left({\displaystyle \frac{({}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}{)}_{\mathrm{Sample}}}{A_{\mathrm{ABS}}}}-\mathrm{1}\right)\cdot \mathrm{1000},\end{array}$$

where (¹⁴C∕¹²C)_Sample denotes the ¹⁴C :¹²C ratio of the sample, and A_ABS denotes the ¹⁴C :¹²C ratio of the standard. A_ABS is set to be $\mathrm{1.176}\times {\mathrm{10}}^{-\mathrm{12}}$ (Karlen et al., 1965; Stuiver, 1980).

The ¹⁴C input can then be calculated as

$$\begin{array}{}\text{(17)}& {}^{\mathrm{14}}i=A_{\mathrm{ABS}}\cdot \left(\mathrm{1}+{\displaystyle \frac{{\mathrm{\Delta }}^{\mathrm{14}}{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ATM}}\left(t\right)}{\mathrm{1000}}}\right){\cdot }^{\mathrm{12}}i,\end{array}$$

where ¹⁴i is ¹⁴C input (Mg C ha⁻¹ yr⁻¹) from plants, and ${\mathrm{\Delta }}^{\mathrm{14}}{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ATM}}\left(t\right)$ is the atmospheric Δ¹⁴C signal at time t (year).

2.2.3 Soil profile evolution due to erosion

In the model, soil profiles are represented as a series of soil layers with equal depths. Given that C input and SOC decomposition rate are related to soil depth, SOC cycling is simulated in each layer independently. Because erosion and deposition change the depth of soil profiles, the model updates the depth of soil profiles and the carbon content of each soil layer every time step. At the eroding locations, soils are removed from the top layer and the soil profile is truncated by the amount of eroded soil. In the meantime, SOC is also lost with the local C enrichment ratio. To keep the soil layer with fixed thickness, soils and associated SOC from soil layers below are incorporated into the upper soil layer at the erosion rate. At the depositional locations, because the top layer is buried by the deposited sediments at the deposition rate, soils and associated SOC are moved downward. For all the soil profiles, the component pools of each C isotope of every layer are updated by homogeneously mixing the component materials every time step.

2.2.4 Advection and diffusion

The vertical transport of mineral and organic components of soil is a complex phenomenon driven by a number of distinct mechanisms such as bioturbation (Jagercikova et al., 2017; Johnson et al., 2014) and chemical mobilization (Taylor et al., 2012). We use the advection–diffusion equation to model vertical transport:

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\mathrm{d}F(z,t)}{\mathrm{d}t}}={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}K\left(z\right)F\left(z,t\right)}{\mathrm{d}{z}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{d}v\left(z\right)F(z,t)}{\mathrm{d}z}},\end{array}$$

where F(z,t) is the concentration of a soil constituent (such as a C isotope pool or ¹³⁷Cs) at depth z (m) at time t (year), v(z) (m yr⁻¹) is the advection term at depth z (m) and K(z) (m yr⁻¹) is a diffusion-type mixing coefficient at depth z (m).

K and v are both depth-dependent and are represented using a sigmoidal scaling function:

$$\begin{array}{}\text{(19)}& v\left(\mathrm{z}\right)={\displaystyle \frac{v_{\mathrm{0}}}{\mathrm{1}+e^{(v_{\mathrm{d}}\cdot \left(z-\mathrm{ct}\right))}}},\end{array}$$

where v_d (m⁻¹) is the depth attenuation of advection, ct (m) is a constant that is set to 0.15 m and v₀ (m yr⁻¹) is the value of v at the topsoil layer.

$$\begin{array}{}\text{(20)}& K\left(z\right)={\displaystyle \frac{K_{\mathrm{0}}}{\mathrm{1}+e^{(K_{\mathrm{d}}\cdot \left(z-\mathrm{ct}\right))}}}\end{array}$$

K_d (m⁻¹) is the depth attenuation of diffusion, and K₀ (m yr⁻¹) is the value of K at the topsoil layer.

2.2.5 ¹³⁷Cs dynamics

¹³⁷Cs is an artificial nuclear radioactive isotope from nuclear tests and reactor incidents. The ¹³⁷Cs in the soil mainly originates from bomb experiments between 1950 and 1970. It falls from the atmosphere primarily in association with precipitation and is rapidly adsorbed to soils by clay materials. As the fallout is well constrained, ¹³⁷Cs has been widely used for tracing the movement of soil and sediment particles in erosion studies (Ritchie and McHenry, 1990). WATEM_C reads the values of the local ¹³⁷Cs fallout. The model then simulates the redistribution of ¹³⁷Cs by soil erosion and deposition at the land surface associated with soil particles (Eqs. 1–6). The model also simulates the downward movement of ¹³⁷Cs in the soil profile by advection and diffusion (Eq. 18). The decay of ¹³⁷Cs (half-life of 30.23 year) is also represented in the model:

$$\begin{array}{}\text{(21)}& C_{\mathrm{s}}\left(z,t+T\right)=C_{\mathrm{s}}\left(z,t\right)\cdot R_{\mathrm{p}}^T,\end{array}$$

where C_s(z,T) (%) is ¹³⁷Cs content at soil depth z (m) at time t (year), R_p (equal to 0.977) is the fraction of ¹³⁷Cs preserved after decay of 1 year and T (year) is the time step of the model iteration.

2.2.6 Model implementation

In order to make the model applicable at various temporal and spatial resolutions, the time step of model iteration and the vertical resolution of the soil profile were not fixed but modifiable as parts of the model input parameters. Given the long-term temporal iteration in SOC cycling processes and the possible large spatial regions where the model may apply, the model was developed using a computation-efficient language (Pascal). The compiled executable file can then be called in other environments such as R (R Development Core Team, 2011), for which the preparation of the input maps is easier. In our model, the default values of the input parameters were given, and at the same time users are allowed to assign custom values to the input parameters in the R environment when calling the executable file. The description and relevant parameters regarding SOC cycling used in this study are listed in Table 1. For the initialization of the ¹³⁷Cs profile, the model first checked if the beginning year of the model simulation was earlier than the beginning of ¹³⁷Cs fallout. If earlier, the initial ¹³⁷Cs profile was set to be zero; if not earlier, the model was run from the beginning of ¹³⁷Cs fallout to the beginning year of the model simulation with relevant parameter values of soil advection and diffusion to generate the initial ¹³⁷Cs profile. For the initialization of C profiles, the initial profile of each C pool was set to be the equilibrium C profile under specific conditions to ensure that the initial C profile was realistic for the study site, i.e., the C profile that made the C input equal to the C mineralization, which was determined by model parameters such as C inputs, C turnover rates and humification coefficients. The model had input parameters of the initial δ¹³C and Δ¹⁴C values of the top and bottom soil layers, and the profiles of δ¹³C and Δ¹⁴C values were then generated by a liner interpolation based on the values of the top and bottom soil layers. The model was then run for a parameter-defined period to get the profiles to reach equilibrium.

Table 1 Values of parameters for SOC cycling used this study.

2.2.7 Model application

Five model scenarios was tested in this study (Table 2). A set of three scenarios was assumed in order to investigate the effect of advection and diffusion as well as lateral soil redistribution by erosion on the spatial and vertical distribution of SOC and δ¹³C and Δ¹⁴C values at the landscape scale. Scenario 1 is a scenario without advection, diffusion or lateral soil redistribution. Scenario 2 is a scenario with vertical advection and diffusion but without lateral soil redistribution. Scenario 3 is a scenario with advection, diffusion and lateral soil redistribution. In order to investigate the effect of plant type change and the Suess effect on the δ¹³C values of soil profiles, the model was applied in another set of three scenarios. Given that advection and diffusion are common in soils, we used the scenario with only advection and diffusion as the reference scenario, i.e., Scenario 2 defined above. The other two scenarios are Scenario 4 with plant type change and Scenario 5 with the Suess effect.

The model was also evaluated using observed C, δ¹³C and Δ¹⁴C soil profiles. Given that an important objective of this study is to investigate the effects of vertical soil advection and diffusion as well as lateral soil redistribution on the profiles of C, δ¹³C and Δ¹⁴C, the model was optimized for parameters of K₀, v₀ and the soil redistribution rate, while the other model parameters were set to be realistic values. For the study site in Belgium, the soil redistribution rates varied between soil profiles due to their locations on the hillslope, while the other parameters were set to be the same among soil profiles. For the study site in the USA, the soil redistribution rates varied between soil profiles. As noted by Berhe et al. (2008), the grass types vary between slope positions. Also, the top C contents showed great differences between profiles (Fig. 1a). Therefore, the C inputs were set to be different between soil profiles. The other parameters were set to be the same among soil profiles for the study site in the USA. The period of erosion was set to be 100 years. The model calibration was performed by comparing the agreement of the simulations and observations, which included both C contents and C isotopic composition (δ¹³C values were available at the Belgian study site, while Δ¹⁴C data were available at the USA study site). A weight factor was introduced to make sure that both the C content and C isotopic composition played equivalent roles in the model calibration.

$$\begin{array}{}\text{(22)}& \mathrm{MRMSE}=\sqrt{\sum _{j=\mathrm{1}}^{C_j}{\left({\displaystyle \frac{{\mathrm{SC}}_j-{\mathrm{OC}}_j}{{\mathrm{SD}}_{\mathrm{C}}}}\right)}^{\mathrm{2}}+\sum _{j=\mathrm{1}}^{I_j}{\left({\displaystyle \frac{{\mathrm{SI}}_j-{\mathrm{OI}}_j}{{\mathrm{SD}}_{\mathrm{I}}}}\right)}^{\mathrm{2}}}\end{array}$$

MRMSE is the modified root mean square error of the model, SC_j (%) is the simulated C content of sample j, OC_j (%) is the observed C content of sample j, C_j is the number of C content observation, SD_C (%) is the standard deviation of the observed C contents of all the samples, SI_j (‰) is the simulated isotopic composition of sample j, OI_j (‰) is the observed C isotopic composition of sample j, I_j is the number of the observed C isotopic composition of all the samples and SD_I (‰) is the standard deviation of the observed C isotopic composition of all the samples.

Figure 1 Observed and simulated C content and Δ¹⁴C values of stable, erosion and depositional areas at the USA study site. Observed and simulated C contents in the format of profiles (a) and 1 : 1 lines (b); observed and simulated Δ¹⁴C values in the format of profiles (c) and 1 : 1 lines (d).

In order to quantify the effects of C decomposition, vertical soil advection and diffusion, and lateral soil redistribution on the C, δ¹³C and Δ¹⁴C profiles, comparisons were performed between the reference profiles (i.e., profiles under the condition of Scenario 1) and profiles under a given condition of SOC decomposition, vertical soil advection and diffusion, and lateral soil redistribution using root mean square error (RMSE):

$$\begin{array}{}\text{(23)}& \mathrm{RMSE}=\sqrt{\sum _{i=\mathrm{1}}^n({\mathrm{SV}}_i-{\mathrm{OV}}_i{)}^{\mathrm{2}}},\end{array}$$

where SV_i is the simulated value of C content or C isotopic composition of soil layer i, OV_i is the observed value of C content or C isotopic composition of soil layer i, and n is the number of soil layers.

The Fourier amplitude sensitivity test (FAST) (Cukier et al., 1973, 1975) was applied using simulations obtained in a Monte Carlo approach to assess the contribution associated with relevant parameters. The FAST method is based on the analysis of variance (ANOVA) decomposition, which quantifies the relative contribution of only one given parameter to the total variance of the model output. Eight parameters of the model relevant to C decomposition (r₀ and r_e), soil advection and diffusion (K₀, K_d, v₀ and v_d), and soil redistribution (soil redistribution rate and erosion time) were tested in 10 000 Monte Carlo scenarios. The values of these parameters were derived from a random distribution within a realistic range (r₀: 0.5–1.5; r_e: 2.6–4 m⁻¹; K₀: 0.005–0.1 m yr⁻¹; K_d: 0.005–0.015 m⁻¹; v₀: 0.01–0.02 m yr⁻¹; v_d: 0.005–0.015; soil redistribution rate: −1 to 1 mm yr⁻¹; erosion time: 1–100 years) for each Monte Carlo scenario. No correlations among these input parameters were assumed for sample generation. The FAST was performed using the MATLAB package developed by Cannavo (2012).

Table 2 Model scenarios implemented in this study.

3 Results

3.1 Model calibration

The optimal parameter values obtained after model calibration are reported in Tables 3 and 4. The model simultaneously simulated both the observed C content and C isotopic composition profiles well, with the MRMSE being 2.26 and 4.46 for the Belgian and the USA study sites, respectively (Figs. 1 and 2). The model not only reproduced the horizontal difference of the C, δ¹³C and Δ¹⁴C profiles between soil profiles well, but the vertical patterns of these profiles were also well represented by the model, except the model underestimated the Δ¹⁴C values at the topsoil layers (Fig. 1c).

Figure 2 Observed and simulated C content and δ¹³C values of stable, erosion and depositional areas at the Belgian study site. Observed and simulated C contents in the format of profiles (a) and 1 : 1 lines (b); observed and simulated δ¹³C values in the format of profiles (c) and 1 : 1 lines (d).

Table 3 Calibrated optimal parameter values for the Belgian study site. K_d and v_d were set to be 0.01 m⁻¹ in the model calibration. The parameter values for SOC cycling are listed in Table 1. SRR indicates the soil redistribution rate.

Table 4 Calibrated optimal parameter values for the USA study site. K_d and v_d were set to be 0.01 m⁻¹ in the model calibration. The parameter values for SOC cycling are listed in Table 1. SRR indicates the soil redistribution rate.

3.2 SOC

Our model is able to reproduce the general pattern of the SOC profile for decreasing SOC content with depth in all scenarios despite rates of advection, diffusion, erosion or deposition (Fig. 3). In Scenario 2, higher rates of soil advection and diffusion result in more SOC transferred to depth, and therefore the difference in SOC content between top layers and bottom layers is smaller under the condition of higher soil advection and diffusion rates compared to SOC profiles with lower advection and diffusion rates (Fig. 3b). In Scenario 3, eroding soil profiles contain less SOC compared to the stable soil profiles free of erosion and deposition, while soil profiles at the depositional area are enriched in SOC compared to the stable soil profile (Fig. 3c).

Figure 3 The simulated C content profiles in (a) Scenario 1, (b) Scenario 2 and (c) Scenario 3. See Sect. 2.2.7 and Table 2 for descriptions of the scenarios. In (b), K₀ (m yr⁻¹) is the diffusion coefficient at the topsoil layer and v₀ (m yr⁻¹) is the advection term at the topsoil layer (Eq. 18). In (c), E indicates the soil redistribution rates, with negative values indicating erosion and positive values indicating deposition. K₀ and v₀ were set to be 0.09 m yr⁻¹ and 0.018 m yr⁻¹, respectively. In (b) and (c), both K_d and v_d (depth attenuation of diffusion and advection) were set to be 0.01 m⁻¹.

3.3 δ¹³C values

In Scenario 1, the δ¹³C profile shows no variation with depth (Fig. 4a). In Scenario 2, the δ¹³C profile decreases with depth (Fig. 4b). The δ¹³C values of the soil profile with higher soil advection and diffusion rates are more negative than that with lower soil advection and diffusion rates (Fig. 4b). In Scenario 3, the δ¹³C values of the eroding profile are less negative than those of the stable soil profile, while soil profiles at the depositional area have more negative δ¹³C values compared to the stable soil profile (Fig. 4c). Our simulation shows that δ¹³C values increase significantly when the vegetation is converted from C₃ vegetation to C₄ vegetation (Fig. 5). When the Suess effect is considered, the δ¹³C values are lower than in scenarios that do not consider the Suess effect (Fig. 5).

Figure 4 The simulated δ¹³C profiles in (a) Scenario 1, (b) Scenario 2 and (c) Scenario 3. See Sect. 2.2.7 and Table 2 for descriptions of the scenarios. In (b), K₀ (m yr⁻¹) is the diffusion coefficient at the topsoil layer and v₀ (m yr⁻¹) is the advection term at the topsoil layer (Eq. 18). In (c), E indicates the soil redistribution rates, with negative values indicating erosion and positive values indicating deposition. K₀ and v₀ were set to be 0.09 and 0.018 m yr⁻¹, respectively. In (b) and (c), both K_d and v_d (depth attenuation of diffusion and advection) were set to be 0.01 m⁻¹.

Figure 5 Effects of plant type change and the Suess effect on the δ¹³C profiles. In the reference C₃ plant scenario, the δ¹³C value of C inputs was set to be −26 ‰; in the Suess effect scenario, the δ¹³C value of C inputs decreased from −26 ‰ to −28.5 ‰ gradually. In the scenario of conversion from C₃ plants to C₄ plants, the δ¹³C value of C inputs was set to be −13 ‰ after vegetation change.

3.4 Δ¹⁴C values

Our model is able to reproduce the general pattern of decreasing Δ¹⁴C values with depth in all scenarios despite rates of soil advection, diffusion, erosion or deposition (Fig. 6). In Scenario 2, soil profiles with higher rates of advection and diffusion have higher Δ¹⁴C values compared to profiles with lower vertical transfer rates (Fig. 6b). In Scenario 3, eroding soil profiles have lower Δ¹⁴C values compared to the stable soil profiles, while soil profiles at the depositional area are enriched in ¹⁴C compared to the stable soil profile (Fig. 6c).

Figure 6 The simulated Δ¹⁴C profiles in (a) Scenario 1, (b) Scenario 2 and (c) Scenario 3. See Sect. 2.2.7 and Table 2 for descriptions of the scenarios. In (b), K₀ (m yr⁻¹) is the diffusion coefficient at the topsoil layer and v₀ (m yr⁻¹) is the advection term at the topsoil layer (Eq. 18). In (c), E indicates the soil redistribution rates, with negative values indicating erosion and positive values indicating deposition. K₀ and v₀ were set to be 0.09 and 0.018 m yr⁻¹, respectively. In (b) and (c), both K_d and v_d (depth attenuation of diffusion and advection) were set to be 0.01 m⁻¹.

3.5 Factors controlling C, δ¹³C and Δ¹⁴C profiles

C decomposition played the primary role in controlling the C, δ¹³C and Δ¹⁴C profiles, with the parameter r₀ accounting for the major variance of the difference between reference profiles and Monte Carlo scenario profiles (Fig. 7). Soil advection and diffusion played a secondary role in controlling the C, δ¹³C and Δ¹⁴C profiles relative to C decomposition, with parameters K₀ and r₀ contributing more to the variance of the difference between reference profiles and Monte Carlo scenario profiles than parameters K_d and v_d. Similarly, soil redistribution played a relatively secondary role compared to C decomposition, with the soil redistribution rate contributing more to the variance of the difference between reference profiles and Monte Carlo scenario profiles than the erosion time. For the C profile, the difference between reference profiles and Monte Carlo scenario profiles was mainly caused by parameter r₀, while for the δ¹³C profile, both r₀ and r_e played important roles. However, for the Δ¹⁴C profiles, r₀, K₀, v₀ and the soil redistribution rate all accounted for more than 15 % of the variance of the difference between reference profiles and Monte Carlo scenario profiles.

Figure 7 The matrix of the proportion of variance of the difference between reference profiles and Monte Carlo scenario profiles caused by model parameters as indicated by the FAST coefficients. SRR indicates the soil redistribution rate, and ET indicates erosion time.

3.6 Spatial variability of soil properties

The model is able to generate a reasonable pattern of soil redistribution, with erosion occurring in upland areas and deposition occurring in footslope areas or valleys (Fig. 8b). Soil redistribution results in higher ¹³⁷Cs inventories in depositional areas than eroding areas (Fig. 8c). The model is also able to generate the spatial variability of SOC stock and properties induced by erosion. The depositional area is enriched in SOC compared with the eroding area (Fig. 8d and e). SOC in the depositional area has lower δ¹³C values (Fig. 8f and g) and higher Δ¹⁴C values (Fig. 8h and i) compared to the eroding area.

Figure 8 Model simulations of erosion and erosion-induced spatial variability of SOC stock and isotopic compositions. (a) DEM (digital elevation model) of the field, (b) erosion and deposition rates (positive values indicate deposition and negative values indicate erosion), (c) ¹³⁷Cs inventory, (d) C stock of topsoil (0–25 cm), (e) C stock of subsoil (26–50 cm), (f) δ¹³C values of topsoil (0–25 cm), (g) δ¹³C values of subsoil (26–50 cm), (h) Δ¹⁴C values of topsoil (0–25 cm), and (i) Δ¹⁴C values of subsoil (26–50 cm).

4 Discussion

In Scenario 1, the shape of the SOC profile is determined by the vertical patterns of SOC input and decomposition rates, both of which decrease with depth. The fact that the basic shape of the SOC profile can be well represented in Scenario 1 shows that the pattern of C input and decomposition rates is the primary controlling factor on the SOC profile, while other factors such as advection, diffusion, erosion or deposition are relatively secondary (Fig. 3a). It is natural that higher rates of advection and diffusion would result in more SOC transferred to deep layers (Fig. 3b). Given that it is less favorable for SOC to be mineralized in deep layers, the transferred SOC by advection and diffusion to depth would be better preserved. Simulations in Scenario 2 show that SOC stock in the top 1 m under the condition of high advection and diffusion rates (K₀=0.09, v₀=0.018) is ca. 14 % higher than that under the condition of low advection and diffusion rates (K₀=0.05, v₀=0.01). Our model can not only reproduce the vertical pattern of SOC distribution in the soil profile, but it can also reproduce the spatial variability of SOC stock due to soil redistribution. The simulations under Scenario 3 are consistent with observations that soil erosion results in spatial variability of SOC stock (Van Oost et al., 2005; VandenBygaart et al., 2012; Yoo et al., 2005). This model assumes the same C input at both the eroding and depositional areas, and therefore in eroding areas, this C input in combination with the decreased heterotrophic respiration rate caused by the decreased SOC stock by erosion result in the replacement of lost SOC at the eroding areas (Harden et al., 1999). Also, SOC buried in the depositional area is partially mineralized over a long period (Van Oost et al., 2012; Wang et al., 2015b). Although offset by the two processes discussed above, soil erosion results in observations that eroding areas are depleted of SOC compared to depositional areas (Van Oost et al., 2005; VandenBygaart et al., 2012; Yoo et al., 2005). In Scenario 1, each soil layer is independent from other soil layers; i.e., there is no mass flux between soil layers due to the neglection of advection, diffusion and soil redistribution. In this case, each soil layer has its C input and decomposition rates, which results in the vertical decrease in both ¹²C and ¹³C with soil depth. The δ¹³C value of each soil layer is therefore determined by the discrimination ratio between ¹³C and ¹²C. If this discrimination ratio is the same between soil layers as implemented in this model, the equilibrium δ¹³C profile would be vertically constant (Fig. 4a). Due to the fact that the condition of no soil advection and diffusion is not realistic, a vertically constant δ¹³C profile with depth is rarely reported. When vertical advection and diffusion are considered as in Scenario 2, the transferred SOC from upper layers is isotopically heavier due to degradation compared to fresh input from plants. This results in an increase in δ¹³C values with soil depth (Fig. 4b). Our simulation shows that vertical soil advection and diffusion can be one of the main causes of the widely observed increase in δ¹³C profiles with depth (Fig. 2). The effect of soil advection and diffusion on the vertical variation of δ¹³C values is more profound in soil profiles with high soil advection and diffusion rates. Because erosion and deposition will truncate or bury the original δ¹³C profiles, this results in the eroding soil profile having higher δ¹³C values compared to the stable soil profile, while the soil profiles at the depositional sites will have lower δ¹³C values in comparison to the stable soil profile (Fig. 4c). This is consistent with the observations made in the croplands in Belgium (Fig. 2). Also, this discrepancy will be more distinct when the erosion or deposition rates become higher. Our simulation shows that soil redistribution by erosion can also cause spatial variability of δ¹³C values on eroding land.

Our model is able to reproduce the widely observed decrease in Δ¹⁴C values with depth in soil profiles (Fig. 4). Also, the model can capture the signal of bomb carbon, with Δ¹⁴C values at the surface layer being positive. In Scenario 1, with no mass fluxes between soil layers, the Δ¹⁴C value is mainly a metric for the turnover rate or residence time of SOC in each layer. The simulated vertical decrease in Δ¹⁴C values is attributed to the vertical variation of environmental conditions that become less favorable for C mineralization. As discussed for δ¹³C profiles, at the same depth the soil profile with low soil advection and diffusion rates contains more degraded and old SOC than the profile with high soil advection and diffusion rates, and therefore the soil profile with low soil advection and diffusion rates has more negative Δ¹⁴C values (Fig. 6b). Similar to δ¹³C profiles, erosion and deposition also have a truncation or burial effect on the Δ¹⁴C profiles, and this results in differences in Δ¹⁴C values between disturbed soil profiles and stable soil profiles at the same soil depth. Therefore, the eroding soil profiles have more negative Δ¹⁴C values compared to the stable soil profiles, while the profiles at the depositional sites have less negative Δ¹⁴C values than the stable soil profiles (Fig. 6c). Our simulation is consistent with observations from an eroding hillslope in northern California by Berhe et al. (2008) (Fig. 1). The causes of more negative Δ¹⁴C values in eroding soil profiles are mainly attributed to the exposure of old SOC from depth, while the observed less negative Δ¹⁴C values in depositional profiles are due to the burial of young SOC from eroding areas.

WATEM_C focuses on the catchment scale, which allows it to account for processes of both erosion and deposition. It is a spatially distributed model with parcel maps denoting various land use types. Also, it allows accounting for soil conservation measurements, which enables the model to investigate anthropogenic effects (such as land use and management) on erosion and SOC cycling. Compared to previous models, the model presented here is more comprehensive. It includes SOC cycling processes and the redistribution of soil and associated SOC by erosion. It is a three-pool C model that discriminates C isotopes (¹²C, ¹³C and ¹⁴C). Thus, it can not only give a three-dimensional representation of C, but also C properties such as δ¹³C and Δ¹⁴C values (Fig. 8). Our model calibration results show that vertical soil advection and diffusion as well as lateral soil redistribution can explain the vertical pattern of C, δ¹³C and Δ¹⁴C profiles as well as their spatial variabilities. FAST shows that C content and C isotopic composition at a given soil depth have different sensitivities to factors such as C decomposition rate, vertical soil advection and diffusion rates, and lateral soil redistribution rates (Fig. 7). The C content is directly related to the C decomposition rate, and thus it is mainly controlled by in situ C decomposition rather than vertical soil advection and diffusion or lateral soil redistribution. The effect of C decomposition on δ¹³C and Δ¹⁴C values is not so dominating as on C content, and vertical and lateral soil redistribution also play important roles in determining the δ¹³C and Δ¹⁴C profiles. The default values of most of the parameters was set in the executable file generated in Pascal, but they can be assigned to custom values before running the executable file in R. This allows the model to be applied in various scenarios of different erosion rates, advection and diffusion rates, or vegetation types by setting relevant parameter values. The model is programmed in a computational efficiency language (Pascal), which makes it suitable to include more C pools and isotopes. Also, the vertical resolution of the soil profile and the temporal resolution of the model iteration is set to be flexible in our model. Users can modify these parameters based on their requirements or circumstances. The arrangement that the model can be called in R makes it easier to prepare various input maps and to proceed with the output of the model. However, it requires users to have experience coding in R. The model is designed to simulate only one period with temporally varying inputs on ¹³⁷Cs fallout as well as ¹³C and ¹⁴C input. For cases of temporal variations such as C input or erosion caused by land use change, the current version of the model is not able to represent these processes.

5 Conclusions

This paper presents a model (WATEM_C) that is capable of simulating SOC dynamics on an eroding landscape. It allows tracking the redistribution of soils and associated ¹³⁷Cs and SOC within the catchment. The model captures the soil profile evolution due to erosion and deposition. The SOC dynamics were simulated using a three-pool C cycling model. All three C isotopes (¹²C, ¹³C and ¹⁴C) are considered in the model and are discriminated with different cycling rates. The model uses a flexible time step and vertical solution of the soil profile. It gives a three-dimensional representation of soil properties such as ¹³⁷Cs activity, SOC stock, δ¹³C values and Δ¹⁴C values. Model calibration shows that the model is able to reproduce the observed spatial pattern of SOC stock: eroding soil profiles are depleted of SOC compared to the stable soil profile, while the depositional soil profile is more enriched in SOC than the stable soil profile. Our simulation is consistent with the observation that the δ¹³C values of the eroding profile are less negative than those of the stable soil profile, while soil profiles at the depositional area have more negative δ¹³C values compared to the stable soil profile. Our model reproduces the observation that eroding soil profiles have lower Δ¹⁴C values compared to stable soil profiles, while soil profiles at the depositional area are enriched in ¹⁴C compared to the stable soil profile. The fact that the spatial patterns of these SOC metrics can be reproduced using the same C cycling processes indicates that physical soil redistribution is the main cause of these spatial variabilities and that our model captures the most important processes and mechanisms in SOC cycling on an eroding landscape. FAST shows that C content is mainly controlled by in situ C decomposition, while δ¹³C and Δ¹⁴C are also to a large extent affected by processes of vertical soil advection and diffusion as well as lateral soil redistribution. We envisage WATEM_C to be a useful tool in simulating the SOC cycling in eroding landscapes, with wide coverage of various soil properties and flexible choices of resolution options and scenario settings.