Application of Semi-Empirical Ventilation Models in A Mediterranean Greenhouse with Opposing Thermal and Wind E ﬀ ects. Use of Non-Constant C d (Pressure Drop Coe ﬃ cient Through the Vents) and C w (Wind E ﬀ ect Coe ﬃ cient)

: The present work analyses the natural ventilation of a multi-span greenhouse with one roof vent and two side vents by means of sonic anemometry. Opening the roof vent to windward, one side vent to leeward, and the other side vents to windward (this last vent obstructed by another greenhouse), causes opposing thermal G T (m 3 s − 1 ) and wind e ﬀ ects G w (m 3 s − 1 ), as outside air entering the greenhouse through the roof vent circulates downward, contrary to natural convection due to the thermal e ﬀ ect. In our case, the ventilation rate R M (h − 1 ) in a naturally ventilated greenhouse ﬁts a second order polynomial with wind velocity u o ( R M = 0.37 u o 2 + 0.03 u o + 0.75; R 2 = 0.99). The opposing wind and thermal e ﬀ ects mean that ventilation models based on Bernoulli’s equation must be modiﬁed in order to add or subtract their e ﬀ ects accordingly—Model 1, in which the ﬂow is driven by the sum of two independent pressure ﬁelds G M 1 = (cid:113)(cid:12)(cid:12)(cid:12) G 2 T ± G 2 w (cid:12)(cid:12)(cid:12) , or Model 2, in which the ﬂow is driven by the sum of two independent ﬂuxes G M 2 = | G T ± G w | . A linear relationship has been obtained, which allows us to estimate the discharge coe ﬃ cient of the side vents ( C dVS ) and roof vent ( C dWR ) as a function of u o [ C dVS = 0.028 u o + 0.028 ( R 2 = 0.92); C dWR = 0.036 u o + 0.040 ( R 2 = 0.96)]. The wind e ﬀ ect coe ﬃ cient C w was determined by applying models M1 and M2 proved not to remain constant for the di ﬀ erent experiments, but varied according to the ratio u o / ∆ T io0.5 or δ [ C wM1 = exp( − 2.693 + 1.160 / δ ) ( R 2 = 0.94); C wM2 = exp( − 2.128 + 1.264 / δ ) ( R 2 = 0.98)].


Introduction
Natural ventilation is perhaps the main means of climate control in greenhouses [1,2], particularly in regions such as the Mediterranean, where new technologies have not been widely incorporated [3]. During most of the year, good management of natural ventilation may prove sufficient to maintain suitable levels of temperature, humidity, and CO 2 concentration inside the greenhouse [1]. It is important to have both quantitative and qualitative knowledge of how natural ventilation functions in greenhouses in order to correctly design and use the vents [4].
The earliest studies on the circulation of air in greenhouses date back to the mid-20th century [5,6]. Since then, many authors have shown interest in studying and understanding natural ventilation than 1, while Bot [21] put this limit at 0.3.
The relationship between wind and thermal effects must be known in order to correctly design greenhouse ventilation. Few authors mention this fact, and most works apply models that take these effects as being complementary to one another. This is not always the case, however, as depending on the arrangement of ventilation openings, wind can assist or oppose the thermal force in natural ventilation. It is most important to bear this in mind when applying semi-empirical models of ventilation [22].
To study the natural ventilation in greenhouses, the gas tracer method can allow us to quantify ventilation rate correctly, but not characterise the real airflow inside the greenhouse [1]. One disadvantage of the tracer gas method is that the position of the sensors can affects the results obtained, as Van Buggenhout et al. [23] observed, with errors of up to 86% according to the position of the sensors.
With sonic anemometers, we can measure the different components of the air velocity vector at the greenhouse vents, determining where the air enters or exits the greenhouse through each vent. This methodology allows us to study when the wind and thermal effects are complementary or opposing. The objective of the present work is to study natural ventilation in a greenhouse with two side vents and one roof vent with the use of sonic anemometry. The roof vent was windward, the windward side vent was protected from the wind by another greenhouse located in close proximity, and the leeward side vent was free of obstacles. Under such conditions, wind and thermal effects maybe fail to contribute to natural ventilation in the same manner. Rather, they can have opposing effects.
In this work, different ventilation models have been applied in a greenhouse with opposition of the wind and thermal effects. The models have been modified to take into account the opposition of both effects. Identifying which model best fits the experimental values was done by statistical analysis. In the western sector of the greenhouse, the vent surface area was 1.05 × 17.5 m 2 for each of the two-side vents (2.1 m in height from the ground to the mid-point of the vent) and 0.97 × 17.5 m 2 for the roof vent (6.2 m in height from the ground to the mid-point of the vent), and 4.1 m of vertical distance between the mid-point of the side and the roof vents. The ventilation surface SV was 11.2% of the greenhouse base area (SV/SA = 0.112). To prevent insects entering, all the greenhouse vents were fitted with insect-proof screens of 13 × 30 threads cm −2 (0.39 porosity; 164.6 µm pore width; 593.3 µm pore height; 165.5 µm thread diameter). Measurement tests were carried out under prevailing winds from the southwest (Figure 1), popularly known as "Poniente wind" (in the province of Almería). The outside climatic conditions remained relatively stable over the ten measurement tests ( Table 1). The greenhouse contained a tomato crop (Solanum lycopersicum L. var. cerasiforme Hort., cv. Salomee) with an average height of approximately 0.85 m for the first test and 1.88 m for the final one, and a leaf area index (m 2 leaf per m 2 ground.) of about 0.37 and 2.75, respectively. The measurements were performed from 7 April 2009 to 2 July 2009. Airflow was measured at 21 points (tests 1 to 4) and 12 points (tests 5 to 10), evenly distributed at the side vents (Figure 2a,b). To measure air velocity at the side vents, two 3D sonic anemometers were used (Figure 2c), recording data over three minutes at each point [19].
Given the difficulty of placing the sonic anemometers in the roof vent, it was divided into three equal surfaces and air velocity measurements were taken continuously at the centre of each one, with two (measurement tests 1 to 4) or one (measurement tests 5 to 10) 2D sonic anemometers ( Figure  2a,b). At the roof vent, these measurements were taken continuously using six (tests 1 to 4) or three (tests 5 to 10) 2D sonic anemometers (Figure 2c). In the western sector of the greenhouse, the vent surface area was 1.05 × 17.5 m 2 for each of the two-side vents (2.1 m in height from the ground to the mid-point of the vent) and 0.97 × 17.5 m 2 for the roof vent (6.2 m in height from the ground to the mid-point of the vent), and 4.1 m of vertical distance between the mid-point of the side and the roof vents. The ventilation surface S V was 11.2% of the greenhouse base area (S V /S A = 0.112). To prevent insects entering, all the greenhouse vents were fitted with insect-proof screens of 13 × 30 threads cm −2 (0.39 porosity; 164.6 µm pore width; 593.3 µm pore height; 165.5 µm thread diameter). Measurement tests were carried out under prevailing winds from the southwest (Figure 1), popularly known as "Poniente wind" (in the province of Almería). The outside climatic conditions remained relatively stable over the ten measurement tests ( Table 1). The greenhouse contained a tomato crop (Solanum lycopersicum L. var. cerasiforme Hort., cv. Salomee) with an average height of approximately 0.85 m for the first test and 1.88 m for the final one, and a leaf area index (m 2 leaf per m 2 ground.) of about 0.37 and 2.75, respectively. The measurements were performed from 7 April 2009 to 2 July 2009. Airflow was measured at 21 points (tests 1 to 4) and 12 points (tests 5 to 10), evenly distributed at the side vents (Figure 2a,b). To measure air velocity at the side vents, two 3D sonic anemometers were used (Figure 2c), recording data over three minutes at each point [19].
Given the difficulty of placing the sonic anemometers in the roof vent, it was divided into three equal surfaces and air velocity measurements were taken continuously at the centre of each one, with two (measurement tests 1 to 4) or one (measurement tests 5 to 10) 2D sonic anemometers (Figure 2a,b). At the roof vent, these measurements were taken continuously using six (tests 1 to 4) or three (tests 5 to 10) 2D sonic anemometers (Figure 2c).   The distribution of measurement points used in the present study was very similar to the one made in Espinoza et al. [24], in the same experimental greenhouse. For lateral vents, the surface corresponding to each point was 0.9 m 2 (tests with 21 points) and 1.5 m 2 (tests with 12 points). These are similar values to those used by other authors. Boulard et al. [25] used 2.6 m 2 per point in a tunnel greenhouse with one roof vent, Teitel et al. [18] used 1.1 m 2 per point in a mono-span greenhouse with two side vent openings, and Molina-Aiz et al. [19] used 2.1 m 2 per point in a five-span Almería-type greenhouse. At the roof vent, the mean surface corresponding to each section was 5.8 m 2 . Teitel et al. [26] used a larger surface area for each measurement point, with 8.5 m 2 per point in a four-span greenhouse with three roof vents.

Equipment and Instrumentation
The three components of air velocity and air temperature were measured with two 3D sonic anemometers. The horizontal components of air velocity (x and y) were measured with six 2D sonic anemometers ( Table 2). The data measured by the sonic anemometers were recorded by two CR3000 Microloggers (Campbell Scientific Spain S.L.). The data registration frequency was 10 Hz [27] for 3D sonic anemometers and 1 Hz for 2D sonic anemometers, respectively.
Outside climatic conditions were recorded by a meteorological station at a frequency of 0.5 Hz (Figure 1). It included a measurement box with a Pt1000 temperature sensor and a capacitive humidity sensor (BUTRON II). Wind speed and direction were recorded by a cup anemometer and a vane. Solar radiation was measured using a Kipp Solari sensor. Air temperature and humidity inside the greenhouse were measured using six autonomous dataloggers HOBO Pro Temp-RH U23-001 ( Table 2). The dataloggers were placed in a vertical profile under the ridge of the three greenhouse spans at heights of 1 and 2 m (Figure 1). These autonomous dataloggers were protected against direct solar radiation with a passive solar radiation open shield.

Ventilation Models
In general, we can consider that the natural ventilation in a greenhouse is the result of combining two different fluxes: one generated by wind and the other generated by buoyancy forces. We can combine these fluxes and obtain different semi-empirical models based on Bernoulli's equation [4,7,19,28]. These models (used by Molina-Aiz et al. [19]) can be modified by applied using a different discharge coefficient for the side vents (C dVS ) and for the roof vent (C dVR ): (1) Model 1. In this model, the flow is driven by the sum of two independent pressure fields. In a greenhouse with side and roof openings, G is calculated as the vector sum of the free component of the flux induced by buoyancy forces G T and the flux induced by wind forces G w : (3) Model 3. In this model, only the wind effect is considered, neglecting the thermal effect: The height of the meteorological station and its distance from the greenhouse condition the results obtained when applying the models. This aspect is important to keep in mind when applying the models described in the literature.

Anemometric Measurement of Volumetric Flow Rate
To describe the air circulation through the greenhouse vents, the mean and turbulent volumetric flow rates can be calculated considering only the component of air velocity perpendicular to the vent, u x [25]: where u x,j are the time average value of air velocity perpendicular to the vent, |u' x,j | are the time average value of the absolute values of the turbulent component, and S V,j are elementary surface corresponding to each measurement point. With only two 3D sonic anemometers measuring at the lateral vents, at different intervals of 3 min (Figure 2), the external conditions can change or fluctuate during the test. This problem can be overcome by correcting (scaling with the wind speed) the air velocities measured at each position j at the lateral vents [19]. The corrected air velocity u * x,j (t) has been calculated by multiplying measured values of u x,j (t), at time t and at each point j, by the ratio between the average wind speed u o for the overall test period, and the measured values of u o (t) at time t [19]:

Estimation of the Pressure Drop/discharge Coefficient of the Greenhouse Openings
The individual pressure drop coefficient of each vent C d was determined following the methodology outlined by Molina-Aiz et al. [19]. Calculation of C d requires knowing the specific permeability (K p = 1.851·10 −9 ) and inertial factor (Y = 0.155) of the insect-proof screens installed in the vents; these values were determined in wind-tunnel tests. A detailed description of the wind-tunnel experiment was provided by Valera et al. [29] or López et al. [30]. To determine K p and Y, the thickness of the screen is required (e = 391.7 µm), which was obtained by measuring a transversal section of the net with an optical measurement unit equipped with a video system (TESA-VISIO 300, TESA SA, Switzerland; with a resolution of 0.05 µm and precision of ±7 µm).

Statistical Analysis
We carried out regression analyses to compare the different variables for statistically significant relationships (p-value < 0.05) using Statgraphics ® Centurion 18 v18.1 (Statgraphics Technologies, Inc., The Plains, VA, USA). In order to analyse the fit of the average volumetric flow rate for the greenhouse G M,sim (simulated by models M1, M2 and M3) with the experimentally observed values of G M,obs , different statistics were used, as well as the determining coefficient R 2 . Two of the most commonly used statistics based on the deviation of G M,sim from G M,obs are RMSD (root mean squared deviation) and MD (or bias) [31]: RMSD represents the average distance between the simulated and observed values. MD corresponds to the mean value of the differences between the simulated and observed values, though in this case negative differences are compensated by positive ones, which may lead to erroneous interpretation. Both statistics represent different aspects of the deviation from simulated values, but the relationship between them has not been well defined [31]. The Nash-Sutcliffe efficiency (NSE) is a normalised statistic indicating the relative magnitude of the residual variance of the model ("noise") with respect to the variance of the observed values ("information") [32]. NSE indicates how well the plot of the observed versus simulated values fits the 1:1 line [33]: where G M,obs is the average value of all the G M,obs values. NSE can take values of between −∞ and 1.0, the latter being the optimal value. Values between 0.0 and 1.0 are generally viewed as acceptable levels of performance [33]. Percent bias PBIAS represents the mean trend of simulated values to greater or lower than their respective observed values [34]. The optimal value of PBIAS is 0.0, and values close to 0 indicate high precision of the model. Positive values indicate that the model provides values that are Agronomy 2019, 9, 736 9 of 23 lower than those observed (model underestimation bias), while negative ones indicate the contrary (model overestimation bias) [34]. PBIAS can be expressed as a percentage: RMSD-observations standard deviation ratio (RSR): this statistic [33] was developed based on the recommendations of [35]. RSR standardises the values of RMSD with standard deviation from the observed values. It combines both an error index and the additional information recommended by Legates and McCabe [36]: RSR takes values from 0 to a large positive value, and when RMSD is equal to zero, the model is considered perfect. Small values of RSR indicate a better performance of the simulation model [33].

Airflow Inside the Greenhouse
Given the particular situation of the greenhouse and the location of the vents, in conditions of natural ventilation with prevailing southwest wind (SW) or "Poniente winds", the eolic and thermal effects oppose each other. The former causes air to enter through the windward roof vent and to leave through the leeward and windward side vents (the windward side vent was protected from the wind by another greenhouse located in close proximity). The thermal effect, on the other hand, causes warm air to rise and leave through the roof vent, which favours the entrance of air through the side vents. Figure 3a,b shows the entrance and exit of air through the roof vent (see the polar histograms), demonstrating the negative interaction of the wind and thermal effects. This flow pattern with opposition of the wind effect and thermal effect in the roof vents, in the same naturally-ventilated three-span Mediterranean greenhouse, was described in more detail, including measurements among the crop lines, in López et al. [37]. And, similar flow patterns, in the same greenhouse and with the same methodology, was described in Espinoza et al. [24], opening the lateral vents combined with two and three roof vents.
The greenhouse natural ventilation rate is affected by the thermal effect due to the inside-outside temperature difference, ∆T io . For greenhouses with side and roof vents, Kittas et al. [4] established that the thermal effect is important when the ratio u o / ∆T io 0.5 < 1, while Bot [21] put this limit at 0.3.
In tests 1, 2, 3, 7, and 10, with u o / ∆T io 0.5 ≥ 1.5, hardly any air left the greenhouse through the roof vent, indicating that the wind effect clearly prevailed over the thermal effect, and the air left in quite a uniform fashion through the side vents ( Figure 4). In the other tests, with u o / ∆T io 0.5 < 1.5, the opposing wind and thermal effects gave rise to less uniformity of air flow at the ventilation surfaces (Figure 3a,b), with alternating positive (entrance) and negative (exit) airflow at the roof vent (see polar histograms). In some of the tests, there was a certain degree of discrepancy between the wind direction and the airflow direction entering the greenhouse through the roof vent. This is likely due to the location of the meteorological station, to the north of the greenhouse, which therefore recorded the wind characteristics once it had passed through the experimental greenhouse.
demonstrating the negative interaction of the wind and thermal effects. This flow pattern with opposition of the wind effect and thermal effect in the roof vents, in the same naturally-ventilated three-span Mediterranean greenhouse, was described in more detail, including measurements among the crop lines, in López et al. [37]. And, similar flow patterns, in the same greenhouse and with the same methodology, was described in Espinoza et al. [24], opening the lateral vents combined with two and three roof vents. The greenhouse natural ventilation rate is affected by the thermal effect due to the inside-outside temperature difference, ∆Tio. For greenhouses with side and roof vents, Kittas et al. [4] established that the thermal effect is important when the ratio uo/∆Tio 0.5 <1, while Bot [21] put this limit at 0.3. In tests 1, 2, 3, 7, and 10, with uo/∆Tio 0.5 ≥1.5, hardly any air left the greenhouse through the roof vent, indicating that the wind effect clearly prevailed over the thermal effect, and the air left in quite a uniform fashion through the side vents ( Figure 4).   The greenhouse natural ventilation rate is affected by the thermal effect due to the inside-outside temperature difference, ∆Tio. For greenhouses with side and roof vents, Kittas et al. [4] established that the thermal effect is important when the ratio uo/∆Tio 0.5 <1, while Bot [21] put this limit at 0.3. In tests 1, 2, 3, 7, and 10, with uo/∆Tio 0.5 ≥1.5, hardly any air left the greenhouse through the roof vent, indicating that the wind effect clearly prevailed over the thermal effect, and the air left in quite a uniform fashion through the side vents ( Figure 4).

Evaluation of the Mean and Turbulent Ventilation Flows
The longitudinal component u x at the side vents is corrected according to Equation 6 to compensate for the change in wind speed during the tests. This did not prove necessary for the roof vent, as the measurements were taken continuously (Table 3).
By applying Equations (4) and (5), we determined the volumetric flow rates at the three greenhouse vents. Knowing the volumetric flow rate of the greenhouse, G M allows the ventilation rate R M to be calculated (Table 4 and Figure 5a).   (8)) and turbulent component G' (Equation (9) Test    The precision of the mean values of G M can be assessed by comparing them with the inflow and outflow volumes measured at the different ventilation surfaces. To verify the degree of satisfaction of the Mass Conservation Law in the greenhouse [1,28], the error in the calculation of the ventilation volumes has been estimated as follows: If we calculate the airflow without correcting the inside air velocity with the outside wind speed, we obtain a mean error for all the assays of E G = 17.2%, but on correcting it, the mean error of the assays is reduced to E G = 15.1%. Following a similar methodology in an Almería greenhouse with two side and two roof vents, Molina-Aiz et al. [19] obtained errors of E G between 3.0% and 37.0% (calculating airflows with u x * ) and between 8.4% and 70.5% (with u x ). In a multi-tunnel greenhouse with a roof vent, other authors obtained errors of 2.2% and 2.6% [1], and of 31.6% [28]. The accuracy of the method of calculating the ventilation volumetric flow rates will depend mainly on the stability of the wind conditions (intensity and direction) and on the influence of the thermal effect on greenhouse ventilation. When wind conditions are not stable or when the dominance of the thermal effect or the wind effect is not clear, a greater degree of error should be expected. Maximum error was recorded in the assays for which the ratio u o / ∆T io 0.5 was close to 1.5 (Figure 5b).
The greenhouse ventilation rate R M under conditions of natural ventilation (Table 4) (Figure 5a). According to this fit, for air velocity below this value, there would be no exchange of air with the outside, and so this cannot be considered a valid fit.

Determining the Discharge Coefficient of the Vents C d
The discharge coefficients (Table 5) have been obtained for each of the vents, following the methodology outlined by Molina-Aiz et al. [19]. The windward roof vent had the highest average discharge coefficient (C dWR = 0.161 ± 0.054), followed by the windward side vent (C dWS = 0.143 ± 0.023), and finally the leeward side vent (C dLS = 0.100 ± 0.062); the average discharge coefficient for both side vents was C dVS = 0.121 ± 0.042. The discharge coefficients of side vents C dVS and of roof vent C dWS can be expressed as a function of air velocity [C dVS = 0.028·u o + 0.028 (R 2 = 0.92 and p-value = 0.0000); C dWR = 0.036·u o + 0.040 (R 2 = 0.96 and p-value = 0.0000)] ( Figure 6).

Determining the Wind Effect Coefficient Cw
Once the discharge coefficient of the greenhouse vents, CdVS and CdWR (Table 5), and the average volumetric flow rate, GM, are known (Table 4), we can calculate the wind effect coefficients Cw for each of the experiments (Table 6).
This dimensionless wind effect coefficient expresses the relationship of the field of velocities measured on a reference level (the meteorological station) and on another level close to the greenhouse vents [2].

Determining the Wind Effect Coefficient C w
Once the discharge coefficient of the greenhouse vents, C dVS and C dWR (Table 5), and the average volumetric flow rate, G M , are known (Table 4), we can calculate the wind effect coefficients C w for each of the experiments (Table 6). Table 6. Wind effect coefficients C w obtained for the experimental greenhouse according to the original ventilation models M1, M2, and M3 (Equations (1), (2) and (3)) and the modified models M1 and M2 (Equations (13) and (14)). E v , coefficient of the effectiveness of the opening E v = C dM C w 0.5 based on model M3. This dimensionless wind effect coefficient expresses the relationship of the field of velocities measured on a reference level (the meteorological station) and on another level close to the greenhouse vents [2].

Original
Equations (1)-(3) have been used to obtain C w in accordance with the three models of ventilation, M1, M2, and M3, respectively, derived from Bernoulli's equation [4,7,10]. Ventilation models M1 and M2 consider that the wind and thermal effects complement each other. However, this does not always occur, and Li and Delsante [22] considered two situations: "fully assisting and fully opposing". According to these authors, the latter can occur when there are only two ventilation openings. Although the greenhouse in the present study has three ventilation openings, the ventilation models present a better fit considering the fully opposing scenario (subtracting the wind and thermal effect) than considering the fully assisting scenario (adding the wind and thermal effect). It is necessary to modify models M1 and M2 so as to indicate the relationship between the two effects: Equation (13) is similar to that proposed by Li and Delsante [22]. Table 6 shows the wind coefficients obtained with the original models M1 and M2 (Equations (1) and (2)) and with the modified models (Equations (13) and (14)). When we do not consider the appropriate relationship between the wind and thermal effects, very low wind effect coefficients are obtained, and even negative coefficients are found in some of the experiments for model M1 (Table 6).
To optimise the proposed ventilation models, it would be necessary to evaluate the real relationship between the two effects in each specific situation. The two simplest scenarios are those studied by Li and Delsante [22], "fully assisting and fully opposing", and we understand that the latter occurs in our greenhouse. In commercial greenhouses with more than three ventilation surfaces, an intermediate situation can occur, making it difficult to determine an appropriate model.
In Espinoza et al. [24], the flow patterns were determined in the east sector of the experimental greenhouse (Figure 1), opening 2 lateral vents combined with 2 or 3 roof vents (making a total of 4 or 5 vents). In the windward roof vent, opposition of the wind effect and thermal effect was observed, but in the leeward roof vent, these two effects were added. By combining windward and leeward roof vents in the greenhouse, an intermittent situation between "fully assisting and fully opposing" would take place. Table 6 and Figure 7 illustrate how the coefficient C w obtained by model M3 (which ignores the thermal effect) maintains relatively similar values for the different tests (between 0.134 and 0.218). This means that a mean value of C w can be considered for model M3 to predict the ventilation rate according to wind velocity, as has been done in previous research works [10,12,13,15,[17][18][19]. However, the same does not apply to models M1 and M2, in which the thermal effect is taken into account in the equations (and with a wind effect opposition). For model M1, the coefficient C w , obtained by applying Equation (13), takes values of between 0.111 and 0.433; while for model M2, applying Equation (14), it takes values of between 0.184 and 0.824 (Figure 7). This variation in C w when applying models M1 and M2 was also observed in an Almería-type greenhouse with two side vents and two roll-up roof vents [38]. The values of C w for models M1 and M2 have been seen to vary according to the ratio u o / ∆T io 0.5 (Figure 7).
Low values of the ratio u o / ∆T io 0.5 would indicate that the thermal effect has a major bearing on natural ventilation of the greenhouse. In this case, since the thermal and wind effects are opposing, on applying models M1 and M2, higher values of C w are obtained, possibly to compensate for the strong opposition of the two effects. High values of the ratio u o / ∆T io 0.5 would indicate that the thermal effect has less bearing on natural ventilation and that the wind effect predominates. In this case, by applying models M1 and M2, the C w values obtained are close to those obtained by model M3 (Figure 7), which ignores the thermal effect. It has been observed that the C w coefficient that corresponds to each test can be calculated  Low values of the ratio uo/∆Tio 0.5 would indicate that the thermal effect has a major bearing on natural ventilation of the greenhouse. In this case, since the thermal and wind effects are opposing, on applying models M1 and M2, higher values of Cw are obtained, possibly to compensate for the strong opposition of the two effects.
High values of the ratio uo/∆Tio 0.5 would indicate that the thermal effect has less bearing on natural ventilation and that the wind effect predominates. In this case, by applying models M1 and M2, the Cw values obtained are close to those obtained by model M3 (Figure 7), which ignores the thermal effect. It has been observed that the Cw coefficient that corresponds to each test can be calculated according to the ratio uo/∆Tio 0.5 or δ [Cw,M1 = exp(−2.693 + 1.160/δ) (R 2 = 0.94 and p-value = 0.0000); Cw,M2 = exp(−2.128+1.264/δ) (R 2 = 0.98 and p-value = 0.0000)] (Figure 7). Kittas et al. [2] observed how this coefficient increased at low air velocities in a tunnel greenhouse with continuous side openings.
Although models M1, M2, and M3 can be used in their current form to predict the greenhouse ventilation rate, as is explained below, more work is required to ascertain the dependence of Cw on the ratio uo/∆Tio 0.5 for models M1 and M2.
Using the original models M1 and M2, the mean values of Cw are 0.013 ± 0.079 (M1) and 0.014 ± 0.018 (M2) ( Table 6), which are much lower than those obtained by other authors (Table 7). On applying the modified ventilation models, these values are 0.186 ± 0.098 (M1) and 0.358 ± 0.192 (M2), which are closer to the values found by other researchers in multi-span type greenhouses and are considerably higher than those found by Molina-Aiz et al. [19] in an Almería-type greenhouse ( Table  7). Table 6. Wind effect coefficients Cw obtained for the experimental greenhouse according to the original ventilation models M1, M2, and M3 (Equations (1), (2) and (3)) and the modified models M1 and M2 (Equations (13) and (14)). Ev, coefficient of the effectiveness of the opening Ev = CdM Cw 0.5 based on model M3.  Although models M1, M2, and M3 can be used in their current form to predict the greenhouse ventilation rate, as is explained below, more work is required to ascertain the dependence of C w on the ratio u o / ∆T io 0.5 for models M1 and M2.

Cw
Using the original models M1 and M2, the mean values of C w are 0.013 ± 0.079 (M1) and 0.014 ± 0.018 (M2) ( Table 6), which are much lower than those obtained by other authors (Table 7). On applying the modified ventilation models, these values are 0.186 ± 0.098 (M1) and 0.358 ± 0.192 (M2), which are closer to the values found by other researchers in multi-span type greenhouses and are considerably higher than those found by Molina-Aiz et al. [19] in an Almería-type greenhouse ( Table 7).
The use of sonic anemometry has allowed us to observe the negative interaction between wind and thermal effects in the experimental greenhouse. Other techniques, such as the tracer gas method, allow correct quantification of the ventilation rates in the greenhouse, but they provide no information on the characteristics of the airflow [1], nor do they allow us to determine the way in which the two effects interact. Had the tracer gas method been used, applying the original ventilation methods would give much lower wind effect coefficients than those obtained taking into account the correct relationship between the wind and thermal effects. The wind effect coefficients obtained by other authors using the tracer gas method (Table 7) may be affected by a negative interaction between the two effects.
Since M3 does not consider the thermal effect, in general, lower wind effect coefficients are obtained (Table 6) than with models M1 and M2. Model M3 includes the thermal effect in the wind effect coefficients, and as the relationship between both effects is negative, lower wind effect coefficients are obtained than when the two effects are treated separately. Model M3 assumes that the volumetric caudal of air in a greenhouse is proportional to: the surface area of the vents S V , the air velocity u o , and the coefficient of effectiveness of the openings E v = C d ·C w 0.5 . The latter is frequently used in the literature to characterise the wind effect in the ventilation of greenhouses [19]. The present study has used a different expression in model M3 to that used by Molina-Aiz et al. [19]; rather than the coefficient C d and the total vent surface area S V , this study opts for the coefficients corresponding to the side vents C dVS , to the roof vent C dWR , and their corresponding surfaces. In our case, the coefficient E v of the vents has been obtained by calculating C d as the weighted average value between C dVS and C dWR , according to the surface of each vent. The coefficients of effectiveness obtained in the present work with M3 (Table 6), in some tests, are lower than those obtained by other authors in greenhouses with side and roof vents, and with insect-proof screens (Table 7), possibly due to the negative interaction between the wind and thermal effects in our greenhouse. ϕ, porosity of the insect-proof screens; S A , greenhouse base surface area; ue, air velocity. a , wind effect coefficient Cw calculated using C d = 0.65 (the discharge coefficient of the openings was not indicated by these authors). * Method: Gas (N 2 O tracer gas), H 2 O (Mass balance on water vapour), 3D (trisonic anemometer), Om (omnidirectional hot-ball anemometer).

Fitting the Semi-Empirical Models to Experimental Data
To assess the fit of the three models with the experimental values, the flow rates were obtained for each experiment applying the three models. These ventilation flow rates have been estimated by applying the ventilation models in four ways ( Figure 8): (i) first, using the mean wind coefficient C w (0.186 for M1, 0.358 for M2 and 0.177 for M3) and the particular discharge coefficients, C dVS and C dWR , for each test (Table 5); (ii) second, using the mean wind coefficient C w (idem) and the mean discharge coefficients C dVS and C dWR (equal to 0.121 and 0.161, respectively); (iii) third, using the mean wind coefficient C w (idem) and the fits C dVS = 0.028·u o + 0.028 and C dWR = 0.036·u o + 0.040; (iv) fourth, also using the fits C w,M1 = exp(−2.693 + 1.160/δ) and C w,M2 = exp(−2.128 + 1.264/δ).
Based solely on the determining coefficient R 2 , in the former case, the three models show a very good fit with the experimental data (with R 2 = 0.992 for M1 and M3, R 2 = 0.991 for M2).
In the second case, model M3 best fits the experimental data with R 2 = 0.952, followed by M1 with R 2 = 0.925, and finally M2 with R 2 = 0.919. Applying the models with the fits C dVS (u o ) and C dWR (u o ) improves the fit of models M1 and M2 (R 2 = 0.986) and of M3 (R 2 = 0.987). Applying the models with the fits C w,M1 (δ), C w,M2 (δ), C dVS (u o ), and C dWR (u o ), the values of R 2 go down slightly for models M1 and M2 (R 2 = 0.973 for M1 and R 2 = 0.980 for M2). In Molina-Aiz et al. [19], the models that obtain a better fit to the experimental data are Model M1 (R 2 = 0.984), considering the resulting pressure distribution as the sum of the pressure fields due to stack and wind effects, and the most simplified M3 (R 2 = 0.985). In the case of Molina-Aiz et al. [19], the relationship between wind and thermal effects was "fully assisting", with two side vents (one to windward and one to leeward) and two roll-up roof vents; the air entered through the side vents to exit through the roof vents.
In the second case, model M3 best fits the experimental data with R 2 = 0.952, followed by M1 with R 2 = 0.925, and finally M2 with R 2 = 0.919. Applying the models with the fits CdVS(uo) and CdWR(uo) improves the fit of models M1 and M2 (R 2 = 0.986) and of M3 (R 2 = 0.987). Applying the models with Analysis based solely on the coefficient R 2 may give the impression that the behaviour of the three models is very good, given the high R 2 values obtained in all cases. On comparing the simulated and observed G M values one by one (Figure 8), the simulated values are seen to vary considerably from the observed ones in certain cases. Table 8 presents the different statistics obtained for the three models, M1, M2, and M3, and the four different applications of the coefficients C w , C dVS , and C dWR (Applications i, ii, iii and iv in Table 8). With model M3, higher values of G M are obtained than those observed experimentally, irrespective of how the coefficients are applied (Figure 8), always obtaining values of PBIAS<0 (model overestimation bias). Models M1 and M2 present values of PBIAS<0 (model overestimation bias) for the first three applications of the coefficients (i, ii, and iii in Table 8) and values of PBIAS>0 (model underestimation bias) on applying the models using the fits C w,M1 (δ), C w,M2 (δ), C dVS (u o ), and C dWR (u o ). For the remaining statistics obtained, no great differences are observed between the different applications of the coefficients C w , C dVS , and C dWR (Applications i, ii, iii and iv in Table 8). We consider the best application of the models uses the fits obtained for C dVS and C dWR as a function of u o (for models M1, M2, and M3), and the fits obtained for C w as a function of the ratio u o / ∆T io 0.5 (for models M1 and M2). In this way, the values of the coefficients applied are closest to the values obtained experimentally.
In summary, model M3 considers a linear relationship between the wind velocity and the ventilation volumetric flow rates G M , a statistical device which improves the quality of the fit, and which therefore must be used and interpreted with great care [39]. It has also been observed that the ventilation rate fits better a second order polynomial than a linear relationship with air velocity. We therefore recommend using the modified model M1, since it includes both the wind and thermal effects and it is the one with the greatest physical foundation, as it combines the pressure fields rather than the flow rates [4,7]. For model M1, the lowest values of RMSD and RSR are obtained when applying the model using the fits C w,M1 (δ), C w,M2 (δ), C dVS (u o ), and C dWR (u o ) (Application iv in Table 8), rather than Applications i, ii, and iii. Moreover, before applying the model, it is essential to determine whether the interaction between the wind and thermal effects is positive or negative. Table 8. Statistical parameters obtained for the average volumetric flow rate for the greenhouse G M,sim (simulated with the different models M1, M2 and M3) with the experimentally observed values G M,obs . R 2 , determining coefficient; RMSD, root mean squared error; MD, bias; NSE, the Nash-Sutcliffe efficiency; PBIAS, percent bias; RSR, RMSD-observations standard deviation ratio. Applying the models with mean values of C w and the particular discharge coefficients C dVS and C dWR for each test (Application i), with mean values of C w , C dVS and C dWR (Application ii); with mean values of C w and the fits C dVS (u o ) and C dWR (u o ) (Application iii); with fits C w,M1 (δ), C w,M2 (δ), C dVS (u o ) and C dWR (u o ) (Application iv). In view of the results obtained in this work, we recommend to experimentally determine the coefficients C d and C w coefficients before applying the semi-empirical models M1, M2, or M3 in a greenhouse. Special care must be taken when they are not determined experimentally and when the coefficients C d and C w coefficients of the bibliography have to be selected. The coefficient C d depends on the geometry of the vent and on the characteristics of the insect-proof screen installed in the vent [19]; a coefficient obtained for a vent with similar characteristics must be selected. But as noted, the coefficient C d varies with air velocity through the vents, and this should also be taken into account. In the case of the coefficient C w , it has been observed that in the current form of the M1 and M2 models, this coefficient varies depending on the interaction of the thermal effect and thermal effect (ratio u o / ∆T io 0.5 ), which complicates the task of selecting coefficient C w of the bibliography that adapts well to the greenhouse where you want to apply these semi-empirical models.

Application i
3.3.4. Estimating the Contribution of the Thermal Effect in the Natural Ventilation of the Greenhouse By applying models M1 and M2, calculating the free component of the flux induced by buoyancy forces G T and the flux induced by wind forces G w , we estimated the contribution of the thermal effect to the natural ventilation of the greenhouse (using fits C w,M1 (δ), C w,M2 (δ), C dVS (u o ), and C dWR (u o )).
The contribution of the thermal effect in all the experiments (Table 9) is sufficient for it to be considered in any study of natural ventilation of greenhouses in regions with a similar climate to the Mediterranean one, in which the temperature gradients between inside and outside the greenhouse are considerable. Indeed, in test 1, with a wind velocity of 6.86 m s −1 and ratio u o / ∆T io 0.5 = 2.98, the thermal effect would be equivalent to 36% or 27% of the wind effect applying models M1 and M2, respectively ( Table 9). The limits established by other authors should be revised. For instance, the value of 2 m s −1 , above which the thermal effect may be ignored, as suggested by Boulard and Baille [7], Kittas et al. [39] and Fatnassi et al. [13], does not fit the results of the present study. Table 9. Ventilation flow rates generated by the wind effect (G w ) and the thermal effect (G T ) estimated by applying ventilation models M1 and M2 using the fits C wM2 (δ), C dVS (u o ), and C dWR (u o ).

Model M1 Model M2
Test

Estimating the Reduction in the Ventilation Rate Caused by the Insect-Proof Screens
The reduction in greenhouse ventilation due to the insect-proof screens may be considered proportional to the reduction in air velocity through the vents [40]. The reduction in the ventilation rate caused by the insect-proof screens (ϕ = 0.39) can be estimated by applying the ventilation models using the discharge coefficients of the openings calculated without the screens (C dVS = 0.657, C dWR = 0.712); determined following the methodology outlined by Molina-Aiz et al. [19]. Applying the mean value of C w (equal at 0.177) for model M3, and the fits C w,M1 (δ) and C w,M2 (δ) for models M1 and M2, the insect-proof screens (ϕ = 0.39) have been estimated to cause a mean reduction of 70% (M1), 61% (M2), and 85% (M3). Natural ventilation is extremely important for optimal plant growth during the summer in Mediterranean countries. For most of the year, a good system of natural ventilation will allow growers to maintain suitable microclimate conditions inside the greenhouse for the crops [41]. However, Valera et al.
[3] point out that the 14.4% average ventilation surface (S V /S A ) in Almería's greenhouses is well below the minimum recommended value of 30% [42,43]. Indeed, it is also considerably lower than the 25% recommended in the Andalusian Regulations on Integrated Production of Protected Horticultural Produce [44]. Von Zabeltitz [45] recommends between 18% and 25% for the Mediterranean Basin, while Kittas et al. [41] quote a value of between 15% and 30%. It should be noted that all of the above recommendations refer to greenhouses that are not equipped with insect-proof screens over the vent openings. The value of ventilation rate R M required for a temperature rise of 5 • C between outside-inside air temperatures vary from 0.02 to 0.09 m 3 s -1 m -2 (15 to 65 h -1 for a greenhouse with an average height of 5 m) depending on the solar radiation and the crop transpiration [46], with an optimum value of 45-60 h -1 [47].
In this work, the experiments were carried out for a ventilation surface, S V /S A , of 0.112. By applying the ventilation models as follows: with the mean value of C w = 0.177 for model M3; the fits C w,M1 (δ) and C w,M2 (δ) for models M1 and M2; and the fits C dVS (u o ) and C dWR (u o ) for the three models; it has been estimated that for the vent arrangement in the experimental greenhouse, a ventilation surface S V /S A of 0.62 (M1), 0.54 (M2), and 0.36 (M3) would be necessary to reach a value of 45 h −1 (under conditions of southwest wind and u o = 5 m s −1 ).

Conclusions
Sonic anemometry has allowed us to identify the entrance and exit vents of greenhouse ventilation air, thus allowing us to establish natural ventilation flow patterns for greenhouses. According to most of the works reviewed in the literature, the ventilation rate of a greenhouse with natural ventilation fits a linear relationship with air velocity. In the present case, it was observed that this relationship fits better a second order polynomial, R M = 0.37 u o 2 + 0.03 u o + 0.75 (R 2 = 0.99 and p-value = 0.0000).
Opening the roof vent to windward, one side vent to leeward, and the other side vents to windward (this last vent obstructed by another greenhouse) brings about the opposition of the thermal and wind effects in natural ventilation of the greenhouse, as it drives the air entering through it downwards, opposing natural convection due to the thermal effect. The modification of the ventilation models obtained from Bernouilli's equation in order to add or subtract the airflows due to wind and thermal effects improves their precision.
G(M1) = G 2 T ± G 2 w and G(M2) = |G T ± G w | A linear relationship has been obtained, which allows us to estimate the discharge coefficient of the side vents (C dVS ) and roof vent (C dWR ) as a function of wind velocity [C dVS = 0.028·u o + 0.028 (R 2 = 0.92); C dWR = 0.036·u o + 0.040 (R 2 = 0.96)]. The wind effect coefficient determined by applying models M1 and M2 has not proved constant for the different tests carried out, but rather it varies depending on the ratio u o / ∆T io 0.5 or δ [C w,M1 = exp(−2.693 + 1.160/δ) (R 2 = 0.94); C w,M2 = exp(−2.128+1.264/δ) (R 2 = 0.98]. The fits for C d and C w allow us to estimate the ventilation flow rate from the outside and inside temperatures and from the wind velocity using the ventilation models. The contribution of the thermal effect to the natural ventilation of the greenhouse has been quantified based on models M1 and M2. The opposition of the thermal effect accounted for at least 36% or 27% of the wind effect (applying models M1 and M2, respectively). The thermal effect must therefore always be taken into account in studies on the natural ventilation of greenhouses. Applying models M1, M2, and M3, it has been estimated that the presence of insect-proof screens on the vents (ϕ = 39.0%) can reduce the ventilation flow rate of the greenhouse by 70%, 61%, and 85%, respectively. Likewise, by applying these models, it has been ascertained that the surface area of ventilation with respect to the surface area of the soil, S V /S A , should be 0. 36