Temporal Sensitivity Analysis of the MONICA Model: Application of Two Global Approaches to Analyze the Dynamics of Parameter Sensitivity

Abstract

Sensitivity analysis (SA) is often applied to evaluate the behavior of ecological models in which the integrated soil and crop processes often vary over time. In this study, the time dependence of the parameter sensitivity of a process-based agro-ecosystem model was analyzed for various sites and model outputs. We applied the Morris screening and extended FAST methods by calculating daily sensitivity measures. By analyzing the daily elementary effects using the Morris method, we were able to identify more sensitive parameters compared with the original approach. The temporal extension of the extended FAST method revealed changes in parameter sensitivity during the simulation time. In addition to the dynamic parameter sensitivity, we noticed different relationships between parameter sensitivity and simulation time. The temporal SA performed in this study improves our understanding of the investigated model’s behavior and demonstrates the importance of analyzing the sensitivity of ecological models over the entire simulation time.

1. Introduction

Complex crop growth models are widely applied to assess the productivity and sustainability of agricultural systems, as well as the effects of climate change on these systems [1,2,3,4]. In addition to simulating crop growth and development, these models also describe biophysical and biological processes in the soil-crop-atmosphere system. Because of the complex nature of the represented processes, crop models contain a large number of biophysical and physiological parameters. Moreover, crop models can be made to describe the characteristics of different crops or cultivars by using distinct parameters to adjust crop development, phenology, partitioning, etc. However, the determination of these parameter values is often difficult [5] because their acquisition via field observations is costly and time consuming. As an alternative approach, parameter values are often estimated using optimization methods, simulated annealing algorithms, genetic algorithms, or Bayesian approaches [6,7,8,9]. However, such methods can only be applied to estimate a small number of model parameters because their computational costs are proportional to the number of analyzed parameters.

Sensitivity analysis (SA) is helpful for identifying the parameters and processes that exert the most influence on the outputs of a model. SA is also applied to improve the understanding of a model’s structure and behavior. Furthermore, SA provides information about how a model behaves under different conditions [10], which can be used to guide model simplification and reduction efforts [11]. SA is an important tool for analyzing ecological models, especially for richly-parameterized crop growth models. Many studies have shown that only a small number of model parameters are responsible for most of a model’s output variability [12,13,14]. Thus, the identification of the most influential model parameters may facilitate model calibration, for example, allowing non-influential parameters to be excluded during parameter estimation.

In many studies in which SA is applied to analyze ecological models, the model responses are assessed at only a single point in time in the simulations [12,13,15,16]. Such studies neglect the fact that process-based ecological models involve time-dependent processes. Because of this time dependence, the influence of the parameters may change during the simulation time. In crop models, some parameters may be influential during earlier developmental stages but decline in influence as plant development advances, or vice versa. Moreover, the climate or soil information that is used to initialize the model simulation has an effect on parameter sensitivity [17,18,19]. Because SA at a single simulation period provides an unrepresentative view of parameter sensitivity, temporal sensitivity analysis is necessary to provide deeper insight into a model’s structure and behavior [18,20,21].

The objective of this study is to highlight the importance of temporal sensitivity analysis in obtaining an improved understanding of model behavior by demonstrating how this approach provides additional information on the sensitivity of model parameters for the whole simulation compared to the widely-used snapshot SA.

We used the agro-ecosystem model MONICA (Model for Nitrogen and Carbon in Agro-ecosystems, http://monica.agrosystem-models.com) [12,22] to analyze the time-dependent parameter sensitivity for the example of winter wheat. Because the MONICA model includes nearly 200 parameters, a two-step approach to the temporal SA was applied. First, we applied a parameter screening method to identify not only model parameters that were generally relevant, but also parameters that might be relevant only for a specific period during the model simulation. We subsequently used a global variance-based method to analyze the main and total effects of the selected parameters according to the results of the parameter screening. The calculation of the sensitivity indices was extended to analyze the changes in sensitivity that occurred throughout the simulation time. We then compared temporal SA results with the results of the original methods.

2. Materials and Methods

2.1. The Crop Growth Model

The MONICA model [12,22] is a process-based agro-ecosystem simulation model for simulating crop growth, water, and nitrogen dynamics for practical applications. A generic crop model, MONICA can be adapted to various crops by using specific model parameters that describe their physiology and development. The simulation time step is one day.

Daily net dry matter production via photosynthesis and respiration is driven by global radiation and temperature. Crop development is calculated by means of a thermal sum (degree-days) and modified, when appropriate, for each stage to account for day length and vernalization [23]. Crop growth is limited by water, nitrogen (N), heat, and oxygen stress. The processes of water and N uptake are calculated from the potential evaporation and crop N status, which depend on the simulated root distribution, as well as the availability of water and N in different soil layers [23]. The impact of N shortage on crop growth is calculated based on the concept of a critical N concentration in plants as a function of crop biomass [24].

The implemented crop growth algorithms are influenced by changes in the CO${}_2$ concentration in the air, which affects the crop’s photosynthesis rate, stomata resistance, and transpiration. The algorithms for determining the moisture concentration in the soil follow the capacity approach, enhanced with modifications introduced by [25]. This approach considers the phenomenon of capillary rise from groundwater. Evapotranspiration is calculated using the reference evaporation provided by the Penman–Monteith equation [26], adapted by means of crop-specific factors ($K_c$).

MONICA simulates the carbon cycle by calculating the population dynamics of microbes in the soil. The organic matter algorithms include the processes of mineralization, nitrification, and denitrification. The calculation of organic matter turn-over in the soil is based on the routines used in the DAISY model [27]. The soil C dynamics are described by three pairs (with slow or rapid decomposition) of conceptual pools (soil organic matter, soil microbial biomass, and added organic matter). The decomposition rate coefficients depend on both temperature and moisture and reflect the environmental conditions of the simulated site.

2.2. Study Sites and Management Settings

We used climate, soil, and management information from an energy cropping experiment (EVA) conducted in Germany. The objective of this experiment was to assess the suitability of various crop rotation schemes for a range of site conditions [28,29,30]. Five experimental sites (Figure 1) were chosen to study the effects of climate and soil on the parameter sensitivity of the MONICA model. At the beginning of the project, an extensive investigation was conducted to describe the soil profiles, based on soil samples drawn from depths of up to 2 m with 6–8 replications. The site-specific characteristic are given in

Table 1. Characteristics of the experimental sites. Soil type is specified according to the FAO classification system. C${}_{org}$ refers to the first 30 cm of the top soil layer. Total annual precipitation and mean temperature are averaged between 1971 and 2000.

Name	Ascha	Dornburg	Ettlingen	Gülzow	Werlte
Geographical location	${48}^{\circ }59$′ N	${50}^{\circ }48$′ N	${48}^{\circ }50$′ N	${53}^{\circ }42$′ N	${52}^{\circ }51$′ N
	${12}^{\circ }65$′ E	${11}^{\circ }35$′ E	$8^{\circ }41$′ E	${12}^{\circ }54$′ E	$7^{\circ }41$′ E
Height above sea level	430 m	250 m	170 m	10 m	32 m
Soil type	Stagnic Cambisol	Luvisol	Regosol-Luvisol	Planosol	Stagnic Cambisol
Soil texture	Loamy sand	Silty clayey loam	Loamy silt	Sandy loam	Loamy sand
Available water capacity	117 mm	189–215 mm	199 mm	120 mm	105 mm
C${}_{org}$	1.3%	1.0%	0.76%	0.7%	1.3%
pH value	6.4	7.2	7.3	6.4	5.2
Total annual precipitation	807 mm	474 mm	742 mm	559 mm	797 mm
Mean temperature	7.5 ${}^{\circ }$C	8.8 ${}^{\circ }$C	11.1 ${}^{\circ }$C	8.4 ${}^{\circ }$C	9.57 ${}^{\circ }$C

. The soil textures differed strongly among the sites, ranging from sandy soils in Gülzow and Werlte to loamy soils in Dornburg. The average annual precipitation varied between 474 mm in Dornburg and 807 mm in Ascha. Of the sites, Ettlingen has the highest mean air temperature (11.1 ${}^{\circ }$C), whereas Ascha has the lowest (7.5 ${}^{\circ }$C).

For the SA, we chose an artificial management setting derived from experiments performed as part of the energy cropping project [29]. The cultivation of winter wheat starting in 2005 was repeated over a four-year simulation period. The management steps and dates (e.g., ploughing, sowing, fertilization, and harvesting) were made consistent at each study site to facilitate the comparison of time-dependent sensitivity indices (

Table 2. Artificial crop management setting used for the sensitivity analysis, in which CAN was used as a calcium-ammonium-nitrate-based nitrogen fertilizer. The management steps were repeated for each cultivation year, starting in 2005. The dates of the management steps were made consistent for each study site to facilitate the analysis of the temporal SA results.

Ploughing	Sowing	Fertilization			Harvesting
3 October	4 October	31 March	70 kg N	CAN	29 July
		30 April	30 kg N	CAN
		22 May	60 kg N	CAN

). Because of the use of the same dates for sowing and harvesting, as well as the same amounts of fertilizers instead of applying individual settings for each site, the simulation results were more susceptible to stress situations, e.g., water, nitrogen, or heat stress.

2.3. Parameter Selection and Model Outputs

One hundred eighty-six model parameters were selected for inclusion in the SA. The model parameters were assumed to be independently and uniformly distributed because of a lack of information about the prior probability distributions for each parameter. The range of parameters was limited to 30% on either side of its reference value [12,18,19,31].

Six different model outputs were selected for this study, including the grain yield (grainYield), the above-ground biomass (AGB), the N content in above-ground biomass (N${}_{AGB}$), the actual evapotranspiration (ET${}_a$), the soil moisture (Moist), and the soil mineral nitrogen (N${}_{min}$), each at soil depths of 0–90 cm. The parameter sensitivity of the model was analyzed for each output separately.

2.4. Parameter Screening

The parameter screening method proposed by [32] and improved by [33] and [34] was applied. We used this approach to identify the most influential model parameters and to screen out non-influential parameters. With k as the number of input parameters, a k-dimensional p-level grid ($\mathsf{\Omega }$) was created by dividing the range of parameters into p discrete levels. The model was evaluated for r trajectories within $\mathsf{\Omega }$. The starting point of a trajectory was selected randomly. For each trajectory, every single parameter was changed separately, whereas the new point of this trajectory was an element of $\mathsf{\Omega }$. The elementary effect (EE, Equation (1)) of the i^th parameter $x_i$ with respect to the variation $\Delta$ is defined as:

$$E{E}_i=\frac{f(x_1,\dots ,x_{i-1},x_i+\Delta ,x_{i+1},\dots ,x_k)-f(x)}{\Delta }$$

(1)

where $x\in \mathsf{\Omega }$ and $x+\Delta$ is still in $\mathsf{\Omega }$, $\Delta$ as a multiple of $1/(p-1)$. EE are scaled due to different parameter ranges [34] and calculated for different trajectories within the input space. The number of model evaluations n required for the Morris method is defined by $n=r(k+1)$, where r is the number of trajectories and k is the number of input parameters. For each input parameter, an EE distribution is generated. Morris proposed two sensitivity indices, $\mu$ and $\sigma$, which represent the calculated mean and standard deviation of this distribution. A high $\mu$ indicates an input with an overall influence on the output, whereas a high $\sigma$ indicates either an input with a non-linear effect on the output or an input that interacts with other parameters. In this study, the absolute values of the EEs were used to calculate ${\mu }^{*}$, as proposed by [33], to avoid the cancellation of single EEs with opposite signs. Because the order of importance of the model parameters was irrelevant, we simply used a threshold value ${\tau }_o$ to identify a set ${\mathsf{\Psi }}_o$ of relevant model parameters:

$${\mathsf{\Psi }}_o=\{x_i|{\mu }^{*}(x_i)\ge {\tau }_o\}$$

(2)

where ${\mu }^{*}(x_i)$ is the ${\mu }^{*}$ of the $x_i$. The choice of the threshold value depends on the objective of the performed SA. A low threshold value of ${\tau }_o=0.2$ can be used to discard non-influential parameters, whereas a high threshold value should be selected to identify parameters that are highly sensitive. In this study, we selected a low threshold value to identify any parameters that somehow affect the output and to make sure that no sensitive parameters are erroneously excluded from further investigation.

The original parameter screening approach provides only a measure of the sensitivity associated with each input parameter at a single point in time (e.g., the above-ground biomass at the harvest date). When this screening method is used to identify the most influential model parameters, some parameters might be classified as negligible even though they may significantly affect some output of the model during a short period of the simulation. We therefore adopted the parameter screening method by calculating the EEs and the respective ${\mu }_t^{*}$ values at each simulation time step t to identify all model parameters that affected the output throughout the entire simulation time. As a result, we obtained a set ${\mathsf{\Phi }}_t$ as follows:

$${\mathsf{\Phi }}_t=\{{\mu }_t^{*}(x_i)\}$$

(3)

where ${\mu }_t^{*}(x_i)$ is the ${\mu }^{*}$ of the $x_i$ at simulation time t and $t\in \{1,\dots ,d_{max}\}$. For the analysis of the four-year simulation period for winter wheat, the simulation times t are defined as $t=\{1,\dots ,d_{max}\}$, where $d_{max}$ is the maximum number of simulation days. Following our approach based on the original screening method, we identified a set of relevant model parameters ${\mathsf{\Psi }}_t$ based on ${\mathsf{\Phi }}_t$ and a threshold value ${\tau }_t$:

$${\mathsf{\Psi }}_t=\{x_i|{\mu }_t^{*}(x_i)\ge {\tau }_t\}$$

(4)

A high threshold value ${\tau }_t=0.8$ has been selected to identify only parameters that are highly influential on the output at some point during simulation time. Parameters with only minor influence during the whole simulation time should be excluded from further investigation.

In this study, the parameter range for the selected parameters was divided into $p=20$ levels and analyzed for $\Delta =5/19$. Among 500 randomly-generated Morris trajectories, we chose r = 40 trajectories with the highest spread to guarantee high coverage of the input space [12].

2.5. The Extended FAST Method

The extended FAST (Fourier amplitude sensitivity test) method [35], a variance-based SA method, is an extension of the classical FAST method [36]. The main features of the classical FAST method are that the input parameter space is sampled using a periodic transformation function [14] and the output variance is separated into the variances associated with the different inputs. Examples of different sampling functions are given in [14]. The transformation function proposed by [35] was used in this paper.

$$x_i=\left(0.5+\frac{1}{\pi }arcsin[sin({\omega }_{i}s+{\phi }_i)]\right)·K_i$$

(5)

where $x_i$ is the i^th parameter, s is the sampling range with $s\in [-\pi ,\pi ]$, ${\phi }_i$ is a random phase shift parameter marking the starting point of the search curve where ${\phi }_i\in [0,2\pi ]$, ${\omega }_i$ ($i=1,\dots ,k$) is the individual assigned integer frequency of parameter $x_i$, and $K_i$ is a scaling factor to scale the value of the transformation function that lies between zero and one to the appropriate parameter range. ${\omega }_i$ ($i=1,\dots ,k$) is the individual frequency of the $i^{\mathrm{th}}$ parameter, and for each parameter, a different integer frequency $\omega$ must be selected, in a way that no single frequency is a multiple of another. Since the transformation function is periodic, a sampling range of $2\pi$ is sufficient for decomposition. Based on the transformation equation and the different frequencies $\omega$ for the input parameters, a set of N samples is generated for the model simulation. The model output can be decomposed into a Fourier series, whereas the variance of the model output is assigned to the inputs. The first-order sensitivity indices ($S_i$) describe the main effects of the parameters by quantifying how the variance of each input contributes to the total output variance. The total sensitivity index $S_{T_i}$ for an input includes the variance of that input, but it also accounts for the variance created by its interactions with other parameters. The difference between the total and first-order sensitivity indices represents the effect of interactions with other parameters. A detailed description of the extended FAST method applied here can be found in [12].

We chose different sets of independent frequencies based on an algorithm proposed by [14]. We set the maximum frequency for the parameter of interest to ${\omega }_{max}=4096$; hence, the maximum frequency for the other parameters was 512. We considered a maximum harmonic of $M=4$, which resulted in a minimum sampling size for each individual output of $N=2M·{\omega }_{max}+1$ = 32,769 [35].

Following our approach in the extended parameter screening, we calculated the variance-based sensitivity measures at each simulation time step t. In this study, we focus on the direct influence of parameters on the output, which is expressed by the first-order sensitivity index. For each input parameter, we calculated the daily main effects $S_i(t)$ for $t\in \{1,\dots ,d_{max}\}$ to analyze the change in the parameter sensitivity throughout the simulation time.

2.6. Top-Down Concordance Coefficients

We used top-down concordance coefficients (TDCCs, [37,38,39]) based on Savage scores [40] to investigate the temporal dynamics of the parameter sensitivities. For each simulation day, we calculated a parameter ranking based on the main effect of parameters. Then, we calculated TDCCs to compare the parameter rankings for different simulation days. TDCCs are characterized by an emphasis on the agreement among important model parameters, whereas disagreements between less important ones are de-emphasized. A high TDCC close to 1.0 indicates that influential parameters are ranked similar, whereas a low TDCC close to zero means that the ranking order of the influential parameters differs considerably.

2.7. Implementation of the SA

The SA of the MONICA model was implemented following the approach of [12]. The SA algorithms, i.e., the screening method and the extended FAST method, were coded in Python [41]. Python/SWIG [42] was used to create a new software interface for the MONICA model that would allow direct access to the model’s data structures, variables, and functions [12]. The MONICA parameter values were thus directly accessed and modified within the Python scripts. The SA of the MONICA model was executed in parallel on a high-performance computer (HPC) using MPI for Python [43,44]. Thirty-two cores were used to run the independent SA calculations, resulting in a speed-up factor of 27.

3. Results

3.1. Parameter Screening

For each study site, we screened the MONICA model for the most relevant parameters based on the original and extended Morris approaches. We compared the numbers of relevant parameters identified by both approaches for all study sites. We analyzed how the number of relevant parameters changed for both approaches by varying the threshold values ${\tau }_o$ and ${\tau }_t$ between 0.2 and 0.9. As an example, Figure 2 depicts the evolution of the number of parameters for $AGB$. Regardless of the threshold value used, the number of relevant model parameters identified by the extended Morris approach was significantly higher than that for the original approach.

Applying both approaches, our objective was to identify all model parameters that were highly influential at some during simulation time. Choosing ${\tau }_o=0.2$ for the original approach, we could identify all model parameters that affected the output at some point during simulation time. With a ${\tau }_t=0.8$ for the extended approach, only parameters were selected that were highly influential at some point during simulation time. When analyzing ${\mathsf{\Psi }}_o$ and ${\mathsf{\Psi }}_t$, we found that the numbers of identified parameters were comparable, and the majority of the identified model parameters were identified by both approaches (Figure 3). However, some parameters were identified by the new approach (${\mathsf{\Psi }}_t$) that were not identified by the original approach, whereas some parameters were also found to be relevant using the original approach that were not included in ${\mathsf{\Psi }}_t$.

Table 3. Relevant model parameters as identified using the time-dependent Morris approach with ${\tau }_t=0.8$. A parameter was marked as relevant if its ${\mu }^{*}(t)$ was higher than the threshold value ${\tau }_t$ in one of the investigated sites.

Parameter	Description	Unit	Nominal	grainYield	AGB	N${}_{\mathit{AGB}}$	ET${}_{\mathit{a}}$	Moist	N${}_{\mathit{min}}$
AOMSlowUtilEff	Added organic matter slow utilization efficiency		0.4		x	x			x
assimPartitioningStage1${}_{leaf}$	Portion of assimilates assigned for leaf growth in Stage 1	%	0.5		x	x	x	x	x
assimPartitioningStage2${}_{leaf}$	Portion of assimilates assigned for leaf growth in Stage 2	%	0.2				x
assimPartitioningStage2${}_{shoot}$	Portion of assimilates assigned for shoot growth in Stage 2	%	0.6				x
baseTemp${}_2$	Base temperature for assimilation in Stage 2	${}^{\circ }$C	1	x			x	x
baseTemp${}_3$	Base temperature for assimilation in stage 3	${}^{\circ }$C	1	x
baseTemp${}_5$	Base temperature for assimilation in Stage 5	${}^{\circ }$C	9				x	x	x
beginSensPhaseHeatStress	Temperature sum marking the start of the sensitive phase for heat stress	${}^{\circ }$C	620	x	x				x
CNRatioSMB	C-to-N ratio of the soil microbial biomass	−	6.7		x	x			x
cropHeightP1	Factor for crop height	−	6				x	x
cropHeightP2	Reduction factor for crop height	−	0.5				x	x
daylengthReq${}_2$	Day length required for maximum growth in Stage 2	h	20	x		x		x	x
daylengthReq${}_3$	Day length required for maximum growth in stage 3	h	20	x			x
daylengthReq${}_4$	Day length required for maximum growth in Stage 4	h	20				x
denit3	Denitrification coefficient		0.9						x
endSensPhaseHeatStress	Temperature sum marking the end of the sensitive phase for heat stress	${}^{\circ }$C	740	x	x				x
initRootingDepth	Initial root depth of the crop	m	0.1			x
LT50Cultivar	Threshold temperature below which 50% of the crop dies from frost injury	${}^{\circ }$C	−24	x	x	x	x	x	x
luxuryNCoeff	Coefficient describing the maximum N concentration relative to the critical N concentration in the crop tissue	−	1.3			x			x
maintRespP1	Q${}_{10}$ factor for maintenance respiration	−	0.08		x				x
maxAssimRate	Maximum assimilation rate per leaf area	kg CO${}_2$ ha${}^{-1}$	52		x	x	x		x
maxCropHeight	Maximum crop height	m	0.83				x	x
maxCropNDemand	Maximum amount of soil mineral N to be taken up by the crop	$\mathrm{kg}{\mathrm{m}}^{-2}$	6			x			x
minAvailableN	Mineral N concentration in the soil that is not available for crop N uptake	$\mathrm{kg}{\mathrm{m}}^{-3}$	0.0008		x	x			x
minTempAssim	Minimum temperature required for assimilation	${}^{\circ }$C	4		x	x	x	x	x
NConcAGB	Default N concentration in above-ground biomass	kg kg${}^{-1}$	0.06		x	x		x
nConcPN	Shape factor for the critical N curve	−	1.6			x			x
nitrRateCoeffStand	Nitrification rate default coefficient	d${}^{-1}$	0.2						x
organMaintResp${}_{shoot}$	Maintenance respiration factor for shoots	kg CO${}_2$ kg DM${}^{-1}$	0.15	x	x	x		x	x
partSMBSlowToSOMFast	Portion of the soil microbial biomass that is added to the fast soil organic matter pool		0.6						x
partSOMToSMBSlow	Portion of the soil organic matter that is added to the slow soil microbial biomass pool		0.015						x
referenceAlbedo	FAO reference albedo for green grass	−	0.23			x	x	x	x
referenceLAI	Leaf area index of the reference crop	m${}^2$kg${}^{-1}$	1.44				x	x
referenceMaxAssimRate	Maximum assimilation rate of the reference crop	kg CO${}_{2}4$ ha${}^{-1}$	30					x
rootFormFactor	Factor describing the root mass distribution pattern with respect to depth	−	3					x
rootPenRate	Vertical root growth rate	m ${}^{\circ }$C${}^{-1}$ d${}^{-1}$	0.0011		x	x
SMBSlowMaintRateStand	Maintenance rate for slowly-reproducing soil microbial biomass	d${}^{-1}$	0.001						x
SOMSlowUtilEff	Microbial utilization efficiency for the slowly-decomposing soil organic matter pool		0.4						x
specificLeafArea${}_1$	Specific leaf area for calculating the leaf area index for stage 1	m${}^2$ kg${}^{-1}$	0.002		x	x	x	x	x
specificLeafArea${}_2$	Specific leaf area for calculating the leaf area index for stage 2	m${}^2$ kg${}^{-1}$	0.0018				x
stageAtMaxHeight	Stage of maximal crop height	−	3					x
stageKcFactor${}_1$	$K_c$ factor for Stage 1	−	0.4				x	x
stageKcFactor${}_2$	$K_c$ factor for Stage 2	−	0.7	x		x	x	x	x
stageKcFactor${}_3$	$K_c$ factor for Stage 3	−	1.1				x	x
stageKcFactor${}_4$	$K_c$ factor for stage 4	−	1.1				x	x
stageKcFactor${}_5$	$K_c$ factor for Stage 5	−	0.8				x	x
stageTempSum${}_1$	Temperature sum for Stage 1	${}^{\circ }$C	148	x	x	x	x	x	x
stageTempSum${}_2$	Temperature sum for Stage 2	${}^{\circ }$C	284	x		x	x		x
stageTempSum${}_3$	Temperature sum for Stage 3	${}^{\circ }$C	380	x			x	x
stageTempSum${}_4$	Temperature sum for Stage 4	${}^{\circ }$C	180				x
stageTempSum${}_5$	Temperature sum for Stage 5	${}^{\circ }$C	420				x	x	x
stomataCondAlpha	Stomata conductivity parameter	−	40				x	x
vernReq${}_2$	Temperature sum required for optimum vernalization in Stage 2	${}^{\circ }$C	50	x	x	x	x	x	x
Number of relevant parameters				13	16	21	29	27	29

lists the model parameters identified with ${\tau }_t=0.8$ as calculated using the time-dependent Morris approach for all analyzed model outputs. Including the results of different study sites, a parameter was marked as relevant if at any site ${\mu }^{*}$ was higher then the threshold value ${\tau }_t$. It shows that each model output is affected by a different set of model parameters. By comparing the results obtained using the extended approach for each study site, we found that each site was affected by different model parameters, whereas the identified parameters differed both in their number and composition. For example, the number (n) of relevant parameters for AGB differed. For Ascha and Guelzow, $n=8$ parameters were selected, for Dornburg $n=10$, and for Ettlingen and Werlte, $n=12$. The parameter $maxAssimRate$ was sensitive for all sites, but $minAvailableN$ affected AGB only in Ascha, Ettlingen, or Werlte. However, because the screening method was only used as a preliminary means of identifying the most relevant model parameters for each output, we did not further analyze the differences among the sites using these screening results; instead, as described in Section 3.2.4, this comparison was performed using the more reliable sensitivity measures calculated using the time-dependent extended FAST method.

3.2. Extended FAST Method

3.2.1. Original Approach

After identifying the most relevant model parameters based on the results of the extended screening method, we calculated the first-order ($S_i$) and total ($S_{T_i}$) sensitivity indices using the extended FAST approach. The computation of both first-order and total sensitivity indices allows one to distinguish between main and interaction effects. The resulting estimates for the $S_{T_i}$ averaged over all sites are shown in Figure 4.

The parameter sensitivity differed slightly among the analyzed sites. For Dornburg, the main variations in grainYield and AGB were caused by different sensitivity indices. For the heat stress parameters, the $S_{T_i}$ and $S_i$ values were slightly lower compared with those at the other locations, where both sensitivity indices were quite similar. For Dornburg, the sensitivity indices for daylengthReq${}_3$, LT50Cultivar, stageTempSum${}_2$, stageTempSum${}_3$, and vernReq${}_2$ were higher compared with the other sites. Similar results were found for N${}_{AGB}$ and N${}_{min}$, for which parameters related to the crop phenology were associated with considerably higher sensitivities in Dornburg compared with the other locations. Analyzing the parameter sensitivity indices for ET${}_a$ at the different locations revealed that although the orders of the parameters based on the sensitivity measures were the same, Ettlingen had the highest $S_{T_i}$ and $S_i$ values, followed by Gülzow and Werlte. The parameter sensitivities for Moist also differed among the analyzed sites, but no generally applicable order of the important parameters based on $S_{T_i}$ and $S_i$ was identified in this case.

3.2.2. Time-Dependent Main Effects

In addition to applying the original approach, we calculated the main effects of the model parameters for each time step of the simulation. Figure 5 displays the results of the temporal sensitivity analysis aggregated over all sites for a selection of the relevant model parameters. A complete overview of the calculated daily main effects for each analyzed parameter differentiated among the different sites is provided in the Appendix A in Figure A1, Figure A2, Figure A3, Figure A4, Figure A5 and Figure A6.

For biomass-related outputs such as grainYield, AGB, and $N_{AGB}$, the time-dependent extended FAST method identified more sensitive parameters compared with the original approach. The temporal analysis of the main effects enabled the identification of not only parameters that were relevant at harvest time, but also ones that were relevant during earlier stages of crop development. assimPartitioningStage1${}_{leaf}$ and specLeafArea${}_1$ strongly affected AGB at the beginning of crop growth, but their influence decreased strongly over time until harvest. For N${}_{AGB}$, the influence of NConcAGB was very high during the first stage of development, but decreased once the flowering stage began.

Similar results were found for the soil-related outputs ET${}_a$, Moist and N${}_{min}$, for which the temporal sensitivity analysis enabled the detection of parameters that were influential at various times throughout the simulation, in contrast to the original approach, which analyzed the parameter sensitivities at only a single point in time. For ET${}_a$ and Moist, the sensitivities to the crop’s transpiration parameters stageKcFactor${}_1$–stageKcFactor${}_4$ were very high in their respective developmental stages, but decreased as harvest time approached.

3.2.3. Temporal Sensitivity Patterns

For each simulation day, we generated a parameter ranking based on the daily main effects. We used the TDCCs to compare the parameter rankings for different simulation days averaged over all sites (Figure 6). For all crop-specific outputs (grainYield, AGB, and N${}_{AGB}$), we compared the parameter sensitivity rankings solely for the cultivation period, because these outputs were not calculated for the fallow time. grainYield was calculated only for the last three developmental stages of the winter wheat. For this short period, we found an annual pattern of similar parameter sensitivities in which the sensitivities were more comparable at the beginning of this period and then decreased. For AGB, regardless of the simulation year, the agreement among parameter rankings was higher in earlier developmental stages at all sites. With the onset of the flowering phase, the differences in the parameter rankings increased. This general pattern could be found in every study site. However, analyzing the TDCC of AGB for each study site (Figure 7) showed that there are some periods where this pattern was only weakly pronounced. In Ettlingen, for example, the parameter ranking of 2006 and 2007 during spring differed highly. This indicates that due to differences in the site-specific climate conditions of both years, the parameter sensitivity for AGB differed. 2007 was characterized by a long spring drought with no precipitation compared to a relative wet spring in 2006. For N${}_{AGB}$, we found few agreements between parameter rankings for the early developmental stages in a single simulation year. For a short period at the end of crop growth, the parameter rankings for the different simulation years were quite similar. Analyzing the TDCC of N${}_{AGB}$ showed less concordance for 2007 compared to the other years; again, presumably due to the exceptional weather conditions in 2007.

Fewer seasonal patterns were found for the soil-specific model output. We identified a small annual pattern of similar parameter sensitivities for ET${}_a$ and Moist, especially at the beginning of crop growth, but also during further crop development. This means that sensitivity ranking did not differ for these periods. Comparing these time periods with the temporal main effects in Figure 5j–l, we found that $E{T}_a$ and $Moist$ were mostly affected by the parameters $stageKcFacto{r}_1$–$stageKcFacto{r}_4$ regardless of the analyzed study site. For both outputs, the TDCC showed high consistency in parameter ranking in the early simulation days ending with the sowing of the crop. However, the TDCC also showed that the parameter sensitivity of the first growing season beginning with sowing in October 2005 differed strongly compared to the following growing seasons. Only a few agreements between the parameter sensitivities for different simulation days were found for N${}_{min}$. Regardless of the simulation year, a high TDCC was found for a short period directly before crop harvest. In this period, N${}_{min}$ was mainly affected by parameters $luxuryNCoeff$ and $NConcPN$. Additionally, the TDCC for N${}_{min}$ showed that the parameter rankings during crop growth for the last two cultivation periods were more similar than those for the first cultivation period.

3.2.4. Influence of the Different Sites

Regardless of the selected output, the parameter sensitivities differed among the analyzed sites (see Figure 8 and Figure A1, Figure A2, Figure A3, Figure A4, Figure A5 and Figure A6 in the Appendix A). Only minor differences in parameter sensitivity were found for grainYield. Site-specific differences for AGB were caused mainly by the parameters assimPartitioningStage1${}_{leaf}$, maxAssimRate, specLeafArea${}_1$, and minTempAssim. For N${}_{AGB}$, the variations were caused by NConcAGB and minAvailableN, especially during the first cultivation period. Regardless of the cultivation period, the parameter sensitivities of the remaining relevant parameters consistently differed among the different sites. Substantial variations in the sensitivities to several parameters were found for ET${}_a$, Moist, and N${}_{min}$, originating from site-specific differences in the most influential parameters.

4. Discussion

The previously-performed SA of the MONICA model [12] increased the understanding of the model and served as input for a later model calibration, but it yielded only a static snapshot of the parameter sensitivity at harvest time. The time-dependent SA presented in this study provides a deeper understanding of the model behavior and more detailed information about the influence of the model parameters throughout the simulation time. Many studies have used snapshot SAs to obtain general information on model behavior [12,15,16,33,39,45,46,47,48]; however, only a few studies have addressed the fact that parameter sensitivity is a function of time, and changes throughout the simulation time [18,19,33,49,50,51]. Ref. [20,21] present an alternative approach to account for the dynamics for parameter sensitivity complementing the approach applied in this study. This study, however, is the first to analyze the effect of the simulation time on the parameter sensitivity based on elementary effects calculated using the Morris screening design. The general approach used in this study can be easily applied for time-dependent SAs of other ecological models.

The results of the original and adopted Morris approaches strongly depend on the chosen threshold values ${\tau }_o$ and ${\tau }_t$ that are used to identify the relevant model parameters. When the threshold value is lower, a greater number of relevant model parameters are identified. We compared the identified parameters in ${\mathsf{\Psi }}_o$ with the parameters in ${\mathsf{\Psi }}_t$. We used a low ${\tau }_o$ to exclude parameters with negligible influence from ${\mathsf{\Psi }}_o$, respectively identify influential parameters. For ${\mathsf{\Psi }}_t$, we chose a high threshold value ${\tau }_t=0.8$ to identify model parameters that had a significant influence at any point during the simulation time. With all outputs considered, the numbers of parameters in both sets were comparable, but not identical. The majority of parameters were identified in both sets, but some parameters were found to be relevant using one approach, but negligible according to the other approach. Similarly to [12,18,19,45,52], we used the Morris approach to screen out less influential parameters. The selection of a high threshold value allowed us to identify model parameters that had a strong effect at some point during the simulation time. Because the results of the original approach represent the aggregated parameter sensitivity based on a single point in time, we found the set of parameters identified using the temporal extension for the objective of this study to be more reliable and used the parameters in ${\mathsf{\Psi }}_t$ for further temporal sensitivity analysis.

The original extended FAST method provides a measure of the sensitivity at only a single point throughout the simulation time, in our case at the harvesting of the winter wheat. The daily analysis of the main effects performed in this study showed that more parameters than those identified by the original approach affected the selected model outputs. With the new approach, we were also able to analyze the relationships between parameter sensitivity and crop growth using daily sensitivity values. Several parameters were highly influential during the early stages of crop development, but the corresponding parameter sensitivities declined during crop growth. For example, stageTempSum${}_1$ had a minor effect on AGB (S${}_i$ = 0.0) according to the original approach. However, upon analyzing the parameter sensitivity using the daily main effect values, we found that this parameter had a major impact on AGB during the initial stages of crop development (S${}_i>0.8$) and then steadily decreased in influence until harvest time. Similar observations apply to other parameters (specLeafArea${}_1$ for AGB, NConcAGB for N${}_{AGB}$, stageKcFactor${}_1$ for Moist), whose associated sensitivities decreased during crop development or with increasing simulation time. Additionally, we identified parameters that had no effect at the beginning of crop growth, but whose influence increased during crop development, e.g., endSensPhaseHeatStress for grainYield and NConcPN for N${}_{AGB}$ or N${}_{min}$. Other parameters showed a peak-shaped sensitivity evolution; the model’s sensitivity to these parameters was low during early stages of crop growth, then increased to a maximum, and finally steadily decreased once again as the simulation progressed further (stageKcFactor${}_1$–stageKcFactor${}_4$ for ET${}_a$ or Moist). referenceAlbedo was the only parameter that had a major effect during fallow time, but only a minor influence on ET${}_a$ and Moist during crop growth.

Often, SA at a single simulation time is used to identify relevant model parameters that are later used for model calibration [12,53]. The time-dependent SA results presented here demonstrate that such snapshot SAs may lead to incorrect interpretations of the parameters’ influences. Calibration studies relying on snapshot SAs may be insufficient and lack the necessary information on how the parameter sensitivity changes over the simulation time. Additionally, the analyses of parameter sensitivities for different sites indicate the strong influence of soil and climate conditions on parameter sensitivity [18,19]. A temporal sensitivity analysis can improve the calibration process for a model, depending on the scope of the calibration and the data available to be used for calibration. Temporal sensitivity analysis enables the identification of influential parameters for different simulation periods. A crop model calibration can benefit from such results if measured data regarding the model outputs are available for different times throughout the model simulation. Then, a set of parameters that are specific to a certain simulation time for which measured data exist can be selectively adopted to improve the results of the model calibration. If experimental data are available for only one specific time, then a snapshot sensitivity analysis should be sufficient for model calibration.

The TDCCs based on the parameter rankings for different simulation days revealed some temporal patterns of similarity in parameter sensitivities. For each output, we found an annual pattern of similar parameter sensitivities for some phase of crop development, although the climate conditions differed for each cultivation year. For the crop-specific model outputs, more pronounced seasonal sensitivity patterns could be identified, which means that the calculated parameter rankings were similar for longer simulation periods, as well as for different simulation years. Such a pattern could be expected as some relevant model parameters are only active during a specific developmental stage. Parameter rankings for soil-related model outputs such as Moist or N${}_{min}$ were less homogeneous and were more affected by site-specific conditions of different simulation years. Similar parameter rankings were only identified for small simulation periods during crop development caused by the high influence of only a few parameters for this period. However, the temporal TDCCs also demonstrated the site-specific effect of soil and climate conditions on the parameter sensitivity, as the parameter rankings from one year differed highly compared to the other simulation years.

Differences in the soil and climate conditions also affected the parameter sensitivities [13,15,18], as indicated by the fact that the values and shapes of the sensitivity results differed for the different analyzed study sites. ET${}_a$ and Moist showed the highest differences in parameter sensitivity. Because both of these outputs strongly rely on the water storage capacity of the soil and the available water (precipitation and capillary rise from groundwater), both are highly sensitive to different site and climate conditions. N${}_{AGB}$ and N${}_{min}$ were affected by a similar set of parameters (e.g., minAvailableN, luxuryNCoeff, NConcPN), although the form of the relationship differed slightly. The parameter sensitivities for the different analyzed sites differed greatly, especially with regard to the most influential parameters. Both outputs depend on the available nitrogen (N) stored in the soil, which is influenced by the soil type at the site. Additionally, mineralization processes of the soil microbial biomass (SMB) that transform added organic matter (e.g., crop residues or organic fertilizer) into plant-available nitrogen are highly dependent on temperature. Both of these factors contributed to the site-specific differences we detected in the parameters to which N${}_{AGB}$ and N${}_{min}$ were sensitive. Because the simulation time over which grainYield was calculated was rather short, only minor site-specific differences were identified in this output.

The temporal sensitivity analysis in this study was based on the calculation of first-order sensitivity indices. We used main effects (1) for analyzing the course of sensitivity during simulation time and (2) for generating the parameters’ ranking on which the calculation of TDCCs was based. As Figure 4 shows, the interaction effects of parameters are quite high. In this study, we deliberately focused on the main effects to quantify the direct effect of model parameters on the output. Using the total effect would not deliver any information on the direct contribution of the model parameters or on the contribution of the interaction with other parameters. For a sound interpretation of total effects, it is still necessary to know the main effect of model parameters.

In this study, we deliberately focused on the main effects to quantify the direct effect of model parameters towards the output. By using the total effect, there would be a high uncertainty for how much this influence was caused directly by the model parameters or how much it was caused by the interaction with other parameters. For a correct interpretation of the total effects, it is still necessary to know the main effect of the model parameters.

In this study, SA was performed at five different sites to demonstrate the influence of different soil and climate conditions on parameter sensitivity. We focused on the general effects of the different sites on the time-dependent SA results. Consequently, this study lacks a detailed analysis of how the parameter sensitivity is affected by single driving variables (e.g., soil texture, air temperature, precipitation). Future SA studies performed to further improve the understanding of the model behavior should consider both the temporal dynamics of parameter sensitivity and the effects of different sites and climate conditions. Such a study using a more comprehensive SA scheme should therefore concentrate on the specific effects of individual driving variables on the time-dependent parameter sensitivity.

5. Conclusions

The temporal dynamics of parameter sensitivity and the use of site-specific data will typically hamper a comprehensive sensitivity analysis of an (agro-)ecological model. Our comparison of the original and time-dependent screening approaches revealed that the original approach, which aggregates the parameter sensitivity over a certain simulation time, will fail to identify parameters to which the model is highly sensitive for only a limited period during the simulation. The time-dependent screening approach permits more reliable identification of the relevant model parameters, which helps to yield a better understanding of the behavior of process-based models and to better assist in their improvement. It is evident that calibration studies in which SA is applied to identify relevant model parameters must include an analysis of the temporal dynamics of parameter sensitivity; otherwise, the results of a snapshot SA may lead to incorrect interpretations of the parameters’ influences.

Most process-based agro-ecosystem models include processes that are only active during a certain period of the simulation, and consequently, their parameters are similarly only relevant during a certain period. A model calibration can benefit from a temporal SA because it enables a detailed identification of the influential model parameters that acknowledges their limited temporal activity. Unlike a snapshot SA, a temporal SA also permits the analysis of the relationships between model parameters and outputs, which supports a more thorough understanding and evaluation of model behavior. This is particularly important for understanding the influence of site conditions on the simulation and considering parameter estimation with respect to a particular environment. A temporal SA may provide hints regarding how to prioritize the experimental observations of crucially sensitive parameters to improve the general logic of the experiment-algorithm-simulation model-prediction chain.

Although it still seems challenging to quantify the influence of site conditions to obtain a more objective interpretation of parameter sensitivity, this study already serves as a step toward reducing the uncertainty of model simulations conducted with limited available data, as demonstrated by [54] for the improvement of a model algorithm that had previously been identified in a variable-focused SA as being a source of disagreement between models [55].