*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.*

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

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

*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.*

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

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

_{FC}(volumetric water content at field capacity);

_{SAT}(volumetric water content at saturation).

Funding agencies/Institutions | Project / Grant |

Spanish Ministerio de Economía y Competitividad, INIA | RTA 2011-00136-C04-01 |

The improvement in some agricultural practices, as irrigation and fertilization, to reduce some of the environmental problems they produce, can be undertaken using crop numerical models (

Consequently, the estimation of the uncertain parameters from experimental data is an important task, as model predictions depend on the accuracy of the parameter estimates (

A local sensitivity analysis is based mainly on the estimation of the local derivative of the model output regarding the model input parameters. This derivative 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 information only at the point where the derivatives are calculated 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 Hypercube-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 (

Experimental data were obtained from field experiments 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 ^{th}
, 2012, with a density of plantation of 23,420 plants ha
^{-1}
. The crop was harvested on February 20
^{th}
, 2013 and yield was 41.9 t ha
^{-1}
. The crop residues were incorporated to soil on March 21
^{st}
2013. In Paterna-2 plot transplantation was on September 16
^{th}
, 2013, with a plant density of 21,140 plants ha
^{-1}
. The cauliflower was harvested on February 15
^{th}
, 2014 and yield was 33.7 t ha
^{-1}
. In this plot, the crop residues were incorporated to soil on March 4
^{th}
2014. The cultivar used in both plots was Triomphant, which is a long cycle type (about 150 days).

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 29
^{th}
, 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 ammonium) 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 is a process-based, one-dimensional model that simulates water and solute movement, and related chemical and biological processes, in the unsaturated soil (

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 (
^{3}
cm
^{-3}
), matric potential (
^{-1}
) are required in the water flow model, which are obtained using Campbell's equations (

where
_{SAT}
is the volumetric water content at saturation (cm
^{3}
cm
^{-3}
),
_{s}
is the hydraulic conductivity at saturation (mm d
^{-1}
), and

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

LEACHM model contemplates different processes related to the N dynamics in the soil (

being
_{i}
the carbon content (mg kg
^{-1}
) and
_{mi}
the mineralization rate (d
^{-1}
) in each compartment. The relation among the three carbon pools is given by the efficiency factor
_{e}
and the humification fraction
_{h}
, (

trification and volatilization. Three organic compartments (litter, humus and manure) and three mineral compartments (urea, nitrate and ammonium) define

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

being
_{nit}
the nitrification rate (d
^{-1}
) and
_{max}
the maximum relation NO
_{3}
^{-}
/NH
_{4}
^{+}
. 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
^{-1}
.

EU-Rotate_N is a decision support system for N management in crop rotations developed under the framework of a European Union project (

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 (

Nitrogen mineralization from soil organic matter is based on the approach implemented in the DAISY model (

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, A
_{0}
and B
_{0}
(

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
^{m}
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

An orthogonal factorial design is used in the simulations, and an ANOVA analysis is performed separately on each output variable (

where
_{i}
is the sum of the squares associated with the main effect of parameter m
_{i}
,
_{ij}
is the interaction between the parameters m
_{i}
_{j}
, and
_{1...m}
is the interaction between up to the

and the total sensitivity indices are defined as:

taking into account the interaction of factor

To reduce the dimension of the problem, a principal components analysis (PCA) of matrix
_{c}
is the matrix
^{T}
_{c}
_{c}
. If
_{1}
_{c}
are the eigenvalues of
^{T}
_{c}
_{c}
, and

_{c}

When an orthogonal factorial design is used, ANOVA-based sensitivity analysis can be applied to each principal component
_{l}
obtaining first order sensitivity indices and total sensitivity indices for each one of the principal components considered,
_{c}
, being
_{c}
_{i}
, is defined for each parameter involved in the sensitivity analysis.
_{i}
is equal to the weighted sum of sensitivity indices of the
_{c}
principal components with weights proportional to the inertia associated with the
_{c}
principal components. The computation of the generalized sensitivity index is implemented in R language in the library

The simulation period for both models has been from September 3
^{rd}
of 2012 until May 30
^{th}
of 2013 for Paterna-1 plot and from September 10
^{th}
of 2013 until May 5
^{th}
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
_{s}
values was estimated using the expressions proposed in

For the EU-Rotate_N model, the volumetric water content at field capacity (
_{FC}
), at saturation (
_{SAT}
) and the drainage coefficient (

Once the most important parameters related with the water dynamics were obtained by determining the
_{i}
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
_{vij}
(cm
^{3}
cm
^{-3}
) is the
^{3}
cm
^{-3}
) is the

In the LEACHM model, most of the selected parameters, shown in _{4}
^{+}
and molecular diffusion coefficient have also been taken into account.

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 (

With respect to the parameters related with cauliflower growth, the target dry matter value was the maximum value measured in the field (14.5 t ha
^{-1}
). The N critical curve coefficients A
_{0}
and B
_{0}
were modified with respect to the default values for cauliflower in the
_{0}
was 3.26 and for B
_{0}
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 (

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
_{vij}
is the
_{3}
^{-}
content in the soil in the
^{-1}
) and
_{vij}
is the
_{3}
^{-}
content in the
^{-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 (

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 _{FC}
and
_{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 parameters in the optimization process did not improve the simulations.

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 _{w}
, obtained was slightly lower in EU-Rotate_N (0.221) than in LEACHM (0.258).

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 (_{FC}
was the most influencing parameter in water drainage and, therefore, the soil water content is mainly determined by this value.

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

GSI indices obtained for N parameters calibration in both models are shown in

Again, the parameters with the four largest values of the total GSI were chosen for calibration of the N module. _{N}
. The other variables were set to the mean value of their range of variation, shown in _{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

45 cm). RRMSE values for the whole simulation were high (0.7 in both models), and RMSE was about 50 kg N-NO
_{3}
^{-}
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

The accumulated nitrate leaching simulated and measured at 45 cm depth is shown in _{3}
^{-}
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.

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.

The RMSE ranged from 50 kg N-NO
_{3}
^{-}
ha
^{-1}
for LEACHM model to 83 kg N-NO
_{3}
^{-}
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

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

With respect to
_{s}
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
_{s}
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

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
_{4}
^{+}
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-NO
_{3}
^{-}
ha
^{−1}
. The model underestimated soil nitrate content, mainly after the incorporation of crop residues.

EU-Rotate_N model uses a simple approach to calculate soil water redistribution in the soil profile and requires few soil hydraulic parameters. GSI sensitivity analysis showed that
_{FC}
and

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

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 _{3}
^{-}
ha
^{−1}
. Nitrate leaching simulation was better than that obtained with the LEACHM model, but still provided an underestimation.

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-NO
_{3}
^{-}
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 (

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 incorporation to soil. To avoid overparameterization, a global sensitivity analysis adapted to discrete time model outputs was used to determine the most influential 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.