Next Article in Journal
Discrete Integral and Discrete Derivative on Graphs and Switch Problem of Trees
Next Article in Special Issue
Splines Parameterization of Planar Domains by Physics-Informed Neural Networks
Previous Article in Journal
RETRACTED: An Evolutionary Technique for Building Neural Network Models for Predicting Metal Prices
Previous Article in Special Issue
A Dynamical Model with Time Delay for Risk Contagion
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Novel Fractional-Order RothC Model

by
Vsevolod Bohaienko
1,
Fasma Diele
1,*,
Carmela Marangi
1,
Cristiano Tamborrino
2,
Sebastian Aleksandrowicz
3 and
Edyta Woźniak
3
1
Istituto per Applicazioni del Calcolo ‘M.Picone’, National Research Council (CNR), Via Amendola 122/D, 70126 Bari, Italy
2
Istituto di Nanotecnologie ‘NANOTEC’, National Research Council (CNR), Via Monteroni 122/D, 73100 Lecce, Italy
3
Centrum Badan Kosmicznych PAN, Bartycka 18a, 00-716 Warszawa, Poland
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(7), 1677; https://doi.org/10.3390/math11071677
Submission received: 1 March 2023 / Revised: 28 March 2023 / Accepted: 29 March 2023 / Published: 31 March 2023

Abstract

:
A new fractional q-order variation of the RothC model for the dynamics of soil organic carbon is introduced. A computational method based on the discretization of the analytic solution along with the finite-difference technique are suggested and the stability results for the latter are given. The accuracy of the scheme, in terms of the temporal step size h, is confirmed through numerical testing of a constructed analytic solution. The effectiveness of the proposed discrete method is compared with that of the classical discrete RothC model. Results from real-world experiments show that, by adjusting the fractional order q and the multiplier term ζ ( t , q ) , a better match between simulated and actual data can be achieved compared to the traditional integer-order model.

1. Introduction

Soil carbon dynamics models describe the processes that govern the exchange of carbon between the soil and the atmosphere. They help in understanding the role of soils in the global carbon cycle and the impact of land use and management practices on soil carbon stocks [1]. One of the most widely used models in this field is the Rothamsted Carbon Model [2,3] (RothC), which is a process-based model that simulates soil organic matter dynamics in response to changes in land use and management. RothC considers a number of important factors, including the decomposition of organic matter, changes in soil moisture and temperature, and the effect of ploughing or other agricultural practices on soil carbon capture and storage. By providing a comprehensive representation of soil carbon dynamics, RothC has become a valuable tool in evaluating the impact of land use and management strategies on soil carbon sequestration and in making predictions about trends in soil carbon stocks under the effect of rising temperatures due to global climate change [4,5].
The RothC model analyzes the dynamics of soil carbon by splitting it into five compartments, four active with different chemical degradation characteristics, and one inactive that is resistant to decomposition. The original monthly time-stepping procedure was recently improved by the non-standard approaches proposed in [6]. As in real soil systems, the processes are continuous in time; the following formulation of the RothC model was given in [7]:
d c d t = ρ ( t ) A c + b ( t )
where c ( t ) = [ c d p m ( t ) , c r p m ( t ) , c b i o ( t ) , c h u m ( t ) ] T indicates the concentration of different compartments: decomposable plant material (DPM), resistant plant material (RPM), microbial biomass (BIO) and humified organic matter (HUM). The rate modifier ρ ( t ) depends on land cover, soil, and climatic variables. The model (1) is given in terms of the matrix A, defined as
A = k d p m 0 0 0 0 k r p m 0 0 α k d p m α k r p m ( α 1 ) k b i o α k h u m β k d p m β k r p m β k b i o ( β 1 ) k h u m .
where k = [ k d p m , k r p m , k b i o , k h u m ] are the decomposition rate constants, while α and β represent the fractions of metabolized carbon incorporated into the compartments c b i o ( t ) and c h u m ( t ) ; the remaining part 1 α β is the fraction of metabolized carbon lost by the system in the form of C O 2 .
At time t, the carbon level entering the system is represented by the vector
b ( t ) : = g ( t ) h ( g ) + f ( t ) h ( f )
where
h ( g ) : = [ γ , 1 γ , 0 , 0 ] T , h ( f ) : = [ η , η , 0 , 1 2 η ] T .
The quantities g ( t ) h ( g ) and f ( t ) h ( f ) represent, respectively, the entry of plant residues and farmyard manure (FYM) with the fraction inputs 0 γ 1 , 0 η 1 / 2 , which sum up to 1. DPM and RPM are suitable compartments for the carbon produced by vegetable waste to penetrate the terrain, while DPM, RPM, and HUM are suitable compartments for the carbon quantity of FYM. To ensure the positivity of the RothC solution, we assume 0 < α + β 1 and ρ ( t ) > 0 , when g and f are both positive [6].
Soils, as was previously supposed in, e.g., [8,9], can be considered as media of fractal structure. Their properties under such an assumption are dependent on scale, and memory effects can be pronounced in the behavior of the processes within them. One such class of processes are chemical reaction processes, to which the decomposition of carbon compounds belongs. To increase the accuracy of such processes’ description, integro-differential models of fractional order are applied in, e.g., [10,11]. They include the so-called fractional-order derivatives that generalize the concept of derivatives and can describe non-local interactions [12]. Regarding the considered RothC model that simulates the dynamics of four species’ concentrations without taking spatial effects into account, memory effects, which can be caused, e.g., by the transition of carbon compounds into and out of the forms that do not exhibit decomposition, similarly to the case of adsorption processes [13], are usually modeled by fractional-order derivatives with respect to the time variable. From another point of view, as the RothC model describes decomposition processes in a generalized manner, some not considered additional influences can lower the accuracy of simulation. Generalization of the model to the fractional-order case can be considered as a procedure that adds degrees of freedom to it, allowing better compliance with measurements through the fitting of such additional parameters and orders of derivatives. As a case of the application of such an approach, see, e.g., [14].
Thus, compared to the classical integer-order RothC model, we propose its generalization that takes into account the fractal properties of soil. In this context, considering that, to the best of our knowledge, the behavior of solutions for fractional-order generalizations of the RothC model is by now an open problem, we propose the statement of such a model with the Caputo fractional derivative with respect to the time variable and a numerical technique for the solution of an initial value problem for such a model; we study the peculiarities of such a problem’s solutions and apply the proposed technique to real-world data. The usage of the Caputo fractional derivative is here substantiated by the need to describe known initial values of SOC compartments. Other Caputo-type derivatives with singular kernels could also be used.
The paper is organized as follows. An introduction to fractional derivatives and their components is provided in Section 2, along with theoretical results about the existence, uniqueness, and positivity of the solution. Section 3 provides discretizations of the fractional-order model for numerical simulations and performs a stability analysis of the numerical scheme. Section 4 gives real-world examples of how the fractional-order continuous RothC model can be used to solve classical and novel scenarios including two case studies in Ukraine and Poland. The results show the effectiveness of the fractional model compared to the classical integer-order model. In Section 5, the paper concludes with a summary of the advantages of the fractional-order approach with respect to the integer-order one and with a brief discussion of the implications of our findings.

2. The Novel Fractional-Order Continuous RothC Model

Fractional differential generalization of the RothC model (1) can be obtained by replacing the integer-order derivatives in its left-hand side with one of the fractional-order derivatives. Introducing the Caputo derivative [15,16] for a scalar function y : [ 0 , ) R in the form
t 0 D t q y ( t ) = 1 Γ ( 1 q ) t 0 t y ( t ) ( t τ ) q d τ , 0 < q 1
where Γ ( x ) : = 0 + t x 1 e t d t is the Euler Gamma function, we make use of the fact that the continuous RothC model can be deduced from the mass conservation law and the equations of chemical reaction kinetics. Considering the fractional differential form of the mass conservation law [17], we obtain the fractional RothC model
t 0 D t q c ( t ) = ρ ¯ ( t ) A c ( t ) + b ¯ ( t ) , t t 0 c ( t 0 ) = c 0
where c 0 0 denotes the vector of the initial concentrations and
t 0 D t q c ( t ) : = t 0 D t q c d p m ( t ) , t 0 D t q c r p m ( t ) , t 0 D t q c b i o ( t ) , t 0 D t q c h u m ( t ) T .
As the processes described by the model (3) occur in a single medium, we consider the orders of fractional derivatives in each equation to be equal to q. The novel functions ρ ¯ ( t ) and b ¯ ( t ) appearing in (3) are necessary to transform the ordinary media dimensions into the fractal media one (for an example of such approaches’ application, see, e.g., [18]). They are obtained by dividing the rate modifier ρ ( t ) and the input functions f ( t ) , g ( t ) by the term ζ ( t , q ) , which has dimension t i m e q 1 to ensure the dimensional correctness of Equation (3):
ρ ¯ ( t ) : = ρ ( t ) / ζ ( t , q ) , b ¯ ( t ) : = b ( t ) / ζ ( t , q ) .
Accordingly, for q = 1 , we recover the classical-order derivative by setting ζ ( t , 1 ) equal to the dimensionless value 1. Examples of ζ functions are given in [18]. The simplest choice ζ ( t , q ) = 1 means that the rates in a medium, where memory effects are observed, are equal to the rates in an ordinary medium. Another choice that can be used for the fractional-order RothC model is ζ ( t , q ) = t q 1 , t t 0 > 0 .

2.1. Existence, Uniqueness, Positively Invariant Set, Equilibria

To prove the existence and uniqueness of the system in (3)’s solution, we use the results for the general system of fractional ordinary differential equations (FODE) presented in [19]. Application of Theorem 3.1 in [19] to the system (3) yields that its solution c ( t ) exists if f ( t , c ( t ) ) : = ρ ¯ ( t ) A c ( t ) + b ¯ ( t ) , with f : [ 0 , ) × R 4 R 4 , is Lebesgue measurable with respect to t, continuous with respect to c ( t ) , and there exist two positive constants ω ¯ , λ ¯ > 0 such that
f ( t , c ( t ) ω ¯ + λ ¯ c ( t )
where · is the supremum norm. The latest inequality holds as
f ( t , c ( t ) = | | b ¯ ( t ) + ρ ¯ ( t ) A c ( t ) b ¯ ( t ) + ρ ¯ ( t ) A c ( t ) = ω ¯ + λ ¯ c ( t ) , ω ¯ : = b ¯ ( t ) , λ ¯ : = ρ ¯ ( t ) A .
According to Remark 3.2 in [19], the solution is also unique if the Jacobian f ( t , c ( t ) ) c is continuous. The model (3) has the same positively invariant and globally attractive set as the set for the integer-order model given in Proposition 1 and Theorem 1 in [6]. The equilibrium state of the model (3) is obviously the same as given in Proposition 2 and Theorem 3 in [6] for the case of the integer-order RothC model.

2.2. Analytic Solution

Applying in a straightforward manner Theorem 5 in [20] to the FODE (3) itself, we have its analytic solution in the form
c ( t ) = Ψ ( t , t 0 ) c 0 + t 0 t Φ ( t , τ ) b ¯ ( τ ) d τ , Φ ( t , τ ) : = m N τ J t m q [ ρ ¯ ( t ) A ] , Ψ ( t , t 0 ) : = m N t 0 J ¯ t m q [ ρ ¯ ( t ) A ] ,
where, for all the continuous matrix functions M ( t ) , the Riemann–Liouville fractional integral is defined as
s J t q [ M ( t ) ] : = 1 Γ ( q ) s t ( t τ ) q 1 M ( τ ) d τ , s < t
while both s J t m q [ M ( t ) ] and s J ¯ t m q [ M ( t ) ] are recursively defined as
s J t 0 q [ M ( t ) ] : = ( t s ) q 1 Γ ( q ) I , s J t ( m + 1 ) q [ M ( t ) ] : = s J t q M ( t ) s J t m q [ M ( t ) ] , m = 0 , 1 , ,
and
s J ¯ t 0 q [ M ( t ) ] : = I , s J ¯ t ( m + 1 ) q [ M ( t ) ] : = s J t q M ( t ) s J ¯ t m q [ M ( t ) ] , m = 0 , 1 , ,
with I the identity matrix. Let us also note that, for q = 1 ,
s J ¯ t m 1 [ M ( t ) ] = s J t m 1 [ M ( t ) ] = 1 m ! s t M ( τ ) d τ m ,
therefore, Φ ( t , t 0 ) [ ρ ¯ ( t ) A ] = Ψ ( t , t 0 ) [ ρ ¯ ( t ) A ] = e ξ ( t ) A , and Equation (4) transforms into the well-known solution of the integer-order model.

2.3. Soil Organic Carbon Content

Following the definition in [6], let us consider the SOC indicator as s o c ( t ) : = e T c ( t ) , with e = [ 1 , 1 , 1 , 1 ] T . Then, from (3), in the case of ζ ( t , q ) > 0 , we have
t 0 D t q s o c ( t ) = g ¯ ( t ) + f ¯ ( t ) ρ ¯ ( t ) δ k T c ( t ) B μ δ k m i n s o c ( t ) ,
where g ¯ = g / ζ ( t , q ) , f ¯ = f / ζ ( t , q ) and k m i n = min k d p m , k r p m , k b i o , k h u m . Using the solution given in Theorem 2 in [20] of the equation t 0 D t q f ( t ) = B A · f ( t ) that corresponds to the inequality (5), we have
s o c ( t ) B μ δ k m i n B μ δ k m i n s o c ( t 0 ) E q ( μ δ k m i n ( t t 0 ) q )
where E q ( z ) : = n = 0 + z n Γ ( q n + 1 ) is the Mittag–Leffler function. By noticing that the Mittag–Leffler function generalizes the exponential function being equal to it when q = 1 , the above-mentioned statement about the positively invariant and globally attractive set of the model (3) follows.
In the case of no C O 2 release ( α + β = 1 ), the soil carbon content dynamics have the form
s o c ( t ) = s o c ( t 0 ) + 1 Γ ( q ) t 0 t ( t τ ) q 1 ζ ( t , q ) ( f ( t ) + g ( t ) ) d τ .
These dynamics are equal to the dynamics in the integer-order case when q = 1 .
When f ( t ) + g ( t ) ζ ( t , q ) is constant, s o c ( t ) grows non-linearly as ( t t 0 ) q in the fractional-order model, compared with linear growth in the integer-order one. Such behavior is similar to the inverse power of the time dynamics of water content in unsaturated soil considered as a fractal medium, as reported in [8].

3. Discretizations of the Fractional-Order Model

Let us consider a computational procedure based on the discretization of the analytic solution (4). In the case of a constant ρ ¯ ( t ) , according to [20],
Φ ( t , τ ) [ ρ ¯ A ] = ( t τ ) q 1 E q , q ( ρ ¯ A ( t τ ) q )
where E a , b M = i = 0 M i Γ ( a · i + b ) is the matrix two-parameter Mittag–Leffler function. For Ψ ( t , τ ) , we have Ψ ( t , τ ) = E q , 1 ( ρ ¯ A ( t τ ) q ) . Thus, for constant ρ ¯ ( t ) , the solution (4) can be expressed as
c ( t ) = E q , 1 ( ρ ¯ A ( t t 0 ) q ) c 0 + t 0 t ( t τ ) q 1 E q , q ( ρ ¯ A ( t τ ) q ) b ¯ ( τ ) d τ .
Introducing a grid ω h = { t n = t 0 + n h ; n = 0 , . . . , N } and considering b ¯ ( t ) = b ¯ ( t n + 1 ) , t n t < t n + 1 , we obtain
c ( t N ) E q , 1 ( ρ ¯ A ( N h ) q ) c 0 + n = 0 N 1 b ¯ ( t n + 1 ) t n t n + 1 ( t N τ ) q 1 E q , q ( ρ ¯ A ( t N τ ) q ) d τ .
As an indefinite integral
( t τ ) q 1 E q , q ( ρ ¯ A ( t τ ) q ) d τ = ( t τ ) q E q , q + 1 ( ρ ¯ A ( t τ ) q )
we have
c ( t N ) E q , 1 ( ρ ¯ A ( N h ) q ) c 0 + n = 0 N 1 b ¯ ( t n + 1 ) ( t N t n ) q 1 E q , q ( ρ ¯ A ( t N t n ) q ) ( t N t n + 1 ) q 1 E q , q ( ρ ¯ A ( t N t n + 1 ) q )
Computation according to (8) is complicated by the need to compute the matrix Mittag–Leffler function with a growing argument value, which leads to fast data type overflow [21]. This fact, accompanying large truncation errors, and the increasing number of terms in the Mittag–Leffler function’s representation needed to compute it with acceptable accuracy, lead to the high computational complexity of this scheme, making it problematic when used for long-term simulations that are needed for SOC level predictions. For non-constant ρ ¯ ( t ) , the operators Φ ( t , τ ) and Ψ ( t , τ ) cannot be simplified. This introduces additional computational complexities.
These computational problems are grounded in the exact solution’s discretization procedure itself, so we consider the methods that, instead, approximate a system of equations producing less difficulties in their application. Such methods include numerous finite difference (see, e.g., [22,23]) and spectral (see, e.g., [24]) ones.
Thus, solely for the purpose of obtaining and studying numerical solutions and not considering the issue of the methods’ performance comparison, we further use the Crank–Nicholson finite difference scheme with L1 approximation of the fractional derivatives of O ( h 2 q ) [25] accuracy order in approximating the fractional-order RothC model (3). Denoting with c ( n , q ) the approximation of c ¯ ( t n ) at t n = t 0 + n h for n = 0 , . . . , N h : = ( t N t 0 ) / h , the scheme starts with c ( 0 , q ) = c ¯ ( t 0 ) and advances at t n + 1 according to
1 Γ ( 2 q ) h q I 1 2 ρ ¯ ( t n + 1 ) A c ( n + 1 , q ) = 1 2 ρ ¯ ( t n ) A c ( n , q ) + b ¯ ( t n + 1 ) + b ¯ ( t n ) + + 1 Γ ( 2 q ) h q c ( n , q ) k = 0 n 1 ( c ( k + 1 , q ) c ( k , q ) ) ( n k + 1 ) 1 q ( n k ) 1 q .
Let us denote
Δ q = 1 Γ ( 2 q ) h q , Δ I = Δ q I , d k n = Δ q ( n k + 1 ) 1 q ( n k ) 1 q , M n ( 1 ) = Δ I 1 2 ρ ¯ ( t n + 1 / 2 ) A , M n ( 2 ) = Δ I + 1 2 ρ ¯ ( t n + 1 / 2 ) A .
Then, the scheme (9) in the notations (10) can be written in a more compact form:
c ( n + 1 , q ) = M n ( 1 ) 1 M n ( 2 ) c ( n , q ) + 1 2 b ¯ ( t n + 1 ) + b ¯ ( t n ) + k = 0 n 1 Δ q d k n M n ( 1 ) 1 c ( k + 1 , q ) c ( k , q ) .
Finally, let us note that the scheme (9) reduces to the classical Crank–Nicholson implicit scheme [26] when q = 1 :
1 h I 1 2 ρ ( t n + 1 ) A c ( n + 1 , 1 ) = 1 2 ρ ( t n ) A c ( n , 1 ) + b ( t n + 1 ) + b ( t n ) + c ( n , 1 ) h .

3.1. Stability Analysis of the Numerical Scheme

Let c ( t n ) be the exact solution of the problem at t = t n with unperturbed initial condition c 0 and let c ˜ ( n , q ) be the numerical solution of the problem at t n corresponding to the perturbed initial conditions c ˜ 0 = c 0 + e 0 induced by the error e 0 . We consider the scheme stable on the step n with respect to the initial conditions if the 2-norm of the error e n = c ( t n ) c ˜ ( n , q ) satisfies | | e n | | 2 C | | e 0 | | 2 , where C is a constant independent of n.
Considering all inputs bounded, from the linearity of (11), the error e n satisfies
e n + 1 = M n ( 1 ) 1 M n ( 2 ) e n + k = 0 n 1 Δ q d k n M n ( 1 ) 1 e k + 1 e k .
Let us first consider the case of n = 0 . Here, we have
e 1 = M n ( 1 ) 1 M n ( 2 ) e 0 , | | e 1 | | 2 = | | M n ( 1 ) 1 M n ( 2 ) e 0 | | 2 ,
| | e 1 2 M n ( 1 ) 1 2 M n ( 2 ) 2 e 0 2 .
The values of ( M n ( 1 ) ) 1 2 and | | M n ( 2 ) | | 2 can be obtained from the eigenvalues λ n ( 1 ) of the matrix M n ( 1 ) given by
λ n , 1 ( 1 ) = Δ q + 1 2 ρ ¯ ( t n + 1 / 2 ) k d p m , λ n , 2 ( 1 ) = Δ q + 1 2 ρ ¯ ( t n + 1 / 2 ) k r p m , λ n , 3 ( 1 ) = Δ q + 1 4 ρ ¯ ( t n + 1 / 2 ) K 1 1 / 2 + K 2 , λ n , 4 ( 1 ) = Δ q + 1 4 ρ ¯ ( t n + 1 / 2 ) K 1 1 / 2 + K 2 ,
and the eigenvalues λ n ( 2 ) of the matrices M n ( 2 ) that have the form
λ n , 1 ( 2 ) = Δ q 1 2 ρ ¯ ( t n + 1 / 2 ) k d p m , λ n , 2 ( 2 ) = Δ q 1 2 ρ ¯ ( t n + 1 / 2 ) k r p m , λ n , 3 ( 2 ) = Δ q 1 4 ρ ¯ ( t n + 1 / 2 ) K 1 1 / 2 + K 2 , λ n , 4 ( 2 ) = Δ q 1 4 ρ ¯ ( t n + 1 / 2 ) K 1 1 / 2 + K 2
where K 1 = ( k b i o k h u m ) 2 ( α 1 ) 2 + 4 α β k b i o k h u m , K 2 = ( k b i o + k h u m ) ( 1 α ) .
From the definition and the properties of the 2-norm, we have, for K 1 1 / 2 K 2 ,
M n ( 2 ) 2 = max i λ n , i ( 2 ) Δ q , ( M n ( 1 ) ) 1 2 = 1 min i λ n , i ( 1 ) 1 Δ q .
Inequalities (13) yield e 1 2 C 1 e 0 | | 2 with C 1 1 .
Proceeding with the induction, suppose that e n 2 C n e 0 2 , with C n 1 . Then,
e n + 1 2 max i λ n , i ( 2 ) min i λ n , i ( 1 ) C n e 0 2 + k = 0 n 1 Δ q d k n ( M n ( 1 ) ) 1 e k + 1 e k 2 .
The maximal value of d k n is achieved for k = n 1 , d k n 2 1 q 1 . Thus, we have
k = 0 n 1 Δ q d k n ( M n ( 1 ) ) 1 e k + 1 e k 2 Δ q 2 1 q 1 ( M n ( 1 ) ) 1 2 e n e 0 2 2 1 q 1 1 C n e 0 2
and further
e n + 1 2 C n + 1 e 0 2 , C n + 1 = max i λ n , i ( 2 ) min i λ n , i ( 1 ) C n + 2 1 q 1 1 C n 1 .
The inequalities (13) hold for K 1 1 / 2 K 2 that yields α β ( 1 α ) 2 . Combining the latter inequality with the restriction 0 < α + β 1 , we obtain α 1 2 . Thus, the considered scheme is stable with respect to the initial conditions when α 1 2 .

3.2. Convergence Testing on a Constructed Analytical Solution

Computational experiments to test the accuracy and convergence of the finite difference scheme were conducted on a constructed analytical solution
c ( t ) = c 0 + ( t t 0 ) r + ( t t 0 ) 2 s
that satisfies Equation (3) for
b ¯ ( t ) = ( t t 0 ) 1 q Γ ( 2 q ) r + 2 ( t t 0 ) 2 q Γ ( 3 q ) s ρ ¯ ( t ) A c ( t ) , ρ ¯ ( t ) = ρ ( t ) / ζ ( t , q ) .
In our experiments, we choose the parameters α , β , γ , η and the decomposition rates k as reported in [6]. We choose ζ ( t , q ) = t q 1 and ρ ( t ) = sin ( ρ 0 ( t t 0 ) + ρ 1 ) + ρ 2 with ρ 0 = 0.1 , ρ 1 = 0.5 , ρ 2 = 1.5 . The initial value is set at c 0 = [ 0.1533 , 4.4852 , 0.6671 , 25.857 ] T , with r = 10 3 · [ 1 , 2 , 1 , 40 ] T , s = 10 4 · [ 0.1 , 1 , 0.05 , 2 ] T , t 0 = 1852 , and the ending time t N = 2000 .
Denote with c i ( t n ) the i-th entry of the theoretical solution c ( t n ) and with c i ( n , q ) the corresponding entry of its numerical approximated solution c ( n , q ) c ( t n ) at t n = t 0 + n h for n = 0 , . . . , N h : = ( t N t 0 ) / h , obtained with the fractional q order scheme (9). The obtained values of average absolute errors
E ( h , q ) = 1 N h n = 0 N h i = 1 4 c i ( n , q ) c i ( t n ) 2 1 2
are given in Table 1. The obtained data confirm the O ( h 2 q ) -order accuracy of the scheme.
In all presented cases, the solutions of the fractional-order model are qualitatively asymptotic periodical (the theoretical proof of this statement in a quantitative case is beyond the scope of this study), with peaks shifted to the lower values of t when compared with the integer-order local behavior (see Figure 1). These shifts are more pronounced with the change in the form of the curves for ζ ( t , q ) = t q 1 (Figure 2). The peak values decrease with the decrease in q for ζ ( t , q ) = 1 and, on the contrary, increase for ζ ( t , q ) = t q 1 .
The influence of the fractional derivatives’ order on c d p m is here slight, with less than 1% difference between the solutions with q = 1 and q = 0.4 for t t 0 > 6 . More pronounced peak shifts and scalings are observed for c r p m and c b i o , with a higher level of changes for ζ ( t , q ) = 1 . The highest level of changes is observed for c h u m .
The average values of concentrations decrease with the decrease in q for ζ ( t , q ) = 1 and, on the contrary, increase for ζ ( t , q ) = t q 1 (see Figure 1 and Figure 2).
The above-described qualitative behavior of the fractional-order model’s solutions can be taken into account while studying real-world datasets for the purpose of selecting the proper form of the factor ζ ( t , q ) and the range of the changes in the order q used for fitting.

4. Real-World Applications

4.1. The Case of “Hoosfield Spring Barley”

To test the new fractional-order model and perform a comparison with the classical model, we consider the data of the well-known case of “Hoosfield spring barley” [3,27]. The Hoosfield experiment took place between 1852 and 2000. Throughout the period, spring barley was cultivated, with an interruption in the years 1912, 1933, 1943, and 1967, due to the control of weeds. The initial SOC content (in t C ha 1 units) in the soil at t 0 = 1852 was measured as c d p m ( t 0 ) = 0.1533 , c r p m ( t 0 ) = 4.4852 , c b i o ( t 0 ) = 0.6671 , and c h u m ( t 0 ) = 25.857 . The values of the function ρ calculated for each month are listed in Table 2, while the other parameters α , β , γ , η , and k have been set as in [6].
Three different scenarios are simulated with the fractional-order RothC model discretized with the first-order scheme—in particular, Scenario 1, named Unmanured Treatment; Scenario 2, named Farmyard Manure, and Scenario 3, named Mixed Treatment, as described in [6]. Six annual real observed SOC values in the soil (1882, 1913, 1946, 1975, 1982, 1987) are compared with the results obtained for several values of the order 0 < q < 1 and with the classical RothC model corresponding to q = 1 : the above-described three scenarios were modeled using the model (3) for q = 1 (classical case) and q = 0.98 , 0.96 , 0.94 , 0.92 with ζ ( t , q ) = 1 , and ζ ( t , q ) = t q 1 discretized by the scheme (9).
The average absolute errors of modeling with a time step equal to 1 month computed similarly to (15) are given in Table 3. The only case wherein the fractional-order model gave a lower error compared with the integer-order one was the case of q = 0.98 and ζ ( t , q ) = t q 1 —here, for Scenario 1, the error did not change significantly; for Scenario 2, it was 19% lower, and for Scenario 3, it was 6% higher. The dynamics of SOC indicator s o c ( t ) in this case are shown in Figure 3. From the obtained solutions, it can be seen that for ζ ( t , q ) = t q 1 and small deviations of q from 1, the fractional-order model changes primarily the form of the curves within yearly cycles, with little influence on the long-term dynamics. Thus, more observations are needed to benefit from the usage of such a model in this case.

4.2. Solution Scheme Comparison

To compare the performance of the considered Crank–Nicholson scheme with the ones used in the literature, we model Scenario 1 in the integer-order case using the scheme (9) for q = 1 with a monthly time step, the non-standard scheme presented in [6], and the basic scheme described in [2].
We study the average absolute differences between the solutions computed similarly to (15) for the whole period and changes in their accumulated values during the simulation (Table 4).
The monthly stepping Crank–Nicholson scheme produced a solution that differed from the other solutions by 2% of the average value of S O C . The greatest differences in terms of the percentage of the average value were observed for the DPM (decomposable plant material) indicator. The average absolute difference between the solutions by the basic and the non-standard models comprised 26% of the average DPM value. It is 9-times more than in the Crank–Nicholson scheme. Let us also note that with the increase in the simulation time, the accumulated average absolute difference in the S O C indicator increases when comparing the basic and the non-standard schemes. A less pronounced increase was observed when comparing the non-standard scheme with the Crank–Nicholson scheme. This means that, despite being close, the solutions produced by these three methods have a tendency to yield an increase in differences with the increase in the simulation period.
The simulation times for the basic and the non-standard schemes implemented in Python were close and equal to 0.37 s (AMD Ryzen 3 5300U CPU was used). In the case of the Crank–Nicholson scheme, it was 0.23 s. In the case of the fractional-order model for q = 0.98 , the differences with the integer-order case (non-standard scheme) comprised 8% of the average S O C value for ζ ( t , q ) = t q 1 and 2% for ζ ( t , q ) = 1 . The execution time in this case significantly increased to 4.9 s. Such a significant execution time makes it efficient to use multi-threading and compiled instead of interpreted languages to implement fractional-order simulation procedures.

4.3. Ukrainian and Polish Scenarios

The proposed model was also applied to simulate SOC dynamics on the basis of data on the agrochemical and physical parameters of the soils in Ukraine and Poland. We selected two sites. For Ukraine, we chose the space point called “Sajivka”, located within an agricultural field. This field is characterized by Chernozem soil formed on loess with high content of carbonates under steppe vegetation. Chernozems are very fertile soils and very high agricultural yields can be obtained. A humid continental climate zone characterizes the site for which data on humus content, clay content, and soil density are available for the years from 2015 to 2020 (https://ismld.com.ua/ (accessed on 27 November 2022)). A weather dataset containing daily maximal, minimal, and average temperatures along with the amount of precipitation was obtained from ftp://ftp.ncdc.noaa.gov (accessed on 27 November 2022) for the weather station located in Kryvyi Rih city.
For Poland, the control point was selected from the database of the Polish Chief Inspectorate for Environmental Protection, which contains measurements of soil properties on a regular basis, every 5 years, since 1995, for a network of 216 control points located on arable land throughout the country (https://www.gios.gov.pl/chemizm_gleb/index.php?mod=monit (accessed on 20 December 2022.)). The site, called “Chrusciel”, is characterized by a fertile Luvisol (IUSS Working Group WRB (2022). “World Reference Base for Soil Resources, 4th edition”. IUSS, Vienna) originating from sandy loam. It is considered of high agricultural usefulness, with complexes (Saturnin Zawadzki (red): Gleboznawstwo. PWRiL, 1999. ISBN 83-09-01703-0.) suitable for the cultivation of sugar beet, wheat, red clover, alfalfa, winter rape, and faba bean. The point is under the climatic influence of the Baltic sea and is characterized by the following conditions: annual mean temperature 8.5 C, mean annual precipitation 710 mm, and growing season length 215 days. The choice of the site was based on its proximity to a meteorological station, Elblag-Milejewo (http://www.meteomanz.com (accessed on 20 December 2022)). Figure 4 presents the spatial distribution of the two points and Table 5 contains their main characteristics.

4.4. Ukrainian and Polish Numerical Experiments

For both scenarios, as only data on humus content are available, we assume that in the initial moment of time, the SOC concentrations of all other compartments but the humidified one are equal to 0. Moreover, as the introduction of manure is not common in large-scale crop production, we consider f ( t ) = 0 . The representation of the density function g ( t ) is given according to [4] in the form g ( t ) = P ( t ) g ^ ( t ) , where P ( t ) is the integral yearly value of g ( t ) and g ^ ( t ) represents the distribution of g ( t ) within the year. The value of g ^ ( t ) , due to the lack of actual data, was taken equal to g ^ r ( a ) ( t ) from Table 1 in [4] ( 1 / 6 for April, May, June; 1 / 2 for July; 0 for other months). Net Primary Production (NPP) estimates (from the MOD17A3HGF.006 database of the NASA Earth Observation System (EOS) program; temporal range from 2015 to 2020. Application for Extracting and Exploring Analysis Ready Samples (AppEEARS)) in the considered points have been used to evaluate the annual values of P ( t ) from the initial value P ( t 0 ) , according to the formula P ( t n + 1 ) = P ( t 0 ) N P P ( t n + 1 ) / N P P ( t 0 ) as in [4]. Given the hypothesis that, at the initial moment of time, the observed values of SOC compounds correspond to their values at a stable equilibrium [4], the value of P ( t 0 ) reported in Table 5 was calculated according to Equation (11) in [4], taking into account the considered assumptions
P ( t 0 ) = 12 ρ ( 0 ) ( r ) δ k h u m c h u m ( t 0 )
The value of ρ ( 0 ) ( r ) = k c ( r ) with k c ( r ) calculated as in [4]
k c ( r ) = 0.6 N p 30 e x ( r ) 1 + e x ( r ) , x ( r ) = 30 r 1 r
with the number of months per year of bare soil set at N b = 4 for the Ukraine point from April to July, while N b = 3 for the Polish site, from August to October. The value of the DPM/RPM ratio r was equal to 1.44 , which corresponds to agricultural crops.
We assume a sufficient moisture supply for the site “Sajivka” as it is located within the lands of irrigation systems, and we assume the respective field to be irrigated. The irrigation system is not taken into account by the RothC model; therefore, it is necessary to add the amount of water due to the irrigation of the crops to the count of the monthly rainfall that occurred during the irrigation period. Therefore, we assume as a soil water deficit the amount of water added with irrigation, as in [28]. Consequently, the soil moisture deficit there should not be significant. For the “Crushiel” scenario, potential evapotranspiration (PET) assessments obtained using GLEAM remote sensing PET products (actual evaporation) [29] were used to evaluate the accumulated soil moisture function.
To assess the performance of different models and the errors between the real observed and predicted values, we use the modeling efficiency index (EF) and the root-mean-square error (RMSE) defined as
E F = i = 1 n R i R ¯ 2 i = 1 n R i P i 2 i = 1 n R i R ¯ 2 R M S E = 1 n i = 1 n R i P i 2
where R are the real SOC field observations and P are the SOC model simulations, n is the number of observations, and R ¯ is the mean. The EF values from 0 to 1 indicate that the model results are better than those using the mean values of observations and values below 0 are worse. The best-performing models are those with EF close to 1. For the RMSE, the values close to zero indicate that the simulated solution describes the data well. Some of the obtained numerical results are shown in Figure 5 and statistical performance is given in Table 6.
For the point "Sajivka", the solution that best matches the observed dynamics was obtained using the fractional-order model for q = 0.97 and ζ ( t , q ) = t q 1 , which reaches the best EF value 0.851173 and RMSE value 1.201889 . We point out that the evaluation of the statistical indices in this case was carried out taking into consideration all the observed measured values of SOC, excluding the year 2018 as it represents an outlier. The use of the fractional-order model with ζ ( t , q ) = 1 yields slight but worse changes compared with the integer-order model in the studied cases with a 5-year simulation period.
For the point “Chrusciel”, the solution that best matches the observed dynamics was obtained using the fractional-order model for q = 0.95 and ζ ( t , q ) = 1 . Here, the use of the fractional-order model with ζ ( t , q ) = t q 1 yields worse changes compared with the integer-order model.
Simulation of the SOC dynamics for the Poland and Ukraine sites showed that the use of the fractional-order model with ζ ( t , q ) = t q 1 seems to correct the simulation dynamics in the cases when, with the best-known values of models’ parameters, the SOC content is overestimated by the integer-order model. Conversely, the choice of a fractional model with ζ ( t , q ) = 1 corrects the simulation dynamics in the cases when the SOC content is underestimated by the model predictions. In general, the use of fractional-order models introduces two added degrees of freedom, the q parameter and function ζ ( t , q ) , which allows us to correct the SOC trend predicted by the classical model. The obtained results in terms of the statistical indicators EF and RMSE confirm the improved performance of the fractional model when compared with the integer-order model.

5. Conclusions

In this study, we formulated the fractional-order generalization of the continuous RothC model and studied the qualitative properties of its numerical solutions. The introduction of the Caputo derivative with respect to the time variable into the model reflects the case of fractal-structured soils in which carbon compartment-related reactions occur. As such an introduction requires the agreement of coefficients’ dimensions, we consider two possible forms of corresponding correction multipliers, obtaining two models with different behavior of the solutions.
The respective solutions were obtained by a finite difference scheme whose conditional stability was investigated. Because of the non-constant coefficients of the model, the usage of a fully numerical method is here justified by the high computational complexity of the exact solution’s discretization, which is the scheme that is widely and effectively applied for integer-order models and fractional-order models with constant coefficients.
The results of the fractional-order model’s application to several real-world scenarios with once-per-year measured SOC values show that it can be used to increase the accuracy of simulation by fitting the order of the fractional derivative. In the case of the simplest correction multiplier ζ ( t , q ) = 1 , significant differences compared to the integer-order model were observed for the carbon compartments’ dynamics within the year. Thus, further studies with a higher frequency of SOC measurements are needed to identify the limits of such a fractional-order model’s usability.
For the Ukranian and Polish real-world scenarios, in the case of the multiplier in the form ζ ( t , q ) = t q 1 , a decrease in SOC content compared with the integer-order model was observed. Conversely, an increase in SOC content compared with the integer-order model was obtained for ζ ( t , q ) = 1 . As such a decrease can also be simulated by the changes in other models’ parameters, particularly the rate modifier ρ ( t ) , additional studies should be performed to distinguish real-world cases in which the influence of each of these parameters prevails.
As a future application, we will investigate the benefits of introducing the fractional-order derivatives for those compartmental models that explicitly consider the microbes’ influence on soil organic matter decomposition and on the long-term stabilization of carbon in soils [30].

Author Contributions

Methodology, V.B. and F.D.; Software, V.B. and C.T.; Validation, V.B., C.T., S.A. and E.W.; Writing—original draft preparation, V.B., F.D. and C.T.; Project administration, C.M.; Funding acquisition, C.M. All authors have read and agreed to the published version of the manuscript.

Funding

Funding was provided for the project by the European Union’s Horizon 2020 research and innovation programme, project ”EOTiST” grant agreement No. 952111. Cristiano Tamborrino has been supported by LifeWatch Italy through the project LifeWatchPLUS (CIR-01-00028). Vsevolod Bohaienko has been supported by grant IAC-BA UA-SC N.01-2022, a specific support measure for the Ukrainian scientific communities.

Data Availability Statement

The developed codes and the datasets used for this work can be released under specific request.

Acknowledgments

This study was carried out within the framework of the H2020 project “Earth Observation Training in Science and Technology” at the Space Research Centre of the Polish Academy of Sciences—”EOTiST”, grant agreement No. 952111. The authors Fasma Diele and Cristiano Tamborrino are members of GNCS Indam. Finally, the authors would like to thank Cosimo Grippa for his valuable technical support.

Conflicts of Interest

The authors declare no competing interests.

References

  1. Diele, F.; Luiso, I.; Marangi, C.; Martiradonna, A. SOC-reactivity analysis for a newly defined class of two-dimensional soil organic carbon dynamics. Appl. Math. Model. 2023, 118, 1–21. [Google Scholar] [CrossRef]
  2. Coleman, K.; Jenkinson, D.; Crocker, G.; Grace, P.; Klir, J.; Körschens, M.; Poulton, P.; Richter, D. Simulating trends in soil organic carbon in long-term experiments using RothC-26.3. Geoderma 1997, 81, 29–44. [Google Scholar] [CrossRef]
  3. Coleman, K.; Jenkinson, D.S. ROTHC-26.3: A Model for the Turnover of Carbon in Soil: Model Description and Users Guide: K. Coleman and DS Jenkinson; IACR: Santa Barbara, CA, USA, 1995. [Google Scholar]
  4. Diele, F.; Luiso, I.; Marangi, C.; Martiradonna, A.; Woźniak, E. Evaluating the impact of increasing temperatures on changes in Soil Organic Carbon stocks: Sensitivity analysis and non-standard discrete approximation. Comput. Geosci. 2022, 26, 1345–1366. [Google Scholar] [CrossRef]
  5. Marangi, C.; Diele, F.; Luiso, I.; Martiradonna, A.; Wozniak, E. SOC indicator of land-degradation: Responses of continuous and non-standard discrete RothC models to environmental changes. In Proceedings of the EGU General Assembly Conference Abstracts, Vienna, Austria, 23–27 May 2022; p. EGU22-5894. [Google Scholar]
  6. Diele, F.; Marangi, C.; Martiradonna, A. Non-standard discrete RothC models for soil carbon dynamics. Axioms 2021, 10, 56. [Google Scholar] [CrossRef]
  7. Parshotam, A. The Rothamsted soil-carbon turnover model—Discrete to continuous form. Ecol. Model. 1996, 86, 283–289. [Google Scholar] [CrossRef]
  8. Pachepsky, Y.; Timlin, D. Water transport in soils as in fractal media. J. Hydrol. 1998, 204, 98–107. [Google Scholar] [CrossRef] [Green Version]
  9. Kong, B.; Dai, C.X.; Hu, H.; Xia, J.; He, S.H. The Fractal Characteristics of Soft Soil under Cyclic Loading Based on SEM. Fractal Fract. 2022, 6, 423. [Google Scholar] [CrossRef]
  10. Singh, A.; Das, S.; Ong, S. Study and analysis of nonlinear (2+1)-dimensional solute transport equation in porous media. Math. Comput. Simul. 2022, 192, 491–500. [Google Scholar] [CrossRef]
  11. Priya, P.; Sabarmathi, A. Caputo Fractal Fractional Order Derivative of Soil Pollution Model Due to Industrial and Agrochemical. Int. J. Appl. Comput. Math 2022, 8, 250. [Google Scholar] [CrossRef]
  12. Tanriverdi, T.; Baskonus, H.M.; Mahmud, A.A.; Muhamad, K.A. Explicit solution of fractional order atmosphere-soil-land plant carbon cycle system. Ecol. Complex. 2021, 48, 100966. [Google Scholar] [CrossRef]
  13. Vieru, D.; Fetecau, C.; Ahmed, N.; Ali Shah, N. A generalized kinetic model of the advection-dispersion process in a sorbing medium. Math. Model. Nat. Phenom. 2021, 16, 39. [Google Scholar] [CrossRef]
  14. Krasnoschok, M.; Pereverzyev, S.; Siryk, S.; Vasylyeva, N. Determination of the fractional order in quasilinear subdiffusion equations. Fract. Calc. Appl. Anal. 2020, 23, 695–722. [Google Scholar] [CrossRef]
  15. Podlubny, I. Fractional Differential Equations; Academic Press: New York, NY, USA, 1999. [Google Scholar]
  16. Samko, S.; Kilbas, A.; Marichev, O. Fractional Integrals and Derivatives; Gordon and Breach Science Publishers: New York, NY, USA, 1993. [Google Scholar]
  17. Wheatcraft, S.; Meerschaert, M. Fractional conservation of mass. Adv. Water Resour. 2008, 31, 1377–1381. [Google Scholar] [CrossRef]
  18. Gomez-Aguilar, J.; Miranda-Hernandez, M.; Lopez-Lopez, M.; Alvarado-Martinez, V.; Baleanu, D. Modeling and simulation of the fractional space-time diffusion equation. Commun. Nonlinear Sci. Numer. Simul. 2016, 30, 115–127. [Google Scholar] [CrossRef]
  19. Wei, L. Global existence theory and chaos control of fractional differential equations. J. Math. Anal. Appl. 2007, 332, 709–726. [Google Scholar]
  20. Matychyn, I. Analytical Solution of Linear Fractional Systems with Variable Coefficients Involving Riemann–Liouville and Caputo Derivatives. Symmetry 2019, 11, 1366. [Google Scholar] [CrossRef] [Green Version]
  21. Garrappa, R.; Popolizio, M. Computing the matrix Mittag-Leffler function with applications to fractional calculus. J. Sci. Comput. 2018, 77, 129–153. [Google Scholar] [CrossRef] [Green Version]
  22. Li, C.; Zeng, F. The Finite Difference Methods for Fractional Ordinary Differential Equations. Numer. Funct. Anal. Optim. 2013, 34, 149–179. [Google Scholar] [CrossRef]
  23. Garrappa, R. Numerical solution of fractional differential equations: A survey and a software tutorial. Mathematics 2018, 6, 16. [Google Scholar] [CrossRef] [Green Version]
  24. Zayernouri, M.; Karniadakis, G. Exponentially accurate spectral and spectral element methods for fractional ODEs. J. Comput. Phys. 2014, 257, 460–480. [Google Scholar] [CrossRef]
  25. Langlands, T.; Henry, B. The accuracy and stability of an implicit solution method for the fractional diffusion equation. J. Comput. Phys. 2005, 205, 719–736. [Google Scholar] [CrossRef]
  26. Hundsdorfer, W. Unconditional convergence of some Crank-Nicolson LOD methods for initial-boundary value problems. Math. Comput. 1992, 58, 35–53. [Google Scholar]
  27. Rothamsted Experimental Station, Great Britain. Rothamsted: Guide to the Classical Field Experiments; AFRC Institute of Arable Crops Research: Harpenden, UK, 1991. [Google Scholar]
  28. Di Bene, C.; Marchetti, A.; Francaviglia, R.; Farina, R. Soil organic carbon dynamics in typical durum wheat-based crop rotations of Southern Italy. Ital. J. Agron. 2016, 11, 209–216. [Google Scholar] [CrossRef] [Green Version]
  29. Martens, B.; Miralles, D.; Lievens, H.; van der Schalie, R.; de Jeu, R.; Fernández-Prieto, D.; Beck, H.; Dorigo, W.; Verhoest, N. GLEAM v3: Satellite-based land evaporation and root-zone soil moisture. Geosci. Model Dev. 2017, 10, 1903–1925. [Google Scholar] [CrossRef] [Green Version]
  30. Wieder, W.R.; Allison, S.D.; Davidson, E.A.; Georgiou, K.; Hararuk, O.; He, Y.; Hopkins, F.; Luo, Y.; Smith, M.J.; Sulman, B.; et al. Explicitly representing soil microbial processes in Earth system models. Glob. Biogeochem. Cycles 2015, 29, 1782–1800. [Google Scholar] [CrossRef] [Green Version]
Figure 1. From the top-left cell following the left-to-right direction: the numerical solutions c d p m , c r p m , c b i o , c h u m with respect to t t 0 for h = 0.02 in the case of periodic inputs ( ζ ( t , q ) = 1 ).
Figure 1. From the top-left cell following the left-to-right direction: the numerical solutions c d p m , c r p m , c b i o , c h u m with respect to t t 0 for h = 0.02 in the case of periodic inputs ( ζ ( t , q ) = 1 ).
Mathematics 11 01677 g001
Figure 2. From the top-left cell following the left-to-right direction: the numerical solutions c d p m , c r p m , c b i o , c h u m with respect to t t 0 for h = 0.02 in the case of periodic inputs ( ζ ( t , q ) = t q 1 ).
Figure 2. From the top-left cell following the left-to-right direction: the numerical solutions c d p m , c r p m , c b i o , c h u m with respect to t t 0 for h = 0.02 in the case of periodic inputs ( ζ ( t , q ) = t q 1 ).
Mathematics 11 01677 g002
Figure 3. Dynamics of SOC indicator in Scenario 3 for q = 1 and q = 0.98 . On the left: ζ ( t , q ) = t q 1 . On the right: ζ ( t , q ) = t q 1 (default), ζ ( t , q ) = 1 .
Figure 3. Dynamics of SOC indicator in Scenario 3 for q = 1 and q = 0.98 . On the left: ζ ( t , q ) = t q 1 . On the right: ζ ( t , q ) = t q 1 (default), ζ ( t , q ) = 1 .
Mathematics 11 01677 g003
Figure 4. Locations of points with SOC data in Poland and Ukraine.
Figure 4. Locations of points with SOC data in Poland and Ukraine.
Mathematics 11 01677 g004
Figure 5. The values of SOC obtained for the integer-order and the fractional-order models (“corr = 0” means ζ ( t , q ) = 1 , “corr = 1” means ζ ( t , q ) = t q 1 ): left for the point “Sajivka”, right for the point “Chrusciel”.
Figure 5. The values of SOC obtained for the integer-order and the fractional-order models (“corr = 0” means ζ ( t , q ) = 1 , “corr = 1” means ζ ( t , q ) = t q 1 ): left for the point “Sajivka”, right for the point “Chrusciel”.
Mathematics 11 01677 g005
Table 1. Average absolute errors for the problem with a constructed analytical solution.
Table 1. Average absolute errors for the problem with a constructed analytical solution.
h q = 0.99 0.80.60.40.20.01
0.080.001270.0007170.0005770.0004130.0002660.000176
0.040.0006340.0003560.0002880.0002060.000133 8.82 · 10 5
0.020.0003150.0001770.0001440.000103 6.66 · 10 5 4.41 · 10 5
0.010.000157 8.82 · 10 5 7.21 · 10 5 5.16 · 10 5 3.33 · 10 5 2.20 · 10 5
Table 2. Hoosfield spring barley experiment. Values of ρ per month, in crop and fallow years, estimated from clay content of 23.4%, water-related data in [3], and soil cover factors.
Table 2. Hoosfield spring barley experiment. Values of ρ per month, in crop and fallow years, estimated from clay content of 23.4%, water-related data in [3], and soil cover factors.
Months123456789101112
ρ ( t n )
Crop years
0.35610.37230.50680.44710.74730.77790.24910.41510.65701.12770.60920.4594
ρ ( t n )
Fallow years
0.35610.37230.50680.74511.24541.29650.41510.41510.65701.12770.60920.4594
Table 3. Average absolute errors of modeling the data of the Hoosfield spring barley experiment.
Table 3. Average absolute errors of modeling the data of the Hoosfield spring barley experiment.
ScenarioInteger-Order ModelFractional-Order Model ( q = 0.98 )
ζ ( t , q ) = 1 ζ ( t , q ) = t q 1
16.847.266.98
295.3215976.97
3185192197
Table 4. Average absolute differences in percent of average indicator values and execution times, s, for integer- and fractional-order schemes. C-N(1)—Crank–Nicholson scheme; NS—non-standard scheme; S—standard scheme; F(1)—Crank–Nicholson scheme for the fractional-order model ( q = 0.98 ).
Table 4. Average absolute differences in percent of average indicator values and execution times, s, for integer- and fractional-order schemes. C-N(1)—Crank–Nicholson scheme; NS—non-standard scheme; S—standard scheme; F(1)—Crank–Nicholson scheme for the fractional-order model ( q = 0.98 ).
Comparison betweenSOCDRMDPMBIOHUMTime, s
C-N(1)S1.98211.877.685.090.160.23
C-N(1)NS2.04229.897.715.660.17
SNS0.1825.970.441.110.020.37
NS 0.37
F(1) ( ζ = 1 )NS2.06211.197.815.140.764.9
F(1) ( ζ = t q 1 )NS8.17232.4918.2216.977.66
Table 5. Summary of the datasets from Ukraine.
Table 5. Summary of the datasets from Ukraine.
PointClay Content, %Humified Layer Depth, cmMean Temperature, °C P ( t 0 ) Initial SOC, t C / ha
Sajivka55.04010.571.6296138.8
Chrusciel11.83308.461.1161687.36
Table 6. Statistical model indicators RMSE and EF, evaluated for the comparison of the fractional-order and RothC models in the Sajivka and Chrusciel sites.
Table 6. Statistical model indicators RMSE and EF, evaluated for the comparison of the fractional-order and RothC models in the Sajivka and Chrusciel sites.
SajivkaChrusiel
q = 0.97 , c o r r = 1 q = 0.97 , c o r r = 0 q = 1 q = 0.99 , c o r r = 1 q = 0.99 , c o r r = 0 q = 1
EF0.8511730.7983650.8086090.8364870.9222870.900954
RMSE1.2018891.3989681.3629673.0442742.0987152.369329
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Bohaienko, V.; Diele, F.; Marangi, C.; Tamborrino, C.; Aleksandrowicz, S.; Woźniak, E. A Novel Fractional-Order RothC Model. Mathematics 2023, 11, 1677. https://doi.org/10.3390/math11071677

AMA Style

Bohaienko V, Diele F, Marangi C, Tamborrino C, Aleksandrowicz S, Woźniak E. A Novel Fractional-Order RothC Model. Mathematics. 2023; 11(7):1677. https://doi.org/10.3390/math11071677

Chicago/Turabian Style

Bohaienko, Vsevolod, Fasma Diele, Carmela Marangi, Cristiano Tamborrino, Sebastian Aleksandrowicz, and Edyta Woźniak. 2023. "A Novel Fractional-Order RothC Model" Mathematics 11, no. 7: 1677. https://doi.org/10.3390/math11071677

APA Style

Bohaienko, V., Diele, F., Marangi, C., Tamborrino, C., Aleksandrowicz, S., & Woźniak, E. (2023). A Novel Fractional-Order RothC Model. Mathematics, 11(7), 1677. https://doi.org/10.3390/math11071677

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop