Simulation of Water Vapor and Near Infrared Radiation to Predict Vapor Pressure Deﬁcit in a Greenhouse Using CFD

: Vapor pressure deﬁcit (VPD) can be used as an indicator to schedule greenhouse irrigation. VPD can be estimated as a function of relative humidity (RH) and temperature (T). The objective of this work was to analyze spatial variation in VPD as an indicator of water stress inﬂuenced by concentration of water vapor and intensity of near infrared (NIR). The study was carried out in an empty three-span sawtooth greenhouse with natural ventilation under the local climate in Montecillo, Mexico; these ﬁndings established a base value to analyze greenhouse ﬁeld conditions prior to the inﬂuence from a crop. The experimental phase consisted of recording data (3 February 2019–24 February 2019) on temperature, humidity, solar radiation, and wind speed, which were used for developing a model in computational ﬂuid dynamics (CFD). Then, this model was used to estimate VPD, considering changes in mass fraction of water vapor and the intensity of NIR. Scenarios with 50, 70, and 90% external RH were evaluated. It was found that without a crop, temperature was not affected by the variation in the mass fraction of water vapor and the intensity of NIR in the simulated scenarios, each of which generated a thermal gradient within the range of 4 ◦ C. When considering the scenario of 90% external RH, we found the best VPD range along the greenhouse (2–3 kPa) that would be a favorable ﬁeld condition for crops. Differences between VPD with and without a crop can be used to estimate the water quantity needs for crop growth based on the climate variables examined in this study, where higher VPD values require more water for irrigation.


Introduction
In greenhouses, the ventilation system favors the distribution of temperature inside. With natural ventilation, air exchange sometimes becomes insufficient due to outside wind speeds equal to or less than 1 m s −1 , causing an increase in temperature and reducing the concentration of CO 2 and relative humidity [% (RH)] [1][2][3]. In 2017, Cemek et al. [4] concluded that ventilation efficiency should be further analyzed with the aim of establishing the best environmental conditions favoring plant growth parameters. Wind speed and direction are factors impacting the environmental comfort of the internal climate of a greenhouse, specifically in terms of the exchange of heat and mass between the inside and outside [5][6][7]. The ventilation system of a greenhouse may eliminate excess heat and maintain acceptable CO 2 and water vapor levels for crops [8,9]. Strategies of natural ventilation are diversified according to the season of the year; for instance, during winter, natural ventilation allows the elimination of excess humidity, while in summer it sometimes provides insufficient cooling, requiring the use of auxiliary mechanical ventilation systems [10]. greenhouse [48][49][50]. VPD is the difference between the amount of water in air as vapor and the amount of moisture the air can retain when it is saturated with water, so it can be considered a useful way to measure the microclimate of a greenhouse; the value of VPD can be estimated through the equation of Rosenberg et al. [51], VPD = 0.61078exp [(17.269 · T)/(T + 237.3)] · [1 + RH/100], where T is the air temperature.
The objective of this work was to develop and evaluate a computational model to estimate VPD values as an indicator of water stress, due to the relationship between water quantity or vapor concentration (RH) and NIR in the greenhouse environment. The idea of simulating an empty greenhouse assumed that the exchange of mass and energy between the exterior and interior of the greenhouse only occurs through natural ventilation, and there is no generation of water vapor inside the greenhouse. Predicting VPD values without cultivation facilitates an easier understanding of the requirements to achieve recommended ranges.

Description of Experimental Site
The greenhouse used in this study is located at Campus Montecillo, Colegio de Postgraduados in the State of Mexico (19 • 27 45.5" N, 98 • 54 12.2" W) with an altitude of 2239 masl, oriented from north to south on a plain where monthly mean temperatures range from 5 to 25 • C. It has an area of 1050 m 2 ; it is a three-span sawtooth type covered with translucent polyethylene on the floor and sides, polycarbonate on the roof, and antiinsect mesh on three side windows and three roof windows. The greenhouse relies on a natural ventilation system with the roof and side windows fully open.
For data acquisition, an eddy covariance system was installed at a distance of 2.4 m, windward, equipped with a sonic anemometer (3D, Campbell Scientific Ltd., Antony, France), an EC150 gas analyzer (Campbell Scientific Inc., Logan, UT, USA), a net radiometer (NR-LITE2 Kipp and Zone BV, Delft, The Netherlands), an LI-200 pyranometer (LI-COR Inc., Lincoln, NE, USA), an LI-190R Quantum sensor (Li-Cor Inc., Lincoln, NE, USA), a temperature and humidity sensor (HMP155A, Campbell Scientific, Logan, UT, USA), and four sensors for recording soil temperature (TCAV thermocouple, Campbell Scientific Inc., Logan, UT, USA). Seven temperature sensors (DS18B20, Maxim Integrated, San Jose, CA, USA) and three humidity sensors (DHT22, Zhengxinyuan Electronics Co., Ltd., Shenzhen, China) were installed in distributed locations in the three spans at a height of 3 m ( Figure 1). Processes 2021, 9,1587 3 of 14 cladding material [44], time of day, season of the year [45,46], the type of ventilation, and its configuration [47]. However, scarce or no information is available to analyze the effect of the intensity of NIR and the concentration of water vapor on the air temperature of a greenhouse [48][49][50]. VPD is the difference between the amount of water in air as vapor and the amount of moisture the air can retain when it is saturated with water, so it can be considered a useful way to measure the microclimate of a greenhouse; the value of VPD can be estimated through the equation of Rosenberg et al. [51], VPD = 0.61078exp [(17.269 · T)/(T + 237. 3)] · [1 + RH/100], where T is the air temperature. The objective of this work was to develop and evaluate a computational model to estimate VPD values as an indicator of water stress, due to the relationship between water quantity or vapor concentration (RH) and NIR in the greenhouse environment. The idea of simulating an empty greenhouse assumed that the exchange of mass and energy between the exterior and interior of the greenhouse only occurs through natural ventilation, and there is no generation of water vapor inside the greenhouse. Predicting VPD values without cultivation facilitates an easier understanding of the requirements to achieve recommended ranges.

Description of Experimental Site
The greenhouse used in this study is located at Campus Montecillo, Colegio de Postgraduados in the State of Mexico (19°27′45.5" N, 98°54′12.2" W) with an altitude of 2239 masl, oriented from north to south on a plain where monthly mean temperatures range from 5 to 25 °C. It has an area of 1050 m 2 ; it is a three-span sawtooth type covered with translucent polyethylene on the floor and sides, polycarbonate on the roof, and anti-insect mesh on three side windows and three roof windows. The greenhouse relies on a natural ventilation system with the roof and side windows fully open.
For data acquisition, an eddy covariance system was installed at a distance of 2.4 m, windward, equipped with a sonic anemometer (3D, Campbell Scientific Ltd., Antony, France), an EC150 gas analyzer (Campbell Scientific Inc., Logan, UT, USA), a net radiometer (NR-LITE2 Kipp and Zone BV, Delft, The Netherlands), an LI-200 pyranometer (LI-COR Inc., Lincoln, NE, USA), an LI-190R Quantum sensor (Li-Cor Inc., Lincoln, NE, USA), a temperature and humidity sensor (HMP155A, Campbell Scientific, Logan, UT, USA), and four sensors for recording soil temperature (TCAV thermocouple, Campbell Scientific Inc., Logan, UT, USA). Seven temperature sensors (DS18B20, Maxim Integrated, San Jose, CA, USA) and three humidity sensors (DHT22, Zhengxinyuan Electronics Co., Ltd., Shenzhen, China) were installed in distributed locations in the three spans at a height of 3 m ( Figure 1).  Data logging and storage were performed every 5 s for temperature and humidity, 100 ms for the anemometer, and 20 s for the net radiometer. The monitoring period was during the month of February 2019, the month when there commonly is a RH gradient that oscillates between 30-100% inside the greenhouse between the day and night periods, respectively. During data collection, there was no crop, thus the humidity inside the greenhouse came from outside. Information storage was carried out with a data logger system (CR3000 Micrologger, Campbell Scientific, Inc., Logan, UT, USA) for the eddy covariance system, and an Arduino Data Logger Shield for the DS18B20 and DHT22 sensors. Outside the greenhouse there was a meteorological station (Vantage Pro2 Plus, Davis Instruments, Hayward, CA, USA) used to collect data on temperature, humidity, solar radiation, and wind direction and speed-measurements that were used as the initial boundary conditions of the computational model. The stored information was downloaded and processed in Microsoft Excel to determine the input environmental conditions of the computational model.

Computational Model
The computational model was developed and simulated in ANSYS ® Fluent ® . The geometry was developed in ANSYS ® Workbench using the SpaceClaim tool. The mesh of the model domain was built in Meshing. The properties of the greenhouse materials were incorporated into the model following Tables 1 and 2 [8,19,44,52]. Mesh sensitivity analysis of the model was performed to obtain reliable predictions in the same way as in [53]. Mesh was built with structured elements, using an orthogonal quality of 0.98 and distortion of 1.935 × 10 −0.02 . Table 1. Properties of cladding materials.

Material
Density  Table 2. Optical properties of materials to thermal radiation.

Materials Emissivity (
Data logging and storage were performed every 5 s for temperature and humidity, 100 ms for the anemometer, and 20 s for the net radiometer. The monitoring period was during the month of February 2019, the month when there commonly is a RH gradient that oscillates between 30-100% inside the greenhouse between the day and night periods, respectively. During data collection, there was no crop, thus the humidity inside the greenhouse came from outside. Information storage was carried out with a data logger system (CR3000 Micrologger, Campbell Scientific, Inc., Logan, UT, USA) for the eddy covariance system, and an Arduino Data Logger Shield for the DS18B20 and DHT22 sensors. Outside the greenhouse there was a meteorological station (Vantage Pro2 Plus, Davis Instruments, Hayward, CA, USA) used to collect data on temperature, humidity, solar radiation, and wind direction and speed-measurements that were used as the initial boundary conditions of the computational model. The stored information was downloaded and processed in Microsoft Excel to determine the input environmental conditions of the computational model.

Computational Model
The computational model was developed and simulated in ANSYS ® Fluent ® . The geometry was developed in ANSYS ® Workbench using the SpaceClaim tool. The mesh of the model domain was built in Meshing. The properties of the greenhouse materials were incorporated into the model following Tables 1 and 2 [8,19,44,52]. Mesh sensitivity analysis of the model was performed to obtain reliable predictions in the same way as in [53]. Mesh was built with structured elements, using an orthogonal quality of 0.98 and distortion of 1.935 × 10 −0.02 .

Material
Density The conditions of the model and of the variables applied to the air flow to solve the transport equations discretized in algebraic equations and calculated by numerical methods are defined in Table 3. For the species model, the mass fraction of water vapor was used to calculate the relative humidity distribution. The conditions of the model and of the variables applied to the air flow to solve the transport equations discretized in algebraic equations and calculated by numerical methods are defined in Table 3. For the species model, the mass fraction of water vapor was used to calculate the relative humidity distribution.

Fundamental Equations of Modeled Flow
The general transport equation that can be solved numerically by CFD to explicitly calculate the velocity and temperature field associated with a flow in three dimensions at steady state can be described in the following terms [19,51,54]: where, Φ represents the concentration of the dimensionless term, that is, momentum, mass (air and water vapor mass fraction), and energy. U, V, and W are the velocity vector components in the three dimensions (m s −1 ), Γ is the diffusive coefficient (kg m −1 s −1 ), and S Φ is the source term. Turbulence inside the greenhouse was considered using the k-e model [18,55] that gave rise to two supplementary equations and species. The choice of the model resulted from the comparison of different turbulence models in greenhouses [38].
For energy, the source term is a function of NIR absorption. Therefore, the equation can be described in general terms as follows [38,41]: where, c p is the specific heat of air (J kg −1 K −1 ), ρ is the air density (kg m −3 ), T is the air temperature (K or • C), u i is the velocity component (m s −1 ), λ f is the thermal conductivity of air (W m −1 K −1 ), µ t is the turbulent viscosity of air (kg m −1 s −1 ), Pr t is the turbulent Prandtl number, and S r is the source term.

Radiative Modeling Equation
The discrete ordinate model (DO) solves the radiative transfer equation (RTE) for a finite number of angles, each associated with a → s vector direction in the global Cartesian system [17,19,33,38,56]. The luminance distribution I λ → r , → s was obtained through the solution of the RTE; the activation of the RTE submodel was used to couple the DO model and divide the near infrared NIR spectral region (0.76-1.1 µm) [33,37,57].
where, I λ is the radiation intensity for the wavelength (W m −2 sr −1 ), φ is the scattering phase function, Ω' is the solid angle, along the radiation propagation direction, → s I is the direction of the scattering vector, σ s is the scattering coefficient (m −1 ), a λ is the spectral absorption coefficient (m −1 ), I bλ is the monochromatic luminance of a black body (W m −3 sr −1 ), and n is the refractive index.
The radiation parameter was considered as constant in the band length. The net radiation flux absorbed by a control volume was then defined by the expression described by Nebbali et al. [38]: where S r is the source term of the energy Equation (1), allowing a coupling between the energy and radiation transfer equations.

Evaluation of the Computational Model
The computational model was evaluated by means of a statistical analysis of the average of the temperature data, RH inside the greenhouse, and the point data of wind speed and NIR, simulated, and experimental. An analysis of variance (ANOVA) with a significance level of 0.05 was used to parametrically contrast the difference of the variables. The data of the experimental models used for the evaluation were obtained between 01:30 p.m. and 03:30 p.m. (3 February 2019); this time was considered to be the solar middle day, where a change of position occurs at the zenith.

Simulation Scenarios
The simulated scenarios were based on a hypothesis regarding variation in the mass fraction of water vapor concentration in the air derived from temporal climatic variations. The change in concentration of the water vapor mass fraction in the computational model was performed on the basis of the RH range recorded in the air from February 3 to 24. The computational model was focused on the highest NIR intensity during the day, which was used for the simulation scenarios. In the computational model, only the exterior RH concentration was artificially changed to three scenarios for simulating: (a) 50, (b) 70, and (c) 90%, which were equivalent to a water vapor mass fraction of 0.013, 0.018, and 0.023 (g/g), respectively. The change in percent of RH to grams of water vapor per grams of air allowed estimating the effect of moisture concentration on the thermal gradient, assuming that water vapor is a greenhouse gas and impacts the temperature and VPD variables of a greenhouse. Subsequently, two more simulations were performed to modify the NIR intensity (b1 = 50 W m −2 and b2 = 100 W m −2 ), for which scenario b was considered because it represents the RH of a typical day in February for the site of study. The simulation of the scenario incorporating cultivation was performed using only the porous medium approach, where there is an effect on pressure drop, which the Forchheimer equation can express. The environmental analysis of the simulations inside the greenhouse was carried out at a height between 1.5 and 2 m, a height considered representative for growing tomatoes, and at which environmental changes affect the quality of the crop.

Validation of the Computational Model
To evaluate the feasibility of using the computational model, we used a statistical indicator that measures the correlation between the simulated temperature, humidity, wind speed, and NIR values and those recorded by the sensors on 9 February 2019. The variation of the experimental temperature data from 2:00 p.m. to 2:30 p.m. shows that the values were concentrated, indicating that the average value of the temperature recorded by the sensors was very probable, and the recorded data were reliable (Figure 2).  Table 4 shows that the mean measurements of temperature, wind speed, RH, and NIR, registered and simulated, are values close to each other. To measure this correlation, an ANOVA statistical analysis was performed with a significance level of 0.05, considering that if the significance level was less than the probability value p (level of significance < p value) there was no significant difference between the registered and simulated variables. (Table 5). Therefore, the evaluated model could be used in the simulations.

Analysis of the Results of Simulated Scenarios
The environmental analysis of the simulated scenarios was performed on the basis of the concentration of water vapor in the air, temperature, NIR, and vapor pressure deficit (VPD)-indicators used to manage the greenhouse environment and maintain suitable conditions for growing tomatoes within it. The analysis of the influence of concentration of the water vapor on the temperature was limited to measurements at the center of each span; in contrast, simulated data were computed as averages along the greenhouse, away from the predominant air inlet in the greenhouse. In each simulated scenario, a dataset of 100 data points was obtained longitudinally at a height of 1.5 m and was used for the  Table 4 shows that the mean measurements of temperature, wind speed, RH, and NIR, registered and simulated, are values close to each other. To measure this correlation, an ANOVA statistical analysis was performed with a significance level of 0.05, considering that if the significance level was less than the probability value p (level of significance < p value) there was no significant difference between the registered and simulated variables. (Table 5). Therefore, the evaluated model could be used in the simulations.

Analysis of the Results of Simulated Scenarios
The environmental analysis of the simulated scenarios was performed on the basis of the concentration of water vapor in the air, temperature, NIR, and vapor pressure deficit (VPD)-indicators used to manage the greenhouse environment and maintain suitable conditions for growing tomatoes within it. The analysis of the influence of concentration of the water vapor on the temperature was limited to measurements at the center of each span; in contrast, simulated data were computed as averages along the greenhouse, away from the predominant air inlet in the greenhouse. In each simulated scenario, a dataset of 100 data points was obtained longitudinally at a height of 1.5 m and was used for the environmental analysis shown in Figure 3. These distances were selected for the analysis Processes 2021, 9, 1587 8 of 14 because the central span was characterized by greater thermal and humidity gradients in the air between the windward and leeward zones. The results in Figure 3 show a decrease in temperature in relation to the increase in water vapor concentration in the simulated scenarios. The thermal gradient range in the greenhouse (length) was around 2 • C in the crop zone (1.5 m height), being remarkable in the first half of the greenhouse. In the area near the outlet, the thermal gradient range was negligible (0.15 • C) between scenarios a, b, and c. environmental analysis shown in Figure 3. These distances were selected for the analysis because the central span was characterized by greater thermal and humidity gradients in the air between the windward and leeward zones. The results in Figure 3 show a decrease in temperature in relation to the increase in water vapor concentration in the simulated scenarios. The thermal gradient range in the greenhouse (length) was around 2 °C in the crop zone (1.5 m height), being remarkable in the first half of the greenhouse. In the area near the outlet, the thermal gradient range was negligible (0.15 °C) between scenarios a, b, and c.
On the other hand, for the VPD along the greenhouse length, there was a 2 kPa difference between scenarios, and this remained constant. This finding was derived from the influence on temperature (T) and relative humidity (RH) from external local climate conditions. In addition, for crop simulations, enhanced environmental conditions favoring crop growth were observed inside the greenhouse. We concluded that the most significant change occurred in the first half of the greenhouse. In the area near the outlet, the increase in water vapor concentration did not influence the increase in temperature. To calculate the VPD in the greenhouse, the equation was programmed with a userdefined function (UDF) in CFD. The VPD value was estimated as the difference between the saturation vapor pressure value and the actual vapor pressure of the air (VPD = PVsat − Pvair), a method recommended by Howell and Dusek, and FAO [58,59]. The greenhouse analysis was performed in all the simulated scenarios as an average of 10 data points uniformly distributed within a height range of 1.9-2 m at the inlet, center, and outlet of the central span. The results in Table 6 show the spatial variation in the water vapor concentration and thermal gradient at three points of observation along the length of the greenhouse.
The temperature increase was generated from the absorption of NIR by water vapor; however, when the mass fraction of water vapor remains constant in the air at the center and outlet areas of the greenhouse, the RH gradient is minimal equal to 0.78%. Therefore, the thermal effect generated by the energy exchange between water vapor and NIR is negligible. In all scenarios, there is a reduction between the outside RH and that of the greenhouse with a gradient that ranges between 45 and 50% caused by the low speed of the outside wind (0.64 m s −1 ) and the anti-insect mesh in the windows. On the other hand, for the VPD along the greenhouse length, there was a 2 kPa difference between scenarios, and this remained constant. This finding was derived from the influence on temperature (T) and relative humidity (RH) from external local climate conditions. In addition, for crop simulations, enhanced environmental conditions favoring crop growth were observed inside the greenhouse. We concluded that the most significant change occurred in the first half of the greenhouse. In the area near the outlet, the increase in water vapor concentration did not influence the increase in temperature.
To calculate the VPD in the greenhouse, the equation was programmed with a user-defined function (UDF) in CFD. The VPD value was estimated as the difference between the saturation vapor pressure value and the actual vapor pressure of the air (VPD = PVsat − Pvair), a method recommended by Howell and Dusek, and FAO [58,59]. The greenhouse analysis was performed in all the simulated scenarios as an average of 10 data points uniformly distributed within a height range of 1.9-2 m at the inlet, center, and outlet of the central span. The results in Table 6 show the spatial variation in the water vapor concentration and thermal gradient at three points of observation along the length of the greenhouse. Table 6. Average of ten data of the variables: T, RH, and VPD simulated at a height from 1.9 to 2 m at the inlet, center, and outlet area of the central span of the greenhouse. The temperature increase was generated from the absorption of NIR by water vapor; however, when the mass fraction of water vapor remains constant in the air at the center and outlet areas of the greenhouse, the RH gradient is minimal equal to 0.78%. Therefore, the thermal effect generated by the energy exchange between water vapor and NIR is negligible. In all scenarios, there is a reduction between the outside RH and that of the greenhouse with a gradient that ranges between 45 and 50% caused by the low speed of the outside wind (0.64 m s −1 ) and the anti-insect mesh in the windows.

Position Scenario T ( • C) RH (%) VPD (kPa)
The increase in temperature inside the greenhouse, given the increase in water vapor concentration in the air, was not significant, ranging from 0.15 to 0.28 • C (Table 6). For this reason, it was inferred that the temperature remained constant under the simulated conditions and was not influenced by the variation in the mass fraction of water vapor.
The simulation of scenario b was carried out with a change in the intensity emitted by the NIR (50 and 100 W m −2 ) to analyze the effect it had on temperature. The results of the simulations (Table 7) show that when there was a change in the intensity of NIR, the highest thermal gradient (0.23 • C) occurred at the inlet and outlet areas of span 2 of the greenhouse, where the wind speeds were below 0.01 m s −1 and air movement was generated by the convective effect caused by temperature differences. This was also demonstrated by Sun et al., and Zeroual et al. [3,7], where outside wind speeds close to 1 m s −1 were insufficient to renew the air inside the greenhouse, affecting the exchange of heat and mass between the outside and the inside. Derived from the change in the concentration of the mass fraction of water vapor and the NIR intensity in the simulated scenarios, the highest thermal gradient occurred in the inlet area (0.29 • C), where the RH of the greenhouse was higher due to the air exchange between outside and inside. Therefore, it could be assumed that the generated temperature gradient was negligible, as in most of the greenhouse it oscillated between 0.08-0.16 • C. Figure 4 shows the distribution of RH and VPD inside the greenhouse in the simulated scenarios at a height of 1.5 m. The results indicate that scenario c presented better environmental conditions for tomato cultivation with lower VPD (kPa), temperature T ( • C), and higher RH (%) inside the greenhouse. Nevertheless, there was a 53% reduction in RH between the outside and inside of the greenhouse, causing this scenario to be outside the desirable optimum for tomato (50-60%) [60]. The RH deficit inside the greenhouse caused an increase in the VPD, with values over 1.3 kPa for all simulated scenarios, values outside the range for tomato (0.3-1.25 kPa) according to Grange et al.,Körner et al. and Shtienberg [15,61,62]. The VPD values in Figure 4 show that under the conditions of the simulated scenarios, the crop inside the greenhouse would undergo a high rate of transpiration, causing the closure of its stomata and impeding its growth [14].
Processes 2021, 9,1587 10 of 14 outside the range for tomato (0.3-1.25 kPa) according to Grange et al.,Körner et al. and Shtienberg [15,61,62]. The VPD values in Figure 4 show that under the conditions of the simulated scenarios, the crop inside the greenhouse would undergo a high rate of transpiration, causing the closure of its stomata and impeding its growth [14].  Table 6.
Scenario b was simulated with the mean monthly wind speed (1.59 m s −1 ) during the month of February in Montecillo, Mexico according to data from the National Meteorological Service (SMN). The purpose of the simulation under these conditions was to predict the effects on environmental conditions that can occur when wind speed increases. These changes were simulated, and the results in Figure 5 show enhanced environmental conditions for tomato crop cultivation compared to the scenario b of Figure 4 with decreases in temperature and VPD that oscillated between 4-6 °C and 0.9-1.63 kPa, respectively. Likewise, it presents an increase in RH of between 8.3-15.4% throughout the greenhouse. However, these generated changes were not enough to maintain a RH between 50-60%, a temperature below 30 °C, and a VPD between 0.3-1.25 kPa throughout the greenhouse.  Table 6.
Scenario b was simulated with the mean monthly wind speed (1.59 m s −1 ) during the month of February in Montecillo, Mexico according to data from the National Meteorological Service (SMN). The purpose of the simulation under these conditions was to predict the effects on environmental conditions that can occur when wind speed increases. These changes were simulated, and the results in Figure 5 show enhanced environmental conditions for tomato crop cultivation compared to the scenario b of Figure 4 with decreases in temperature and VPD that oscillated between 4-6 • C and 0.9-1.63 kPa, respectively. Likewise, it presents an increase in RH of between 8.3-15.4% throughout the greenhouse. However, these generated changes were not enough to maintain a RH between 50-60%, a temperature below 30 • C, and a VPD between 0.3-1.25 kPa throughout the greenhouse.
The results of the simulations carried out in this work also confirm what was described by Chu, Lan, Tasi, Wu, and Yang; and Tong et al. [46,63], where lengths greater than six times the height of the greenhouse present problems in the ventilation rate due to the obstruction generated by the anti-insect mesh and the length of the greenhouse. Therefore, the use of fans and extractors is recommended to avoid stagnant areas in the center of greenhouse buildings when natural ventilation is insufficient to improve the ventilation rate [9,10,64]. The results of the simulations carried out in this work also confirm what was described by Chu, Lan, Tasi, Wu, and Yang; and Tong et al. [46,63], where lengths greater than six times the height of the greenhouse present problems in the ventilation rate due to the obstruction generated by the anti-insect mesh and the length of the greenhouse. Therefore, the use of fans and extractors is recommended to avoid stagnant areas in the center of greenhouse buildings when natural ventilation is insufficient to improve the ventilation rate [9,10,64].

Conclusions
The greenhouse air inlet area was the location with the highest energy exchange in all of the simulated scenarios, with a temperature increase of 0.29 °C generated by the absorption of NIR and by the exchange of air between the outside and inside. The greenhouse scenarios did not present favorable environmental conditions for tomato cultivation. The temperature exceeded 30 °C with an average thermal gradient between the outside and inside of up to 12 °C, causing an RH deficit and an increase in VPD to values higher than 2 kPa, which consequently would cause excess perspiration in tomato and an expected stomatal closure. The temperature gradient between the outside and inside of the greenhouse was the result of the low wind speed (0.64 m s −1 ) and the anti-insect mesh, causing air movement within the greenhouse to be generated by the difference between soil and air temperatures. Environmental conditions improved when the scenario was simulated with the monthly mean wind speed for the month of February (1.59 m s −1 ) in Texcoco; however, this was not enough to maintain favorable environmental conditions. From this, it is inferred that natural ventilation is insufficient to reduce the thermal gradient and increase the RH to acceptable values for tomato. Therefore, the use of fans and extractors is recommended to increase the exchange of heat and mass between the outside and inside and help to distribute the RH and temperature evenly in the stagnant area of the greenhouse (in the center of the spans).

Conclusions
The greenhouse air inlet area was the location with the highest energy exchange in all of the simulated scenarios, with a temperature increase of 0.29 • C generated by the absorption of NIR and by the exchange of air between the outside and inside. The greenhouse scenarios did not present favorable environmental conditions for tomato cultivation. The temperature exceeded 30 • C with an average thermal gradient between the outside and inside of up to 12 • C, causing an RH deficit and an increase in VPD to values higher than 2 kPa, which consequently would cause excess perspiration in tomato and an expected stomatal closure. The temperature gradient between the outside and inside of the greenhouse was the result of the low wind speed (0.64 m s −1 ) and the anti-insect mesh, causing air movement within the greenhouse to be generated by the difference between soil and air temperatures. Environmental conditions improved when the scenario was simulated with the monthly mean wind speed for the month of February (1.59 m s −1 ) in Texcoco; however, this was not enough to maintain favorable environmental conditions. From this, it is inferred that natural ventilation is insufficient to reduce the thermal gradient and increase the RH to acceptable values for tomato. Therefore, the use of fans and extractors is recommended to increase the exchange of heat and mass between the outside and inside and help to distribute the RH and temperature evenly in the stagnant area of the greenhouse (in the center of the spans).