On the Existence of Solutions of a Two-Layer Green Roof Mathematical Model

: The aim of this article is to ﬁll part of the existing gap between the mathematical modeling of a green roof and its computational treatment, focusing on the mathematical analysis. We ﬁrst introduce a two-dimensional mathematical model of the thermal behavior of an extensive green roof based on previous models and secondly we analyze such a system of partial differential equations. The model is based on an energy balance for buildings with vegetation cover and it is presented for general shapes of roofs. The model considers a vegetable layer and the substratum and the energy exchange between them. The unknowns of the problem are the temperature of each layer described by a coupled system of two partial differential equations of parabolic type. The equation modeling the evolution of the temperature of the substratum also considers the change of phase of water described by a maximal monotone graph. The main result of the article is the proof of the existence of solutions of the system which is given in detail by using a regularization of the maximal monotone graph. Appropriate estimates are obtained to pass to the limit in a weak formulation of the problem. The result goes one step further from modeling to validate future numerical results.


Introduction
A green roof is a roof that contains a soil (growing media) and vegetation layer as its outermost surface. Green roofs are used in buildings since ancient times, they have multiple benefits, among them, the reduction of the effect of the urban heat island, the building energy savings, the storm water reduction, aesthetic effect and the acoustic benefits (see [1][2][3] and references therein for more details).
Two main types of green roof exist: extensive green roof and intensive green roof. In the extensive green roofs, the depth of the growing media is less than the depth in the intensive green roofs. The choice of the vegetation is also different. Sedum, small grasses, herbs and flowering herbaceous plants which do not need permanent irrigation system are the typical vegetation in an extensive green roof.
During the last few decades, green roofs have been studied from different points of view. Several mathematical models of heat transfer or heat and mass transfer for green roofs have been proposed (see [4] for a review of this type of models). Two of the most representative green roof mathematical models were proposed by del Barrio [5] and Sailor [6]. These models are based on an energy balance. In [5], a set of relevant parameters in the design of green roofs were described, among them, the leaf area index (LAI) or soil density and thickness. The energy balance mathematical model considered in [6] was validated with experimental data of a green roof installed in Florida.
Moreover, several energy balance models describing the thermal behaviour of green roofs have appeared (e.g., [2]). Most of these energy balance models are one-dimensional models, that is, they consider the substratum and the canopy as homogeneous. See [4] for a review of green roof mathematical models.
In [3] experimental results about the influence of the choice of the vegetation species, the weather conditions and the watering on the cooling potential in prefabricated green roofs were obtained. Experimental data were used in some other works to validate green roof mathematical models (see e.g., [1,6]).
Most studies agree that evapotranspiration is one of the main factors that affect the thermal behaviour of green roofs (see [4,7] and the references therein). Evapotranspiration is a combination of soil evaporation and plant transpiration. Comparing the energy balance of a green roof to the energy balance of a conventional roof, we notice that evapotranspiration is not in the second one. Evapotranspiration effect plays an important role in Biosphere models (see [8]).
The purpose of our work is to present a mathematical model of the energy balance of a green roof which includes the shape of the green roof and its mathematical treatment. A system of two coupled nonlinear parabolic equations describing the temperature in a vegetation layer and the temperature in the substrate has been considered. The model includes the main feedback mechanisms of the energy balance than previous models consider, like evapotranspiration. The proposed model in this work considers that the properties of the substratum depend on the temperature, in fact, the equation modeling the evolution of the temperature of the substratum also consider the change of phase of water described by a maximal monotone graph, denoted by γ. In consequence, the heat capacity of the substratum is not constant, it is a function depending on the temperature. This is the main difficulty in the mathematical treatment of the problem. In the proof of the existence of solutions we have used a regularization of the maximal monotone graph to avoid the lack of regularity produced by γ. The space domain is a surface M of R 3 , then the obtained model is a two dimensional model which allow us to include the shape of the green roof in the space domain; formulation of partial differential equations on Riemannian manifolds is required. Although most of the green roofs are flat we can also find green roofs with different shapes, spherical, cylindrical, or non-standard shape adapted to the terrain. We mention two examples of buildings whose green roofs are not contained in a plane: California Academy of Sciences by the architect Renzo Piano and la Maison Vague (Wave House) by the architect Patrick Nadeau. The model will allow us to analyze the response of the temperature to different choices of shape of the green roof, choice of the vegetation or the influence of the local climate conditions, for example.
The article is organized as follows. In Section 2, we introduce the mathematical model and the required assumptions and simplifications. In Section 3, we study the system of PDEs from a mathematical point of view. After introducing the notion of weak formulation of the equations we proceed to prove the existence of solutions. Finally, Sections 4 and 5 are devoted to discussion and conclusions respectively.
Notice that due to the shape of γ, we do not expect a regular solution T s in the classical sense, i.e., T s ∈ C 1,2 t,x ((0, T) × M) because the variation of T s with respect to the heat presents a discontinuity at T s = 273 K due to the change of phase. Since, the existence of a measurable and continuous temperature is expected, as real experiments show, a new notion of solution has to be introduced. In such direction, it is natural to present the notion of weak solution and weak formulation in appropriate functional spaces. Following the original concept of distributions and generalized functions studied by L. Schwartz and the extensive literature existing from the 40s to our days, we introduce the notion of weak solution in Definition 1, Section 3. The weak formulation is an alternative to validate the model, that can not be interpreted in terms of classical solutions and it should be understood in terms of distributions. The model is based on an energy balance, hence the temperature T s is linked to energy. In particular, T s is related to the latent heat of fusion, which is implicitly defined by γ and can only be interpreted in terms of distributions.
To obtain the existence of solutions, we first regularize the problem by using a regularization of γ to obtain uniform bounds of M |∇T s | 2 dA independent of . Such estimates guarantee the pass to the limit in the weak formulation of the approximated problem to the limit one.

A Green Roof Mathematical Model
In this section we model the energy balance of a green roof in a given geometry, which is not necessarily flat. The model is based on an energy balance and it includes the interactions between these two layers: vegetation and substratum. The model is based on the proposed ones by del Barrio [5] and Sailor [6]. In our case, the space domain considered is a two-dimensional surface, M.
The model describes the evolution of the temperature in two different layers: the vegetable layer and the substratum layer, denoted by T v and T s respectively (temperatures are given in Kelvin).
The substratum is a porous medium. It is a mixture of solid (minerals and organic material), liquid (water) and gas (air and water steam). See De Vries [9] for thermal properties of soils. We consider the graph γ(u) in order to include the change of phase in the model and so, the heat capacity of the substratum depends on the temperature.
In the vegetable layer, the energy balance depends on the following mechanisms solar radiation absorbed by the vegetation, -long wave radiation exchange between vegetation and the exterior, -radiative exchange between the two layers (vegetation and ground), -the sensible flux between vegetation and the air surrounding the vegetation, -evapotranspiration, -conduction of the heat in every layer.
Following Sailor [6], the energy balance in the vegetable layer is given by the function where • σ v is the fractional vegetation coverage. • R av is the absorbed solar radiation (short wave) by the vegetal layer given by where t is the time and x ∈ M. Q is the solar constant, S is the isolation function (depending on the orientation of the surface). β v is the coalbedo in the vegetal layer, i.e., 1 − β v is the albedo function, that is, the fraction which is reflected.
s is the emissivity of the soil layer.
LAI is the leaf area index (see Frankenstein et al. [10]) defined as follows where A 0 and A 1 are positive constants in the following ranges If we consider small scales of time, the parameter LAI can be consider as constant.

•
H v is the sensible heat flux between the vegetal layer and the air (see Frankenstein et al. [10] and references therein). Following [10] we define H v as follows where e 0 windless exchange coefficient for sensible heat (2.0 W/m 2 ).
ρ a f is the air density in the foliage (kg/m 3 ) near the atmosphere-foliage interface. It is a given function and for simplicity we assume homogeneous in space, i.e., .
-C f bulk transfer coefficient; -T a f is the air temperature in the foliage approximated by and T a is the air temperature measured at the shelter (see Frankenstein et al. [10] for more details).
-W a f wind speed at the air/foliage interface (m/s). -c p,a specific heat of air at constant pressure.
• L v is the latent heat flux expressed in terms of the unknown temperature. We notice that L v depends on the stomatal resistance coefficient of the leaves of the vegetal layer.
where l v is the latent heat of vaporization, it is the amount of energy required to convert a unit mass of water to steam. It is measured in units of J/kg and is inversely proportional to the temperature. From Henderson-Sellers [11] it is estimated as follows (see also Sailor [6]) ρ a f , C f , W a f have been defined previously. -r is given byr and r 0 is a given constant (see [10] for more details). -q as and q f ,sat are the mixing ratio of the air within the canopy and the saturation mixing ratio at the foliage surface temperature. The explicit expression can be found in [6]. We assume that it is a given function.
For simplicity we assume The energy balance in the substratum is given by F 2 (see Sailor [6]) where • R as is the absorbed solar radiation by the substratum layer, and β s represents the coalbedo function of the substratum layer. • H s is the sensible heat flux where e 0 have been already defined. - The bulk transfer coefficient for sensible heat C g h is given as follows (see also [10]) In the model, C hg is assumed as a given function.
-T a f air temperature in the foliage approximated by (see Frankenstein et al. [10] and references therein) and T a is the air temperature measured at the shelter (known data).
-W a f wind speed at the air/foliage interface (m/s) -c p,a specific heat of air at constant pressure. Then Notice that C hg is assumed as a given function, therefore H s becomes • L s is the latent heat flux. Following Sailor [6] we have where -C e,s is the bulk transfer coefficient, it is analogous to C f , assumed constant. -l s is the latent heat of vaporization at the ground surface temperature, assumed constant.

-
W a, f has been previously defined as the wind speed at the air/foliage interface and it is a given datum. ρ a,s is the air density near the soil surface. It is a known datum.
q as and q s are the mixing ratio at the foliage-atmosphere interface and the mixing ratio at the ground surface. We assume they are known data. The explicit expression is presented in [6,10].
Then, L s is assumed a known function, i.e., Notice that F 1 and F 2 depend on the type of vegetation, in particular, they depend on two functions which characterizes the plants used in the green roof: Leaf Area Index (LAI) and Fractional vegetation coverage (σ v ). If we consider small scales of time, these functions (LAI and σ v ) can be included in the model as constant. We notice that they are two of the relevant parameters of the model.
One of the characteristics of the energy balance of the green roof is the solar shading produced by the foliage and the cooling by evapotranspiration (these two effects are not in the energy balance of the conventional roof). We consider some relevant parameters related with vegetation: leaf area index (LAI), foliage fractional coverage (σ v ), vegetation stomatal resistance (r s ) and the vegetation coalbedo β v . These parameters could depend on the vegetation type and the seasons. The analysis of the sensitivity of the model front fluctuations of these parameters would be a useful future research.
The model also considers the heat conduction in every layer, whose coefficients are given by k 1 , k 2 respectively. Considering the heat conduction and the energy balance, under suitable boundary and initial conditions, we arrive to following system of partial differential equations Here, the heat capacity in every layer could depend on the moisture.
γ is defined as the graph for γ 1 and γ 2 positive numbers. These constants could depend on the moisture. See the climate model studied in [12] for more details. See also [13]. Notice that F 1 and F 2 can be expressed as follows and

Well-Posed Problem
We start the mathematical treatment of the obtained model by analyzing the existence of solutions under the following assumptions: (H M ) M is a two-dimensional connected oriented Riemannian C ∞ manifold with compact closure and regular boundary ∂M.  (7) and (8), resp. whose coefficients a i = a i (t, x) (i = 2, .., 7) and b j = b j (t, x) (j = 2, .., 4) are uniformly bounded and a 0 , a 1 and b 1 are positive constants. (H k ) The heat conduction coefficients k 1 and k 2 satisfy and there exists a constants k 0 > 0 such that (H g ) The boundary conditions g v and g s satisfy and g v ≥ 0, g s ≥ 0.
(H 0 ) The initial data T 0v , T 0s ∈ H 1 (M) ∩ L ∞ (M), and satisfy the boundary condition Notice that divergence and gradient in the diffusion terms are understood in the sense of the Riemannian metric of M. In particular if k 1 ≡ 1 then the diffusion operator div(∇ M T v ) is the Laplace-Beltrami operator of T v (see [14]). For simplicity in the notation we denote ∇ M by ∇. See also [15] where an energy balance model on a Riemannian manifold was considered. For maximal monotone operators properties see [16]. Now, following [14] we define the functional spaces on manifolds,   See [14,15,17] for more details about functional spaces defined on manifolds. We now introduce the notion of weak solution.

Definition 1. We say that
We remark that although γ is multivalued, in order to solve the problem we introduce the following approximation of γ, denoted by γ for > 0 in the following way Notice that γ ∈ C 1 and there exists a positive constant c( ) such that |γ | < c( ) a.e. We denote the inverse of γ by ϕ , i.e., ϕ satisfies ϕ (γ (s)) = s.
We first obtain the following a-priori estimates.
Proof. We denote by P v and P s the following functions where k 0 > 0 satisfies, for |s| ≤ k 0 + 1 and P s (s, T s ) < 0, for any T s > k 0 , P v (s, T s ) > 0 for any T s < −k 0 . Since In similar way we obtain We introduce the truncation function τ defined by and consider the problem (10) where F 1 and F 2 are replaced byF 1 andF 2 Notice that, thanks to Mean Value Theorem, and their symmetric once we multiply by −(−T v − k) + and −(−T v − k) + . We multiply by (T v − k) + and −(−T v − k) + in the first equation of (10) and by (T s − k) + and −(−T s − k) + in the second equation to obtain, after integration by parts Then, Gronwall's Lemma and the choice of the truncation function end the proof.
Proof. We multiply the first equation in (10) by T v − w 1 and integrate by parts over M to get In view of Lemma 1 we have that After integration we obtain In the same way, we multiply by T s − w 2 in the second equation in (10). Then, after integration over M we have The subscript t represents the time derivative. We can express the first term by using for Φ : R → R defined as follows After integration over (0, T) and thanks to previous lemmas and assumption H g , we get, for which ends the proof.
Proof of Theorem 1. The proof is divided into two different steps. In the first step we prove the existence and uniqueness of solution to the approximated problem (10). In the second step we prove the convergence of the approximated solution to (10) to the solution of (5).
Step 1. We consider the equation satisfied by (T v , Γ ), with boundary conditions and initial data We construct a fixed point argument in the following way. Let with boundary conditions (15) and initial data (16).
Thanks to Lemmas 1-3 and Aubin-Lions-Simon Lemma, there exists a subsequence, (T v ,T s ) such thatT strongly in L p (0, T : L p (M)) and C(0, T : L p (M)) for any p < ∞ and weakly in L p (0, T : H 1 (M)) for any p < ∞.
Since any Z 1 ∈ γ(T s ) and Z 2 ∈ γ(T s − ) are uniformly bounded in L ∞ (M T ), there exist γ and γ such that In view of the boundedness in L ∞ (M T ) of γ and γ we have that there exists a subsequence of (denoted with the same sub index) such that and γ γ * (23) weakly in L p (M T ) for any p < ∞.
We now take limits in the weak formulation of problem (10) to conclude the result.

Discussion
A two-layers model has been proposed to analyze the thermal behaviour of green roofs. The model is based on an energy balance on the vegetation layer and the soil layer. The model incorporates the evapotranspiration, that is, the combination of soil evaporation and plants transpiration, this effect is one of the differences between a green roof and a conventional roof. The modeling has considered some simplifications in the reaction terms F 1 and F 2 , the case where the reaction terms are more general functions will be considered in the future including non-polynomial growth terms. The shape of the green roof is defined in a two dimensional spatial domain which leads to formulate the model as a system of partial differential equations on manifolds. One of the difficulties of the mathematical analysis of the model is the multivalued term γ(T s ) in the balance on the substratum. The multivalued term and the coefficient C v present a dependence of the humidity w, assumed known. Further development must consider w as unknown of a third equation, which will allow us to decrease the data collection.
We have proved the existence of solutions (temperature in every layer) in a suitable functional space and we have obtained some estimates of them. This mathematical treatment open new questions about this mathematical model. The model will allow to analyze the impact of fluctuations of the parameters, among them, the leaf area index (LAI), fractional vegetation coverage (σ v ).

Conclusions
In this work, we have proposed and analyzed a mathematical model of the thermal behavior of the green roof. The system consists of two coupled parabolic equations describing the temperature in a vegetation layer and the temperature in the substrate. Two features of this nonlinear model are remarkable in comparison with other previous models in the study of green roofs: the space domain which is a surface (and so, a two dimensional space variable) and the term γ to describe the effect of the change of phases in the substratum. The other terms in the energy balance model were expressed as polynomials with power less or equal to four.
The second result of the article is the mathematical proof of the existence of solutions which is presented in detail. Due to the shape of γ, classical solutions are not expected if the temperature falls bellow 273 K, therefore weak solutions have to be introduced.
Future work of this two-layer energy balance model will be the analysis of the relevant parameters (e.g., LAI), the inclusion of non-polynomial reaction terms and the numerical approximation of the solutions. Funding: The research of the first two authors is partially supported by the project MTM2017-83391-P Ministerio de Ciencia e Innovación, Spain. Second author is also supported by project MTM2017-85449-P, Ministerio de Ciencia e Innovación, Spain.

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