Sample size, range of parameters and time-dependent effects on global sensitivity analysis in sugarcane modelling

Abstract

Process-based crop models (PBCM) for sugarcane present many genotype parameters compared to other crops, making it harder to calibrate. Global sensitivity analysis (GSA) has thus become an important tool for understanding, calibrating and further developing PBCMs. This work used a recently updated Agronomic Modular Simulator for Sugarcane (SAMUCA) to simulate crop growth and development, conducted with two treatments: with green cane trash blanket (GCTB) and under bare soil (BARE). Using the extended Fourier Amplitude Sensitivity (eFAST) algorithm, GSA was performed on the 24 genotype parameters of SAMUCA. The objective was to determine the sample size (SZ): how many samplings are necessary to quantify the sensitivity indices. Additionally, we aimed to assess the influence of parameter range variation and identify which genotype parameters explain the highest variance in simulations of the SAMUCA model under BARE and GCTB conditions. The results showed that SZ greatly affected the convergence and sensitivity indices, and the SZ required here needed to be >2049 for the analysis to cover all variables. Two sets of parameter ranges were used for analysis (the first used maximum and integer values of each parameter reported in the literature; the second applied 25% perturbation to the previously calibrated values). The results indicated that the parameter range affected the parameters' order of importance. Furthermore, we identified that at different phenological stages of sugarcane development, distinct parameters were responsible for explaining the most variance of the output. However, there was no difference among ratoons or interference in the results of BARE or GCTB.

Introduction

In agriculture, process-based crop models (PBCM) represent the state of the art for simulating crop growth and development (Jones et al., 2017; Marin et al., 2017). When properly calibrated, they are commonly used to simulate crop growth and development under different conditions, thus being able to test hypothetical management, climate and soil scenarios (Faivre et al., 2009). Scientists and decision-makers have used modelling as a tool to address issues related to the sugar and bioenergy sectors, including climate change (Marin et al., 2013; Singels et al., 2013; Jones et al., 2015), plant breeding (Hoffman et al., 2018), risk analysis (Everingham et al., 2002) and crop forecast (Everingham et al., 2016).

Sugarcane is a crucial crop for world bioenergy (Raza et al., 2019), and several authors have studied sugarcane modelling (Keating et al., 1999; Singels and Bezuidenhout, 2002; Inman-Bamber and Smith, 2005; Thorburn et al., 2005; Marin and Jones, 2014; Valade et al., 2014; Jones and Singels, 2018; Vianna et al., 2020). Singels (2013) presented a detailed review of the main sugarcane models in the literature, highlighting sugarcane as one of the crops with highest need to be represented in PBCM given its specific farming systems and logistic requirements. To represent its physiological complexity, sugarcane PBCM have many genotype parameters compared to other crops, such as maize and wheat. For instance, the sugarcane models DSSAT/CASUPRO (Villegas et al., 2005) and DSSAT/CANEGRO (Jones and Singels, 2008) have 42 and 22 cultivar coefficients, respectively, while DSSAT/CERES-MAIZE (Jones et al., 1986) and DSSAT/CERES-WHEAT (Ritchie and Otter, 1985) have only respectively six and seven parameters to be calibrated.

According to Sinclair and Seligman (1996), the development of different PBCM by more research groups allows to improve the understanding of processes. In this sense, Marin and Jones (2014) developed the Agronomic Modular Simulator for Sugarcane (SAMUCA) focusing on the specific characteristics of sugarcane farming systems in Brazil. Recently, the SAMUCA model was improved by reducing the uncertainties around the soil water balance, heat flux and physiological mechanisms such as carbon partition, photosynthesis, tillering and root growth (Vianna et al., 2020).

As any other PBCM, SAMUCA represents a simplification of the real system and requires several parameters whose determination is a problem for practical operational applications (Makowski et al., 2002). Most parameters are acquired through field observations, which are expensive and time-consuming, and the acquisition of certain parameters is difficult.

Furthermore, for reliable simulations, accurate parameter estimation is required (Guérif and Duke, 2000; Wallach et al., 2019a), and so several parameter estimation algorithms have been developed, part of which are based on Bayesian approaches (He et al., 2010, 2009; Sreelash et al., 2012; Marin et al., 2017; Sheng et al., 2019; Zhang et al., 2020). To some extent, these methods solved the problem of difficult-to-acquire parameters, and they are quite efficient but applicable to a small number of parameters only (Varella et al., 2010). Still, the inclusion of many parameters in a PBCM raises a dilemma related to the difficulty in simultaneously estimating all unknown parameters and ensuring at the same time that their biophysical meaning is coherent.

In practice, it is well known that only part of the parameters is usually responsible for most of the model uncertainty, while most of them have only minor influence (Varella et al., 2010; Li et al., 2019; Zhang et al., 2020). The parameter sensitivity analysis (SA) method can identify the most important parameters for a given model output variable, which allows users to focus on the most important model parameters during the calibration process. Furthermore, based on the SA, the balance and robustness of the model can be analysed for future improvement, model development and applications (Fraedrich and Goldberg, 2000; Confalonieri, 2010; Chu-Agor et al., 2011; Hirabayashi et al., 2011).

The SA can be divided into two groups: the local sensitivity analysis (LSA) and the global sensitivity analysis (GSA) (Wallach et al., 2019b). The LSA consists of changing a single parameter at a time, while the other parameters are kept at their reference values; in other words, this method is based on the local derivatives of the model's output concerning the variation of a single parameter, which indicates how strong the output changes around the reference parameter values (Saltelli et al., 1999). The GSA allows you to evaluate the entire uncertainty range of parameters, considering changes in all parameters along with their range, as well as the interactions among parameters (Saltelli et al., 1999).

The GSA methods can also be classified into three groups: screening, regression and variance; all of these follow the Bayesian sampling principle (Wei, 2013). The most used screening method is the Morris method, which permits to define the most important model parameters and it is often considered a qualitative method (Morris, 1991; Dejonge et al., 2012). Regression methods, such as the Partial Rank Correlation Coefficient (PRCC), provide the correlation between the model's output and the selected parameters and they are mandatory when the parameter and the model outputs have a monotonic relationship (Marino et al., 2008; Krishnan and Aggarwal, 2018). Methods based on variance are the most used one, and the three mains are: Sobol (2001), Fourier Amplitude Sensitivity Test (FAST) (Cukier et al., 1978) and the extended Fourier Amplitude Sensitivity Test (eFAST) (Saltelli et al., 1999); they provide what parameters cause the greater variability for the model's output, and they usually demand a high computational cost. To define which method is the most suitable, some properties of model must be known (linearity, prior distribution of parameters and monotonic), furthermore considering the number of parameters to be evaluated and the computational cost (Iooss and Lemaître, 2015). In this sense, Sobol and eFAST are the methods most applicable to any type of PBCM, but while Sobol is very computationally expensive, eFAST integrates the merits of FAST and Sobol's algorithms, representing a method with high efficiency and precision, and the ability to adequately compute interaction effects among parameters (Iooss and Lemaître, 2015).

The GSA has several aspects that can affect sensitivity indices and their uncertainty, regardless of the method adopted. In general, the most important uncertainty sources of GSA are: (i) sample size (SZ), (ii) range of parameters and (iii) complexity of the model (Xu and Gertner, 2011; Gan et al., 2014; Song et al., 2015). To our knowledge, there are no studies in the literature investigating effects in GSA caused by SZ and parameter range on eFAST method in sugarcane models. In the case of parameter SZ, studies that utilized the eFAST method were based on the evidence presented by Wang et al. (2013), where the recommended SZ is greater than 1024. However, in studies that did not follow this recommendation, a very large SZ was adopted without clear criteria (Tan et al., 2016). By adopting an SZ like that of Wang et al. (2013), the model characteristics are ignored, as the SZ is influenced by various factors such as the number of parameters, the analysed output variable, and the methods used to estimate variables within the model. Additionally, choosing a very large SZ requires excessively long computational time, which can render the analysis infeasible in some cases.

The parameter range is another source of uncertainty in PBCM when using Bayesian approaches for parameter estimation (Makowski et al., 2006), consequently affecting the sensitivity indices in a GSA. Depending on the range of parameters, it is possible to generate calibrations that do not represent the desired genotype (He et al., 2010, 2009; Sexton et al., 2016; Marin et al., 2017) or correctly quantify the uncertainty (Soetaert and Petzoldt, 2010; Dzotsi et al., 2013; Gan et al., 2014; Zhang et al., 2020; Pereira et al., 2021). Wang et al. (2013) compared a range of parameters measured for maize and a relative range of 10% in relation to a reference calibration and found important differences in the GSA results. Li et al. (2019) evaluated different relative ranges, from 10 to 50%, in relation to a reference value and concluded that the most important parameters for the 10% range diverged from those obtained using the 50% range. Many recent PBCM studies (Vazquez-Cruz et al., 2014; Tan et al., 2016; Jin et al., 2018; Li et al., 2019) have adopted relative parameter ranges to apply a GSA, which can result in serious methodological errors. According to Homma and Saltelli (1996), the GSA principle is to identify the parameters that cause the greatest uncertainty in the model, and this is not possible when all parameters are disturbed to create relative parameters range. Finally, there is still the model complexity, which is a source of considerable uncertainty and a complicated issue to be considered in the GSA (Razavi and Gupta, 2015).

The SAMUCA crop model has been relatively little evaluated compared to well-established models such as DSSAT/CANEGRO and APSIM-Sugar (Thorburn et al., 2005; Marin and Jones, 2014; Marin et al., 2015; Sexton et al., 2017). Pereira et al. (2021) used the PRCC to perform a GSA to identify the most important parameters for SAMUCA and then to explore the model uncertainty. However, some issues remained unclear due to the following limitations: (i) the GSA was performed only for the end-of-season output values; (ii) the PRCC method is limited when the parameter relation is not monotonic, and only four out of the 24 genotype parameters are monotonic in SAMUCA, which means that some parameter responses might be neglected during the GSA; (iii) authors did not consider the effect of the range of parameters on GSA results. The first two limitations may result in the omission of important parameters because they affect certain variables at a different time of simulation (Lamboni et al., 2009; Dejonge et al., 2012), and the third has never been investigated for sugarcane crop models, and even for other crops the range of parameter effect on GSA was little studied (Wang et al., 2013).

Considering the aspects mentioned above, this paper aimed to: (i) determine the optimal SZ for the eFAST method; (ii) investigate whether there is a difference between the ranges of parameters used in GSA; (iii) identify which parameters are responsible for the greatest uncertainty in the SAMUCA model.

Material and methods

SAMUCA crop model

The SAMUCA model is a process-based model that simulates the growth and development of sugarcane. Originating from the argument of Sinclair and Seligman (Sinclair and Seligman, 1996), which emphasizes the importance of research groups developing their own models to incorporate region-specific processes, this approach enhances the understanding of the uncertainties and limitations of the model. Thus, the SAMUCA model was mainly developed by Marin and Jones (2014) and Marin et al. (2017), using a large database for different Brazilian production conditions. Subsequently, Vianna et al. (2020) improved the model structure by decreasing the uncertainty in the soil water balance and including the effect of straw cover on sugarcane growth and development, modifying routines of soil moisture and the flow of water and heat from soil, compared to the previous versions. The most recent version is also included in the DSSAT platform v.4.8 (https://dssat.net/).

The standalone version of SAMUCA model previously presented by Marin and Jones (2014) simulates fibre and sucrose biomass partitioning at the internode level. The new version keeps the same implementation; however, leaves and internodes are grouped into phytomer units (Vianna et al., 2020). The phytomers, which are associated with tillers that have emerged in the simulation area. The root system is treated as a single set of organs and develops as thermal age progresses. Plant development is calculated separately for roots, tillers, leaves and internodes, considering the accumulation of thermal units. Soil temperature is used for underground processes such as root development and the development of tiller phytomers that have not yet emerged. Air temperature, on the other hand, is used to calculate the development of above-ground phytomers. The phytomer appearance rate for each tiller is simulated based on the plastochron interval, which is a measure of thermal time required for the development of each phytomer. During the initial phase of sugarcane development, the tiller appearance rate is calculated based on soil temperature and the emergence interval of new tillers, known as the tiller emergence rate. The senescence phase of tillering is triggered when the amount of light transmitted through the canopy reaches a specific threshold. In this phase, younger tillers gradually senesce over thermal time until the number of tillers reaches the expected final population, determined by the parameter popmat. The leaf senescence rate is controlled by the parameter maxgl, which defines the maximum number of leaves per tiller. When the number of leaves on a tiller exceeds this value, the oldest leaf is shed as it has completed its life cycle. Soil temperature is used to calculate thermal age and temperature-dependent processes until the organ in question (such as below-ground internodes, roots and tillering) emerges at the soil surface. After emergence, all new leaves use air temperature to calculate accumulated thermal units and contribute to the formation of the crop canopy, expressed by the leaf area index. As for internodes, the switch from using soil temperature to air temperature for the development of new internodes occurs when the stalk emerges at the soil surface, controlled by the parameter n_lf_when_stk_emerg. Internodes are considered ‘appeared’ at the top of the tiller after a minimum number of leaves have emerged for their formation, while older leaves are attached to already differentiated internodes, enclosed in leaf sheaths. At the apex of the tiller, a set of developing leaves surrounds the apical meristem of the bud, at the centre of the leaf whorl. These leaves will become visible at the tip of the tiller and eventually develop into internodes.

Data and management

The crop growth simulation scenario was based on a field experiment conducted in the College of Agriculture ‘Luiz de Queiroz’, Piracicaba, São Paulo (Lat: 22°41′55″S, Lon: 47°38′34″W, Alt: 540 m). The sugarcane cultivar was the RB86-7515, a widely used genotype in Brazil (ca. 30% of Brazil's planted area). It was planted on 16 October 2012, with a row spacing of 1.4 m and depth of 0.2 m. The BARE treatment was conducted over four consecutive years, while the GCTB treatment commenced in the first ratoon (October 2013) and was carried out for 3 years. Agricultural practices were adopted to represent high yield farming systems and to ensure the crop was free from pests, diseases and nutritional stress. The climate is characterized by hot and humid summers with dry winters (Cwa-Koppen classification), and the soil classified as Typic Hapludox. The experiment was irrigated by a centre-pivot, based on monitoring the soil moisture by Frequency Domain Reflectometry (FDR) and the evapotranspiration by Bowen Ratio Method (BRM) in both treatments (Table 1).

Table 1. Description of seasons, planting and harvesting dates, duration of the season in days, treatments and variable measurements of the field experiment in Piracicaba, Brazil

Green cane trash blanket (GCTB), stalk dry mass (SDM; t/ha) and stalk fresh mass (SFM; t/ha) of leaf area index (LAI; m²/m²), sucrose concentration in fresh matter (POL; %) and tiller population (TIL; #/m²).

As the eFAST method requires a high computational time, we divided our study into two steps. In the first step (STp1), we simulated the first ratoon of sugarcane with GCTB and BARE, testing different SZ and two sets of parameter ranges (Table 2). The main objective of STp1 was to define what minimum SZ is needed to obtain a reliable GSA and then apply them to the second step (STp2) of the study. In STp1 we ran the GSA using end-of-season values as reference of the variables: stalk dry mass (SDM), stalk fresh mass (SFM), leaf area index (LAI), sucrose concentration in the fresh matter (POL) and tiller population (TIL). In the STp2, after determining the SZ where the standard deviation of the sensitivity indices tends to zero and there is no risk of altering the parameter importance order, we performed a long simulation considering the different ratoons and ran the GSA in function of the daily values simulated in the whole season of each variable (SDM, SFM, POL, LAI and TIL).

Table 2. Description of the processes performed in the first step (STp1) and the second step (STp2); the simulated season, the sample size (SZ) evaluated, soil cover type, number of repetitions (NR) and parameter range set (PRS)

The details about the SZ are presented in the section on sample size, and the details about PRS1 and PRS2 are provided in Table 3.

Sensitivity analysis

Extended Fourier Amplitude Sensitivity Test

The eFAST is an algorithm that combines two GSA methods: the FAST and the Sobol (Saltelli et al., 1999, 2010), which in turn, use the model output variance principle. While FAST can scan the entire parameter space and obtain quantitative sensitivity measures in terms of the main sensitivity index (S _i) of each parameter to output variance, the Sobol calculates the total sensitivity index (ST _i) and provides an indication of the overall effect of a given parameter, considering all possible interactions of that parameter with others (Sobol, 2001). The eFAST is based on the decomposition of the model's output variance, determining which fraction of the variance explained by each parameter. This variation is quantified using the statistical notion of variance (analogous to ANOVA):

(1)$$\sigma ^2 = \mathop \sum \limits_{i = 1}^N \displaystyle{{( {y_i-\bar{y}} ) ^2} \over {( {N-1} ) }}$$

where N is number of models runs, y _i is i _th model output and $\bar{y}\;$sample mean. Partitioning of variance in eFAST works by varying different parameters at different frequencies, encoding the identity of parameters in the frequency of their variation. In recent years, due to these advantageous properties, eFAST has become more popular in hydrological, ecological and agronomy modelling (Varella et al., 2010; Reusser et al., 2011; Xing et al., 2017; Li et al., 2019). We implemented the eFAST method for the SAMUCA model in the sensitivity R-package available at: https://cran.r-project.org/web/packages/sensitivity/index.html.

The main sensitivity index (S _i) of a given parameter (i) is calculated as the variance at a particular parameter's unique frequency (and harmonics of that frequency) divided by total variance (VAR _t). First, variance (VAR _i) is calculated from the Fourier coefficients at the frequency of interest (j):

(2)$$VAR_i = 2( {A_j^2 + B_j^2 } ) $$

(3)$$A_j = \mathop \smallint \limits_{-\pi }^\pi f( s ) cos( {\,js} ) \;ds$$

(4)$$B_j = \mathop \smallint \limits_{-\pi }^\pi f( s ) sin( {\,js} ) \;ds$$

(5)$$S_i = \displaystyle{{VAR_i} \over {VAR_t}}$$

where s is a scalar variable within the range −∞ < s <   + ∞; A _j and B _j are the Fourier coefficients (or Fourier amplitude) over the domain of integer frequencies j ∈{ − ∞, …, − 1, 0, 1, …, ∞}. S _i represents the fraction of the output variance of the model explained by the input variation of a given parameter. ST _i is calculated as the remaining variance after the complementary set contribution is removed. Thus, to estimate ST _i for the given parameter i, the eFAST algorithm first calculates the sensitivity indices except for parameter i using the identification frequencies.

(6)$$ST_i = \displaystyle{{VAR_t-VAR_{{-}i}} \over {VAR_t}}$$

where VAR _−i is the sum of all the variance terms that do not include the parameters i.

S _i and ST _i must vary between 0 and 1, and effects are greater the indices reach values closer to 1, whereas values close to 0 indicate negligible effects. ST _i considers both S _i and the interactions between the parameters, such interactions be evaluated by the difference between the ST _i and the S _i. The two sensitivity indices S _i and ST _i are equal if the effect of the parameter i 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 additive into parameter i. We considered only the parameters that had S _i > 0.05 and ST _i > 0.1 as significant and relevant in the GSA analysis, as suggested by Dejonge et al. (2012) and Xing et al. (2017). This allowed us to identify the parameters that individually explain at least 5% of the model's variance through the S _i index. Additionally, using the ST _i index, we could identify the parameter that, when combined with others, explains more than 10% of the variance, indicating that the contribution of the other parameters explains at least an additional 5% of the variance.

Sample size

The SZ refers to the number of samples generated within the range of each parameter. In other words, for a given parameter x, SZ samples are generated while respecting the limits of variation for that parameter. This SZ is a key factor in determining the number of model evaluations (Ne), where the number of evaluations is obtained by multiplying the number of evaluated parameters (Np) by the SZ (Ne = Np × SZ). In some cases, a large number of model evaluations are required, which can restrict the method's use. Therefore, the relationship between SZ and the convergence of the sensitivity measure is of utmost importance. The SZ is determined by the sampling frequency and the exploration curve, as described by Saltelli et al. (1999). To investigate the influence of SZ on the convergence of sensitivity indices, an SA was performed with different SZ. Seven SZ cases were used: 65, 129, 257, 513, 1025, 2049 and 4097, as proposed by Wang et al. (2013). For each SZ, we conducted ten repetitions considering the 1st ratoon for the BARE and GCTB treatments. We employed two criteria to determine the SZ. The first criterion involved calculating the sum of S _i values and checking if this sum converged to 1. If it did not meet this criterion, the SZ would be discarded. The second criterion involved calculating the mean and standard deviation of S _i values based on the ten repetitions. In this way, we identified the smallest SZ where the standard deviation of S _i was sufficiently small and the sum of S _i approached 1, classifying it as the most suitable for the variable analysis. Additionally, the smallest SZ value that simultaneously met these two criteria for all variables was defined as the most appropriate for the SAMUCA model. In this analysis, we considered two ranges of parameter variation, which will be described in the next item (Parameters range set).

Parameters range set

The GSA is affected by the uncertainty range of the parameters (Wang et al., 2013), so we investigated this factor for the GSA of the SAMUCA model by constructing two ranges of genotype parameters. The first parameter range set (PRS1) was constructed based on the literature, containing the maximum and minimum values of each genotype parameter, regardless of the sugarcane cultivar (Pereira et al., 2021). For the second parameter range set (PRS2), we considered the studies of Wang et al. (2013) and Jin et al. (2018), and concluded that the order of importance of the parameters converges to above 10% disturbance. Thus, in order not to use an excessively small disturbance that would cause inconsistency in the GSA results, and to avoid an excessively large disturbance to generate parameter values outside their genotype reality, we choose to cause a ±25% perturbation in the values calibrated by Vianna et al. (2020). In Table 3, we present the description of the parameters and their respective values for PRS1 and PRS2.

Table 3. Description of parameters, calibrated values for genotype RB867515 (μ) by Vianna et al. (2020), first set of parameters (PRS1) extracted from Pereira et al. (2021), based on the literature, and second set of parameters (PRS2) with ±25% perturbation relative to μ; CV is the coefficient of variation with respect to μ; Max and Min values are range used for random parameters uniform distribution

Results

Sample size

Figures 1 and 2 presented the S _i values for the main parameter obtained from the PRS1 intervals (Fig. 1) and PRS2 intervals (Fig. 2), in relation to different SZ in the BARE and GCTB treatments. We observed a dispersion for small SZ. The SZ of 65 and 127 showed high variability both in the sensitivity indices and in the order of the main parameters explaining the model variance. For the lowest SZ evaluated (65 and 127), it was not possible to obtain sensitivity indices and accurately quantify the order of parameter importance due to the divergence of the sensitivity indices' sum, S _i and ST _i from 1. In these cases, the eFAST method was unable to quantify the sensitivity indices for the PRS1 (Table 4) and PRS2 (Table 5) sets.

Figure 1. Evolution of sensitivity index of the most important parameter for BARE and GCTB with increasing SZ for variables mass stalk dry (a) and fresh (b), sucrose content (c), tiller population (d) and leaf area index (e), for PRS1; red line is the average of the ten simulations for each SZ and in blue we have the max and min S _i of each SZ.

Figure 2. Evolution of sensitivity index of the most important parameter for BARE and GCTB with increasing SZ for variables mass stalk dry (a) and fresh (b), sucrose content (c), tiller population (d) and leaf area index (e), for PRS2; red line is the average of the ten simulations for each SZ and in blue we have the max and min S _i of each SZ.

Table 4. Order of importance of the genotype parameters for the PRS1 interval, considering the main sensitivity index (S _i) and treatment bare soil (BARE) and green cane trash blanket (GCTB)

It was not possible to define an order of importance for the parameters for the sample size (SZ) values that do not appear in the table; stalk dry mass (SDM), stalk fresh mass (SFM), leaf area index (LAI), sucrose concentration in the fresh matter (POL) and tiller population (TIL).

Table 5. Order of importance of the genotype parameters for the PRS2 interval, considering the main sensitivity index (S _i) and treatment bare soil (BARE) and green cane trash blanket (GCTB)

It was not possible to define an order of importance for the parameters for the sample size (SZ) values that do not appear in the table; stalk dry mass (SDM), stalk fresh mass (SFM), leaf area index (LAI), sucrose concentration in the fresh matter (POL) and tiller population (TIL).

We noticed that the SZ varied according to the analysed variable; however, we did not observe any effect of the GCTB or BARE treatments on the SZ (Tables 4 and 5). Tables 4 and 5 present the order of significant parameters for each SZ, considering different parameter sets and treatments. For the SDM variable, the required SZ for the PRS1 set was 1025 (Table 4), which was smaller than the 2049 obtained for the PRS2 set (Table 5). In the case of the SFM, POL and LAI variables, the same SZ were found for both the PRS1 and PRS2 sets: 1025 for SFM and 513 for POL and LAI (Tables 4 and 5). The TIL variable achieved the sensitivity indices with the smallest SZ in both parameter sets. For the PRS1 set, the SZ was 257 (Table 4), while for the PRS2 set, it was 513 (Table 5). We observed a divergence among the variables, and the minimum SZ required to calculate the sensitivity indices using the eFAST method varied depending on the variable. Consequently, the lowest SZ at which it was possible to calculate the sensitivity indices for all variables was 2049 (Tables 4 and 5).

Crop features

We consider two key points to evaluate; (i) the order of importance of the parameters and (ii) the values of S _i and ST _i in both parameter ranges (PRS1 and PRS2). We noticed that GCTB affected the order of the main parameter only in PRS2, while S _i and ST _i values were slightly affected between GCTB and BARE treatments.

In PRS1, the results between BARE and GCTB were similar since the main parameter for all variables was maintained regardless of the presence or absence of GCTB (Fig. 3). For example, the main parameter for SDM and SFM was the plastochron, and it was responsible for explaining 29 and 31.6% of the variance in the BARE and GCTB treatments, respectively; for SFM the explained variance was 27.9 and 30.4% for BARE and GCTB, respectively (Tables 6 and 7). The GCTB only influenced the order of importance of the parameters of the variable SDM. In the BARE treatment the parameter n_lf_when_stk_emerg was the 6th parameter (explaining 5.4% of the variance) and in the treatment with GCTB it was the 3rd (explaining 12.7% of the variance) (Tables 6 and 7). Furthermore, in the SFM variable, the plastochron parameter appears to have greater relevance in the presence of GCTB compared to BARE, as the percentage explained only by it increased, decreasing the values of the parameters mid_tt_it_growth and max_it_dw (Tables 6 and 7).

Figure 3. Main sensitivity index (S _i) and total sensitivity index (ST _i) of the first parameter range set apply to bare soil (PRS1-BARE) and green cane trash blanket (PRS1-GCTB) treatments, as well as for the second parameter range set applying to bare soil (PRS2-BARE) and green cane trash blanket (PRS2-GCTB) treatments; the sensitivity analysis was calculated for the end-of-season value of each output variable.

Table 6. Relative value of the model output variance (σ ²) explained individually by each parameter, and the variance sum ($\mathop \sum \nolimits^ {\sigma \,^2}\;)$ the parameters; we only considered the parameters that presented S _i > 0.05 and the sample size (SZ) of 2049 in treatment bare soil (BARE) to first set of parameters (PRS1) and second set of parameters (PRS2); stalk dry mass (SDM), stalk fresh mass (SFM), leaf area index (LAI), sucrose concentration in the fresh matter (POL) and tiller population (TIL)

Table 7. Relative value of the model output variance (σ ²) explained individually by each parameter, and the variance sum ($\mathop \sum \nolimits^ {\sigma \,^2}\;)$ of the parameters; we only considered the parameters that presented S _i >0.05 and the sample size (SZ) of 2049 in treatment green cane trash blanket (GCTB) to first set of parameters (PRS1) and second set of parameters (PRS2); stalk dry mass (SDM), stalk fresh mass (SFM), leaf area index (LAI), sucrose concentration in the fresh matter (POL) and tiller population (TIL)

In PRS2, the presence of GCTB affected the results regarding the analysis of the variables SDM, SFM and LAI, while there was no influence of GCTB on TIL and POL (Fig. 3). The main parameter for SDM in the BARE treatment was popmat (explaining 33.9% of the variance; Table 6), while in the GCTB treatment, the main parameter was plastochron (explaining 29.0% of the variance; Table 7). Similarly, in the SFM variable, the main parameter was popmat in the BARE treatment (explaining 26.6% of the variance; Table 6) and plastochron in the GCTB treatment (explaining 30.5% of the variance; Table 7). In relation to LAI, the presence of GCTB led to the inclusion of the parameter mla among the significant parameters, going from three parameters in BARE (maxgl, popmat and mid_tt_lf_growth; Table 6) to four parameters in the GCTB treatment (maxgl, popmat, mid_tt_lf_growth and mla; Table 7).

Effects of parameter range

For SDM, considering the SZ of 2049, we observed a reduction in the number of significant parameters between PRS1 and PRS2. In PRS1 there were six significant parameters (plastochron, max_it_dw, n_lf_it_form, eff, popmat and n_lf_when_stk_emerg), which together were responsible for explaining 81.6% (Table 6) and 82.2% (Table 7) of the variance in the BARE and GCTB treatments, respectively. In the case of PRS2, only four parameters (popmat, plastochron, max_it_dw and eff) were significant, responsible for explaining 70% (Table 6) and 70.1% (Table 7) of the variance in the BARE and GCTB treatments, respectively.

Among the output variables analysed, SFM was the one with the highest number of significant parameters, from six to seven parameters. In PRS1 there were six significant parameters in the BARE treatment (plastochron, mid_tt_it_growth, max_it_dw, end_tt_it_growth, eff and n_lf_it_form), which explained 96.8% (Table 6) of the variance, while in the GCTB treatment there were seven parameters (plastochron, max_it_dw, mid_tt_it_, n_lf_when_stk_emerg, end_tt_it_growth, eff and n_lf_it_form), which explain 93.3% of the variance (Table 7). In the case of PRS2, there were six significant parameters in both treatments, responsible for explaining around 93.1% (Table 6) and 95% (Table 7) of the variance in the BARE and GCTB treatment, respectively. Two results should be highlighted: (i) in the BARE treatment, popmat was not a significant parameter in PRS1, while in PRS2 it was the most important parameter (Table 6); (ii) on the other hand, in the GCTB treatment, popmat remained irrelevant in PRS1 and was the second most important parameter in PRS2; the most important parameter was plastochron (Table 7).

In the case of the POL, there was no change in the group of significant parameters or in their order, with the parameters mid_tt_it_growth and end_tt_it_growth being the only significant ones, regardless of the treatment and the PRS1 and PRS2 set (Tables 6 and 7). However, we noticed that mid_tt_it_growth in PRS1 explained 65.8% (Table 6) and 64.2% (Table 7) of the variance, while in PRS2 it was 39.2% (Tables 6 and 7). For the TIL variable, there was no difference between PRS1 and PRS2, being popmat responsible for explaining 99.8% of the variance regardless of treatment.

The variable LAI in PRS1 presented only two significant of treatment (BARE or GCTB); the two main parameters were popmat and mla (Tables 6 and 7). However, in PRS2 the number of significant parameters increased, to 3 (maxgl, popmat and mid_tt_lf_growth; Table 6) in BARE and to 4 (maxgl, popmat, mid_tt_lf_growth and mla; Table 7). The biggest discrepancy is in the variance explained by the significant parameters, where in PRS1 they were 33.9% (Table 6) and 34.6% (Table 7), and in PRS2 they were 55.2% (Table 6) and 60% (Table 7).

Time-dependent effects on global sensitivity analysis

Having defined the appropriate SZ (2049), we performed an analysis considering a temporal variation in the different stages of the sugarcane crop season (plant cane to 3rd ratoon). We observed that the GSA over time provided two important results: (i) the influence of the sugarcane season did not affect the order of importance of parameters. The order of parameters was not changed compared to ratoons; for example, for SDM we noticed that the same parameters have the same degree of importance in all ratoons (Fig. 4). The same pattern was repeated for all variables and their respective parameters and (ii) at different times of sugarcane growth, there were different parameters responsible for explaining the greater variance of the model.

Figure 4. Main sensitivity index (S _i) and total sensitivity index (ST _i) calculated for first parameter set (PRS1) and second parameter set (PRS2) for bare soil (BARE) and green cane trash blanket (GCTB) treatments across plant cane, 1st, 2nd and 3rd ratoons for stalk dry mass (SDM). The red dashed line separates the stages from plant cane to the 3rd ratoons; ‘days after planting’ represents the time in days from plant cane to the end of 3rd ratoons.

For the variable SDM, during the different growth phases, there was a change in the order of the main parameter. For example, during the crop establishment and development phase, the plastochron parameter was responsible for explaining more than 50% of the SDM variance, becoming the main parameter in this period when using the PRS2 set. On the other hand, at the end of the maturation phase, the plastochron was the second most important parameter, behind the popmat with the PRS2 set (Fig. 4 – PRS2 BARE and GCTB). In PRS1, during the emergence phase until mid-establishment, the main parameter was n_lf_it_form, later the plastochron became the main parameter until the end of the growth phase (Fig. 4 – PRS1 BARE and GCTB). Furthermore, in short periods and isolated in the growth phase, the parameters n_lf_it_form, end_tt_lf_growth and n_lf_when_stk_emerg were significant using the PRS2 set (Fig. 4 – PRS2 BARE and GCTB); these three parameters were not reported in Tables 6 and 7 as they were not significant at the end of the growth phase. This pattern was not observed in PRS1 (Fig. 4 – PRS1 BARE and GCTB), that is, all significant parameters presented in Tables 6 and 7 appeared in the analysis throughout the growth phases.

For SFM, plastochron had a greater impact on the simulation regardless of the parameter set (PRS1 and PRS2) or treatment (Fig. 5 – BARE and GCTB). Among the sets of parameters, PRS1 and PRS2, there was the inclusion of the popmat parameter in PRS1 in the BARE treatment; however, its influence was minimal (5%). For PRS2, the popmat remained relevant during most parts of the growth phase, but in our view not enough to have more impact on the SFM simulations. The parameter end_tt_lf_growth showed significance at specific moments during the tillering phase for the PRS2 set, in both the BARE and GCTB treatments (Fig. 5 – BARE and GCTB), a situation not observed in the PRS1 set. Furthermore, even within the PRS2 set, its relevance was more pronounced in the GCTB treatment, as it was evident in all three ratoons.

Figure 5. Main sensitivity index (S _i) and total sensitivity index (ST _i) calculated for first parameter set (PRS1) and second parameter set (PRS2) for bare soil (BARE) and green cane trash blanket (GCTB) treatments across plant cane, 1st, 2nd and 3rd ratoons for stalk fresh mass (SFM). The red dashed line separates the stages from plant cane to the 3rd ratoons; ‘days after planting’ represents the time in days from plant cane to the end of 3rd ratoons.

For the variable POL, we had the inclusion of several parameters depending on the of the analysed growth phase (Fig. 6). In this case, in PRS1, in addition to mid_tt_it_growth and end_tt_it_growth, there was the inclusion of plastochron, n_lf_when_stk_emerg and n_lf_it_form (Fig. 6 – PRS1 BARE and GCTB). In PRS2, the parameters plastochron, n_lf_when_stk_emerg, n_lf_it_form and end_tt_lf_growth were included (Fig. 6 – PRS2 BARE and GCTB). For PRS1, the parameter end_tt_it_growth was not significant in the 1st ratoon (Fig. 6 – PRS1 BARE and GCTB), evidencing a combination of the range of parameters with the weather conditions for that season. Yet, for PRS2, end_tt_lf_growth was significant, while in PRS1 the same was not observed.

Figure 6. The main sensitivity index (S _i) and total sensitivity index (ST _i) calculated for first parameter set (PRS1) and second parameter set (PRS2) for bare soil (BARE) and green cane trash blanket (GCTB) treatments across plant cane, 1st, 2nd and 3rd ratoons for sucrose content (POL). The red dashed line separates the stages from plant cane to the 3rd ratoons; ‘days after planting’ represents the time in days from plant cane to the end of 3^rd ratoons.

The TIL variable had two main parameters, from the beginning to the middle of the season, that was tillochron and from the middle to the end of the season popmat (Fig. 7). In the first half of the season, the tillochron influence was around 90%, while in the second half the popmat explained more than 90% of the variance. In the transition period between these two phases of the crop season, that takes place between the end of the establishment and half of the development of the crop, we observed that there was a punctual influence of some other parameters, such as mla and plastochron. However, during this transition, there were strong indications of the effect of climate on GSA, as these two parameters (mla and plastochron) were not in all ratoons. To corroborate this result, we did not observe relevant influences of the parameter range (PRS1 or PRS2) and the evaluated treatment (BARE or GCTB).

Figure 7. The main sensitivity index (S _i) and total sensitivity index (ST _i) calculated for first parameter set (PRS1) and second parameter set (PRS2) for bare soil (BARE) and green cane trash blanket (GCTB) treatments across plant cane, 1st, 2nd and 3rd ratoons for tiller population (TIL). The red dashed line separates the stages from plant cane to the 3rd ratoons; ‘days after planting’ represents the time in days from plant cane to the end of 3rd ratoons.

The LAI variable, contrary to what was seen in sections ‘Crop features’ and ‘Effects of parameter range’, with a maximum number of significant parameters of four at the end of growth (Tables 6 and 7), was influenced by 15 parameters at different growth phases (Fig. 8). However, the high number of parameters can be attributed to isolated events. These events occurred at specific moments during the growth phase but did not follow a consistent pattern or a sequential period of influence. In PRS1, we observed that tillochron, plastochron, mla and init_leaf_area parameters had greater relevance, as they explained, individually, more than 40% of the variance of LAI at different moments of emergence phase and of establishment phase (Fig. 8 – PRS1 BARE and GCTB). In the establishment and maturity stages, the variance of the LAI limited to the mla and popmat parameters with the PRS1 set and, in our view, the mla parameter was the most important parameter as it was significant in all growth phases (Fig. 8 – PRS1 BARE and GCTB). We arrived at this result because the mla parameter explains between 20 and 50% of the LAI variance throughout the entire growth phase, while the popmat parameter contributes approximately 10–30% of the variance and has influence over a shorter time interval. Thus, even though the variance explained by popmat and mla may be similar at the end of the cycle, the impact of mla is greater and more consistent throughout the entire growth phase. In PRS2, the emergence and establishment phases were dominated by plastochron and mla, and in the development and maturation phases, the parameters maxgl , popmat and mid_tt_lf_growth were responsible for explaining most of the variance (Fig. 8 – PRS2 BARE and GCTB). In our view, considering the PRS2 set, the most important parameter was maxgl because it explained more than 30% of the variance in almost all growth phases. On the other hand, popmat and mid_tt_lf_growth explained 10–20% of the variance, and over a shorter period compared to maxgl. This result showed the uncertainty present in the range of the chosen parameters, since maxgl was not significant at growth phases in the PRS1 set (Fig. 8 – PRS1 BARE and GCTB).

Figure 8. The main sensitivity index (S _i) and total sensitivity index (ST _i) calculated for first parameter set (PRS1) and second parameter set (PRS2) for bare soil (BARE) and green cane trash blanket (GCTB) treatments across plant cane, 1st, 2nd and 3rd ratoons for leaf area index (LAI). The red dashed line separates the stages from plant cane to the 3rd ratoons; ‘days after planting’ represents the time in days from plant cane to the end of 3rd ratoons.

Discussion

The SZ can vary depending on the variable of interest. For example, the variable POL required a minimum size of 513 for convergence (Tables 4 and 5), while the variables SDM and SFM required an SZ of at least 2049 (Tables 4 and 5). Additionally, we observed that for the TIL variable, the SZ was 257 for PRS1 and 513 for PRS2 (Tables 4 and 5). This difference implies that computational time can be reduced depending on the output variable. Therefore, the high computational cost, which is one of the major limitations of GSA (Jeuffroy et al., 2006; Marino et al., 2008; Gilardelli et al., 2018), can be optimized through previous studies exploring the SZ. Thus, the most suitable SZ for conducting SA in the SAMUCA model was 2049, as it encompasses all the analysed variables. This result is in line with the findings reported by Wang et al. (2013), where an SZ of 2049 yielded the most stable sensitivity indices.

As far as we know, Wang et al. (2013) is the main study that used the eFAST method in a PBCM and evaluated different SZ, being considered the main reference source for other GSA studies. Another study, conducted by Ma et al. (2017), also determined the most suitable SZ for their model and found lower values compared to those observed by Wang et al. (2013) and obtained in our study. This highlights the uniqueness of GSA, where the adopted boundary conditions and the analysed model have an impact on the convergence of sensitivity indices. Our study has two distinct characteristics compared to the works of Ma et al. (2017) and Wang et al. (2013): we evaluated the impact of SZ on different variables and considered two sets of parameter ranges. We observed that, in addition to the output variable, the SZ differed depending on the applied parameter sets (PRS1 and PRS2; Tables 4 and 5). Although there are several studies (Xu and Gertner, 2011; Gan et al., 2014; Song et al., 2015) that consider the effect of SZ and parameter range on sensitivity indices, there are currently no available studies on PBCM that consider the relationship between SZ in different analysis variables and the parameter range. Therefore, we can state that the SZ in GSA may vary depending on (i) the analysed output variable and (ii) the parameter ranges.

As in many studies on the influence of the environment on sensitivity indices (Dejonge et al., 2012; Sexton et al., 2017; Zhang et al., 2020), there was an effect of GCTB on SAMUCA sensitivity indices, but it was not sufficient to change the order of importance of the parameters and it did not influence the parameter range (PRS1 and PRS2) (Figs 4–8). For example, the SFM had the plastochron parameter as the main parameter in PRS1 in both treatments (BARE and GCTB), and in PRS2, the main parameter was plastochorn in BARE treatment and popmat in GCTB treatment. To consider that a GCTB was sufficient to change the order of importance of the parameters, an alteration between the parameters should be observed in both sets of parameters (PRS1 and PRS2; Tables 6 and 7). In the case of GCTB, it has a direct influence on the soil heat flux and on the soil water dynamics (Pereira et al., 2021; Vianna et al., 2020), which would hypothetically refer to the lower effect of GCTB on GSA for the variables LAI, POL, SDM, SFM and TIL; if we would evaluate soil temperature or evapotranspiration, we would possibly have observed a greater impact of GCTB, as it is directly related to these variables. We assume this based on the results of Dejonge et al. (2012), who identified the radiation use efficiency as the most important parameter for yield, both for irrigated and rainfed treatments. However, for crop evapotranspiration, in the irrigated environment, the main parameters were related to the crop, while in the rainfed environment the main parameters were related to the soil (Dejonge et al., 2012). Thus, the effect of management on the SA is dependent on the variable of interest and did not affect all simulation variables in our study.

The range of parameters was the main source of uncertainty for the SA, changing the order of parameters for the variables SDM, SFM and LAI. Li et al. (2019) have already identified that parameter range affected the order and variance explained by each parameter. Many studies have applied parameter ranges with relative values from GSA studies in different PBCM (Vazquez-Cruz et al., 2014; Tan et al., 2016; Jin et al., 2018; Li et al., 2019), which, in our view, may not be the most appropriate approach, based on the influence of the parameter range. According to Santelli et al. (1999), the principle of the GSA is to quantify the model uncertainty based on perturbations on the environmental variables. In this sense, the use of relative values in the parameter ranges without experimental foundations may not represent well the influence of the parameters of a PBCM. For example, for the variable LAI, mla and maxgl were observed as the main parameters in PRS1 and PRS2, respectively. The CV of these parameters was 29% (mla) and 9% (maxgl), in PRS1, and 25% for both parameters in PRS2; for this configuration of PRS1, maxgl was not significant, but it was the most important parameter in PRS2. Yet, this can be related to the correlation that exists between the parameters, neglected in many Bayesians and PBCM applications (Marin et al., 2017; Pereira et al., 2021), which demonstrated that any uncertainty analysis using relative values seems to be inadequate.

In SA, it is necessary to use a time interval for one evaluation, as changes in the order of importance of main parameters occurred during the growth of the crop. For example, considering our final SFM yield value, the parameters popmat and plastochron accounted for between 21.8 and 30.4% of the variation, respectively (Table 7). However, the popmat is less relevant than the plastochron, as the plastochron influenced a longer period of the season, being responsible for explaining up to 60% of the variance in the first 150 days of each ratoon in the PRS2 set (Fig. 5, BARE and GCTB). For the TIL variable, this was even more evident, having distinct parameters and explaining over 90% of the variation in different times. In addition, the occurrence of significant different parameters between ratoon years indicated a possible influence of climate variation (Fig. 7). Several GSA studies have demonstrated the impact of seasonality on sensitivity indices, regardless of the method employed (Vazquez-Cruz et al., 2014; Tan et al., 2016; Sexton et al., 2017; Xing et al., 2017; Li et al., 2019; Attia et al., 2021). In our study, we conducted GSA using a daily time interval from plant cane to the 3rd ratoon and observed that the ranking of parameter importance between ratoons remained consistent (Figs 4–8). Over a short span of 4 years, we did not detect any impact of climate variability on the sensitivity indices. However, there is evidence that climate variability influenced the sensitivity index, as noted by Anderson et al., 2014; Attia et al., 2021; Gilardelli et al., 2018; Sexton et al., 2017. Therefore, for sugarcane models, it is crucial to incorporate time intervals in GSA, and future studies should consider the influence of climate in GSA.

Conclusion

The results showed that the SZ and parameter range were important for GSA, and that an SZ of at least 2049 was required for the sensitivity indices to converge regardless of the variable. However, some variables needed smaller sizes, such as the case of TIL that predicted an SZ of 257. The range of parameters must be carefully investigated, and we demonstrate that the use of relative values without biophysical basis to determine the parameter ranges is inappropriate for the uncertainty analysis, and measured thresholds should always be used, even if from different genotypes, to determine the model's response across the full range of parameters. We did not identify the influence of GCTB and the different ratoons on the order of importance of the parameters, they only slightly affected the values (S _i and ST _i) of the sensitivity indices.