Inﬂuence of the Height in a Colombian Multi-Tunnel Greenhouse on Natural Ventilation and Thermal Behavior: Modeling Approach

: The dimensions of a passive greenhouse are one of the decisions made by producers or builders based on characteristics of the available land and the economic cost of building the structure per unit of covered area. In few cases, the design criteria are reviewed and the dimensions are established based on the type of crop and local climate conditions. One of the dimensions that is generally exposed to greater manipulation is the height above the gutter and the general height of the structure, since a greenhouse with a lower height has a lower economic cost. This has led some countries in the tropical region to build greenhouses that, due to their architectural characteristics, have inadequate microclimatic conditions for agricultural production. The objective of this study was to analyze the effect on air ﬂows and thermal distribution generated by the increase of the height over gutter of a Colombian multi-tunnel greenhouse using a successfully two-dimensional computational ﬂuid dynamics (CFD) model. The simulated numerical results showed that increasing the height of the greenhouse allows obtaining temperature reductions from 0.1 to 11.7 ◦ C depending on the ventilation conﬁguration used and the external wind speed. Likewise, it was identiﬁed that the combined side and roof ventilation conﬁguration (RS) allows obtaining higher renovation indexes (RI) in values between 144 and 449% with respect to the side ventilation (S) and roof ventilation (R) conﬁgurations. Finally, the numerical results were successfully ﬁtted within the surface regression models responses.


Introduction
In Colombia, plasticulture and protected agriculture have promoted a method of agricultural production that generates higher crop yields per unit of available land area, compared to open field production [1]. The use of passive greenhouses, where the microclimate is managed by means of natural ventilation, has become more widespread. This type of greenhouse is also quite common in other regions of the world [2]. This cultivation method allows, among others, the intensification of horticultural or ornamental production [3], better management of water resources and fertilizer applications [4,5], partial or total control of the micro-climatic variables that affect crop growth and development [6]. In recent years, this production technology has also become a crop protection tool used by producers in the face of increasingly frequent and severe weather events due to climate change [3,7,8].
The use of greenhouses or protected agriculture structures with plastic roofing worldwide has been estimated at more than 3.5 million hectares, which have been established tially, which in turn generates higher temperature values inside the greenhouse [31][32][33]. For the height dimension, Boulard and Fatnassi [34] using a CFD-2D study on an 8-span arch-type greenhouse found that when the height above the greenhouse gutter increases from 3 to 5 m, the reduction of the thermal gradient between the interior and the exterior is reduced from 3.5 to 2.0 • C. This allowed the authors to conclude that the greenhouse height positively affects the natural ventilation phenomenon and allows to generate better thermal conditions inside the greenhouse [34].
Recently Fatnassi et al. [35] conducted a study in a tunnel type greenhouse with butterfly roof vents, where the thermal behavior of the greenhouse was evaluated numerically under four different gutter heights, 4, 5, 6 and 7 m, respectively. It was found that when the height is increased from 4 to 5 m, the temperature inside the greenhouse decreases by approximately 2 • C, while when it is increased from 6 to 7 m the temperature value only decreases by 0.2 • C. According to these results, the authors concluded that the effect of increasing greenhouse height on air temperature is clearly asymptotic, therefore, the increase in height should be analyzed as a climatic benefit versus economic cost ratio.
In this study, an experimentally validated CFD-2D model has been implemented under the climatic conditions of a region of the savanna of Bogota, Colombia. The CFD model was used to evaluate the effects on the air flow velocity and the interior temperature of a gothic-type greenhouse as the height level increases from 2.5 m to 5.0 m, under three ventilation configurations, side ventilation (R), roof vents (S) and combined ventilation (RS). Finally, the numerical data obtained in each of the developed simulation scenarios were grouped by ventilation configuration and fitted to response surface regression models. This was done to generate a more integrated and simpler method of analysis of the response variables analyzed and their relationship with the microclimate generated

Description of Prototype Proposed for Analysis
The analysis proposed in this research was carried out on a gothic multi-tunnel greenhouse with a 200 µm thick polyethylene cover located in the Savannah of Bogota, Colombia. The greenhouse was composed of a total of six spans, each span had a width length of 9.33 m, for a total greenhouse width length of 55.98 m ( Figure 1). The greenhouse had a gutter height of 4.0 m and a maximum roof height of 8.3 m. It was equipped with ventilation areas on the side walls, with an effective opening of 3.7 m, which equals a lateral ventilation surface with respect to the covered floor area of 13.21%. Likewise, the ventilation surface was complemented with a roof ventilation area of 1.5 m of total opening in each of the spans, therefore, the roof ventilation surface with respect to the covered surface is 16.1%.
Sustainability 2021, 13, x FOR PEER REVIEW 3 of 29 Regarding the greenhouse width, it is known that in Almeria type structures [31], which include the arc type [32] and gothic tunnel type [33], as the more spans that are joined laterally, the greater the width of the greenhouse, the ventilation rate is reduced exponentially, which in turn generates higher temperature values inside the greenhouse [31][32][33]. For the height dimension, Boulard and Fatnassi [34] using a CFD-2D study on an 8-span arch-type greenhouse found that when the height above the greenhouse gutter increases from 3 to 5 m, the reduction of the thermal gradient between the interior and the exterior is reduced from 3.5 to 2.0 °C. This allowed the authors to conclude that the greenhouse height positively affects the natural ventilation phenomenon and allows to generate better thermal conditions inside the greenhouse [34].
Recently Fatnassi et al. [35] conducted a study in a tunnel type greenhouse with butterfly roof vents, where the thermal behavior of the greenhouse was evaluated numerically under four different gutter heights, 4, 5, 6 and 7 m, respectively. It was found that when the height is increased from 4 to 5 m, the temperature inside the greenhouse decreases by approximately 2 °C, while when it is increased from 6 to 7 m the temperature value only decreases by 0.2 °C. According to these results, the authors concluded that the effect of increasing greenhouse height on air temperature is clearly asymptotic, therefore, the increase in height should be analyzed as a climatic benefit versus economic cost ratio.
In this study, an experimentally validated CFD-2D model has been implemented under the climatic conditions of a region of the savanna of Bogota, Colombia. The CFD model was used to evaluate the effects on the air flow velocity and the interior temperature of a gothic-type greenhouse as the height level increases from 2.5 m to 5.0 m, under three ventilation configurations, side ventilation (R), roof vents (S) and combined ventilation (RS). Finally, the numerical data obtained in each of the developed simulation scenarios were grouped by ventilation configuration and fitted to response surface regression models. This was done to generate a more integrated and simpler method of analysis of the response variables analyzed and their relationship with the microclimate generated

Description of Prototype Proposed for Analysis
The analysis proposed in this research was carried out on a gothic multi-tunnel greenhouse with a 200 µ m thick polyethylene cover located in the Savannah of Bogota, Colombia. The greenhouse was composed of a total of six spans, each span had a width length of 9.33 m, for a total greenhouse width length of 55.98 m ( Figure 1). The greenhouse had a gutter height of 4.0 m and a maximum roof height of 8.3 m. It was equipped with ventilation areas on the side walls, with an effective opening of 3.7 m, which equals a lateral ventilation surface with respect to the covered floor area of 13.21%. Likewise, the ventilation surface was complemented with a roof ventilation area of 1.5 m of total opening in each of the spans, therefore, the roof ventilation surface with respect to the covered surface is 16.1%.

Computational Fluid Dynamics Modeling
The numerical calculation of the physical processes involving fluid movement and related variables such as pressure, force, density, temperature and their associated changes in heat and mass transfer phenomena can be solved using computational fluid dynamics [36]. Currently, CFD is one of the most widely used analysis, design or redesign methodologies in different sectors of the economy, therefore, it is possible to find CFD studies approached from mechanical, environmental, aeronautical and automotive engineering [37]. On the other hand, its use in the agricultural sector has not been insignificant and in the last two decades it has been implemented as a tool for process analysis involving natural ventilation processes or energy optimization of greenhouse structures [22,38].
CFD analysis is a resolution technique that solves a set of nonlinear partial equations from algebraic discretization by means of numerical simulation using the finite volume method. A CFD simulation process is composed of three phases of development, preprocessing, processing and post-processing [39]. In the preprocessing phase, the computational domain, and the geometry of the structure on which the natural ventilation process will be analyzed are designed. In this phase the numerical meshing process applied to the whole computational domain including the greenhouse geometry is also generated, these processes must be performed based on criteria of good CFD simulation practices [22,40].
In the processing phase, the numerical solution of the evaluated problem is performed, when the air and energy flow in the model is calculated by solving the governing equations based on the physical laws of conservation of energy, mass and momentum. In this phase, the sub models considered as source terms of the governing equations are selected, such as turbulence, buoyancy or free convection, solar radiation, porous media and phenomena associated with mass transfer [41]. Finally, in the post-processing phase, it is possible to develop a qualitative and quantitative analysis of the numerical results obtained for each scenario evaluated [42].
It should be noted that the CFD-2D model used in this research has been previously successfully validated experimentally, verifying its high predictive capacity for air flow patterns and thermal distribution inside the greenhouse analyzed. The validation results can be verified in the works developed by Villagrán and Bojacá [16,43] and by Villagrán et al. [1]. Therefore, for this research work, details related to the validation process will not be discussed, although relevant aspects of each of the CFD simulation phases will be mentioned.

Pre-Processing
The commercial software ANSYS ICEM CFD (v. 18.0, ANSYS Inc., Canonsburg, PA, USA, EE. UU.) was used, by means of which a large computational domain was created in a two-dimensional configuration. It included the geometry of the cross section of the analyzed greenhouse, as well as the process of meshing the computational domain was performed by means of this software. The two-dimensional CFD studies of natural ventilation are useful and have a capacity to predict the ventilation rate and the thermal distribution, inside a structure when the dominant external wind currents blow perpendicular to the ventilation areas [27,44].
It was determined that the dimensions of the computational domain, establishing as the reference for the maximum height of the structure (H), should have dimensions from the region of the airflow inlet to the windward sidewall of the structure of 15H. From the leeward sidewall to the airflow outlet boundary of 20H and a minimum height from ground level of 10H ( Figure 2). These dimensions are established following the guidelines given in numerical studies of natural ventilation of greenhouses and are the dimensions that allow to adequately describe the physical phenomena associated with the analyzed problem within the atmospheric boundary layer [1,45]. In the meshing process, the computational domain volume was discretized into an unstructured grid of square or rectangular elements ( Figure 3). The size of the numerical grid was defined by a mesh independence test as described in the work done by Villagran et al. [1], at the end of the analysis process, the selected numerical grid was composed of a total of 1,104,369 elements. The quality parameters evaluated were cell size and cell-tocell size variation; it was found that 95.2% of the cells of the mesh were within the highquality interval (0.95-1). The orthogonality criterion was also evaluated, where the minimum value obtained was 0.92, a value considered as high quality [10]. Once the geometry has been constructed, it is possible to define the boundary conditions to be established in the limits of the computational domain and in the geometry of the greenhouse (Figure 1). In this specific case the right side was defined as the air flow inlet boundary, for which a condition determined by a logarithmic wind speed profile was established, this logarithmic profile depends on the climatic characteristics and the type of soil in the study region, factors that have been previously determined in the work developed by Villagrán et al. [1].
The left boundary on the other hand was determined as a pressure and airflow outflow boundary, while the upper boundary of the computational domain was set as a wall boundary condition with an imposed solar radiation flux. The wall-type boundary condition was also established for the lower boundary of the computational domain, the floor, the walls and the roof of the greenhouse, while the ventilation areas were set to the indoor boundary condition depending on the ventilation configuration analyzed. The physical and optical properties of the materials within the computational domain were also defined with the values summarized in Table 1. In the meshing process, the computational domain volume was discretized into an unstructured grid of square or rectangular elements ( Figure 3). The size of the numerical grid was defined by a mesh independence test as described in the work done by Villagran et al. [1], at the end of the analysis process, the selected numerical grid was composed of a total of 1,104,369 elements. The quality parameters evaluated were cell size and cell-to-cell size variation; it was found that 95.2% of the cells of the mesh were within the high-quality interval (0.95-1). The orthogonality criterion was also evaluated, where the minimum value obtained was 0.92, a value considered as high quality [10]. In the meshing process, the computational domain volume was discretized into an unstructured grid of square or rectangular elements ( Figure 3). The size of the numerical grid was defined by a mesh independence test as described in the work done by Villagran et al. [1], at the end of the analysis process, the selected numerical grid was composed of a total of 1,104,369 elements. The quality parameters evaluated were cell size and cell-tocell size variation; it was found that 95.2% of the cells of the mesh were within the highquality interval (0.95-1). The orthogonality criterion was also evaluated, where the minimum value obtained was 0.92, a value considered as high quality [10]. Once the geometry has been constructed, it is possible to define the boundary conditions to be established in the limits of the computational domain and in the geometry of the greenhouse (Figure 1). In this specific case the right side was defined as the air flow inlet boundary, for which a condition determined by a logarithmic wind speed profile was established, this logarithmic profile depends on the climatic characteristics and the type of soil in the study region, factors that have been previously determined in the work developed by Villagrán et al. [1].
The left boundary on the other hand was determined as a pressure and airflow outflow boundary, while the upper boundary of the computational domain was set as a wall boundary condition with an imposed solar radiation flux. The wall-type boundary condition was also established for the lower boundary of the computational domain, the floor, the walls and the roof of the greenhouse, while the ventilation areas were set to the indoor boundary condition depending on the ventilation configuration analyzed. The physical and optical properties of the materials within the computational domain were also defined with the values summarized in Table 1. Once the geometry has been constructed, it is possible to define the boundary conditions to be established in the limits of the computational domain and in the geometry of the greenhouse (Figure 1). In this specific case the right side was defined as the air flow inlet boundary, for which a condition determined by a logarithmic wind speed profile was established, this logarithmic profile depends on the climatic characteristics and the type of soil in the study region, factors that have been previously determined in the work developed by Villagrán et al. [1].
The left boundary on the other hand was determined as a pressure and airflow outflow boundary, while the upper boundary of the computational domain was set as a wall boundary condition with an imposed solar radiation flux. The wall-type boundary condition was also established for the lower boundary of the computational domain, the floor, the walls and the roof of the greenhouse, while the ventilation areas were set to the indoor boundary condition depending on the ventilation configuration analyzed. The physical and optical properties of the materials within the computational domain were also defined with the values summarized in Table 1. Subsequently, six greenhouse models (M1-M6) were created with the only variations being the height under the gutter (H1) and the lateral ventilation dimension (SV), as described in Table 2. The meshing process for models M1 to M6 was performed on the same geometric model already discretized from model M4, simply decreasing the greenhouse height for models M1 to M3 and increasing it for models M5 and M6.  Subsequently, six greenhouse models (M1-M6) were created with the only variations being the height under the gutter (H1) and the lateral ventilation dimension (SV), as described in Table 2. The meshing process for models M1 to M6 was performed on the same geometric model already discretized from model M4, simply decreasing the greenhouse height for models M1 to M3 and increasing it for models M5 and M6.

H1 (m) SV (m) Scheme
Finally, the simulation scenarios to be analyzed were defined, these were built from the ventilation configurations used locally, which are: ventilation by the lateral sides (S), by the roof vents (R) and combined ventilation (RS). After this, the greenhouse model to be evaluated was defined from M1 to M6 and finally the climatic conditions to be used as starting conditions for each simulation were determined. These conditions were known from previous studies developed in the same study region and where it was established that the maximum daytime temperature averaged 22.3 °C, the maximum solar radiation level was 853 Wm-2 and wind conditions can vary between values of 0.2 and 1.6 ms −1 . Therefore, it was determined to perform the evaluations under four wind speeds S1 (0.2 ms −1 ), S2 (0.5 ms −1 ), S3 (1 ms −1 ) and S4 (1.6 ms −1 ). According to the above, 72 simulations were performed and coded as shown in Table 3. Table 3. Coding of the 72 simulated scenarios.
Finally, the simulation scenarios to be analyzed were defined, these were built from the ventilation configurations used locally, which are: ventilation by the lateral sides (S), by the roof vents (R) and combined ventilation (RS). After this, the greenhouse model to be evaluated was defined from M1 to M6 and finally the climatic conditions to be used as starting conditions for each simulation were determined. These conditions were known from previous studies developed in the same study region and where it was established that the maximum daytime temperature averaged 22.3 • C, the maximum solar radiation level was 853 Wm-2 and wind conditions can vary between values of 0.2 and 1.6 ms −1 . Therefore, it was determined to perform the evaluations under four wind speeds S1 (0.2 ms −1 ), S2 (0.5 ms −1 ), S3 (1 ms −1 ) and S4 (1.6 ms −1 ). According to the above, 72 simulations were performed and coded as shown in Table 3. Table 3. Coding of the 72 simulated scenarios.

Processing
For the numerical solution of the Reynolds-averaged Navier-Stokes (RANS) equations, the ANSYS FLUENT processing software was used. The air flow in the greenhouse was considered to be turbulent and with a density change associated with temperature changes that can be modeled by the Boussinesq approximation. The method of resolution selected for the CFD model was semi-implicit for the pressure and velocity bound equations using the SIMPLE algorithm. Finally, residual convergence criteria were established in 10 −6 for the energy equation and in 10 −3 for the other variables such as continuity, turbulence and momentum [41].
The influence of solar radiation intensity on the spatial behavior of temperature was considered using the discrete order (DO) radiation model, which allows modeling the radiation and calculating the convective exchanges that occur in the computational domain, [46].
It was also decided not to include any type of crop, because the purpose of the study is not to analyze the heat and mass transfer flows that occur between the indoor environment and any species of plants. On the contrary, the main objective is to make an analysis of the ventilation phenomenon under the most critical condition that can occur in a real scenario and that is in empty conditions and under the most extreme climatic conditions [37,47]. This is an approach that is still valid and continues to be implemented in a significant number of studies carried out in various countries [48][49][50][51]. In addition, as different types of horticultural, ornamental, aromatic and medicinal crops are grown in the region, the analysis developed in this study can be applied to any type of crop.

Post-Processing
For this phase, where the main objective is to perform the qualitative and quantitative analysis of the evaluated scenarios, ANSYS CFD-Post software was used. Therefore, for each of the 72 simulations developed, two-dimensional plots of temperature and airflow patterns were generated and numerical values of air velocity and temperature in the greenhouse cross section at a height of 1.4 m above ground level were extracted. These data sets were analyzed for each simulation scenario and as a whole using response surface regression models.

Response Surface Modeling
Response surface regression models were fitted to each of the three ventilation scenarios under consideration (roof (R), side (S), and roof-side (RS)). For each scenario, models were fitted to evaluate the response of varying levels of exterior wind speed and greenhouse height on the internal temperature and air velocity. The objectives of applying response surface models on top of the CFD simulated scenarios was threefold: (1) to establish an approximate relationship between the response variables and the exterior wind speed and greenhouse height that can be used for predicting temperature or internal wind speed values for a given combination of the predictors; (2) to determine the statistical significance of the explanatory variables through hypothesis testing; and (3) to determine the optimum combination of the explanatory variables that result in the maximum response over the region under consideration.
Response surface models are an extension of linear models and works similarly; however, extra arguments should be added to consider the response-surface portion of the model [52]. In the present study, second order models including first-order response, two-way interactions and quadratic terms for the explanatory variables were included in the model's formulation. Once calibrated, the models were tested for their goodness of fit through measures such as the coefficient of determination, lack of fit and pure error [53,54].
The output of the models was plotted as perspective plots to depict the combined effect of the experimental factors (external wind speed and greenhouse height) over the dependent variables (internal temperature and air velocity). The models were fitted using the rsm package (v. 2.10.2; [52]) included in the statistical software R (v. 4.0.4; [55]). The airflow patterns under the ventilation configuration through the sidewall (S) areas of the greenhouse can be seen in Figure 4. In general terms, it is observed that the airflow pattern enters the greenhouse through the ventilation area on the windward side and moves horizontally in an airflow pattern that has a height very similar to the height under the gutter of each of the greenhouses evaluated (M1-M6), this horizontal flow leaves the greenhouse again through the ventilation area arranged on the leeward wall. Another factor that showed a direct relationship with height was the uniformity of the air flow velocity behavior along the cross axis of the greenhouse ( Figure 5). In general terms, it is observed that as the height of the structure increases, the velocity behavior curves show less oscillation between the points of high and low air velocity. This can also be verified with the standard deviation (sd) values. In S-M1S4 the sd value was ± 0.27 ms −1 while for S-M6S4 it was ± 0.16 ms −1 . This will undoubtedly promote a more homogeneous microclimate behavior inside the greenhouse, helping to limit the negative impacts on growth and development generated by non-uniform microclimates [58]. Another common feature is the air circulation loops that are generated in the upper part of each of the greenhouse spans, this occurs because due to the buoyancy phenomenon and the associated pressure components, the warm air inside the greenhouse moves vertically towards the roof area and its space in the lower part of each span is occupied by the fresh air that enters from the outside environment [3,56]. It should also be noted that qualitatively it is observed that as the height of the greenhouse increases, the air flow that moves over the nearby area where the crops are grown has a higher velocity; this same behavior is observed in the air flow that leaves the greenhouse through the leeward ventilation.
To quantitatively analyze the airflows, velocity data were extracted from a cross profile at a height of 1.4 m above ground level. In general, the airflows show a behavior for velocity that oscillates between a minimum of 0.25 ± 0.05 ms −1 for S-M1S1 and a maximum of 1.42 ± 0.16 ms −1 for S-M6S4. It was also observed that the air flows inside the greenhouse show a direct relationship with the wind speed outside, therefore, the higher the wind speed, the higher the air velocity inside the greenhouse. This occurs in greenhouses where this type of ventilation is analyzed and where porous insect-proof screens are not contemplated in the ventilation areas [29,57].
On the other hand, it is also observed that as the height of the greenhouse increases, the average air flow velocity inside the greenhouse increases for the wind speed scenarios analyzed (S1 to S4). In summary, the increase in velocity between the extreme simulation scenarios was 42.3, 20.9, 19.4 and 16.2% for M1S1, M1S2, M1S3, M1S4 compared to M6S1, M6S2, M6S3, M6S4, respectively. These velocity increases in the higher greenhouse models occur because the effect of air friction is reduced as the ratio (L/H1) between the greenhouse length (L) and the height above the gutter (H1) decreases, which coincides with what was previously reported in the study developed by Chu et al. [44].
Another factor that showed a direct relationship with height was the uniformity of the air flow velocity behavior along the cross axis of the greenhouse ( Figure 5). In general terms, it is observed that as the height of the structure increases, the velocity behavior curves show less oscillation between the points of high and low air velocity. This can also be verified with the standard deviation (sd) values. In S-M1S4 the sd value was ±0.27 ms −1 while for S-M6S4 it was ± 0.16 ms −1 . This will undoubtedly promote a more homogeneous microclimate behavior inside the greenhouse, helping to limit the negative impacts on growth and development generated by non-uniform microclimates [58]. 3.1.2. Roof Ventilation (R) Figure 6 shows the simulated airflow patterns for this ventilation configuration. In this case, since the side windows are closed, the airflow patterns are forced in and out through the windows in the roof region [16]. For the low velocity scenarios S1 where the thermal component of natural ventilation dominates, a behavior where three flow patterns are generated is observed. In span 1, an air flow was identified that enters through the ventilation area and becomes a convective movement loop between the floor and the roof of this span.
Subsequently, between span 2 and the middle of span 3 there is an airflow pattern 3.1.2. Roof Ventilation (R) Figure 6 shows the simulated airflow patterns for this ventilation configuration. In this case, since the side windows are closed, the airflow patterns are forced in and out through the windows in the roof region [16]. For the low velocity scenarios S1 where the thermal component of natural ventilation dominates, a behavior where three flow patterns are generated is observed. In span 1, an air flow was identified that enters through the ventilation area and becomes a convective movement loop between the floor and the roof of this span.
Subsequently, between span 2 and the middle of span 3 there is an airflow pattern that moves in the direction of the outside wind at ground level and a flow that moves in the opposite direction to the outside wind at the top of the roof of these spans. Finally, through the ventilation area of span 3 an airflow enters and mixes with another airflow entering from the ventilation area of span 4, flows that move towards the leeward wall and exits through the ventilation areas of spans 4 and 5. It should also be noted that the confluence of flows over the central area of span 3 generates a small low velocity loop near the floor region, being this an inadequate movement pattern that may cause the generation of a heat spot over this region [29,33,35].
For the case of the S4 scenario wind speed, the natural ventilation of the greenhouse will depend on the thermal and wind components together [59,60]. In this case, we were able to identify an air flow pattern that enters through span 1 and forms a recirculation pattern between the roof of the span, the floor and the wall of the windward side, where the highest air flow velocity occurs at ground level, which coincides with that reported by Kwon et al. [61]. Likewise, part of the air flow entering through the window of span 1 is mixed with another air flow entering through the window of span 2, which generates an acceleration of the air flow above ground level; the same behavior is repeated with the air flows entering through spans 3 and 4; finally, the air flows leave the greenhouse through spans 4 and 5. In numerical terms, the air flow velocity inside the greenhouse over the area where the crops are grown ranged from a minimum of 0.26 ± 0.07 ms −1 for R-M1S1 to a maximum of 0.92 ± 0.36 ms −1 in R-M1S4. In this specific case no increase in air flow velocity is observed as the greenhouse height increases in the case of the low velocity scenario S1 it is observed that between M1 and M6 there is only an increase of 11.5% but this same value is obtained between M1 compared to M2, M3 and M4.
For the case of S2 the air velocity increased between M1 and M2 by 2.5% but compared to the other greenhouse models it decreased by up to 12.5% with M6 specifically. In the S3 scenarios the air flow velocity decreased up to 20.5% compared M1 with M5 and Figure 6. Simulated air distribution patterns with ventilation configuration R and for the six greenhouse models (M1 to M6) under wind speeds S1 (0.2 ms −1 ) and S4 (1.6 ms −1 ).
In numerical terms, the air flow velocity inside the greenhouse over the area where the crops are grown ranged from a minimum of 0.26 ± 0.07 ms −1 for R-M1S1 to a maximum of 0.92 ± 0.36 ms −1 in R-M1S4. In this specific case no increase in air flow velocity is observed as the greenhouse height increases in the case of the low velocity scenario S1 it is observed that between M1 and M6 there is only an increase of 11.5% but this same value is obtained between M1 compared to M2, M3 and M4.
For the case of S2 the air velocity increased between M1 and M2 by 2.5% but compared to the other greenhouse models it decreased by up to 12.5% with M6 specifically. In the S3 scenarios the air flow velocity decreased up to 20.5% compared M1 with M5 and M6 and finally in S4 the air velocity also decreased up to 15% compared M1 with M6.
Regarding the behavior of the air flow in the cross axis of the greenhouse, it was found that for models M1, M2 and M3 there is a very similar distributed behavior in space, while for M4, M5 and M6 there is a small change in this spatial distribution, which may be caused by the change in greenhouse heights (Figure 7). Likewise, the inside air velocities increase as the wind speed increases, which has already been demonstrated in several research works [22,33,49,62]. On the other hand, there is a quite heterogeneous behavior among each of the six spans, finding that the smallest standard deviation among these scenarios was ± 0.07 ms −1 in the low velocity scenarios S1, while for the high velocity scenarios S4 this standard deviation was up to ± 0.36 ms −1 . The lowest air speed points were the regions near the leeward and windward wall and the highest air speed points were the areas below spans 2 and 4.

Side and Roof Ventilation (RS)
The spatial distribution of the airflow patterns for the combined ventilation configuration can be seen in Figure 8. Qualitatively, more continuous and intense airflows can be observed over the entire cross-sectional area analyzed, both in the crop and canopy regions, which should positively impact the greenhouse renovation rates and directly the thermal performance of the greenhouse [62]. On the other hand, there is a quite heterogeneous behavior among each of the six spans, finding that the smallest standard deviation among these scenarios was ±0.07 ms −1 in the low velocity scenarios S1, while for the high velocity scenarios S4 this standard deviation was up to ± 0.36 ms −1 . The lowest air speed points were the regions near the leeward and windward wall and the highest air speed points were the areas below spans 2 and 4.

Side and Roof Ventilation (RS)
The spatial distribution of the airflow patterns for the combined ventilation configuration can be seen in Figure 8. Qualitatively, more continuous and intense airflows can be observed over the entire cross-sectional area analyzed, both in the crop and canopy regions, which should positively impact the greenhouse renovation rates and directly the thermal performance of the greenhouse [62]. Numerically, the air flows presented a mean velocity that ranged from a minimum value of 0.31 ± 0.09 ms −1 in RS-M1S1 to a maximum of 1.53 ± 0.17 in RS-M6S4. For the case of S1, it is observed that the airflow velocity increases up to 29.1% in M5 and M6 compared to M1. Likewise, in S4 the increase in velocity was 15% in the M6 greenhouse with respect to the value obtained in M1, therefore, the greenhouse height has a significant impact on the airflow velocity. This is mainly due to the fact that, according to the analysis of natural ventilation carried out in previous studies, the greater the distance between the central axis of the side ventilation with respect to the central axis of the roof ventilation, the more dynamic the wind and thermal component is and therefore the better the air flow patterns are obtained [56,63].
On the other hand, the spatial behavior of the air flow velocity inside the greenhouse cross section can be found for some of the simulated scenarios in Figure 9. In general terms, it can be observed that as the height of the greenhouse increases, the less oscillation between the points of higher and lower velocity exists in the spatial distribution curves. At the same time, it should be noted again that these changes in air flow velocity are mainly due to the dynamics of incoming and outgoing air flows and to the interaction of warm and cold air masses that mix and generate buoyancy flows [10,64]. Although these values are less critical in the M3 to M5 greenhouse scenarios under the RC combined ventilation configuration. Numerically, the air flows presented a mean velocity that ranged from a minimum value of 0.31 ± 0.09 ms −1 in RS-M1S1 to a maximum of 1.53 ± 0.17 in RS-M6S4. For the case of S1, it is observed that the airflow velocity increases up to 29.1% in M5 and M6 compared to M1. Likewise, in S4 the increase in velocity was 15% in the M6 greenhouse with respect to the value obtained in M1, therefore, the greenhouse height has a significant impact on the airflow velocity. This is mainly due to the fact that, according to the analysis of natural ventilation carried out in previous studies, the greater the distance between the central axis of the side ventilation with respect to the central axis of the roof ventilation, the more dynamic the wind and thermal component is and therefore the better the air flow patterns are obtained [56,63].
On the other hand, the spatial behavior of the air flow velocity inside the greenhouse cross section can be found for some of the simulated scenarios in Figure 9. In general terms, it can be observed that as the height of the greenhouse increases, the less oscillation between the points of higher and lower velocity exists in the spatial distribution curves. At the same time, it should be noted again that these changes in air flow velocity are mainly due to the dynamics of incoming and outgoing air flows and to the interaction of warm and cold air masses that mix and generate buoyancy flows [10,64]. Although these values are less critical in the M3 to M5 greenhouse scenarios under the RC combined ventilation configuration.
Regarding air flow velocities, in general, under the three configurations evaluated, the values obtained are within the ranges reported in passive greenhouse ventilation studies [5,65]. Although for the low velocity condition S1, 100% of the airflows present velocities lower than the minimum recommended value of 0.5 ms −1 for plant growth inside greenhouses [42,65].  Regarding air flow velocities, in general, under the three configurations evaluat the values obtained are within the ranges reported in passive greenhouse ventilation stu ies [5,65]. Although for the low velocity condition S1, 100% of the airflows present vel ities lower than the minimum recommended value of 0.5 ms −1 for plant growth insi greenhouses [42,65].

Renewal Index (RI)
The renewal index (RI) was calculated using the method of integration of the vo metric flow of air leaving the ventilation areas of the greenhouse, the results obtained c be seen in Figure 10. The RI values ranged from a minimum of 5.89 Vol h −1 obtained in M1S1 to a maximum of 72.3 Vol h −1 obtained in S-M6S4. It is important to highlight th these are the contrasting scenarios, therefore, the minimum RI was obtained in the gree house with the lowest height and ventilation area, evaluated at the lowest wind spe (S1), while the maximum RI was obtained in the greenhouse with the highest height a ventilation area evaluated at the highest wind speed (S4).
For the lowest wind speed (S1), RI values between 5.89 Vol h −1 and 19.1 Vol h −1 w obtained in the side ventilation configuration S for M1 and M6, respectively. Likewise, the ventilation configuration through the vents of the roof region R, RI values of 19.4 a 19.7 Vol h −1 were obtained in M1 and M6, which are values compared to those obtain in S of 329% and 3.68%, respectively. Finally, for the RS combined ventilation configu tion, an RI value was obtained for M1 of 26.5 Vol h −1 which is a value 449% higher co pared to the S ventilation configuration and 36.6% higher compared to the R ventilati

Renewal Index (RI)
The renewal index (RI) was calculated using the method of integration of the volumetric flow of air leaving the ventilation areas of the greenhouse, the results obtained can be seen in Figure 10. The RI values ranged from a minimum of 5.89 Vol h −1 obtained in S-M1S1 to a maximum of 72.3 Vol h −1 obtained in S-M6S4. It is important to highlight that these are the contrasting scenarios, therefore, the minimum RI was obtained in the greenhouse with the lowest height and ventilation area, evaluated at the lowest wind speed (S1), while the maximum RI was obtained in the greenhouse with the highest height and ventilation area evaluated at the highest wind speed (S4).
For the lowest wind speed (S1), RI values between 5. To highlight the scenarios where RI values above or equal to the recommended minimum (RI ≥ 40 Vol h −1 ) are achieved for naturally ventilated agricultural structures [11,66]. Under ventilation configuration S, this was only obtained under greenhouse models M4, M5 and M6 and under external wind speed conditions of 1.5 ms −1 , with this same speed condition for ventilation configuration R adequate RI values are obtained in all greenhouse models (M1-M6).
While for the RS combined ventilation configuration, the RI values are adequate for all greenhouses (M1-M6) under wind speed scenarios higher than 1 ms −1 , the same is true for the M4, M5 and M6 models under wind speeds higher than 0.5 ms −1 . It is therefore relevant for greenhouse growers or greenhouse builders or decision makers to seriously analyze the prevailing wind speed conditions in the study region [62,67].
These results reaffirm some conclusions already determined in previous works related to the natural ventilation of greenhouses, where it was identified that the renewal indexes are dependent on the ventilation configuration used [68][69][70][71]. It was also identified that there is a linear relationship between the renewal index and wind speed [38,56,72,73], the ventilation configuration that allows to obtain the highest renewal index values is the combined configuration of roof and side vents [10,16,74,75].
On the other hand, in the side ventilation configuration, the increase of renewal indexes as a function of the increase of the side ventilation area is relevant and significant only in narrow greenhouses (width < 60 m) or with few attached spans (<6 spans) [76][77][78][79]. Finally, the IR in low external wind speed conditions are more stable and higher in greenhouse structures with relevant ventilation areas in the roof region [14,29,80].

Spatial Distribution of Temperature
The distribution of the thermal behavior for each of the evaluated scenarios can be

Side Ventilation (S)
The distribution of the thermal behavior for each of the evaluated scenarios can be seen in Figure 11. In general terms, irrespective of the outside wind speed, it was found that the cool zones correspond to the ventilation regions where the air enters the greenhouse, while the regions of higher temperature are located in the region where the air flow exits from the interior of the structure. This behavior occurs mainly because the air that enters the greenhouse through the windward window of span 1, as it crosses the cross section, mixes with the warm air and increases its energy level by heat transfer [10,45].
On the other hand, the height of the greenhouse directly influences the magnitude and spatial distribution of the temperature; it was observed that as the height of the greenhouse increases, the magnitude of the temperature and the percentage of the area of the structure with high temperature values decrease. This is an effect generated by the higher renovation index and the higher level of thermal inertia obtained in higher greenhouses [1,35]. Likewise, it is possible to observe the effects of the external wind speed, where for the highest speed scenario S4 it was found that the temperature value inside the greenhouse was lower and also presented a greater homogeneous behavior, coinciding with what was reported Flores-Velázquez [81].
In numerical terms, the average air temperature inside the greenhouse ranged between maximum values of 37.5 ± 8.9 °C for S-M1S1 and a minimum of 23.3 ± 0.6 °C for S-M6S4. For the S1 scenario it can be observed that the temperature decreases 11.7 °C comparing greenhouses M1 and M6 respectively, likewise for the S4 scenario this temperature reduction between M1 and M6 is only 2.5 °C, therefore, the increase in greenhouse height will be more relevant in regions where calms or low wind speeds predominate. Although under these conditions it should also be analyzed up to what level it is convenient from the technical and economic point of view to increase the greenhouse height, since for the M4 and M5 models under this condition no less important reductions of 9.7 and 10.8 °C, Figure 11. Simulated thermal distribution patterns with ventilation configuration S and for the six greenhouse models (M1 to M6) under wind speeds S1 (0.2 ms −1 ) and S4 (1.6 ms −1 ).
In numerical terms, the average air temperature inside the greenhouse ranged between maximum values of 37.5 ± 8.9 • C for S-M1S1 and a minimum of 23.3 ± 0.6 • C for S-M6S4. For the S1 scenario it can be observed that the temperature decreases 11.7 • C comparing greenhouses M1 and M6 respectively, likewise for the S4 scenario this temperature reduction between M1 and M6 is only 2.5 • C, therefore, the increase in greenhouse height will be more relevant in regions where calms or low wind speeds predominate. Although under these conditions it should also be analyzed up to what level it is convenient from the technical and economic point of view to increase the greenhouse height, since for the M4 and M5 models under this condition no less important reductions of 9.7 and 10.8 • C, respectively, would be obtained.
Another relevant criterion to be analyzed is the uniform distribution of temperature in the cross axis of the greenhouse. This is undoubtedly a factor that has received more attention in recent years, since it has a direct influence on the physiological behavior of the plants and on the final yield of the crops [58,82]. Therefore, for this study, the values obtained in the cross section of each of the greenhouses were plotted as a function of wind speed ( Figure 12).
In general, it is again observed the difference that exists between span 1 on the windward side (cool area) and span 9 on the leeward side (warm area) with relevant thermal differentials higher than 15 • C. Likewise, in some critical scenarios, the temperature in some areas of the cross section of the greenhouse exceeds 40 • C, a value that is quite inadequate for the growth and development of any vegetal species [83].

Roof Ventilation (R)
In qualitative terms, it is observed that under ventilation configuration R, lower te peratures are obtained with respect to configuration S ( Figure 13). In general, it is o served that the highest temperature areas are located on the windward and leeward sid of the span, respectively. This behavior is related to the ventilation configuration and the characteristics of the airflows of these regions near the greenhouse walls previou analyzed, a similar behavior was reported in a work on natural ventilation of four typ of multi-tunnel greenhouses by Park et al. [84].

Roof Ventilation (R)
In qualitative terms, it is observed that under ventilation configuration R, lower temperatures are obtained with respect to configuration S ( Figure 13). In general, it is observed that the highest temperature areas are located on the windward and leeward sides of the span, respectively. This behavior is related to the ventilation configuration and to the characteristics of the airflows of these regions near the greenhouse walls previously analyzed, a similar behavior was reported in a work on natural ventilation of four types of multi-tunnel greenhouses by Park et al. [84].
For this configuration it is also observed that there is a positive effect of wind speed and greenhouse height on the magnitude and distribution of the temperature inside the structure. In the most critical case S-M1S1 it is observed that there are relatively important areas of high temperature in spans 1, 3 and 5, while in the S-M6S1 scenario these areas are significantly reduced. In S-M1S4 the high temperature area disappears in span 3 and the hot areas in spans 1 and 6 are reduced to an area very close to the side wall of each side and, finally, in S-M6S4 these high temperature areas disappear completely in the spans already described. In numerical terms the average air temperature ranged between maximum value of 32.7 ± 3.1 °C for R-M1S1, value that is lower by 4.8 °C with respect to S-M1S1, the minimum value obtained was 25.0 ± 1.3 °C in R-M6S4 value that is higher by 1.7 °C with respect to S-M6S4. For the S1 scenario the temperature reduction obtained was 1.3, 2.2, 2.8, 3.6 and 4.3 °C in M2 to M6 compared to M1 respectively, while in S4 these reductions with respect to M1 were 0.5, 0.7, 1.1, 1.1, 1.4, 1.7 °C in M2 to M6, successively, results that continue to confirm the importance of roof ventilation in low wind speed conditions (S1).
For this case, the thermal behavior on the greenhouse cross axis in each scenario evaluated shows a totally different spatial temperature distribution than that observed in the previous ventilation configuration (Figure 14). Due to the air inlet and outlet flows through the roof vents and the airflow distribution patterns discussed in Figure 6, it is possible to observe the temperature variations occurring between span 1, 3 and 6 successively. It can also be seen how this temperature distribution tends to be more uniform and smaller in magnitude as the height of the greenhouse and the outside wind speed increase. In numerical terms the average air temperature ranged between maximum value of 32.7 ± 3.1 • C for R-M1S1, value that is lower by 4.8 • C with respect to S-M1S1, the minimum value obtained was 25.0 ± 1.3 • C in R-M6S4 value that is higher by 1.7 • C with respect to S-M6S4. For the S1 scenario the temperature reduction obtained was 1.3, 2.2, 2.8, 3.6 and 4.3 • C in M2 to M6 compared to M1 respectively, while in S4 these reductions with respect to M1 were 0.5, 0.7, 1.1, 1.1, 1.4, 1.7 • C in M2 to M6, successively, results that continue to confirm the importance of roof ventilation in low wind speed conditions (S1).
For this case, the thermal behavior on the greenhouse cross axis in each scenario evaluated shows a totally different spatial temperature distribution than that observed in the previous ventilation configuration (Figure 14). Due to the air inlet and outlet flows through the roof vents and the airflow distribution patterns discussed in Figure 6, it is possible to observe the temperature variations occurring between span 1, 3 and 6 successively. It can also be seen how this temperature distribution tends to be more uniform and smaller in magnitude as the height of the greenhouse and the outside wind speed increase.
ainability 2021, 13, x FOR PEER REVIEW 20 o Figure 14. Spatial distribution of air temperature inside the greenhouse for the simulated scenar under the R configuration.

Side Ventilation and Roof (RS)
The spatial distribution of the temperature for this configuration allows us to obser qualitatively that the temperature value in each of the simulated scenarios has a low magnitude compared to the S and R configurations ( Figure 15). For these scenarios, i observed that the regions of lower temperature coincide with the ventilation areas wh the air flow enters, both in the lateral sides and in the roof ventilation areas, which co cides with what was analyzed in the work by Villagran et al. [10].

Side Ventilation and Roof (RS)
The spatial distribution of the temperature for this configuration allows us to observe qualitatively that the temperature value in each of the simulated scenarios has a lower magnitude compared to the S and R configurations ( Figure 15). For these scenarios, it is observed that the regions of lower temperature coincide with the ventilation areas where the air flow enters, both in the lateral sides and in the roof ventilation areas, which coincides with what was analyzed in the work by Villagran et al. [10].
For scenario S1, it can be observed that the high temperature regions are in the spans located on the leeward side. This high temperature region gradually disappears as the height of the greenhouse increases. In the case of M1, the highest temperature region occurs in spans 4, 5 and 6, with a heat spot in the middle of the area of span 6. However, in M6, this heat patch disappears from the span 6 and only a small region of higher temperature is observed near the ventilation area on the leeward side of the greenhouse. Similar behavior is observed in scenario S4, although in this case the temperature values appear to be qualitatively lower than those obtained in S1.
Numerically, it was found that the temperature value inside the structure presented a maximum value of 29.4 ± 3.8 • C in RS-M1S1, being this value lower by 8.1 and 3.3 • C compared to S and R, respectively. Therefore, the RS ventilation configuration also allows obtaining a higher cooling efficiency under the same climatic conditions, which is a very relevant factor for the management of the microclimate in greenhouses located in the savanna of Bogota [1].
On the other hand, the minimum mean inside air temperature value was 23.2 ± 0.6 • C in RS-M6S4, a value that is lower by 1.7 and 0.1 • C with respect to those obtained in R and S, respectively. Therefore, again it can be observed that the greatest benefits in terms of thermal distribution are obtained for low wind speed conditions. It is also important to note that the temperature values obtained in most of the scenarios are within the range of 25 and 30 • C, ranges where species of commercial and food interest such as Tomatoes tend to develop adequately. As well as some ornamental species of commercial interest for the international market such as the Rose or Carnation [1,85].
Regarding the spatial distribution, it was again found that as the greenhouse height increases and the external wind speed increases, more stable and uniform temperature values are obtained in the cross section of the greenhouse (Figure 16). The thermal gradients for the same moment inside the greenhouse were reduced from 13.2 • C in M1S1 to 1.94 • C in M6S4, the latter value being a recommended limit to guarantee the homogeneity of a naturally ventilated greenhouse [39,86]. For scenario S1, it can be observed that the high temperature regions are in the spans located on the leeward side. This high temperature region gradually disappears as the height of the greenhouse increases. In the case of M1, the highest temperature region occurs in spans 4, 5 and 6, with a heat spot in the middle of the area of span 6. However, in M6, this heat patch disappears from the span 6 and only a small region of higher temperature is observed near the ventilation area on the leeward side of the greenhouse. Similar behavior is observed in scenario S4, although in this case the temperature values appear to be qualitatively lower than those obtained in S1.
Numerically, it was found that the temperature value inside the structure presented a maximum value of 29.4 ± 3.8 °C in RS-M1S1, being this value lower by 8.1 and 3.3 °C Figure 15. Simulated thermal distribution patterns with the RS ventilation configuration and for the six greenhouse models (M1 to M6) under wind speeds S1 (0.2 ms −1 ) and S4 (1.6 ms −1 ).
Regarding the spatial distribution, it was again found that as the greenhouse heig increases and the external wind speed increases, more stable and uniform temperatu values are obtained in the cross section of the greenhouse (Figure 16). The thermal gra ents for the same moment inside the greenhouse were reduced from 13.2 °C in M1S1 1.94 °C in M6S4, the latter value being a recommended limit to guarantee the homogene of a naturally ventilated greenhouse [39,86].

Response Surface Modeling
Finally, looking for a better and simpler interpretation of the results obtained by C simulation, these numerical results of the evaluated scenarios were adjusted to respon surface models. This methodology allows observing the behavior of the response varia in scenarios not simulated numerically, if these scenarios are within the limiting ranges the simulated boundary conditions The results for the fitted response surface models are summarized in Table 4. T second order model alternative gave the best results for all models. The lack of fit t showed that the models accurately fit the data exception made for the sidewalls scena and with temperature as the response variable. The adjusted R-squared coefficient of d termination indicated varying results whether the outcome variable was the internal velocity or the temperature. Highest R-squared values were obtained for the models w wind speed as the dependent variable as compared to their temperature pairs under same scenario.

Response Surface Modeling
Finally, looking for a better and simpler interpretation of the results obtained by CFD simulation, these numerical results of the evaluated scenarios were adjusted to response surface models. This methodology allows observing the behavior of the response variable in scenarios not simulated numerically, if these scenarios are within the limiting ranges of the simulated boundary conditions The results for the fitted response surface models are summarized in Table 4. The second order model alternative gave the best results for all models. The lack of fit test showed that the models accurately fit the data exception made for the sidewalls scenario and with temperature as the response variable. The adjusted R-squared coefficient of determination indicated varying results whether the outcome variable was the internal air velocity or the temperature. Highest R-squared values were obtained for the models with wind speed as the dependent variable as compared to their temperature pairs under the same scenario.
Despite the varying results for the adjusted R-squared coefficient, all terms included in the models were considered significant. Since most of the results indicated a good fit of the model to the calibration data, we evaluated the effect of varying levels of greenhouse height and external wind speed for each ventilation scenario through response surface plots. The response surface plots indicated a similar behavior for internal temperature and air velocity irrespective of the ventilation scenario ( Figure 17). The individual effects of the predictors showed that increasing the greenhouse height effectively resulted in lower temperatures, while the internal air velocity was slightly affected by the greenhouse height. The major effect of the external wind speed on both, the temperature, and the internal air velocity, is clearly depicted in Figure 17. Increasing the greenhouse height under increasing levels of external wind speed resulted in lower temperature values within the greenhouse. On the other hand, lower greenhouse heights combined with decreased external wind speeds resulted in lower internal air velocity. However, while all scenarios depicted the same trends, differences arise when looking at the scale of variation of the internal temperatures and air velocity. The sidewalls scenario showed that increasing the greenhouse height under increased external wind speeds have a more dramatic effect in terms of lowering the internal temperature of the greenhouse (Figure 17c), compared to roof or combined ventilation (Figure 17a,e).
The option to apply only the roof ventilation decreased the internal wind speed under almost null external wind speed no matter what the greenhouse height is. However, when the external wind speed increases along with the greenhouse height, internal air velocity also increases but the highest internal air velocity is reached with the lowest greenhouse height and the maximum external wind speed (Figure 17b). This results is opposed to the other scenarios where the top internal air velocity was reached when the maximum greenhouse height and external wind speed were evaluated (Figure 17d,f).
The optimal combination of greenhouse height and external wind speed that minimized the temperature and achieved an adequate internal wind speed seemed to be located around the maximum height and external wind speed for all scenarios. These results implied that higher values for both explanatory variables should be considered to find the lowest internal temperature and the highest internal air velocity. However, going beyond and increasing the region of analysis will yield unrealistic results due to greenhouse construction restrictions and maximum local wind speeds.
The results presented here are in line with those reported by Villagran et al. [1], and Bustamante et al. [87], indicating that increased greenhouse heights and the combination of roof and sidewalls ventilation yields improved internal climate conditions, particularly for greenhouse established in high altitude tropical regions.
It is also important to recommend that a future study regarding this greenhouse typology should focus on the structural design of the model evaluated and an economic analysis of the height alternatives proposed. This to be able to dimension the robustness of the structure and an associated cost per m 2 of covered area for each of the 6 models proposed, although in Colombia there are currently no regulations in force for the construction of greenhouses, a work of this magnitude can provide economic and structural tools for decision making. This information can be used for the generation of public policy management documents that can promote and regulate the use of greenhouses conceived with design criteria.

Conclusions
The numerical simulation tool is still an agile and accurate alternative to determine the renewal indexes and thermal distribution in passive greenhouses. It is also analyze unconstructed scenarios, facilitating decision making for farmers and providing an alternative analysis that can contribute to improve the technical and economic sustainability of protected agriculture.
The use of previously validated CFD models allows generating simulated data series of temperature and air velocity inside a greenhouse and associating them to response variables such as the height of the structure evaluated or others. This facilitated the association of the obtained data to other analysis methodologies such as regression with response surface models obtaining acceptable adjustments in the analyzed variables, therefore, it was possible to generate response surface graphs which facilitates the interpretation of the numerical results as a whole for each ventilation configuration analyzed (S, R, RS).
The natural ventilation of a greenhouse is highly dependent on the behavior of the wind speed of the study region. For the three ventilation configurations analyzed (S, R and RS), it was found that the renewal indexes were maximized as a function of the increase in the external wind speed. Although it should be noted that regardless of the ventilation configuration for the lowest wind speed S1 (0.2 ms −1 ), none of the greenhouse models exceeded the recommended minimum value of 40 hourly renovations, although these IR values are less critical in the RC configuration. The optimization of the IR in the RS ventilation configuration allowed obtaining a higher cooling efficiency inside the greenhouse and a uniformity in the spatial distribution of the temperature. Maximum temperature values of 29.4 ± 3.8 • C were obtained in RS-M1S1, being this value lower by 8.1 and 3.3 • C with respect to S-M1S1 and R-M1S1 respectively.
The crop production in greenhouse in the future should be with the efficient use of resources, so that numerical simulation techniques will have the goal of adapting the climatic environment in the design and operation of protected agriculture. In addition to the greenhouse dimensions, the position and size of windows are among the other influential factors to be considered in future studies Funding: The present study was funded by the Servicio Nacional de Aprendizaje (SENA), la Asociación Colombiana de Exportadores de Flores (Asocolflores) y al Centro de Innovación de la Floricultura Colombiana (Ceniflores) through the project "Generación de una herramienta de diseño u optimización de ventilación natural de los invernaderos dedicados a la producción de flores de corte en cuatro subregiones de la Sabana de Bogotá, mediante el uso de herramientas de simulación basadas en la técnica Dinámica de Fluidos Computacional (CFD)". Data Availability Statement: Not applicable.