Global sensitivity analysis for choosing the main soil parameters of a crop model to be determined ()
1. INTRODUCTION
Dynamic crop models are very useful to predict the behavior of crops in their environment and are widely used in a lot of agro-environmental work such as crop monitoring, yield prediction or decision making for cultural practices [1-3]. These models usually have many parameters and their spatial application for agro-environmental predictions is difficult without a good knowledge of these parameters [4-6].
The crop model parameters can be divided in three groups: those related to the agricultural techniques, those related to the genotype of the crop and those related to the soil properties. Generally, agricultural techniques are quite easy to know as they are those used by the farmer. Crop parameters can be determined from literature, or estimated from experimental work or calibrated on a large database [7-9]. The knowledge of the soil parameters is an important issue because the spatial variability of the crop model simulations depends for a large part on the soil parameter values [10] and predictions obtained with the model are not reliable when inaccurate parameter values are used. This knowledge may be especially difficult to acquire because parameter values can greatly vary in space [11,12]. The use of existing soil maps and associated pedotransfer functions can be considered where accurate soil map are available [13]; but in many cases, the spatial accuracy of the map is too limited for accurate applications such as for example precision agriculture [14]. In those cases, these parameters should be determined in another way. Measurements can be made directly with soil sampling analysis at different locations of the study area or indirectly by using electrical geophysical measurements [15,16]. Whatever the technique of measurement used, it is submitted to practical limitations and to time and financial constraints. Another way of gathering quite accurate values on soil parameters consists in estimating them through an inverse modeling approach using a crop model and observations on the crop state variables [11,17,18]. However, the soil parameters may not have the same contribution to the performance of the crop model and do not require the same precision of determination for a given objective: some of them deserve an accurate determination while the others can be fixed at nominal values [19,20]. Considering this aspect, the practical limitations of soil parameter measurements, as well as time and financial constraints should be reduced by considering only a subset of the crop model soil parameters depending on the objective and configuration of the study.
The combination of uncertainty analysis and sensitivity analysis techniques should help in identifying these key parameters. The objective of sensitivity analysis is to study how the variation of selected outputs of a model can be apportioned to different sources of variation [21]. In particular, sensitivity analysis methods can be used to rank uncertain input factors with respect to their effects on the model output variables by calculating quantitative or qualitative indices. Nevertheless, the fact that some factors are detected as important for a given output variable on the basis of sensitivity analysis results is not sufficient to decide that the uncertainties on these factors should be reduced. Indeed, if the variation of the considered output variable induced by the uncertainties on the factors is low, the results of sensitivity analysis on this output variable should not be taken into account. The description and quantification of these variations is the objective of uncertainty analysis.
Some authors [19,22-26] used uncertainty analysis techniques to quantify the uncertainties of a selection of crop models output variables generated by uncertainties on some selections of input parameters. Others authors [27-32] used global sensitivity analysis to evaluate the contribution of the parameters to the variance of the model output variables. In this study we propose a helpful combination of these techniques to identify soil parameters that need particular accuracy for simulating a set of given output variables of interest in spite of the financial and practical interests of such a study.
A variance-based sensitivity analysis method is used in order to rank the soil parameters relatively to their importance on some selected outputs of the crop model STICS (Simulateur multidisciplinaire pour les Cultures Standard) [33] and to select those which deserve an accurate determination by considering also the coefficient of variation of the outputs, that is the variation of the outputs compared to their magnitude. We considered 13 soil parameters and their effects on 5 dynamic output variables of the STICS crop model, at different phenological stages, which are involved in decision making for crop management. Two different crops (winter wheat and sugar beet) growing on different seasons are considered in order to illustrate the impact of soil properties on crop growth. Each crop considered is simulated under different pedological conditions and weather.
2. METHODS
2.1. The STICS Model
The STICS model [33,34] is a nonlinear dynamic crop model simulating various crops. For a given crop, STICS takes into account the weather, the type of soil and the cropping techniques used, and simulates the carbon, water and nitrogen balances of the crop-soil system on a daily time-scale. In this study, winter wheat and sugar beet crops are simulated. The crop is essentially characterized by its above-ground biomass carbon and nitrogen, and leaf area index. The soil is considered as a series of layers where the transfer of water and nitrate is described by a reservoir-type analogy. The main outputs are agronomic variables (yield, grain protein content for wheat) as well as environmental variables (water and nitrate leaching).
The STICS model includes more than 200 parameters. The global sensitivity analysis described in this study only concerns the soil parameters. The values of the parameters related to the crop have been determined either from literature, from experimental works conducted on specific processes included in the model (e.g. mineralization rate, critical nitrogen dilution curve etc.) or from a calibration based on a large experimental database [35]. Cropping techniques and soil parameters ranges are described in Section 2.5.
2.2. The Soil Parameters
Among the available options for simulating the soil system, the simplest was chosen in this study, by considering only the transfers in the microporosity and ignoring those in the macroporosity, the cracks, pebbles, and processes like capillary rise and nitrification. We then considered the soil as a succession of two horizontal layers, the top layer having a thickness fixed at 30 cm.
These different hypotheses made on the soil description lead to consider a set of 13 soil parameters, defined in Table 1. They refer to permanent characteristics and initial conditions. Among the permanent characteristics, clay and organic nitrogen content of the top layer are
Table 1. Definition of the 13 soil parameters and their ranges of variation.
involved mainly in organic matter decomposition processes. Water content at field capacity of both layers affects the water (and nitrogen) movements and storage in the soil reservoir and the thickness of the second layer defines the volume of the reservoir. The initial conditions correspond to the water and nitrogen content, Hinit and NO3init, at the beginning of the simulation, in this case the sowing date.
2.3. Model Output
In this study, the STICS output variables of interest are:
1) The amount of nitrogen absorbed by the plant (QN) and the leaf area index (LAI) at two (for sugar beet) or three (for wheat) different key stages during the season, which are important variables for making a diagnosis on crop growth2) The yield, and the mineral nitrogen content in the soil at harvest (for both crops) plus the grain protein content (for wheat), which are of particular interest for decision making, especially for monitoring nitrogen fertilization.
The different stages of interest and the corresponding variables are displayed for each crop on Table 2. For the wheat, the three key stages concern the maximum leaf growth rate-beginning of stream elongation-(AMF), the maximum leaf areaor booting-(LAX) and the flowering (FLO). For the sugar beet, the two key stages concern the maximum leaf growth rate (AMF) and the maximum leaf area (Summer).
2.4. Sensitivity and Uncertainty Analysis
Among the available methods of sensitivity analysis, variance-based methods are well adapted for non-linear models that need less than 1 minute for a simulation [36]. These methods are widely used in different contexts [27,28,30-32]. Their principle is to evaluate the contribution of the given uncertain factors to the variance of the model output variables selected. We will describe in this section the sensitivity indices that can be computed with these methods and the EFAST variance-based method we have used in this study to compute these indices. Uncertainty analysis is performed here by computing the coefficient of variation of the output variables considered from the simulations realized for the sensitivity analysis.
2.4.1. Sensitivity Indices and Coefficient of Variation
We note further Y an output variable of STICS. Y will represent in turn LAI and QN computed at the different phenological stages and the variables computed at harvest. The total variance of Y, V(Y), is partitioned as follows [37]:
(1)
where V(Y) is the total variance of the output variable Y induced by the 13 soil parameters q,
Table 2. Definition of the variables and the stages of interest.
measures the main effect of the parameter qi, i = 1, ···, 13, and the other terms measure the interaction effects. Decomposition (1) is used to derive two types of sensitivity indices defined by:
(2)
(3)
where V−i is the sum of all the variance terms that do not include the index i.
Si is the first-order sensitivity index of the ith parameter. It computes the fraction of Y variance explained by the uncertainty of parameter qi and represents the main effect of this parameter on the output variable Y. STi is the total sensitivity index of the ith parameter and is the sum of all effects (first and higher order) involving the parameter qi. Si and STi are both in the range (0, 1), low values indicating negligible effects, and values close to 1 huge effects. STi takes into account both Si and the interactions between the ith parameter and the 12 other parameters, interactions which can therefore be assessed by the difference between STi and Si. The interaction terms of a set of parameters represent the fraction of V(Y) induced by these parameters but that is not explained by the sum of their main effects. The two sensitivity indices Si and STi are equal if the effect of the ith parameter on the model output is independent of the values of the other parameters: in this case, there is no interaction between this parameter and the others and the model is said to be additive with respect to qi. Selecting the parameters that have a negligible effect and that can thus be fixed to nominal values is called factor fixing in the literature [38]. Total effects must be considered in this case. Indeed, a factor that has a small main effect but a medium to high total effect cannot be considered as negligible: its effect depends on the value of other uncertain factors and can be important in some cases.
The coefficient of variation of the output variable Y can be calculated by:
(4)
where, s(Y) is the standard deviation of the output variable Y and is the mean of Y induced by the 13 soil parameters q.
2.4.2. Extended FAST
Sobol’s method and Fourier Amplitude Sensitivity Test (FAST) are two of the most widely used methods to compute Si and STi [11]. We have chosen here to use the extended FAST (EFAST) method, which has been proved, in several studies [30,39,40], to be more efficient in terms of number of model evaluations than Sobol’s method. The main difficulty in evaluating the first-order and total sensitivity indices is that they require the computation of high dimensional integrals. The EFAST algorithm performs a judicious deterministic sampling to explore the parameter space which makes it possible to reduce these integrals to one-dimensional ones using Fourier decompositions. The reader interested in a detailed description of EFAST can refer to [40].
We have implemented the EFAST method in the Matlab® software, as well as a specific tool for computing and easily handling numerous STICS simulations. The uncertainties considered for the soil parameters are described in the next section. A preliminary study of the convergence of the sensitivity indices allowed us to set the number of simulations per parameter to 5000, leading to a total number of model runs of 13 × 5000 = 65,000 to compute the main and total effects for all output variables and parameters considered here. One run of the STICS model takes about 1s with a Pentium 4, 2.9 GHz processor.
2.5. Data
In this study, we have considered two crops: winter wheat and sugar beet. This allows to illustrate the difference of sensitivities of different crops to the soil properties. For the same reason, each crop is simulated for two different weathers and two different types of soil depth. The weather data were obtained from the meteorological station of Roupy (49.48˚N, 3.11˚E). The first set of data is chosen to characterize a dry weather (1975-1976) and the second set is chosen to characterize a wet weather (1990-1991). Table 3 shows the amount of rainfall calculated for each season and weather data set. The wheat crop simulated in this study is sown on October 30th while the sugar beet crop is sown on March 30th. The amount of fertilizer provided on wheat varies between 200 kg (shallow soil) and 240 kg (deep soil), while the amount provided on sugar beet varies between 150 kg (shallow soil) and 200 kg (deep soil).
The range of parameter values considered in this study correspond to the soil description of the precision agriculture experimental site in northern France near Laon, Picardie (Chambry 49.35˚N, 3.37˚E) [41]. In this study, the uncertainties of these 13 soil parameters are observed in the literature (for parameters related to albedo, evapotranspiration or drip rainfall) or measured in the experimental site (for the other parameters), and their ranges of variation are displayed on Table 1. Concerning the parameter NO3init two ranges of variation are considered, depending on the crop cultivated just before the one considered: in this study, the wheat is cultivated after sugar beet and the sugar beet is cultivated after a bare soil. The different previous crops used determine the quantity of nitrogen NO3init at the beginning of the corresponding crop season. The two different types of soil depth are defined by their ranges of variation (Table 1) and correspond to a shallow soil and a deep soil. The uncertainties considered in the global sensitivity analysis for the soil parameters are assumed independent and follow uniform distributions. The ranges of variation of the distributions are given in Table 1.
3. RESULTS AND DISCUSSION
Only the main results of the study are presented here for the sake of clarity. These results concern: 1) wheat crop simulated with dry, then wet weather and a shallow soil and 2) sugar beet crop simulated with dry, then wet weather and a deep soil.
3.1. Global Sensitivity Analysis
Figure 1 shows the sensitivity indices calculated for the 13 soil parameters and for each output variable of the wheat crop simulated with a dry weather and a shallow soil. For the early stage the initial water content is dominant because in the considered weather, the rainfall is light in autumn when the wheat is sown (see Table 3): at the stage AMF (Figure 1(a)), Hinit is the only parameter
Table 3. Amount of rainfall (in mm) calculated for each season and weather.
contributing (for more than 90%) to the output variance for both variables LAI and QN. In the later stages, the effects of parameters evolve with the soil volume explored by the roots (first layer, then second one), up to the flowering stage where the root development is maximum: at the stage LAX (Figure 1(b)) and FLO (Figure 1(c)), the effect of Hinit disappears and that of HCC(1) and epc(2) increase, with a dominant position of epc(2). At Harvest (Figure 1(d)), the variables are much sensitive to epc(2) followed by HCC(1), HCC(2) and Hinit for the three variables and albedo for the variable protein of the grain. In those conditions of dry weather and shallow soil, the parameters related to water availability (epc(2), HCC(1), HCC(2) and Hinit) are the main parameters contributing to the variance of the outputs. Those concerned by the turnover of organic nitrogen in the soil are not concerned, because the water stress is the dominant limiting factor and also because the mineralization processes are reduced by dry conditions.
When considering wet conditions (Figure 2), the water stress is not so much a limiting factor: maximum LAI is equal to 3.61 in average, whereas it is equal to 2.57 in dry conditions (see Table 4). The roots grow more rapidly at the beginning of the season and the size of the soil reservoir (via the parameter epc(2)) is important since the AMF stage: the depth of roots is equal to 55.84 in average (3 months after the sowing), whereas it is equal to 45.62 in dry conditions (see Table 4). Moreover, in these conditions, the mineralization of the soil organic matter is increased and the effects of the concerned parameters argi and Norg do so: the cumulative mineral nitrogen arising from humus is equal to 23.95 in average (at stage LAX), whereas it is equal to 18.09 in dry conditions (see Table 4). This does not seem to influence the effects of the different parameters on LAI at stage AMF since they are very similar to these obtained with the dry weather. On the contrary, the sensitivities of the variable QN to the different parameters are very different from the ones obtained with a dry weather: there is no contribution of Hinit but epc(2), HCC(1) and parameters involved in the mineralization process (argi, Norg and NO3init) contribute to the output variance.
This is also the case for both LAI and QN on later stages, with an increasing dominancy of epc(2). At Harvest (Figure 2(d)), the variables are sensitive to the parameters epc(2), HCC(1) and HCC(2) with still a slight sensitivity to argi and Norg for the soil mineral nitrogen
Figure 1. Sensitivity indices of the 13 soil parameters for each model output of the wheat crop simulated with a dry weather and a shallow soil. The outputs (a) correspond to LAI and QN at stage AMF; (b) correspond to LAI and QN at stage LAX; (c) correspond to LAI and QN at stage FLO and (d) correspond to Yld, Prot and Nit at Harvest. First-order indices are in black and interactions in white.
content. The main difference between these results and those presented in Figure 1 lies in the sensitivity to parameters involved in the mineralization process (especially argi and Norg).
When considering wet conditions (Figure 2), the water stress is not so much a limiting factor: maximum LAI is equal to 3.61 in average, whereas it is equal to 2.57 in dry conditions (see Table 4). The roots grow more rapidly at the beginning of the season and the size of the soil reservoir (via the parameter epc(2)) is important since the AMF stage: the depth of roots is equal to 55.84 in average (3 months after the sowing), whereas it is equal to 45.62 in dry conditions (see Table 4). Moreover, in these conditions, the mineralization of the soil organic matter is increased and the effects of the concerned parameters argi and Norg do so: the cumulative mineral nitrogen arising from humus is equal to 23.95 in average (at stage LAX), whereas it is equal to 18.09 in dry conditions (see Table 4). This does not seem to influence the effects of the different parameters on LAI at stage AMF since they are very similar to these obtained with the dry weather. On the contrary, the sensitivities of the variable QN to the different parameters are very different from the ones obtained with a dry weather: there is no contribution of Hinit but epc(2), HCC(1) and parameters involved in the mineralization process (argi, Norg and NO3init) significantly contributes to the variance of this variable. This is also the case for both LAI and QN on later stages, with an increasing dominancy of epc(2). At Harvest (Figure 2(d)), the variables are sensitive to the parameters epc(2), HCC(1) and HCC(2) with still a slight sensitivity to argi and Norg for the soil mineral nitrogen content. The main difference between these results and those presented in Figure 1 lies in the sensitivity to parameters involved in the mineralization process (especially argi and Norg).
Figure 2. Sensitivity indices of the 13 soil parameters for each model output of the wheat crop simulated with a wet weather and a shallow soil. The outputs (a) correspond to LAI and QN at stage AMF; (b) correspond to LAI and QN at stage LAX; (c) correspond to LAI and QN at stage FLO and (d) correspond to Yld, Prot and Nit at Harvest. First-order indices are in black and interactions in white.
Table 4. Ranges of some output variables uncertainties generated by the uncertainties on the soil parameters. The output concerns the value of maximum LAI, the cumulative mineral nitrogen arising from humus Qminh (calculated at the stage LAX or Summer) and the depth of roots Zrac (calculated 3 months after the sowing date).
Figure 3 shows the sensitivity indices calculated for the 13 soil parameters and for each output variable of the sugar beet crop simulated with a deep soil and a dry weather. In this case, the crop grows mainly in summer where it experiences a severe water stress, leading to a value of maximum LAI equal to 1.42 in average (see Table 4). The depth of the second layer (parameter epc(2)) does not have any importance here. This is also the case for wheat crop with a deep soil (results not shown here). Indeed, as the root growth is no longer limited by the