Simulation of a time dependent advection-reaction-diffusion problem using operator splitting and discontinuous Galerkin methods with application to plant root growth

—Motivated by the need of developing numerical tools for the simulation of plant root growth, this article deals with the numerical resolution of the C-Root model. This model describes the dynamics of plant root apices in the soil and it consists in a time dependent advection-reaction-diffusion equation whose unique unknown is the density of apices. The work is focused on the implementation and validation of a suitable numerical method for the resolution of the C-Root model on unstructured meshes. The model is solved using Discontinuous Galerkin (DG) ﬁnite elements combined with an operator splitting technique. After a brief presentation of the numerical method, the implementation of the algorithm is validated in a simple test case, for which an analytic expression of the solution is known. Then, the issue of the positivity preservation is discussed. Finally, the DG-splitting algorithm is applied to a more realistic root system and the results are discussed.


I. INTRODUCTION
The article is devoted to the numerical modeling of plant root growth. This work has been originally motivated by the need of developing numerical tools for the simulation of plant growth dynamics. Due to the difficulty of doing nondestructive observations of the underground part of plants (that allow to do long term studies of the dynamics of tree roots for example), mathematical models are achieving an essential role. Several theoretical and numerical challenges arise in the field of the simulation of the dynamics of plant roots [48], [47], [38], [2], [39]. The mathematical description of plant root is not trivial, due to the presence of many interactions arising in the rhizosphere and also due to the diversity of plant root types. Mathematical models based on the use of partial differential equations are useful tools to simulate the evolution of root densities in space and time [43], [44], [44], [45], [46], [41], [40], [1]. This formalism facilitates the coupling with physical models such as water and nutrient transports [42], [43], [44], [41], [49]. And the computational time for the simulation of such models is not dependent on the number of roots which is useful for applications at large scale. The C-Root model [1] is a generic model of the dynamics of root density growth. This model takes only one unknown which is related to root densities such as the density of apices, root length density or biomass density. It has only three parameters. The model is said to be generic in the sense that it can apply to a wide variety of root system types. The model consists in a single time-dependent advection-reaction-diffusion equation, and one of the challenge is to numerically solve the equation. In [1] and [2] the authors solved the problem with the finite difference method on Cartesian mesh grids combined with an operator splitting technique. Unfortunately, Cartesian mesh grids do not allow easily to mesh complex soil geometries. From the theoretical and computational point of view, Cartesian grids also lead to difficulties for a rigorous study and validation of the model. That is why this article focus on the development and implementation of a suitable numerical method for the resolution of the C-Root model on triangular mesh grids, that allow to mesh complex geometries. However, one of the main difficulties in the C-root model is that the advection and diffusion terms are not always of the same order of magnitude. It depends on the phase of the root system development [2]. As a result, the properties of the equation may vary along the simulation: it can be either close to a hyperbolic problem or close to a parabolic problem.
In a previous work [3], the use of the Discontinuous Galerkin method has been implemented and validated. Indeed, the usual choice of the classical Lagrange finite element method suffers from a lack of stability when the advection term is dominant [4]. For this reason, we implemented a discontinuous Galerkin (DG) method for both the advection and diffusion terms. All the three operators where solved simultaneously using the same time approximation scheme (θ-scheme).
However, as explained in [6], for multibiophysic problems it is not efficient to use the same numerical scheme for the different operators of the system. For example, we may want to use the Euler explicit scheme for the advection term and an Euler implicit scheme for the diffusion. The operator splitting technique [7], [8] is a well known alternative for the resolution of equations having a multi-biophysic behaviour that allows the use of different time schemes for each operator of the equation. The idea of the splitting technique is to split the problem into smaller and simpler parts of the problem so that each part can be solved by an efficient and suitable time scheme. This methods has been used for a wide range of applications dealing with the advection-reaction-diffusion equation [9]. Operator splitting techniques have been extensively used in combination with finite difference methods [10], [2], finite volume methods [11], [12] but also with Continuous Galerkin methods [13], [14], [15], [16], [17]. To the best of my knowledge, only very few articles deal with the use of the operator splitting technique in combination with the discontinuous Galerkin approximation [18], [19], [20], [21]. In this paper, we present a new application of the operator splitting technique combined with discontinuous finite elements.
The paper is structured as follows. In section II, the C-root growth model [1], [2], [3] is briefly described. An analysis is also provided, where I showed the existence and uniqueness of a positive real solution. In section III, the splitting operator technique is introduced and applied to the C-Root model, combined with the use of discontinuous Galerkin approximations. In section IV, the algoritm is implemented and validated using a simple test case for which an analytic expression of the solution is known. As an application, I provide simulations of the development of eucalyptus roots in section V. Finally, the paper ends with a conclusion and further improvements.

A. Modelling root growth with PDE: the C-Root model
The C-Root model [1] was developed to simulate the growth of dense root networks, usually composed of fine roots, with negligible secondary thickening. As presented in [1], the unknown variable u is the number of apices per unit volume, but it can also stand for the density of fine root biomass. The soil is considered as a subdomain of R d (with d = 1, 2 or 3). It is assumed that Ω has smooth boundaries (Lipschitz boundaries) denoted ∂Ω. The C-Root model combines advection, diffusion and reaction, which aggregate the main biological processes involved in root growth, such as primary growth, ramification and root death. The reaction operator gives the quantity of apices (or root biomass) produced in time, whereas advection and diffusion operators spatially distribute the whole apices (or biomass) in the domain.
The reaction operator describes the evolution in time of the root biomass in a given domain. In the C-Root model it is a linear term characterized by the scalar parameter ρ which is the growth rate of the root system. The diffusion corresponds to the spread of the root biomass over space. It is described by the parameter σ which is a d × d matrix that characterizes the growth of the root biomass in any direction exploiting free space in the soil. The advection corresponds to the displacement of the root biomass in a direction and velocity given by v which is a vector in R d .
On the boundaries of Ω, what happens for the quantity being transported is different depending if the growth makes the roots to come inside Ω or to go outside of Ω. If v is going inside Ω (at the inlet boundary) the root biomass u will enter the domain and increase. On edges where v is going out of the domain (outlet boundary) the root biomass u is going to be pushed out of Ω. Since this phenomena is oriented (causality) and the behaviour of the solution is different on inlet and outlet boundaries, we need to specify in the model these parts of the boundaries. Mathematically, it is required to define the inlet boundary with respect to v as The outlet boundary Ω + is given by ∂Ω + = ∂Ω\∂Ω − . The dynamics of the root system is studied between the time t 0 and t 1 with 0 ≤ t 0 < t 1 .
The problem reads as follow: find u such that where g ∈ L 2 (∂Ω) and g in ∈ L 2 (∂Ω − ) are given. And u 0 is the given initial solution.
Problem (2) is known as the time dependent advection-reaction-diffusion problem and belongs to the class of parabolic partial differential equations. This equation is a model problem that often occurs in fluid mechanics but also in many other applications in life sciences (see for instance [22], [23], [24]).
Depending on the boundary conditions, the problem has different meanings. To simplify the presentation we only consider the Neumann boundary condition combined with an inlet boundary condition at the inlet of the domain. The Neumann condition specifies the value of the normal derivative of the solution at the boundary of the domain. The inlet boundary condition specifies the quantity of u convected by v that enters in the domain.

B. The weak problem
Since the goal is to solve the problem on unstructured meshes, the spatial operators are approximated using finite element methods. Within this framework, it is classical to write the problem in a variational form. Let us first introduce some functional spaces [50].

C. The positivity preserving property of the solution
In the framework of our applications to the simulation of root biomass densities one of the crucial property of the problem is the preservation of the positity of the solution along time. For a positive initial solution u 0 , the solution of (3) stays positive.
Proof: I follow [25]. See also [26], [27]. We consider the function u − defined by Let us note that That is we have We verify that u − is an admissible test function in W . Using the following obvious equations By adding the same quantity on both sides of the equation we get Since u satisfy (3) we have One can notice that d t u − , u − + d t u, u − = 0. Then we have with g(t, x) ≥ 0 a.e. in ]t 0 , t 1 [×∂Ω. Now from the coercivity of the bilinear form a we obtain where c is a strictly positive constant. The estimate is then By the Gronwall lemma we have that Since c > 0 and t ≥ t 0 ≥ 0 , we have that e −2ct ≤ 1 , so we obtain Since

A. The operator splitting technique
Here we focus on the implementation of the operator splitting technique. The time inter- At each iteration step we solve the following problems . The bilinear forms a A (t, u, v), a D (t, u, v) and a R (t, u, v) are respectively given by (6), (7) and (8). And (t, v) is the linear form (4). If the operators are commutative, then the splitting error vanishes. Otherwise, if the operators are not commutative, then the splitting error does not vanish and a second order splitting would be required (see [6]). In the following, I present the different schemes related to each operator.

B. The advection step: DG upwind scheme
The advection step consists in solving the following transport problem : Find u such that ∀v ∈ H, for a.e t ∈]t n , t n+1 [, where a A (t, u, v) is the bilinear form (6). For the space approximation of this problem, we implemented the DG upwind method presented below.
Let T h be a regular family of decomposition in triangles of the domain Ω such that The h subscript in T h denotes the size of the mesh cells and it is defined by are the sets of edges belonging to ∂Ω − and ∂Ω + respectively. And E i h is the set of interior edges. Let us consider an element of E i h . We denote by T + and T − the two mesh elements sharing the edge e so that e = ∂T + ∩ ∂T − where the minus and plus superscripts depend on the direction of the advection vector. By convention we suppose that v goes from is the outward normal vector of e in T + (resp. T − ). When it is not necessary to distinguish the orientation of the normal vectors n + e and n − e we denote by n the unitary normal of e.
Let us consider the advection problem on each element K i of the domain : for all K i , i = 1, N we look for u the solution of the equation (12) defined on K i . Similarly to the problem defined on all the domain Ω, we look for a solution u that is in L 2 (K i ) and such that ∇u is in L 2 (K i ) for all K i in T h . Let us introduce the following broken Sobolev space: with v e the trace of v along e. The mean value of u on e is defined by Besides for edges on the boundaries we take h . Let us denote by X the functional space defined such that This space is a Hilbert space equipped with the norm The DG variational formulation of the advection step written on the broken Sobolev space takes the following form: Find u in X such that for a.e where the form a up h (t; u, v) is the approximation of the advection term. It consists in the upwind formulation of the DG method [28]. It reads: The approximated linear form of the the right hand side reads The DG-formulation (14) is consistent and stable, see for example [32]. The discontinuous Galerkin method consists in searching the solution in the approximation space X h defined such that where W k h is given by Let us note that the functions of W k h can be discontinuous from one element of the mesh to the other. Let us note that W k h is embedded in In this basis the approximated solution takes the form: where the ξ i (t) are the degrees of freedom. Let us define X the vector of degrees of freedom: The approximated problem then reduces to find and L up h (t) is the vector of size n defined such that The problem reduces to a linear system of ordinary differential equations. The time approximation is based on a finite difference scheme.
At each iteration step we solve the following problem: Find X N +1 ∈ R n such that where θ is a real parameter taken in [0, 1]. For θ = 0, we have the explicit Euler schema. For θ = 1, it is the implicit Euler schema. For θ = 1/2, it is the Crank-Nicolson schema.

C. The diffusion step
The diffusion step consists in solving the following problem : Find u such that ∀v ∈ H, for a.e t ∈]t n , t n+1 [, where a D (t; u, v) is the bilinear form (7) and (t; v) is the linear form (4). In the setting introduced before, the DG variational formulation of the diffusion step written in the broken Sobolev space takes the following form: The form a ip h (t; u, v) is the approximation of the diffusion term. It consists in the interior penalty formulation (IP) that reads where η is a positive penalization factor. The linear This formulation was introduced in [31] and is known as the non-symmetric interior penalty (NSIP) formulation, see [30], [32]. In matrix form the problem reduces to find where M is defined by (15) and A ip is defined such that Similarly to the advection step, the time approximation of the problem is based on a finite difference scheme of the form (17).

D. The reaction step
The reaction step consists in solving the following problem : Find u such that ∀v ∈ H, for a.e.
where a R (t; u, v) is the bilinear form (8). This problem takes the following matrix form find where we recall that ρ is a constant real parameter. This problem can be solved by an exact scheme (a kind of schemes that provide exact solutions, i.e. a solution equal to the analytical solution). At each iteration we find X N +1 such that ). This scheme is unconditionally stable, meaning that we can choose the time step independently from the space step. It is also positively stable, meaning that if X N ≥ 0 so is X N +1 .

IV. VALIDATION OF THE SPLITTING ALGORITHM WITH A SIMPLE TEST CASE
Problem (3) has been already solved using discontinuous Galerkin elements (DG) [3]. Advection and Diffusion operators were solved simultaneously using the Crank-Nicolson scheme providing stable results. However, even for simple test cases some simulations did not always provide positive numerical solutions. One reason is that the same time approximation scheme is not necessarily suitable for both the advection and for the diffusion. That is why a new operator splitting algorithm has been implemented with a different time scheme for each operator.
The goal of this section is to validate the implementation of the code. To this end I compare the convergence of the approximation with and without the splitting technique. I briefly explore the question of the positivity of the approximated solution.

A. Description of the simple test-case
First let me introduce a simplified test-case for the validation of the splitting algorithm. Set L > 0, and Ω =] − L; L[ 2 . Let v = (v 1 , v 2 ) ∈ R 2 and d ∈ R be a constant and 0 ≤ t 0 ≤ t 1 . Find u such that ∂u ∂t The initial condition and the boundary condition are chosen such that the solution of problem (18) is explicitly given by ∀(x, y, t) ∈ Ω×]t 0 , t 1 [ with κ(x, y, t) where c 0 > 0, a > 0, x 0 and y 0 are real parameters and v 1 and v 2 are the two components of v. Notice that u(x, y, t) > 0 for all (x, y, t) in Ω×]t 0 , t 1 [.

B. Numerical validation and convergence
To validate the implementation of the splitting technique, I ran the previous test case with different mesh sizes and time steps and I computed the global L 2 -errors such that where t k = t 0 + kδt, with k ∈ N + * and t N = t 1 . The flexibility of the splitting technique allows to choose different time schemes for each operator. I consider a θ-scheme with θ = 0 (explicit Euler), θ = 1 (implicit Euler), and θ = 1 2 (Crank-Nicolson) for both the advection step and the diffusion step, and I consider an exact scheme for the reaction step. For the simulations I took the parameters such that v = (0.1, 0) T , σ = 0.01 and ρ = −1. The triangular meshes used for the simulations are identified by h which is the size of the biggest triangle of the mesh.  L = 1/2, the simulations are performed between t 0 = 0 and t 1 = 1 for different values of the time step δt. Fig. 1, page 9, shows the solution at t = t 0 and t = t 1 . The code is implemented in Fortran 90 and it is run under a 64-bit Linux operating system on a 8-core processor Intel R Core TM i7-7820HQ at a frequency of 2.9GHz and with 32 GB of RAM. The sparse matrices resulting from the finite element approximation are inverted using a solver provided by the library MUMPS [51], [52].
According to Fig. 2, page 10, all the three temporal schemes provide results with approximately the same level of accuracy with a spatial convergence rate of 2 computed with the global L 2 - Fig. 2. Convergence of the solution with respect to the mesh size: plot of the total L 2 -error computed between t = 0 and t = 1 with and without the splitting technique for the explicit Euler scheme (θ = 0), the implicit Euler scheme (θ = 1) and the Crank-Nicolson scheme (θ = 1/2) for δt = 5 × 10 −5 . norm. The same order of convergence is obtained when the problem is solved without the splitting technique. Figure 3 on page 10 shows that the Crank-Nicolson scheme (θ = 1/2) converges in δt 2 while the Euler Implicit scheme (θ = 1) converges in δt, with and without the splitting technique. The convergence rate in time has to be computed with a really refined mesh grid (here h ≈ 4.1 × 10 −3 ). Fig. 4. Validation of the test case: plot of the CPU time against the mesh size (h) for the computations performed with a processor Intel R Core TM i7-7820HQ at 2.9 GHz and RAM 32 GB, between t = 0 and t = 1 with and without the splitting technique for the explicit Euler scheme (θ = 0), the implicit Euler scheme (θ = 1) and the Crank-Nicolson scheme (θ = 1/2) for δt = 5 × 10 −5 . Fig. 5. Validation of the test case: plot of the CPU time against the time step (δt) for the computations performed with a processor Intel R Core TM i7-7820HQ at 2.9 GHz and RAM 32 GB, between t = 0 and t = 1 with and without the splitting technique for the explicit Euler scheme (θ = 0), the implicit Euler scheme (θ = 1) and the Crank-Nicolson scheme (θ = 1/2) for h ≈ 4.1 × 10 −3 and δt ranging from 5 × 10 −1 to 5 × 10 −4 . Note that the computations performed here with θ = 0 gave unstable results. It results an additional cost in term of CPU time, since it behaves like 1/h 2 , as shown on figure 4 page 10. For bigger values of h the plot of the errors gave convergence order in time less than 1 and 2 for the implicit Euler scheme and the Crank-Nicolson scheme respectively. As expected, the explicit Euler scheme is conditionally stable, such that, when the CFL condition is fulfilled, the computational time becomes prohibitive. Indeed, it behaves like 1/δt, as shown on figure 5. For instance, the computation with h ≈ 4.1 × 10 −3 and δt = 10 −5 takes more than 30 hours with the device specified above. That is why, in the rest of the paper, we will only focus on implicit Euler and Crank-Nicolson schemes. However I present here additional computations performed with a bigger mesh size (h ≈ 1.6 × 10 −2 ) and smaller time steps chosen such that the CFL condition is fulfilled. The global L 2 -errors and the CPU time are shown on table II. Clearly, the mesh is not fine enough to recover the convergence order in δt, indeed decreasing the time step results only in an increase of the computational time but not in a significant decrease of the errors.
With the Crank-Nicolson scheme (θ = 1/2), for a given mesh size, if δt is too big, there is no threshold of positivity in t k ∈]t 0 , t 1 ] and the computed solution is not non-negative all along the simulation. For θ = 1/2 and θ = 1, still with a given mesh size, if δt is too small, the simulations showed that there is a threshold of positivity t + such that the approximated solution becomes non-negative for t k ≥ t + . The thresholds of positivity slightly depend on the time step and tend to increase when the time step δt is decreased. The computations clearly showed that the positivity thresholds diminish with the mesh size h (see for example [36]).
Altogether, the positivity of the approximated 0.2104 0.0877 0.0221 Implicit Euler scheme (θ = 1) solution is obtained at the expense of the computational cost, but for a given mesh size h computations performed with too small time step can also lead to a loss of positivity for small t k . In [36] (and references therein), Thomée showed that threshold values of t k > 0 may exist such that X(t) > 0 when t > t k .
At this stage, one may wonder how each term of the splitting behaves in terms of positivity preservation. The reaction term is approximated using an exact scheme, so obviously the positivity of the solution is preserved. What about the diffusion and the advection term ?
3) Positivity of the pure advection problem: Here I set σ = 0 and ρ = 0, while keeping all others parameters to the same values as in the first test. Table IX shows that none of the computations performed gave a non negative solutions, even though the minimum value of the dof can be really close to zero for small mesh sizes. Besides, I did Unfortunately, this threshold of negativity is really small compared to the ending time of the computation (t 1 = 1), while the computational time was reaching more than 14 hours (Intel R Core TM i7-7820HQ at 2.9 GHz, RAM 32 GB) for both the Crank-Nicolson and the implicit Euler schemes. In fact it is well known that for the advection term the solution can be polluted by overshoot and undershoot oscillations near a discontinuity or a sharp layer, see [34], [33], [35], [30]. For low order accurate spacial approximations one can prove the positivity preserving property of the scheme [33]. But for high order schemes slopes limiters are often required to guarantee the positivity of the approximated solution. When slope limiters are used, explicit time schemes seem to be suitable for the advection [6]. However, in the next section we will only privilege a numerical scheme that is unconditionally stable, i.e. the Crank-Nicolson scheme, that is a two-order scheme.

SYSTEM GROWTH
In this section, I apply the previous DG-splitting approach to solve numerically the C-Root model. First, I detail the parameters used for the simulations, then, I present and validate the results of the simulations.

A. The C-Root parameters for Eucalyptus root growth
The parameters and operators' coefficients are chosen based on the previous calibration done in [2]. The diffusion coefficient, σ, is build using the following Gaussian function where r(x, y) = (x − x 0 ) 2 + (y − y 0 ) 2 and (x 0 , y 0 ) ∈ Ω =] − L, L[. The function f α,µ (x, y) depends on two real and positive parameters: α, related to the maximum amplitude of f α,µ , and µ, the distance from (x 0 , y 0 ) to the point where the function f α,µ reaches its maximum.
The diffusion tensor is taken such that σ(x, y) = f αd,µd (x, y) 1 0 0 1 , for all (x, y) ∈ Ω, and α d , µ d ∈ R + are given parameters. The advection vector is taken such that v(x, y) = (0, −v 0 ) T , for all (x, y) ∈ Ω, with v 0 a positive constant. The reaction term is constant in space and splited into two contributions: β r and µ r , the branching and mortality rates, respectively. That is The branching rate, β r , is estimated from biological knowledge: it is equal to zero before 9 months and equal to 1/3 after, since no roots die before 9 month. However, for the following simulations we will not distinguish the contribution of β r and µ r , so that the reaction term will only be described by the parameter ρ.

B. Some simulations
For the simulation the initial solution is chosen equal to the following function: with A = 2 · 10 −4 and b = 1. The parameters' values µ r , α d , µ d are estimated using the code described in [2]. I run the simulations from t 0 = 1 to t 1 = 24 months, with L = 13. The simulations are performed for different values of the mesh size. Fig. 6, page 15, shows the solution computed at four different stages of the root system development. One can notice the diffusion of the apices in the soil and also the transport of the apices from the top to the bottom of the soil layer. Since there is no analytic solution, the convergence of the computation is evaluated by measuring the L 2errors with respect to the approximated solution computed with the finest mesh (h ≈ 8.97 × 10 −2 ) and with δt = 10 −3 . The curves of the errors against the mesh size are plotted on figure 7 and clearly show that the DG-splitting algorithm converges with a convergence rate of almost two. However, one can note that the mesh sizes and the time steps chosen for the simulations presented here might not be small enough. The positivity of the solution is not preserved at all times and the full convergence might not be acheived. Unfortunatly, refining the mesh sizes and the time steps can lead to prohibitive computational time as shown on table X. On top of that simulation of root system growth can last for a long period of time, particularly for trees. Finally, this application shows promising results for future simulations of h (≈) δt = 10 −1 δt = 10 −2 δt = 10 −3 1. 44   the root system growth, provided that the computational cost is not limiting. Further simulations requiring much more computational power has to be done to check if the convergence is acheived. This application also point out the difficulties related to the rigorous simulation validation in realistic test-cases of root system growth.

VI. CONCLUSION
In this work, a discontinuous Galerkin approximation method based on unstructured mesh combined with operator splitting has been described, implemented and tested, to solve an advectiondiffusion-reaction equation used to model the growth of root systems. The code has been validated in a simple test case for which an analytic expression of the solution is known. The computations showed that the method convergences with a convergence rate of two in space with P 1 -finite elements. A convergence rate of one and two in time were obtained for respectively the implicit Euler scheme and the Crank-Nicolson scheme both with and without the splitting technique. The computations of those convergence rates required the use of fine mesh grids. For the explicit Euler scheme, such fine mesh computations were not performed since they require really small time steps to fulfill the CFL condition, resulting in huge additional computational cost. Indeed the computational time of the DG-splitting algorithm behaves like 1/δt and 1/h 2 where δt and h are respectively the time step and the mesh size.
Similarly, the positivity of the approximated solution is obtained at the expense of the computational time since it requires meshes of small size and small time steps. In fact, there is a CFLlike condition for positivity that has to be fulfilled to guarantee the positivity of the approximated solution. But for a given mesh size computations performed with too small time step can also lead to a loss of positivity at the beginning of the computation [36]. In that cases, the computations showed that there is a positivity threshold in time after which the solution becomes positive. This positivity threshold clearly appeared to diminish with the mesh size. This behavior is specific to the diffusion term. For the advection term, the computations also showed that the positivity of the solution can be preserved, but only at the beginning of the simulation and it required a really small mesh size and time step leading to huge computational time. Further studies in terms of numerical analysis has to be done in that direction.
I also performed a more realistic simulation of root system growth. The computations showed that the algorithm converged but additional simulations with smaller time steps and mesh sizes might be performed to recover the full convergence order and positivity. Validation of the computation, but above all the computational time appeared to be the major limitations of the root growth simulation based on the C-Root model, particularly when it comes to deal with trees for which the life span is rather a long period of time. Further improvements on the numerical method has to be done so that the scheme preserves the positivity of the approximated solution under acceptable CFL conditions in terms of computational time. However, our work shows promising results for the simulation of the C-Root model which appears to be an appropriate methodology for future improvements, like rootsoil coupling or nonlinear terms arising to handle competition phenomena.

ACKNOWLEDGMENT
The author would like to thank Y. Dumont (CIRAD, University of Pretoria) for discussions and valuable comments about the numerical schemes.