Non-Standard Discrete RothC Models for Soil Carbon Dynamics

Soil Organic Carbon (SOC) is one of the key indicators of land degradation. SOC positively affects soil functions with regard to habitats, biological diversity and soil fertility; therefore, a reduction in the SOC stock of soil results in degradation, and it may also have potential negative effects on soil-derived ecosystem services. Dynamical models, such as the Rothamsted Carbon (RothC) model, may predict the long-term behaviour of soil carbon content and may suggest optimal land use patterns suitable for the achievement of land degradation neutrality as measured in terms of the SOC indicator. In this paper, we compared continuous and discrete versions of the RothC model, especially to achieve long-term solutions. The original discrete formulation of the RothC model was then compared with a novel non-standard integrator that represents an alternative to the exponential Rosenbrock–Euler approach in the literature.


Introduction
The United Nations Convention to Combat Desertification (UNCCD) is an international agreement, established in 1994, that links the environment and development with sustainable land management. The first objective indicated in the UNCCD 2018-2030 Strategic Framework is to improve the conditions of affected ecosystems, combat desertification/land degradation and promote sustainable land management [1]. For each country, the commitment is to achieve no net loss of land-based natural capital by 2030 [2]. No net loss means that the quantity and quality of land-based natural capital are maintained or increased, despite the impacts of global environmental change, whether due to human or natural causes. Land degradation is monitored through the changes of the values of a specific set of consistently measured indicators from their baseline quantities, conventionally identified as their initial values. The deviations from the baseline values of these indicators are the basis for monitoring land degradation.
Soil Organic Carbon (SOC) is one key indicator of land degradation [3]. Monitoring the SOC stocks and the loss of soil organic carbon due to land use changes is fundamental for maintaining the physical, chemical and biological quality of soil [4]. Soil organic carbon positively affects soil functions with regard to habitats, biological diversity, soil fertility, crop production potential, erosion control and water retention. A high SOC content improves the processes of soil formation, nutrient storage, water holding capacity and the absorption of organic or inorganic pollutants. Thus, a reduction of the SOC stock not only indicates soil degradation, but may also have potential negative effects on soil-derived ecosystem services.
Starting from the SOC baseline, predictive spatial modelling can simulate the carbon dynamics, estimate carbon sequestration under the actual land use and evaluate the deviation from the baseline average value of total carbon [5]. Moreover, a dynamical model may determine the optimal potential land use pattern that is suitable to achieve land degradation neutrality in terms of SOC indicator values. Well-validated models, such as Rothamsted Carbon (RothC) [6], CENTURY [7] and MOMOS [8], which take into account the interactions among the climate, pedology, cropping systems and soil and crop management, can be used to predict SOC changes under different management practices and climatic conditions. These models are essentially compartmental, meaning that they represent soil organic matter as a few discrete compartments (generally two to five) characterized by different chemical characteristics of the soil's degradation. The decomposition rates, applied to each compartment, are governed by kinetic and stoichiometric laws and are mainly ruled by the environmental conditions (e.g., soil moisture level, aeration and soil temperature). These models are used in a variety of ways and often for long-term studies. Indeed, being able to compute and predict long-term solutions is extremely valuable for various reasons: it gives a synthetic view of the system in the given agro-climatic conditions; it makes it possible to test if a studied soil has reached equilibrium or not; and it allows envisioning what would be the consequences of specific events on a given soil [9]. In this paper, we compared continuous and discrete versions of the RothC model, especially regarding long-term solutions. Moreover, since the discrete RothC version can be interpreted as a first-order approximation of the continuous model, we introduced a non-standard discrete approximation that can be interpreted as a novel discrete model of soil carbon dynamics. The original discrete formulation of the RothC model was then compared to the novel non-standard integrator, which represents a different approach with respect to the Exponential Rosenbrock-Euler (ERE) discretization [10].

Rothamsted Carbon Model-Continuous Formulation
The SOC indicator is considered to be the result of the equilibrium between the inputs and outputs of the soil system. SOC contained in the organic matter is constantly built up and decomposed and is then released into the atmosphere as CO 2 and recaptured through photosynthesis. Inputs in the soil organic matter decomposition model consist of two major components: living organisms' biomass (mainly plant roots and microbial biomass) and plant and animal residues at various stages of decomposition. Outputs result from the heterotrophic respiration processes when soil organic carbon is used as an energy source by soil organisms and returned to the atmosphere as CO 2 fluxes. The Rothamsted Carbon model [6,11] (RothC) is a model of carbon turnover in non-waterlogged soils [12]. Initially developed for arable soils, it was later expanded to grasslands and forests. It takes into account the effects of temperature, moisture content and soil type. The RothC model divides the soil carbon into four active compartments and one inactive, characterized by different chemical decomposition rates of degradability (see Figure 1). The Inert Organic Matter (IOM) represents the inactive pool, resistant to decomposition, which does not receive carbon (C) inputs. At each time step, incoming plant residues are split between easily Decomposable Plant Material (DPM) and Resistant Plant Material (RPM), depending on the ratio γ 1 − γ , which estimates the decomposability of the particular plant material inputs, which in turn depends on the specific cultivation being considered. The fraction 2 η of the input of Farmyard Manure (FYM), if any, is equally split between the DPM and RPM compartments; the remaining part 1 − 2 η enters in the system directly as Humified organic matter (HUM). Both DPM and RPM decompose to form CO 2 , microbial Biomass (BIO) and more HUM. The fraction α + β of metabolised C incorporated into the sum of compartments BIO+HUM is determined by the clay content of the soil, while the remaining part 1 − α − β is released as CO 2 and lost by the system. The BIO+HUM carbon content is then split into α α + β percent BIO and β α + β percent HUM. Finally, both BIO and HUM decompose to form more CO 2 , BIO and HUM. The four active compartments undergo decomposition as a function of different rate constants, which correspond to the entries of the vector k = [k dpm , k rpm , k bio , k hum ] and of the rate modifier ρ(t), which depends on the clay content of the soil, on climatic variables (rainfall, temperature, open pan evaporation) and land cover.
In real soil systems, processes involved in the RothC model are continuous in time, and thus, in [10], the author proposed the following continuous formulation: where c(t) = [c dpm (t), c rpm (t), c bio (t), c hum (t)] T and c 0 ≥ 0 denotes the vector of the initial concentrations. The matrix A is given by: The vector b(t) represents the carbon amount entering the system at time t. It considers both the input of plant residues g(t) a (g) and the input of FYM f (t) a ( f ) , so that: The entries of vectors a (g) := [γ, 1 − γ, 0, 0] T and a ( f ) := [η, η, 0, 1 − 2 η] T are the fraction inputs 0 ≤ γ ≤ 1, 0 ≤ η ≤ 1/2, which sum up to one; the carbon of plant residues enters the soil only through the DPM and RPM compartments, while the carbon amount of FYM enters through the DPM, RPM and HUM compartments.
It can be shown that, under the hypothesis 0 < α + β ≤ 1 and ρ(t) > 0, the solution of the homogeneous part of the RothC system (1) with non-negative initial data verifies the more general assumptions stated in [13,14] for biochemical systems, which ensure the wellposedness and positivity of the solution. Comparison theorems guarantee the positivity of the solution of the complete RothC system when g(t), f (t) are assumed positive. Definition 1. We define as the SOC indicator of the continuous RothC model (1) the function SOC(t) = c iom + c dpm (t) + c rpm (t) + c bio (t) + c hum (t) for t ≥ t 0 , where c iom denotes the constant carbon content in the compartment IOM. Theorem 1. Let α, β > 0 and δ := 1 − α − β > 0; suppose 0 < g(t) + f (t) ≤ B and ρ(t) are uniformly bounded from below by a constant µ > 0. The set: µ δ k min with k min = min i k(i) is positively invariant and globally attractive for Model (1).
Under the hypothesis of Theorem 1, the set: is globally attractive for the SOC model: Finally, to complete the understanding of the mathematical features of the RothC model, we provide the following theorem: Theorem 2. Suppose g(t) and f (t) are integrable on every finite subinterval of [t 0 , +∞). If Proof. Under the assumption α + β = 1, the unit vector e ∈ ker(A T ). We haveṠOC = e Tċ = ρ e T A c + g e T a (g) + f e T a ( f ) = g + f .

Long-Term Solutions
When the RothC model is used in real application, the first step is to run the model to equilibrium to calculate the required carbon inputs needed to match the initial SOC content measured. Hence, being able to compute and predict the long-term solution are extremely valuable in order to avoid long-run simulations, which can lead to numerical artefacts if numerical tools are not properly used [15].

Theorem 2 indicates that:
• if there is no external carbon input (g = f = 0), then the SOC indicator is a linear invariant for the SOC model (3); • if the carbon input g + f is replaced by its average valueĝ +f , then SOC(t) grows linearly in time; • if g(t) and f (t) are integrable on every finite subinterval of [t 0 , +∞) and the improper integral +∞ t 0 (g(s) + f (s)) ds converges to its value C ∞ , the indicator SOC R (t) increases in time tending to the finite amount SOC * := SOC(t 0 ) + C ∞ (steady-state solution); • if both g(t) and f (t) are integrable on every finite subinterval of [t 0 , +∞) and t t 0 (g(s) + f (s)) ds is a bounded periodic function of period T, then SOC(t) oscillates indefinitely with period T (periodic solution). Theorem 3. Suppose α, β, δ > 0 with α + β < 1 and k(i) > 0, i = 1, . . . 4. For ρ, g, f , not varying with time, the RothC model admits a unique positive globally stable equilibrium, which has the following expression: Consequently, the SOC indicator has SOC * cont = c iom + e T c cont as the equilibrium, which Proof. From the assumptions, it follows that 0 < α, β < 1. The eigenvalues of the matrix A are given by λ 1 = −k dpm < 0, λ 2 = −k rpm < 0, and λ 3,4 are the roots of the second0order polynomial: with discriminant ∆(λ) = (α − 1) 2 k 2 bio + (β − 1) 2 k 2 hum + 2 k bio k hum (α(β + 1) + β − 1) > 0. Consequently, from Descartes' rule of sign, λ 3 , λ 4 < 0. The matrix A turns out to be negative definite and admits an inverse. It is trivial to check that det(A) = δ For constant (positive) values of ρ, g, f , the equilibrium is achieved by solving the linear system ρ A c = −b, i.e., with: where we adopted the convention that k i,j,m := k(i) k(j) k(m) with i, j, m ∈ {1, 2, 3, 4}. Finally, the stability of the equilibrium is guaranteed as A is negative definite.
As concerns the equilibrium of the SOC indicator, first notice that, from (2), it results When climatic and agricultural variables are considered, the functions ρ(t), g(t) and f (t) are chosen to be time varying on a periodical basis, and the system defined by c(t) is expected to tend toward an oscillatory state as t − → +∞. If we introduce ξ(t) := t t 0 ρ(s) ds for all t ≥ t 0 , then the solution of (1) is given by: The study of the eigenvalues of ξ(t) A enables characterizing the solution behaviour. We can first observe that if the eigenvalues of A ξ(t) are negative, c(t) as t → +∞ does not depend on initial conditions c 0 . Secondly, if there exists a periodic solution c(t) with period T, then c(t 0 + T) = c 0 , and the following theorem holds: Theorem 4. Assume that ρ(t), g(t) and f (t) are periodic with period T. If α + β = 1 and k(i) = 0 for all i ∈ [dpm, rpm, bio, hum], then I − e ξ(t 0 +T) A is not singular. Starting from: the RothC model admits a (unique) periodic solution.
Proof. From the periodicity of ρ(t), it follows that: By imposing in (6) that c(t 0 + T) = c 0 , it follows that: Under the assumed hypothesis, A is not singular, and consequently, the matrix I − e ξ(t 0 +T) A has the inverse. Exploiting the periodicity of b(t) inherited by the periodicity of both g(t) and f (t), we have that: This proves the periodicity of c(t).
The condition α + β < 1 is always true for the RothC model due to the definition of α and β (see [6]), and thus, the stock of carbon in each compartment tends towards a periodic solution in large times whatever the input values are.
Although the continuous approach (1) gives a simple explicit solution, the function of both the input variables and the parameters of the model, in real applications, firstorder discretized versions are applied. We analyse the original discrete formulation of the RothC model, the Exponential Rosenbrock-Euler (ERE) version and a novel first-order nonstandard procedure, closer to the classical discrete RothC procedure than the ERE model.

Original Discrete RothC
By denoting with I the identity matrix and with the discrete (monthly) formulation of the RothC Version 26.3 [11] in vectorial form is given by: where c n ≈ c(t n ) and the discrete temporal grid t n+1 = t 0 + n ∆t advances with stepsize ∆t = 1. Discrete RothC has been applied using data from long-term experiments across several ecosystems, climate conditions and Land Use (LU) classes. It has been extensively applied in Europe for SOC modelling, and applications of RothC to a long-term experiment in semi-arid conditions in Italy can be found in [16,17].
In the case when g(t), f (t) and ρ(t) do not depend on time t, then b = g a (g) + f a ( f ) ; by setting F(∆t) = Λ + (I − Λ) e −∆t ρD , one can demonstrate that, whenever k(i) = 0, (I − F(∆t)) has an inverse, and the system yields a steady-state solution.
From the definition of c * cont in (5) and by noticing that the matrix A that defines the continuous problem (1) verifies the relation A = (Λ − I) D, it follows that: As the original discrete RothC model is applied with ∆t = 1, we can estimate the deviation of the equilibrium of the continuous model with respect to the equilibrium c * RothC (1) of the original discrete RothC model. Given the matrix norm · induced by a vector norm, as −ρ D is negative definite, it results that ϕ(−ρ D) . This indicates that the equilibrium c * RothC (1) is, in norm, an overestimation of the theoretical equilibrium c * cont . The relative error depends on the stepsize ∆t according to: As concerns the equilibrium of the SOC indicator, it results that: and consequently, where: a vector with positive entries that verifies lim ∆t →0 w(∆t) = 0.
Since the original discrete RothC model is applied with ∆t = 1, the deviation of the SOC * cont of the continuous model with respect to the equilibrium SOC * RothC (1) of the original discrete RothC model is given by SOC * RothC (1) = SOC * cont + w T (1) c * cont , thus indicating that the evaluation of soil organic carbon by means of the discrete RothC model is an overestimation of the theoretical value SOC * cont . Usually, ρ, g and f vary through time, but it can be assumed that they have a periodic behaviour. Typically, if the agricultural practices are cyclic and if the weather conditions can be considered periodic, then ρ = ρ(t), g = g(t) and f = f (t) will also behave periodically. Assuming that the periodicity of these variables is T = N∆t, one looks for a solution of c such that c 0 = c N . Then, we can write: for n = 0, . . . N − 2 and then impose the periodic condition c N = c 0 : Setting F n := F n (∆t) := Λ + (I − Λ) e −∆t ρ(t n )D , the above relations can be reformulated as:

Exponential Rosenbrock-Euler Model
If we regard the stepping procedure in (8), which defines the original discrete RothC model as a first-order approximation of the solution of the continuous model (1), then different discrete RothC models can be formulated. As an example, the discretization of Equation (1) from t n = t 0 + n ∆t to t n+1 = t n + ∆t by means of the exponential Rosenbrock-Euler model (for non-autonomous systems) [10,18,19] leads to: Notice that A is negative definite, and for z < 0, it results in 0 < ϕ(z) < 1. In [10], it was shown that A = (Λ − I) D and: Of course, the approximated solutions via the Exponential Rosenbrock-Euler (ERE) method (13) differ from the values given by the original discrete RothC model (8); the major consequence is that the constant discrete steady-state solution, which for the ERE method coincides with the continuous equilibrium c * cont , does not depend on ∆t. As concerns its stability, given A negative definite, it is enough to notice that the eigenvalues of e ∆t ρ(t n ) A are positive, but all less than one.
To find periodic solutions via the ERE method, we can generalize the approach followed for the discrete original RothC model as follows. Suppose that T = N ∆t, and let us impose that c 0 = c N : for n = 0, . . . N − 2 and: with b n := ϕ(∆t ρ(t n ) A) b n , for n = 0, . . . , N − 1.
The sequence of states c n , n = 0, . . . N − 1 characterizes the oscillatory state of the carbon stock in each compartment. Of course, the solution differs from the periodic solution provided in (12).

Novel Non-Standard Discrete RothC Model
The ERE procedure described in (13) belongs to the more general class of non-standard finite difference schemes [14,20]. Different first-order non-standard approximations can be used, depending on the function of the stepsize used to advance the first-order procedure. Each of them can be considered as discrete RothC models, alternative to the ERE model.
A discrete model, closer to the original formulation of the RothC model, can be introduced as follows. Consider the discrete formulation of RothC (8). Simple evaluations lead to the equivalent expressions: are the corresponding eigenvectors. A novel Non-Standard (NS) discrete RothC model is given by: Indeed, we can prove that: (17) is a non-standard first-order approximation of the continuous model (1), i.e., ∆t ϕ(∆t ρ(t n ) A) = O(diag(∆t)), for ∆t → 0.
Proof. From the definition of the exponential matrix, we have that: It follows that ∆tϕ(∆t ρ(t n ) A) = O(diag(∆t)), for ∆t → 0.
By comparing the ERE flow (13) with the above non-standard formulation (17), we notice that the main difference is in the replacement of the matrix A with the matrix A, which has a simple representation in Jordan form. This simplifies the evaluation of the matrix function ϕ(∆tρ(t n ) A) because it can be evaluated on the known eigenvalues −k(i) and eigenvectors in (16).
Moreover, the comparison of the original discrete formulation of the RothC model (8) with the novel NS procedure (17) allows us to notice that, different from the ERE method (13), the decomposition dynamic is now treated in the same way as the original discrete RothC model. The time updating procedure related to the carbon amount entering the system at time t + ∆t is now evaluated up to a factor given by ϕ(∆t ρ(t n ) A) , different from the ERE approximation, which advances evaluating ϕ(∆t ρ(t n ) A).
As the ERE method, in the case of constant functions ρ, g, f , the novel method NS has c * cont as the steady-state solution.
Proof. It is enough to prove that the eigenvalues of F n (∆t) = Λ + (I − Λ) e −∆t ρ(t n )D are in modulus less than one. It is easy to see that two eigenvalues are given by e −∆t ρ(t n ) k dpm and e −∆t ρ(t n ) k rpm . The others are the eigenvalues of the sub-matrix: The eigenvalues lie in the union of the two Gershgorin sets: Notice that F 3,4 has positive entries. Set evaluate: Similarly, Hence it results that the eigenvalue of F 3,4 are less than one in modulus.
To search for periodic solutions via the novel non-standard model, we need to solve the linear system: whereb n := ϕ(∆t ρ(t n ) A) b n , for n = 0, . . . , N − 1. As already observed, we have the same coefficient matrix of the system (12), while the knowledge of the explicit Jordan form of the matrix A simplifies the evaluation of the coefficientsb n .
In the next two sections, we compare the behaviour of the different discrete models on the evaluation of steady-state and long-term periodic solutions. Firstly, the comparison is made on a theoretical basis comparing the accuracy of the methods in approximating the solutions of the continuous model; secondly, we tested both the ERE and the novel nonstandard model with respect to the original discrete RothC model on a classical monthly time-scale experiment where real measurements were available.
A freely accessible MATLAB routine named NSRothC [21] was implemented to replicate all the simulations presented in the next sections. The package includes two versions. The first one (contNSRothC) allowed us to run the model with different stepsizes, when the continuous periodic input functions ρ(t) and b(t) have the particular form presented in Section 4. In the second version (monNSRothC), the stepsize was fixed to one month, and the discrete monthly values of input residuals and FYM were required, as well as the monthly values of weather variables (temperature, rainfall and moisture). As an example, data from the Hoosfield spring barley experiment [6] were used.

Numerical Tests
Let us compare long-term solutions obtained by using the original discrete RothC model, the exponential Rosenbrock-Euler method and the novel non-standard first-order scheme. We set parameters α = 0.1, β = 0.12, γ = 0.59 and η = 0.49. The vector of the decomposition rates k = [0.8333, 0.0250, 0.0550, 0.0017] and functions ρ(t), g(t) and f (t) are supposedly expressed on a monthly scale. They vary with the same period T = 12 and have Gaussian distributions according to: Values a ρ = 0.3561, µ (2) ρ = 1.4996 were chosen, while dispersion coefficients were given by σ Figure 2 on the left.

The function ρ(t) in a time-span of three years is shown in
Values The four components of the input function b(t) := g(t) a (g) + f (t) a ( f ) are depicted in Figure 2, on the right.

Steady-State Solution
To compare long-time periodic and steady-state solutions, we consider the average values of ρ(t) and b(t) over one period T: We recall that both the ERE and the NS method have c * cont as the steady-state solution of their discrete flows, and consequently, they provide the value of SOC * cont as the equilibrium for soil organic carbon; differently, the steady-state solution of the discrete RothC model (9) depends on the chosen stepsize ∆t. In Table 1, we report the equilibrium solutions c * RothC (∆t) for decreasing values of the stepsize ∆t. When we proceed with ∆t = 1, we obtain the original monthly time-stepping procedure, while ∆t = 1/30 represents a daily temporal updating. Notice that, as predicted by the theoretical results, c * cont 2 = 62.6516 ≤ c * RothC (1) 2 = 62.695, i.e., the original discrete RothC model overestimates in norm the theoretical equilibrium value.
The reduction of the temporal stepsize of 1/30 causes approximately the same reduction of the absolute errors as .534, with this confirming the first-order accuracy of the approximation of the steady-state solution. In Figure 3a, we report the relative errors in correspondence of halved stepsizes starting from ∆t = 1 together with their bounds evaluated in (10). In the same figure, we report the error on the equilibrium value of SOC normalized respect to c * cont , i.e., where w(∆t) is defined in (11), with respect to the same reduction of the stepsize ∆t.

Periodic Solutions
In these experiments, we compared the periodic solutions provided by the three discrete models, obtained for the functions ρ(t) and b(t) plotted in Figure 2. Different from the evaluation of steady-state solutions, the ERE and the novel NS procedure, as well as the discrete RothC model provide long-term periodic solutions affected by their own numerical errors, also depending on the chosen stepsize ∆t. As shown in Figure 3b, the RothC discrete model provides the worst performance when compared with the other procedures in terms of errors with respect to the reference solution. In Figure 4, we plot the related SOC indicator (1) in a one-period time span obtained with the three discrete procedures applied with reduced stepsizes ∆t = 1, 1/2, 1/4, 1/8. In the same figure, the reference solution obtained in the long run with the ode45 MATLAB function, with tolerance set at machine precision, is plotted. Notice that all the methods need to be applied with a ∆t << 1 in order to have an acceptable level of accuracy when compared to the reference solution. Again, the original discrete RothC model was confirmed as a less accurate integrator of the continuous model (1) with respect to both the ERE and NS models. In our second experiment, we started at t 0 = 1 with c 0 = c * cont in (18), and we proceeded until the final time T f = t 0 + 15 T was reached. We recall that the period was set at T = 12, and we used the monthly, weekly and daily update by setting ∆t = 1, ∆t = 1/4 and ∆t = 1/30, respectively. In Figure 5, we report the values of SOC obtained for the three different procedures. Two main observations have to be made: first, the initial value c * cont corresponding to the equilibrium solution with mean valuesρ andb was higher than the attained value of SOC reference value at T f when ρ(t) and b(t) were not averaged. This indicates that the temporal oscillations of ρ(t) and b(t) around their mean values cannot be neglected when evaluating, through the SOC indicator, the achievement of land degradation neutrality in the fifteen years from 2015-2030. Second, whatever discrete model is chosen, a qualitative long-term accordance between an approximate value of the SOC indicator and its theoretical solution needs, at least, a weekly update procedure.

The Hoosfield Spring Barley Experiment
We illustrate the model and test the three methods using data from the Hoosfield spring barley experiment, one of the classical long-term experiments carried out at the Rothamster Experimental Station [22]. The same dataset was used as an example of the use of the RothC model in [6]. The Hoosfield experiment was conducted from 1852 till 2000. Spring barley has been grown continuously in the whole period, except in the years 1912, 1933, 1943 and 1967, when the experiment was fallowed to control weeds. The initial SOC content in the soil at 1852 was measured as 33.8632 t C ha −1 , split out in the soil compartments in the following way: 2.7 t C ha −1 in IOM, 0.1533 t C ha −1 in DPM, 4.4852 t C ha −1 in RMP, 0.6671 t C ha −1 in BIO and 25.8576 t C ha −1 in HUM.
The ρ function was estimated on a monthly basis from the weather dataset in [6], including average monthly air temperature, monthly open pan evaporation and monthly rainfall. It was also affected by the percentage of clay in the soil (23.4% in Rothamsted soil) and by the monthly soil cover factor (covered or fallow). Assuming that the soil was covered only from April to July for each crop year, the ρ function assumed the values in Table 2. The other parameters of the model were set as α = 0.10, β = 0.12, γ = 0.59 and η = 0.49. Moreover, the monthly decomposition rate vector k = [0.8333, 0.0250, 0.0550, 0.0017] was considered. Three different scenarios were simulated with the three numerical methods analysed in Section 3, and the results were compared with direct observations of SOC quantity in the soil in 1882, 1913, 1946, 1975, 1982 and 1987: Scenario 1. Unmanured treatment: For this scenario, the annual input of plant residuals was calculated to be 1.60 t C ha −1 y −1 , distributed as in Figure 6, in every year except in those that were fallow. For these years, a null plant residuals input was considered. In this scenario, FYM input was zero.

Scenario 2. Farmyard manure:
The second treatment included inputs from both plant residuals and FYM, as scheduled in Figure 7. The annual input of plant residuals was calculated to be 2.8 t C ha −1 y −1 , greater than the one in Scenario 1, due to the farmyard treatment.

Scenario 3. Mixed treatment, farmyard manure till 1871:
In this scenario, farmyard manure was assured every February till 1871 and nothing thereafter. This caused an annual plant residual input equal to 2.8 t C ha −1 y −1 from 1852 till 1876 and equal to 1.60 t C ha −1 y −1 in the following years, except for the fallow ones. Input data are shown in Figure 8.     11 show the modelled data for total soil organic C in the three treatments, together with the measured data. The modelled results for Scenario 3 when the treatment considered FYM for only the first twenty years were considerably lower than the measurements; agreement was closer with the other two treatments (Scenarios 1 and 2). To test the models' ability to predict the achievement of land degradation neutrality, note that all the simulations gave results in agreement with measured data, i.e., the loss of neutrality in 2000 for Scenario 1 and the achievement of neutrality for Scenario 2 with respect to the initial value in 1852. Measurements indicated the achievement of land degradation neutrality also for Scenario 3; in this case, however, the three models failed in predicting for SOC(2000) a value lower than the initial one, i.e., SOC(1852) = 33.80 t C ha −1 .  To assess the performance of the discrete models and compare measured and simulated valued, as in [16], we used two well-known statistical indices: RMSE (Root Mean Squared Error) and EF (modelling Efficiency): where O i and S i are observed and simulated SOC at the ith value, O is the mean of the observed data and N is the number of observations. The closer the RMSE is to zero, the better the simulated solution describes the data. EF can range from −∞ to one, with the best performance at EF = 1. Negative values of EF indicate that the observed mean is a better predictor than the model.
The results of the evaluation of the above indicators related to the three discrete models applied to approximate Hoosfield barley experiments data for the three different scenarios are reported in Table 3. The obtained values confirm that Scenario 3 was the worst case, i.e., the case when simulated data were more different from the measurements. The three discrete models provided similar values for the simulation of the three different treatments, and we could not select a method that clearly outperformed the others. However, from Table 3, we could set the ERE model as the best model for Scenario 1, while the RothC model provided better results for both Scenarios 2 and 3. Comparing the indicators for all scenarios with respect to the approximation of the experimental data, the minimum value of RMSE = 1.1596 was assumed by the ERE model in Scenario 1, while the maximum value of EF = 0.9609 was assumed by the RothC method in Scenario 2.

Conclusions
Soil organic carbon is one of the key indicators of land degradation status. In this paper, we analysed the RothC model, a simple tool developed for predicting the dynamic evolution of the content of soil organic carbon under the effect of weather conditions and land use data. Both the continuous version, based on a linear, non-autonomous differential system, and the original discrete monthly time-stepping procedure were considered. The aim of our study was to compare the qualitative analysis of both continuous and discrete dynamics in approximating steady-state and long-term periodic solutions. We focused on these aspects since they provide the information concerning the possible achievement of land degradation neutrality. We found that the steady-state solutions of the original discrete model were first-order accurate approximations of the steady-state solutions of the continuous differential model. The accuracy of the approximation represents the main weakness of the RothC discrete model when considered as a first-order accurate approximation of its continuous version. We point out that well-established numerical procedures such as, for example, the Exponential Rosenbrock-Euler (ERE) method, have, by construction, the same equilibria of the approximating differential system. Conversely, the discrete RothC model overestimates the steady-state equilibrium of the continuous flow, so that, in order to obtain the correct estimate with first-order accuracy, we have to reduce the stepsize of the updating time procedure. This fact may have a negative effect in real applications. It is necessary to run the RothC model to equilibrium first, in order to align the initial content of soil organic carbon with the measured one [6]. Nevertheless, the good matching of real data and numerical approximations, together with the available implementation in the RothC 26.3 open-source interface [6], justifies the wide use in the literature of the original discrete RothC model and explains the lack of use of more accurate numerical integrators. For this reason, in this paper, our additional objective was to propose a novel non-standard first-order procedure, the so-called NS discrete model, able to approximate the decomposition dynamics in the same way as the original discrete RothC model. This procedure features the same steady-state solution of the continuous model and exhibits a computational cost lower than the one required by the ERE time-stepping procedure. We also provided a simple code that implements both the original discrete RothC model and the monthly time-stepping NS and ERE models in a MATLAB environment. That allows the reproduction of our results and facilitates future comparisons among all the approaches. The discrete models were firstly tested on a hypothetical example as first-order accurate numerical integrators, and secondly, they were tested on the classical Hoosfield barley experiment to evaluate their ability in reproducing observed data. If we consider the discrete RothC model as a difference scheme, which approximates the continuous RothC differential system, numerical tests showed its coarseness with respect to both the ERE and NS procedures. When applied to the Hoosfield barley experiment, the three discrete models provided very similar results. However, the evaluation of the statistical error indicators still revealed the discrete RothC model as the best scheme for approximating real data in both Scenarios 2 and 3, while the ERE model outperformed the discrete RothC and the novel procedure in Scenario 1. In predicting the long-term behaviours, which are useful to establish the achievement of land degradation neutrality, the three discrete models failed in giving results in accordance with theoretical values in our hypothetical test, as well as with real data in the case of Scenario 3 of the Hoosfield spring barley experiment. This motivates our future research to consider more complex SOC dynamics. In particular, we will analyse the MOMOS model [23,24], which puts the microbial community at the centre of all transformation processes in the cycle of organic matter, from assimilation by degradation to loss by mineralization. Further analysis will also consider nonlinear SOC dynamics [25,26] and spatially explicit models described by partial differential equations [27,28] in order to deal with real domains, as protected areas, covered by non-homogeneous land use patterns. Finally, using field measurements and remote sensing information for modelling land degradation will significantly improve the obtained results [29]. The final aim is to select and provide a robust modelling tool for evaluating the trends of SOC stocks in the Alta Murgia National Park, where the achievement of land degradation neutrality is hampered by a combination of anthropic pressures and climate change [30][31][32].