Next Article in Journal
Practices of the Random Variable Proposed in the Chilean Mathematics Curriculum of Secondary Education
Previous Article in Journal
Convex Obstacles from Travelling Times
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Carlos E. Villarreal-Olavarrieta
1,
Néstor García-Chan
2,* and
Miguel E. Vázquez-Méndez
3
1
Departmento de Oceanografía Física, Centro de Investigación Científica y de Educación Superior de Ensenada, Ensenada 22860, Mexico
2
Departamento de Física, CU de Cs Exactas e Ingeniería, Universidad de Guadalajara, Guadalajara 44430, Mexico
3
Departamento de Matemática Aplicada, Universidade de Santiago de Compostela, EPSE, 27002 Lugo, Spain
*
Author to whom correspondence should be addressed.
Mathematics 2021, 9(19), 2431; https://doi.org/10.3390/math9192431
Submission received: 20 May 2021 / Revised: 13 September 2021 / Accepted: 23 September 2021 / Published: 30 September 2021
(This article belongs to the Section Computational and Applied Mathematics)

Abstract

:
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 finite difference method and the canopy shape is defined 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 fitting any tree species. This allows formulating a fully continuous SHAW-derived model, which is numerically solved by a finite element approach of P1 Lagrange type. For this novel approach, several numerical experiments were carried out to understand how the shape of well distinguishable canopies influences heat and water transport.

1. 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, multi-species, and soil coupled model uses a continuous vertical function of LAD to upscale CO2, 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.

2. Materials and Methods

2.1. 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
ρ a c a t T a = z [ ρ a c a K e z T a + H l ( T a , T l ) ] on Ω a × [ 0 , T ]
t p v = z [ K e z p v + E l ( p v , T l ) ] on Ω a × [ 0 , T ]
m c c c z T l = S n + L n H l ( T a , T l ) L v E l ( p v , T l ) on Ω c × [ 0 , T ]
ρ s c s t T s = z ( λ s z T s ) on Ω s × [ 0 , T ]
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
K e = k u ( z d + z H ) / ϕ H z > d k u z H / ϕ H z d
being k the Von Karman constant, z H the canopy height, d 0.77 h c the zero-level displacement, ϕ H is the dimensionless heat stability factor, and u the friction velocity given by
u = u r e f k ln z r e f d z m + ψ m 1
where u r e f is the reference wind speed at a given height z r e 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 ( L A I ) as a vertical profile of the foliage mass, which will be defined later,
H l = 2 ρ a c a L A I T l T a r h
E l = L A I p v s p v r l s + r v p v s p v 0 p v s < p v
The ideal gas law is used to estimate the humidity in the leaf stomata p v s as a function of T l ,
p v s = M v R T l e s ( 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 l s , are all defined as
r h = 7.4 ρ a M a d l u = 7.4 P R T a d l u
r v = 6.8 ρ a M a d l u = 6.8 P R T a d l u
r l s = r l s 0 f S t f T f V P D 1 + Ψ l Ψ c n s
being u the wind profile, d l the leaf characteristic length, P the atmospheric pressure, r l s 0 the unstressed leaf interchange, Ψ l the leaf water potential, Ψ c the leaf critical water potential, and f S t , f T , f V P D 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
u = u k ln z d z m + ψ m z h c u k ln h c d z m + ψ m exp 0.28 L A I t o t 2 h c w l 1 / 3 z h c 1 z < h c
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 L A I t o t 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:
ψ H = 4.7 ζ ζ 0 2 ln 1 + 1 16 ζ 2 ζ < 0
ϕ H = ( 1 16 ζ ) 1 / 2 ζ 0 1 + 6 ζ 1 + ζ ζ < 0
being the atmospheric stability ζ given by
ζ = k z r e f g H ζ ρ a c a T a u 3
with the total exchange of sensible heat formulated as
H ζ = k 2 ρ a c a u ( z ) [ T a ( d + z H ) T a ( z ) ] ln z d z m + ψ m ln z d z H + ψ H
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.

2.2. 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 L A I and —implicitly— of the foliage mass m c . However, L A I and m c are reported in field measurements as a 2D-flat projection, which in this work are called total leaf area index L A I t o t and total canopy biomass m c , t o t , respectively. So a vertical parametrization of both L A I 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]:
f ( v ) = β α v α β 1 exp v c β
F ( v ) = 1 exp v c β
where β is the Weibull module, α is the scale parameter and v is the independent variable. Is very important to note that f ( v ) is the anti-derivative of F ( v ) ( f ( v ) = F ( v ) ). Applying these concepts to the tree canopy, another element to consider is the called leaf area density L A D = L A D ( z ) , a vertical shape function that quantified mass change concerning the canopy height. Assuming that L A D ( z ) is known, we can relate it with different L A I definitions (see [13])
L A I t o t = 0 h c L A D ( z ) d z
L A I a = z h c L A D ( s ) d s
being L A I a the cumulative leaf area index. The shape function L A D ( z ) can be adjusted by the following Weibulls’ probability distribution [14,21,22]
L A D = L A I t o t β c α c 1 z / h c α c β c 1 exp 1 z / h c α c β c
Similarly, the cumulative L A I a is defined as the following Weibulls’ cumulative distribution
L A I a = L A I t o t 1 exp 1 z / h c α c β c
where in both cases L A I t o t is taken as a proportionality constant. Defining v = 1 z / h c , we have that L A D ( v ) = L A I a ( v ) holds, and L A I is computed straightforwardly in an arbitrary canopy layer of length v 1 v 0 in the following way
L A I = v 0 v 1 L A D ( v ) d v = L A I a ( v 1 ) L A I a ( v 0 )
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 L A I in a canopy layer of length Δ z is formulated as
L A I ( z ) = L A I t o t exp 1 z / h c α c β c exp 1 ( z + Δ z ) / h c α c β c
In a totally similar way, the biomass of the foliage m c as a function of z could be formulated as
m c ( z ) = m t o t exp 1 z / h c α c β c exp 1 ( z + Δ z ) / h c α c β c
being m t o t the total biomass in the canopy.

2.3. 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 L A I (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 ( . , 0 ) = T a 0 ( . ) , p v ( . , 0 ) = p v 0 ( . ) , T l ( . , 0 ) = T l 0 ( . ) and T s ( . , 0 ) = T s 0 ( . ) . 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 , r e f ( t ) , p v , r e f ( t ) and T s , r e f ( t ) , the following boundary conditions are imposed for all t 0
ρ a c a K e z T a ( h a , t ) = k a ( T a , r e f ( t ) T a ( h a , t ) )
K e z p v ( h a , t ) = k v ( p v , r e f ( t ) p v ( h a , t ) )
λ z T s ( h s , t ) = k s ( T s , r e f ( t ) T s ( h s , t ) )
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 T a ( h a , t ) = T a , r e f ( t ) , p v ( h a , t ) = p v , r e f ( t ) and T s ( h s , t ) = T s , r e f ( t ) .
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,
ρ a c a K e z T a ( 0 , t ) = H g ( T a , T s )
K e z p v ( 0 , t ) = E g ( p v , T s )
λ s z T s ( 0 , t ) = S n , 0 + L n , 0 H g ( T a , T s ) L v E g ( p v , T s )
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
H g = ρ a c a T s T a r s z = 0
E g = p v g p v r s z = 0 p v g p v 0 p v g < p v
where in turn p v g 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
p v g = M v R T s e s ( T s ) z = 0
r s = 1 u k ln z r e f , s d s + z H , s z H , s + ψ H
where z r e 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].

2.4. 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 = z ^ 0 > z ^ 1 > > z ^ 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 spaces V h = { φ h C ( Ω ¯ a ) : φ h | I i P 1 ( I i ) } and W h = { φ ^ h C ( Ω ¯ s ) : φ ^ h | I ^ i P 1 ( I ^ i ) } such that V V h and W W h , then it is possible formulate the following semi-discrete problem: Find ( T a , p v , T l , T s ) V h × V h × V h × W h and T l C ( Ω c ) such that for φ i V h and φ ^ m W h satisfies the following system equations
ρ a c a t T a , φ i + ρ a c a K e z T a , φ i + k a T a ( h a , t ) φ i ( h a ) + k a 0 T a ( 0 , t ) φ i ( 0 ) = z H l , φ i + k a T a , r e f φ i ( h a ) + k a 0 T s ( 0 , t ) φ i ( 0 )
t p v , φ i + K e z p v , φ i + k v p v ( h a , t ) φ i ( h a ) + k v 0 p v ( 0 , t ) φ i ( 0 ) = z E l , φ i + k v p v , r e f φ i ( h a ) + k v 0 p v g ( 0 , t ) φ i ( 0 )
m c C c t T l = S n + L n H l L v E l
ρ s C s t T s , φ ^ j + λ s z T s , z φ ^ j k s T s ( h s , t ) φ ^ j ( h s ) = ( S n , 0 + L n , 0 H g L v E g ) φ ^ j ( 0 ) k s T s , r e f φ ^ j ( h s )
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
Φ ( z , t ) = j = 0 n x ξ j ( t ) Φ ^ j ( z )
being ξ j ( t ) = Φ ( z j , t ) the time depend state variable and { Φ ^ j ( z ) } j = 0 n x the set of nodal basis functions. Thus, we define the following spacial discretized functions T a ( t ) = { T ˜ a , j ( t ) } j = 1 n a , p v ( t ) = { p ˜ v , j ( t ) } j = 1 n s , T s ( t ) = { T ˜ s , k ( t ) } k = 1 n s , and T l ( t ) = { T ˜ l , j ( t ) } j = n 0 n a . From the initial conditions, we have known the vectors T a ( 0 ) = { T a , i 0 } i = 0 n a , p v ( 0 ) = { p v , i } i = 0 n a , T s ( 0 ) = { T s , m 0 } m = n s 0 , T l ( 0 ) = { T l , i 0 } i = n 0 n c which are the initial conditions for the following system of ordinary differential equations
ρ a c a M d T a d t ( t ) + [ R a + B a ] T a ( t ) = h a ( t ) + b a ( t )
M d p v d t ( t ) + [ R v + B v ] p v ( t ) = e v ( t ) + b v ( t )
C c m c d T l d t ( t ) = S n ( t ) + L n ( t ) H l ( t ) L v E l ( t )
M s d T s d t ( t ) + R s T s ( t ) = b s ( t )
where the contribution of 2 × 2 matrices and 2 × 1 vectors for arbitrary elements I l = [ z l , z l + 1 ] of τ a and I ^ r = [ z ^ r , z ^ r + 1 ] of τ s are given by
( M ) i , l l = φ m , φ l
( R a ) i , j l = ρ a c a K e , j φ j , φ i
( B a ) i , j l = k L a φ j ( h a ) φ i ( h a ) + k a 0 φ j ( 0 ) φ i ( 0 )
( h a ) i l = z H l , φ i
( b a ) i l = k a 0 T s ( 0 , t ) φ i ( 0 ) + k a T a , r e f φ i ( h a )
( R v ) i , j l = K e , j φ j , φ i
( B v ) i , j l = k v φ j ( h a ) φ i ( h a ) + k v 0 φ j ( 0 ) φ i ( 0 )
( b v ) i l = k v 0 p v g ( 0 , t ) φ i ( 0 ) + k v p v , r e f φ i ( h a )
( e v ) i l = z E l , φ i
( m c ) i l = m c ( z i ) , ( S n ) i l ( t ) = S n , i , ( L n ) i l ( t ) = L n , i
( M s ) m , k r = φ ^ k , φ ^ m
( R s ) m , k r = λ s φ ^ k , φ ^ m + k s φ ^ m ( h s ) φ k ( h s )
( b s ) k r ( t ) = ( S g + L g H g L v E g ) φ ^ k ( 0 ) + k s T s , r e f φ ^ k ( h s )
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.

2.4.2. 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 a 0 , p v 0 , T s 0 , T l 0 find the vectors T a n + 1 , p v n + 1 , T s n + 1 , T l n + 1 such that for n = 0 , , N 1 satisfies the following system
ρ a c a M T a n + 1 T a n Δ t + [ R a + B a ] T a n + 1 = h a n + 1 + b a n + 1
M p v n + 1 p v n Δ t + [ R v + B v ] p v n + 1 = e v n + 1 + b v n + 1
C c m c T l n + 1 T l n Δ t = G n + 1 , n = 0 C c m c 3 T l n + 1 4 T l n + T l n 1 2 Δ t = G n + 1 , n = 1 , , N 1
M s T s n + 1 T s n Δ t + R s T s n + 1 = b s n + 1
being G n + 1 = S n n + 1 + L n n + 1 H l n + 1 L v E l n + 1 .
The system above is nonlinear, so numerical methods such as Piccard iteration (Fix-point), 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
u q + 1 = F ( u q )
u ^ q + 1 = β u q + 1 + ( β 1 ) u q
being q the fix-point index, u q = [ T a q , p v q , T l q , T s q ] 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 , u ^ 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]
( u q + 1 u q ) 2 ( u q + 1 ) 2 < ϵ
for a predefined tolerance ϵ .

2.5. 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,
z i = h 0 i n 0 1 2 0 z < h 0 , 0 < i < n 0 z i 1 + h c h 0 n c h 0 z h c , n 0 i < n c + n 0 h c + h a h c i N C n 0 + 1 n a 2 h c < z h a n c + n 0 i < n a + n c + n 0
and for τ s these criteria concentrate the nodes in the interface between the soil and air,
z ı ^ = h s ı ^ n s 1 2 0 z < h s , 0 ı ^ < n s
The stomata resistance r l s 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
r l s , i = r l s 0 , i f ( S n , i S n , m i n ) ( T l , i > T l , m i n ) ( p v s > p v ) r l s 0 S n , m i n S n , i , i f ( S n , i < S n , m i n ) ( T l , i > T l , m i n ) ( p v s > p v ) , i f ( S n , i = 0 ) ( T l , i < T l , m i n ) ( p v s p v )
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 ,
T s ( z ı ^ , 0 ) = T s , r e f T a , r e f z b c z ı ^ + T a , r e f z ı ^ h b c T s , r e f z ı ^ > h b c
being h b c < h s the maximum depth of linearization.

3. Results

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 L A I t o t = 3.25 , total foliage mass m c , t o t = 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 l s 0 = 5 s m3 kg−1. Furthermore, parameters related to turbulent transport and theoretical wind profile for the air column z r e f = 50 m, z H = 0.078 m, d = 0.77 h c and soil resistance, roughness, and soil level displacement z r e 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 , r e f , p v , r e f and T s , r e f , similarly for initial conditions with constants values for T a 0 , p v 0 and T s 0 . Conversely, we define four different shape canopies using the Weibull distribution, holding all of them the same L A I t o t and m c , t o t . 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).
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 , r e 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 r e f = 10 m/s. To demonstrate that, we take lower velocity values u r e f = 5 m/s and u r e 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 , r e f = 293.15 K imposed as boundary condition in the top of the air column and some parameters as the total biomass m t o t .
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,b) wherein the humidity line reaches 2.5 m in the first case, meanwhile only reaches 1.5 m in Case 1, matching with their respective canopies top.
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.

4. 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.

Author Contributions

Conceptualization, C.E.V.-O., N.G.-C. and M.E.V.-M.; methodology C.E.V.-O. and N.G.-C.; investigation C.E.V.-O. and N.G.-C. writing—original draft preparation C.E.V.-O., N.G.-C. and M.E.V.-M. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

First author thanks the support from CONACyT (México) under Grant 766062. Second author thanks the support from Sistema Nacional de Investigadores (México) under Grant SNI-52768, Programa para el Desarrollo Profesional Docente (México) under Grant PRODEP/103.5/ 16/8066 and CONACyT by Ciencia de Frontera (México) under Grant 217556. Also, the authors appreciate the comments of the anonymous reviewers that helped to improve this paper.

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 } j = 1 n c + 1 , { S b , j } j = 1 n c + 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 , s k y I s k y
S b , n c + 1 = τ t , s k y τ d , s k y I s k y
where τ d , s k y is the atmospheric coefficient of diffuse transmission, τ t , s k y is the maximum transmission at the clear sky and I s k y 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
S b , i = τ b , i S b , i + 1
meanwhile the diffuse downward radiation S d , i is given by
S d , i = τ d , i + α l f d , i , + τ l f d , i , 1 τ d , i S d , i + 1 + α l f b , i , + τ l f b , i , 1 τ b , i S b , i + 1 + α l f d , i , + τ l f d , i , 1 τ d , i S u , i 1
and the diffuse upward radiation S u , i is given by
S u , i = τ d , i + α l f d , i , + τ l f d , i , 1 τ d , i S u , i 1 + α l f d , i , + τ l f d , i , 1 τ d , i S d , i + 1 + α l f b , i , + τ l f b , i , 1 τ b , i S b , i + 1
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 )
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 s k y 4
being ϵ a the sky emissivity, σ the Stephan-Boltzman constant, T s k y 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
L d , i = τ d , i + 1 τ d , i 1 ϵ c f d , i , L d , i + 1 + 1 τ d , i 1 ϵ c f d , i , L u , i 1 + 1 τ d , i ϵ c σ T l , i 4
and diffuse upward
L u , i = τ d , i + 1 τ d , i 1 ϵ c f d , i , L u , i 1 + 1 τ d , i 1 ϵ c f d , i , L d , i + 1 + 1 τ d , i ϵ c σ T l , i 4
meanwhile the longwave radiation absorbed by the soil is given by the balance
L u , 0 = 1 ϵ s L d , 1 + ϵ s σ T s , 0 4
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
L n , 0 = L d , 1 L u , 0

Appendix B. Computing Contributions Matrices of the Gradients of El and Hl

We focus on two terms in the semi-discrete problem (22), which was approached by a finite difference scheme, is now listed as a finite element approach. We begin with z E l , φ i , doing formal operations we have
z E l , φ i = z L A I r l s + r v ( p v s p v ) , φ i = z L A I r l s + r v p v s , φ i z L A I r l s + r v p v , φ i
being p v s = p v s ( T l ) a nonlinear relation and where p v admits a representation as lineal combination, then we write
z L A I r l s + r v ( p v s p v ) , φ i = z , j L A I r l s + r v p v s ( T ˜ l , j ) , φ i j = 1 n a η j z L A I r l s + r v φ j , φ i
now, taking an arbitrary element I i = [ z i , z i + 1 ] we compute its contribution to the system matrix
z L A I r l s + r v φ l , φ m m l i = z i z i + 1 z L A I r l s + r v φ l φ m d z = L A I ¯ r l s + r v ¯ z i z i + 1 φ l φ m d z
for m , l = i , i + 1 and where L A I ¯ , r l s + r v ¯ are the average on the element I i . Given the value to sub-indexes m , l the follow contribution matrix is obtained
K v i = L A I ¯ r l s + r v ¯ z i z i + 1 φ i φ i φ i + 1 φ i φ i φ i + 1 φ i + 1 φ i + 1 d z = L A I ¯ r l s + r v ¯ 1 1 1 1
On other part the first term at element I i is computed as
z , j L A I r l s + r v p v s ( T ˜ l , j ) , φ i m i = L A I ¯ r l s + r v ¯ z i z i + 1 z p v s ( T ˜ l ) φ i z p v s ( T ˜ l ) φ i + 1 d z
applying a forward-backward derivative approximation and a simple quadrature rule we have
k v i = L A I ¯ r l s + r v ¯ z i z i + 1 z p v s ( T ˜ l ) φ i z p v s ( T ˜ l ) φ i + 1 d z = L A I ¯ r l s + r v ¯ p v s , i φ i , i + p v s , i + 1 φ i , i + 1 2 Δ z p v s , i φ i + 1 , i + p v s , i + 1 φ i + 1 , i + 1 2 Δ z = L A I ¯ r l s + r v ¯ p v s , i + 1 p v s , i 2 p v s , i + 1 p v s , i 2 = L A I ¯ 2 r l s + r v ¯ ( p v s , i + 1 p v s , i ) 1 1 p v s p v 0 p v s < p v
Respect to the second term z H l , φ i in the semi-discretized model (22), as above, with formal operations we have
z H l , φ i = z 2 ρ a c a L A I r h ( T l , j T a , j ) , φ i = 2 ρ a c a z L A I r h T l , j , φ i z L A I r h T a , j , φ i = 2 ρ a c a j = 1 n a ξ l , j z ( L A I r h φ j ) , φ i ξ a , j z ( L A I r h φ j ) , φ i
evaluating on arbitrary element I i = [ z i , z i + 1 ]
z L A I r h φ l , φ m = z i z i + 1 z L A I r h φ l φ m d z = L A I ¯ r h ¯ z i z i + 1 φ l φ m d z
being r h ¯ is an average resistance in the element. Using a simple quadrature rule and forward-backward derivative approximation
K a i = 2 ρ a c a L A I ¯ r h ¯ z i z i + 1 φ i φ i φ i + 1 φ i φ i φ i + 1 φ i + 1 φ i + 1 d z = ρ a c a L A I ¯ r h ¯ 1 1 1 1
therefore
z 2 ρ a c a L A I r h ( T l , j T a , j ) , φ i = ρ a c a j = 1 n a ξ l , j L A I ¯ r h ¯ 1 1 1 1 ξ a , j L A I ¯ r h ¯ 1 1 1 1

Appendix C. Nomenclature and Values of the SHAW-derived Model Parameters

Table A1. Paremeter values used in the numerical experiences.
Table A1. Paremeter values used in the numerical experiences.
Symbol/ValueDefinitionSymbol/ValueDefinition
t f ( h ) = 24 hSimulation time L v = 2.45 × 10 6 J kg 1 vaporization enthalpy
c a = 1004.5 J ( kg K ) 1 specific heat air ρ s c s = 3.7 × 10 6 J m 3 K 1 soil conductivity
c c = 2200 J ( kg K ) 1 specific heat foliage λ = 53.66 ° latitude
c s = 3.77 × 10 6 J ( kg K ) 1 specific heat soil d = 0.77 h c zero level displacemet
d l = 0.08 m characteristic leaf length d s = 0.01 soil level displacement
d n = 172 day number of the year τ l = 0.35 foliage transitivity
h a = 50 m air column height Ω = 0.9 pile leaves factor
h c = 3 m canopy height τ t = 0.6 air transitivity
h s = 2 m soil layer width m c = 10 kg m 2 canopy biomass
α s = 0.08 soil albedo w l = 0.03 m mean leaf width
α c = 0.513 foliage vertical profile ρ a = 1.21 kg m 3 air density
α l = 0.15 leaf albedo L A I t = 2.25 total LAI
ϵ c = 0.90 canopy emissivity r l s 0 = 5 s m 1 stomata resistance
ϵ s = 0.94 soil emissivity T s k y = 278.15 K sky temperature
u r e f = 10 m s 1 reference velocity T a , r e f = 293.15 K air reference temperature
z r e f = 50 m reference height T a 0 = 293.15 K initial air temperature
z r e f , s = 1 soil reference height T l 0 = 280 K initial leaf temperature
z H , s = 0.001 soil convection roughness T s , r e f = 293.0 K reference soil temperature
n a = 30 mesh elements in air T s 0 = 293.15 K initial soil temperature
n c = 500 mesh elements in canopy λ s = 2.9 J / ( m K ) soil thermal conductivity
n s = 30 mesh elements in soil Δ t = 15 s time step
n 0 = 25 mesh elements in trunk p v , r e f = 0.013 kg / m 3 air reference humidity
S n , m i n = 0.5 W m 2 minimum canopy radiation p v 0 = 0.012 kg / m 3 initial humidity
T l , m i n = 277 K minimum stomata temperature

References

  1. Bonan, G.B.; Patton, E.G.; Finnigan, J.J.; Baldocchi, D.D.; Harman, I.N. Moving beyond the incorrect but useful paradigm: Reevaluating big-leaf and multilayer plant canopies to model biosphere-atmosphere fluxes—A review. Agric. For. Meteorol. 2021, 306, 108435. [Google Scholar] [CrossRef]
  2. Yun, S.H.; Park, C.Y.; Kim, E.S.; Lee, D.K. A Multi-Layer Model for Transpiration of Urban Trees Considering Vertical Structure. Forests 2020, 11, 1164. [Google Scholar] [CrossRef]
  3. Dargahi, M.; Newson, T.; Moore, J.R. A Numerical Approach to Estimate Natural Frequency of Trees with Variable Properties. Forests 2020, 11, 915. [Google Scholar] [CrossRef]
  4. Zhao, W.; Qualls, R.J. A multiple-layer canopy scattering model to simulate shortwave radiation distribution within a homogeneous plant canopy. Water Resour. Res. 2005, 41. [Google Scholar] [CrossRef]
  5. Wu, C.; Chau, K.; Huang, J. Modelling coupled water and heat transport in a soil–mulch–plant–atmosphere continuum (SMPAC) system. Appl. Math. Model. 2007, 31, 152–169. [Google Scholar] [CrossRef] [Green Version]
  6. Banimahd, S.; Zand-Parsa, S. Simulation of evaporation, coupled liquid water, water vapor and heat transport through the soil medium. Agric. Water Manag. 2013, 130, 168–177. [Google Scholar] [CrossRef]
  7. Flerchinger, G.; Pierson, F. Modeling plant canopy effects on variability of soil temperature and water. Agric. For. Meteorol. 1991, 56, 227–246. [Google Scholar] [CrossRef]
  8. Flerchinger, G.N.; Kustas, W.P.; Weltz, M.A. Simulating Surface Energy Fluxes and Radiometric Surface Temperatures for Two Arid Vegetation Communities Using the SHAW Model. J. Appl. Meteorol. 1998, 37, 449–460. [Google Scholar] [CrossRef] [Green Version]
  9. Flerchinger, G.; Yu, Q. Simplified expressions for radiation scattering in canopies with ellipsoidal leaf angle distributions. Agric. For. Meteorol. 2007, 144, 230–235. [Google Scholar] [CrossRef]
  10. Flerchinger, G.; Xiao, W.; Sauer, T.; Yu, Q. Simulation of within-canopy radiation exchange. NJAS-Wagening. J. Life Sci. 2009, 57, 5–15. [Google Scholar] [CrossRef] [Green Version]
  11. Flerchinger, G.N.; Reba, M.L.; Link, T.E.; Marks, D. Modeling temperature and humidity profiles within forest canopies. Agric. For. Meteorol. 2015, 213, 251–262. [Google Scholar] [CrossRef]
  12. Coops, N.C.; Smith, M.L.; Jacobsen, K.L.; Martin, M.; Ollinger, S. Estimation of plant and leaf area index using three techniques in a mature native eucalypt canopy. Austral Ecol. 2004, 29, 332–341. [Google Scholar] [CrossRef]
  13. Mori, S.; Hagihara, A. Crown profile of foliage area characterized with the Weibull distribution in a hinoki (Chamaecyparis obtusa) stand. Trees 1991, 5, 1432–2285. [Google Scholar] [CrossRef]
  14. Launiainen, S.; Katul, G.G.; Lauren, A.; Kolari, P. Coupling boreal forest CO2, H2O and energy flows by a vertically structured forest canopy–Soil model with separate bryophyte layer. Ecol. Model. 2015, 312, 385–405. [Google Scholar] [CrossRef]
  15. Campbell, G.S.; Norman, J. An Introduction to Environmental Biophysics; Springer Science + Business Media: New York, NY, USA, 1988. [Google Scholar] [CrossRef]
  16. Flerchinger, G.N.; Reba, M.L.; Marks, D. Measurement of Surface Energy Fluxes from Two Rangeland Sites and Comparison with a Multilayer Canopy Model. J. Hydrometeorol. 2012, 13, 1038–1051. [Google Scholar] [CrossRef]
  17. Absi, R. Reinvestigating the Parabolic-Shaped Eddy Viscosity Profile for Free Surface Flows. Hydrology 2021, 8, 126. [Google Scholar] [CrossRef]
  18. Harman, I.N.; Finnigan, J.J. Scalar Concentration Profiles in the Canopy and Roughness Sublayer. Bound.-Layer Meteorol. 2008, 129, 323–351. [Google Scholar] [CrossRef]
  19. Liu, R.; Huang, W.; Ren, H.; Yang, G.; Wang, J.; Li, X. Research on FPAR vertical distribution in different variety maize canopy. In Proceedings of the 2011 IEEE International Geoscience and Remote Sensing Symposium, Vancouver, BC, Canada, 24–29 July 2011; pp. 2745–2748. [Google Scholar] [CrossRef]
  20. Lazzari, L. 8-Statistical Analysis of Corrosion Data. In Engineering Tools for Corrosion; Lazzari, L., Ed.; European Federation of Corrosion (EFC) Series; Woodhead Publishing: Kidlington, UK, 2017; pp. 131–148. [Google Scholar] [CrossRef]
  21. Yang, X.; Witcosky, J.J.; Miller, D.R. Vertical Overstory Canopy Architecture of Temperate Deciduous Hardwood Forests in the Eastern United States. For. Sci. 1999, 45, 349–358. [Google Scholar] [CrossRef]
  22. Coops, N.C.; Hilker, T.; Wulder, M.A.; St-Onge, B.; Siggins, A.; Trofymow, J.A. Estimating canopy structure of Douglas-fir forest stands from discrete-return LiDAR. Trees 1999, 21, 349–358. [Google Scholar] [CrossRef] [Green Version]
  23. Oleson, K.W.; Lawrence, D.M.; Bonan, G.B.; Flanner, M.G.; Kluzek, E.; Lawrence, P.J.; Levis, S.; Swenson, S.C.; Thornton, P.E.; Dai, A.; et al. Technical Description of Version 4.0 of the Community Land Model (CLM); Technical Report; National Center for Atmospheric Research: Boulder, CO, USA, 2010. [Google Scholar]
  24. Larson, M.G.; Bengzon, F. The Finite Element Method: Theory, Implementation, and Applications; Springer: Berlin/Heidelberg, Germany, 2013. [Google Scholar] [CrossRef]
  25. Reddy, J.N. An Introduction to Nonlinear Finite Element Analysis with Applications to Heat Transfer, Fluid Mechanics, and Solid Mechanics, 2nd ed.; Oxford University Press: Oxford, UK, 2014. [Google Scholar] [CrossRef]
  26. Nakamura, S. Applied Numerical Methods with Software; Prentice-Hall: Harlow, UK, 1991; ISBN 9780130410474. [Google Scholar]
  27. Hansen, F.V. Surface Roghness Lengths; Technical Report; U.S. Army Research Laboratory: Adelphi, MD, USA, 1994. [Google Scholar]
  28. SHAW Model: USDA ARS. 2021. Available online: https://www.ars.usda.gov/ARSUserFiles/20520500/SHAW/302/SHAWReferences.pdf’ (accessed on 4 February 2021).
  29. Iqbal, M. An Introduction to Solar Radiation; Academic Press Canada: Waterloo, ON, Canada, 1983. [Google Scholar] [CrossRef]
  30. Zhao, W.; Qualls, R.J. Modeling of long-wave and net radiation energy distribution within a homogeneous plant canopy via multiple scattering processes. Water Resour. Res. 2006, 42. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Inputs for the SHAW-derived model. (a) Different canopies shapes derived with the function L A I ( z ) given by (13) and different combinations of the parameters pair ( α c , β c ) . Cases, 1 (solid), 2 (dashed), 3 (dotted), and 4 (dash-dot). (b) Zoom of the Wind-profile for the four canopies shape, the black line is the canopy height. (c) 24 h radiation forcing from the sky at the upper layer in the canopy, Short wave-direct (solid), Short wave-diffusive (dotted), and Long wave-direct (dashed).
Figure 1. Inputs for the SHAW-derived model. (a) Different canopies shapes derived with the function L A I ( z ) given by (13) and different combinations of the parameters pair ( α c , β c ) . Cases, 1 (solid), 2 (dashed), 3 (dotted), and 4 (dash-dot). (b) Zoom of the Wind-profile for the four canopies shape, the black line is the canopy height. (c) 24 h radiation forcing from the sky at the upper layer in the canopy, Short wave-direct (solid), Short wave-diffusive (dotted), and Long wave-direct (dashed).
Mathematics 09 02431 g001
Figure 2. Energy change Δ E at 24 h for all canopy shapes. (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4.
Figure 2. Energy change Δ E at 24 h for all canopy shapes. (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4.
Mathematics 09 02431 g002
Figure 3. Leaf temperature variations at 24 h for all canopy shapes. (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4.
Figure 3. Leaf temperature variations at 24 h for all canopy shapes. (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4.
Mathematics 09 02431 g003
Figure 4. Air column temperatures at 24 h for all canopy shapes. (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4.
Figure 4. Air column temperatures at 24 h for all canopy shapes. (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4.
Mathematics 09 02431 g004
Figure 5. Wind velocity partly explains the scant variability of the temperature in the air column. Temperatures at 24 h for the canopy shape 2 subjected to different wind velocities. (a) Case 2 u r e f = 5 m/s; (b) Case 2 u r e f = 2.5 m/s.
Figure 5. Wind velocity partly explains the scant variability of the temperature in the air column. Temperatures at 24 h for the canopy shape 2 subjected to different wind velocities. (a) Case 2 u r e f = 5 m/s; (b) Case 2 u r e f = 2.5 m/s.
Mathematics 09 02431 g005
Figure 6. Leaf humidity at 24 h for all canopy shapes. (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4.
Figure 6. Leaf humidity at 24 h for all canopy shapes. (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4.
Mathematics 09 02431 g006
Figure 7. Temperature and heat balance in the soil, Cases 1 (solid), 2 (dashed), 3 (dotted), and 4 (dash-dot). In essence, the soil temperature and heat balance is the same for all the canopies, maximum values at midday and energy released at night. However, Cases 1 and 2 represent the extreme values from the rest of them. (a) Soil temperature (averaged in width); (b) Energy balance.
Figure 7. Temperature and heat balance in the soil, Cases 1 (solid), 2 (dashed), 3 (dotted), and 4 (dash-dot). In essence, the soil temperature and heat balance is the same for all the canopies, maximum values at midday and energy released at night. However, Cases 1 and 2 represent the extreme values from the rest of them. (a) Soil temperature (averaged in width); (b) Energy balance.
Mathematics 09 02431 g007
Table 1. Different pair values ( α c , β c ) of the parameters given to the Weibull distribution to define four distinguishable canopy shapes.
Table 1. Different pair values ( α c , β c ) of the parameters given to the Weibull distribution to define four distinguishable canopy shapes.
Case 1 (Top)Case 2 (Low)Case 3 (Middle)Case 4 (Height-Wise)
α c 0.250.750.50.5
β c 3.511.07.53.25
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Villarreal-Olavarrieta, C.E.; García-Chan, N.; Vázquez-Méndez, M.E. Simulation of Heat and Water Transport on Different Tree Canopies: A Finite Element Approach. Mathematics 2021, 9, 2431. https://doi.org/10.3390/math9192431

AMA Style

Villarreal-Olavarrieta CE, García-Chan N, Vázquez-Méndez ME. Simulation of Heat and Water Transport on Different Tree Canopies: A Finite Element Approach. Mathematics. 2021; 9(19):2431. https://doi.org/10.3390/math9192431

Chicago/Turabian Style

Villarreal-Olavarrieta, Carlos E., Néstor García-Chan, and Miguel E. Vázquez-Méndez. 2021. "Simulation of Heat and Water Transport on Different Tree Canopies: A Finite Element Approach" Mathematics 9, no. 19: 2431. https://doi.org/10.3390/math9192431

APA Style

Villarreal-Olavarrieta, C. E., García-Chan, N., & Vázquez-Méndez, M. E. (2021). Simulation of Heat and Water Transport on Different Tree Canopies: A Finite Element Approach. Mathematics, 9(19), 2431. https://doi.org/10.3390/math9192431

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