Performance analysis method for model-based irrigation strategies under uncertainty

Highlights • A Richards equation-based soil water model was combined with an evapotranspiration model to compute drainage and crop water stress resulting from irrigation within a MATLAB programming environment.• Monte Carlo sampling was used to simulate how uncertainty in soil parameters and evapotranspiration propagates into predictions of drainage, water use and crop water stress.• Soil water pressure head was used as a threshold value to indicate crop water stress based on soil profile interactions with the crop root zone.


Method details
In this contribution we present the code used in [3] to study the impact of uncertainty in evapotranspiration and soil parameters on drainage and crop stress predictions. For that, we adopted a model that mimics the real world in terms of the simulation model that describes water movement in soil (the EMMAN3G model). The software package (in Fortran) of the EMMAN3G, as described in [2] , is included. Details of the model can be found in [3] and are briefly repeated here.
Based on the uncertainty in model input and model parameters, a Monte Carlo analysis was developed in order to compute the uncertainty in predictions. This framework was developed in a MATLAB 2018b environment.

Richards equation
The Richards equation describes the change of water content in a soil column based on the water pressure head, hydraulic conductivity, water uptake, and irrigation: Here θ is the volumetric water content [ml cm −3  to the non-linear constitutive relationships between h, θ , and K , as presented next, and thus needs to be solved numerically.

The Mualem-Van Genuchten relationship
The Van Genuchten function ( Eq. (2) ; [1] ) describes the water retention as a function of the pressure head according to Here S is the effective saturation [dimensionless], θ r is the residual water content, θ s is the water content at saturation, and α [cm −1 ], n [dimensionless], and m [dimensionless] are curve shape parameters. The hydraulic conductivity characteristic is given by the Mualem equation ( Eq.
Here K r is the relative hydraulic conductivity [dimensionless], K s is equal to K at saturation [cm d −1 ], and λ is a curve shape parameter [dimensionless].

Root water uptake
Under well-watered conditions it is known that root water uptake is proportional to the root length density distribution. Since the main aim of crop production in greenhouse conditions is to maintain optimal conditions in the root zone, root water uptake was assumed to be proportional to fractional root distribution in the root zone times the potential transpiration. When this yields a pressure head exceeding a certain threshold value, this is regarded as a signal of water limitation.

The De Graaf model
The De Graaf model ( Eq. (4) ; [5] ) is used to calculate the evapotranspiration inside a greenhouse based on the global radiation received, the additional radiation by artificial lighting, the air temperature inside the greenhouse, the additional air temperature supplied by heating pipes, and the ratio between actual crop length and maximum crop length. This is formulated as follows, Here ET [mm] is the evapotranspiration, R [J cm −2 ] is the global radiation outside the greenhouse combined with the radiation from supplementary light, T tube [ °C] is the temperature from heating pipes, T gh [ °C] is the greenhouse indoor air temperature, a [mm cm 2 J −1 ] is an empirical crop factor for the effect of radiation, b [mm °C −1 ] is an empirical crop factor of the heating influence, L is the current crop length, and L max is the maximum crop length. In this study we focused on the final stage of the cropping period so that L = L max .
The ET is then divided in the potential transpiration (demand for crop water uptake) T pot and potential evaporation at the soil surface E pot according to [ mm ] (5) where LAI is the leaf area index. For a nearly full grown crop, as considered here, ET is dominantly assigned to T pot , and, therefore, we used T pot = ET and E pot = 0 .

The De Graaf model
First, the De Graaf model [5] was set up as a function

Starting the model framework and data input
A script is used to start the framework. First, any stored previous data is cleared and data sets of irrigation, temperature, and radiation are stored as single column double arrays. If multiple irrigation strategies are used then for each irrigation, the values are stored in a separate column, under irrigation.mat Next, the evapotranspiration function is run and stored as the variable ET

Simulation settings
An IF statement was used to select if the simulations will include uncertainty or not. The variables et_unc and soil_unc indicate if uncertainty will be present in the study ( > 0) or not ( = 0). After, the number of days the dataset spans are included. Next, the possible clay and sandy soils names were included in cell values named clay_soils and sand_soils . If evapotranspiration is uncertain, the evapotranspiration data is divided into the total number of days the dataset represents and stored as the variable ET_unc . Depending on the type of soil, an IF statement stores the clay or sand type of soil as the variable t ype_soil.
After selecting whether soil or evapotranspiration uncertainty will be included, the file PlantToSoil.dat is read as a ASCII delimited file. Then, the number of different irrigation strategies to run, and the number of sampling repeats using Monte Carlo sampling, are established as variables.

Monte Carlo sampling -evapotranspiration uncertainty
A FOR loop of the number of irrigation strategies is included as irrigations_run . Within the irrigations_run loop, another FOR loop of the number of times the framework will run depending on the number of Monte Carlo samplings as model_run . Within the model_run FOR loop, an IF statement and a FOR statement are added to indicate the Monte Carlo sampling of ET values. The FOR loop ET_unc_run includes a random sampling of ET_unc storing them as the variable ET_sampling to create a new ET variable for the current run simulation.

Monte Carlo sampling -soil uncertainty
An IF statement is added to indicate if soil uncertainty is included in the study when the variable soil_unc is higher than zero. Within the IF statement, a Monte Carlo random sampling of the soil parameters is included from the soil list type_soil and stored as the variable soil . The SoilInfo.dat is opened using the 'fopen' command and giving the ID fid_soil . Data from fid_soil is read using the command 'textscan' and stored as a string variable soil_data. Fid_soil is closed using the 'fclose' command, and the soilinfo.dat file is re-opened using the 'fopen' command with the option to rewrite. The new soil values stored in the variable soil_data are written into SoilInfo.dat using the 'fprintf' command. Finally fid_soil is closed using the 'fclose' command.

Sending the case studies to the EMMAN3G model
New values for evapotranspiration, temperature and irrigation are written on the previous established variable plant_soil . Evaporation and root length are commented out as it was assumed in the study [3] that the crop was fully grown, so evaporation from soil was neglected, and root length was assumed to be a constant of 30 cm. The file PlantToSoil.dat is re-opened with writing rights using the 'fopen' command. The headers for the data are included with the 'fprintf' command. The data from plant_soil is written to PlantToSoil.dat using the command 'dlmwrite' (ASCII delimited file) with a specific precision required by the EMMAN3G model and using '-append' to append the matrix to the file, if not, 'dlmwrite' overwrites the previous included headers.

Running the EMMAN3G model
Before running the EMMAN3G model, the previous output file is deleted. EMMAN3G model is run using the 'system' command. A new EmMan3G_W.csv file with be created by the EMMAN3G model.

Gathering results data
The EMMAN3G_W.csv file is opened and the number of columns are established as the variable columns_number . This is because, depending on the number of Z planes in the LocalFlux.dat file, the number of columns of the csv file will change. The data in the csv file is saved in the variable csv_data using the 'textscan' command. The fid is closed and the data is reshaped into the number of original rows and columns, and is stored in the variable data .

Storing results data as MATLAB variables
The variable data is in strings format and it is required to be as double. A FOR loop is used to store the results from the data variable as cells for the results of drainage, water pressure head of the first soil layer (8 Z-planes) and crop transpiration.
Each Monte Carlo sampling iteration outputs values of drainage, pressure head and transpiration predictions. The values stored in drainage, pressure_head , and transpiration are stored in new variables with the suffix _MC, which changes in which column the data is stored depending on the FOR loop model_run . Finally, the FOR loop model_run is closed after completing the Monte Carlo samplings.
Finally, the Monte Carlo sampling results for drainage, crop transpiration and water pressure head are stored in a cell which changes cell position depending on the irrigation strategy FOR loop irrigations_run . The FOR loop irrigations_run is closed after completing the different irrigation strategies previously established.

Defining crop stress ratio
The crop stress ratio consists of values of the mean water pressure head that are lower than the selected variable threshold. The crop is considered to be under stress when the mean water pressure is below the threshold. Three nested FOR loops were used to run the crop water stress ratio. The first loop runs for the number of irrigation strategies established ( stress_irrigation_run ). The next loop runs The mean pressure head value is obtained from the mean of the bottom 4 Z planes in the first soil layer. The mean value is stored in the variable mean_p_head . Also, the variable mean_p_head stores the mean pressure head value from different irrigation strategies studied in separate rows of the variable. If the values of mean_p_head are below the threshold value then it is stored as a risk event under the variable risk . The crop stress ratio represents the ratio of hours the crop was below the threshold value over the total amount of hours. The ratio values are stored in a cell ( final_risk ), this cell changes cell positions to store the ratio values of the different irrigation strategies within the FOR loop stress_irrigation_run.

Display of results
The mean and standard deviation values of the model outputs: drainage, risk of crop stress, and evapotranspiration was represented graphically. This was done with the FOR loop plot_results , storing the mean and standard deviation results in variables with the first column as the mean and the second column as the standard deviation. See Table 1 for a tabulated data example. Fig. 1. Example of a case study using the tabulated data. The shadowed area represents the uncertainty in the cumulative evapotranspiration. The bars represent the standard deviation of the prediction due to uncertainty in evapotranspiration.
A figure is created, the X-axis is established by including a variable ( total_irrigation ) which includes the different irrigation values used as irrigation strategies. The command 'Yyaxis right' is used to set an errorbar figure for the risk, using the plot_risk variable first column for the mean and the second column for the error. 'Yyaxis left' is used to set another errorbar for the drainage, using plot_drainage variable first column for the mean and the second column for the error. The uncertainty in evapotranspiration was represented as a set of limits between the mean and standard deviation of the plot_ET variable stored as the variable rx1 for the lower limit and variable rx2 for the higher limit. The space between rx1 and rx2 was filled an colored using the 'area' command. This colored area was made transparent using the 'FaceAlpha' and removing the line area with 'LineStyle'. See Fig. 1 as an example of the figure created using the tabulated data.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.