A Comparison of Moment-Independent and Variance-Based Global Sensitivity Analysis Approaches for Wheat Yield Estimation with the Aquacrop-OS Model

: The present work reports the global sensitivity analysis of the Aquacrop Open Source (AOS) model, which is the open-source version of the original Aquacrop model developed by the Food and Agriculture Organization (FAO). Analysis for identifying the most inﬂuential parameters was based on di ﬀ erent strategies of global SA, density-based and variance-based, for the wheat crop in two di ﬀ erent geographical locations and climates. The main objectives were to distinguish the model’s inﬂuential and non-inﬂuential parameters and to examine the yield output sensitivity. We compared two di ﬀ erent methods of global sensitivity analysis: the most commonly used variance-based method, EFAST, and the moment independent density-based PAWN method developed in recent years. We have also identiﬁed non-inﬂuential parameters using Morris screening method, so to provide an idea of the use of non-inﬂuential parameters with a dummy parameter approach. For both the study areas (located in Italy and in China) and climates, a similar set of inﬂuential parameters was found, although with varying sensitivity. When compared with di ﬀ erent probability distribution functions, the probability distribution function of yield was found to be best approximated by a Generalized Extreme Values distribution with Kolmogorov–Smirnov statistic of 0.030 and lowest Anderson–Darling statistic of 0.164, as compared to normal distribution function with Kolmogorov–Smirnov statistic of 0.122 and Anderson–Darling statistic of 4.099. This indicates that yield output is not normally distributed but has a rather skewed distribution function. In this case, a variance-based approach was not the best choice, and the density-based method performed better. The dummy parameter approach avoids to use a threshold as it is a subjective question; it advances the approach to setting up a threshold and gives an optimal way to set up a threshold and use it to distinguish between inﬂuential and non-inﬂuential parameters. The highly sensitive parameters to crop yield were speciﬁcally canopy and phenological development parameters, parameters that govern biomass / yield production and temperature stress parameters rather than root development and water stress parameters.


Introduction
Crop growth simulation models are important tools for ecology, agriculture and environmental research and management. Such models are used to simulate the behavior of soil-crop systems in response to climate and agriculture management [1][2][3]. Agricultural processes, the interaction of the crop with the environment and soil over time, are represented by mathematical equations based on theory and empirical research. These equations and parameterizations inevitably entail simplifications of reality, which leads to uncertainty and inaccuracy of output variables [4]. Regardless of the model used, the identification of the parameters and input variables that most affect its output is a fundamental problem for most environmental applications whenever large uncertainty on their values is expected, such as in regional-scale applications [5][6][7].
In general, crop models contain many parameters and input variables, of which some cannot be measured directly, but can only be inferred by a calibration of the responses observed from the model [8,9]. Manual, one-at-a-time calibration is a very tedious and time-consuming task, which is not feasible for complex models [10]. Numerical optimization or Monte Carlo methods of calibration are widely used [11], but they are prone to errors, e.g., when local minima are present and their efficiency is reduced as the number of parameters increases [12]. Therefore, it is often not possible to include all the parameters in the model calibration [13,14]. In this view, Sensitivity Analysis (SA) is the set of methods that serve the purpose of selecting a limited subset of parameters or input variables, excluding those that are non-influential and have a low impact on the model output. The aim of the SA is to determine how sensitive the outputs of a crop model are with respect to the elements of the model, which are subject to uncertainty or variability [15]. The main purposes of SA are Factor Prioritization (FP) or ranking, Factor Fixing (FF) or screening and Factor Mapping (FM) [16]. This work focuses on ranking and screening.
There are two different strategies for SA: local and global. Local SA focuses on the sensitivity of the factors at particular points in the feature space [17], whereas global SA methods assess the sensitivity in the entire feature space and include higher-order, variance-based methods, e.g., Sobol or Fourier methods [18,19]. An extensive analysis of SA methods on crop models has been performed by different authors [16,20], which shows the inadequacy of local SA methods for crop models due to the complex nature of such models and to the inability to determine the combined effects of the parameters. For these reasons, global SA methods are currently preferred for the SA of the crop models.
Many different global SA methods have been developed [21][22][23][24][25]. Among all these methods, variance-based methods, e.g., Sobol or Fourier, are probably the most widely used for global SA. In general, these methods measure sensitivity to an uncertain parameter or input variable by assessing the contribution of that input to the total variance of the model output variable. The advantage of the variance-based methods is their ability to quantify the contribution of the individual parameters and also quantify the contribution due to interacting effects of the parameters, independently from assumptions on the form of the input-output relation (e.g., linearity and additivity). Variance-based sensitivity indices are also easy to interpret, as they represent the fraction of the output variance caused by the variation of input [26].
Variance-based methods are moment-dependent methods and consider the second-order moment, i.e., variance as a measure of the output uncertainty. It is assumed that this moment is sufficient to describe the output variability [26]. However, it has been recognized that the variance does not adequately represent output uncertainty when the model output is not represented by variance, e.g., output that is highly skewed or multi-modal [25,27,28]. To overcome this limitation (moment-dependent methods), moment-independent global SA measures have been developed [24,25,27]. These methods are also known as density-based methods. They use the entire output distribution to fully characterize the output uncertainty and to quantify the relative influence of the uncertain parameters. The main advantage of these methods, as compared to the variance-based methods, is that they do not use any specific moment of the output distribution to measure the output variability and, therefore, are applicable regardless of its shape (e.g., multi-modal or highly skewed). The moment-independent density-based Global Sensitivity Analysis (GSA) method, PAWN, was introduced by [25]. It measures sensitivity based on the difference between unconditional output distribution and conditional output distribution, which are obtained simultaneously when all the parameters are free to vary and when one of the parameters is fixed. The difference between PAWN and other density-based methods, e.g., the entropy-based method [27] and the δ-sensitivity measure [24], is that the output distribution in PAWN is characterized by the cumulative distribution function (CDF), whereas in other methods the output distribution is characterized by the probability density function (PDF). The approximation of CDFs by using empirical distributions of the data sample is much easier than the approximation of their PDFs. This facilitates the analysis of the robustness and the convergence of the estimated sensitivity indices [25].
The PAWN method was initially tested with a simple conceptual hydrological model with only five parameters. To further investigate its efficiency, it was applied to the complex environmental model -Soil and Water Assessment Tool (SWAT) [29] with 26 parameters and compared with variance-based method Sobol.
The Aquacrop model [30] is a crop water productivity model developed by the Food and Agriculture Organization (FAO) of the United Nations specifically for the purpose of assessing crop response to water and increasingly used by scientists and agronomists [31][32][33][34][35]. The Aquacrop Open Source (AOS) version was developed by [31] and made available as an open-source Matlab version. In principle, this open-source version is the same as the original developed FAO version of the Aquacrop model. Recently, a version of AOS, v6.0a, has been released. It implements all calculation procedures performed by the original FAO Aquacrop v6.0a, with the exception of soil fertility stress and soil salinity stress, as these are yet to be tested extensively for different crop types and environments worldwide [36]. The soil water balance, stress effects, runoff and plant growth are controlled by a large number of parameters (around 100). Even when some of the parameters can be fixed a priori, calibration of the AOS remains quite challenging, given the relatively large number of parameters (54 in our case) that are typically left to be varied simultaneously. Therefore, SA is often applied prior to the calibration process to identify the most influential parameters and the non-influential ones [7,37]. To avoid any confusion with the model, we will hereafter refer to the FAO developed model as Aquacrop and its open-source Matlab version (The MathWorks Inc., Natick, MA, USA) as AOS that is used in this study.
There have been studies of sensitivity analyses with the FAO developed Aquacrop model [7,37,38], but none yet on the open-source version AOS, which nevertheless lends itself more easily to the integration into regional-scale data assimilation schemes than the standard version. The studies of Aquacrop model SA was based on the arbitrarily chosen thresholds for yield output. According to the authors' best knowledge this is the first time the dummy parameter approach [29] is being used to identify the most influential parameters with the recently released AOS model (v6.0).
The main objectives of this study are to (i) assess and compare the application of the Aquacrop model approach with the Extensive Fourier Amplitude Sensitivity Test (EFAST) method and the density-based PAWN method with a high number of parameters (soil and crop) in terms of parameter ranking results; (ii) perform a sensitivity analysis of the AOS model, provided as an open-source Matlab version, that in principle is similar to the original FAO developed Aquacrop model, but differs in some parameters; (iii) compare the sensitivity indices (SI) of the EFAST and density-based PAWN methods for two different geographical locations and climates. As both methods have different backgrounds and rationale, the comparison cannot be made on the basis of the Sensitivity Index (SI); therefore, the comparison is made in terms of parameter ranking as a measure of the effectiveness with multiple years as a measure of robustness. To assess the parameter ranking, we have used the dummy parameter approach proposed by [29], while the SI of the dummy parameter has used as a threshold to identify non-influential parameters.

Sensitivity Analysis
First, the Morris method was used to distinguish between influential and non-influential parameters of the AOS model. The Morris method is designed for an initial screening or the elementary effects that the variation of input factors (parameters and input variables) produces on model outputs. Secondly, the variance-based EFAST method and finally the density-based PAWN method were applied to all the chosen parameters considered for sensitivity analysis (see Table 1), including a dummy parameter to find the higher-order effects and the combined effect of different parameters.

Morris-Based Elementary Effects
The Morris method is designed for the screening of the most influential parameters [39]. It is considered more a qualitative rather than a quantitative method as it does not quantify the sensitivity of a parameter individually, rather it represents the parameter sensitivity in terms of mean (m*) and sigma (s). High m* indicates higher sensitivity and higher s indicates interacting effects with other parameters. The method determines whether the effects are negligible, additive, non-linear or the factors are having interacting effects with other parameters [40]. To sample input space, the experimental design is structured in a group of points, called trajectories. In this work, the number of trajectories has been set to 20, which is a value between the minimum suggested by [40] (i.e., 10) and the value used by [37] (i.e., 25) in order to avoid an excessively high computational cost. It computes two indices for each input: the mean of the effects, which measures the total effect of input over the output, and the standard deviation of the effects, which measures the degree of interactions with the other inputs. The input factor with a mean below a threshold value can be considered as having negligible effects on the model output, thus it can be eliminated from further processing by keeping its value as fixed. Whereas, other authors have defined an arbitrary threshold [7,21]. In this work, we chose to identify influential parameters by comparing them with a dummy parameter used as a threshold [29]. All the parameters having a mean value greater than that of the dummy parameter are considered as influential.

EFAST Method
The EFAST method [4] combines the FAST [41] and the Sobol methods [42]. It calculates sensitivity indices using the total variance of the output of the model and the conditional variances depending on the parameters. The interaction among the factors can be quantified by calculating the main SI (Si) and the index of total sensitivity (Sti), i.e., the sum of the main effect plus the interaction between the variation of the terms of the parameters to all orders. Si and Sti range between 0 and 1, with higher values indicating more important effects [7].
In this study, the input parameters plus the dummy parameter were implemented using the EFAST method to identify influential parameters and to rank them. The number of sample points considered for the chosen 55 parameters, including a dummy variable, is 10,000. The required number of runs are thus 55 × 10,000. The number of higher harmonics M s is 4, and the number of search curves N cs was set to 1.

Density-Based PAWN Method
The PAWN method is a moment-independent global sensitivity analysis method that takes into account the entire model output distribution for the quantification of the sensitivity of the parameters [25], in comparison to the other variance-based methods that consider only the moments, e.g., the variance of the model output. In general, density-based sensitivity indices measure the sensitivity to parameter x i by the distance between unconditional PDF and conditional PDF of the model output [24,27]. Practically, PDFs are usually unknown and must be approximated using a data sample. Authors in [25] pointed out the difficulties and limitations of deriving empirical PDFs and suggested to use CDFs instead of PDFs, as it is much easier than using the approximation of PDFs.
In the PAWN method, sensitivity is measured by estimating the variations that are induced in the output distribution when removing the uncertainty about one or more inputs. The sensitivity to the input parameter is measured using the distance between the unconditional probability distribution and the conditional distributions that are obtained by keeping the parameter fixed using the Kolmogorov-Smirnov (KS) statistics [43,44].
The unconditional cumulative distribution function of the model output y is denoted by F y (y), and the conditional cumulative distribution is denoted by F y|xi (y), where x i is a parameter under consideration.
If F y (y) and F y|xi (y) coincide, it can be concluded that the parameter is non-influential. If the distance between F y (y) and F y|xi (y) is increased, the influence of the parameter increases as well. As a measure between unconditional and conditional CDFs, the KS statistics were used (Equation (1) [25]) As KS depends on the value at which we fix x i , the PAWN index T i considers a statistic (e.g., the maximum or the median) over all possible values of x i , i.e., By definition, T i varies between 0 and 1; the lower the value of T i, the less influential x i , and if T i = 0, then x i has no influence on y.
The total number of model evaluations to compute the sensitivity indices are N u + n × N c × M. M is the number of parameters, N u and N c are the numbers of unconditional random samples and conditional samples and n is the number of sampling points. In this work, we used N u = 300, N c = 200, n = 10 and M = 55.
The density-based PAWN method estimates sensitivity by implementing unconditional and conditional cumulative distribution functions. In the case of unconditional distribution functions, all the parameters are free to vary within their defined range, and in the case of the conditional distribution function, the parameter into consideration is fixed to the nominal value. The parameter with the low sensitivity or called PAWN index has the same overlapping unconditional and conditional distribution functions.
In the present study, the Morris and PAWN methods were implemented with the Matlab (The MathWorks Inc., Natick, MA, USA)-based SAFE toolbox [45] and EFAST based on the SA toolbox Eikos [46]. To figure out the best PDF fit for the yield output and different statistical distributions, tests including KS statistics and Anderson-Darling statistics were carried out with EasyFit Evaluation (Ver. 5.6, http://www.mathwave.com/).

Identifying Non-Influential Parameters by Using a Dummy Parameter
Theoretically, the SI of a non-influential parameter is zero. In the case of the PAWN method, SI of zero means that the conditional and unconditional CDF coincides, i.e., fixing the parameter into consideration will have no influence on the model output distribution. The value of zero in the Morris screening test indicates that the parameter has no influence on the model output. However, numerical algorithms and computations are utilized to implement SI, and small but non-zero values can be obtained for non-influential parameters. To set a threshold to identify non-influential parameters, a dummy parameter approach was proposed by [29]. The SI of this dummy parameter provides an indication of the approximation error of the sensitivity analysis. The parameters having SI greater than the SI of the dummy parameter are considered as the influential parameters. It should be noted that the dummy parameter is not added to the model; its SI is calculated by using the sampled data.

Crop Model: Aquacrop Open Source
Crop growth simulation model Aquacrop was originally developed by the Food and Agriculture Organization (FAO) of the United Nations. As the code of the original model is not accessible, it does not allow, for example, the implementation of data assimilation techniques such as the ensemble Kalman filter, which are increasingly employed for ingesting remote sensing observations into crop models [30]. In this work we have used the open-source Matlab version of the model called Aquacrop Open Source (AOS) developed by [31].
Aquacrop model evolved from the most commonly used empirical approach to assess crop yield response to water [47,48], shown as Equation (3) where Y x and Y a are the maximum and actual yield, ET x and ET a are the maximum and actual evapotranspiration, and k y is the proportionality factor between relative yield loss and relative reduction in evapotranspiration. The separation of ET into Tr and E avoids the confounding effects of the nonproductive consumption use of water (E), which is important especially during incomplete ground cover, and leads to the conceptual equation at the core of the model [48] (Equation (4), derived by [49]). Yield (Y) as a function of biomass (B) and harvest index (HI) is described as Equation (5): where Tr is the crop transpiration (in mm), and WP is the water productivity (kg of biomass per m 2 and per mm of cumulated water transpired over the time period in which the biomass is produced), which is a constant for the given climatic condition. By normalizing it for different climatic conditions, WP becomes a conservative parameter [50]. Reference evapotranspiration (ET 0 ) is provided as a model input: it can easily be estimated using an external calculator [51] that uses the Penman-Monteith equation, following the procedures outlined in FAO Irrigation and Drainage paper no. 56 [52]. The model progressed by separating by evapotranspiration (ET) into crop transpiration (Tr) and soil evaporation (E) using canopy cover (CC). Tr is proportional to the fraction of land surface covered by the canopy, while E is proportional to the surface not covered by the canopy, with some minor adjustments. Tr in the absence of water stress is calculated using the adjusted green canopy cover (CC*) and the crop coefficient of the full-grown crop, prior to senescence (Kcb), as in Equation (6): A distinctive feature of Aquacrop is the use of green canopy cover instead of Leaf Area Index (LAI) to simulate canopy expansion. This allows Aquacrop to be easily used with remote sensing data [48]. Soil evaporation is based on the two-stage Ritchie's approach [53], adopted also in other crop models (e.g., DSSAT crop model). While the yield-transpiration relationship accounts for the biomass accumulation, canopy expansion and crop phenology are simulated following the thermal time concept, which uses growing degree days (GDDs) as an internal clock of crop development as described in [54], with minor modifications. CC development is simulated through a two-stages curve. An exponential equation is used up to half the maximum canopy cover (Equation (7)): where CC is the fractional canopy cover of soil at time t, CC 0 is the initial CC (at time 0) and CGC is the canopy growth coefficient in the percentage of existing CC at time t. The second half of canopy development follows an exponential decay equation (Equation (8)): where CC x is the maximum canopy cover for optimal conditions. An empirical canopy decline coefficient (CDC) is used to characterize the CC fractional reduction overtime after maturity has been reached (Equation (9)): The soil-water balance is simulated first by using the runoff curve number (CN) approach of the United States Department of Agriculture (USDA) Soil Conservation Service to separate runoff and infiltration, and then an exponential drainage function is used to simulate downward flow when the soil water content is above field capacity (Equation (10)): where θ i is the actual soil water content, θ SAT is the soil water content at saturation θ FC is the soil water content at field capacity, ∆t is the time step (daily) and τ is a coefficient to account for drainage characteristics, derived from the saturated hydraulic conductivity (Ksat). The effect of water-limiting conditions on crop growth is calculated with the ratio of actual available water to total available water TAW (θ SAT -θ FC ) with an approach evolved from the FAO Irrigation and Drainage paper no. 56 [52]. Water stress coefficients (Ks) act as modifiers (with values ranging from 0 to 1) of target processes and are used to modulate transpiration, canopy expansion, canopy senescence and pollination when the readily available water (RAW) is depleted. RAW is defined as a fraction of TAW by the p_up parameters. Each of the Ks coefficients is defined on the basis of the fraction of actual to total available water, using three parameters: the upper threshold (p_up, when water stress begins: Ks = 1), the lower threshold (e, when water stress is complete: Ks = 0) and the shape factor of the curve (fshape_w).
Extensive descriptions of Aquacrop processes are presented in [30,55]. AOS does not differ in the calculation procedures but is provided as an open-source Matlab version of Aquacrop [31,56]. Test simulations were conducted to verify the suitability of AOS 5.0a to reproduce accurately the calculations and outputs from the original FAO Aquacrop model, considering a range of crop types and environments: wheat production in Tunisia, rice and wheat production in Hyderabad, India, and potato production in Brussels, Belgium. The maximum RMSE and minimum r 2 for crop yields across all simulations were 0.045 t ha −1 and 0.993, respectively, which confirms the ability of AOS to reproduce successfully Aquacrop calculations and outputs [56].
Recently, a version of AOS, AOS v6.0a, has been released. AOS v6.0a implements all calculation procedures performed by the original FAO Aquacrop v6.0a model, with the exception of soil fertility stress and soil salinity stress, as these are yet to be tested extensively for different crop types and environments worldwide [36].
There is a total of 89 parameters in the AOS model, as mentioned in the reference manual [57]. Both Aquacrop v6.0a and AOS software include crop files with default parameter values of biomass and canopy growth, thermal time-based phenological stages, transpiration coefficients and temperature and water stresses. In AOS each of these parameters is assigned a specific abbreviation and description in a clearer way with respect to Aquacrop v6.0a. The abbreviations and descriptions of the 54 AOS parameters considered in this work for the SA are listed in Table 1, along with their Aquacrop v6.0a counterparts, as presented in two of the most recent and extensive global SA [58,59]. For clarity, when referring to the previous study, all the parameter abbreviations used in this paper will be written using the AOS notation.

Setting up the Parameters for AOS Model
In this study, the EFAST and PAWN methods of SA have been applied to the AOS model and compared in terms of parameter ranking. We used all the parameters that were not fixed as mentioned in the AOS manual [57] to be analyzed for SA. The winter wheat grain yield was considered as the output. AOS parameters are reported in Table 1. The range of variation of the selected 54 parameters (reported in Table 1) was determined based on the AOS manual [57] and default values provided with the open-source AOS model database for the wheat crop and also on the previous studies of sensitivity analysis, carried out with the original FAO Aquacrop model [7,37,38]. Since there was no prior information on parameter frequency distributions, they were assumed to have an a priori uniform distribution. It should be noted that [37], in their study, confirms the parameter ranges are more important than their distributions in SA studies.

Sites Description
In the scenarios applied for the SA test, we used the climate and soil data from two experimental sites. The first site was located in Maccarese (41.833 • N, 12.217 • E; 8 m a.s.l), western coast of Central Italy, near the Fiumicino Airport, Rome. This area has a typical Mediterranean climate with wet autumn and dry hot summer falling into the Koppen-Geiger class Csa (temperate, dry summer, hot summer) [60]. The average yearly maximum temperature is 27.4 • C and the average minimum temperature is 5 • C. The average annual precipitation is 812 mm [61]. The soil in Maccarese is Cutanic Luvisol, with a prevailing sandy clay loam texture, becoming more clayey toward the north-east of the site.
The other site was located in Xiaotangshan (40.167 • N, 116.4 • E; 36 m a.s.l), Chanping district, Beijing, China. It is characterized by a continental climate, with a cold, dry winter and a hot, wet summer, belonging to the Koppen-Geiger [60] class Dwa (cold, dry winter, hot summer), i.e., monsoon-influenced, hot-summer humid continental climate. The soil in the experimental area was loam. The mean annual precipitation is approximately 602 mm. The mean temperature is approximately −4.5 • C in winter and 26.2 • C in summer. There are on average 180 frost-free days per year. The main crop is winter wheat. Table 2 shows sowing and harvest dates used for model simulations at the different test sites [38], which were set on the basis of local agronomic practice. Soil information for both the test sites is presented in Table 3.

Weather Data
Sensitivity analysis was carried out for four different wheat growth seasons (2014-2018) for the Maccarese site and for three seasons (2015-2018) for the Xiaotangshan site. The meteorological data used in this study to drive the model simulations for the SA were obtained from the Regione Lazio Agrometeorological Service, Italy (http://www.arsial.it/portalearsial/agrometeo/index.asp) and the National Meteorological Information Center of the China Meteorological Administration (http://cdc. cma.gov.cn/). The daily minimum and maximum temperatures, relative humidities, wind speed, rainfall data and daily solar radiation for the period were provided. The daily reference evapotranspiration was calculated based on the FAO Penman-Monteith method as described in [51] and the ET 0 calculator [62]. The minimum and maximum temperature, rainfall and reference evapotranspiration for both test sites are shown in Figure 1. For both sites, the years considered in this study allowed to include wet and dry conditions, which were assessed by examining the Standardised Precipitation-Evapotranspiration Index SPEI (https://spei.csic.es/index.html). In particular, the spring of 2017 was severely dry in Maccarese and also in Xiaotangshan, although for the latter the drought was limited to the month of May as shown in the Supplementary Materials ( Figure S1).  Figure 2 shows the sensitivity of the AOS model parameters obtained by applying the Morris method for the two different locations, climates and years. In Figure 2, the average value for all the years is plotted, and the error bars indicate the variation (standard deviation) of the average value.  Figure 2 shows the sensitivity of the AOS model parameters obtained by applying the Morris method for the two different locations, climates and years. In Figure 2, the average value for all the years is plotted, and the error bars indicate the variation (standard deviation) of the average value.

Morris Method
A large magnitude of µ* indicates the importance of the model parameter. A high standard deviation of σ implies either a non-linear relationship or interaction with other parameters. The parameters CGC and HIstart were found to be of high importance for both the study sites regardless of the climate. The number of relevant parameters differs slightly for each year and study site.  Table 1.
A large magnitude of μ* indicates the importance of the model parameter. A high standard deviation of σ implies either a non-linear relationship or interaction with other parameters. The parameters CGC and HIstart were found to be of high importance for both the study sites regardless of the climate. The number of relevant parameters differs slightly for each year and study site. When considering all years of data and study sites, 29 out of 55 (including dummy) were found important for yield output. These 29 parameters are plotted in Figure 2. For clarity purposes, only the parameters with the mean value of or greater than 1 t ha −1 are labeled. Figure 3 shows the EFAST indices plots for all the parameters and all the years for the Maccarese test site. Both the first-order (Si) and the total-order (Sti) indices are calculated, as proposed in [19], to distinguish between the main effect and the interaction effect.  Table 1. Figure 3 shows the EFAST indices plots for all the parameters and all the years for the Maccarese test site. Both the first-order (S i ) and the total-order (St i ) indices are calculated, as proposed in [19], to distinguish between the main effect and the interaction effect.

Variance-Based Extended Fourier Amplitude Sensitivity Test (EFAST)
For both study sites and all the years of weather data, the most relevant parameter was CGC, with the exception of 2016-2017 for both test sites; thus, this parameter is responsible for most of the variance of the output (Figures 3 and 4). The 2016-2017 season was particularly dry in Maccarese, with a three-month accumulated standardized precipitation-evapotranspiration index (SPEI) below −1 from February to July, indicating drought conditions ( Figure S1). During that season, the parameter that resulted to have the highest sensitivity was HIstart. For Xiaotangshan test site ( Figure S2 Table 1. A large magnitude of μ* indicates the importance of the model parameter. A high standard deviation of σ implies either a non-linear relationship or interaction with other parameters. The parameters CGC and HIstart were found to be of high importance for both the study sites regardless of the climate. The number of relevant parameters differs slightly for each year and study site. When considering all years of data and study sites, 29 out of 55 (including dummy) were found important for yield output. These 29 parameters are plotted in Figure 2. For clarity purposes, only the parameters with the mean value of or greater than 1 t ha −1 are labeled. Figure 3 shows the EFAST indices plots for all the parameters and all the years for the Maccarese test site. Both the first-order (Si) and the total-order (Sti) indices are calculated, as proposed in [19], to distinguish between the main effect and the interaction effect.  . First-order and total-order SI for different years for Maccarese. Only the parameters with a total sensitivity index higher than the sensitivity index of the dummy parameter are shown. The details on the parameters can be found in Table 1.

Variance-Based Extended Fourier Amplitude Sensitivity Test (EFAST)
For both study sites and all the years of weather data, the most relevant parameter was CGC, with the exception of 2016-2017 for both test sites; thus, this parameter is responsible for most of the variance of the output ( Figure 3; Figure 4). The 2016-2017 season was particularly dry in Maccarese, with a three-month accumulated standardized precipitation-evapotranspiration index (SPEI) below −1 from February to July, indicating drought conditions ( Figure S1). During that season, the parameter that resulted to have the highest sensitivity was HIstart. For Xiaotangshan test site ( Figure  S2  First-order and total-order SI for different years for Maccarese. Only the parameters with a total sensitivity index higher than the sensitivity index of the dummy parameter are shown. The details on the parameters can be found in Table 1. The average indices for all the years are plotted in Figure 4, error bars indicating the extreme values for the total-order sensitivity indices. As expected, a small, but non-zero, value was estimated for the dummy parameter SI. For easy visualization and interpretation, the sensitivity indices are sorted in descending order. The parameters with SI values less than the SI of the dummy parameter were, thus, considered as non-influential parameters. The average SI for the Xiaotangshan location is shown in Figure 4.  Figure 5 and Figure 6 show an example of two different obtained parameters with different and similar CDFs, illustrating the behavior of non-influential and influential parameters using the PAWN method. In Figure 5 it can be seen that for the parameter CGC the unconditional CDF and conditional CDF are different, and the value of the KS statistics is above the critical value at 95% confidence interval. Therefore, the model output has a high sensitivity to the parameter, which can be thus considered highly influential. On the contrary, Figure 6 shows an example of another parameter, CDC, with similar unconditional and conditional CDFs; the value of the KS statistic is below the critical value, and thus the parameter is non-influential and has a low SI. The results of EFAST for winter wheat were quite similar for both study sites. Although, there were minor differences in the specific value of the St i , and the general importance of the analyzed parameters was identical for both study sites. In the end, we used the average of the total sensitivity indices across all the years and study sites to find a similar set of sensitive parameters that could be adapted when calibrating AOS for regional studies (Figure 4).

Density-Based PAWN Method
The average indices for all the years are plotted in Figure 4, error bars indicating the extreme values for the total-order sensitivity indices. As expected, a small, but non-zero, value was estimated for the dummy parameter SI. For easy visualization and interpretation, the sensitivity indices are sorted in descending order. The parameters with SI values less than the SI of the dummy parameter were, thus, considered as non-influential parameters. The average SI for the Xiaotangshan location is shown in Figure 4.

Figures 5 and 6 show an example of two different obtained parameters with different and similar
CDFs, illustrating the behavior of non-influential and influential parameters using the PAWN method. In Figure 5 it can be seen that for the parameter CGC the unconditional CDF and conditional CDF are different, and the value of the KS statistics is above the critical value at 95% confidence interval. Therefore, the model output has a high sensitivity to the parameter, which can be thus considered highly influential. On the contrary, Figure 6 shows an example of another parameter, CDC, with similar unconditional and conditional CDFs; the value of the KS statistic is below the critical value, and thus the parameter is non-influential and has a low SI.  The PAWN index (Ti) for the different years is presented in Figure 7 for Maccarese and in Figure  S3 (Supplementary Material) for Xiaotangshan. The average values of the indices for all the years are plotted in Figure 8, in which error bars represent the range of values for the PAWN sensitivity indices across the years. In comparison to EFAST, the sensitivity index appeared to be more stable. For both study sites, considering all the years of weather data, the most relevant parameter was CGC, also identified by EFAST method, and is also responsible for most of the variance of the model output. The second most influential parameter was HIstart. However, differently from EFAST, this parameter did not have a sensitivity index higher than CGC even in the water-stressed conditions of 2017. Again, as expected, a small but non-zero value was estimated for the SI of the dummy parameter. The parameters with the SI less than the SI of the dummy parameter are thus non-influential parameters. The average SI for the Xiaotangshan location is shown in Figure 8.  The PAWN index (Ti) for the different years is presented in Figure 7 for Maccarese and in Figure  S3 (Supplementary Material) for Xiaotangshan. The average values of the indices for all the years are plotted in Figure 8, in which error bars represent the range of values for the PAWN sensitivity indices across the years. In comparison to EFAST, the sensitivity index appeared to be more stable. For both study sites, considering all the years of weather data, the most relevant parameter was CGC, also identified by EFAST method, and is also responsible for most of the variance of the model output. The second most influential parameter was HIstart. However, differently from EFAST, this parameter did not have a sensitivity index higher than CGC even in the water-stressed conditions of 2017. Again, as expected, a small but non-zero value was estimated for the SI of the dummy parameter. The parameters with the SI less than the SI of the dummy parameter are thus non-influential parameters. The average SI for the Xiaotangshan location is shown in Figure 8. The PAWN index (T i ) for the different years is presented in Figure 7 for Maccarese and in Figure S3 (Supplementary Material) for Xiaotangshan. The average values of the indices for all the years are plotted in Figure 8, in which error bars represent the range of values for the PAWN sensitivity indices across the years. In comparison to EFAST, the sensitivity index appeared to be more stable. For both study sites, considering all the years of weather data, the most relevant parameter was CGC, also identified by EFAST method, and is also responsible for most of the variance of the model output. The second most influential parameter was HIstart. However, differently from EFAST, this parameter did not have a sensitivity index higher than CGC even in the water-stressed conditions of 2017. Again, as expected, a small but non-zero value was estimated for the SI of the dummy parameter. The parameters with the SI less than the SI of the dummy parameter are thus non-influential parameters. The average SI for the Xiaotangshan location is shown in Figure 8.  Only the parameters with a total sensitivity index higher than the sensitivity index of the dummy parameter are shown. The details on the parameters can be found in Table 1. Identifying sensitive parameters using the PAWN method is suitable in the cases where the output distribution is not normal, and it follows skewed or multi-modal distribution [25]. An example of the output distribution is compared with 54 different statistical distributions. The top six best-fitted PDFs are mentioned in Table 4. The comparison is made on the basis of KS statistics and Anderson-Darling Statistics, with low values representing the best distribution. As can be seen, the Generalized Extreme Distribution has the lowest statistic value, with an Anderson-Darling statistic of 0.16401 and is second lowest in KS statistics. It is the best-suited PDF for the data, the reason for the choice of the PAWN method and is feasible to implement for identifying the most influential set of parameters for the local calibration. The Generalized Extreme Distribution is plotted in Figure 9, and statistics including the p-value and critical values at alpha (α) are presented in Table 5. The sample was from conditional CDF of the PAWN method, consisting of 200 samples, Nc = 200 (Section 2.1.3). It shows the statistics value and rank, indicating the raking of the distribution when compared with 54 PDFs (Table 4). Table 4. Best probability distribution with Kolmogorov-Smirnov (KS) and Anderson-Darling statistics, ranking the best from top to bottom. Superscript in red indicates the rank of the KS and Anderson-Darling statistics.

KS Statistics
Anderson-Darling Statistics Gen. Extreme Values 0.03253 (2) 0.16401 (1) Johnson SB 0.03047 (1) 0.19169 (2) Gumbel Max 0.05197 (3) 0.4357 (3) Logistic 0.12206 (4) 3.7302 (4) Hypersecant 0.1276 (5) 4.0165 (5) Normal 0.12242 (6) 4.0999 (6)  Identifying sensitive parameters using the PAWN method is suitable in the cases where the output distribution is not normal, and it follows skewed or multi-modal distribution [25]. An example of the output distribution is compared with 54 different statistical distributions. The top six best-fitted PDFs are mentioned in Table 4. The comparison is made on the basis of KS statistics and Anderson-Darling Statistics, with low values representing the best distribution. As can be seen, the Generalized Extreme Distribution has the lowest statistic value, with an Anderson-Darling statistic of 0.16401 and is second lowest in KS statistics. It is the best-suited PDF for the data, the reason for the choice of the PAWN method and is feasible to implement for identifying the most influential set of parameters for the local calibration. The Generalized Extreme Distribution is plotted in Figure 9, and statistics including the p-value and critical values at alpha (α) are presented in Table 5. The sample was from conditional CDF of the PAWN method, consisting of 200 samples, N c = 200 (Section 2.1.3). It shows the statistics value and rank, indicating the raking of the distribution when compared with 54 PDFs (Table 4). Table 4. Best probability distribution with Kolmogorov-Smirnov (KS) and Anderson-Darling statistics, ranking the best from top to bottom. Superscript in red indicates the rank of the KS and Anderson-Darling statistics.

Elementary Effects Using the Morris Method
The results obtained in this study, performed using a contrasting range of climatic scenarios of winter wheat growing areas in Maccarese, Italy (Mediterranean climate) and Xiaotangshan, China (continental climate), allowed to obtain essential information on the sensitivity of the AOS model, especially required for their application within regional-scale studies [63][64][65]. It is well known that the results of the sensitivity analysis studies depend on the boundary conditions chosen [66]. In this work, these conditions are climate datasets, actual data from both sites and the range of variation of the parameters [7,37,38]. The Morris method was used for initial screening of the parameters and to identify the influential and non-influential parameters. The influential parameters were found to be 29 for the Maccarese test site and 28 for the Xiaotangshan test site out of the total 54 model parameters considered in this study. These parameters are mentioned in Table 6, and the description of the parameters is reported in Table 1. A study [37] on the SA of the Aquacrop model used the mean value of 0.25 as a threshold, and parameters below 0.25 t ha −1 were classified as non-influential parameters. It was found that the parameters (translated in AOS notation) exc, dHI0, p_up1, Tmin_up, Tmax_up, p_up4, SxTopQ and SxBotQ were non-influential (<0.25 t ha −1 ), and the other parameters PlantPop, Flowering, HI0, dHI_pre, a_HI, Kcb, p_up2 and fshape_w2 had low but non-negligible impacts (0.25-1.0 t ha −1 ) for all the climatic scenarios implemented. In this study, the CGC parameter was found to be highly influential in all the conditions for both test sites. Another study carried out by [7] found almost the same parameter set of the non-influential parameters as by [37] in all the climatic conditions, with the exception of the soil curve number (CN) and the shape factor for water stress coefficient inducing early senescence (fshape_w3). It should be noted that the threshold considered in the study was 0.1 t ha −1 i.e., less than that used by [7]. The tabular comparison of the sensitive parameters is presented in Table 6.

Variance-Based Extended Fourier Amplitude Sensitivity Test (EFAST)
The Morris method is only used for initial screening, and it does not quantify the higher-order interaction between the parameters. This can be assessed by using the variance-based EFAST method and the density-based PAWN method for GSA. The comparison between these algorithms is not meaningful if compared on the basis of sensitivity indices as they have different rationales and backgrounds. This work applies a comparison based on the ranking of the parameters in both methods. So, to achieve this goal, all the input parameters, i.e., 55 including a dummy parameter, are used to carry out SA in both the methods. Figure 3 shows both the first-order and the total-order sensitivity indices for the EFAST methods for all the years considered in this study. The sensitive parameters are, in general, common in all the climate years considered, but with different values of the SI. These common parameters are Emergence, Senescence, Maturity, HIstart, CGC, Kcb, WP, WPy and HI0, Tmin_up, GDD_up and fshape_b. The average sensitivity indices plots for all the years identify the non-influential parameters having a sensitivity less than the sensitivity of the dummy parameter. As can be seen from Figure 4, dHI0, p_up3, p_lo2, Tmax_lo, p_up1, PlantPop, Tmin_lo, Zmin, fshape_w2, fshape_w3, fshape_w4 and SxBotQ had low sensitivities in comparison to the dummy parameter. These parameters are clearly identified as non-influential parameters.
Similarly, for the Xiaotangshan test site, the most influential parameters are reported in Table 6. The highly influential parameters for both the locations were CGC, HIstart, WP, HI0, Kcb, Emergence, Senescence, WPy, GDD_up, Maturity, fshape_b, Tmin_up, GDD_lo and CDC (Figure 4). These parameters are in common for both the test sites, although the ranking is not the same. Authors in [37] identified that for winter wheat in a temperate maritime climate, many parameters had an equivalent effect on the output variance, among which parameters describing canopy development were SeedSize, CCx, CGC, Senescence, WP, GDD_up and Zmax. In a previous study [7], it was found that the parameters describing the phenological cycle of the crop and the water stress were among the most highly influential. The ranking of the parameters was different from that of [37], due to the climatic conditions employed. Also, the total-order effects were different than those of the first-order effects, indicating more complex higher-order interacting effects rather than individual effects. It should be noted that both these studies applied EFAST on the results of the Morris method, that is not the case of our study. The EFAST method for SA was directly applied to all the parameters in [38] under different irrigation treatments. In the presence of rainfall, CGC, CCx, Emergence, PlantPop, SeedSize, Kcb, Maturity, Zmax, p_up3, MaxRooting, GDD_up and WP were the most sensitive to yield. In normal irrigation, CGC, CCx, Emergence, PlantPop, SeedSize, Kcb, Maturity, Zmax and p_up3 were among the highly influential, and in over-irrigation, CGC, CCx, Emergence, PlantPop, SeedSize, Kcb, Maturity, Zmax, p_up3 and CDC were the most sensitive parameters. The translations of these parameters of Aquacrop and AOS model are mentioned in Table 1.

Density-Based PAWN Method
In the PAWN method, the sensitivity of the model output (yield) to the parameter x i is measured by the distance between the conditional CDF of output obtained by fixing x i , to a specific value and the unconditional CDF. Since all the variability related to x i , both due to direct and interaction effects, is removed in the conditional CDF, the comparison illustrates the total effect of x i . As a distance measure, the maximum absolute difference between two CDFs, i.e., the Kolmogorov-Smirnov statistic, is considered. The PAWN SI for the different years for Maccarese location is presented in Figure 7. The sensitive parameters are similar for all different years, but with a varying degree of the PAWN index. The average index for all the years is plotted in Figure 8. Table 6 summarizes the set of sensitive parameters based on the PAWN index.
The PAWN method is simpler and less computationally demanding in comparison to the EFAST method. Out of all the methods tested in this study, the Morris method was the simplest and required the lowest computational time, but the disadvantage with the method is that it does not quantify higher-order effects. The PAWN method, despite being simple, has two problems: first is the reliability of the CDF produced, which may become unstable as can be observed for both study sites. Fifteen parameters were identified as most influential for Maccarese study site, whereas 46 parameters appeared to have sensitivity higher than the sensitivity of the dummy parameter for the Xiaotangshan study site. However, the parameters with the highest sensitivity were in common for both the study sites. The other drawback of PAWN is the use of KS statistics to compute SI that leads to slow convergence.

Sensitive Parameters in Relation to Model Processes and Agronomic Conditions
AOS requires calibration of parameters to properly simulate crop growth, especially in new locations: SA can be used to identify the parameters that have the highest impact on crop yield, before proceeding to calibration. The interpretation of the SA results, in light of previously published literature and of the knowledge of model calculations and simulated processes, can provide useful insights for advanced as well as basic use of the model by researchers. Aquacrop parameters are divided into two main groups: conservative and non-conservative parameters. Conservative crop parameters do not change with time, cultivar, location nor water regime and in principle should require no adjustment to local conditions. If simulations do not match observed values, the first thing to look at is the fertility regime before changing the parameters. Non-conservative parameters, on the other hand, are dependent on the cultivar (e.g., phenology parameters) and conditions.
In our case, both in Maccarese and in Xiaotangshan, both methods (EFAST and PAWN) found CGC, HIstart, WP, HI0 and Kcb to be the first five most sensitive parameters (Table 6), even though the order changed slightly between locations and methods. CGC, WP, HI0 and Kcb are conservative parameters that showed a high impact on crop yield, while HIstart is a phenological non-conservative parameter. These parameters are mainly related to crop transpiration, which is directly used by AOS to calculate biomass accumulation (AOS is a water-driven model). CGC describes how fast the canopy cover (CC) increases (Equation (7)). CC and Kcb are directly involved in the estimation of crop transpiration, which in turns determines biomass production through water productivity (WP). Grain yield is also highly influenced by the harvest index (HI0) and by the thermal time to start yield formation (HIstart). High fluctuations of yield estimations can therefore be expected when adjusting these parameters, and model users should be careful in this regard. Moreover, these parameters are also used to mimic the effects of fertility stress on crop growth in a semiquantitative way, since nutrient cycling is not directly incorporated into the model [59]. Testing parameter sensitivity under stress conditions is important for calibration purposes [67].
Other phenological parameters (e.g., Emergence, Senescence and Maturity) followed up in the SA with both locations and methods. Authors in [68] also included Emergence and Maturity as sensitive parameters for winter wheat cultivation in China, using the EFAST method. Authors in [7] found among the most influential parameters Senescence, HIstart and Maturity in winter wheat cultivated in a location in Central Italy (in the same region of Maccarese) and two locations in China (with Xiaotangshan being one of them). Similarly, authors in found Emergency, Senescence and Maturity to be among the highest influential parameters for winter wheat in Belgium. Authors in [38] also reported HIstart and Maturity among influential parameters for wheat cultivated near Beijing (using the EFAST method).
Water stress parameters (p_up, p_lo and fshape_p) are often the target of SA and calibration of Aquacrop. However, datasets for water-deficient conditions are usually needed to properly estimate the impact of these parameters on crop yield [69]. In this study, winter wheat was rainfed and not irrigated. Indeed, in the two locations considered, rainfall is frequent and evapotranspiration is low during the growing period (Figure 1), so it is not surprising that water stress parameters had lower sensitivity with respect to the canopy and phenological parameters for both locations and methods (Table 4). There are several examples that support this argument. Authors in [37] found water stress parameters to have a higher sensitivity order for winter wheat under drier conditions and more sandy soils. Authors in [33] simulated satisfactorily barley growth in southern Italy under rainfed conditions without the need to adjust water stress parameters for each year, while other crop growth and phenological parameters were specifically changed over the years. Authors in [70] obtained good simulations of winter wheat yield in Iran under moderate water stress without including water stress parameters in their SA and calibration. The SA carried out by [38] highlighted that the parameters with higher-order sensitivity to winter wheat yield in Beijing were mainly related to phenology, biomass accumulation and canopy growth, with p_up3 (the threshold for water stress-related senescence) being the only water stress-related parameter included. In another recent global SA [59], maize water stress parameters were less sensitive than crop growth and phenology parameters, under all the irrigation treatments (designed to obtain a range of water stress intensities, from no to moderate stress). These experiences confirmed our results and highlighted that water stress parameters do not have high impact on winter wheat yield when water-limiting conditions are moderate.

Conclusions
In this paper, we compared the application of two global SA techniques, the variance-based EFAST and the moment-independent density-based PAWN method, to the analysis of 54 parameters of the AOS model, a crop growth simulation model widely used for canopy cover, biomass and yield simulations. With respect to the widely used variance-based sensitivity indices, the main advantage of PAWN is that it can be applied whatever the output distribution may be. The comparison was performed in terms of parameter ranking for different climates and geographical locations. In all the methods, highly influential parameters were common, and the ranking of these parameters was similar in EFAST and PAWN.
The PAWN method could be considered more adequate in our study, since the frequency distribution of the output variable was right-skewed. The use of a dummy parameter approach as a viable option to set a threshold value for parameter screening was demonstrated for all the methods considered. The results reported in this study provide insight into identifying the most influential parameters of the AOS model for two different geographic locations with two different climates.
This set of most influential parameters can be used for model simplification and calibration of the most sensitive parameters. Simplification of the model is necessary for regional-scale studies and particularly for remote sensing data assimilation. Both EFAST and PAWN methods identified the most influential parameters for all different years for both the test sites, but with a varied sensitivity and ranking. According to the PAWN method, the highly influential parameters belonged to crop parameters specifically, to the canopy and phenological development (CGC, HIstart, Senescence, maturity), biomass/yield production (WP, WPy) and temperature stress (Tmin_up, Tmin_lo, GDD_up) rather than root development and water stress.
This study thus concludes that variance-based and density-based global SA methods are both useful, with the condition that the variance-based method should be preferable where a variance is an adequate proxy of the model output, i.e., where the model output distribution is normal. Although this work is limited to compare these two (EFAST and PAWN) methods on parameter ranking, for two different climates with seven years data overall, to verify the robustness of the comparison, it can be extended on the basis of efficiency, convergence or computational cost with other geographic locations, climates and crops.