An Evapotranspiration Evolution Model as a Function of Meteorological Variables: A CFD Model Approach

The study of meteorological variables and evapotranspiration in open spaces using the three spatial dimensions represents a technical challenge since the high computational resources required only enable the problem to be addressed on a very small scale. This research sets forth a three-dimensional computational fluid dynamics numerical model, characterized by its simplicity, which allows problems to be addressed over large areas (scale of kilometres). Similarly, the corresponding design and software developments carried out allow for a more dynamic introduction of meteorological and evapotranspiration boundary conditions. In the numerical domain created, the Reynolds-Averaged Navier–Stokes equations are solved, supplemented by a multispecies model (to distinguish the movements of dry air, evapotranspiration and air humidity) and one of solar radiation. The numerical model was applied to a semi-arid area in southern Spain, obtaining the three-dimensional special evolution of evapotranspiration, temperature, air humidity and wind velocity, specifically concerning its variation in horizontal and vertical planes of the three-dimensional domain, as well as vertical profiles at discrete points.


Introduction
Climate change not only impacts countries' economies, but it also affects human health. The Paris Agreement [1] focuses attention on the greatest challenge for the future and provides a catalyst for research lines directed towards climate and meteorology.
The global climate crisis has been closely linked to issues such as global changes in the water cycle, deterioration of urban air quality or desertification. Therefore, the study of the relationship between meteorological variables and climate change, especially their spatial distribution in large-scale open areas, is essential.
So far there have been several studies which have analysed the spatial evolution of meteorological variables in indoor or outdoor spaces.
In the case of indoor spaces (such as greenhouses) the existence of boundaries and their reduced dimensions simplify the study (simulation and validation) of the evolution of the meteorological variables. Such simplicity has allowed for the undertaking of different studies using Computational Fluid Dynamics (CFD) models with reduced computational resources. Investigations have also been held in greenhouses, using three-dimensional (3D) CFD models combined with those of radiation, to analyse the distribution of both temperature and humidity [2] or the CO 2 distribution [3] and other cases with 2D CFD models to analyse the vertical distribution of humidity and temperature [4] or even the humidity condensation in a greenhouse roof [5].
Another research line in indoor spaces is focused on evaluating human thermal comfort using CFD models, studying indoor temperature and air flow distributions. For example, using a 2D model in a house ventilated with a solar chimney [6]; with a 3D model of test rooms to evaluate the efficiency of radiant cooling systems [7]; or a 3D model of a house which enables simulation over long periods of time, demonstrating its accuracy compared with real data [8].
In the case of outdoor spaces, the investigations are focused on different scenarios using 3D CFD models to study distributions of airflow, temperature, humidity, radiation, and evapotranspiration (ET). For example, using a model of: A house with a tree located nearby [9]; street canyons [10]; multiple street canyons [11]; open spaces like a stadium [12] or a dam [13].
All these models are located on a micro scale due to the high computational cost of achieving great precision in large-scale models. While there are simulations with high spatial resolution on large scales, it is nevertheless necessary to use supercomputing centres [14]. Therefore, when studying large areas, it is common to use remote sensing, i.e., focusing on evaluating ET with data obtained from remote sensors [15][16][17], or hydrological balances using multiple sensors [18,19]. However, this approach (based on remote sensors) provides 2D information at surface level, and there are cases in which the 3D air movement is required in order to provide an in-depth study of the distribution of meteorological variables, i.e., urban forest interface regions [20] or semi-arid ones [21,22].
The aim of this research is to outline a methodology based on a 3D numerical CFD model, which in turn enables us to secure the values of the meteorological variables of air temperature, radiation, air humidity, wind characteristics as well as the humidity associated with the mass fraction of ET (named evapotranspiration rate, ETR) in large-scale outdoor 3D spaces (created from the zone topography), minimizing the computational resources needed. The model also allows for the introduction of the ET flow from the ground-the sum of evaporation from the land surface plus transpiration from plants-at the different points where it is discretized, known as the ET profile. The methodology and model proposed have been applied as a practical example in a zone corresponding to the semi-arid region of Cabo de Gata in Almería (Spain). The data obtained highlight the importance of topography when analyzing the characteristics of the meteorological variables and ETR distribution in the 3D domain.

Geometric Model
Firstly, a geometric model (with a variable extension depending on the computational resources available) from the selected study area referred to as a "study zone" is created. The structure of the study zone will be made up of the topography and kilometres of atmosphere. The model is separated into three different volumes by four different planes defined as Figure 1: Top of the atmosphere (a horizontal plane); land surface (created using data from a Digital Terrain Model and a GIS software tool); soil limit; and ET frontier, the latter two are acquired by copying the land surface.
Secondly, the complete model, which represents a virtual wind tunnel with three zones, is defined ( Figure 2a): The study zone in the centre and an interface zone both united, as well as an exterior domain. The circular shape of the interface zone allows both zones to rotate with respect to the exterior domain and create the angle necessary to simulate different wind directions without geometric modifications.
The wind velocity is included as a uniform profile on one side of the domain and develops as it passes through the exterior domain and the interface zone until being fully developed by the time it reaches the study zone.
An additional function of the interface zone is to adapt the elevation of the base of the plane of the exterior domain to fit the land surface of the study zone with the use of a slope (Figure 2b).

Theoretical Model
The numerical model used to study the fluid-dynamic behaviour of the domain employs CFD methodology based on Reynolds-Averaged Navier-Stokes equations (RANSs) for the three-dimensional flow complemented with the multispecies model to consider air, ET and air humidity-and a radiation model to include solar radiation.
Once the geometry is developed, the meshing of the volume of the completed model is carried out. Hexahedral cells with base dimensions between 100 to 200 m are used with a boundary layer mesh to adapt to the terrain.
Volume conditions must also be conveniently chosen. In the study zone, the air and water vapour mixture volume will have the condition of fluid, the ET volume will be considered as water vapour source (to simulate the existing vapour flow because of ET flow from the ground) and the soil volume will have the condition of solid. Additionally, the interface zone and the exterior domain will have the fluid condition ( Figure 3).
In general, multispecies models are used to study the mixing and transport of chemical species by solving, for each component, the conservation equation that describes its convection, diffusion and reaction sources. In this specific application, only the species mixing is considered without any chemical reaction. To predict the local mass fraction of a species i, the general form of the conservation equation is: where Y i is the local fraction, J i is the diffusion flux (for turbulent flows), ρ is the density, v is velocity vector, R i is the net rate of production by chemical reaction, S i is the rate of production by addition from the dispersed phase plus any user-defined source, Sc t is the turbulent Schmidt number, µ t is the turbulent viscosity and D i is the turbulent diffusivity (all units are in the International System). The Discrete Ordinates (DO) radiation model [23] is used to consider the radiation that reaches the land surface. This model was selected due to its ability to calculate through transparent materials and the possibility of defining various radiation bands. Additionally, the intensity of solar radiation, as well as the position and inclination of the sun each day, can be used as input.
The DO model solves the radiation transfer for a finite number of discrete solid angles, each associated with a vector direction, fixed in the global Cartesian system. The accuracy of the angular discretization is set following two guidelines: The number of angles and the number of pixels. This model transforms the general equation of heat transmission into a transport equation for radiation intensity in spatial coordinates.
where r is the position vector, s is the direction vector, s is the dispersion direction vector, s is the length of the path, a is the absorption coefficient, n is the refractive index, σ s is the dispersion coefficient, σ s is the Stefan-Boltzmann constant (5.672 × 10 −8 W/m 2 · K 4 ), I is the radiation intensity that depends on the position and the direction, T is the temperature, Φ is the phase function and Ω is the solid angle (all units are in the International System).

Simulation Methodology
The numerical calculations were completed using the commercial code Ansys Fluent 18.2. This code was used to solve the unsteady Reynolds-Averaged Navier-Stokes equations (URANS) by the finite volume method, converting them from differential equations into their numerical analogues (Eulerian method). The Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm was used to resolve the coupling between pressure and velocity fields. The discretizations of the spatial and temporal derivatives for the velocities and concentration in the equations were carried out by means of second-order schemes. The pressure staggering option (PRESTO!) discretization scheme was used for the pressure. Turbulence effects in the fluid flow were incorporated by means of the K-epsilon-RNG model [24], including buoyancy effects. Zheng [25] discovered that this turbulence model produced the best results for cases with features similar to the current test case.

Boundary Conditions
Different boundary conditions have been considered to simulate in the same model different processes that affect the evolution of meteorological variables.
A set of boundary conditions are considered for the exterior domain to allow for a means to include meteorological variability (Figure 4). One surface is set as a velocity inlet, aimed to introduce a wind velocity uniform profile and its humidity and temperature. The opposite surface is set as a pressure outlet. The other two lateral surfaces are set as symmetries (or wall without roughness) and the same for the ceiling. The floor is defined as a rough wall. An interface condition will allow the cylinder, where the study area is included, to rotate with respect to the domain, offering the chance to select a precise angle for wind inlet.
Furthermore, a set of boundary conditions for the study zone ( Figure 5) are applied to the four limiting surfaces. The upper one, called "top of the atmosphere", is defined as a wall without roughness, constant temperature and a semi-transparent material (such as glass). Moreover, in this wall, the radiation is included with an orthogonal beam using Φ and Ω values of 0.53 [26] to define the angle, so that direct and diffuse irradiation is included as input magnitude. The ET frontier is defined as an interior surface (virtual wall which is 100% permeable for flow). The land surface is a wall with friction, and a temperature dependant of the air temperature and solar radiation. The material of the land surface is opaque. The soil limit is considered also an opaque wall with frictionless and constant temperature conditions. The value of this temperature will be that of the ground at a depth of 10 m. The value used can be the measured or calculated by a thermal gradient law from the land surface temperature. Finally, lateral surfaces of the study area are all interiors (virtual wall which is 100% permeable for flow).

Study Area
The aim of this section is to provide an example of the possibilities that this methodology offers in order to study meteorological variables' behaviour in large 3D areas.
The study was carried out in a land area of Cabo de Gata, Almería, Spain, (Figure 6), in an area where databases of meteorological variables can be accessed and used as model inputs. It is a semi-arid region near the coast that includes an elevated zone with a hill of around 1000 m of elevation and a valley with altitudes close to sea level (0 m) (Figure 7).

Geometric and Numerical Model
The land surface of the study zone measures 23,600 × 24,800 sq. m and has been obtained using data from IGN (National Geographic Center) with a spatial resolution of 5 m. Before meshing the complete study zone, a pre-processing is carried out which simplifies it before transforming it into a .stl file that will be imported by the meshing software. The mesh consists of 2,497,360 cells with hexahedral structure (Figure 8). Tests were made with a larger number of cells (3.5 million) without appreciable variations in the results. To reduce the number of cells, the boundary layer extends from 2 m in height to 200 m. The exterior domain is carried out with tetrahedral cells. The size of the tetrahedrons diminishes while descending from high-order to low-order generations.  Before meshing the volumes, a boundary layer mesh was built in the study zone using structured meshing to obtain a better description of the boundary layer in the numerical calculations. The minimum volume cell is 3.17 × 10 4 m 3 . In addition, a quality analysis of the mesh yielded very satisfactory results, indicating a magnitude of the equitize skew below 0.4 for 100% of the cells in the mesh. This parameter shows the shape of the cells formed; values near 0 indicate more regular cells, which are far more likely to achieve satisfactory results.

Input Data
The meteorological data, wind speed (m/s) and direction, air humidity (%) and temperature ( • C) were drawn from the 30 min logs of a weather station from June 2006 to December 2014. The seasonal values used as model inputs are calculated on the basis of the mean of air humidity and temperature while the wind velocity and direction are obtained using compass roses per each season (Figure 9).  A representative value of radiation for each season was obtained from a study of a nearby zone (Tabernas desert) performed by [27]. The soil limit temperature, at 10 m depth, is estimated by using a typical value gradient at shallow levels of about 3 • C decreasing at every 100 m of depth [28] considering the ground surface at the same temperature as the air. Table 1 represents all meteorological data used as inputs. The average values for each season of ET profile (kg/m 3 s) used in the model as a vapour source input were obtained from the data of ET flow (mm/day) with a spatial resolution of 1100 m of GeoTIFF images provided by [21]. For this purpose, a Raster package implemented in the R software [29] was applied Figure 10.

Simulation Methodology
The numerical calculations were performed using Ansys Fluent 18.2 software. The SIM-PLE algorithm is chosen to resolve coupling between pressure and velocity fields. A secondorder scheme for spatial and temporal discretization was carried out, PRESTO! scheme for pressure and k-Epsilon-RNG for turbulence model.
To represent real conditions, it is necessary to perform two simulations. First, an initial simulation in a transitory state to place the study zone in the correct position to reproduce wind direction. Secondly, in a steady state, introducing all the inputs for each season. This part involves the greatest computational effort.

Results
Some results are shown as an example of those that can be obtained in the study zone giving the necessary information to evaluate ETR and the meteorological variables evolution in a 3D domain. As an example, four types of results are shown: (a) ET vertical flow; (b) maps of ETR, wind velocity (both closely to ground level) and ground temperature; (c) ETR evolution in zones of a vertical cross section; (d) vertical profiles of wind speed, ETR and temperature. Those profiles were obtained for points A, B, C and D. Point P represents the location of an ET flux tower ( Figure 11).
The convergence criterion to conclude the iterative computational process was set at values lower than 10 −5 for the normalised residuals of each governing equation. With these values the time required for the simulation of each season was about 36 working hours in parallel mode on a computer with an AMD-Ryzen9 processor of eight cores. The ET vertical flow (ET v f ) represents the ET transported by turbulent diffusion, which is several orders of magnitude lower that the ET transported in the horizontal layers by convection. It can be calculated using the results of the model at a certain height, such as ET v f = ETR · d · V z , where ETR is the evapotranspiration rate (p.u.), d is the density of the species mixture (kg/m 3 ) and V z is the vertical component of the wind velocity (m/s).
Comparing the values of ET v f (kg/m 2 s) registered in the ET flux tower at 3.5 m height with the mean values calculated using the model, it was found that they reliably represent ET movements, especially considering the accuracy with the low magnitude of ET v f (less than 0.1 kg/m 2 s). Figure 12 shows this comparison for the season of Spring as an example. Additionally, in all cases, it is noticeable that wind direction justifies the ETR concentrations more specifically in zones exposed to the incoming wind where the ETR is transported by convection out of the study zone (Figures 13 and 15), while the ETR is more concentrated in places sheltered from the wind (because of the hill effect) which for the same reason have ground temperatures slightly higher than the surrounding zones.
The model allows us to observe the evolution of the variables studied not only on horizontal planes but also on vertical ones. As an example, Figure 16 represents the variation of the magnitude of the ETR during the Spring season in two zones of a vertical cross section. Figures 17 and 18 represent profiles of wind velocities, air temperature and ETR during Summer and Winter of points A-D.
In the Summer (Figure 17) points A-D have similar parabolic profiles for the three variables, as they are facing the wind direction with no obstacles. The ETR and ground temperatures are slightly higher at point C (at the different heights), indicating that it is in a zone a little more protected from the incoming wind.
In the case of the Winter (Figure 18) all profiles (wind, ETR and temperature) of points A and B are distorted parabolic as they are in a zone highly protected from the incident wind, contrary to points C and D with parabolic profiles as there is no obstacle that affects the wind profile. Moreover, in this case, due to the previously mentioned protection effect, the wind velocity values of A and B are lower than those of C and D, which leads to higher ground temperatures and ETR at different heights.

Conclusions
The methodology and 3D CFD model put forward provide the means to simulate meteorological variable evolution as well as ETR variations in outdoor spaces. Its design facilitates a simple introduction of the meteorological variables air temperature and humidity, wind velocity and angle, solar radiation and ET profile, as model inputs. Furthermore, it ensures the adequate development of the wind profile and that the change in wind velocity direction can be easily reproduced. All the above together with the meshing employed create the conditions required to apply the model to large land extensions and thereby obtain satisfactory results with reduced computational resources.
The model was applied to a land area of Cabo de Gata (Almería, Spain) of 23,600 × 24,800 m, as this is an area with similar meteorological conditions on different days of each season of the year, and therefore, using average data per season, representative results can be obtained.
Some relevant results were obtained that demonstrate the usefulness of the model, such as the influence of wind direction and topography in the ET transported by convection and land temperature or the unerring accuracy of the model concerning calculating ET vertical fluxes.