Simulation of Heat and Water Transport on Different Tree Canopies: A Finite Element Approach

: Heat and water transport modeling is a widely explored topic in micro-meteorology, agriculture, and forestry. One of the most popular models is the Simultaneous Heat and Water (SHAW) model, which includes partial differential equations (PDEs) for air-soil temperature and humidity, but with a priori discretized PDE for the foliage temperature in each canopy layer; it is solved using the ﬁnite difference method and the canopy shape is deﬁned as a simple rule of proportionality of total quantities such as the total leaf area index. This work proposes a novel canopy shape characterization based on Weibull distribution, providing a continuous vertical shape function capable of ﬁtting any tree species. This allows formulating a fully continuous SHAW-derived model, which is numerically solved by a ﬁnite element approach of P 1 Lagrange type. For this novel approach, several numerical experiments were carried out to understand how the shape of well distinguishable canopies inﬂuences heat and water transport.


Introduction
Forest and plant canopies have vertically varying profiles of different physical quantities related to complex transport processes. Variables as leaf temperature, water potential, wind velocity, longwave and shortwave radiation, photosynthesis, air and soil temperatures, but also the fluxes between soil-canopy-air have vertical profiles which change in function of the vertical distribution of the canopy biomass-mainly its leaves-(see the review by [1] and references therein). Such variability of physical properties deeply affects evapotranspiration and heat fluxes within the soil-vegetation-atmosphere system (SVAT).
To study complex processes within the canopy, mathematical models with a single homogeneous biomass layer -the so-called "big leaf" models-were used as a first approach. Later models were improved by including soil-vegetation-atmosphere heat and mass fluxes, and finally, multi-layer models arise as a more realistic choice as the vertical canopy shape is considered. However, the multi-layer models sometimes are not necessarily better; thus, model selection is subject to the scope of study [1]. This allows the formulation of different models, single or multi-layered, to study the different phenomena in forest and plants canopies, some of which are harmful to human life quality. Several examples could be enunciated: an evaluation of the radiation that reaches the soil and tree transpiration as means to mitigate the urban heat island (UHI) [2], a study of trees resistance to wind loads, and efficiency evaluation of trees as windbreakers [3], formulation of a multi-layer model to estimate the radiation distribution inside the canopy [4], use of a high-resolution model to study the water stress in a forest mulch layering [5], estimation of heat and water fluxes in the soil layer [6].
One of the well-known models is the Simultaneous Heat and Water (SHAW) model formulated by Flerchinger et al. [7] as a PDEs' model which has been updated several times [8][9][10], improving with it resolution associated with radiation, mass, sensible and latent heat fluxes, and recently [11] to model vertical profiles of temperature and humidity in a given tree canopy.
The SHAW model significantly depends on the leaf area index (LAI, a parameter that quantifies the canopy biomass) to model the leaf-air interchange of sensible heat and humidity. Frequently LAI is estimated merely as the total surface of leaves projected on a horizontal plane, known as total LAI, and holds no information about the vertical biomass distribution through the canopy [12]. The SHAW model requires a known value of LAI in each canopy layer, turning LAI in a priori discretized function; for example, in [11], this requirement is fulfilled with in situ measurements. This mandates a semi-discretized formulation of the SHAW model, leading to a finite difference scheme for its numerical resolution. Therefore, a continuous vertical function of LAI is needed for a continuous formulation of the SHAW model. The comments above lead us to focus on combining a continuous vertical canopy biomass profile with a multi-layer model to better approach vertical profiles of modeled state variables. As antecedents, in [13], the Weibull distribution is used successfully in curve fitting the leaf area density LAD-a vertical function of biomass canopy-of some tree species. On the other hand, in [14] a multi-layer, multispecies, and soil coupled model uses a continuous vertical function of LAD to upscale CO 2 , water, heat, and momentum exchange from the canopy, validating the methodology with data from a tree canopy in Finland.
The accurate results shown in [14] suggest using a similar methodology to study -from a numerical point of view-the influence of the vertical canopy profile on the humidity and heat transport and air, leaf and soil temperatures. Therefore, in this paper, a SHAW-derived model is formulated in a fully continuous form (Section 2.1) using a continuous vertical profile of LAI drawn from the Weibull distribution (Section 2.2). After a detailed description of the initial and boundary conditions (Section 2.3), a combination of numerical techniques is proposed to solve the model (Section 2.4). It combines a finite element approach of P 1 Lagrange type for space discretization (a novelty in this kind of models), a finite difference scheme for time discretization, and a relaxed fix-point algorithm to address the resulting non-linear system. This method is applied in several numerical experiments where the resulting vertical profiles of leaf-air-soil temperatures, and humidity and heat fluxes corresponding to four different-shaped canopies showed notable variations (Section 3). Final remarks, discussions, and future work are presented in Section 4.

The Mathematical Model
In the model an air column with height h a has within a tree canopy with height h c and trunk height h 0 (then 0 < h 0 < h c < h a ). Furthermore, a soil layer with depth h s is considered. Then we define the following sub-domains, the air column Ω a = (0, h a ), the foliage layer of the canopy Ω c = (h 0 , h c ) and the soil layer Ω s = (0, h s ). All of them are such that Ω c ⊂ Ω a and Ω s ∩ Ω a = ∅. Now, we look for the state functions of air temperature T a (z, t), air humidity p v (z, t), leaf temperature T l (z, t), and soil temperature T s (z, t) such that satisfies the following nonlinear system of PDEs The mathematical model (1a)-(1d) is similar to the called Simultaneous Heat and Water (SHAW) model proposed by Flerchinger et al. [7,8] but not include a PDE related to soil humidity and Equation (1c) is continuous instead of discrete. The system incorporates air-canopy interchanges of sensible heat and humidity, absorbed shortwave and longwave net radiations, and assumes a well-mixed air column by turbulent transport. Consequently, the model needs complex empirical formulas, nonlinear coefficients, and correlations to define interchanges rates, resistances, friction velocity, and vertical wind profile. All of them are written down below as they were formulated in [7,8,11,15].
The turbulent transport coefficient K e is given by the widely use K-theory being k the Von Karman constant, z H the canopy height, d ≈ 0.77h c the zero-level displacement, φ H is the dimensionless heat stability factor, and u * the friction velocity given by where u re f is the reference wind speed at a given height z re f . This formulations of K e and u * are simple and commonly accepted in Micrometerology but has limitations. It draws constant transport coefficients below the zero displacement height d; the stability factor φ H is calculated by sensible heat fluxes above the canopy; the effect of a sparse canopy and beyond the canopy could generate an undesirable damping effect [16,17]. Alternatives to address those limitations could be a non-parabolic approach for the free-surface layer given by a semi-theoretical coefficient of exponential type [17], the extended K-theory [18] or the called L-theory based on a Lagrangian point view of transport [16].
The interchange of sensible heat H l and humidity E l between air and canopy depends on the amount of foliage quantified by the leaf index area (LAI) as a vertical profile of the foliage mass, which will be defined later, The ideal gas law is used to estimate the humidity in the leaf stomata p vs as a function of T l , where in turn M v is the molecular mass of water, R the universal gas constant, and e s (T l ) is the water vapor pressure. Meanwhile, the interchange resistances to convective heat between air and leaf r h , humidity between air and leaf r v , and the evaporation in the leaf stomata r ls , are all defined as being u the wind profile, d l the leaf characteristic length, P the atmospheric pressure, r ls0 the unstressed leaf interchange, Ψ l the leaf water potential, Ψ c the leaf critical water potential, and f St , f T , f VPD are corrections on the stomata resistance to solar radiation, temperature and water vapor deficit, respectively. The wind profile has a logarithm profile above the canopy but exponentially decline in its interior being u * the friction velocity, z m is the momentum transfer roughness, ψ m is the diabatic momentum stability factor, w l is the mean leaf width and LAI tot is the total Leaf Area Index of the canopy. Again, the velocity profile (7) is classical and agree with the comments to Equation (2); there are alternatives, for example, in the free-surface beyond the canopy z > h c , a non-parabolic profile proposed recently in [17]. Other relations needed for computing parameters related to atmospheric stability are: being the atmospheric stability ζ given by ζ = kz re f gH ζ ρ a c a T a u 3 * with the total exchange of sensible heat formulated as where ψ m = 0.6ψ H for unstable conditions and ψ m = ψ H for stable conditions. Finally, L n and S n in (1c) represent the net longwave and shortwave radiations absorbed within the canopy, computed with a multi-layer model [4,9,15,19] (see Appendix A). To understand the effect of the canopy shape in the models' state variables, ideal conditions must be considered. Therefore, the soil layer has a water concentration constant along time, justifying eliminate the PDE related to soil humidity in the model. Moreover, the water vapor on leaf stomata is assumed to be saturated, implying that Equation (5) holds.

The Weibull Distribution and Canopy Shape Characterization
In the above subsection, the vertical interchange of heat and humidity between air and canopy depends on the amount of foliage quantified by LAI and -implicitly-of the foliage mass m c . However, LAI and m c are reported in field measurements as a 2D-flat projection, which in this work are called total leaf area index LAI tot and total canopy biomass m c,tot , respectively. So a vertical parametrization of both LAI and m c is needed. Therefore this subsection defines them as the most straightforward continuous vertical functions using the well-known Weibull distribution, being so one of the main contributions of this work.
The Weibull distribution is a two-parameter distribution used widely in biology, meteorology, sciences materials, and others. It can be described by its probability density function f (v) and its cumulative distribution function F(v) as described in [20]: where β is the Weibull module, α is the scale parameter and v is the independent variable.
Applying these concepts to the tree canopy, another element to consider is the called leaf area density LAD = LAD(z), a vertical shape function that quantified mass change concerning the canopy height. Assuming that LAD(z) is known, we can relate it with different LAI definitions (see [13]) being LAI a the cumulative leaf area index. The shape function LAD(z) can be adjusted by the following Weibulls' probability distribution [14,21,22] Similarly, the cumulative LAI a is defined as the following Weibulls' cumulative distribution where in both cases LAI tot is taken as a proportionality constant. Defining v = 1 − z/h c , we have that LAD(v) = LAI a (v) holds, and LAI is computed straightforwardly in an arbitrary canopy layer of length v 1 − v 0 in the following way taking v 1 = 1 − (z + ∆z)/h c and v 0 = 1 − z/h c as relative positions in the plant canopy, an useful expression of LAI in a canopy layer of length ∆z is formulated as In a totally similar way, the biomass of the foliage m c as a function of z could be formulated as being m tot the total biomass in the canopy.

Initial and Boundary Conditions
In this subsection, the initial boundary conditions will be formulated. For this purpose, it must be noted that the model domain has three physical borders: soil (z = h s ), interface soil-air (z = 0) and air column height (z = h a ). There are no boundaries between the canopy and air; instead, there are interchanges of heat and water, which can be computed using the continuous function LAI (Equation (16)). Therefore, assuming that in an initial time the air temperature, leaf temperature, humidity and soil temperature are known the following initial conditions are imposed to the state variables of the model T a (., . Concerning the boundary conditions, a heat and humidity flux is assumed in the upper part of the air column (z = h a ), and the same in the lower part of the soil layer (z = h s ). Then, with knowing reference time functions T a,re f (t), p v,re f (t) and T s,re f (t), the following boundary conditions are imposed for all t ≥ 0 observing that for large value of the interchanges resistances k a , k v , k s → ∞, the previous conditions are transformed into the Dirichlet type conditions Finally, in the border between the soil layer and the air column (z = 0), interchange of heat and water and soil heat balance are considered. Therefore we impose the following conditions, where H g and E g are sensible heat and vaporization fluxes between air and soil, L v is the heat of vaporization, S n,0 and L n,0 are the net shortwave and longwave radiation absorbed by the soil. These last quantities must consider the process, which begins when the radiation reaches the upper part of the canopy and downward until reaching the soil. Therefore, a multi-layer model is used to approach them (see Appendix A). To compute the air-soil fluxes H g and E g similar formulations those above (in the canopy) are used, then where in turn p vg is the water steam concentration in the soil surface and r s is the soil resistance to convection [16], computed both by the following equations where z re f ,s , d s , z H,s are the reference height, zero displacement height and heat transfer roughness. As before, there are alternatives to the logarithm profile to estimate r s , for example, which is derived from the extended K-Theory and introduced in [23].

Numerical Solution 2.4.1. Finite Element Approach
On the part of Equations (1a), (1b), (1d) the usual finite element treatment will be applied (see [24,25]); however, it must be observed that Equation (1c) has no spatial gradient, so requires an accurate time discretization instead, using finite differences [26]. To avoid tedious calculations and large formulations, a straightforward process is displayed in this subsection. Let the Sobolev spaces V = H 1 (Ω a ) and W = H 1 (Ω s ) and its L 2 (Ω) inner product denoted by ., . . Furthermore, the spacial domain Ω a ∪ Ω s is meshed using two meshes τ a and τ s , given by the partitions h a = z n a > z n a −1 > · · · > z 0 = 0 and 0 =ẑ 0 >ẑ 1 > · · · >ẑ n s = h s with n a + 1 and n s + 1 nodes respectively. It is important consider the existence of a canopy mesh τ c with n c + 1 nodes covering the foliage layer Ω c = (h 0 , h c ) and such that τ c ⊂ τ a . Once the meshes are defined let be the discrete for each i = 0, . . . , n a ; j = 0, . . . , n s and being k a 0 = ρ a c a r s , k v 0 = 1 r s . Using the function basis for V h and W h is it possible to define over τ a and τ s the functions T a , p v , T s , T l as the following lineal combination being ξ j (t) = Φ(z j , t) the time depend state variable and {Φ j (z)} n x j=0 the set of nodal basis functions. Thus, we define the following spacial discretized functions T a (t) = {T a,j (t)} n a j=1 , p v (t) = {p v,j (t)} n s j=1 , T s (t) = {T s,k (t)} n s k=1 , and T l (t) = {T l,j (t)} n a j=n 0 . From the initial conditions, we have known the vectors T a (0 which are the initial conditions for the following system of ordinary differential equations where the contribution of 2 × 2 matrices and 2 × 1 vectors for arbitrary elements I l = [z l , z l+1 ] of τ a andÎ r = [ẑ r ,ẑ r+1 ] of τ s are given by for i, j = l, l + 1, and k, m = r, r + 1. The vectors h a and e v depend on nonlinear relations of temperature and humidity. Taking advantage of the 1D mesh and by simplicity, they are computed using a finite difference approach. A finite element approach is also possible, and it is formulated in Appendix B.

Time Discretization and Relaxed Fixed Point Scheme
As was the case in spatial discretization, the treatment of the system equations will not be the same. Therefore, Equations (24a), (24b), (24d) were time discretized using an Euler Backward scheme. Meanwhile, in Equation (24c), a bi-level type discretization with an Euler-Backward scheme was applied for the first level and a three-point derivative approach on the second level and higher, as is described in [26]. Therefore our problem is now: Given the initial vectors T 0 a , p 0 v , T 0 s , T 0 l find the vectors T n+1 a , p n+1 v , T n+1 s , T n+1 l such that for n = 0, . . . , N − 1 satisfies the following system The system above is nonlinear, so numerical methods such as Piccard iteration (Fixpoint), Newton-Galerkin, or others are needed to approximate the system solution. A first choice is the Piccard method due to its simplicity (it is derivative-free) and when it converges, it does to the solution. However, the complexity and different scales of the empirical formulas, nonlinear coefficients, and correlations could make that the fix-point fails. Anticipating this problem, a more robust variant of the Piccard method consisting of a bi-level iteration of type predictor-corrector was applied q s ] the state variables vector, F the application derived at isolating the linear state variable in each equation of the system, β a relaxation parameter such that 0 < β ≤ 1,û q+1 the relaxed fix-point, u q+1 the predicted fix-point and u q the initial guess. The idea of this bi-level scheme is made β → 0 as at least one state variable has slow convergence in the fixed-point iterations, imposing with it the same convergence velocity for all the state variables in the sense of the following error criteria [25] ( for a predefined tolerance .

Mesh Configuration, Adjustments, and Sensible Analysis
Several numerical experiments were conducted to evaluate the stability, convergence, and parameter sensitivity of the solution. Those experiences showed the need for adjustments in the sub-domains meshes and some parameters. Therefore this subsection is devoted to showing the mesh configuration and solution stabilization criteria.
The abrupt changes between air-foliage-soil suggest large heat and humidity gradients, so to forestall instability, a standard strategy is using no regular meshes. Therefore we maximize node concentration on the border of the different sub-domains. To this end, the following criteria are used to concentrate the nodes of τ a in the interface air-soil and the upper part of the canopy, h c < z ≤ h a n c + n 0 ≤ i < n a + n c + n 0 (29) and for τ s these criteria concentrate the nodes in the interface between the soil and air, The stomata resistance r ls is stressed by solar radiation, temperature, and humidity. We made an adjustment substituted (6c) by a discrete function which is based on a minimum value of radiation and heat to stabilize the interchange between leaf and air Finally, for the initial condition for T s is imposed a piecewise function which is lineal close to the interface soil-air and constant in-depth until z = h s , being h bc < h s the maximum depth of linearization.

Numerical Experiments, Comparative between Different Artificial Canopies
We first define the values of the fixed parameters that characterize all the canopies considered in the numerical experiments and others parameters mentioned through this paper (see Table A1 of Appendix C). However, we emphasizing here the follow parameter values, air column height h a = 50 m, canopy height h c = 3 m, total leaf area index LAI tot = 3.25, total foliage mass m c,tot = 10 kg m −2 , leaves orientation coefficient x f = 1, length d l = 0.08 m, width w l = 0.03 m and stomata resistance r ls0 = 5 s m 3 kg −1 . Furthermore, parameters related to turbulent transport and theoretical wind profile for the air column z re f = 50 m, z H = 0.078 m, d = 0.77 h c and soil resistance, roughness, and soil level displacement z re f ,s = 1 m, z H,s = 0.078 m, d s = 0.001 m which corresponds to grass [27]. In the boundary we impose Dirichlet type conditions with constant values of T a,re f , p v,re f and T s,re f , similarly for initial conditions with constants values for T 0 a , p 0 v and T 0 s . Conversely, we define four different shape canopies using the Weibull distribution, holding all of them the same LAI tot and m c,tot . To do it, we use different combinations of the pair (α c , β c ) (see Table 1 and Figure 1a). The intention is to define four distinguishable canopy shapes with maximum foliage concentration at the top, middle and low parts of the canopy. Additionally, a height-wise foliage concentration in the canopy is added as the last case (see Figure 1a).

Case 1 (Top) Case 2 (Low) Case 3 (Middle) Case 4 (Height-Wise)
α c 0.25 0.75 0.5 0.5 β c 3.5 11.0 7.5 3.25 Between the inputs of the model, we have the vertical wind profile computed by Equation (7) and displayed in Figure 1b, including the line h c = 3 m, and the radiation forcing on the top of the canopy and shown in Figure 1c. Concerning the vertical profile of the wind, hold the decreasing values within the canopy and the exponential beyond it, noting that the interchange resistance coefficients (6a) and (6b) depends on wind speed. For radiation, diffusive and direct shortwave radiations reach their maximum at midday and are zero at night time. Meanwhile, the longwave radiation is set to be constant throughout the day, as it is dependent of the sky temperature.
Before showing the spatial-time evolution of the state variables, Figure 2 shows an approach of the thermal energy change in the canopy computing as ∆E = m c c c (T l − T a,re f ). As expected, the maximum energy change is located where is the maximum foliage mass position from the soil and where the wind velocity is minimum, being Case 2, which both conditions hold (Figure 2b). In contrast, the minimum energy change is a consequence of a foliage mass upper located from the soil and faster wind. In Case 1, both conditions hold with almost zero change (Figure 2a). Additionally, in all canopies (Figure 2), the foliage mass quantity on its limits is so slight that the thermal energy change is practically zero. This leads us to expect maximum leaf temperatures and quick humidity saturation at canopies' ends. However, these extremes values of temperature and humidity will be mitigated within the canopies. This is confirmed by the results explained and shown below. We show in the following figures the space-time evolution of the state variables of the SHAW-derived model for all canopies. Numerical solutions have been computed in the mesh defined in Section 3 and a time interval of 24 h. As initial conditions to leaf-air-soil temperatures and leaf humidity were taken the outputs of a previous models' run of four days.
For the leaf temperature, Cases 1, 2, and 3 are all characterized by maximum temperatures on canopy ends due to minimum foliage mass and forcing by radiation from sky and soil (see Figure 3a-c). In particular, Case 4 shows notably more forcing on the low end of the canopy, in consequence, showing notably more significant temperature change there (due to a small foliage mass and a very short distance to the soil that radiates longwave energy and reflects part of the short wave radiation due to its albedo). Inside the canopy, the foliage mass allows moderate temperatures due to a significant net heat capacity. Now, the temperature of the air column is shown in Figure 4, taking only the part of the air column that matches the canopy height of 3 m. It is important to note that all canopies have the same temperature profile: maximum values close to the soil and decreasing through the canopy. We highlighted Case 2 (Figure 4a), which has a lower located foliage mass and higher energy change ∆E with air. In consequence, this canopy presents the maximum temperatures. Conforming the foliage mass is located in upper positions, the temperatures shown lower values due to a significant distance from the soil and less energy change capacity. Concerning Case 4, a well-mixed air temperature column (Figure 4d) is shown due to its height-wise foliage distribution. Nevertheless, all canopies show a scant variation (≈0.4 K). This could be partly due to a high reference wind velocity u re f = 10 m/s. To demonstrate that, we take lower velocity values u re f = 5 m/s and u re f = 2.5 m/s for Case 2, and results are shown in Figure 5a,b where a major variation (≈1 K) is observed. There are others factors than potentially influence too, such as the constant temperature T a,re f = 293.15 K imposed as boundary condition in the top of the air column and some parameters as the total biomass m tot . About humidity, it is greatly influenced by the saturated soil (because of the more significant moisture flux in the soil-air inter-phase) is apparent the effect of the canopy shape, trapping, and adding (due to leaf stoma) humidity, avoiding its distribution along with the canopy. This can be observed clearly in opposites Cases 1 and 2 (Figure 6a Finally, the soil temperature is computed from surface to a depth of −2 m but only presents a notable variation in a thick layer. Due to this, is shown only an average soil temperature. Then in Figure 7a is shown the average leaf temperature focuses our attention in the midday, wherein Case 2 (low foliage mass location) reaches the maximum temperature, opposite to its minimum corresponding to Case 1 (top foliage mass location). This indicates that to an upper location of the foliage, less soil temperature (due to a higher sensible heat soil cooling caused by steeper air and soil temperature difference). Of course, the model generates other interesting outputs: the net shortwave and longwave radiations in both soil and canopy. Energy fluxes at the soil surface (Equation (19c)) are shown in Figure 7b for all canopies. We expect that the canopy absorbs energy during the day and releases it at night. The shape canopy made a clear difference in magnitude of released-absorbed energy with all canopies having the same exposed behavior. We point out that at midday, a difference of 20 Wm −2 between Cases 1 and 2, being the maximum and minimum values of the four cases.

Discussions
In this work, a SHAW-derived model was formulated and solved numerically with the finite element method, being an alternative to finite difference schemes used by Flerchinger et al. [7][8][9][10][11] and other authors (see reference list [28]). This new approach allows using irregular meshes and potentially scaling the model to a higher dimension which was possible thanks to the continuous vertical function LAI proposed in Section 2.2. Moreover, numerical solutions were derived by combining relaxed fix-point iterations (to address nonlinearity) and a time scheme of predicted-corrector type (to yield stability).
Typically, theoretical and experimental studies of heat and humidity fluxes in canopies are made on a fixed foliage profile either take it as a simple geometry shape -commonly triangular or rectangular-or derivate it by curve fitting of in situ measures. In this work, were compared different canopies characterized exclusively by their foliage profile or shape and analyzed their heat and water transport dynamics. Numerical results suggested that the canopy shape notably influences the vertical profile of leaf, air and soil temperatures, humidity and heat fluxes. Therefore, this work could be a first step in selecting a canopy tree to maximize its benefits. Examples could be several: windbreak, soil cooling or heating, reduced air temperature, maximized air humidity, covering the soil from radiation, and others. Furthermore, the simplicity of the Weibull distribution to characterize a canopy shape with only two parameters opens the possibility of formulating optimization problems related to minimizing harmful phenomenons.
We emphasize that the model could be sensitive to other canopy parameters: leaf width, orientation and diameter, stomata resistance, and specific heat capacity. Moreover, latitude, day of the year, and meteorological parameters are significant too. Therefore, changing these parameters would imply changes in the model outputs. However, address a model sensibility analysis would overextend this paper leaving it to future work.
A value-added of this work is a complete description of the mathematical methodology to obtain the numerical solution of the model, this could be attractive for researchers to reproduce results, update the model, or as another example of how to approach a nonlinear system of PDEs.
Finally, a comparison between the SHAW model and the SHAW-derived model (Equations (1a)-(1d)) could be interesting. However, there are essential differences between them: turbulent transport (K-theory), soil humidity, shape canopy characterization and numerical approach, which have been pointed out through the paper. Indeed, in [16] the substitution of the K-theory by the L-theory made a better approach of the modeled variables to the data. This motivates us to a model update for reach goals such as validation or model comparison. However, we do not consider this observations as a step backward but a direction for future work.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. A Multi-Layer Model to Estimate the Radiation Extinction through the Canopy
Solar radiation is divided into long and shortwave radiation. Moreover, shortwave radiation is subdivided into diffusive or direct; longwave radiation is always diffusive. Moreover, the radiation could be downward or upward, reflected or absorbed. These radiation patterns are essential in the model, so this appendix shows how the amount of solar radiation in the canopy was estimated.
Let be the discrete functions {S d,j } n c +1 j=1 , {S b,j } n c +1 j=1 the diffusive and direct radiation on each layer of the canopy, so in its upper part (j = n c + 1) S d,n c +1 = τ d,sky I sky (A1) where τ d,sky is the atmospheric coefficient of diffuse transmission, τ t,sky is the maximum transmission at the clear sky and I sky is the radiation reaching the upper section of the canopy. The detailed process to compute these parameters could be viewed in [29].
Once the radiation penetrates the canopy, an estimation of the absorbed radiation is needed. This is done by a multi-layer model [4,10,30], which divided the canopy into a finite number of layers, considering that one part of the radiation could pass to the following layer. Another part is reflected turn in diffusive downward and diffusive upward radiations. Therefore, let us consider the case of the short wave radiation, then for an arbitrary layer i, the direct radiation that passes between leaves is meanwhile the diffuse downward radiation S d,i is given by where f b,i,↓↑ and f d,i,↓↑ are, respectively, the upward fraction of direct and diffusive radiation, f b,i,↓↓ and f d,i,↓↓ are, respectively, the downward fraction of direct and diffuse radiation, τ b,i and τ d,i are the percentage of direct and diffuse radiation which pass trough the gap between leaf; meanwhile, τ l is the percentage of radiation which pass trough the leaves, and α l is the leaf albedo (see [4,9,15,19]). Finally, the upwards reflected radiation by the soil is computed as the reflected sum of the direct and diffuse radiations that reach the soil, S u,0 = α s (S d,1 + S b,1 ) and with this last equation, is completed the multi-layer model (A3)-(A6) for the shortwave radiation.
Once the system (A3)-(A6) is solved, the net short wave radiation absorbed by the canopy and soil is computed by the following conservative equations S n,i = (S b,i+1 + S d,i+1 + S u,i−1 ) − (S b,i + S d,i + S u,i ) (A7a) S n,0 = S b,1 + S d,1 − S u,0 With completely similar reasoning, the multi-layer model for the diffuse longwave radiation could be written. Therefore, the longwave radiation that reaches the upper part of the canopy (j = n c + 1) is computed with the Stephan-Boltzman law L d,n c +1 = a σT 4 sky (A8) being a the sky emissivity, σ the Stephan-Boltzman constant, T sky the sky temperature in Kelvin and L d,n c +1 the longwave radiation reaching the canopy.
On an arbitrary layer i of the canopy, the longwave radiation is divided in diffuse downward being T s,0 the ground surface temperature and s the soil emissivity. This last equation completes the multi-layer model by the longwave radiation in the canopy (A9)-(A11).
Once the system (A9)-(A11) is solved, the net radiation absorbed by the canopy and the soil is given by the following conservative equations L n,i = L d,i+1 + L u,i−1 − (L d,i + L u,i ) (A12a) L n,0 = L d,1 − L u,0 (A12b)