IntroductionTop

Material and methodsTop

Greenhouse condition and data acquisition

Two experiments were carried out in greenhouses located at the University of Chapingo, Mexico (19°29’ N, 98°53’ W; 2240 m a.s.l.). During the autumn-winter and spring-summer seasons, tomato (Solanum lycopersicum L) cultivar "CID F1" was grown in a hydroponic system. Plastic bags of 10 L capacity were used, which were f illed with volcanic sand (Tuff) as substrate. Plants were distributed with a density of 3.5 plants m-2 for both crop seasons. In the first experiment, tomato seeds were sown on July 18, 2015, and the seedlings were transplanted on August 21, 2015, in an 8 × 8 m chapel type glasshouse. For the second experiment, the seeds were sown on March 24, 2016, and the seedlings were transplanted on April 24, 2016, in 35 × 35 cm 12-L polyethylene bag pots in an 8 × 15 m plastic greenhouse with natural ventilation. Both experiments were fertilized with Steiner nutrient solution (Martínez-Ruiz et al., 2019; 2020). A weather station (HOBO, Onset Computer Corporation, Bourne, MA, USA) was installed inside the greenhouse. Temperature and relative humidity were measured with an S-THB-M008 model sensor (Onset) placed at a height of 1.5 m. Global Solar Radiation was measured with an S-LIB-M003 sensor (HOBO, Onset) located at 3.5 m above the ground. Both sensors were connected to a U-30-NRC model data logger (HOBO), and the data were recorded every minute, and subsequently the data were processed to obtain average data at hourly intervals.

In each experiment, three plants were randomly chosen and harvested every ten days to measure the DMP, LAI, Nup, Pup, Kup, Caup, and Mgup. The plants were dried for 72 h at 70 °C in a convection drying oven (Binder, ED-400 model). Nutrients in the stems, leaves and fruits were determined. Then 0.5 g dry matter was subjected to wet digestion with mixture of 5 mL of sulfuric acid and perchloric acid at a 4:1 ratio, both acids with 99% purity. Samples were digested until fully mineralized at a temperature of 250 ºC for approximately six hours; subsequently, 2 mL of hydrogen peroxide at 30% were added to the mineralized samples, and they were adjusted to 50 mL with deionized water. Nitrogen was determined by the Kjeldahl method (Sáez-Plaza et al., 2013). Phosphorus concentration was determined by the yellow molybdovanadate method, and K, Ca and Mg concentrations were measured by atomic absorption spectrophotometry (Oliveira et al., 2010). LAI was estimated by a nondestructive method which consisted of randomly taking four plants to measure the width and length of their leaves, and total leaf area was measured using a plant canopy analyzer (LAI-3100, LI-COR, Lincoln, NE, USA). From those measurements, nonlinear regression models were fitted to calculate the crop LAI. ETc was measured every minute by a weighing lysimeter located in the central row of each greenhouse. The device included an electronic balance Sartorius QA model (scale capacity=120 kg, resolution ±0.5 g equipped with a tray holding four plants. The weight loss measured was assumed to be equal to the ETc (Martinez-Ruiz et al., 2012).

 

Model description

The dynamic HORTSYST model assumes no water and nutrient limitations (Martínez-Ruiz et al., 2019; 2020), and it simulates PTI (MJ m-2 d-1), DMP (g m-2), Nup (g m-2), Pup (g m-2), Kup,(g m-2), Caup (g m-2), and Mgup (g m-2) as the state variables, while ETc (kg m-2) and LAI (m2 m-2) were considered as output variables. Table S1 [suppl.] lists the mathematical equations of the seven state and two output variables. Fig. 1 shows the general structure of the model using a Forrester diagram, where the inputs, parameters, state variables and outputs of the model are drawn.

The model structure was based on the VegSyst model (Gallardo et al., 2011; 2014; 2016; Giménez et al., 2013; Granados et al., 2013). The input variables of the model are hourly measurements of air temperature (°C), relative humidity (%), and global solar radiation (Wm-2). The models in the light use efficiency approach (Kang et al., 2008; Lemaire et al., 2008; De Reffye et al., 2009) allow calculation of daily ∆DMP (Eq. 9, Table S1 [suppl.]), as a function of the photosynthetically active radiation (PAR) (Eq. 6), crop characteristics such as LAI (Eq. 8), and the parameter of radiation use efficiency (RUE, g MJ-1) as has been used by Shibu et al. (2010) and Soltani & Sinclair (2012).

The fraction of light intercepted (fi-PAR) is the fraction of global solar radiation used for the photosynthesis process that enters through the canopy of a crop and is characterized by the LAI. The extinction coefficient (dimensionless k parameter) is related to leaf size and leaf orientation. Leaf area (LA) was modeled as a function of PTI using the Michaelis-Menten equation and it was multiplied by the planting density d in order to calculate the LAI. For this purpose, the normalized thermal time (TT, ºC), def ined as the ratio of the growth rate and the conditions of actual and optimum temperature ( Dai et al., 2006), was calculated with Eq. (4). Then, the daily ΔPTI (Eq. 7) was calculated as the product of normalized thermal time with the fraction of light intercepted (fi-PAR) and PAR radiation, and the accumulation of PTI was computed by Eq. (1) (Xu et al., 2010).

Daily ΔNup, ΔPup, ΔKup, ΔCaup and ΔMgup were calculated as the product of simulated DMP and nutrient concentration (Eq. 11). Previously the %N, %P, %K, %Ca and %Mg concentrations were determined by a power equation (Tei et al., 2002) (Eq. 10). Then their accumulated values were computed by Eqs. (12-16). Finally, the ETc was simulated hourly by Eq. (18), using the data recorded for global solar radiation, vapor pressure deficit, the fraction of light intercepted, and LAI as shown in Eq. (17). The daily transpiration was accumulated for the 24-h period using Eq. (3).

 

Monte Carlo uncertainty method

Monte Carlo simulation is a statistical technique for stochastic modeling and analysis of error propagation in calculations. Its aim is to trace out the structure of the probability distributions of the model output variables. These distributions are mapped by quantifying the deterministic results (realizations) for a large number of unbiased random draws (Matott et al., 2009) from the individual distribution function of input data and model parameters (Monod et al., 2006).

The uncertainty analysis consisted of the following four steps (Monod et al., 2006): a) for the spring-summer season, uniform PDFs were selected for HORTSYST model parameters (Table 1). Other PDFs are possible but the only information available for model parameters were their nominal values. Nominal values for all parameters were taken from literature. The lower and upper limits of the uncertainty intervals were defined with 10% and 20% of parameters variation around their nominal values as listed in Table 1; b) Latin hypercube sampling (LHS) was applied to choose model parameters values for the generation of N=10,000 scenarios; c) the output variables predicted by the model were calculated for all N scenarios, running N model simulations, using the climatic data measured in the greenhouse (Figs. 2a, 2b, 2c); and d) for the predicted variables (PTI, LAI, DMP, Nup, Pup, Kup, Caup, Mgup and ETc), the following statistical indicators were calculated: minimum, maximum, mean, coefficient of variation (CV), skewness, and kurtosis, as well as the histograms.

 

Figure 1.  EForrester´s relational diagram for the HORTSYST model of a greenhouse tomato crop: inputs, outputs, state variables, and parameters of the crop model. State variables are represented by rectangles, rate variables by valves, parameters with a horizontal line, input variables with a circle and a horizonal line, and auxiliary variables with ellipses. Flows of material are represented by normal arrows and information flows with dashed lines.

 

 

 

The generalized likelihood uncertainty estimation (GLUE) uncertainty method

This method is based on the Monte Carlo simulation, in which parameter sets may be sampled from some probability distribution function (PDF). The most commonly used PDF is a uniform distribution. Parameters values are also sampled from those PDF. Each parameter set is used to produce a model output; the acceptability of each model run is then assessed using a goodness-of-fit criterion which compares the predicted and observed values over some calibration period. As part of the GLUE procedure (Beven & Freer, 2001; Makowski et al., 2002; Stedinger et al., 2008; Beven & Binley, 2014), several likelihood functions can be used such as RMSE (root mean square error), inverse error variance, efficiency index, among others.

The GLUE procedure was applied to the HORTSYST model using LHS with 2,000 samples. Thus, 2,000 simulations were run using 10% and 20% of parameters variations around the nominal values of the 24 parameters (listed in Table 1), for the autumn-winter crop cycle. In order to confirm and compare the results that will be obtained with Frequentist Method what is considered the standard method for running a uncertainty analysis. The Matlab toolbox sensitivity analysis for everybody (SAFE) was used to implement the GLUE uncertainty analysis method since SAFE contains a module that allows not only per-forming uncertainty analyses but also it includes visua-lization tools such as scatter plots. This toolbox is freely available from the authors for non-commercial research and educational uses (Pianosi et al., 2015).

 

Table 1. HORTSYST model parameters with 10% and 20% of the variation of their nominal value, used for uncertainty simulation under the experimental condition for the spring-summer and autumn-winter crop cycles.

[1] Martinez-Ruiz et al. (2020). [2] Martinez-Ruiz et al. (2012). [3] Challa & Bakker (1999)

 

 

 

 

Figs. 2a, 2b and 2c show the daily averaged measu-rements of the integration of global solar radiation, air temperature, and relative humidity inside the greenhouses for the autumn-winter and spring-summer seasons. These data were fed as input variables into the HORTSYST mo-del. Input variable values represent two different weather conditions since averaged values of global radiation were 3.99 and 10.59 MJ m-2 d-1, respectively, although avera-ged values of air temperature were 18.3°C and 17.8°C for both seasons. However, relative humidity was also somewhat contrasting not because of its averaged values of 78.6% and 76.8%, but rather for its minimum values of 62.5% and 29.5% compared to its maximum values of 93.4% and 93.2%, respectively. Thus, the variation in re-lative humidity during spring-summer (63%) was roughly twice that for autumn-winter (31%).

Model output uncertainty with Monte Carlo method

HORTSYST model behavior regarding PTI, LAI and DMP is shown in Figs. 3a, 3c, 3e, respectively. Measu-red values are only reported for LAI and DMP, becau-se PTI is computed during the simulations and has not been measured. The corresponding histograms are also reported (Figs. 3b, 3d, 3f). These results were obtained using the following input variables: global solar radia-tion, temperature, and relative humidity recorded over the spring-summer cultivation period (Figs. 2a, 2b, 2c) and the 10,000 scenarios generated by the Monte Carlo procedure. The model simulation results for ETc, Nup and Pup (Figs. 4a, 4c, 4e) and their histograms (Fig. 4b, 4d, 4f) are presented. Finally, model predictions for Kup, Caup, and Mgup (Figs. 5a, 5c, 5e) together with their respective histograms (Figs. 5b, 5d, 5f) are shown.

Descriptive statistical measures calculated for all HORTSYST model predicted variables are summarized (Table 2). It is worthwhile noting that those quantities were calculated at the end of the cultivation period. In general, the uncertainty of the model predictions was increased with larger uncertainty intervals, as expected. The CV values of all variables were larger for 20% than 10% of the uncertainty variation in the model parameters. Furthermore, all predicted variables have CV values lower than 15% in the case of a 10% variation in the parameter values, which means the model is highly reliable and robust. Low CV values mean a reduction of average values for all predicted variables. Measured values at the end of the cultivation cycle were 6.86 m2 m-2 for LAI, 1304 g m-2 for DMP, and 291.69 kg m-2 for ETc. When these values were compared with HORTSYST model predictions, namely averaged values of predicted variables in Table 2 with 10% parameter variations, the deviations (difference between observed and estimated values) were: -1.7% for LAI, -1.7% for DMP and -0.3% for ETc. In the case of macronutrients, measured values and their corresponding deviations were: 27.4 g m-2 (1.6%) for Nup, 8.59 g m-2 (-6.9%) for Pup, 68.76 g m-2 (6.1%) for Kup, 20.42 g m-2 (2.9%) for Caup, and 7.63 g m-2 (3.3%) for Mgup. This means the average predicted values of LAI and DMP were slightly over-estimated and the predicted average values of ETc were close to the measured ones, whereas they were under-estimated for Nup, Kup, Caup and Mgup, and only Pup was over-estimated with respect to the observed data.

When the simulation used larger uncertainty intervals for model parameters (parameters varied 20% around their nominal values), the errors calculated from residuals between estimated (averaged values of predicted variables in Table 2 with 20% parameter variations) and measured values were as follows: 0.1% for LAI, 0.7% for DMP, 1.5% for ETc, 3.2% for Nup, -5.0% for Pup, 8.2% for Kup, 5.0% for Caup, and 3.9% for Mgup. Except for Pup, all the average values predicted by the model were under-estimated. Remarkably, the averaged values of LAI, DMP, and Pup were closer to the measured data with larger uncertainty. The average error of Nup, Pup, Kup, Caup and Mgup was higher with a larger uncertainty interval for model parameters. In the case of ETc, the average was more underestimated when the uncertainty included in the parameters was larger.

The differences between maximum and minimum values, for 10% variation of the parameters (Table 2), were 255.48 MJ m-2 d-1 for PTI, 4.74 m2 m-2 for LAI, 1014.9 g m-2 for DMP, 297.66 kg m-2 for ETc, 22.94 g m-2 for Nup, 8.09 g m-2 for Pup, 66.84 g m-2 for Kup, 18.27 g m-2 for Caup, and 7.04 g m-2 for Mgup. In the case of 20% variation (Table 2) these values were 255.48 MJ m-2 d-1 for PTI, 9.3 m2 m-2 for LAI, 1858.7 g m-2 for DMP, 579.43 kg m-2 for ETc, 50.7 g m-2 for Nup, 16.03 g m-2 for Pup, 121.94 g m-2 for Kup, 32.62 g m-2 for Caup, and 13.96 g m-2 for Mgup. This shows that the intervals of the predicted values increased more than two-fold with a large variation of 20%. Larger differences between maximum and minimum values of all HORTSYST predicted variables were observed with larger uncertainty intervals for model parameters.

The skewness values were positive for all predicted variables, which means that data are more spread out to the right of the distribution, which was observed in the corresponding histograms for 10% parameter variation (Figs. 3b, 3d and 5f; Figs. 4b, 4d and 4f; Figs. 5b, 5d and 5f). All skewness values were remarkably close to zero, which means all predicted variables fit a normal distribution very well, but with greater variation (20% of uncertainty) of the parameters, more asymmetric distributions are expected. In fact, skewness values are all greater for all variables in the case of 20% parameter variations than for 10% ones. Kurtosis values of predicted variables (Table 2) slightly deviate for both uncertainty intervals. This means that for both situations the behavior of predicted variables is remarkably close to a normal distribution (Figs. 3b, 3d and 3f; Fig. 4b, 4d and 4f; Fig. 5b, 5d and 5f).

 

Figure 2.  Daily averaged values of: a) the integration of global solar radiation, b) the air temperature and c) relative humidity, measured inside the greenhouses located in Chapingo, Mexico during autumn-winter, 2015, and spring-summer, 2016. DAT: days after transplant.

 

 

Figure 3.  HORTSYST model predicted variables: PTI= photo-thermal time (a), LAI = leaf area index (c), DMP = dry matter production (e) and their corresponding histograms (b), (d) and (f) calculated with 10% parameter variation, using Latin hypercube sampling and data collected during the spring-summer season. Measured LAI and DMP are indicated by circles.

 

 

Figure 4. HORTSYST model predicted variables: ETc = crop transpiration (a), Nup = nitrogen uptake (c), Pup = phosphorus uptake (e), and their corresponding histrograms (b), (d), (f), calculated with 10% parameter variation, using Latin hypercube sampling and data collected during the spring-summer season. Measured ETc, Nup and Pup are represented by circles.

 

 

Figure 5. HORTSYST model predicted variables: Kup = potassium up-take (a), Caup = calcium uptake (c), Mgup = magnesium uptake and their corresponding histograms, calculated with 10% parameter variation, using Latin hypercube sampling and data collected during the spring-summer season. Measured Kup, Caup and Mgup are represented by circles.

 

Table 2. Statistical summary for HORTSYST model predicted variables: photo-thermal time (PTI), crop leaf area index (LAI), dry matter production (DMP) and crop transpiration (ETc), nitrogen (Nup), phosphorus (Pup), potassium (Kup), calcium (Caup) and magnesium uptake (Mgup) with 10% and 20% model parameter variations around their nominal values.

 

 

Model output uncertainty with GLUE method

Fig. 6 shows HORTSYST model predictions with the GLUE procedure with a 95% confidence interval around measured data when 10% model parameter variation was used. Fig. 7 shows the GLUE uncertainty method outcomes in the case of 20 % parameter variation. It is apparent that with a smaller uncertainty interval for model parameters (10% variation around their nominal values), there was lesser uncertainty in all predicted variables (Fig. 6) than when a larger uncertainty interval (20% variation) was used (Fig. 7).

Although scatter plots between each HORTSYST parameter and the predicted variables generated by the GLUE uncertainty method were constructed, only those plots where minimum RMSE values can be identified are shown in Figs. S1-S3 [suppl.]. The best RUE parameter value was between 4.0 and 5.5 MJ m-2 d-1 (Fig. S1c [suppl.]), which corresponds to the smallest values of RMSE. The best value for parameter c1 was between 2.5 and 3.3 m2 (Fig. S2a [suppl.]), and for c2 it was between 60 and 85 (Fig. S2b [suppl.]). In the case of Bd this was between 15 and 35 W m-2 kPa-1 (Fig. S1d [suppl.]). Figs. S1e and S1f [suppl.] show the parameter values of Nup; a (from 6.0 to 7.5 g m-2) and b (from -0.2 to -0.15).

According to the RMSE values shown on the scatter plots, the best value for parameter a2 was between 0.55 and 0.65 g m-2 (Fig. S2a [suppl.]), for b2 between -0.08 and -0.05 (Fig. S2b [suppl.]), for a3 between 2.5 and 3.5 g m-2 (Fig. S2c [suppl.]), for b3 between 0.08 and 0.12 (Fig. S2d [suppl.]), for a4 between 2.5 and -3.5 g m-2 (Fig. S2e [suppl.]), for b4 between -0.13 and -0.06 (Fig. S2f [suppl.]), for a5 between 1.7 and 2.3 g m-2 (Fig. S2g [suppl.]), and for b5 between -0.14 and -0.07 (Fig. S2h [suppl.]). The parameter b (the slope of nutrient concentration curve) for all macronutrient concentrations were negative, except for Kup. The scatter plots (Fig. S3 [suppl.]) between the plant density values (d) and the RMSE of HORTSYST predicted variables show that the best value for this parameter was between 3 and 4 plants m-2. Therefore, plant density plays an important role in the model’s general performance and it had a major effect on LAI and ETc. A density higher than three plants m-2 yielded better fitness values for Nup, Pup, Kup, Caup, Mgup and DMP (Fig. S3 [suppl.]), whereas the model’s performance was inferior with a density lower than three plants m-2.

 

 

Figure 6. H HORTSYST model predicted variables with the GLUE uncertainty method with 10% parameter variations. LAI: leaf area index (a), DMP: dry matter production (b), ETc: crop transpiration (c), Nup: nitrogen uptake (d), Pup: phosphorus uptake (e), Kup: potassium uptake (f), Caup: calcium uptake (g), Mgup: magnesium uptake (h). DAT: days after transplant. Measured variables are represented by circles. Continuous lines show 95 % confidence intervals.

 

 

 

Figure 7. HORTSYST model predicted variables with the GLUE uncertainty method with 20% parameter variations. LAI: leaf area index (a), DMP: dry matter production (b), ETc: crop transpiration (c), Nup: nitrogen uptake (d), Pup: phosphorus uptake (e), Kup: potassium uptake (f), Caup: calcium uptake (g), Mgup: magnesium uptake (h). DAT: days after transplant. Measured variables are represented by circles. Continuous lines show 95 % confidence intervals.

 

 

 

According to Monte Carlo uncertainty analysis, with small uncertainty intervals for model parameters the uncertainty of HORTSYST model predictions, quantified by the CV, was in decreasing order: ETc > Mgup > Kup > Nup > Pup > PTI > Caup > DMP > LAI. The CV for all variables ranged from 10% to 14%. Although similar behavior was observed for larger uncertainty intervals, the range of uncertainty variation (CV values ranged from 22% to 30%) was also larger as expected. In both cases, LAI and DMP were predicted more accurately than Mgup, Kup, Nup, Pup, PTI, Caup, and ETc. Considering the error of the mean value of Monte Carlo simulations against the measured data for output variables, it is apparent that the HORTSYST model had good predictive ability, even though the parameter values included more uncertainty. However, it was observed that with larger uncertainty the accuracy of LAI, DMP and Pup estimations was increased. This is because the uncertainty analysis used nominal values of model parameters (values taken from the literature) instead of parameter values estimated by model calibration. Nominal parameter values can be near or far away from the optimal parameter values. Further work is needed to compare model uncertainty quantification before and after parameter estimation by using an optimization procedure.

The GLUE uncertainty analysis procedure confirmed that HORTSYST predictions are reliable, according to the calculated 95% confidence intervals for both model parameter uncertainty intervals (10% and 20% parameter variations around nominal values). Furthermore, this method provides parameter ranges for fitting model predictions to measured data. The intervals for optimal parameter values were in agreement with those obtained by model calibration (Martínez-Ruiz et al., 2019; 2020). Also, the value of the parameters RUE, required for DMP, parameters a and b, which are needed for Nup, were found in the ranges reported by Gallardo et al. (2014; 2016) for the VegSyst model.

According to Pathak et al. (2012), estimation of the GLUE and Monte Carlo uncertainty methods is based on the assumption that model parameters are independent. However, it is unlikely that this is the case since there are other sources of uncertainty, such as the input variables of the HORTSYST model (solar radiation, air temperature and humidity), initial conditions of state variables and the equations that are part of the model structure. Therefore, this work could be further extended by including other uncertain factors.

Based on the results of this work, the Monte Carlo and GLUE uncertainty analysis methods are necessary and complementary approaches since the former does not directly use any measured data of predicted variables in contrast to GLUE in which the use of that information is compulsory.

On the other hand, Lopez et al. (2018) found better fit to the measurement of HORTSYST model against VegSyst (Gallardo et al., 2016), the RMSE and the mean absolute error resulted three times lower. The improvement in quality of the prediction of DMP by HORTSYST model can be explained by the good modeling of LAI and the introduction of PTI as state variable. For the quality of the predictions of the HORTYSYS model it is possible to use it in the development of a decision support system (DSS) for irrigation management and nutrition in greenhouses, like those implemented with the VegSyst model developed by Gallardo et al. (2014) for N management an irrigation in greenhouse crops in Mediterranean climates. Elia & Conversa (2015) applied the GesCoN-DSS for management of fertigation in open field vegetable crops and Pérez et al. (2017) implemented a cFertigUAL-DSS based on a model for the control of the fertilization dose.

The results obtained in this study indicate that uncertainty analysis using both the Monte Carlo and GLUE methods can help in quantifying uncertainties in HORTSYST model predictions. Due to the small uncertainty associated with the model outputs, the model provides acceptable predictions (and it is reliable) when the nominal values of its parameters are varied between 10% and 20% under two different growing conditions (crop season). However, more research work is needed to determine whether HORTSYST can be applied to irrigation and nutrient supply management in greenhouse tomatoes grown under soilless culture. For example, future studies could consider various crop densities, different irrigation and nutrient levels, and temperature-driven stress conditions.