Sensitivity analysis and parameterization of two agricultural models in cauliflower crops

Antonio Lidón (Lidón, A)

Universitat Politècnica de València, Instituto Universitario de Ingeniería del Agua y del Medio Ambiente. Camí de Vera s/n, 46022 València, Spain.

Damián Ginestar (Ginestar, D)

Universitat Politècnica de València, Instituto de Matemática Multidisciplinar. Camí de Vera sn, 46022 València, Spain.

Sofía Carlos (Carlos, S)

Universitat Politècnica de València, Dept. Ingeniería Química y Nuclear. Camí de Vera s/n, 46022 València, Spain.

Carlos Sánchez de Oleo (Sánchez de Oleo, C)

Universidad Autónoma de Santo Domingo, Escuela de Física, Laboratorio de Investigación en Matemática Aplicada. San Juan de la Manguana, Dominican Republic.

Claudia Jaramillo (Jaramillo, C)

Universitat Politècnica de València, Dept. Producción Vegetal. Camí de Vera sn, 46022 València, Spain.

Carlos Ramos (Ramos, C)

Instituto Valenciano de Investigaciones Agrarias. Apartado oficial 46113. Moncada (Valencia), Spain.



Aim of study: The development of a procedure to calibrate the LEACHM and EU-Rotate_N models for simulating water and nitrogen dynamics in cauliflower crops.

Area of study: Calibration was performed using experimental data obtained from measurements in a cauliflower crop sited in Valencia (Spain) region.

Material and methods: A procedure based on generalized sensitivity indices for time-dependent outputs was used to determine the most influencing model parameters, in order to reduce the number of parameters to be calibrated and to avoid overparameterization. The most influencing parameters were introduced in an optimization process that uses the experimental measurements of soil water and nitrate content to determine its optimal value and obtain calibrated models.

Main results: After this analysis, the most important hydraulic parameters found were the coefficients of Campbell’s equation for the LEACHM model and the soil water content at field capacity and drainage coefficient for the EU-Rotate_N model. For the N cycle, the most influencing parameters were those related with the nitrification, humus mineralization rate and residue decomposition for both models. Both calibrated models provided good simulation of soil water content with an error between 5-7%. However, larger errors in soil-nitrate content simulation were found, mainly in the period corresponding to the crop residues incorporation. The prediction of the calibrated models in a different plot gave error values of about 7-9% for soil water content, but for soil nitrate content errors computed were 34% and 58%.

Research highlights: After calibration, both models can be used to optimize the farmer water management and fertilization practices in horticultural crops, although in the N case further studies should be performed.

Additional keywords: soil water content; soil nitrogen; global sensitivity analysis; model calibration; Brassica oleracea.

Abbreviations used: DC (drainage coefficient); ETo (potential evapotranspiration); GSI (generalized sensitivity index); Ks (hydraulic conductivity at saturation); LEACHM (Leaching Estimation And Chemistry Model); NR (nitrification rate coefficient); NSE (Nash and Sutcliffe coefficient of efficiency); PCA (principal components analysis); RMSE (root mean square error); RRMSE (relative root mean square error); θFC (volumetric water content at field capacity); θSAT (volumetric water content at saturation).

Authors' contributions: The analysis was conceived and designed by AL, DG. The numerical computations were conducted by SC and CSO. The experimental design, sampling and analysis were conducted by CR, AL and CJ. The paper was written by AL, DG and SC. All authors read, reviewed and approved the final manuscript.

Citation: Lidón, A; Ginestar, D; Carlos, S; Sánchez de Oleo, C; Jaramillo, C; Ramos, C (2019). Sensitivity analysis and parameterization of two agricultural models in cauliflower crops. Spanish Journal of Agricultural Research, Volume 17, Issue 4, e1106.

Received: 19 Jun 2019. Accepted:

Copyright © 2019 INIA. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International (CC-by 4.0) License.

Funding agencies/Institutions Project / Grant
Spanish Ministerio de Economía y Competitividad, INIA RTA 2011-00136-C04-01

Competing interests: The authors have declared that no competing interests exist.

Correspondence should be addressed to Antonio Lidón:





Material and methods





The improvement in some agricultural practi­ces, as irrigation and fertilization, to reduce some of the environmental problems they produce, can be undertaken using crop numerical models (Makowski et al., 2006; Cannavo et al., 2008). Such models take into account different processes to simulate the soil water and N dynamics and plant growth for agricultural systems. These models, after calibration, allow the estimation of nitrate leaching, soil mineral N and soil water content for different crops under different irrigation, rainfall and fertilization conditions, being an inexpensive and fast technique to evaluate the effects of various agricultural management practices on N fluxes to groundwater and atmosphere (Kersebaum et al., 2007; Cannavo et al., 2008). However, these models require a large number of uncertain or unknown parameters and input variables, what represents the major source of inaccuracy on the model predictions (Lamboni, 2009; Stella et al., 2014). Some of these parameters can be obtained from previous studies but, in general, it is necessary to cali­brate the models using experimental measurements, obtained from field experiments that, usually, cannot be continuously performed and are difficult to obtain and expensive.

Consequently, the estimation of the uncertain para­meters from experimental data is an important task, as model predictions depend on the accuracy of the parameter estimates (Makowski et al., 2006). More­over, the models can present overparameterization and overfitting problems, as a large number of parame­ters are required, and few experimental measurements are available (Kersebaum et al., 2007). Thus, it is interesting to select the most influencing parameters, by carrying out a sensitivity analysis, and use them to calibrate the model, setting the rest of parameters to a nominal value (Monod et al., 2006). There are different types of sensitivity analysis for the input parameters of a simulation model that can be performed depending on the intended objectives (Saltelli et al., 2008).

A local sensitivity analysis is based mainly on the estimation of the local derivative of the model output regarding the model input parameters. This deriva­tive indicates the speed at which the output grows or decreases locally with respect to the values of the input parameters. The local sensitivity analysis provides infor­mation only at the point where the derivatives are cal­culated without taking into account the rest of the inputs variation range. Other type of sensitivity analysis is called global sensitivity analysis and assigns the uncertainty of the model output to the variation of the input parameters, allowing to calculate the sensitivities about the full variation range of the input parameters. Global sensitivity methods are much more expensive from the computational point of view, therefore, its use is limited by the computational cost of the evaluation of the model to study.

Global sensitivity analyses, such as the Latin Hyper­cube-One Factor at a Time (LH-OAT) and the Fourier Amplitude Sensitivity Test (FAST), which are based on a scalar output variable, have been used to calibrate two agricultural computational models (Sánchez de Oleo, 2016). Nevertheless, in agricultural systems the responses provided by the models are time-dependent variables, since they depend on the different seasons of the crop, management practices and weather con­ditions. Thus, the classical sensitivity analysis for scalar variables is difficult to be applied. A global sensitivity analysis for discrete-time model outputs, which are vectorial outputs, has been proposed by Lamboni (2009). This methodology is based on a complete factorial design together with a principal components analysis, and it is the one used here to establish an importance ranking for the parameters to be calibrated in the model. The set of most important parameters is used to calibrate two agricultural models using soil water content and soil N content measures in a given plot of a cauliflower crop. LEACHM (Leaching Estimation And Chemistry Model) and EU-Rotate_N are the agricultural models selected in this study. The first one, is a popular model for water and N transport in soil, and the other one is specifically developed as a decision support system for water and N management in crop rotations. Both models have been widely used in different crops to simulate water and N dynamics in the soil-plant system (Fang et al., 2008; Doltra & Muñoz, 2010; Lidón et al., 2013; Soto et al., 2014; Suarez-Rey et al., 2016). In this work, the use of the time-dependent output global sensitivity analysis to select the most influ­encing parameters of LEACHM and EU- Rotate_N models was assessed. Finally, the predictive capacity of both calibrated models for another cauliflower plot, not used in the models calibration, was evaluated.

Material and methodsTop

Field experiments and measurements

Experimental data were obtained from field expe­riments conducted in two commercial plots (called Paterna-1 and Paterna-2) in Valencia province (39º 29’ 32’’ N and 0º 26’ 21’’ W, 14 meters above sea level), Spain. The soil in both plots is a deep alluvial one with clay loam texture and without stones. Main soil properties for Paterna-1 and Paterna-2 plots are shown in Table 1. In Paterna-1 plot, cauliflower (Brassica oleracea L. var. botrytis) plants were transplanted on September 17th, 2012, with a density of plantation of 23,420 plants ha-1. The crop was harvested on February 20th, 2013 and yield was 41.9 t ha-1. The crop residues were incorporated to soil on March 21st 2013. In Paterna-2 plot transplantation was on September 16th, 2013, with a plant density of 21,140 plants ha-1. The cauliflower was harvested on February 15th, 2014 and yield was 33.7 t ha-1. In this plot, the crop residues were incorporated to soil on March 4th 2014. The cultivar used in both plots was Triomphant, which is a long cycle type (about 150 days).

Table 1. Soil physical and chemical properties of the experimental plots Paterna-1 and Paterna-2, at three different depths (0-15, 15-30 and 30-45 cm).

In both plots, irrigation was by furrows and the water applied was measured at each irrigation event, and its content of nitrate, ammonium and chloride analyzed. In Paterna-1, six irrigations were applied being the total water depth of 698 mm. In Paterna-2, there were nine irrigations with a total amount of water applied of 815 mm. The N fertilizer in the experimental area of Paterna-1 plot was 100 kg N ha-1 and was supplied on October 29th, 2012, while in the area of Paterna-2 plot, 100 kg N ha-1 were applied on October 13th, 2013. In both cases, the fertilizer used was ammonium sulphate.

The meteorological data were collected from Manises weather station located less than 2 km away from the two plots. In the annual period, between September of 2012 and May of 2013, the precipitation was of 527 mm

and the potential evapotranspiration (ETo) was of

770 mm, while between September 2013 and May 2014 total rainfall was only 99 mm and ETo was 670 mm. The monthly potential evapotranspiration ranged from 39.2 mm in December 2013 to 178.6 mm in July of the same year, which is between 1.3 and 5.8 mm day-1 of average daily evapotranspiration.

In each plot, three soil samples were taken at three depths (0-15, 15-30 and 30-45 cm) with a periodicity of 15-21 days during the cultivation period and after the incorporation of the crop residues, until the beginning of the next crop at the middle of May. For each soil sample, moisture, mineral N content (nitrate and ammo­nium) and chloride content were determined. The gravimetric water content was determined by drying the soil at 105ºC, and the nitrate and ammonium content, by extraction in 2 mol L-1 KCl and subsequent colorimetric determination by flow injection analysis. Crop yield and N plant uptake were also measured. During the grow season, five cauliflower plants were sampled at seven different dates, for measuring the total N content and dry matter content. For each plant sample, fresh above-ground plant material was weighted and a representative sample oven-dried at 65ºC until constant weight was achieved. The N content in the plant sample was determined with an elemental analyzer. A soil chloride balance was applied to estimate water drainage for every period between two soil-sampling dates. Nitrate leaching was calculated as the product of the drainage and the average nitrate concentration in the soil water at 45 cm depth.

LEACHM model

LEACHM is a process-based, one-dimensional model that simulates water and solute movement, and related chemical and biological processes, in the unsaturated soil (Wagenet & Hutson, 1989). The model describes the one-dimensional water flow in the unsaturated zone using Richards’ equation and the solutes transport using the convection-dispersion equation. The main pro­­cesses described in the N module are mineralization, nitrification, denitrification and volatilization.

Water flow equation is solved for each soil layer

(5 cm in this case) and each water flow interval with a periodicity of 0.01 day. Equations relating volumetric water content (θ, cm3 cm-3), matric potential (h, kPa) and hydraulic conductivity (K, mm d-1) are required in the water flow model, which are obtained using Campbell’s equations (Campbell, 1974):

where θSAT is the volumetric water content at saturation (cm3 cm-3), a and b are constants. Parameter a can be interpreted as an air entry value and is given in kPa. Parameter b is a dimensionless fitting parameter. Ks is the hydraulic conductivity at saturation (mm d-1), and p is a pore interaction parameter, normally set to 1. In addition, this model uses the wet-end modification of Hutson & Cass (1987).

Input data for the LEACHM model include soil physical and chemical properties for the different soil layers as well as weather and crop data. The soil physical properties include bulk density, hydraulic conductivity and water retention curve parameters. Daily potential evapotranspiration is calculated as 1/7th of the weekly ETo (obtained with Penman Montheith method) multiplied by a crop factor. Transpiration potential (T) is obtained by multiplying the daily ETo by the he crop cover fraction, and potential evaporation is obtained as the difference between daily ETo and T. For each time step, LEACHM calculates potential ET assuming that potential evapotranspiration flux density (mm d-1) varies in a sinusoidal way throughout the day. The potential evaporation flux density during a time step is compared to the maximum possible evaporative flux density, using hydraulic parameters of the soil layers (matric potential, conductivity, and air-dry potential). This whole calculation process was based upon the methods of Childs & Hanks (1975). Water uptake by plant roots was calculated following Nimah & Hanks (1973).

LEACHM model contemplates different processes related to the N dynamics in the soil (Johnsson et al., 1987). The N cycle is linked to the carbon cy­cle, modelled using compartments: litter, humus and manure. The decomposition of the carbon pools is car­ried out by the soil microbial biomass according to a first order kinetics:

being Ci the carbon content (mg kg-1) and kmi the mineralization rate (d-1) in each compartment. The relation among the three carbon pools is given by the efficiency factor fe and the humification fraction fh, (Wagenet & Hutson, 1989). Regarding the N cy­cle, the main transformation processes considered are mineralization, immobilization, nitrification, deni­­-

tri­fication and volatilization. Three organic compart­ments (litter, humus and manure) and three mineral compartments (urea, nitrate and ammonium) define

N cycle in LEACHM model. The organic N trans­formations follow a first order kinetics, as for organic carbon. The nitrification process is described by the equation:

being knit the nitrification rate (d-1) and rmax the maximum relation NO3-/NH4+. Denitrification also follows a first order kinetics with respect to the nitrate content. The crop growth is taken into account considering the crop cover and root distribution as a function of time and depth. The N plant uptake takes place only when there is transpiration and is a function of the root density and the concentration of mineral N in the soil solution. Furthermore, this model uses a dimensionless curve relating the total extracted N with the growth period, which requires data on the maximum N crop uptake when there is not limiting nutrient levels in the soil. The potential daily uptake of N was calculated for each day according to Watt & Hanks procedure as described in Wagenet & Hutson (1989). In our study, the maximum N uptake for cauliflower crop was set to 250 kg ha-1.

EU-Rotate_N model

EU-Rotate_N is a decision support system for N management in crop rotations developed under the framework of a European Union project (Rahn et al., 2010). The model consists of a number of modules that simulate plant growth both, below and above ground, N mineralization from the soil and crop residues and subsequent N uptake, simulating the flow of water and N in the soil, into the plant, evapotranspiration and leaching. This model operates on a daily basis, utilizing soil properties data and crop residues, fertilizer and weather data. The model is bi-dimensional, and the unit cell in which soil is divided was 50 × 50 mm. The lower soil boundary was at 2000 mm depth.

EU-Rotate_N uses a “tipping bucket” approach in the soil profile to simulate water dynamics. The basic properties of the soil layers are provided by the user and include soil moisture at permanent wilting point (θWP), at field capacity (θFC), and at saturation (θSAT). These hydraulic properties control water availability to the plant and allow drainage calculation. Crop evapotranspiration was calculated using the FAO approach (Allen et al., 1998). Potential evapotranspiration was calculated from the climate parameters of weather data file or can be given as an input by the user. The values of the basal crop coefficients and the length of the different development periods are included in a model database for most vegetable crops, and can be easily adjusted to local growing practices, if required. Water infiltration and redistribution in the soil follow a capacitance approach using a drainage coefficient (DC) that allows water transfer between layers above soil moisture at field capacity to be progressively controlled (in more than one day) and more or less rapidly depending on soil type (Ritchie, 1998). Water not infiltrating in the daily time step is stored on the surface and infiltrates the following day.

Nitrogen mineralization from soil organic matter is based on the approach implemented in the DAISY model (Hansen et al., 1991), which considers three organic compartments (soil humus, microbial biomass and added organic matter) with materials that decompose at two speeds (slow and fast). In addition, decomposition rate constants are affected by temperature and soil mois­ture. Nitrification follows Michaelis-Menten kine­tics. The denitrification process follows a similar equa­tion but, in this case, it depends on soil nitrate concentration. EU-Rotate_N also considers ammonia volatilization after the manure incorporation and the other organic amendments, according to the ALFAM model (Søgaard et al., 2002).

EU-Rotate_N incorporates a root growth module, which calculates the vertical and horizontal extension of the root system from the length of the plant roots and their distribution. This process depends on the soil temperature. The crop growth is a function of the daily concentration of N in the aerial biomass and of the maximum growth potential, which is defined by the user through the dry matter value predicted in the harvest, and two specific crop parameters, A0 and B0 (Rahn et al., 2010). N plant uptake was calculated according to the daily crop demand. Crop growth is given by dry matter increment according to Greenwood’s model (Greenwood et al., 1991). This model has a target dry matter approach for crop growth, and the user has to provide the potential aboveground crop dry matter attainable for the given cultivar, soil, climate and cultivation practices in that region.

PCA-based multivariate global sensitivity analysis

This sensitivity analysis is based on a complete factorial design, taking into account that there is no variability in the model output for a given value of the input parameters, that is, that we work with a deterministic model.

The output of a discrete time-dependent model is denoted here as y(t,m), where m is a vector of m values for the input parameters and t is a discrete time, which takes values from 1 to tf. It is assumed that each parameter, or factor, can take only a finite number of values inside its range of variation, called levels. A complete factorial design for a model with m parameters and p levels considers n = pm possible scenarios to be computed. For this reason, when the number of parameters is large, complete factorial designs are performed using two or three levels at most. In this work, two levels for each parameter are considered, the upper and lower values in its rage of variation.

If the values of the model parameters are changed, according to a factorial design, the output can be stored in a matrix n × tf where n is the number of scenarios to be considered:

An orthogonal factorial design is used in the simu­lations, and an ANOVA analysis is performed sepa­rately on each output variable (Monod et al., 2006; Lamboni, 2009). The response variability is given by:

where µ is the average of the outputs, SSi is the sum of the squares associated with the main effect of parameter mi, SSij is the interaction between the parameters mi and mj, and SS1...m is the interaction between up to the m factors. The first order sensitivity indices (Monod et al., 2006) are defined as:

and the total sensitivity indices are defined as:

taking into account the interaction of factor i with the other factors.

To reduce the dimension of the problem, a principal components analysis (PCA) of matrix M is performed. In this way, Mc is the matrix M which each column centered around its mean (subtracting to each element the mean of the values of its column) and possibly normalized (data are divided by their standard deviation). The PCA decomposition is based on the eigenvalues and eigenvectors of MTc Mc. If λ1, … , λc are the eigenvalues of MTc Mc, and L the matrix of normalized eigenvectors, then the principal components matrix is

H = McL (8)

When an orthogonal factorial design is used, ANOVA-based sensitivity analysis can be applied to each prin­cipal component hl obtaining first order sensitivity indices and total sensitivity indices for each one of the principal components considered, l = 1, 2, …, Pc, being Pc<<tf (Lamboni, 2009). In practice, only the first principal components carry useful information on the model output and higher-order interactions can be neglected. To summarize the information given by the PCA based sensitivity indices, a generalized sensitivity index, GSIi, is defined for each parameter involved in the sensitivity analysis. GSIi is equal to the weighted sum of sensitivity indices of the Pc principal components with weights proportional to the inertia associated with the Pc principal components. The computation of the generalized sensitivity index is implemented in R language in the library multisensi (Lamboni, 2009; Lamboni et al., 2011), which is the library used to perform the computations presented here.

Models calibration methodology

The simulation period for both models has been from September 3rd of 2012 until May 30th of 2013 for Paterna-1 plot and from September 10th of 2013 until May 5th of 2014 for Paterna-2 plot. The analysis has been carried out in two stages; first, the most important parameters of the soil water dynamics are calibrated. Then, setting the value of the parameters of the water model to its optimal value, the parameters related to the N cycle are calibrated. For calibration, data from Paterna-1 plot were used, since more data for the irrigation and N fertilization were available, and this provided a wider range of water and N contents.

As in Jung et al. (2010), the parameters selected for the sensitivity analysis, in LEACHM model, have been the constants a and b of the Campbell’s equation, and the saturated hydraulic conductivity. Their maximum and minimum values, shown in Table 2, are similar to the ones proposed by Jung et al. (2010). Nevertheless, the range of variation for the Ks values was estimated using the expressions proposed in Saxton et al. (1986).

Table 2. Range of values used for the hydraulic parameters in the sensitivity analysis for LEACHM and EU-Rotate_N models for the three soil layers considered (0-15, 15-30 and 30-45 cm).

For the EU-Rotate_N model, the volumetric water content at field capacity (θFC), at saturation (θSAT) and the drainage coefficient (DC) for all the soil layers have been selected for the sensitivity analysis, as the most important parameters in irrigated crops (Lidón et al., 2011). Maximum and minimum values for these parameters have been obtained from Saxton et al. (1986), using the texture, organic matter and bulk density of each soil layer. For the drainage coefficient a value of 0.5 is recommended (Lidón et al., 2011), but a range of variation between 0.1 and 1 has been considered (see Table 2).

Once the most important parameters related with the water dynamics were obtained by determining the GSIi indices, they were used in an optimization process to obtain their best value by minimizing the normalized difference between simulated and measured values of the soil water content (εw), using the expression:

where nlayers is the number of soil layers considered, nmeas is the number of measures available in each layer, θvij (cm3 cm-3) is the j-th measure of the soil water content in the i-th layer and θ* (cm3 cm-3) is the j-th value of the soil water content in the i-th layer computed by the model. To obtain the optimal value of the parameters the Nelder-Mead simplex algorithm as described in Lagarias et al. (1998) has been used. An analogous process has been followed to find the most influential parameters of N dynamics.

In the LEACHM model, most of the selected para­meters, shown in Table 3, are related with N trans­formations in the soil (Sogbedji et al., 2006; Jung et al., 2010; and Lidón et al., 2013). Particularly, the rate constant for mineralization of humus, residues and organic manure, nitrification and denitrification were considered. According to Jung et al. (2010), other chemical properties as adsorption coefficient for NH4+ and molecular diffusion coefficient have also been taken into account.

Table 3. Range of values used for the parameters of nitrogen cycle in the sensitivity analysis for LEACHM and EU-Rotate_N models.

Although EU-Rotate_N model simulates the N cycle in a different way, for the sensitivity analysis the parameters selected for this model are those associated with similar N transformations as the ones considered for LEACHM. Therefore, the parameters selected were the rate constant of decomposition of the soil organic matter, the efficiency factor, and the nitrification rate constant. Other parameters associated with the crop residues decomposition, such as carbon and N content, were established using laboratory measurements (Jara­millo, 2016).

With respect to the parameters related with cau­liflower growth, the target dry matter value was the maximum value measured in the field (14.5 t ha-1). The N critical curve coefficients A0 and B0 were modified with respect to the default values for cauliflower in the croptable file of the model. The value used for A0 was 3.26 and for B0 was 0.728, the dry matter content (9%) and the harvest index (0.34) were modified to better fit the measured values for cauliflower growth in the Valencia area (Ramos, 2014).

As it has been done in the water module, to calibrate the most important parameters to the N cycle, the following error function, εN, was minimized:

where Nvij is the j-th measure of the N-NO3- content in the soil in the i-th layer (kg ha-1) and N*vij is the j-th value of the N-NO3- content in the i-th layer (kg ha-1) computed.

To evaluate the agreement between simulated and measured values, several statistical indices were used: the root mean square error (RMSE), where a RMSE=0 indicates a perfect fit, and the relative root mean square error (RRMSE) to estimate the mean error. For quantifying how good the fit of the models the Nash and Sutcliffe coefficient of efficiency (NSE) has been used; a NSE=1 denotes a perfect fit while a negative value indicates that the measured mean value is a better predictor (Ritter & Muñoz-Carpena, 2013).


Calibration of water dynamics modules

For the sensitivity analysis, the GSI was computed for each one of the parameters of the models using a factorial sampling with two levels associated with the minimum and maximum values of the parameters shown in Table 2. The most sensitive parameter for soil water content in the LEACHM model was the pa­rameter b of Campbell’s equation in each of the three soil layers (see Fig. 1), which is the expected behavior since this parameter is an exponent in Campbell’s model. Total GSI values for this parameter in the three layers ranged from 0.23 to 0.28. The a parameter of Campbell’s equation had total GSI values for the three soil layers ranging from 0.02 to 0.06, whereas those corresponding to the saturated hydraulic conductivity varied between 0.003 and 0.01, indicating a very small incidence of this parameter on the soil moisture in these conditions. In the EU-Rotate_N model, the most sensitive parameters corresponded to θFC and DC, while the values related to volumetric water at saturation were not very sensitive. Total GSI values for field capacity in the different soil layers ranged from 0.12 to 0.28, while those related to the drainage coefficient ranged from 0.11 to 0.22. For the θSAT the total GSI varied between 0.003 and 0.007. The GSI provides an importance order of the model parameters and, to select the number of parameters to be considered in the optimization process, a sensitivity analysis has to be performed. As discussed in Sanchez de Oleo (2016), for this application the four most important parameters provided the optimal value of the error function, that is, to consider additional para­meters in the optimization process did not improve the simulations.

Figure 1. Global sensitivity indices (GSI) for hydraulic and nitrogen parameters of the LEACHM and EU-Rotate_N models (parameter codes in Tables 2 and 3).

The four parameters with the largest GSI values of the water module of each model were calibrated, and the obtained results, together with the initial values are shown in Table 4. The other hydraulic parameters were set to the mean value of their range of variation shown in Table 2. The calibration process modified the values of all the parameters selected in both LEACHM and EU-Rotate_N, and the error, εw, obtained was slightly lower in EU-Rotate_N (0.221) than in LEACHM (0.258).

Table 4. Initial and optimum values of the hydraulic and nitrogen parameters for LEACHM and EU-Rotate_N models, together with the error obtained.

Once calibrated, a good fit was observed between simulated and measured soil water content for both models. The RMSE was 8.3 mm for LEACHM model and 6.4 mm for EU-Rotate_N, with RRMSE values of 0.07 and 0.05, respectively (Table 5). The NSE index was higher for the EU-Rotate_N model (0.52) than for the LEACHM model (0.18). Fig. 2a shows the soil water content in the considered soil profile (0-45cm) for both models, together with experimental data. LEACHM showed higher values of soil water content after irrigation events or important rainfall, which is related to its physical soil character. However, in the compartmental EU-Rotate_N model, θFC was the most influencing parameter in water drainage and, therefore, the soil water content is mainly determined by this value.

Table 5. Statistical indices for the comparison between observed and simulated values for the calibration period (Paterna-1) and validation period (Paterna-2) in both models at 45 cm.

Figure 2. Paterna-1 plot measurements together with LEACHM and EU-Rotate N results for a) soil water content, b) drainage, c) nitrate content and d) nitrate leaching at 45 cm depth.

Since the evaluated parameters not only had an effect on soil water content but also on water movement through the soil profile, and therefore on the nitrate leaching, the accumulated drainage at 45 cm, simulated with both models, together with the drainage calculated by chloride balance, are shown in Fig. 2b. The EU-Rotate_N model predicted at the end of the period an accumulated drainage of 825 mm, higher than the value estimated by the chloride balance, which was 673 mm. EU-Rotate_N overestimated accumulated drainage throughout the simulation. LEACHM simu­lated a total cumulative drainage value of 528 mm, and underestimated drainage mainly after 100 days of simulation. Finally, in Fig. 3a, soil water content measures are plotted against the simulated data, sho­wing a similar behavior of both agricultural models.

Figure 3. Soil water content measures against the simulation with the two models for a) Paterna-1 and b) Paterna-2. Soil nitrate content measures against the simulation for c) Paterna-1 and d) Paterna-2.

Calibration of nitrogen modules

GSI indices obtained for N parameters calibration in both models are shown in Fig. 1. For LEACHM only three parameters out of the thirteen evaluated had total GSI values above 0.1: the nitrification rate coefficient, NR1 (0.64), the synthesis efficiency factor, SE (0.27) and the mineralization rate of the top layer soil, HMR1 (0.14). In the EU-Rotate_N model, GSI for the nitrification rate coefficient, NR, was 0.8, and for all other parameters it was lower than 0.07.

Again, the parameters with the four largest values of the total GSI were chosen for calibration of the N module. Table 4 shows their initial and optima values together with the error εN. The other variables were set to the mean value of their range of variation, shown in Table 3. After the calibration process the error obtained was εN =1.86 for the LEACHM model and εN =2.10 for EU-Rotate_N model.

Nitrate content in the soil profile (0-45cm) was simulated with both calibrated models, and results are shown in Fig. 2c. The agreement between measured and simulated values was reasonably good during the crop period, but after the incorporation of the crop residues both models clearly underestimated the soil mineral N content. This suggests that crop residue decomposition was not well simulated in any model, and this was not due to the use of a wrong value for the parameters, since the GSI indices for these parameters were quite low (Fig. 1). Table 5 shows the agreement indices between measured and simulated values of nitrate content in soil profile (0-

45 cm). RRMSE values for the whole simulation were high (0.7 in both models), and RMSE was about 50 kg N-NO3- ha−1 for both models. A negative value of the NSE index was obtained with both models. This indicates that the calibration process for the soil N dynamics was not satisfactory with both models and further analyses would be necessary. Finally, in Fig. 3c, soil nitrate content measures are plotted against the simulated data. Underestimation of the nitrate content is observed for both LEACHM and EU-Rotate_N models.

The accumulated nitrate leaching simulated and mea­sured at 45 cm depth is shown in Fig. 2d. N-NO3- leaching at 45 cm was underestimated throughout the period by both models, with a percentage of variation from the total estimated leaching of 48% and 24% for LEACHM and EU-Rotate_N, respectively.

Predictive ability of calibrated models

The two models with calibrated parameters using data from Paterna-1 plot were used to evaluate their ability to predict soil water and nitrate content in a different experimental plot (Paterna-2) located about 300 m away from the first one. Fig. 4 shows the measured and simulated soil water content down to 45 cm depth in Paterna-2. Table 5 shows the statistical indices for the soil water content in soil profile (0-45 cm) obtained in Paterna-2. The calibration performed with water parameters allows a reasonable prediction of the soil water content in the profile at the new plot. A good fit was observed between measured and simulated soil water content, with RRMSE values about 0.09 for LEACHM and EU-Rotate_N, (see validation period in Table 5) which are slightly higher than those obtained in the calibration period. Nevertheless, the NSE in­dex decreased with respect to the calibration period, being slightly negative for the EU-Rotate_N model. Furthermore, both models were able to predict well the accumulated drainage in the Paterna-2 plot (see Fig. 4b), especially the EU-Rotate_N model. Soil water content measured against the simulated values for Paterna-2 (Fig. 3b) show that, qualitatively, both models have a similar behavior to predict the soil water content. Figs. 4c and 4d show the measured and simulated soil nitrate content down to 45 cm with both models and nitrate leaching in plot Paterna-2, and Table 5 shows the statistical indices of agreement for the nitrate content in the soil profile (0-45 cm). Soil nitrate content was clearly underestimated by both models, although LEACHM seemed to behave better than EU-Rotate_N, with RRMSE values of 0.34 and 0.58, respectively.

Figure 4. Paterna-2 plot measurements together with LEACHM and EU-Rotate_N results for a) soil water content, b) drainage, c) nitrate content and d) nitrate leaching at 45 cm depth.

The RMSE ranged from 50 kg N-NO3- ha-1 for LEACHM model to 83 kg N-NO3- ha-1 for EU-Rotate_N model. The NSE indices increased with respect to the calibration period. Again, the main differences between measured and simulated values occurred after the incorporation of crop residues to soil, as observed in Fig. 4c and Fig. 3d. The cumulative nitrate leaching in the simulation period is underestimated by 27% and 84% by LEACHM and EU-Rotate_N, respectively. Since both models predicted water drainage reasonably well, the underestimation of nitrate leaching must be related to inadequate simulation of the soil N proces­ses in both models.


Calibration of LEACHM model

Calibration of LEACHM hydraulic parameters, based on the GSI sensitivity analysis, allowed a good simulation of soil water content with a RRMSE of 7% and a RMSE of 8 mm. As it was expected, parameter b of the Campbell’s equation was the most influen­tial parameter in soil water dynamics in a vegetable crop, similar results have been obtained by other authors (Jung et al., 2010; Jabro et al., 2011). This must be due to the influence of this parameter in the functional relations of the volumetric water content with matric potential and hydraulic conductivity, given by equation (1). This shows the importance of performing a good measurement and/or calibration of this parameter.

With respect to a and Ks parameters, since their influence is maximum close to saturation, it would be convenient to take soil water content measurements soon after an irrigation or a major rain event. The soil samples were not obtained with large values of soil moisture, so it was not possible to evaluate the real influence of these parameters. However, simulated drainage, was similar to the one obtained by chloride balance. This indicates that the mean values used for the Ks were adequate.

The greatest discrepancies between measured and simulated soil water content values occurred with the incorporation of crop residues. This required using a rotovator that produced a change in the soil bulk density and therefore in soil porosity, not taken into account in the simulation. The effect of the bulk density on the soil water and N dynamics has been evaluated in Mahmood & Tillman (2015), which recommend accurate measurements of bulk density. To introduce a possible change of the bulk density along the simulation period would improve the model results.

With respect to N module, the GSI sensitivity analysis clearly showed the nitrification rate, the synthesis efficiency factor and the mineralization rate of the residues as the most influential parameters. These parameters have been also identified as relevant in the model evaluation by other authors, as Schmied et al. (2000), Jung et al. (2010) and Mahmood & Tillman (2015). These authors also show that the adsorption coefficient of N-NH4+ is important at all conditions evaluated. In our approach, this last coefficient was considered in the N module calibration of LEACHM model, since the associated GSI index had a large value.

After calibration, the results for the soil nitrate content evolution resulted in a RRMSE of 0.72 and a RMSE of 49 kg N-NO3- ha−1. The model underestimated soil nitrate content, mainly after the incorporation of crop residues.

Calibration of EU-Rotate_N model

EU-Rotate_N model uses a simple approach to cal­culate soil water redistribution in the soil profile and requires few soil hydraulic parameters. GSI sensitivity analysis showed that θFC and DC are the most influential parameters on the soil water dynamics in the cauliflower crop. These results agree with other studies done in vegetable crops (e.g. Rahn et al., 2010; Suarez-Rey et al., 2019). Different results were observed for the sensitivity analysis of the EU-Rotate_N model (Doltra & Muñoz, 2010). Calibration of these parameters re­sulted in a good fit between the simulated water con­tent and the experimental measured values, giving a RRMSE of 5% and a RSME about 10 mm.

Although soil moisture at field capacity can be estimated by the pedotransfer functions, the drainage coefficient could be determined by the model according to the field capacity and saturation water contents, but this value has been included in the calibration process, as it has been done in Soto et al. (2014). However, the drainage was overestimated in the calibration period. This can be due to the daily time scale used by this model, which may not be adequate for some circumstances as furrow irrigation, where the water movement takes place in a smaller time scale.

The lack of water content measurements at or near saturation could be conditioning the non-sensitivity of volumetric moisture at saturation parameter in the simulations. It would be advisable in future works to have moisture measures that cover the entire range of soil moisture that is usually produced in agricultural soils under furrow irrigation.

The nitrification rate coefficient was practically the only parameter affecting the soil nitrate content according to the GSI sensitivity analysis. This is quite different from what was found for LEACHM model (see Fig. 1 and Table 4). The other parameters used in the calibration were related with the crop residues and soil organic matter mineralization. These results show that the fast residues decomposition rate and the slow soil organic matter decomposition rate are very important to model nitrate leaching and soil N content using EU-Rotate_N, as remarked in Doltra & Muñoz (2010), Nendel et al. (2013) and Soto et al. (2014). The optimization of these parameters gave a good fit of the soil nitrate content to the measured data during the crop period, but after crop residues incorporation there was an underprediction. The global RRMSE error was 0.73, which is similar to the one obtained with LEACHM and the mean difference was of 30 kg N-NO3- ha−1. Nitrate leaching simulation was better than that obtained with the LEACHM model, but still provided an underestimation.

Models predictions

Using the calibrated parameters, the prediction for soil water content in a different plot (Paterna-2) was quite reasonable for both models, giving RRMSE values of 8.7% and 7% for LEACHM and EU-Rotate_N, respectively. However, simulated values of soil N-NO3- content for the Paterna-2 plot gave RRMSE values of 0.34 and 0.58 for LEACHM and EU-Rotate_N, respectively, obtaining a slight better fit for LEACHM. Although the obtained values for RMSE are high, this range of values have been also reported in other works where the soil N content is simulated (Sun et al., 2013; Soto et al., 2014). Both models underestimate the nitrate content in all the simulation period, what indicates that more field data are needed to improve the N cycle calibration, as it is more complex than the water dynamics. In addition, a revision of the N dynamics models implemented in the codes would be necessary to improve the prediction capabilities of the models for horticultural crops when organic amendments are used and crop residues are incorporated.

In summary, LEACHM and EU-Rotate_N models have been compared for simulating soil water and N dynamics in a cauliflower crop. Using data from a two-year experiment, these models were calibrated with soil water and nitrate content evolution during a growing crop period and after crop residue incor­poration to soil. To avoid overparameterization, a global sensitivity analysis adapted to discrete time model outputs was used to determine the most influ­ential parameters. After this, the values of the four most influencing parameters were determined by means of an inverse optimization process. From the results obtained, we can conclude that these simulation models are able to predict the soil water content with better accuracy than the nitrate content. Water drainage is difficult to be measured in the commercial fields, but the simulated values with both calibrated models were in a reasonable agreement with estimates obtained using a chloride balance approach. Since processes associated with the N cycle are more complex, the calibration of the N module in both models was not successful and predictions of the soil nitrate content did not match the measurements, especially after the crop residues incorporation. Nitrate leaching at 45 cm depth, was better simulated by EU-Rotate N model in the calibration period, whereas LEACHM obtained better results in the prediction period. The obtained results show that for a better simulation of agricultural systems where different organic inputs are present, the implemented nitrogen dynamics models should be revised.


Allen RG, Pereira LS, Raes D, Smith M, 1998. Crop evapotranspiration. Guidelines for computing crop water requirements. FAO Irrig Drain Paper nº 56, Roma.

Campbell G, 1974. A simple method for determining unsaturated conductivity from moisture retention data. Soil Sci 117: 311-314.

Cannavo P, Recous S, Parnaudeau V, Reau R, 2008. Modeling N dynamics to assess environmental impacts of cropped soils. Adv Agron 97: 131-174.

Childs SW, Hanks RJ, 1975. Model of soil salinity effects on crop growth. Soil Sci Soc Am Proc 39: 617-622.

Doltra J, Muñoz P, 2010. Simulation of nitrogen leaching from a fertigated crop rotation in a Mediterranean climate using the EU-Rotate N and Hydrus-2D models. Agric Water Manag 97 (2): 277-285.

Fang Q, Ma L, Yu Q, Malone RW, Saseendran SA, Ahuja LR, 2008. Modeling nitrogen and water management effects in a wheat-maize double-cropping system. J Environ Qual 37 (6): 2232-2242.

Greenwood DJ, Gastal F, Lemaire G, Draycott A, Millard P, Neeteson JJ, 1991. Growth rate and % N of field grown crops: Theory and experiments. Ann Bot 67: 181-190.

Hansen S, Jensen HE, Nielsen NE, Svendsen H, 1991. Simulation of nitrogen dynamics and biomass production in winter wheat using the Danish simulation model DAISY. Fert Res 27: 245-259.

Hutson JL, Cass A, 1987. A retentivity function for use in soil-water simulation models. J Soil Sci 38: 105-113.

Hutson JL, Wagenet RJ, 1991. Simulating nitrogen dynamics in soils using a deterministic model. Soil Use Mana­ge 7: 74-78.

Jabro JD, Hutson JL, Jabro AD, 2011. Parameterizing LEACHM model for simulating water drainage fluxes and nitrate leaching losses. In: Methods of introducing system models into agricultural research; Ahuja LR & Ma L (eds.). pp: 95-115. Am Soc Agron, Crop Sci Soc Am, Wisconsin, USA.

Jaramillo C, 2016. Mineralización de la gallinaza y de los restos de cosecha en el suelo. Aplicación al cultivo de coliflor en la huerta de Valencia. Doctoral thesis, Univ. Politécnica, Valencia, Spain.

Johnsson H, Bergstrom L, Jansson PE, Paustian K, 1987. Simulated nitrogen dynamics and losses in a layered agricultural soil. Agric Ecosyst Environ 18 (4): 333-356.

Jung YW, Oh DS, Kim M, Park JW, 2010. Calibration of LEACHN model using LH-OAT sensitivity analysis. Nutr Cycl Agroecosys 87: 261-275.

Kersebaum KC, Hecker JM, Mirschel W, Wegehenkel M, 2007. Modelling water and nutrient dynamics in soilcrop systems: a comparison of simulation models applied on common data sets. In: Modelling water and nutrient dynamics in soil-crop systems; Kersebaum KC et al. (eds.). pp: 1-17. Springer, Dordrecht.

Lagarias JC, Reeds JA, Wright MH, Wright PE, 1998. Convergence properties of the Nelder-Mead simplex method in low dimensions. SIAM J Optimiz 9: 112-147.

Lamboni M, 2009. Multivariate global sensitivity analysis for dynamic crop models. Field Crop Res 113: 312-320.

Lamboni M, Monod H, Makowki D, 2011. Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models. Reliab Eng Syst Safe 96: 450-459.

Lidón A, Lado L, Berbegall F, Ramos C, 2011. Influencia de la calibración de los parámetros hidráulicos del modelo EU-Rotate_N en el balance de N en un cultivo de col china. Actas Hort 61: 45-51.

Lidón A, Ramos C, Ginestar D, Contreras W, 2013. Assessment of LEACHN and a simple compartmental model to simulate nitrogen dynamics in citrus orchards. Agric Water Manag 121: 42-53.

Mahmood B, Tillman R, 2015. Sensitivity of LEACHN model predictions to changes in Nitrogen transformation rate constants and bulk density. T ASABE 58: 353-366.

Makowski D, Hillier J, Wallach B, Andrieu B, Jeuffroy MH, 2006. Parameter estimation for crop models. In: Working with dynamic crop models; Wallach D et al. (eds.). pp: 101-149. Elsevier, Amsterdam.

Monod H, Naud C, Makowski D, 2006. Uncertainty and sensitivity analysis for crop models. In: Working with dynamic crop models; Wallach D et al. (eds.). pp: 55-99. Elsevier, Amsterdam.

Nendel C, Venezia A, Piro F, Ren T, Lillywhite RD, Rahn CR, 2013. The performance of the EU-Rotate_N model in predicting the growth and nitrogen uptake of rotations of field vegetable crops in a Mediterranean environ­ment. J Agric Sci 151: 538-555.

Nimah MN, Hanks RJ, 1973. Model for estimating soil water, plant, and atmospheric interrelations: I. Description and sensitivity. Soil Sci Soc Am J 37(4): 522-527.

Rahn C, Zhang K, Lillywhite R, Ramos C, Doltra J, De Paz JM, Riley H, Fink M, Nendel C, Thorup K, Pedersen A, Piro F, Venezia A, Firth C, Schmutz U, Rayns F, Strohmeyer K, 2010. EU-Rotate_N - a decision support system - to predict environmental and economic consequences of the management of nitrogen fertiliser in crop rotations. Eur J Hortic Sci 75: 20-32.

Ramos C, 2014. Integración de medidas de suelo, planta y modelos de simulación para el manejo eficiente del nitrógeno en los cultivos hortícolas. Project RTA2011-00136-C04-01, Final Report, INIA, Madrid, Spain.

Ritchie JT, 1998. Soil water balance and plant water stress. In: Understanding options for agricultural production; Tsuji GY et al. (eds.). pp: 41-54. Kluwer Acad Press, Dordrecht.

Ritter A, Muñoz-Carpena R, 2013. Performance evaluation of hydrological models: Statistical significance for reducing subjectivity in goodness-of-fit assessments. J Hydrol 480: 33-45.

Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, Saisana M, Tarantola S, 2008. Global sensitivity analysis: the primer. John Wiley & Sons Ltd, Chichester, England. 312 pp.

Sánchez de Oleo C, 2016. Estimación de parámetros en modelos de transporte de agua y nitrógeno en el suelo. Doctoral thesis, Univ. Politécnica, Valencia, Spain.

Saxton KE, Rawls WJ, Romberger JS, Papendick RI, 1986. Estimating generalized soil-water characteristics from texture. Soil Sci Soc Am J 50 (4): 1031-1036.

Schmied B, Abbaspour K, Schulin R, 2000. Inverse estimation of parameters in a nitrogen model using field data. Soil Sci Soc Am J 64: 533-542.

Sogbedji JM, Van Es HM, Melkonian J, Schindelbeck R, 2006. Evaluation of the PNM model for simulating drain flow nitrate-N concentration under manure-fertilized maize. Plant Soil 282: 343-360.

Søgaard HT, Sommer SG, Hutchings NJ, Huijsmans JFM, Bussink DW, Nicholson F, 2002. Ammonia volatilization from field-applied animal slurry: the ALFAM model. Atmos Environ 36: 3309-3319.

Soto F, Gallardo M, Giménez C, Peña-Fleitas T, Thompson RB, 2014. Simulation of tomato growth, water and N dynamics using the EU-Rotate_N model in Mediterranean greenhouses with drip irrigation and fertigation. Agric Water Manag 132: 46-59.

Stella T, Frasso N, Negrini G, Bregaglio S, Cappelli G, Acutis M, Confalonieri R, 2014. Model simplification and development via reuse, sensitivity analysis and composition: a case study in crop modelling. Environ Modell Softw 59: 44-58.

Suarez-Rey EM, Romero-Gámez M, Giménez C, Thompson RB, Gallardo M, 2016. Use of EU-Rotate_N and CropSyst models to predict yield, growth and water and N dynamics of fertigated leafy vegetables in a Mediterranean climate and to determine N fertilizer requirements. Agric Syst 149: 150-164.

Suárez-Rey EM, Gallardo M, Romero-Gámez M, Giménez C, Rueda FJ, 2019. Sensitivity and uncertainty analysis in agro-hydrological modelling of drip fertigated lettuce crops under Mediterranean conditions. Comput Electron Agr 162: 630-650.

Sun Y, Hu K, Fan Z, Wei Y, Lin S, Wang J, 2013. Simulating the fate of nitrogen and optimizing water and nitrogen management of greenhouse tomato in North China using the EU-Rotate_N model. Agric Water Manag 128: 72-84.

Wagenet RJ, Hutson JL, 1989. LEACHM: Leaching Estimation and Chemistry Model: A process based model of water and solute movement, transformations, plant uptake and chemical reactions in the unsaturated zone, vers. 2. Water Resour Inst, Cornell Univ, Ithaca, NY.