Unsteady Coupled Moisture and Heat Energy Transport through an Exterior Wall Covered with Vegetation

: A mathematical model that governs unsteady coupled moisture and heat energy transport through an exterior wall covered with vegetation is described. The unknown temperature and moisture content of the plants and canopy air are represented by a system of nonlinear ordinary differential equations (ODEs). The transport of moisture and heat through the support structure, which includes insulation and soil layers, is deﬁned in a series of nonlinear partial differential equations (PDEs). After setting out the model, this article presents and discusses a set of numerical applications. First, a simpliﬁed system consisting of a brick wall covered with climbing vegetation is used to study the role of individual variables (e.g., wind speed, minimum stomatal internal leaf resistance, leaf area index, and short-wave extinction coefﬁcient) on the hygrothermal behaviour of the green wall. Thereafter, more complex green wall systems comprising a bare concrete wall, mortar, cork-based insulation (ICB), soil and vegetation are used to evaluate the inﬂuence of the thermal insulation and substrate layers on the heat ﬂux distribution over time at the interior surface of the wall, and on the evolution of the relative humidity, water content, and temperature throughout the cross section of the green wall. The numerical experiments proved that vegetation can effectively reduce exterior facade surface temperatures, heat ﬂux through the building envelope and daily temperature ﬂuctuations.


Introduction
It is increasingly being acknowledged that urban sustainability requires a comfortable and healthy environment. In this context, vegetation plays an important part in meeting this expectation, dealing with many challenges related to biodiversity loss, air pollution, energy performance and climate change [1,2].
In addition, buildings are expected to be increasingly efficient, in particular with respect to heat exchange through the walls [3,4]. In recent years, designers have promoted the idea of incorporating plants into building envelopes to address ecological [5][6][7], aesthetic [8] and energy aspects [9][10][11]. Additional advantages of these solutions include their ability to reduce stormwater runoff [12][13][14], filter out harmful particulate matter and greenhouse gases [15,16], restore biodiversity [17,18], and mitigate the urban heat island effect [19][20][21]. Note, however, that whereas green roofs have been actively promoted from the standpoint of building engineering, green walls have received much less attention.
Green walls can involve a large variety of plants ranging from vines to shrubs, planted directly on the facade or supported on trellises and wires. Some green walls also have a layer of soil and are installed as an integral part of the building envelope [8,22]. There are two major categories of green wall technologies: Those where plants are rooted in the ground (indirect greening systems) and those where plants are rooted in artificial substrates or compost (direct greening systems). The indirect greening system uses the facade to train plants to grow upwards, while in the direct greening system, the plants are directly attached to the facade [23].
Although the thermal benefits provided by green walls are difficult to ascertain for certain, since many factors are involved in the phenomena that are involved in the overall behaviour (e.g., building design, plant species and local climate conditions), a green facade are more likely to affect the building envelope performance during relatively warm weather conditions, when the indoor temperature is lower than the facade's surface temperature and the heat flux is channelled indoors. A number of authors have tried to show the efficacy of green walls for insulating building facades and blocking solar radiation. For example, Lee and Jim [24] studied the shading effect of a green facade in Hong Kong, estimating an average daily energy saving of 0.226 kWh/m 2 . Perini et al. [25] and Coma et al. [26] also studied the thermal performance of green vertical systems in the Mediterranean climate; they found energy savings ranging from about 25% to almost 60% in the cooling season. Vox et al. [27] and Mazzali et al. [28] studied a number of living walls installed in Italy and found differences in the outer surface ranging from 9 • C to 20 • C between plant covered and bare walls. Arranz et al. [29] compared the thermal performance of south-and west-facing green walls in Madrid, concluding that the maximum temperature variation occurred for the south-facing wall, in both summer and winter.
The mechanisms underlying heat and mass transfer in green walls are complex. Broadly speaking, plants reflect and transmit long-wave radiation, while absorbing the short-wave radiation during photosynthesis. However, the capacity to absorb radiation depends on many variables (e.g., water content, plant properties, and radiation wavelengths). Plants also exchange heat with surrounding air by transpiration and convection, and with the substrate by conduction. Wind speed and moisture content of the air and plant species alike play an important role in these processes [11,[30][31][32].
While many studies report on the hygrothermal behaviour of green walls, very few of them are carried out under environmental conditions that can be easily recreated, making rationalization of results very difficult. Thus, this paper sets out to present a numerical model that can simulate the dynamic heat and moisture transport behaviour of real-world green walls that consist of the canopy and the structural support. The numerical tool governing the transport phenomena in green roofs previously proposed by the authors [11] has been modified to be able to simulate the transport of heat and moisture in green vertical walls. The green wall model estimates the effects of plants/foliage on lowering the surface temperature of the facade and the heat/moisture flux through the outer walls of buildings. A canopy layer is a very complex physical-biological system involving heat and moisture transport; this means that it is extremely hard to develop a detailed mathematical model. We have taken a canopy to be a single homogeneous layer comprising a system of concentrated parameters with one value for leaf temperature and one value for the temperature and moisture content of the canopy air [11,23,33,34]. The canopy is therefore represented by a system of nonlinear ordinary differential equations (ODEs) where the temperature of the plants that is not known, and the temperature and moisture content of the canopy air where it is also not known.
The structural support, including the bare wall, insulation layer and soil, was regarded as an unsaturated homogeneous porous medium for which the heat and moisture transport are interlinked and coupled. It is governed quantitatively by a set of coupled nonlinear partial differential equations (PDEs). Many researchers have studied coupled heat and moisture transfer in a rigid porous medium under temperature and moisture gradients [35][36][37][38].
We used the finite difference model [39] to solve the ODEs and the PDEs were discretized using the boundary element model [40]. Numerical solutions of the green wall model are strongly related to the plant parameters that define the thermal energy transport, mass transport and storage behaviour, empirical relations involved in modelling heat and moisture fluxes, crop transpiration, choosing the relevant different stems of the plants, and characteristics of the canopy, etc. In general, the empiricism applied in this contribution was taken from [33,[41][42][43]. The values of the parameters vary significantly and have a clear impact on the numerical solutions.
First, we set out the governing transport equations for the canopy and for the porous solid wall. It goes on to describe the finite difference and the boundary element methods used to solve the canopy and porous solid transport equations, respectively. Based on this model, three numerical examples of increasing complexity are discussed: a two-layer greenery system (canopy-brick wall), a four-layer greenery system (canopy-ICB-mortar-concrete wall), and a five-layer greenery system (canopy-substrate-ICB-mortar-concrete wall). Results are further discussed in terms of the effect of individual variables (e.g., wind speed, minimum stomatal internal leaf resistance, leaf area index, and short-wave extinction coefficient) on the hygrothermal behaviour of the green wall. The role of the thermal insulation and substrate in the heat flux distribution at the interior surface of the wall, as well as in the relative humidity, water content, and temperature throughout the cross section of the green wall is also discussed. The main conclusions are summarized at the end.

Mathematical Model for the Coupled Moisture and Heat Transport in a Canopy
The canopy layer (c) consists of the plants and their leaves (p) and the canopy air (ca) within the leaf cover. Defining the heat and mass transport phenomenon in the canopy layer is very difficult, making an exact mathematical physical description very nearly impossible. In general, the canopy is treated as a homogeneous layer of concentrated parameters. It is bounded on one side by the facade surface (sw) and on the other by an ideal air ambient surrounding (a).
The mathematical model governing the non-isothermal heat energy and moisture conservation in the canopy layer is given by a set of three coupled conservation ordinarydifferential equations [32][33][34]41,42,[44][45][46][47]. These are the heat energy conservation equation for the plant and the heat energy and moisture mass conservation equations for the canopy air, formulated for the temperature of the plant T p (t), temperature of the canopy air T ca (t) and the canopy air vapour pressure p v,ca (t) as the primitive driving potentials: These equations describe the heat energy accumulation within the plant, and the heat energy and moisture accumulation within the canopy air, respectively, due to physically different heat and mass fluxes, e.g., convective (conv), shortwave radiation (rs), longwave radiation (rl) and transpiration (trans), or basically due to different sources and sinks of heat energy and mass in the greenery system. The quantities c p = (ρc p ) p , δ p , and LAI in the plant accumulation term indicate the effective specific heat per unit volume, the average leaf thickness, and the leaf area index, which is the ratio of the total leaf top surface area to the facade surface. The quantities c ca = (ρc p ) ca , δ c , ρ ca , and θ p are the specific heat per unit volume of the canopy air, the average canopy thickness, the canopy air mass density and θ p = dω ca /dp v,ca , respectively, where ω ca stands for the specific humidity of the canopy air. The subscripts of the heat energy and mass fluxes refer to heat or mass exchange between two spatial positions, e.g., the subscript (ca − p) in the convective flux expression q conv ca−p denotes the heat energy exchange between the canopy air (ca) and the plant (p), and the same applies to the other flux expressions in Equations (1)-(3).
A canopy is such a complex living, physical-biological system that it is not possible to make any kind of analytical prediction of the heat and mass fluxes involved in Equations (1)- (3). The characteristics of the greenery system are the intrinsic inhomogene-ity of the foliage layer and the turbulent nature of the air stream within and around a canopy [34]. The only reasonable approach to estimate heat and mass fluxes inside the canopy layer, between the canopy layer and the ambient, and between the canopy layer and the facade surface, is to use empirical correlations based on the Newton constitutive heat and mass transfer model and linearized Stefan-Boltzmann radiation law.
The radiation heat flux terms in Equation (1), e.g., the short wave solar radiation absorbed by the plant q rs p , longwave radiation heat exchange with the sky q rl sky−p , with the ground/surrounding surfaces q rl gro−p , and with the facade surface q rl sw−p , such that q rl p = q rl sky−p + q rl gro−p + q rl sw−p , can be estimated by using empirical rheological models: where q sol represents the solar radiation at the top of the canopy and the quantities ρ r p , ρ r ∞ and ρ r sw are, respectively, the shortwave reflectance of the leaves, a dense canopy, and a facade surface; σ f is a fractional vegetation coverage, whilst T sky , T gro and T sw are the sky temperature, the temperature of the ground, assumed to be the same as the outside air temperature T a , and the exterior facade surface temperature. The linearized longwave radiation heat transfer coefficients are given as [11,32,34,47]: The quantities τ s and τ l are the shortwave and longwave transmission coefficients of the greenery system calculated from LAI via an extinction law in a turbid medium, characterized by a shortwave k s and longwave k l radiation extinction coefficients, k s ≈ 0.74k l , depending on the geometric characteristics of the foliage [32]: τ s = exp(−k s LAI) and τ l = exp(−k l LAI).
The shortwave extinction coefficient k s indicates a drop in the absorbed solar radiation in the canopy. It ranges between k s = 0-1, where lower values, k s = 0.3-0.5, relate to leaves with leaf angles less than 45 • and higher values, k s = 0.7-1, relate to leaves with leaf angles greater than 45 • [34]. When the leaves are perpendicular to the wall, k s ≈ 0, and when they are parallel to the wall, k s ≈ 1. The effective emissivity ε sw−p between the plant and the exterior facade surface is expressed as: where the quantities ε p and ε sw are the leaves and exterior facade surface emissivity, respectively, and σ is the Stefan-Boltzmann constant. The geometric relationships for longwave radiative exchange between radiation sources and radiation receivers are described by the view factor F (0-1). In general, the view factor indicates what proportion of the longwave radiation leaving an object with a particular shape is intercepted by an object with a similar or different shape. So the view factors to the sky F sky and to the ground F gro can be expressed as: where θ is the tilt angle of the facade surface in relation to the ground, e.g., θ = 90 • for a vertical wall and θ = 0 • for a roof. For the vertical green wall, with a tilt angle θ = 90 • , the view factors are equal F sky = F gro = 1/2 [34].
The convective and transpiration heat and mass flux terms in Equations (1)-(3), e.g., the sum of the heat fluxes for the plant q p = q conv ca−p + q conv a−p − q trans p−ca , the sum of the heat fluxes for the canopy air q ca = q conv sw−ca + q conv a−ca + q conv p−ca , and the sum of the moisture mass fluxes for the canopy air j v,ca = j trans v,p−ca + j conv v,sw−ca + j conv v,a−ca , can be estimated using an empirical constitutive Newton type transfer model: where p v,p is the vapour pressure at the leaf tissues equal to saturated vapour pressure (2) in the equations accounts for the lower and upper leaf surfaces. The heat transfer coefficients α ca−p , α a−p = α a−ca , α sw−ca and α trans p−ca can be formulated as follows: where the quantities r e , r i , r ea , r sw−ca and r s express heat transfer aerodynamic resistances, e.g., the external single leaf resistance, internal/stomatal single leaf resistance, resistance to the surrounding environment, the facade surface to canopy air resistance, and additional moisture-dependent soil surface resistance, respectively, whilst c a = (ρc p ) a and γ refer to the specific heat per unit volume of the ambient air and thermodynamic psychometric constant. The vapour mass transfer coefficients β sw−ca , β a−ca and β trans p−ca can be formulated as: where the quantity h e (T) is the specific latent enthalpy of water evaporation.

Empirical Correlations for Heat/Mass Transfer
The transport phenomena in the greenery system can be estimated in terms of nondimensional criterion numbers, e.g., Reynolds number Re, Grashoff number Gr and Prandtl number Pr, based on some kind of similarity between the fluid flow along the leaf surfaces in the greenery system and the simplified geometrical cases of developed external forced and free flow along the flat plates and vertical walls. Empirical corelations for estimating the magnitude order of resistance to heat and mass transfer were found in the literature [33,34,41,43,47,48] and are given by the expressions indicated below.
In general, according to the ratio of non-dimensional criterion numbers Gr/Re 2 the heat/mass transfer regime can be determined. For example, the forced convection occurs for Gr/Re 2 << 1, mixed convection for Gr/Re 2 ≈ 1, and the free convection occurs for Gr/Re 2 >> 1, where both non-dimensional Reynolds Re and Grashoff Gr numbers are based on the characteristic length of the leaf L ch = (L lea f · W lea f ) 0.5 , where the quantities L lea f and W lea f stand for the leaf length and width. The wind speed v is an average for the height of the plant usually estimated from the wind profile. Transport properties of the air, e.g., heat conductivity λ, mass density ρ, specific heat c p and kinematic viscosity ν, should all be based on the mean boundary layer temperature, e.g., T m = (T p + T a )/2.
Heat transfer coefficients between the ambient air and the foliage/plant α a−p = c a /r ea and between ambient air and the canopy air α a−ca = c a /r ea can be estimated according to the following empirical correlations, e.g., for the natural convection with the values (C = 0.63, n = 1/4) for the laminar boundary layer (1 × 10 4 < PrGr < 1 × 10 9 ), while in the case of turbulent boundary layer (1 × 10 9 < PrGr < 1 × 10 13 ) the values for the constants are (C = 0.15, n = 1/3). When the forced convection is dominant over the heat exchange the following empirical relation can be used: where the values for the constants C, m and n are valid in the case of the laminar flow (C = 0.664, m = 1/3, n = 1/2), while the values (C = 0.037, m = 0.43, n = 0.8) are valid for the turbulent flow regime Re > 1 × 10 5 [33,49]. Mixed convection occurs when the forced and natural convection overlap, in which case the heat transfer coefficient can be estimated by using formula [50]: where Nu f c is the Nusselt number of the forced convection and Nu nc is that of the natural convection, while the exponent n is generally n = 3. The heat/mass exchange between the leaves and the canopy air through the boundary layer formed on the leaf surface is expressed through the aerodynamic resistance r e , which depends on wind speed and surface roughness. The aerodynamic resistance of a single leaf surface (one side) r e is defined as: where the Nusselt number is estimated from Equations (17) and (18). Transpiration defines the mass vapour flux or latent heat flux from a plant to the canopy air. Some leaves transpire through both sides, while some leaves transpire through only one side. An acceptable phenomenological model predicting the canopy single leaf stomatal resistance of a fullyvegetated green wall r i can be adopted: where r i,min is the plant's characteristic minimum stomatal resistance for the particular crop. F i represents functions greater than unity, which describe the relative increase of stomatal resistance if any of the applicable variables of the surrounding environment (shortwave irradiation q sol , leaf temperature T p , leaf-to-air vapour pressure p v,p − p v,ca , carbon dioxide of the ambient air CO 2 ) is affecting the vapour transfer rate. These functions are chosen: where q sol is the mean flux density per unit leaf area, T m = 24.5 • C is the temperature at which resistance is minimal, CO 2 = 200 vpm is the volume concentration for which resistance is minimal, and the empirical coefficients are Typical minimum stomatal resistance values ranges from 80 to 160 s/m [34]. It is evident that the estimation of the stomatal resistance r i is very complex. At night, r i dramatically increases, e.g., described by Equation (22), thus transpiration is reduced to a minimum value (5-15% of daytime values) [47]. The wall surface to canopy air aerodynamic resistance r sw−ca is defined as Equation (20): where the Nusselt number is estimated from Equations (17) and (18). Additional soil surface resistance r s to water vapour flow can be given as: where the minimum surface resistance r s1 = 10 sm −1 , the minimum soil volume moisture content Θ min = 0.15, while the diffusion coefficient α s = 35.63.

Mathematical Model for the Coupled Moisture and Heat Transport through an Unsaturated Porous Solid Wall
A system of two coupled partial-differential equations governing the general case of non-isothermal moisture and heat energy transport through an unsaturated porous solid wall can be formulated as [11,35,36], i.e., written for the continuous field functions relative humidity ϕ(r j , t) and temperature T(r j , t) as the primitive diffusion driving potentials where θ = dW/dϕ is the slope of the sorption isotherm W = W(ϕ), and the quantities W, p v , ρ l and g represent the mass moisture content, vapour pressure, liquid water mass density and the gravity acceleration. The sink termṁ RU defines the mass of water removed from a unit volume of soil per unit of time by plant water uptake [51][52][53][54]. The transport coefficients D ϕ and D T are given as where the transport coefficients δ p and D l stand for the vapour and liquid permeability of the porous solid, respectively, whilst the quantities p s (T) and R w are the saturation pressure and the water vapour gas constant.
The quantity h lat = [h e + (c pv − c pl )T] in Equation (26) denotes specific latent enthalpy, h e is the specific latent enthalpy of evaporation or condensation and the amounts c pv , c pl denote the specific heat per unit mass of water vapour and of liquid water, respectively. The effective transport coefficients of the porous solid, e.g., the effective specific capacity per unit volume c e f f = (ρc p ) e f f and the effective thermal conductivity λ e f f are moisture content dependent variables.

Canopy-Facade Surface Mathematical Closure Model
The canopy-facade surface mathematical closure model represents the coupling of the mathematical models governing the heat and moisture transport in the greenery system and in the porous solid wall, appearing in the form of mixed or Cauchy type boundary conditions on the facade surface.
The heat and vapour mass balances on the facade surface or the continuity requirements of the heat and vapour fluxes at the canopy-facade interface can be formulated as [11]: where the normal total heat flux q sw = ( q · n) sw = [( q sens + q lat ) · n] sw and moisture flux j v,sw = ( j v · n) sw at the facade surface, defined with the unit normal vector n sw , flowing into the porous solid wall must be equal to the total heat and mass inflow from the surroundings, e.g., canopy, sky, and ground or canopy air, respectively, to the facade surface. The continuity requirements of the driving potentials at the canopy-facade surface interface imply: T sw (t) = T(t, sw) and p v,sw (t) = p v (t, sw). The convective heat flux q conv ca−sw , shortwave radiation heat flux received at the facade surface q rs sw , e.g., the fraction which is transmitted through the greenery system minus the part reflected from the facade surface, and longwave radiation flux term q rl sw , can be formulated as: where the heat transfer coefficient α ca−sw = α sw−ca and Σ f = 1 − σ f , whilst the linearized longwave radiation heat transfer coefficients can be given as and the following equality is valid, α rl p−sw = α rl sw−p .

Computational Solution Numerical Models
The analytical solution of systems of coupled nonlinear ordinary differential equations (ODEs) and partial differential equations (PDEs) given by Equations (1)-(3) and Equations (25) and (26), respectively, governing the heat and moisture transport through the green wall, is not foreseen, therefore, it is necessary to use approximate numerical methods to estimate heat fluxes through the green wall, e.g., finite difference method (FDM), finite volume model (FVM), finite element method (FEM) and boundary element method (BEM).
The concept of coupling the finite difference method (FDM) to solve the set of ordinary differential equations (ODEs) governing the transport phenomena through the greenery system and the boundary element method (BEM) to solve the set of partial differential equations (PDEs) modelling the transport phenomena through the porous solid, was successfully developed and applied in a simulation of heat and moisture circumstances through the green roofs [11]. Therefore, this solution approach is briefly summarized with some modification due to additional flux terms.
In order to generalize the finite difference solution approach, the set of ordinary differential equations (ODEs) formulated for the primitive variables, e.g., temperature and vapour pressure for the canopy, can be formulated generally in a form of conservation ordinary differential equation written for the field function u(t): where the first term denotes the accretion of an extensive property available per unit surface area, such as heat energy, moisture content, whilst the second term takes into account (n) flux terms in Equations (1)-(3) of the appropriate extensive property, where the amounts α(u) and u a denote a nonlinear transfer coefficient and respective surrounding ambient potential. The Euler first order two-time level scheme or three-time level second order asymmetric difference formula can be applied on a time-axis (v) [39], such as: where ∆t = t v+1 − t v being the time increment, yielding the following explicit finite difference approximation for the potential corresponding to Equation (34) where, for instance, A = 1/∆t and A = 3/2∆t for the first and second finite difference scheme, respectively. Note that the FDM is used only to solve the ODEs governing the transport phenomena through the canopy, using an explicit finite difference approximation.
Although the canopy physics is very complex with a lot of empiricism, the solution of the final discretized equations does not represent more significant difficulties than the numerical solution of the porous media. In fact, both problems are coupled through the Cauchy-type boundary conditions. This approach of coupling the two models proved to be stable and accurate.
The system of partial differential equations (PDEs) governing the heat energy and moisture transport through the porous solid Equations (25) and (26) can be unified as a single generic transport equation for the field function u r j , t : where the notation L[ · ] represents the parabolic diffusion linear operator, and the terms b j r j , t and b r j , t represent sources due to the material nonlinearity and production of the conservative field function, respectively, with the following corresponding integral representation written for a time increment ∆t = t 2 − t 1 [11,40,55]: where q = q j n j = ∂T/∂n is the field function normal flux, and the notation u stands for the parabolic diffusion fundamental solution. Following the basic concept of boundary element solution strategy, the boundary Γ and the domain Ω in Equation (38) are discretized into a series of boundary elements and internal cells. Furthermore, the variation of all functions depends on interpolation functions over space and time [11], e.g., high order elements were applied with quadratic variation over space and the linear variation over time of all functions. When dealing with nonlinear transport problems, and to cut the storage and CPU time requirements, the subdomain technique/macro element model should be used [55], which results in sparse and diagonal block-banded influence matrices.

Two-Layer Greenery System (Canopy-Brick Wall Model)
This example represents a single brick wall without thermal insulation, where the vegetation climbs along its surface, with the roots on the ground. This is the case of several old buildings built without hygrothermal care.
This green wall example, presented in Figure 1, considers the heat and moisture transport in a coupled brick wall-canopy structure 0.45 m thick (δ = δ b + δ c , with δ b = 0.2 m and δ c = 0.25 m) and 2.2 m in height (L w ). The system was exposed to periodic fluctuations of the prevailing ambient parameters, such as the solar heat flux Q sol (t) and the temperature T a (t), given by Equations (40) and (41), while the ambient vapour pressure p v,a = p v,a (T a,o , ϕ a,o ) was kept constant.
while on the right/outdoor side the boundary conditions are estimated by the empirical correlations given by Equations (17)- (24). The time-dependent analysis involved running the simulation from the initial state with a time-step value ∆t = 1800 s. The solar heat flux was assumed to be zero at night and to increase gradually during the day, as represented by this sinusoidal expression: where ω = 2π/τ, τ = 24 h, q o = 600 W/m 2 , whilst the external ambient temperature T a variation is given by the relation: with the values T a,o = 22 The solid porous brick wall was reasonably thick, δ b = 0.2 m, to enable study of the sensitivity of different canopy parameters with an influence on the hygrothermal behaviour of a green wall, with special emphasis on the facade surface temperature and indoor heat flux. The sensitivity of the presented greenery system/canopy was analysed by changing the values of one influential parameter at a time while keeping other influential parameters constant. The effect of the following influential parameters was studied: wind speed, v The objective was to calculate the moisture and temperature distribution in the greenery system and in the brick wall over time sequence of t = 240 h (the results are only plotted from t = 120 h to t = 240 h to allow an easy interpretation of the graphics). A uniform dense mesh of M = 100 × 1 macro-elements was used to model the coupled heat and moisture transport phenomenon through the brick layer. The convergence tests concern the space and time discretization of the solution domain, to check the validity of the proposed BEM numerical model to solve PDEs. They were conducted rigorously and are presented in previous works on the transport phenomena through multilayer porous walls [11,36,37]. The model was tested on several well documented benchmark examples to demonstrate its accuracy and stability.
The required data defining the reference greenery system are summarised in Table 1, while the hygrothermal properties for the porous solid brick layer are taken from the literature [37,52].
Using the presented numerical model, the sensitivity of natural, mixed and forced convection heat exchange from the bare brick and vegetated brick wall surface and the ambient to different wind velocities, v    For the tested modelled range in wind velocities, the maximum bare facade surface temperatures, T sw , were 56.0 • C, 55.9 • C, 52.6 • C, 49.8 • C, 47.6 • C, and 45.9 • C, whilst the corresponding indoor heat flux values were in the range from q int = 54.9 W/m 2 to q int = 36.8 W/m 2 (see Table 2). The numerical simulation results for the vegetated facade are depicted in Figures 4 and 5. It is evident that the foliage had beneficial effects on the hygrothermal behaviour of the facade by lowering the maximum temperatures, T sw (35.7 • C, 32.5 • C, 31.5 • C, 31.2 • C, 31.1 • C, 31.0 • C), and consequently the corresponding heat load values on the indoor wall were significantly reduced in the range from q int = 21.6 W/m 2 to q int = 14.2 W/m 2 (see Table 2).  Numerous canopy parameters have significant impacts on the hydrothermal performance of the vegetated facades, including minimum stomatal resistance (r i,min ), leaf area index (LAI), shortwave extinction coefficient (k s ), etc.
The cooling effect of the transpiration heat flux from the leaves was evaluated by changing the values of a minimum stomatal resistance r i,min in the range (80 s/m, 120 s/m, 160 s/m, 200 s/m, 240 s/m, 280 s/m), shown in Figures 6 and 7. The biggest difference in facade surface temperature for the maximum and minimum resistance, r i,min , is approximately ∆T sw ≈ 1.8 • C, with the corresponding approximate difference in indoor heat flux ∆q int ≈ 3.2 W/m 2 .  Besides the wind velocity and the ambient temperature, shortwave solar radiation is one of the most crucial weather parameters that affects the behaviour of green walls. A canopy absorbs a fraction of incident solar radiation and transmits the rest through the greenery system. The transmittance of the solar radiation through the canopy exponentially decreases with foliage density, given by Equation (9), where the attenuating effect is expressed through an extinction/attenuation coefficient, k s , and a leaf area index, LAI. Figures 8-11 show the impact of various LAI and k s values on the hygrothermal behaviour of the vegetated wall. Increasing the leaf area index LAI values and the extinction coefficient k s values, the facade surface temperature and the heat flux through the structure porous wall decrease significantly. A dense canopy with a LAI > 3 and with k s > 0.7 (leaves parallel to the wall) is particularly effective in reducing the facade surface temperature T sw and heat flux load q int .

Four-Layer Greenery System (Canopy-ICB-Mortar-Concrete Wall Model)
This example represents a single concrete wall, initially built without thermal insulation, but later coated with a layer of insulation after the interposition of a mortar water-proofing layer. As can be seen in the example above, the vegetation climbs along it with the roots on the ground.
The realistic direct green wall system, a schematic representation of which is depicted in Figure 12, considers the heat and moisture transport in a four-layer canopy-ICB-mortarconcrete structure, where the thickness of the individual layers are: canopy, δ c (0.20 m); ICB, δ ICB (0.00 m, 0.05 m, 0.10 m, 0.15 m, and 0.20 m); mortar, δ m (0.01 m); and concrete, δ co (0.20 m). This system was exposed to the periodic fluctuation of the prevailing weather conditions given by Equations (40) and (41). The time-dependent simulation was performed by running the computation from the initial thermodynamic state with a time-step value ∆t = 1800 s. The simulation time was 240 h (however, the results are only plotted from t = 120 h to t = 240 h to allow an easy interpretation of the plots). Two non-uniform density macro-element meshes of M = 45 × 1 and M = 75 × 1 were applied to numerically solve a highly nonlinear coupled transport of heat and moisture through the multilayered system. The first (coarser) one adequately captured the transport phenomena in the systems with the canopy, while for the cases without the canopy, only the finer mesh was good enough, since the nonlinearity in moisture transport became more severe.
The influence of the canopy on the multilayered system with two different thicknesses of the medium density expanded cork (ICB), namely δ ICB = 0.10 m and δ ICB = 0.15 m, is depicted in Figure 13. It is evident that the canopy reduces the heat load on the interior side in both cases of insulation thickness, although the increasing thickness of the ICB insulation material has a greater impact. We can say that in the multilayer greenery system the influence of the insulation as well as the canopy on reducing the thermal load of the interior is significant. However, it is also clearly shown that the canopy acts mainly as the solar radiation barrier, namely by increasing the shortwave radiation extinction coefficient from the value k s = 0.6 to the value k s = 0.9, such that the reduction of heat load on the interior is most noticeable.  Figure 14 illustrates how the driving potentials, such as relative humidity ϕ(t, x), water content W(t, x), and temperature T(t, x), change with time (t = 216, 222, 228 and 234 h) throughout the multilayer wall cross section. The bare and the vegetated facades were analysed for the δ ICB = 0.10 m thick ICB insulation layer. During sun exposure, the bare facade exterior surface temperature rose as high as T ≈ 54 • C, whilst the maximum ambient temperature was T = 33 • C; the solar heat radiation is absorbed and stored by the exterior ICB surface. Due to very low heat diffusivity of the ICB material, the propagation of heat from the surface to the interior is very slow, resulting in a high surface temperature. Large surface temperature fluctuations result at the same time in large fluctuations of surface relative humidity. The slow moisture diffusion process means the fluctuations are limited to the wall surface. In fact, the exterior surface temperature of the vegetated facade is much lower, T ≈ 45 • C for K s = 0.6 and T ≈ 34 • C for K s = 0.9, suggesting that the choice of suitable plants for the canopy is of crucial importance. It can be also noted that the computed difference in the surface night temperature between the bare and vegetated facade is negligible. Heat convection is the dominant heat transfer mechanism and the ambient temperature controlled the facade surface temperature. The longwave heat radiation from the bare facade is slightly more intense, which resulted in a lower night surface temperature.  Figures 15 and 16 give the time variations of the heat load at the inner surface of the wall and temperature at the external surface. It is evident that the system without the ICB layer leads to the highest heat load of the interior. Simultaneously, the daytime external surface temperature is the lowest due to the high heat diffusivity coefficients of mortar and concrete.

Five-Layer Greenery System (Canopy-Substrate-ICB-Mortar-Concrete System)
As before, this example represents a single concrete wall, initially built without thermal insulation, but later given an insulating layer, after the placement of a mortar waterproofing layer. However, the vegetation now grows on a substrate layer superimposed to the insulating layer.
This green wall, represented in Figure 17, considers the heat and moisture transport in a five-layer system composed of canopy-substrate-ICB-mortar-concrete: canopy, δ c (0.20 m); substrate, δ s (0.10 m); ICB, δ ICB (0.00 m, 0.05 m, 0.10 m, 0.15 m, and 0.20 m); mortar, δ m (0.01 m); and concrete, δ co (0.20 m). The system was exposed to the periodic fluctuation of the prevailing ambient parameters. Boundary and initial conditions were taken to be the same as those previously prescribed by Equations (39)- (41), with the exception that the initial conditions for the substrate relative humidity were equated to value ϕ = 0.95. The substrate/clay loam soil thermohydraulic properties were taken from the literature [33]. Numerous time-dependent analyses of transport phenomena in a multilayer structure were performed by running the simulations from the initial state with a time-step value ∆t = 1800 s; the simulation time was 240 h (as before, the results are only plotted from t = 120 h to t = 240 h to allow an easy interpretation of the plots). Non-uniform density macro-element meshes of M = 75 × 1, M = 100 × 1 and M = 115 × 1 were applied to numerically model a highly nonlinear coupled transport of heat and moisture through the canopy, substrate, ICB, mortar and concrete layer; to cope with the severe nonlinearity of the outer layer/substrate 50 macro-elements were used. Figures 18 and 19 present the time variations of the heat flux at the interior surface of the wall and temperature at the external surface as a function of ICB thickness; the thickness of the substrate layer was kept constant at δ s = 0.10 m. Comparing these results with those shown in Figures 15 and 16, it can be concluded that a substrate layer considerably reduces heat flux through the multilayer wall and the facade surface temperature. Overall, the maximum surface temperature of the vegetated facade is T sw ≈ 38 • C, regardless of ICB insulation layer thickness. This is because of the favourable thermohydraulic properties of the soil in terms of heat diffusion when compared to ICB.

Conclusions
In this contribution, the numerical simulation tool governing the transport phenomena in green roofs [11] was modified in such a way as to be able to estimate the transport of heat and moisture in vegetated vertical walls. The proposed numerical solution procedure deals with the highly nonlinear coupled heat and moisture transport through a multilayer structure composed of a canopy and a support structure. The relevant physical-mathematical models governing the transfer of heat and moisture through the canopy and porous solid media have been discussed. The required empiricism to close the mathematical models has been considered, too. The respective boundary and interface conditions, and the coupling between the canopy and the facade surface were then formulated.
A standard first and second order finite difference FDM scheme was used to numerically solve the ODEs governing the transport phenomenon through the canopy layer, whilst the PDEs that govern the transport phenomenon through the porous solid material layers are solved with the boundary element numerical model (BEM). The finite difference and singular integral representations of the respective ODEs and PDEs were considered. The basic idea of transforming the integral equations to the boundary element numerical model has also been presented briefly. We have considered a very accurate numerical scheme based on quadratic approximations of all relevant field functions over the space domain, and linear over time.
The presented numerical simulation tool can be used to quantify the dependence of the heat load on various influential parameters of the vegetated systems, such as the parameters defined for the canopy and the structural parameters of the multilayer wall.
Several test multilayer wall cases were analysed to prove the capability and the sensitivity of the developed model.
Based on the analysis of the simulation results presented in this contribution and results taken from the literature, some conclusions can be established. The local environmental conditions, i.e., wind velocity, temperature and shortwave solar radiation, are the most critical influential parameters. The density of the canopy (LAI) and the orientation of the leaves also have a dominant influence on the transmittance of solar radiation, and thus on the facade surface temperature and heat flux load. The cooling effect of transpiration, governed by the stomatal internal leaf resistance, is noticeable as well. Additionally, the influence of different parameters is more pronounced during the day than at night.
Ultimately, the developed model is a compromise engineering tool to assess very complex transport phenomena through the greenery system based on the extended empiricism, which should be further tested and improved. Funding: This research has been supported by the Project EvaporCork (POCI-01-0247-FEDER-033357/ALG-01-0247-FEDER-033357) funded by the Operational Program for Competitiveness and Internationalization (POCI) of Portugal 2020 and the European Regional Development Fund (FEDER).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are partly available on request from the corresponding author.