Modeling and numerical simulation of a vegetation cover

: The aim of this work is to introduce a mathematical model representing the evolution of 1 the temperature in a vegetation cover and the ground underneath it. Vegetation, and its interaction 2 with soil, plays a very important role in the protection of soil surface from the action of sun 3 and precipitations. A reduction in the vegetated mass increase the risk of desertiﬁcation, soil 4 erosion or surface runoffs which which can give rise to soil loss and sediment retention. These 5 processes can favour climate change and global warming, which are major concerns nowadays. The mathematical model presented takes into account the main processes involved in vegetation cover and the interaction with the soil, among which, we can mention the Leaf Area Index, which is a 8 dimensionless quantity deﬁned as the one-sided green leaf area per unit ground surface area, or 9 albedo and co-albedo which are clearly inﬂuenced by the vegetation. It is also considered a nonlinear 10 heat capacity in the soil which incorporates the latent heat of fusion, when the phase change takes place. The numerical technique used to solve the mathematical model is based on a ﬁnite volume 12 scheme with Weighted Essentially Non Oscillatory technique for spatial reconstruction and the third 13 order Runge-Kutta Total Variation Diminishing numerical scheme is used for time integration. Some 14 numerical examples are solved to obtain the distribution of temperature both in the vegetation cover 15 and the soil. 16 time integration and a seventh order WENO approach for spatial 241 reconstruction. These techniques are very efﬁcient for solving this type of problems. 242

environments is given in [8]. 48 Soil temperature is relevant on soil properties and plant growth. In [2], and references therein, 49 interesting aspects on the influence of soil temperature are described in detail, such as the influence of 50 soil temperature on bioactivity, soil micro and macro-organisms and decomposition of organic matter.

51
In addition, the authors describe other characteristics where soil temperature is very important. For 52 instance, related to chemical properties of the soil, important phenomena are: cation exchange capacity 53 (CEC), that is the total capacity of a soil to hold exchangeable cations, available phosphorous and soil 54 pH. Concerning the effect of the temperature on soil physical properties (see [2]) a dehydration process 55 may happen, affecting soil structure. Furthermore thermal transformation of iron and aluminium 56 oxides may produce changes in clay particles. Also, an increase in soil temperature can reduce water 57 viscosity allowing it to percolate through the soil. Moreover, higher evaporation rates may arise, 58 restricting the movement of water through the soil. Concerning soil aeration, carbon dioxide amount 59 in the soil air can be affected by an increase in temperature. Finally, according to [2], soil temperature 60 also has an influence on plant growth, such as nutrient and water uptake by plants or roots growth.

61
In this work we introduce and solve numerically a coupled model vegetation cover-soil, in order 62 to analyze the thermal effects involved and the distribution of temperatures both in the vegetation 63 cover and the soil. The model is based on ideas from a green roof model, proposed in [9,10], and on a 64 climate model, see [11,12]. The numerical resolution is performed by means of a finite volume (FV) The diffusion equation representing the evolution of the temperature in the vegetation cover is given as where T v = T v (x, t) is the temperature in the vegetation cover, T s = T s (x, z, t) is the temperature 77 of the soil, K H0 is the heat conductivity in the vegetation cover, K V is the conductivity of the soil in 78 vertical direction, C v the heat capacity of the vegetation cover, x is the spatial coordinate, t is time. The 79 term K V ∂T s ∂n which considers the normal derivative of the soil temperature, represents the influence 80 of the soil temperature onto the upper boundary where the vegetation is located. This term has been 81 introduced to account for conduction phenomena from the ground to the vegetation cover.

82
This kind of dynamic and diffusive boundary condition appears also in many contexts, among 83 them, we can mention a model arising in Biology ([13,14]) for the analysis of inflammatory process in 84 the artery wall and a climate model including surface-ocean interaction ([12]).

85
In (1) the source term is given by where H v represents the sensible heat flux, which stands for the energy loss by the vegetation cover by heat transfer to the surrounding air and, according to [9,15] being e 0 the windless exchange coefficient which depends just on temperature difference, LAI is the 86 leaf area index, which is a dimensionless quantity which is defined as the one-sided green leaf area per 87 unit ground surface area. The density of the foliage is represented by ρ a f and the units are kgm −3 . C a 88 is the specific heat of air at constant pressure, C f stands for the bulk transfer coefficient, given by the The latent heat flux is given by where l v is the latent heat of vaporization which, according to , is given by As for the equation for the evolution of temperature in 2D coordinates, it is given by where K H and K V are the heat conductivities of the soil in x and z directions respectively. The coordinate x is along the surface and z represents depth. The equation (6) considers the water phase change, which is described by γ(T s ). We notice that, in this model, heat capacity is not a constant, but a function of temperature. In this work it is introduced a similar expression to the one used in global climate models for latent heat of fusion, see [13], which is also similar to an expression utilised in green roof modelling for heat capacity, see [9], reading whereT s = T s − 273. When solving the model (6) the heat capacity γ(T s ) is obtained. According to the 99 expression given in (7), the temperature of the soil, T s , can be computed applying a solver of nonlinear conducted in the next sections, but it is included in this one for the sake of the explanation of nonlinear 113 solver. We recall that this is an automatic process, in the sense that regula falsi or bisection technique 114 only intervenes when Newton-Raphson does not converges.

116
In order to consider that the interaction between the vegetation cover and the soil takes place within the mixed layer, the right hand side of equation (6) is given by where, as aforementioned, Dmix is the depth of the mixed layer, where the influence of the vegetation on the soil takes place. In (6) f s,1 and f s,2 are defined making use of Stefan-Boltzmann's law for radiation and are expressed as In the expression (9) it is introduced E = R as + s I ir where R as is the radiation absorbed by the soil 117 given by R as = Q|(sin( πt 12 ) + cos( πt 12 )|β s where β s is the coalbedo of the soil. In (8) the following quantities are introduced and where is the resistance 119 to moisture exchange and r s is the stomatal resistance. 120 We also consider the Bulk Richardson number given by . According to this value it is defined the so-called atmospheric stability factor as We where κ is the universal von Kármán constant, which is frequently used in turbulence modelling, as 122 it happens in boundary-layer meteorology to account for fluxes of momentum, heat and moisture 123 from the atmosphere to the land surface. We also have the values The values taken in this work for the different parameters involved in this model for the numerical 126 simulations are displayed in Table 1.

128
In this section it is briefly described the numerical technique followed, which is developed in the 129 finite volume framework. The main reasons behind using finite volume schemes are: 130 • They are conservative by construction in the sense that physical variables are conserved.

131
• Due to their discontinuity nature, since it makes use of cell averages of the solution, 132 discontinuities are well resolved.

133
In this work a semi-discrete finite volume scheme is built for the vegetation cover, equation (1), and for the soil, equation (6). Regarding the vegetation cover (1), the 1D domain [−l, l] is discretized in N x control volumes. The equation (1) being the spatial cell average of the temperature T v (x, t) in the control volume S i at time t, is the right interface numerical flux at time t, and is the spatial average of the source term F v in the control volume S i at the particular time t.

134
Considering the soil equation (6) , integration on the control volume gives where is the cell average of the soil temperature inside the control volume I ij , while the value F i+ 1 2 ,j is the right intercell numerical flux in x−direction and G i,j+ 1 2 is the intercell numerical flux in z−direction at time t, and are the spatial average of physical fluxes over cell faces at time t and is the spatial average of the source term F s (T s ) over the control volume I ij .

136
The numerical solution may be evolved in time by means of different procedures. In this work it is performed a third order TVD Runge-Kutta method, as introduced in [18] and described in [12] for global climate modelling. The expressions of this method read where η refers to the temperature of the vegetation cover, T v or to γ(T s ) for the evolution of the 137 temperature in the soil. Furthermore, the operator Λ(·) is l i (·) in (13) for the vegetation cover whereas 138 the operator L ij (·) in (17) for the soil part.

139
The main features of this process are:

140
• High-order reconstruction of fluxes at cell interfaces, achieved using Weighted Essentially Non 141 Oscillatory (WENO) reconstruction.

142
• Third-order evolution in time, using a RK3-TVD approach.

143
For each particular time step, we first solve the evolution of the temperature in the 1D vegetation 144 cover, T v , given the equation (1)  In the 2D dimension-by-dimension WENO procedure we consider the control volume I i,j , where cell averages of the solution are computed, and introduce the following set of one-dimensional reconstruction stencils where L and R represent left and right spatial extension of the stencil. The process considered here 168 follows the strategy put forward in different references such as [12,26,27], where three candidate 169 stencils are used for odd order while four possible stencils are taken into account for the even order 170 case. See the aforementioned references for details on this process.

171
We introduce the mapping [ The way to perform, which is described in detail in several references, 173 see, for instance [26], is briefly explained here.

174
We start by carrying out reconstruction in x−direction. For each control volume the reconstruction polynomials are with the basis interpolation functions ψ p , which in this work are Lagrange ones, although different types, such as Legendre ones, could also be used perfectly well. The valuesŵ n,s ij,p are the coefficients of interpolation. Now we get 1 ∆x e x e+1/2 x e−1/2 M ∑ p=0 ψ p (ξ(x))ŵ n,s ij,p dx =ū n ej , ∀I ej ∈ S s,x ij , applying integral conservation. If we carry out now a data-dependent nonlinear combination of the coefficients of the polynomials for each particular stencil it is obtained. Applying now a nonlinear combination, typical of WENO approach, of the coefficients it is achieved with ω s =ω s ∑ kωk , withω s = λ s (σ s + ) ν , and introduced to avoid division by zero. In this work it is used = 10 −20 and the parameter ν is taken as ν = 3. The oscillation indicators are with the oscillation matrix The next step consists of reconstruction in z-direction, performing in a similar way to the previous case, but applying conservation to each particular degree of freedomŵ n ij,p . Then we have Integral conservation in z-direction is carried out for each particular degree of freedom in x-direction, 175 for all the control volumes of the stencil in z-direction, that is where S s,z ij are the control volumes in z−direction.
Eventually it is obtained for the control volume I ij as the final reconstruction polynomial.

179
In the present work we are using a WENO7 approach, where four cells are used for each stencil 180 (r = 4, M = 3) and the developed scheme is seventh order accurate in space. This way to proceed is 181 computationally less expensive than the fully two-dimensional reconstruction. The mathematical model under study is a coupled model formed by a vegetation cover and soil underneath it, describing the evolution of temperature. The main processes and equations are described in Section 2. The system (33) represents the full coupled model.
In the problem (33) C is a constant and it is considered the spatial domain (−l, l) × (−H, 0) as local 184 reference coordinates. The numerical simulation is conducted applying the finite volume numerical 185 scheme RK3TVD with WENO technique for spatial reconstruction as described in Section 3. As 186 previously detailed, the nonlinearity appearing due to the heat capacity, which incorporates the latent heat of fusion, is solved by a Newton-Raphson's numerical scheme combined with regula 188 falsi technique. The spatial mesh is formed by 50 × 50 control volumes. The size of the time step is 189 computed as ∆t = min(CFL × ∆x 2 K H0 , CFL × ∆x 2 K H , CFL × ∆z 2 K V ) where CFL = 0.35 for stability reasons.

190
Concerning the heat conduction we take K H0 = 0.01, K H = 0.03 and K V = 0.03.

191
In Figure   The results displayed represent the spatial distribution of the temperature T s (0, z, 0.5), that is 213 cutting the 2D plots with x = 0.

214
The last example considered takes into account a colder situation than the previous ones, taking 215 C = −40 in the initial condition of the system (33) and, as ambient temperature, T a = 265 + sin(t).

216
Since the temperature ranges from values smaller than 273K to values higher than 273K the effect of 217 the latent heat of fusion, that is present in the heat capacity of the soil. This is clearly visible in the 218 top-left panel of Figure 6, where when the temperature in the soil is around 273K the temperature 219 remains approximately constant since the energy goes to the phase change, instead of to increase 220 the temperature. This effect has also been studied in climate models, such as reported in [11]. In 221 the top-right panel, the effect of latent heat of fusion is not considered. This means that it is used   Figure 4. temperature in the soil (T s ) for output times=0.05, 0.1, 0.5, 1, 3, 5 when the ambient temperature is T a = 290 + sin(t) and C = 6 in the initial condition of (33). The depth is represented by the coordinate z.
γ(T s ) = T s in equation (33) and the plot shows that the variation of temperature is smoother than 223 when latent heat of fusion is taken into consideration.

224
The bottom panel represents a cut with x = 0 of the plots above for a more clear comparison of 225 the cases with and without latent heat of fusion. T a =265+sin(t) T a =273+sin(t) T a =280+sin(t) T a =290+sin(t) T a =293+sin(t) T a =308+sin(t) T a =313+sin(t) Figure 5. Soil temperature distribution with depth (T s (0., z, 0.5)).  Figure 6. Temperature in the soil (T s ) for output time=0.5 when the ambient temperature is T a = 265 + sin(t) and C = −40 in the initial condition of (33). The depth is represented by the coordinate z. Top: temperature in the soil considering latent heat (left panel) and without latent heat (right panel). Bottom: evolution of temperature in the soil (T s ) cutting with x = 0. the function γ(T s ), where T s is the temperature in the soil. Therefore, heat capacity is not a constant, 234 but a function of temperature of the soil. This nonlinearity is resolved by means of a Newton-Raphson's 235 technique combined with a regula falsi one (which acts automatically when Newton-Raphson's scheme 236 fails to converge).

237
The main interactions between the vegetation and the soil take place in a thin region close to the 238 surface, denoted here as mixed layer.

239
The mathematical model is solved using a finite volume numerical scheme, with a third order 240 TVD Runge-Kutta scheme for time integration and a seventh order WENO approach for spatial 241 reconstruction. These techniques are very efficient for solving this type of problems.
The results obtained show that the main temperature variations take place in the vicinity of the 243 surface, namely the mixed layer, whereas as depth increases the temperature tends to a constant value, 244 as expected. This conclusion agrees with the real situation as reported in references on this topic.

245
Further research will include the consideration of water content and moisture in the soil, as well 246 as taking into account different types of plants. In addition, it will be interesting to introduce the effect 247 of precipitations, that may affect the integrity of the vegetation cover.  The following abbreviations are used in this manuscript: