Contribution to the Sustainability of Agricultural Production in Greenhouses Built on Slope Soils: A Numerical Study of the Microclimatic Behavior of a Typical Colombian Structure

: The use of covered structures is an alternative increasingly used by farmers to increase crop yields per unit area compared to open ﬁeld production. In Latin American countries such as Colombia, productive areas are located in with predominantly hillside soil conditions. In the last two decades, farmers have introduced cover structures adapted to these soil conditions, structures for which the behavior of factors that directly affect plant growth and development, such as microclimate, are still unknown. Therefore, in this research work, a CFD-3D model successfully validated with experimental data of temperature and air velocity was implemented. The numerical model was used to determine the behavior of air ﬂow patterns and temperature distribution inside a Colombian passive greenhouse during daytime hours. The results showed that the slope of the terrain affects the behavior of the air ﬂow patterns, generating thermal gradients inside the greenhouse with values between 1.26 and 16.93 ◦ C for the hours evaluated. It was also found that the highest indoor temperature values at the same time were located in the highest region of the terrain. Based on the results of this study, future researches on how to optimize the microclimatic conditions of this type of sustainable productive system can be carried out.


Introduction
Cultivation systems in greenhouses are becoming one of the most widely used alternatives for the sustainable intensification of food production such as fruit and vegetables in many countries [1,2]. These agricultural production systems allow the increase of crop yields per unit area and optimization of water management per unit of product, compared to open field crops [3]. In addition, the use of these structures allows food production in regions where climatic conditions are adverse for agricultural production, ensuring the supply of fresh food and increasing the levels of food security for these territories.
In tropical countries located in the Central American and Caribbean region, the most commonly used greenhouses are passive greenhouses, which depend entirely on natural ventilation for microclimate management [4,5]. Therefore, they are low-cost and low-tech structures that do not have the equipment for heating, cooling and humidification and are therefore tentative structures for low-income producers [6]. Likewise, since these structures depend on natural ventilation, they are considered to be environmentally friendly, as the environmental load from energy consumption due to the use of fossil fuels is almost non-existent [7][8][9].
In passive greenhouses, natural ventilation is the phenomenon in charge of managing the microclimate conditions in the interior environment of each structure. Therefore, by means of air flows generated by thermal or wind pressure or a combination of these two driving forces, it is possible to regulate the thermal and humidity excesses generated in the interior environment of a greenhouse [10,11]. These same air flows allow the exchange of air with the outside environment, becoming the only means of carbonic enrichment of the air inside the greenhouse [12].
The study of natural ventilation in greenhouse structures or any other structure for agricultural use is not a simple activity to perform experimentally. Although the development of climate monitoring equipment allows the study of air flows through sonic anemometry, these sensors only allow the determination of speed and direction of air flow at a spatial point for a given time [6,7,13]. Therefore, to characterize the airflow patterns, many anemometers would be needed, which in logistical and economic terms is quite laborious.
Other study alternatives include indirect methods such as tracer gas, experiments with colored smoke streams, 3D imaging or wind tunnel tests with scale models and the results of some of these investigations can be reviewed in the work developed by Akrami et al. [14]. However, there is another methodology, based on modeling and simulation that allows studying the air flow patterns inside a greenhouse and its effects on the spatial temperature distributions. This methodology is computational fluid dynamics (CFD), which is considered an agile and robust simulation technique [15][16][17].
The CFD study has allowed the study of almost all existing passive greenhouse models worldwide, and has also been implemented as a tool for the optimization of climate management in greenhouses with active climate control [18][19][20]. There are currently a considerable number of studies applied to agricultural investors that have used CFD simulation to analyze thermal distribution [10,21,22], moisture distribution [5,21,23], efficiency of heating systems [24][25][26], the operation of cooling systems [27][28][29] and the microclimate generated as a function of different ventilation configurations [7,12,[30][31][32]. In the specific case of Colombia, the use of CFD and other modelling techniques to determine the behavior of the main greenhouse structures used has been carried out for the specific conditions of the Bogotá savannah and for a topography of flat soils, which is where the ornamental and cut flower industry is generally developed [33][34][35].
In Colombia, horticultural and fruit production, especially that of small producers, is located in the Andean region. In this natural region there are approximately 14 million hectares with soils suitable for agriculture, although approximately 75% of these areas are located in hillside regions [36]. In these regions, for about a decade, different roof structures have been established, mainly greenhouses and mesh houses, where their microclimatic behavior is still unknown [2]. The microclimate variables inside the greenhouse, such as temperature and relative humidity, are directly related to processes such as transpiration and photosynthesis; processes that affect the growth and development of the plants and consequently the final production obtained [37]. Therefore, knowledge of microclimatic behavior allows the establishment of cultural management practices such as irrigation and fertilization [2].
There are limited studies analyzing the greenhouse conditions in hillsides such as that by Rojas Rishor [38], who found temperature and relative humidity gradients inside a tomato greenhouse built in a hillside area of Costa Rica. For this same greenhouse in the study developed by Flores-Velasquez et al. [39], temperature and relative humidity gradients of 3.14 • C and 11.25% in the slope direction, and 0.63 • C and 6.04% in the cross section, respectively, were reported. Finally, in a study developed by Kuroyanagi et al. [40]. For a greenhouse located in Japan, thermal gradients of up to 2.5 • C were reported inside the structure, gradients that were located at the top of the slope of the land where the greenhouse was built.
Accordingly, the objective of this research was to determine the behavior of air flows and thermal distribution for a Colombian tunnel greenhouse located on a hillside terrain. For this purpose, a 3D numerical simulation model was implemented, adopting the dominant hourly climatic conditions of the experimental site.

Description of the Greenhouse and Experimental Site
The greenhouse evaluated in this study is of the multitunnel type with a plastic cover, composed of five adjoining spans, with the width of each span being 6.8 m. Therefore, the total width of the greenhouse is 34 m (X-axis). The longitudinal dimension of the greenhouse is 30 m (Z-axis), the greenhouse has a minimum height above the gutter of 3.5 m (H min ) and a maximum height (H max ) above the ridge of the zenithal ventilation of 5.8 m (Figure 1). The greenhouse has lateral ventilation areas on both sides, with a maximum opening of 2.3 m, so that the ventilation area in the lateral region is 138 m 2 (13.5% of the covered surface area). The ventilation surface is complemented by a ventilation area on each of the façades of 2.2 m total opening, giving a total ventilation surface on the façades of 149.6 m 2 (14.6% of the covered surface area), and finally, each span has a fixed roof ventilation of 0.66 m opening, with the total roof ventilation area being 99 m 2 (9.7% of the covered surface area) distributed over five spans. The greenhouse was built on a site that had a transversal slope (Sx) of 7.35% and a longitudinal slope (Sz) of 50% ( Figure 1).
The experimental site where the greenhouse was located corresponds to the municipality of Filandia (75 • 42 57.39 W, 4 • 41 11.7 N and altitude: 1586 m), in the department of Quindío. This region has a climatic behavior characterized by humid tropical climate conditions, where the average multiannual temperature value for a 35-year period is 20.95 • C and the minimum and maximum mean values are 15.4 • C and 26.6 • C respectively, while the mean annual value of accumulated precipitation is 2060 mm.

CFD Numerical Model
The numerical solution is developed from the spatial discretization of the nonlinear partial Navier-Stokes equations. The transport phenomena by free convection can be described by Equation (1), which is considered the general transport equation for a fluid in steady state in a three-dimensional field of action.

∂(u∅) ∂x
where y, x and z represent the coordinates in Cartesian space, u, v and w are the components of the velocity vector, ∇ 2 is the Laplacian operator, Γ is the diffusion coefficient, ∅ represents the concentration of the transported quantity in dimensionless form (momentum, mass and energy) and S ∅ is the source term. The turbulent nature of the air flow was simulated by using the empirical turbulence model k-ε standard, a model that has been extensively validated in studies of natural ventilation of agricultural buildings, showing a good fit of the real airflow behavior with respect to the simulated process [7]. The buoyancy effects of air caused by gravity and air density changes were added to the momentum equation as a source term by the Boussinesq approximation, and this model is represented by Equation (2) [41].
where ρ and ρ 0 are the density at a given time and the reference density respectively; β is the coefficient of thermal expansion of air; T and T 0 are the temperature and the reference temperature. The radiation model selected was the discrete ordinate (DO) model with angular discretization. This model allows the calculation of radiation and convective exchanges occurring in the computational domain, treating the greenhouse canopy as a semitransparent medium [42]. Equation (3) describes the DO radiation model.
where I λ is the intensity of radiation at one wavelength; ⇒ σ s and a λ are the spectral scattering and absorption coefficients; n is the refractive index; ∇ is the divergence operator and Φ, T and Ω are the phase function, local temperature ( • C) and solid angle, respectively. The non-grey model was activated by dividing the radiative spectrum into two specific wavelengths, the first being the solar radiation length (0.4-2.4 µm) and the second the long wavelength (2.4-100 µm). Likewise, in order to simplify the resolution of the 3D CFD model, no crops were included, in order to speed up the numerical calculation and establish the behavior of air flow and temperature under the worst possible scenario, which are conditions where there are no plants present, a scenario where a large part of the incident radiation on the interior of the greenhouse is converted into heat, which generates an increase in the temperature of the interior air [2].

Discretisation of the Computational Domain
The computational domain that included the greenhouse was configured in terms of its dimensions following the recommendations reported in the work developed by Franke et al. [43] and implemented by Kim et al. [44]. Where minimum lengths are suggested from the edges of the domain to the evaluated structure, these dimensions are established according to the maximum height of the greenhouse (H max ), as shown in Figure 2. Therefore, the overall dimensions of the computational domain for the x, y and z axes were 208 m, 58 m and 204 m respectively. This computational domain was divided into an unstructured numerical grid with a size determined from a test of independence on the size of the numerical grid. The results obtained in this test will be discussed later in this paper. The quality of the selected numerical grid was determined from the orthogonality parameter, finding that 95.2% of the grid elements had a value higher than 0.93, considered to be of high quality [45]. It should be noted that both the size of the numerical grid and the quality of its elements should be verified by the two methodologies discussed above, since these factors are highly determinant of the quality and accuracy of the numerical results obtained [11].
As for the solver settings, a semi-implicit method was used for the coupling of the pressure-velocity equations using the SIMPLE algorithm. For the momentum, energy, and pressure terms, a second order discretization scheme was used, which allows more accurate results to be obtained compared to the first order discretization schemes [46]. Finally, convergence criteria were established in 10 −6 for the energy equation and in 10 −3 for the other variables of the model [1]. For the radiation model, we defined the following discretization conditions; 3 theta and phi divisions, 3 theta and phi pixels and 10 energy iterations per radiation iteration.

Boundary Conditions
In the computational domain, the physical and optical properties of the materials described in Table 1 were used, taking the reference values used in the work carried out by Senhaji et al. [47]. On the other hand, in these 13 simulations a power absorption coefficient (αλ) was established to obtain a transmittance of zero for long wavelengths and 0.8 for short wavelengths. This can be determined by using Equation (4), which relates the absorptivity (α) and the thickness of the medium (e). The boundary conditions established were symmetrical for the computational domain boundaries parallel to the air flow. In the case of the greenhouse floor and computational domain, it was simulated as a mixed boundary condition for radiative and convective heat transfer and as an opaque medium for diffuse radiation. On the other hand, the greenhouse roof was simulated as a semitransparent material with a finite thickness of 0.2 mm and a two-sided wall boundary condition and a coupled thermal condition. The greenhouse walls were simulated as an opaque and coupled medium.
In the case of the computational domain roof region, a wall boundary condition was established with a solar radiation flux. A logarithmic input velocity profile was considered at the input of the computational domain. The profile was linked to the main CFD module using the user-defined function (UDF), as explained in the work developed by Villagran et al. [7], using the Equation (5): where V 2 is the calculated wind velocity in h 2 height; V 1 is the reported wind velocity in h 1 height and z 0 is the roughness coefficient length, the surface roughness is a nominal value that depends on the type of coverage and roughness of the soil material, the type of landscape and the spacing between obstacles, for this case was considered a value of 0.03 m in relation to class 1 of roughness according to Bañuelos-Ruedas et al. [48]. In the case of the air flow outlet limit, a pressure outlet limit was established, and, finally, the ventilation areas of the greenhouse were set with an interior condition limit.
Likewise, for the recording of the temperature inside the structure over the longitudinal middle region (Z-axis = 15 m) and over a cross section of the structure and a height of 1.5 m above ground level were installed 7 thermocouples (range:−(40)-70 • C, resolution:  With the climatic information collected during the experimental period, the average values of temperature, solar radiation, wind speed and direction were calculated on an hourly scale. These values were used as initial conditions to perform 13 steady-state simulations for each of the hours of the day between 6:00 and 18:00 h, as shown in Table 2. The value of the wind speed outside corresponds to the value of the reference speed (V 1 ) in Equation (5).
To evaluate the validity of the numerical model, once the hourly simulations had been carried out, the temperature and wind speed data were extracted for each sampling point inside the greenhouse. In a first analysis, the validity or rejection of the numerical simulation model was checked through a statistical analysis by means of a hypothesis test for the difference of measurements with homogeneous variances. In a subsequent step, the simulated data were compared with the measured data through goodness-of-fit parameters such as mean absolute error (MAE) Equation (6) and root mean square error (RMSE) Equation (7).
where Dm j and Ds j are the measured and simulated wind speed and temperature values at point j respectively and M is the number of samples. Once the CFD model was validated, we proceeded to the post-processing phase where air flow behavior curves were extracted from the numerical solution. In this case, inside the greenhouse in a plane (z, x) at a height (y) of 1.5 m above the ground level, a total set of 18,951 spatial data were extracted for temperature and wind speed; data that will be analyzed in the results section.
In this post-processing phase, the data for the calculation of the ventilation rate are also extracted. This process is performed by the conventional method that consists of verifying the greenhouse outflows at each ventilation opening. Therefore, the greenhouse air exchange rate (GAER) per unit time in minutes is calculated by Equation (8).
where V i is the air velocity of the outflow (ms −1 ), A vs is the area of the ventilation with exhaust air flow (m 2 ) and V G is the volume of the greenhouse (m 3 ).

Grid Independence Test
For the grid independence test, the mean values calculated for hour 10 were selected as initial conditions (Table 2). With these values, eight simulations were performed, where the variable of change was the number of elements of the numerical grid, for which a total of eight grids of different sizes were available. These grids varied between a minimum of 846,819 elements for grid 1 and a maximum of 20,735,202 elements for grid 7.
For this case, grid number 4 was selected, which is made up of a total number of 6,717,534 elements. This selection was made because quantitatively it was possible to verify that from this grid size the solutions obtained for temperature and velocity are independent of the number of elements that make up the numerical grid, as can be seen graphically in Figure 4. This independence test should always be performed in CFD studies since its results can be used to define the numerical grid size that guarantees the independence and accuracy of the solutions at the lowest possible computational cost [49,50].

Data and Validation of CFD Model
The air velocity values obtained under experimental and simulated conditions at each of the three measurement points, for each of the 13 scenarios evaluated (6:00-18:00 h) can be found in Table 3. The same values obtained for the seven sampling points were obtained for temperature and are summarized in Table 4. These values were the data used to perform the validation process of the CFD model. Table 5 shows the results obtained for both air velocity and temperature from the homogeneity of variances tests and the respective mean comparison test between measured and simulated data. In general, it was found that there were no significant differences between the variances of the simulated and measured data for temperature and air velocity at a 95% interval, and likewise, no statistically significant differences were observed between these sets of measured and simulated data at a 95% confidence interval for the difference in means. In all these statistical tests, the null hypothesis (H 0 ) is accepted.      Figure 5 shows the behavior of the measured and simulated data for the temperature variable, over the experimentally studied cross section. In general terms, it can be observed that the data sets for the spatial points sampled show very similar behavior in terms of magnitude and trend in the totality of the hours evaluated. This analysis was complemented with the numerical analysis based on the goodness-of-fit criteria selected, the results of which are shown in Table 6.
For the temperature variable, it was found that the MAE value oscillated between minimum and maximum values of 0.124 • C and 0.346 • C for 6:00 and 13:00 h respectively, while the RMSE values oscillated between 0.140 and 0.397 • C for these same hours. These values found are within the ranges reported in previous studies where CFD models were implemented to determine airflow patterns and thermal distribution in naturally ventilated greenhouses [34,51].  In the case of the air velocity variable inside the greenhouse evaluated, the spatial distribution of the measured and simulated values for each hour was plotted, finding that the experimentally sampled values showed high coincidence with the simulated values ( Figure 6). Likewise, the results of the quantitative analysis of the measured and simulated data sets through MAE and RMSE (Table 6) showed values lower than 0.107 ms −1 and 0.128 ms −1 respectively, values that can be considered highly acceptable for these types of modelling and simulation study [7,12]. Therefore, under this quantitative and qualitative analysis of the temperature and wind velocity data, it can be concluded that the numerical model has a high capacity for the prediction of the air flow patterns and the spatial distribution of temperature in the greenhouse studied.

Qualitative and Quantitative Characteristics of Airflow Patterns
The airflow patterns inside the greenhouse for some of the simulated scenarios are shown in Figure 7. For the 07:00 and 09:00 h where the wind comes from the south side, air flow patterns were observed entering the greenhouse through the side vents and the ventilation area of the façade exposed to windward, with these air flows moving in the longitudinal direction of the greenhouse ascending the slope of the land, until they leave the interior of the structure through the ventilation area located on the north façade of the greenhouse. It is also observed that in spans 3, 4 and 5 the airflow patterns enter with a low velocity of approximately 0.23 ms −1 , and then accelerate to a velocity of 0.67 ms −1 as these airflows ascend the slope of the terrain. This acceleration of the airflow may be influenced by the convective effect caused by changes in density as a function of the change in temperature of the air inside the structure [47].
It is also observed that the fixed vents in the roof region for this case function as air inlet areas over a region approximately in the lower 3/4 of the longitudinal section and an air outlet area approximately over the upper 1/4 of the longitudinal section in all the spans. This may indicate that the longitudinal slope generates on the fixed windows of the greenhouse roof regions positive pressure for the air inlet zones and negative pressure for the window section where there is an air outlet flow [52,53].
For the 11:00 and 13:00 h, a change in the air flow pattern inside the greenhouse is observed, mainly influenced by the direction of the dominant outside wind, which in this case comes from the west side of the greenhouse. Therefore, the airflows enter the greenhouse through the windward side ventilation with average velocities of 1.07 ms −1 and 1.19 ms −1 for each case. Between spans 1 and 2 these airflows presented a deceleration and again from span 3 to span 5 they accelerated until they reached and exceeded the entry velocity values (Figure 7).
The air flows inside the greenhouse moved in two identifiable patterns (Figure 7, hour 11:00 and 13:00). The first moved close to the lower façade of the greenhouse and as it crossed the cross section it moved in the opposite direction to the longitudinal slope until it exited completely through the leeward side vent. The other flow moved through the greenhouse in a counter-slope direction and exited the structure through the façade located on the north side (Figure 7, hour 11:00 and 13:00). In these scenarios a flow pattern was again identified in the fixed roof windows similar to the one already discussed in the two previous scenarios.
In the case of 15:00 h, where the wind comes from the south-southwest side, the airflow entered through the lateral side of the west side and through the windward façade. This air flow then adopted a very similar behavior to that at 13:00 h with the only difference being that in this case a higher velocity (1.29 ms −1 ) was observed in the entry of the air flow into the structure and in turn a higher velocity in the fixed roof ventilation areas (between 1.18 and 1.39 ms −1 ).
Finally, for the 17:00 h scenario, the same characteristic pattern already discussed for the 07:00 and 09:00 h scenarios are observed. Within these results it can be observed that the variables have already been discussed in previous works where it was defined that the air flow patterns are influenced by the shape, height and size of the greenhouse, the geometry of the roof, the arrangement of windows, and the speed and direction of the outside air flow [6,11,54]. The topography of the terrain where the greenhouse is built must be added, since this factor also influences the behavior of the air flows, both in their direction of displacement and in their velocity value.
The presence of longitudinal and transverse slopes on the ground influences the air flow patterns since these slopes favor a behavior where air flows lose velocity once they enter the greenhouse and in turn the air flow pattern tends to be directed vertically towards the roof ventilation areas while moving in a counter-slope direction. This type of flow pattern had already been reported in the study developed by Taloub et al. [55], who concluded that any type of slope on the floor of the greenhouse, no matter how small, destabilizes the horizontal airflow and promotes the generation of vertical airflows and recirculatory movement patterns with low velocity zones near the leeward vents.

Air Velocity Inside the Greenhouse and Calculated Ventilation Rates
The average air velocity values and their standard deviation for each of the evaluated scenarios are shown in Table 4. In general terms, the average air flow velocity in the greenhouse ranged from a minimum of 0.312 ± 0.08 ms −1 at 06:00 h to a maximum of 1.336 ± 0.29 ms −1 at 15:00 h. The average velocities found for this greenhouse were higher than those reported for other greenhouses studied under the conjugated ventilation configuration with side and roof ventilation areas, where the value reported in several studies ranged between 0.5 and 0.7 ms −1 [56][57][58]. This behavior may be highly influenced by the slopes of the land where the greenhouse was built, and also by the lack of insect netting in the ventilation areas, since in Colombia the use of insect netting is not a common practice.
On the other hand, the behavior of the normalized velocity, which represents the relationship between the average inside air velocity and the outside wind speed, showed a minimum and maximum of 83.1% and 146.8%, respectively (Table 7). These values of normalized velocity, contrary to other studies of natural ventilation, were relatively high, which is influenced by the acceleration of the air flow inside the greenhouse. It was also found that for the period between 15:00 and 18:00 h, the value of the normalized velocity was higher than 100%, which means that the air flow velocity inside the structure was higher than the external wind flow. The values of the highest normalized velocity agreed with the time when the temperature of the outside environment started to decrease, therefore it can be mentioned that this acceleration of the air flow inside the greenhouse may be strongly associated with the thermal effect of natural ventilation by buoyancy; an effect that generally produces rapid changes in temperature and velocity inside a greenhouse [59]. This buoyancy phenomenon is caused by convective movements generated between the soil and the greenhouse cover, since at that time the soil is the surface of higher temperature due to energy storage throughout the day, while the cover usually cools rapidly to the level of the outside ambient temperature [34,60].
This type of behavior analyzed serves as a basis for generating a recommendation for management and microclimatic optimization for night-time hours in these structures that do not have opening and closing systems for side and roof vents. A useful alternative to prevent the accelerated loss of stored heat during the hours of high temperature and solar radiation would be the use of a mechanism to close the ventilation areas just around 15:00 h, which is the time when the ambient temperature begins to decrease. Therefore, a mechanism to close each of the side and roof vents would increase the hermeticity of the greenhouse and prevent thermal inversion phenomena during nighttime hours [50,61,62].
Greenhouse air exchange rates (GAERs) in volumes per minute (Vmin −1 ) were calculated for each of the simulated hours, with these values ranging between a minimum and maximum value of 0.28 Vmin −1 and 1.29 Vmin −1 (Figure 8). These data show that only for a period of 3 h (12:00-15:00) was the GAER value higher than the minimum recommended value for passive greenhouses, which should be at least 1 volume renewed per minute (Vmin −1 ) [59]. Therefore, it can be deduced that the greenhouse evaluated presents deficient ventilation rates for 69.7% of the hours of the day-time period, and so it is recommended that future studies address design alternatives in the structure to maximize this index, since its influence on the microclimate behavior of the structure is very relevant [63]. In general terms, it can be observed that as the day went by 6:00-14:00 h and at the same time there was an increase in the outdoor wind speed, there was an increase in the value of the air exchange rate of the greenhouse (Figure 8). Therefore, for the period between 6:00 and 14:00 h, the GAER value increased from 0.28 to 1.23 Vmin −1 , which represents an increase of 331% between 6:00 and 14:00 respectively. For this same period, the wind speed had an increase of 303%, increasing from a value of 0.32 ms −1 at 6:00 to a value of 1.29 ms −1 at 14:00 ( Figure 8).
The above reaffirms what has already been reported in other studies of natural ventilation, where it was identified that there is a direct relationship between the ventilation rate and the outside air flow velocity [8,64]. Figure 9 shows the linear relationship between GAER (y) and the external wind speed (x), which can be represented mathematically with the equation Y = 1.0291 (X) − 0.1099, which in this case is a relationship with an acceptable coefficient of determination R 2 of 0.8085.

Qualitative and Quantitative Characteristics of the Spatial Distribution of Temperature
The spatial distribution of temperature for each of the simulated scenarios can be seen in Figure 10. The behavior found shows the thermal gain that occurs in the structure via the greenhouse effect between 6:00 and 14:00 h, where the volume of air enclosed inside the structure is heated by the combined effect of higher values of solar radiation and the increase in temperature in the outside environment, which is a typical behavior of structures at the tropical level [65]. Subsequently, from 15:00 to 18:00 h, a totally opposite phenomenon is observed, where an accelerated cooling of the air volume occurs due to the reduction of the temperature and radiation levels in the external environment. This heat loss is more accelerated in this structure since it does not have enclosure systems that limit the air flow through the lateral and roof ventilation areas.
Another characteristic observed is the spatial heterogeneity in the value of the temperature presented in this type of greenhouse. For the first hours of the day some zones of higher temperature are observed just below the areas covered by spans 2, 3 and 4, and this same behavior is observed for the hour 17:00 (Figure 10), while for 11:00, 13:00 and 15:00 a zone of higher temperature is observed in the region below spans 4 and 5.
Another aspect to be highlighted is the influence of the wind direction on the location of the high temperature zones. For the cases where the wind direction comes from the south, the region of higher temperature is concentrated from the middle longitudinal zone of the greenhouse towards the leeward façade, hours 7:00, 9:00 and 17:00 ( Figure 10). For the hours 11:00, 13:00 and 15:00, however, where the outside wind comes from the west, the warm regions are located on the leeward side and on the area near the north façade just below spans 3, 4 and 5 ( Figure 10).
Qualitatively, it can be observed on the spatial distribution patterns of temperature for each evaluated scenario that tends to generate a zone with higher temperature values just above the regions of higher altitude in the greenhouse, these are areas that in our case are located near the north facade and the top of the east side of the structure. These regions of higher temperature located in the upper slopes of the terrain were reported in a Costa Rica greenhouse evaluated in the work developed by Rojas Rishor [38], who also reported these same heat patches even with the outside airflow moving downslope. Therefore, it is important to recommend that developers building on slopes should monitor and establish cooling strategies in the higher altitude regions of the terrain.

Temporal Behavior of the Temperature Inside the Greenhouse
Temperature is one of the main microclimate variables directly involved in the growth and development of most crops of commercial and food interest, and this variable also affects the photosynthetic activity of the plant and the quality of the final product. For this reason it is one of the most studied microclimate variables inside greenhouses [63,66].
The qualitative analysis of the spatial temperature patterns was performed by calculating the average values of the temperature inside the greenhouse (T inside ) and the thermal gradient generated inside the structure (∆T inside ). This gradient represents the difference between the maximum temperature (T max ) and the minimum temperature (T min ) in the indoor environment for each of the hours (Table 8). The values of Ti nside ranged between a minimum value of 16.98 ± 0.28 • C for the 6:00 h and a maximum value of 27.92 ± 1.75 • C for the 14:00 h. Between these times it was observed that the greenhouse temperature was gradually increasing, and this type of behavior is very characteristic of greenhouses located in the tropics. Once the highest T inside value was achieved at about 14:00 h, this value decreased again to a value of 19.70 ± 0.19 • C. Therefore, the greenhouse presented a thermal gain of 10.94 • C in a period of 9 h and an energy loss of 8.22 • C in only 4 h, demonstrating numerically what has already been analyzed qualitatively in the previous section in reference to the accelerated energy loss of this greenhouse.
The values of T inside found for this structure are within the ranges reported as adequate for most of the species of food/fruit interest such as tomato, cucumber and pepper, additionally, at no time of the day does the temperature exceed the maximum recommended for these types of crop, which is around 35 • C [5,67]. Another factor related to the behavior of temperature is its spatial distribution. In recent years, the interest of researchers has focused on determining and studying the heterogeneous distribution of temperature inside greenhouses, since it is a factor that affects the quantity and quality of the final products harvested [68][69][70].
Regarding the evaluation of the spatial distribution of the temperature, it can be observed that there are behaviors that can be considered as homogeneous for the hours 6:00 and 17:00 where the ∆T inside value showed values lower than 2 • C [70]. For the remaining hours, a high thermal heterogeneity is observed between 9:00 and 16:00 h, where the ∆T inside value reached a magnitude of up to 16.93 • C at 13:00 h. This value may be considered too high and inadequate, even more so if we take into account that it is a greenhouse composed of only five spans and that additionally it does not have porous insect-proof mesh protection in the ventilation areas.
This behavior is not ideal for this type of structure since the plants are subjected to different internal temperature conditions for the same external climatic conditions, which generates non-uniform growth and production of crops. The above has a particular aggravating factor, and that is the fact that for these greenhouses, cultural tasks such as irrigation and fertilization were calculated and applied in the same way and quantity for all plants.
On the other hand, the thermal differential between the outdoor and indoor environments (TD = T outside -T inside ) was calculated, where mean values ranging from 0.28 • C for 6:00 h to a maximum value of 2.61 • C for 13:00 h were obtained ( Figure 11). The figure shows that both the TD thermal gradient and the climatic heterogeneity are dependent on the level of solar radiation; therefore, the higher the radiation levels, the greater the thermal gradients and spatial variability of the temperature inside the greenhouse. These results are consistent with those reported in previous studies of the microclimatic behavior of Colombian greenhouses [65,70,71]. Figure 11. Temporal behavior of the thermal differential between the outside and inside temperatures of the greenhouse and hourly behavior of solar radiation.
Finally, it is important to mention that the results obtained in this research are a relevant technical input that will allow future microclimatic optimization strategies to manage the temperature conditions inside roof structures built on slope soil. Thermal heterogeneity inside this type of structure is an undesirable condition that puts at risk the sustainability of agricultural production in the future year, considering the increase of temperatures due to climate change as a serious threat in the tropical region [1,2]. Therefore, it is necessary that future studies consider the search for alternatives that allow the management of radiation and temperature in such structures. Currently there are several studies where the use of passive methods has been used to optimize or increase the degree of management of the microclimate in passive greenhouses [1,2]. These passive methods include the management of the radiation level inside the structure through the use of shading nets. The use of shading nets reduces the temperature inside the greenhouse, improves the distribution of climatic parameters such as temperature and relative humidity, and additionally helps to reduce the water consumption of the plants and improve the quality of the final products [72].
It is also possible to contemplate the use of solar panels on the region of the roof where the greatest thermal heterogeneity is generated. This would reduce the temperature value and additionally obtain an energy benefit through the collection of solar energy that can be used for other daily tasks on the farm. This has also showed in several studies that this alternative is economically and technically feasible, although it must always be accompanied by studies that allow establishing the optimal level of shade that does not negatively affect the yield and quality of the crops [73,74].
Another alternative that can be analyzed is the use of bedrock storage systems and water bags as proposed by Gourdo et al. [75] and Bazgaou et al. [76], who reported temperature reductions in a plastic-covered greenhouse between 3 and 5 • C in the daytime hours of the day period. Finally, the CFD model validated in this research can also be used to simulate unbuilt scenarios that modify the current structure and ventilation configuration of this greenhouse, trying to improve the movement of air flow patterns and the spatial distribution of temperature following the methodological approach developed in several studies of natural ventilation in greenhouses that seek to optimize the generated microclimate [5,7,34,77,78].

Conclusions
The experimental validity of the 3D CFD model analyzed by means of trend curves and numerical parameters of goodness-of-fit between measured and simulated data exhibited high agreement. Therefore, the numerical model implemented in this research showed a high level of prediction of the air flow patterns and the spatial distribution of temperature generated inside a Colombian greenhouse established on slope soil.
This study identified that the greenhouse evaluated under the current ventilation configuration does not comply with the minimum ventilation rates recommended for a naturally ventilated greenhouse, for 69.7% of the hours evaluated. This generates a spatial distribution of temperature that reaches differences of up to 16.93 • C between the highest and lowest temperatures. These conditions are considered highly heterogeneous and will certainly affect the growth of the crops in this type of greenhouse.
Although the behavior of natural ventilation and the spatial distribution obtained for this greenhouse showed direct relationships with wind speed and the level of solar radiation of the external environment, it was observed that the slope of the terrain affects the displacement and speed of air flows, which generates high temperature zones on the higher regions of the terrain.
Finally, this study can be the baseline for future research aimed at evaluating alternatives to optimize the microclimatic performance of this type of greenhouse built on sloping terrain, since there is currently limited available information for greenhouses under these construction conditions. The validated CFD model can also be implemented for the study of the microclimatic behavior of other types of greenhouses built in Latin American or Caribbean countries, where protected crops are increasingly implemented in soils with slope conditions.