Microclimatic Evaluation of Five Types of Colombian Greenhouses Using Geostatistical Techniques

In Colombia, the second-largest exporter of cut flowers worldwide and one of the South American countries with the largest area of crops under cover, passive or naturally ventilated greenhouses predominate. Locally, there are several types of greenhouses that differ in architecture, size, height, shape of roof and ventilation surfaces, of which many characteristics of the microclimate generated in their interior environment are unknown. This generates productive limitations that in some way may be limiting the yield, quality and health of the final products harvested; in addition, Colombian producers do not have the ability to monitor the microclimate of their farms, much less to correlate microclimate data with data on crop production and yield. Therefore, there is a need for the Colombian grower to know the most relevant microclimate characteristics generated in the main greenhouses used locally. The objective of this work was to carry out a microclimatic characterization of the five most used types of greenhouses in Colombia. The main results allowed determining that in these structures, there are conditions of high humidity and low vapor pressure for several hours of the day, which affects the physiological processes of growth and development of the plants. It was also identified that for each type of greenhouse, depending on the level of radiation, there is a significant microclimatic heterogeneity that may be the cause of the heterogeneity in plant growth, which is a common characteristic observed by the technical cultivation personnel. Therefore, it can be concluded that it is urgent to propose microclimatic optimization strategies to help ensure the sustainability of the most important production systems in the country.


Introduction
In crops grown under protective environment structures, microclimate management is undoubtedly one of the most relevant factors affecting crop production in terms of both quality and quantity [1,2]. Worldwide, there are more than 5.2 million hectares covered with different types of protective agriculture structures [3,4]. These structures present different technological levels, from a low technological level, such as that found in the small tunnels or structures typically used in tropical countries of Latin America and Asia [5], to an extremely high level, such as that found in the high-tech greenhouses found mainly in northern European countries, where it is possible to fully manage the variables that affect crop growth and development [6].
In Colombia, there is an important area dedicated to cultivation under cover, mainly for the production of cut flowers for the international market; locally, different types of cover structures are built, mostly with plastic covers and made of low-cost materials, A microclimatic optimization of Colombian greenhouse structures aiming at increasing the yields of the main crops under cover must be preceded by a characterization of the microclimate behavior of the most-used greenhouse models at the local level. This information will allow us to establish when and how much we should manage and control variables such as temperature and air humidity. Therefore, the objective of this work was to characterize the microclimate behavior of the five most-used greenhouses in Colombia. The methodology consisted of measuring variables such as temperature, humidity, vapor pressure deficit (VPD), internal solar radiation level and CO 2 concentration inside the greenhouse on an hourly scale. Finally, using a grid of 40 temperature and humidity sensors and implementing geostatistical techniques, a spatial variability analysis was carried out to determine the climatic heterogeneity present in this type of greenhouse.

Description and Location of the Greenhouses
For this study, 5 types of greenhouses built in the savanna of Bogota were selected, among which there were already well-developed ornamental crops of Rosa and Carnation. Each of the analyzed greenhouse models presents variations in their shape, size, surface and arrangement of the ventilation areas. The main characteristics of each greenhouse are described below.

Traditional Greenhouse (TG)
The traditional Colombian greenhouse (TG) is one of the most widely used models in the country for flower production, as well as for vegetable and aromatic production [8]. The greenhouse evaluated was located in the municipality of Guasca, department of Cundinamarca (4 • 50 57.60 N, 73 • 52 58.19 N and altitude: 2664 m.s.a.l).
The structure of the TG greenhouse is made of wood, and its cover is made of commercial polyethylene of 0.2 mm thickness, the management of the microclimate is performed through the management of the ventilation areas arranged on the 4 sides of the greenhouse and in the region of the roof, the technical specifications and some architectural details of TG can be consulted in Table 1. Table 1. General characteristics of the TG greenhouse.

Greenhouse Schematic Description
Sensors 2022, 22, x FOR PEER REVIEW A microclimatic optimization of Colombian greenhouse structures aim ing the yields of the main crops under cover must be preceded by a charact microclimate behavior of the most-used greenhouse models at the local le mation will allow us to establish when and how much we should mana variables such as temperature and air humidity. Therefore, the objective o to characterize the microclimate behavior of the five most-used greenhous The methodology consisted of measuring variables such as temperature, h pressure deficit (VPD), internal solar radiation level and CO2 concentra greenhouse on an hourly scale. Finally, using a grid of 40 temperature and sors and implementing geostatistical techniques, a spatial variability analy out to determine the climatic heterogeneity present in this type of greenho

Description and Location of the Greenhouses
For this study, 5 types of greenhouses built in the savanna of Bogota among which there were already well-developed ornamental crops of Ro tion. Each of the analyzed greenhouse models presents variations in their s face and arrangement of the ventilation areas. The main characteristics of ea are described below.

Traditional Greenhouse (TG)
The traditional Colombian greenhouse (TG) is one of the most widel in the country for flower production, as well as for vegetable and aromatic The greenhouse evaluated was located in the municipality of Guasca, depa dinamarca (4°50′57.60′′ N, 73°52′58.19′′ N and altitude: 2664 m.s.a.l).
The structure of the TG greenhouse is made of wood, and its cover is mercial polyethylene of 0.2 mm thickness, the management of the micro formed through the management of the ventilation areas arranged on th greenhouse and in the region of the roof, the technical specifications and tural details of TG can be consulted in Table 1.  Another commonly used greenhouse is the multi-tunnel (MG) type; the model evaluated was located in the municipality of El Rosal, department of Cundinamarca (4 • 50 56.27 N, 74 • 14 16.72 N and altitude: 2592 m.s.a.l). This type of greenhouse is manufactured with a metallic structure, and the cover is made of commercial plastic of 0.2 mm thickness, the management of the microclimate is performed through the management of the ventilation areas arranged on the sides and a region of the roof, the main characteristics of MG can be reviewed in Table 2. Table 2. General characteristics of the MG greenhouse.

Greenhouse Schematic Description
Sensors 2022, 22, x FOR PEER REVIEW (4°50′56.27′′ N, 74°14′16.72′′ N and altitude: 2592 m.s.a.l). This type of greenh ufactured with a metallic structure, and the cover is made of commercial pla thickness, the management of the microclimate is performed through the m the ventilation areas arranged on the sides and a region of the roof, the mai tics of MG can be reviewed in Table 2.  To a lesser proportion, we find greenhouses of the so-called hanging ty type of structure is made of wood and with a commercial plastic cover of 0 ness. The structural form of HG presents a structural variation with respec tional TG greenhouse in order to obtain higher thermal inertia and a lower ing, based on a larger volume of air enclosed in the roof area and a small structural elements with respect to TG. In this study, two models of hanging were evaluated, HGT1, which had a larger covered area than HGT2; bot HGT2 had ventilation surfaces with a similar opening on all four sides of th and in the roof area. On the other hand, another marked difference betwee HGT2 is related to the arrangement and orientation of the pivot roof windo the case of HGT1, are all oriented in the same direction and arranged in the of each span, while in HGT2, the pivot windows have alternate orientations arrangements in each span. The greenhouses evaluated were located in the of Nemocón, department of Cundinamarca. (5°03′26.44′′ N, 73°55′47.34′′ N 2570 m.s.a.l) and some general characteristics are summarized in Tables  To a lesser proportion, we find greenhouses of the so-called hanging type (HG), this type of structure is made of wood and with a commercial plastic cover of 0.2 mm thickness. The structural form of HG presents a structural variation with respect to the traditional TG greenhouse in order to obtain higher thermal inertia and a lower level of shading, based on a larger volume of air enclosed in the roof area and a smaller number of structural elements with respect to TG. In this study, two models of hanging greenhouses were evaluated, HGT1, which had a larger covered area than HGT2; both HGT1 and HGT2 had ventilation surfaces with a similar opening on all four sides of the greenhouse and in the roof area. On the other hand, another marked difference between HGT1 and HGT2 is related to the arrangement and orientation of the pivot roof windows, which, in the case of HGT1, are all oriented in the same direction and arranged in the same region of each span, while in HGT2, the pivot windows have alternate orientations and different arrangements in each span. The greenhouses evaluated were located in the municipality of Nemocón, department of Cundinamarca. ( Finally, a spatial greenhouse (SG) was evaluated. This type of structure is built with concrete pillars, elements that function as columns, while the roof region is made of intertwined steel cables on which the plastic covers are fastened, and a ventilation area is generated. This type of greenhouse is equipped with ventilation surfaces on the 4 sides and in a region of the roof of each span. The model evaluated was located in the municipality of Soacha, department of Cundinamarca. (4 • 35 35.53 N, 74 • 12 32.70 N and altitude: 2546 m.s.a.l) and its main characteristics are shown in Table 5. Table 3. General characteristics of the HGT1 greenhouse.   Finally, a spatial greenhouse (SG) was evaluated. This type of structure is concrete pillars, elements that function as columns, while the roof region is mad twined steel cables on which the plastic covers are fastened, and a ventilation a erated. This type of greenhouse is equipped with ventilation surfaces on the 4 in a region of the roof of each span. The model evaluated was located in the mu of Soacha, department of Cundinamarca. (4°35′35.53′′ N, 74°12′32.70′′ N and alti m.s.a.l) and its main characteristics are shown in Table 5.  Table 4. General characteristics of the HGT2 greenhouse.   Finally, a spatial greenhouse (SG) was evaluated. This type of structure is concrete pillars, elements that function as columns, while the roof region is mad twined steel cables on which the plastic covers are fastened, and a ventilation a erated. This type of greenhouse is equipped with ventilation surfaces on the 4 in a region of the roof of each span. The model evaluated was located in the mu of Soacha, department of Cundinamarca. (4°35′35.53′′ N, 74°12′32.70′′ N and alti m.s.a.l) and its main characteristics are shown in Table 5.  Table 6 lists all of the meteorological measuring equipment used both in the outdoor environment and inside each of the greenhouses. Inside each of the greenhouses, it was decided to monitor the values of temperature, relative humidity, vapor pressure deficit (VPD), radiation level and CO 2 concentration level since these are the main microclimate variables that affect the growth and development of crops inside protected agriculture structures [24,25]. Generally, the presence of the translucent cover and the geometry of the greenhouse cause the temperature and humidity inside the greenhouse to increase, the radiation level to decrease, and the CO 2 level to fluctuate more in the time scale, so it is expected that these variables will affect the quality and quantity of the harvested product [26]. Finally, it is also important to mention that, in terms of microclimate and crop growth simulation, the variables studied in this research are the most predicted from any simulation methodology since they are directly related to agricultural production and energy and resource efficiency inside a greenhouse [27]. Table 5. General characteristics of the SG greenhouse.    Table 6 lists all of the meteorological measuring equipment used both in the outdoor environment and inside each of the greenhouses. Inside each of the greenhouses, it was decided to monitor the values of temperature, relative humidity, vapor pressure deficit (VPD), radiation level and CO2 concentration level since these are the main microclimate variables that affect the growth and development of crops inside protected agriculture structures [24,25]. Generally, the presence of the translucent cover and the geometry of the greenhouse cause the temperature and humidity inside the greenhouse to increase, the radiation level to decrease, and the CO2 level to fluctuate more in the time scale, so it is expected that these variables will affect the quality and quantity of the harvested product [26]. Finally, it is also important to mention that, in terms of microclimate and crop growth simulation, the variables studied in this research are the most predicted from any simulation methodology since they are directly related to agricultural production and energy and resource efficiency inside a greenhouse [27]. The pyranometer and the CO2 sensor inside the greenhouse were placed in the central zone of each structure, while the 40 thermocouples were evenly distributed according to the covered area of each greenhouse; the 40 sensors were included inside a ventilated capsule that was placed 1.6 m above the ground level just above the zone where the plant canopy grows. Each of the greenhouses was monitored for 40 calendar days (Table 7), and  The pyranometer and the CO 2 sensor inside the greenhouse were placed in the central zone of each structure, while the 40 thermocouples were evenly distributed according to the covered area of each greenhouse; the 40 sensors were included inside a ventilated capsule that was placed 1.6 m above the ground level just above the zone where the plant canopy grows. Each of the greenhouses was monitored for 40 calendar days (Table 7), and the frequency of data recording was established every ten minutes; these data recording periods are comparable since climatic conditions in Colombia do not vary greatly throughout the months of the year. The vapor pressure deficit variable was determined by means of psychrometric Equations (1)-(3) that relate air temperature and relative humidity, as described in the work developed by Huang et al. [28].

Greenhouse Schematic Description
where (e s ) is the saturation vapor pressure and (e a ) is the actual vapor pressure.

Geostatistical Analysis
The analysis of the collected data using geostatistical techniques was carried out in three phases, the first of which is known as exploratory analysis, the second phase as structural analysis and the third and last phase is the one that allows us to generate the spatial prediction curves of the variables of interest and is known as the prediction phase [29]. To perform the geostatistical analysis, the datasets were organized, grouped and averaged according to 5 intervals categorized according to the outdoor solar radiation level (ESRL). These radiation levels were determined from exploratory analyses of the data, and an attempt was made to group the data into radiation intervals representative of the study region, as performed by Bojacá et al. [12] (Table 8). In these types of structures, the spatial variability of variables such as temperature, relative humidity and VPD is dependent on the level of radiation; therefore, more specific analyses should focus on what happens in the daytime period, the period between 6 and 18 hours of the day (ERL 2-ERL 5), while at night, the spatial variability of the variables is less evident due to the lack of solar radiation in the nighttime period between 18 and 6 h of the following day (ERL 1).

Exploratory Analysis Phase
This is performed with the objective of verifying that the datasets comply with some geostatistical assumptions, such as stationarity and seasonality. Usually, the graphs are generated to inspect the stationarity of the study variables [30]. Additionally, it is verified if the datasets are seasonal; that is, if they present a periodic behavior, for it is assumed that the degree of association of the data depends on the spatial separation and not on the location of the data in space [31]. Finally, basic statistics, such as the Shapiro-Wilk test, were calculated to determine the normal distribution of the data, which, together with the initial graphical analysis, allow us to rule out the existence of any trend in space [12]. In this work, these calculations were performed on an hourly scale for temperature, relative humidity and vapor pressure deficit. On the other hand, radiation and CO 2 concentration data were also averaged both indoors and outdoors, although for these two variables, the geostatistical analysis is not complete since there is only one data recording point for each variable inside the greenhouses.

Structural Analysis Phase
In this phase, the empirical semivariograms associated with each dataset of each variable collected must be constructed, which is nothing more than the graphical relationship that exists between the spatial separation distances and the semivariance of the 3 variables studied [32]; this semivariogram was constructed by means of the equation where Mp d is the number of pairs at distance d, d is the increase, W(xi) is the experimental values and xi is the locations where the values are measured W(xi). Empirical semivariograms should be fitted to theoretical semivariogram models; in this case, we used and evaluated the fit for 5 theoretical models (Gaussian, circular, spherical, Cauchy and Exponential). These models have the following parameters: range (a), sill (c 1 ) and nugget (c 0 ) . The nugget represents both sampling error and spatial variability for the minimum distance interval, the sill represents the upper asymptote of the semivariogram and corresponds to the maximum variability explained by the semivariogram, and the range represents a measure of the maximum distance of influence [12]. The theoretical semivariogram model to be used in the prediction phase was selected according to the Bayesian information criterion (BIC) and the Akaike information criterion (AIC). The model with the lowest values for these criteria will be the model that directly presents the best fit to the experimental data.

Prediction Phase
The information from the two previous phases is used to predict the behavior of the variables of interest at unsampled points based on experimentally sampled values and the calculation of their spatial continuity structure [33]. Likewise, the predicted values are used to create the spatial distribution curves in a two-dimensional plane at a height of 1.6 m above ground level. The selected prediction method is ordinary kriging, which is considered to be an optimal linear predictor, unbiased and with minimum variance. The prediction is calculated with the following equation where λ i is the weight assigned to each value of the variables in the observed positions W(x 0 ) based on the parameters of the previously adjusted theoretical semivariogram.

Exploratory Analysis
For the analysis of each of the variables, the data were averaged on an hourly scale for the total number of sensors used both indoors and outdoors of the greenhouse and for the total number of days of the measurement.
3.1.1. Temperature Figure 1 shows the average temperature of all sensors on an hourly scale inside the greenhouse and the hourly average value of the outside sensor. In general terms, the typical behavior of this type of structure in tropical climates is observed for all the greenhouses evaluated, with the lowest temperature conditions occurring at night and the highest between noon and 15:00 hours, just after the moments of maximum external solar radiation [23]. The mean temperature (T m ) values for the 24 h inside the greenhouses were 15.1 ± 3.6, 14.0 ± 3.5, 15.9 ± 5.1, 15.7 ± 3.5 and 15.9 ± 3.1 • C for TG, MG, SG, HGT1 and HGT2, respectively, which are 2.0, 0.9, 1.6, 1.4 and 1.8 • C higher than the outside ambient conditions for each of these greenhouses.
On the other hand, the maximum mean values (T mh ) occur just between 11 and 14 h for TG, MG, SG, HGT1 and HGT2, where values of 21.4, 19.4, 24.3, 21.5 and 21.1 • C, respectively, were obtained. These T mh values, with the exception of SG, are not within the ranges established as optimal daily temperature values, which should range between 24 and 30 • C for the optimal growth and development of a great variety of ornamental and vegetable species cultivated in greenhouses [4]. For the case of the minimum mean values (Tml), these values were 11.3, 9.8, 10.1, 11.8 and 12.7 • C, representing thermal differentials between the indoor and outdoor environment of 1.1, 0.03, −0.7, 0.07 and 1.7 • C, for TG, MG, SG, HGT1 and HGT2, respectively. It is important to highlight that the values in all the greenhouses are below the minimum value of optimum temperature for species such as Rosa and Carnation, which, according to the literature, should range between 14 and 16 • C [34]. spectively, were obtained. These Tmh values, with the exception of SG, are not within the ranges established as optimal daily temperature values, which should range between 24 and 30 °C for the optimal growth and development of a great variety of ornamental and vegetable species cultivated in greenhouses [4]. For the case of the minimum mean values (Tml), these values were 11.3, 9.8, 10.1, 11.8 and 12.7 °C, representing thermal differentials between the indoor and outdoor environment of 1.1, 0.03, −0.7, 0.07 and 1.7 °C, for TG, MG, SG, HGT1 and HGT2, respectively. It is important to highlight that the values in all the greenhouses are below the minimum value of optimum temperature for species such as Rosa and Carnation, which, according to the literature, should range between 14 and 16 °C [34]. Finally, it is also important to mention that for MG and SG, we can observe, in Figure  1 and in the thermal differentials, values that show that in the interior of these greenhouses, in the hours of the night and the dawn of 20 and 5 hours, the phenomenon of thermal inversion occurs, a phenomenon where the greenhouse temperature is lower than the outside environment [35]. This phenomenon undoubtedly affects crop growth rates and is caused by the ineffectiveness of these greenhouses in retaining the loss of thermal radiation that occurs from the greenhouse floor to the outside environment, as well as the poor hermeticity of these local structures [36]. Therefore, for these conditions, it is imperative to establish climate control strategies to optimize the climate under nocturnal conditions in these types of greenhouses used in Colombia, as mentioned in previous studies [36,37].

Relative Humidity
In the case of relative humidity (RH), the hourly behavior can be observed for each greenhouse in Figure 2. This variable is undoubtedly the most difficult to control within Finally, it is also important to mention that for MG and SG, we can observe, in Figure 1 and in the thermal differentials, values that show that in the interior of these greenhouses, in the hours of the night and the dawn of 20 and 5 hours, the phenomenon of thermal inversion occurs, a phenomenon where the greenhouse temperature is lower than the outside environment [35]. This phenomenon undoubtedly affects crop growth rates and is caused by the ineffectiveness of these greenhouses in retaining the loss of thermal radiation that occurs from the greenhouse floor to the outside environment, as well as the poor hermeticity of these local structures [36]. Therefore, for these conditions, it is imperative to establish climate control strategies to optimize the climate under nocturnal conditions in these types of greenhouses used in Colombia, as mentioned in previous studies [36,37].

Relative Humidity
In the case of relative humidity (RH), the hourly behavior can be observed for each greenhouse in Figure 2. This variable is undoubtedly the most difficult to control within the greenhouse microclimate management since its value will depend on the psychrometric condition of the air and the mass contributions that occur via transpiration and water evaporation phenomena inside the greenhouse [4,38]. It can be observed that the humidity values in the greenhouses always show values higher than 67% for the average hourly value; the maximum humidity conditions inside the greenhouses occur during the night hours and generally exceed 90%, while the minimum values occur between 11 and 14 h and range between 67% and 75%.
The RH value for the 24 h of the day was 93 ± 6%, 87 ± 7%, 84 ± 11%, 88 ± 11% and 87 ± 12% inside TG, MG, SG, HGT1 and HGT2, respectively, values differing by 6%, 2%, 6%, -1.4% and -0.5% with respect to the outside environment for each greenhouse. Therefore, it is evident that for TG, MG and SG, the RH conditions are higher than the outdoor environment, and the RH values inside the greenhouses are above the RH values recommended for crops under cover, which range between 60 and 80% [39]. These conditions of high humidity inside the greenhouses generate a favorable environment for ornamental or horticultural species to be infested and affected by pathogens that generate fungal diseases, which have a high capacity to cause damage to plants. Therefore, productive losses and environmental damage associated with the application of chemical control substances for diseases, such as gray mold and downy mildew, caused by Botrytis cinérea and Chromista Peronospora sparsa, occur [40,41].
hours and range between 67% and 75%.
The RH value for the 24 h of the day was 93 ± 6%, 87 ± 7%, 84 ± 11%, 88 ± 11% and 87 ± 12% inside TG, MG, SG, HGT1 and HGT2, respectively, values differing by 6%, 2%, 6%, -1.4% and -0.5% with respect to the outside environment for each greenhouse. Therefore, it is evident that for TG, MG and SG, the RH conditions are higher than the outdoor environment, and the RH values inside the greenhouses are above the RH values recommended for crops under cover, which range between 60 and 80% [39]. These conditions of high humidity inside the greenhouses generate a favorable environment for ornamental or horticultural species to be infested and affected by pathogens that generate fungal diseases, which have a high capacity to cause damage to plants. Therefore, productive losses and environmental damage associated with the application of chemical control substances for diseases, such as gray mold and downy mildew, caused by Botrytis cinérea and Chromista Peronospora sparsa, occur [40,41].

Vapor Pressure Deficit (VPD)
One of the variables of interest for decision-making in climate and irrigation management inside the greenhouse is the VPD, a variable that relates the difference between the water vapor pressure at saturation and the actual water vapor pressure at a pressure from the greenhouse temperature [39]. VPD is a variable dependent on the temperature and RH of the greenhouse environment; its importance for climate management lies in the fact that this variable is the driving force of plant transpiration [4]. Its behavior in each of the greenhouses can be observed in Figure 3; it is observed that the maximum values of VPD occur during daytime hours, precisely in the periods where the temperature increases and the RH decreases, while the lowest values of VPD occur during the nighttime

Vapor Pressure Deficit (VPD)
One of the variables of interest for decision-making in climate and irrigation management inside the greenhouse is the VPD, a variable that relates the difference between the water vapor pressure at saturation and the actual water vapor pressure at a pressure from the greenhouse temperature [39]. VPD is a variable dependent on the temperature and RH of the greenhouse environment; its importance for climate management lies in the fact that this variable is the driving force of plant transpiration [4]. Its behavior in each of the greenhouses can be observed in Figure 3; it is observed that the maximum values of VPD occur during daytime hours, precisely in the periods where the temperature increases and the RH decreases, while the lowest values of VPD occur during the nighttime period, where the cold and saturated environments predominate in each of the greenhouses.
The mean values of VPD for the 24 h of the day were 0.12 ± 0.12, 0.24 ± 0.17, 0.36 ± 0.34, 0.23 ± 0.24 and 0.24 ± 0.26 kPa, in TG, MG, SG, HGT1 and HGT2, respectively, values that are lower by 0.08, 0.02, and 0.03 kPa for TG, MG and SG with respect to the outside environment, while HGT1 and HGT2 values are higher with respect to the outside environment value by 0.04 and 0.02 kPa, behaviors that are directly related to the thermal and humidity conditions present in each greenhouse. The maximum mean VPD values in the greenhouses were 0.34, 0.5, 0.98, 0.70 and 0.71 kPa for TG, MG, SG, HGT1 and HGT2, respectively. These VPD values, except for TG, were during the 11 and 15 hours within the recommended ranges for plant growth, where VPD should range between 0.45 and 1.25 kPa, although the ideal should be at 0.9 kPa [42]. These values present in greenhouses undoubtedly limit the productive potential of crops, as VPD affects photosynthetic activity and plant nutrition, causing deficiencies due to limited sap flow movement [43,44].
In the case of the minimum mean VPD values, these values were 0.001, 0.06, 0.05, 0.01 and 0.008 kPa for TG, MG, SG, HGT1 and HGT2, respectively. These values, so close to 0, generally occur during the entire night period and are quite inadequate for plant growth, generating a favorable environment for the appearance and survival of fungal diseases [45,46]. Therefore, it is necessary to analyze these situations and look for alternatives to improve the microclimate conditions related to VPD present in the greenhouses evaluated.

Solar Radiation
Solar radiation is the energy engine of the greenhouse and is the variable responsible for generating the greenhouse effect inside the structures and its respective increase in temperature; also, part of this energy is captured by the soil to be liberated at night [47]. The hourly behavior of solar radiation can be seen in Figure 4. In general, it is observed that the radiation value increases from 6 to 13 hours and then decreases again between 13 and 18 hours, this being a typical behavior of the value during the daytime period in tropical regions. In the case of the minimum mean VPD values, these values were 0.001, 0.06, 0.05, 0.01 and 0.008 kPa for TG, MG, SG, HGT1 and HGT2, respectively. These values, so close to 0, generally occur during the entire night period and are quite inadequate for plant growth, generating a favorable environment for the appearance and survival of fungal diseases [45,46]. Therefore, it is necessary to analyze these situations and look for alternatives to improve the microclimate conditions related to VPD present in the greenhouses evaluated.

Solar Radiation
Solar radiation is the energy engine of the greenhouse and is the variable responsible for generating the greenhouse effect inside the structures and its respective increase in temperature; also, part of this energy is captured by the soil to be liberated at night [47]. The hourly behavior of solar radiation can be seen in Figure 4. In general, it is observed that the radiation value increases from 6 to 13 hours and then decreases again between 13 and 18 hours, this being a typical behavior of the value during the daytime period in tropical regions.
The mean radiation values obtained during daylight hours were 180 ± 142, 235 ± 159, 224 ± 177, 145 ± 115 and 144 ± 108 Wm −2 for TG, MG, SG, HGT1 and HGT2, respectively. These values represent a transmission percentage of 75%, 87%, 76%, 74% and 71% for each of the greenhouses mentioned above. It should be noted that the roof that generates the highest radiation transmission is the tunnel type of the MG greenhouse; unlike the flat roofs of the other greenhouses, this percentage of transmission depends on the roof material, which was similar for all greenhouses in material and age of use, while the geometric shape of the roofs was different and this generates a different percentage of transmission [48].

Carbon Dioxide (CO 2 ) Concentration
Finally, the concentration of carbon dioxide (CO 2 ) is another microclimatic variable that directly affects the photosynthetic assimilation of plants and, therefore, is a variable that must be monitored inside greenhouses [49]. In the case of this study, it should be mentioned that the greenhouses studied do not have carbon enrichment systems, and there is no local information on the behavior of this variable in these types of greenhouses. Figure 5 shows the hourly trend of the CO 2 concentration in each of the greenhouses, the mean values of CO 2 concentration in the greenhouses were 381. 4  The mean radiation values obtained during daylight hours were 180 ± 142, 235 ± 159, 224 ± 177, 145 ± 115 and 144 ± 108 Wm -2 for TG, MG, SG, HGT1 and HGT2, respectively. These values represent a transmission percentage of 75%, 87%, 76%, 74% and 71% for each of the greenhouses mentioned above. It should be noted that the roof that generates the highest radiation transmission is the tunnel type of the MG greenhouse; unlike the flat roofs of the other greenhouses, this percentage of transmission depends on the roof material, which was similar for all greenhouses in material and age of use, while the geometric shape of the roofs was different and this generates a different percentage of transmission [48].

Carbon Dioxide (CO2) Concentration
Finally, the concentration of carbon dioxide (CO2) is another microclimatic variable that directly affects the photosynthetic assimilation of plants and, therefore, is a variable that must be monitored inside greenhouses [49]. In the case of this study, it should be mentioned that the greenhouses studied do not have carbon enrichment systems, and there is no local information on the behavior of this variable in these types of greenhouses. This is due to the fact that plants are consuming CO2 in their photosynthesis process at a faster rate than the CO2 renewal rate of the air via natural ventilation; therefore, this factor may be indicating that the ventilation rates in greenhouses are not adequate [50]. The relationship between the level of CO2 concentration in the greenhouse air and its ventilation rates has been demonstrated [51], and in this study, it can be verified since it was found that the greenhouses with the largest ventilation surfaces (HGT1 and HGT2) are the ones with CO2 curves in the diurnal period that are very close for both outside and inside environments.

Testing of Geostatistical Assumptions
In order to rule out any kind of trend in the datasets collected and, additionally, to verify that these same datasets fit a normal distribution, the graphs shown in Figure 6 are constructed. In this case, trends can be observed for data collected in MG and ESRL 5 for temperature ( Figure 6A), in ESRL 3 for relative humidity ( Figure 6B) and in ESRL 4 for VPD ( Figure 6C). These curves show that for both the width and length of the greenhouse, most of the values of the three variables analyzed are randomly distributed without any trend; the temperature ranges between 17.8 and 22.0 °C, the relative humidity between This is due to the fact that plants are consuming CO 2 in their photosynthesis process at a faster rate than the CO 2 renewal rate of the air via natural ventilation; therefore, this factor may be indicating that the ventilation rates in greenhouses are not adequate [50]. The relationship between the level of CO 2 concentration in the greenhouse air and its ventilation rates has been demonstrated [51], and in this study, it can be verified since it was found that the greenhouses with the largest ventilation surfaces (HGT1 and HGT2) are the ones with CO 2 curves in the diurnal period that are very close for both outside and inside environments.

Testing of Geostatistical Assumptions
In order to rule out any kind of trend in the datasets collected and, additionally, to verify that these same datasets fit a normal distribution, the graphs shown in Figure 6 are constructed. In this case, trends can be observed for data collected in MG and ESRL 5 for temperature ( Figure 6A), in ESRL 3 for relative humidity ( Figure 6B) and in ESRL 4 for VPD ( Figure 6C). These curves show that for both the width and length of the greenhouse, most of the values of the three variables analyzed are randomly distributed without any trend; the temperature ranges between 17.8 and 22.0 • C, the relative humidity between 70% and 96% and the VPD between 0.01 and 0.86 kPa.
It is also possible to infer, in principle, that from the bell-shaped histograms of the three variables, there is an adjustment to the normality of the dataset. This was confirmed by calculating the Shapiro-Wilk test, where it is established that to guarantee the normality of the data, the p-value must be greater than the defined alpha (0.05), and for these datasets, p-values of 0.48, 0.39 and 0.33 were obtained for temperature, relative humidity and VPD, respectively. This normality of the datasets is always desirable in this type of work since it guarantees adequate predictions of the study variables by means of the kriging method [52]. This same procedure was performed for all 75 datasets collected in this study, obtaining the same satisfactory results as for the 3 datasets discussed above. Table 9 shows the values obtained for 75 datasets corresponding to the analysis of the fit of the experimental data to theoretical semivariogram models. In general terms, it can be observed that as ESRL 5 radiation levels increase, there is a greater variation in temperature, humidity and VPD inside the greenhouses; on the other hand, when there is a lower level of ESRL 1 radiation, these microclimatic variations inside the greenhouses are less significant. For the case of temperature, it was found that for 19 datasets, the best fitting theoretical semivariogram model (lowest AIC and BIC value) was the circular model, followed by the spherical and Gaussian models each with 3 datasets, respectively. In the case of relative humidity, in 20 datasets, the best theoretical semivariogram model was the circular model, and in 5 datasets, the Gaussian model; finally, for VPD, 24 datasets showed a better fit for the circular semivariogram model and only 1 dataset under the Gaussian model.

Structural Analysis
The values of the nugget or nugget effect (c 0 ) for temperature ranged from a maximum value of 0.37 in MG and ESRL 5 to a minimum value of 0.0 obtained in 11 of the 25 datasets, in the case of relative humidity; these values of c 0 ranged from a maximum value of 6.9, obtained in SG and ESRL 5, to a minimum value of 0.0 obtained in 15 of the 25 datasets. Finally, in VPD, the value of the nugget effect ranged between a maximum of 1.88, obtained in HGT2 and ESRL 2, and a minimum of 0.0, obtained in 21 of the 25 datasets for this variable, all these values of c 0 can be considered low, which allows us to deduce that the spatial dependence due to the error of the experimental measurement is not relevant [12]. For the sill parameter (c 1 ), the values for temperature ranged from a minimum value of 0.09, obtained in TG and ESRL 1, to a maximum of 3.1, in SG and ESRL 5. In the case of relative humidity, the values of c 1 ranged from a minimum value of 0.0, obtained in HGT1 and ESRL 1, to a maximum value of 114.1 in SG and ESRL 5. For VPD, the minimum value was 0.0, obtained in several datasets, and the maximum value was 0.22, obtained in SG and ESRL 5. In general terms, it was found for the three variables that as the level of radiation increases (ESRL 1 to ESRL 5), the values of c 1 increase, which allows us to affirm that as the level of radiation increases, a greater spatial variation of temperature, relative humidity and VPD is found inside the greenhouses [41].  Table 9 shows the values obtained for 75 datasets corresponding to the analysis of the fit of the experimental data to theoretical semivariogram models. In general terms, it can be observed that as ESRL 5 radiation levels increase, there is a greater variation in temperature, humidity and VPD inside the greenhouses; on the other hand, when there is a lower level of ESRL 1 radiation, these microclimatic variations inside the greenhouses are less significant. For the case of temperature, it was found that for 19 datasets, the best fitting theoretical semivariogram model (lowest AIC and BIC value) was the circular  Table 9. Model parameters of the selected theoretical semivariograms and their quality of fit measures (Akaike Information Criteria; AIC, and Bayesian; BIC) for the temperature, relative humidity and VPD datasets of the five greenhouses evaluated.

Temperatura ( • C)
Humedad Relativa (%) VPD (kPa) For the ratio (c 0 /c 1 ), it was found that the values for temperature ranged from minimum values of 0.0% obtained in several datasets to a maximum of 18.6% obtained in MG and ESRL 4; for humidity the minimum value was 0.0%, obtained in several datasets, and a maximum value of 22.1% obtained in TG and ESRL 5, while for VPD, the minimum value was 0.0% in most of the datasets and a maximum of 4.5% in SG and ESRL 4. These values, being less than 25%, can be considered adequate and are indicative of a strong spatial dependence between the data [53]. These values also allow us to conclude that the number of sensors used in each of the greenhouses is sufficient for the prediction of temperature, relative humidity and VPD values in the unsampled points of each greenhouse [12].

Prediction of the Spatial Behavior of the Microclimate
The prediction of the variables analyzed for non-sampled points inside the greenhouse was carried out using the ordinary kriging method. These estimations are a function of the results of the adjustment of the theoretical semivariograms, so their validity is given by the quality of the adjustment of the theoretical semivariograms.
For the analysis of spatial behavior, we selected the contrasting ESRLs (ESRL 1 and ESRL 5), which will allow us to analyze the scenarios of greater and lesser microclimatic homogeneity. The other scenarios (ESRL 2, ESRL 3 and ESRL 4) are included graphically in Appendix A of this document.

Spatial variability of microclimate in MG
For the case of the aggregated data for the night hours (ESRL 1), it was found that for the MG greenhouse, the value of Tm was 11.1 ± 1.1 • C (Figure 7). It was observed that inside the greenhouse, the temperature in the largest amount of covered area ranges between 10.4 and 11.3 • C, while the maximum magnitude of the thermal differential inside the greenhouse (∆ TG ) was 4.1 • C between small areas located in the central region and some points of the front sides. In the case of RH (Figure 7), the mean value was 93 ± 1.6 • C, with a very important area between values of 94% and 97%, while the maximum humidity differential (∆ RHG ) was 11% between some perimeter points of MG. For the VPD variable, the mean value was 0.09 ± 0.03 kPa, while the VPD differential (∆ VPDG ) was 0.16 kPa between areas located in the center and the perimeter of MG.

Prediction of the Spatial Behavior of the Microclimate
The prediction of the variables analyzed for non-sampled points inside the greenhouse was carried out using the ordinary kriging method. These estimations are a function of the results of the adjustment of the theoretical semivariograms, so their validity is given by the quality of the adjustment of the theoretical semivariograms.
For the analysis of spatial behavior, we selected the contrasting ESRLs (ESRL 1 and ESRL 5), which will allow us to analyze the scenarios of greater and lesser microclimatic homogeneity. The other scenarios (ESRL 2, ESRL 3 and ESRL 4) are included graphically in Appendix 1 of this document.

Spatial variability of microclimate in MG
For the case of the aggregated data for the night hours (ESRL 1), it was found that for the MG greenhouse, the value of Tm was 11.1 ± 1.1 °C (Figure 7). It was observed that inside the greenhouse, the temperature in the largest amount of covered area ranges between 10.4 and 11.3 °C, while the maximum magnitude of the thermal differential inside the greenhouse (ΔTG) was 4.1 °C between small areas located in the central region and some points of the front sides. In the case of RH (Figure 7), the mean value was 93 ± 1.6 °C, with a very important area between values of 94% and 97%, while the maximum humidity differential (ΔRHG) was 11% between some perimeter points of MG. For the VPD variable, the mean value was 0.09 ± 0.03 kPa, while the VPD differential (ΔVPDG) was 0.16 kPa between areas located in the center and the perimeter of MG. The spatial behavior of temperature in ESRL 5 shows a greater climatic heterogeneity in MG (Figure 8). The Tm value was 19.2 ± 1.8 °C, while the ΔTG value was 6.2 °C. On the other hand, the mean RH value was 76 ± 11% with a ΔRHG value of 42%, finally, the mean VPD value was 0.52 ± 0.29 kPa and the ΔVPDG value was 1.05 kPa. Figure 8. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in MG and for ESRL 5.

Spatial Variability of Microclimate in TG
For TG of the night hours (ESRL 1), it was found that Tm presents a value of 11.5 ± The spatial behavior of temperature in ESRL 5 shows a greater climatic heterogeneity in MG (Figure 8). The Tm value was 19.2 ± 1.8 • C, while the ∆ TG value was 6.2 • C. On the other hand, the mean RH value was 76 ± 11% with a ∆ RHG value of 42%, finally, the mean VPD value was 0.52 ± 0.29 kPa and the ∆ VPDG value was 1.05 kPa.

Spatial Variability of Microclimate in TG
For TG of the night hours (ESRL 1), it was found that Tm presents a value of 11.5 ± 0.9 • C, while the value of ∆ TG was 0.8 • C, which can be considered a highly homogeneous value. In the case of RH and VPD, the mean values were 99 ± 0.2% and 0.006 ± 0.04 kPa, respectively, while the values of ∆ RHG and ∆ VPDG were 0.8% and 0.01 kPa (Figure 9). For ESRL 5, the values of T m , RH and VPD were found to be 23.1 ± 2.3 • C, 89 ± 7% and 0.41 ± 0.31 kPa, while the values of ∆ TG , ∆ RHG and ∆ VPDG were found to be 6.5 • C, 10% and 0.54 kPa (Figure 10). In the case of temperature, it was observed that there were three quite differentiated areas in magnitude of the energy value in the lateral, central and frontal zones of the TG, while for humidity and VPD, areas with different values were observed in different regions of the TG. The spatial behavior of temperature in ESRL 5 shows a greater climatic heterogeneity in MG (Figure 8). The Tm value was 19.2 ± 1.8 °C, while the ΔTG value was 6.2 °C. On the other hand, the mean RH value was 76 ± 11% with a ΔRHG value of 42%, finally, the mean VPD value was 0.52 ± 0.29 kPa and the ΔVPDG value was 1.05 kPa.

Spatial Variability of Microclimate in TG
For TG of the night hours (ESRL 1), it was found that Tm presents a value of 11.5 ± 0.9 °C, while the value of ΔTG was 0.8 °C, which can be considered a highly homogeneous value. In the case of RH and VPD, the mean values were 99 ± 0.2% and 0.006 ± 0.04 kPa,  respectively, while the values of ΔRHG and ΔVPDG were 0.8% and 0.01 kPa (Figure 9). For ESRL 5, the values of Tm, RH and VPD were found to be 23.1 ± 2.3 °C, 89 ± 7% and 0.41 ± 0.31 kPa, while the values of ΔTG, ΔRHG and ΔVPDG were found to be 6.5 °C, 10% and 0.54 kPa (Figure 10). In the case of temperature, it was observed that there were three quite differentiated areas in magnitude of the energy value in the lateral, central and frontal zones of the TG, while for humidity and VPD, areas with different values were observed in different regions of the TG.

Spatial Variability of Microclimate in SG
For SG, in regards to the spatial distribution of the microclimate in ESRL 1, it was found that some patches appear with a differentiated condition for temperature, relative humidity and VPD; these patches extend in the longitudinal direction of the greenhouse and appear in the same positions distributed along the width of the structure (Figure 11). For this condition, the Tm, RH and VPD values were found to be 11.6 ± 0.4 °C, 96 ± 3% and 0.07 ± 0.02 kPa, respectively, while the ΔTG, ΔRHG and ΔVPDG values were 1.4 °C, 15.1% and 0.13 kPa, although these extreme differential values only appear in small areas of SG.
For the case of ESRL 5, it was mainly observed that the areas of differing magnitude in the microclimate variables expand within SG, increasing the climatic heterogeneity of the greenhouse (Figure 12). The Tm, RH and VPD values were found to be 24.1 ± 3.9 °C, 72 ± 11% and 0.73 ± 0.35 kPa, respectively, while the ΔTG, ΔRHG and ΔVPDG values were 7.2 °C, 39% and 1.6 kPa, with extreme values appearing just above the central zone in SG.  respectively, while the values of ΔRHG and ΔVPDG were 0.8% and 0.01 kPa (Figure 9). For ESRL 5, the values of Tm, RH and VPD were found to be 23.1 ± 2.3 °C, 89 ± 7% and 0.41 ± 0.31 kPa, while the values of ΔTG, ΔRHG and ΔVPDG were found to be 6.5 °C, 10% and 0.54 kPa (Figure 10). In the case of temperature, it was observed that there were three quite differentiated areas in magnitude of the energy value in the lateral, central and frontal zones of the TG, while for humidity and VPD, areas with different values were observed in different regions of the TG.

Spatial Variability of Microclimate in SG
For SG, in regards to the spatial distribution of the microclimate in ESRL 1, it was found that some patches appear with a differentiated condition for temperature, relative humidity and VPD; these patches extend in the longitudinal direction of the greenhouse and appear in the same positions distributed along the width of the structure (Figure 11). For this condition, the Tm, RH and VPD values were found to be 11.6 ± 0.4 °C, 96 ± 3% and 0.07 ± 0.02 kPa, respectively, while the ΔTG, ΔRHG and ΔVPDG values were 1.4 °C, 15.1% and 0.13 kPa, although these extreme differential values only appear in small areas of SG.
For the case of ESRL 5, it was mainly observed that the areas of differing magnitude in the microclimate variables expand within SG, increasing the climatic heterogeneity of the greenhouse (Figure 12). The Tm, RH and VPD values were found to be 24.1 ± 3.9 °C, 72 ± 11% and 0.73 ± 0.35 kPa, respectively, while the ΔTG, ΔRHG and ΔVPDG values were 7.2 °C, 39% and 1.6 kPa, with extreme values appearing just above the central zone in SG.

Spatial Variability of Microclimate in SG
For SG, in regards to the spatial distribution of the microclimate in ESRL 1, it was found that some patches appear with a differentiated condition for temperature, relative humidity and VPD; these patches extend in the longitudinal direction of the greenhouse and appear in the same positions distributed along the width of the structure (Figure 11). For this condition, the T m , RH and VPD values were found to be 11.6 ± 0.4 • C, 96 ± 3% and 0.07 ± 0.02 kPa, respectively, while the ∆ TG , ∆ RHG and ∆ VPDG values were 1.4 • C, 15.1% and 0.13 kPa, although these extreme differential values only appear in small areas of SG. Sensors 2022, 22, x FOR PEER REVIEW 19 of 31 Figure 11. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in SG and for ESRL 1.

Figure 12.
Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in SG and for ESRL 5.

Spatial Variability of Microclimate in HGT1
The spatial behavior of the microclimate for the nighttime hours can be seen in Figure  13. Regions with a differential value in temperature, humidity and VPD were identified in various areas of the interior of HGT1, although it should be noted that the magnitude of the differentials appears to be low. The values of Tm, RH and VPD were 12.5 ± 0.12 °C, 99 ± 0.2% and 0.0018 ± 0.0011 kPa, respectively, while the values of ΔTG, ΔRHG and ΔVPDG were 0.25 °C, 0.3% and 0.0021 kPa. According to these values, it can be considered that the microclimatic condition in this greenhouse is highly homogeneous [54]. For the case of ESRL 5, it was mainly observed that the areas of differing magnitude in the microclimate variables expand within SG, increasing the climatic heterogeneity of the greenhouse (Figure 12). The T m , RH and VPD values were found to be 24.1 ± 3.9 • C, 72 ± 11% and 0.73 ± 0.35 kPa, respectively, while the ∆ TG , ∆ RHG and ∆ VPDG values were 7.2 • C, 39% and 1.6 kPa, with extreme values appearing just above the central zone in SG.

Spatial Variability of Microclimate in HGT1
The spatial behavior of the microclimate for the nighttime hours can be seen in Figure  13. Regions with a differential value in temperature, humidity and VPD were identified in various areas of the interior of HGT1, although it should be noted that the magnitude of the differentials appears to be low. The values of Tm, RH and VPD were 12.5 ± 0.12 °C, 99 ± 0.2% and 0.0018 ± 0.0011 kPa, respectively, while the values of ΔTG, ΔRHG and ΔVPDG were 0.25 °C, 0.3% and 0.0021 kPa. According to these values, it can be considered that the microclimatic condition in this greenhouse is highly homogeneous [54].

Spatial Variability of Microclimate in HGT1
The spatial behavior of the microclimate for the nighttime hours can be seen in Figure 13. Regions with a differential value in temperature, humidity and VPD were identified in various areas of the interior of HGT1, although it should be noted that the magnitude of the differentials appears to be low. The values of T m , RH and VPD were 12.5 ± 0.12 • C, 99 ± 0.2% and 0.0018 ± 0.0011 kPa, respectively, while the values of ∆ TG , ∆ RHG and ∆ VPDG were 0.25 • C, 0.3% and 0.0021 kPa. According to these values, it can be considered that the microclimatic condition in this greenhouse is highly homogeneous [54]. For the case of ESRL 5, it is again observed that under this level of radiation inside the greenhouse, there is greater microclimatic heterogeneity ( Figure 14). Under this condition, the values of Tm, RH and VPD were 22.4 ± 1.6 °C, 89 ± 5% and 0.28 ± 0.06 kPa, respectively, while the values of ΔTG, ΔRHG and ΔVPDG were 2.7 °C, 3.8% and 0.15 kPa. Figure 14. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT1 and for ESRL 5. For the case of ESRL 5, it is again observed that under this level of radiation inside the greenhouse, there is greater microclimatic heterogeneity ( Figure 14). Under this condition, the values of T m , RH and VPD were 22.4 ± 1.6 • C, 89 ± 5% and 0.28 ± 0.06 kPa, respectively, while the values of ∆ TG , ∆ RHG and ∆ VPDG were 2.7 • C, 3.8% and 0.15 kPa.

Spatial Variability of Microclimate in HGT2
In this greenhouse for the night hours, it was found that the values of Tm, RH and VPD were 13.2 ± 0.17 • C, 99 ± 1% and 0.09 ± 0.05 kPa, respectively, while the values of ∆TG, ∆RHG and ∆VPDG were 0.3 • C, 1% and 0.03 kPa. Further, this greenhouse, under this condition can be considered to present a homogeneous microclimate ( Figure 15). In the case of ESRL 5 for HGT2, some differentiated areas of microclimate were found, although to a lower percentage than in the greenhouses analyzed above (Figure 16). For this case, the values of T m , RH and VPD were 22.1 ± 0.7 • C, 86 ± 5% and 0.34 ± 0.12 kPa, respectively, while the values of ∆ TG , ∆ RHG and ∆ VPDG were 2.8 • C, 12% and 0.34 kPa.
Finally, these results obtained from the spatial variability analysis allowed a close relationship between the level of solar radiation and thermal heterogeneity inside the greenhouses to be identified; the higher the level of radiation, the greater the magnitude of thermal gradients inside the greenhouses and the greater the area of the internal surface of the structures, results that had been reported by Bojacá et al. [12]. Even this thermal heterogeneity in such naturally ventilated greenhouses has been reported by computational fluid simulation analysis [55]. On the other hand, it was also a common characteristic to find environments with quite high humidity contents and with VPD values very close to 0, both during the day and at night. This is undoubtedly a factor that may be limiting the yield of the crops grown in these greenhouses due to the negative implications that high-humidity environments have on plant growth. Figure 13. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT1 and for ESRL 1.
For the case of ESRL 5, it is again observed that under this level of radiation inside the greenhouse, there is greater microclimatic heterogeneity ( Figure 14). Under this condition, the values of Tm, RH and VPD were 22.4 ± 1.6 °C, 89 ± 5% and 0.28 ± 0.06 kPa, respectively, while the values of ΔTG, ΔRHG and ΔVPDG were 2.7 °C, 3.8% and 0.15 kPa. Figure 14. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT1 and for ESRL 5.  In this greenhouse for the night hours, it was found that the values of Tm, RH and VPD were 13.2 ± 0.17 °C, 99 ± 1% and 0.09 ± 0.05 kPa, respectively, while the values of ΔTG, ΔRHG and ΔVPDG were 0.3 °C, 1% and 0.03 kPa. Further, this greenhouse, under this condition can be considered to present a homogeneous microclimate ( Figure 15). In the case of ESRL 5 for HGT2, some differentiated areas of microclimate were found, although to a lower percentage than in the greenhouses analyzed above (Figure 16). For this case, the values of Tm, RH and VPD were 22.1 ± 0.7 °C, 86 ± 5% and 0.34 ± 0.12 kPa, respectively, while the values of ΔTG, ΔRHG and ΔVPDG were 2.8 °C, 12% and 0.34 kPa.  Finally, these results obtained from the spatial variability analysis allowed a close relationship between the level of solar radiation and thermal heterogeneity inside the greenhouses to be identified; the higher the level of radiation, the greater the magnitude of thermal gradients inside the greenhouses and the greater the area of the internal surface of the structures, results that had been reported by Bojacá et al. [12]. Even this thermal heterogeneity in such naturally ventilated greenhouses has been reported by computational fluid simulation analysis [55]. On the other hand, it was also a common character- Figure 15. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT2 and for ESRL 1. In this greenhouse for the night hours, it was found that the values of Tm, RH and VPD were 13.2 ± 0.17 °C, 99 ± 1% and 0.09 ± 0.05 kPa, respectively, while the values of ΔTG, ΔRHG and ΔVPDG were 0.3 °C, 1% and 0.03 kPa. Further, this greenhouse, under this condition can be considered to present a homogeneous microclimate ( Figure 15). In the case of ESRL 5 for HGT2, some differentiated areas of microclimate were found, although to a lower percentage than in the greenhouses analyzed above (Figure 16). For this case, the values of Tm, RH and VPD were 22.1 ± 0.7 °C, 86 ± 5% and 0.34 ± 0.12 kPa, respectively, while the values of ΔTG, ΔRHG and ΔVPDG were 2.8 °C, 12% and 0.34 kPa.  Finally, these results obtained from the spatial variability analysis allowed a close relationship between the level of solar radiation and thermal heterogeneity inside the greenhouses to be identified; the higher the level of radiation, the greater the magnitude of thermal gradients inside the greenhouses and the greater the area of the internal surface of the structures, results that had been reported by Bojacá et al. [12]. Even this thermal heterogeneity in such naturally ventilated greenhouses has been reported by computational fluid simulation analysis [55]. On the other hand, it was also a common character- Figure 16. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT2 and for ESRL 5. In the same way, it is important to remember that this type of behavior of nonhomogeneous microclimates generates differentiated yields in the crops established in these greenhouse structures [56,57]. Even this climatic heterogeneity affects plant breeding research experiments conducted in small greenhouses, causing noise and uncertainties in growth data collected at sampling points inside the greenhouse [58]. Therefore, at the local level, the technical managers of the farms observe accelerated production cycles in some areas within the same greenhouse, as well as low yields and poor-quality products in other regions of the structure [3]. This is even more critical because the management of cultural, agronomic, irrigation and fertilization tasks are established for a greenhouse, assuming that the microclimate within the whole area is homogeneous [59]. Finally, in accordance with the work developed by Fatnassi et al. [60], this type of micro-climatic behavior allows the establishment of pests and diseases that generate large economic losses and can generate a high rate of reports of pest interceptions at the control points located at the ports of export, which are monitored by the sanitary organizations of each country, affecting the phytosanitary status of the exporting country.
In terms of sustainability of the production systems, it is important to mention that all the disadvantages discussed in this research work, together with the scientific evidence that has proven the negative environmental impact loads generated by the existing production systems in places where large concentrations of greenhouses are established, have led to a significant reduction in the environmental impact of the production systems [50,61]. They provide a great opportunity for the search for sustainable strategies in the local context to be implemented in Colombian protected agriculture. It is possible to propose circularity strategies, strategies for the microclimatic optimization of the greenhouses analyzed. This is aimed at improving the yields of production systems, increasing the efficiency of irrigation and fertilization resources and promoting a lower use of chemical products for pest and disease control due to lower incidence and damage events [62].
Regarding the strategies for air conditioning, it should be mentioned that these should be based on the architectural redesign of the structures, increasing the ventilation surfaces above 40% of the covered area, seeking to maximize the renovation rates to values suitable for humid tropical conditions, guaranteeing renovation rates above 60 vol h −1 [63]. Strategies based on the use of photovoltaic agriculture or other renewable energy sources to be implemented for heating, cooling and dehumidification are another possibility widely implemented in recent years [64]. At the local level, there is also a lack of implementation and experimentation with passive technologies, such as other types of plastic coverings with special additives, such as anti-drip films to limit the effect of condensation, the use of shading nets to manage temperature and internal radiation at times of high solar radiation. For nighttime hours, possible solutions include the use of energy-saving thermal screens to limit thermal inversion phenomena, as well as the use of double roofs and green or plastic padding. Finally, feasibility studies and experimentation with passive thermal storage systems for heating purposes are also necessary.

Conclusions
In this work, a microclimatic characterization of the five main types of greenhouses used in Colombia was carried out. As relevant results, it can be mentioned that in each of the structures evaluated, it was found that there is a direct relationship between the conditions of temperature, relative humidity and DPV generated inside the greenhouse and the level of external solar radiation.
It was found that by using 40 temperature and relative humidity sensors uniformly distributed in each of the structures and installed at the same level above ground level, it is possible to predict the spatial variability of temperature, relative humidity and VPD in each of the 5 greenhouses evaluated by using geostatistical techniques at unsampled points. Likewise, this spatial variability analysis allowed the heterogeneous microclimate conditions generated in this type of passive greenhouse to be identified; the mean nighttime temperature conditions varied for the same greenhouse interior environment from a minimum of 0.25 • C in HGT1 to a maximum of 4.1 • C in MG, while for daytime conditions, this value varied between a minimum of 2.7 • C in HGT1 and a maximum of 7.2 • C in SG.
In terms of relative humidity and VPD, it was found that there is also heterogeneity in the spatial behavior inside each greenhouse evaluated. In addition, it was identified that the dominant average conditions inside the greenhouses are high humidity and inadequate DPV values; this behavior can be a productive limiting factor of the crops and the main cause of loss of quality of the harvested products due to fungal diseases and attacks.
Due to these results, it is important to evaluate climate optimization strategies at the local level to improve the microclimatic behavior of the greenhouse structures used in Colombia; this will undoubtedly contribute to increasing the sustainability of the most relevant production systems in the country. Likewise, with the rise of agriculture 4.0 and the internet of things, climate monitoring practices and the use of microclimatic information collected for decision-making in agronomic crop management practices should be promoted among Colombian producers.
Finally, it is important to mention that methodologically, this work did not include cultivation variables and measurements related to plant growth, neither physical nor physiological variables. This was due to time constraints for the development of the microclimatic characterization of each of the greenhouses. However, in the future, complementary experimental studies can be proposed to monitor the growth and productivity of the main crops of ornamental or horticultural interest; under these types of fluctuating and heterogeneous protected environments. Likewise, these future studies will allow us to know the responses that plants can offer in the short and long term to environmental conditions since crops have acclimatization mechanisms that can influence physiological or agronomic responses. The above does not detract from the scientific importance of the research work carried out since the microclimate management of greenhouse crops requires the greatest possible knowledge and understanding on the part of decision-makers and producers themselves; therefore, identifying the relationships that exist between the microclimate generated in a type of structure and the local climatic conditions is key to the management of this type of greenhouse. 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)".

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement: Not applicable.
Acknowledgments: To the Universidad Jorge Tadeo Lozano for the support for the development of the doctoral research project called "Development of a methodological proposal for the design, evaluation and microclimatic optimization of passive greenhouses through the use of numerical simulation tools".

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
In Figures A1-A3 Figure A1. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in MG and for ESRL 2.  In Figures A4-A6, the spatial variability curves for the TG greenhouse in ESRL 2, ESRL 3 and ESRL 4 are presented. Figure A4. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in TG and for ESRL 2.    In Figures A4-A6, the spatial variability curves for the TG greenhouse in ESRL 2, ESRL 3 and ESRL 4 are presented. Figure A4. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in TG and for ESRL 2. Figure A2. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in MG and for ESRL 3.   In Figures A4-A6, the spatial variability curves for the TG greenhouse in ESRL 2, ESRL 3 and ESRL 4 are presented. Figure A4. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in TG and for ESRL 2. In Figures A4-A6, the spatial variability curves for the TG greenhouse in ESRL 2, ESRL 3 and ESRL 4 are presented.   In Figures A4-A6, the spatial variability curves for the TG greenhouse in ESRL 2, ESRL 3 and ESRL 4 are presented. Figure A4. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in TG and for ESRL 2. Figure A4. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in TG and for ESRL 2.  In Figures A7-A9, the spatial variability curves for the SG greenhouse in ESRL 2, ESRL 3 and ESRL 4 are presented. Figure A7. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in SG and for ESRL 2. Figure A5. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in TG and for ESRL 3.  In Figures A7-A9, the spatial variability curves for the SG greenhouse in ESRL 2, ESRL 3 and ESRL 4 are presented. Figure A7. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in SG and for ESRL 2. In Figures A7-A9, the spatial variability curves for the SG greenhouse in ESRL 2, ESRL 3 and ESRL 4 are presented. Figure A5. Spatial distribution predicted by kriging, for temperature, relative humidity and in TG and for ESRL 3. In Figures A7-A9, the spatial variability curves for the SG greenhouse in E ESRL 3 and ESRL 4 are presented. Figure A7. Spatial distribution predicted by kriging, for temperature, relative humidity and in SG and for ESRL 2. Figure A7. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in SG and for ESRL 2.  In Figures A10-A12, the spatial variability curves for greenhous ESRL3 and ESRL4 are shown. In Figures A13-A15, the spatial variabilit house HGT2 in ESRL2, ESRL3 and ESRL4 are shown.   In Figures A10-A12, the spatial variability curves for greenhouse HGT1 in ESRL3 and ESRL4 are shown. In Figures A13-A15, the spatial variability curves for house HGT2 in ESRL2, ESRL3 and ESRL4 are shown. Figure A9. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in SG and for ESRL 4.
In Figures A10-A12, the spatial variability curves for greenhouse HGT1 in ESRL2, ESRL3 and ESRL4 are shown. In Figures A13-A15, the spatial variability curves for greenhouse HGT2 in ESRL2, ESRL3 and ESRL4 are shown. . Figure A11. Spatial distribution predicted by kriging, for temperature, relative humidity and in HGT1 and for ESRL 3. Figure A10. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT1 and for ESRL 2.
Sensors 2022, 22, x FOR PEER REVIEW 27 Figure A10. Spatial distribution predicted by kriging, for temperature, relative humidity and V in HGT1 and for ESRL 2.
. Figure A11. Spatial distribution predicted by kriging, for temperature, relative humidity and V in HGT1 and for ESRL 3. Figure A11. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT1 and for ESRL 3.  Figure A13. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT2 and for ESRL 2. Figure A14. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT2 and for ESRL 3.       Figure A13. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT2 and for ESRL 2. Figure A14. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT2 and for ESRL 3. Figure A14. Spatial distribution predicted by kriging, for temperature, relative humidity and VPD in HGT2 and for ESRL 3.