Sensitivity and Uncertainty Analyses of Flux-based Ecosystem Model towards Improvement of Forest GPP Simulation

Abstract

An ecosystem model serves as an important tool to understand the carbon cycle in the forest ecosystem. However, the sensitivities of parameters and uncertainties of the model outputs are not clearly understood. Parameter sensitivity analysis (SA) and uncertainty analysis (UA) play a crucial role in the improvement of forest gross primary productivity GPP simulation. This study presents a global SA based on an extended Fourier amplitude sensitivity test (EFAST) method to quantify the sensitivities of 16 parameters in the Flux-based ecosystem model (FBEM). To systematically evaluate the parameters’ sensitivities, various parameter ranges, different model outputs, temporal variations of parameters sensitivity index (SI) were comprehensively explored via three experiments. Based on the numerical experiments of SA, the UA experiments were designed and performed for parameter estimation based on a Markov chain Monte Carlo (MCMC) method. The ratio of internal CO₂ to air CO₂ ($f_C{}_i)$, canopy quantum efficiency of photon conversion $({\alpha }_q)$, maximum carboxylation rate at 25 $°\mathrm{C}$ ($V_m^{25})$ were the most sensitive parameters for the GPP. It was also indicated that ${\mathsf{\alpha }}_{\mathrm{q}},E_{\mathrm{V}}{}_m^{}\mathrm{and}{Q}_{10}$ were influenced by temperature throughout the entire growth stage. The result of parameter estimation of only using four sensitive parameters (RMSE = 1.657) is very close to that using all the parameters (RMSE = 1.496). The results of SA suggest that sensitive parameters, such as $f_{ci},{\mathsf{\alpha }}_{\mathrm{q}},{\mathrm{E}}_{\mathrm{V}}{}_m^{},V_m^{25}$ strongly influence on the forest GPP simulation, and the temporal characteristics of the parameters’ SI on GPP and NEE were changed in different growth. The sensitive parameters were a major source of uncertainty and parameter estimation based on the parameter SA could lead to desirable results without introducing too great uncertainties.

1. Introduction

Ecosystem models are valuable tools that describe and explain the processes and variable dynamics of photosynthesis and respiration in a forest ecosystem [1,2]. They also play a key role in assessing the carbon equilibrium and the response of a forest ecosystem to global change at local and global scales [3,4]. Existing ecosystem models integrate many ecosystem processes and are coupled into the earth system model (ESM) [5,6,7]. Terrestrial gross primary productivity (GPP) was the largest component flux of the global carbon cycles [8,9]. There had uncertainty in the interannual variability of GPP simulation of the majority of ecosystem models [10]. The more complicated the models are, the more parameters they may have. This may introduce larger uncertainty caused by the parameter sensitivity, thus, it is indispensable to quantify the uncertainty and sensitivity of the parameters in forest ecosystem models. Under such background, various sensitivity analyses (SAs) were proposed and applied to ascertain the corresponding responses in the output variables when input parameters alter within their valid ranges [11,12]. Indeed, various SA algorithms have been widely expanded to ecological models [13,14], SCOPE model [15], hydrological models [16,17], remote sensing models [18] and crop growth models [12]. Thus, performing an SA is a feasible way to characterize and reduce uncertainties in ecosystem models, and to improve their performance [19,20].

The global and local SA algorithms can identify the influence of parameters on variations of the outputs in the model and are therefore key to understand the model performance in response to variations in environmental factors [12,21]. Especially, global SA could be applied to non-linear models, but local SA could not [22,23]. The global SA considers non-linear responses and parameter interactions [24], providing a comprehensive identification of parameter sensitivity on model outputs. For example, the Sobol’ [25], Fourier amplitude sensitivity test (FAST), and extended FAST (EFAST) are variance-based SA methods, and they can quantify the sensitivities of model parameters, including first-order sensitivity index (or main sensitivity index, MSI), total sensitivity index (TSI) and interactions [26].

SA helps parameter uncertainty analysis (UA) and can identify the major uncertainties when combined with parameter optimization. Note that parameter optimization is an essential means of calibrating ecosystem models and increasing the accuracy of their predictions [27,28]. However, the parameter calibration of complex process models is time-consuming and inaccurate owing to a lack of results in SA in general [29,30,31,32,33,34,35].

The flux-based ecosystem model (FBEM) was also widely used for parameter estimation or data assimilation studies [28,36,37,38]. However, the model has as many as 16 parameters, and how these parameters govern the model behavior is still not clear. Thus, it is necessary to identify the sensitivities of the parameters, informing understanding of the uncertainties in the model. Specifically, we need to estimate the parameter sensitivities in FBEM and their variations in the growing season, as well as the relationship between parameters SI and environmental factors. In this study, the FBEM and EFAST algorithm were combined for parameters SA experiment. Three experiments were designed to systematically evaluate the parameters’ sensitivities. Various parameter ranges, different model outputs and temporal variations of parameters SIs were carefully considered.

The objective of this study is to evaluate the sensitivity index (SI) of all the parameters on the GPP of FBEM using the EFAST method and to improve the parameter estimation process based on the results of the SA for quantizing uncertainty. Four scientific questions are addressed: (1) how the parameter ranges influence their sensitivities on forest GPP simulation? (2) which parameter(s) predominates the model behavior and leads to uncertainty of the model, i.e., which parameters were most sensitive for GPP and net ecosystem carbon exchange (NEE) simulation? (3) how the parameter SI variated during the growing season, and how the environmental factors influence the variation of SI? and (4) did the uncertainties caused by parameters mainly come from sensitive parameters? To answer these questions, we designed numeric experiments: (1) to analyze the influence of parameter ranges, models outputs on parameters SIs, (2) to identify the most influential parameters on the output in the FBEM model, (3) to further analyze the relationship between sensitivity index, uncertainty and environmental factors based on the results of the SA.

2. Materials and Methods

2.1. Flux-based Ecosystem Model and Parameters

In this study, FBEM is used to conduct ecosystem modeling. FBEM is originally designed to compute net CO₂ ecosystem exchange with a parameter estimation module [38,39]. FBEM covers two major carbon cycling processes: canopy-level photosynthesis (A_c) and ecosystem respiration (ER). In the model, the two processes are mainly regulated by four environmental variables: leaf area index (LAI), air temperature (T_a), relative humidity (RH) and photosynthetically active radiation (PAR). Leaf-level photosynthesis (A) is calculated by electron transport rates (J_e) of light and the rates of carboxylation enzymes (J_c), and it is finally scaled up to A_c in FBEM. Details of each part of the model and its equations are listed in

Table 1. Symbols, definitions, units, initial values, ranges and sources of 16 parameters of the Flux-based ecosystem model (FBEM) used in a sensitivity analysis.

Parameter	Definition	Unit	Value	Range		Reference
				Minimum	Maximum
${\mathit{\alpha }}_{\mathit{q}}$	Canopy quantum efficiency of photon conversion	mol mol−1 photon	0.28	0	0.5	[40]
${\mathit{K}}_{\mathit{c}}^{25}$	Michaelis–Menten constant for carboxylation	$\mathsf{\mu }$mol mol−1	460	50	600	[40]
${\mathit{E}}_{\mathit{K}}{}_{\mathit{c}}^{}$	Activation energy of ${\mathrm{K}}_{\mathrm{c}}^{25}$	Jmol−1	59,356	30,000	150,000	[40]
${{\mathit{E}}_{\mathit{K}}}_{\mathit{o}}$	Activation energy of ${\mathrm{K}}_{\mathrm{o}}^{25}$	Jmol−1	35,948	10,000	60,000	[40]
${\mathit{K}}_{\mathit{o}}^{25}$	Michaelis–Menten constant for oxygenation	mol mol−1	0.33	0.2	0.5	[40]
${\mathit{E}}_{\mathit{V}}{}_{\mathit{m}}^{}$	Activation energy of ${\mathrm{V}}_{\mathrm{m}}^{25}$	Jmol−1	58,520	10,000	100,000	[40]
${\mathit{\Gamma }}_{\ast }^{25}$	CO2 compensation point without dark respiration	μmol mol−1	42.5	10	200	[40]
${\mathit{r}}_{{\mathit{J}}_{\mathit{m}{\mathit{V}}_{\mathit{m}}}}$	Ratio of ${\mathrm{J}}_{\mathrm{m}}$ to ${\mathrm{V}}_{\mathrm{m}}^{25}$ at 25 ℃	-	1.79	1	5	[40]
${\mathit{ER}}_{\mathbf{0}}$	Whole ecosystem respiration at 0 °C	μmol CO2 m−2 s−1	2.5	1	5	[41]
${\mathit{Q}}_{\mathbf{10}}$	Temperature dependency of ecosystem respiration	-	2	1	3	[41]
${\mathit{V}}_{\mathit{m}}^{\mathbf{25}}$	Maximum carboxylation rate at 25 °C	μmol CO2 m−2 s−1	29	10	300	[40]
${\mathit{f}}_{\mathit{C}}{}_{\mathit{i}}$	Ratio of internal CO2 to air CO2	-	0.87	0.5	0.9	[40]
${\mathit{K}}_{\mathit{n}}$	Canopy extinction coefficient for light	-	0.8	0.7	0.9	[40]
${\mathit{E}}_{{\mathit{\Gamma }}_{\ast }^{\mathbf{25}}}$	Activation energy of CO2 compensation point at 25 °C	J mol−1	60,000	30,000	100,000	[40]
${\mathit{g}}_{\mathit{l}}$	Empirical coefficient in Leuning model	-	1657	100	2000	[42]
${\mathit{D}}_{\mathbf{0}}$	Empirical coefficient in Leuning model	kPa	2.74	1	10	[42]

.

2.1.1. Leaf-level photosynthesis

Leaf-level photosynthesis (A) is based on the model developed by Farquhar [43]. For C3 plants, gross leaf CO₂ uptake (A, $\mathsf{\mu }$mol CO₂ m⁻² s⁻¹) is calculated as follows:

$$A=\min \left\{J_C,J_e\right\}$$

(1)

where $J_e$ and $J_C$ represent light electron transport ($\mathsf{\mu }$mol CO₂ m⁻² s⁻¹) and the rates limited by carboxylation enzymes ($\mathsf{\mu }$mol CO₂ m⁻² s⁻¹), respectively.

$$J_C=V_m\times \frac{C_i-{\mathsf{{\rm T}}}_{\ast }}{C_i+K_c\times (1+\frac{O_x}{K_o})}$$

(2)

$$C_i=f_{C_i }\times C_a$$

(3)

$${\mathrm{V}}_{\mathrm{m}}{=\mathrm{V}}_{\mathrm{m}}^{25}\times \exp \frac{{\mathrm{E}}_{{\mathrm{v}}_{\mathrm{m}}}{\times (\mathrm{T}}_{\mathrm{k}}-298)}{{\mathrm{R}\times \mathrm{T}}_{\mathrm{k}}\times 298}$$

(4)

where $C_i$ is the leaf’s internal CO₂ concentration ($\mathsf{\mu }$mol CO₂ mol⁻¹ air), $C_a$ is ambient CO₂ concentration, $f_{C_i }$ the ratio of the leaf’s internal CO₂ content to its concentration in the ambient air, O_x is oxygen concentration in the air (0.21 mol O₂ mol⁻¹ air), and $V_m$ is the maximum carboxylation rate ($\mathsf{\mu }$mol CO₂ m⁻² s⁻¹), $E_{v_m}$ is the activation energy, $V_m^{25}$ is the maximum carboxylation rate at 25 $°\mathrm{C}$, and R is the universal gas constant (8.314 J K⁻¹mol⁻¹), respectively. The CO₂ compensation point without dark respiration is represented as ${\Gamma }_{\ast }$($\mathsf{\mu }$mol CO² m⁻¹). It is modified by using Arrhenius’s equation.

$${\Gamma }_{\ast }{=\Gamma }_{\ast }^{25}\times \exp \left(\frac{E_{{\Gamma }_{\ast }^{25}}\times \left(T_k-298\right)}{{R\times T}_k\times 298}\right)$$

(5)

where, ${\Gamma }_{\ast }^{25}$ is the CO₂ compensation point without dark respiration at 25 °C and $E_{{\Gamma }_{\ast }^{25}}$ describes the temperature dependence of ${\Gamma }_{\ast }$. Two Michaelis–Menten constants have a temperature dependence based on Arrhenius’s equation similar to $V_m$. $K_c$, the Michaelis–Menten constant for carboxylation ($\mathsf{\mu }$mol mol⁻¹), is calculated by using the following function:

$$K_c{=K}_c^{25}\times \exp (\frac{E_{K_C}{\times (T}_k-298)}{{R\times T}_k\times 298})$$

(6)

where, $E_{K_C}$ is activation energy and $K_c^{25}$ is the Michaelis–Menten constant for carboxylation at 25 $°\mathrm{C}$. Ko, the Michaelis–Menten constant for oxygenation ($\mathsf{\mu }$mol mol⁻¹), is computed by using the following function:

$$K_O=K_O^{25}\times \exp \left(\frac{E_{K_O} \times \left(T_k-298\right)}{R\times T_k\times 298}\right)$$

(7)

where, $E_{K_O}$ is activation energy and $K_O^{25}$ is the Michaelis–Menten constant for oxygenation at 25 $°\mathrm{C}$. The light electron transport processes (J_e with the unit of $\mathsf{\mu }$mol CO₂ m⁻² s⁻¹) are described by using the following function:

$$J_e=\frac{{\alpha }_q\times I\times J_m}{\sqrt{J_m^2+{\alpha }_q^2\times I^2}}\times \frac{C_i-{\Gamma }_{\ast }}{4\times \left(C_i+2{\Gamma }_{\ast }\right)}$$

(8)

where, I is the absorbed PAR ($\mathsf{\mu }$mol m⁻² s⁻¹), ${\alpha }_q$ is the quantum efficiency of photon capture (mol mol⁻¹ photon), and $J_m$ is the maximum electron transport rate ($\mathsf{\mu }$mol CO² m⁻² s⁻¹), respectively. $J_m$ depends on temperature, and is computed by using the following function:

$$J_m{=r}_{J_{{mV}_m}}{\times V}_m^{25}\times \exp (\frac{E_{V_m}{\times (T}_k-298)}{{R\times T}_k\times 298})$$

(9)

where, $r_{J_{m{V}_m}}$ is the ratio of $J_m$ to $V_m^{25}$ at 25 $°\mathrm{C}$.

2.1.2. Stomatal conductance

Stomatal conductance (G_s) is coupled with leaf-level photosynthesis by the Leuning model [44,45]. The flux in carbon at the top layer of the leaf (A_n) is estimated by using the following equations:

$${\mathrm{A}}_{\mathrm{n}}=G_s\times \left(C_a-C_i\right)$$

(10)

$$G_s{=\mathrm{g}}_{\mathrm{l}}\times \frac{A}{{(C}_a{-\Gamma }_{\ast })\times (1+\frac{D}{D_0})}$$

(11)

where, ${\mathrm{g}}_{\mathrm{l}}$ and D₀ (kPa) are empirical coefficients, and D is a deficit in vapor pressure (kPa) calculated by air temperature (T_k) and RH (in %):

$$\ln e_{\mathrm{s}}=21.382-\frac{5347.5}{T_k}$$

(12)

$$D=0{.1\times e}_{\mathrm{s}}\times 1-RH$$

(13)

where, $e_{\mathrm{s}}$ is saturation vapor pressure (mbar), RH is relative humidity.

2.1.3. Canopy-level Photosynthesis

An approach developed by Sellers [46] is used to scale up leaf-level photosynthesis to canopy-level photosynthesis. It describes the relationship between canopy photosynthesis (A_c) and the carbon influx at the top leaf layer.

$$A_c=A_n\times \frac{1-\exp \left(-k_n\times LAI\right)}{k_n}$$

(14)

where $k_n$ is the light extinction coefficient, and $A_c$ equals to GPP.

2.1.4. Ecosystem Respiration

Ecosystem respiration ($ER$) is modeled as a function of temperature (T_a, $°\mathrm{C}$) by using the widely-used van’t Hoff equation:

$${ER=ER}_0\times Q_{10}^{T_a/10}$$

(15)

where, ${ER}_0$ is ecosystem respiration at 0 °C and Q₁₀ is the relative increase ($ER$/${ER}_0$) in respiration per 10 $°\mathrm{C}$ in temperature. NEE is calculated by using the following function:

$$\mathrm{NEE}=\mathrm{GPP}-ER$$

(16)

In total, 16 parameters dominate model behaviors (Table 1).

2.2. Data

Driving data, used in this experiment, were collected flux tower site (68.740° W, 45.204° N) in Howland forest from 2002~2006, Maine, USA, and the vegetation was mainly evergreen needle leaf forest (ENF), and the climate was warm summer continental. The vegetation of the old evergreen needle leaf forest was red spruce (Pinus rubens Sarg), eastern hemlock (Tsuga canadensis (L.) Carr.) and other conifers. The mean annual temperature (MAT) is 5.27 °C and the mean annual precipitation (MAP) is 1070 mm [37,47,48]. The flux dataset contained daily ecological and environmental data from 2000 to 2006, such as photosynthetically active radiation (PAR), air temperature (T_a), and relative humidity (RH). In addition, the flux and ecological data, such as net ecosystem carbon exchange (NEE), gross primary productivity (GPP), leaf area index (LAI), and ecosystem respiration (ER) were collected from the AmeriFlux (http://ameriflux.lbl.gov/sites). The observations of Howland forest flux site are comprehensive and of high quality, so a large number of ecological models, model-data fusion and remote sensing studies selected the site as a case [28,37,48,49]. This study also selected this flux site for parameter sensitivity and uncertainty research of the forest ecosystem model.

2.3. Extended Fourier Amplitude Sensitivity Test (EFAST)

The EFAST is a well-recognized global SA algorithm that is especially widely applied in the sensitivity analysis for nonlinear models [11], such as the Biome-BGC model [50], crop growth model [12], canopy reflectance model [51], microwave remote sensing models [18]. The algorithm utilizes the main sensitivity index (MSI) and total sensitivity index (TSI) to quantify the sensitivity of the model outputs on various model inputs. The EFAST has also been extensively applied in evaluating the sensitivity of ecosystem models [12,13,52]. Since the algorithm has been well documented and applied, we briefly introduced how to combine the algorithm and FBEM to analyze the sensitivity of the parameters on the model output. In this study, a total of 16 parameters to the FBEM were comprehensively evaluated by computing their MSI and TSI:

$$\mathrm{MSI}=\frac{{\widehat{Var}}_i\left(Y\right)}{\widehat{Var}\left(Y\right)}$$

(17)

$$\mathrm{TSI}=\frac{1-{\widehat{Va\mathrm{r}}}_{-i}\left(Y\right)}{\widehat{Var}\left(Y\right)}$$

(18)

where, $\widehat{Var}\left(Y\right)$ is the estimated conditional variance of the output Y refers to GPP and NEE in this study), and ${\widehat{Var}}_i\left(Y\right)$ is the variance of ith parameter as a factor, ${\widehat{Va\mathrm{r}}}_{-i}\left(Y\right)$ is the variance of all except for the ith parameter as a factor, MSI is supposed to the contribution of each input parameter to the $\widehat{Var}\left(Y\right).$ TSI is the total amount of the MSI and all interaction effects involving the parameter, and additional information is available in [11,23].

2.4. Numerical Experiments for Sensitivity Analysis

2.4.1. Parameters’ SI Variation with Parameter Range

Determining the range and distribution of each parameter was significant in a global SA [53] because they may affect SIs and importance rankings of the parameters. Thus, the purpose of this experiment was to test the impact of parameter ranges on their SIs and importance rankings. The ranges of parameter variation were used to set the upper and lower limits of the 16 parameters of the EBFM, as derived from the references [28,54,55] (Table 1), and GPP was the major model output in this experiment.

One set range was determined by ±10% perturbation and ±30% perturbation of the parameter values described in Table 1, and another was between the minimum and maximum (described in Table 1) to improve the parameter estimation based on the SA. Uniform distribution was in this experiment. Moreover, the foremost objective of the models was the estimation of the GPP, which was therefore considered the model output to assess the influence of the range of parametric variation on the SA.

2.4.2. Comparison of Parameter SIs for GPP and NEE

The above-mentioned experiments were designed to analyze the parameter SIs on GPP. However, for different model outputs, e.g., GPP, NEE and ER, the same parameter may present significantly different SI on them. Therefore, analyzing the roles of the parameters for various state variables was useful. To this end, this numerical experiment was designed to compare the SIs of parameters on different model outputs, which are GPP and NEE.

2.4.3. Temporal Characteristics of Parameter Sensitivity for GPP

The main eco-physiological processes of vegetation vary in each growth stage. For example, in the initial growth stage, vegetation growth was dominant, whereas reproductive growth was dominant after anthesis [2]. The dominant biological processes differ throughout the entire growth period, which results in differences among the dominant parameters. To this end, this experiment was designed to analyze the temporal characteristics of the parameter SIs.

The GPP and NEE were the basis of the model outputs and were thus selected as the target of the analysis of the temporal characteristics of the parameters. In the numerical experiment, minimum and maximum values described in Table 1 were set as the lower and upper boundary of parameter ranges, respectively, and the TSI was used as an evaluation criterion.

2.4.4. Results of Uncertainty Analysis based on Parameter Sensitivity Analysis and Parameter Estimation

We selected the most sensitive parameters from all input parameters through sensitivity analysis. We then determined whether this could enhance our understanding of model uncertainty and help optimize the configuration of the model and whether replacing this model’s parameters with another quickly yields the most appropriate parameter configuration.

Two experiments were carried out to compare the results: (1) Experiment 1: All 16 input parameters were selected and optimized by the Markov chain Monte Carlo (MCMC) method. (2) Experiment 2: Four input parameters with a high SI were used to optimize the model. The results of the experiments were compared to reduce errors in the model and make the optimization more efficient.

The MCMC method can be used to estimate model parameters by incorporating observations into the model. According to Bayes theorem, the posterior probability density functions (PDFs) of the model parameters (p) can be calculated from prior knowledge and information generated by comparing the model with the observed values. The method can be expressed as [36,56]:

$$p\left(\theta |Z\right)=\frac{p\left(\theta |Z\right)p\left(\theta \right)}{p\left(Z\right)}$$

(19)

where $p\left(\theta |Z\right)$ is the posterior distribution of parameters $\theta$ given observations Z, $p\left(\theta \right)$ is a set of uniform distributions over the ranges specified in Table 1, and p(Z) is the probability distribution function of the observations. $p\left(\theta |Z\right)$ is a likelihood function.

We used the initial ranges of model parameters proposed by [55], and details of the MCMC process are described in the references [36,56].

3. Results

3.1. Parameters’ SI variation with Parameter Range

In this experiment, the model output was GPP and three-parameter ranges were set as listed in Table 1: $\pm$10% and $\pm$30% perturbation of the values and between the minimum and maximum values. It was found that parameter sensitivities varied with the ranges (Figure 1). For the range of $\pm$10% perturbation of value, only $f_{ci}$ and $V_m^{25}$ showed the highest SIs on the GPP (Figure 1a). The two parameters’ MSIs all exceeded 0.1, and their TSIs were 0.67 and 0.22 respectively. For this experiment, however, the other 14 parameters had a relatively weak influence on the GPP. These observations indicated that in the given range, $V_{\mathrm{m}}^{25}$ and $f_{ci}$ dominated model behavior. For the range of $\pm$30% perturbation of value, Parameters $f_{ci}\mathrm{and}$ ${\alpha }_q$, for which the TSIs were 0.54 and 0.43, were the most sensitive on GPP (Figure 1b). When the ranges were determined by the minimum and maximum values, the sensitivities parameters were as shown in Figure 1c. Parameters $f_{ci}\mathrm{and}$ ${\alpha }_q$, for which the TSIs were 0.51 and 0.49, were the most sensitive on GPP, and the TSIs of $V_m^{25}\mathrm{and} E_V{}_m^{}$ were 0.19 and 0.16, as well those of the other parameters.

By comparing the two cases, it can be seen that the ranges of the parameter significantly influenced their sensitivities and importance order [18]. The range selected in the first experiment was narrower and limiting some important but small parameters, such as ${\alpha }_q$, so the range was given in the second experiment was suggested for FBEM. Overall, in the two cases, the common findings were that the most sensitive parameters were $f_{ci}$, ${\alpha }_{\mathrm{q}},V_m^{25},\mathrm{and}{E}_{\mathrm{V}}{}_m^{}$, which was similar to the results of the earlier studies [57,58]. $V_{\mathrm{m}}^{25}$ was the key parameter for estimating photosynthesis and respiration, its inaccurate estimation may lead to non-negligible uncertainty in GPP estimation [59]. In addition, $f_{ci}$, ${\alpha }_q$ were a highly sensitive parameter, and relative experiments were also focused on them [60].

3.2. Parameter Sensitivity for GPP and NEE

For different model outputs, the importance of a given parameter may vary. This experiment examined differences in parameters’ SIs for NEE and GPP, and the results were shown in Figure 2.

As the model output was the GPP, the most sensitive parameters were ${\alpha }_{\mathrm{q}}$, $f_{ci}$, $V_{\mathrm{m}}^{25}$, and $E_{\mathrm{V}}{}_m^{}$. However, if the model output changed to NEE, ecosystem respiration ($Q_{10})$, whose TSI was 0.40, was the most sensitive parameters and was followed by the temperature dependence of ecosystem respiration at 0 °C (${ER}_0)$ and $E_{\mathrm{V}}{}_m^{}$. Both $f_{ci}$ and ${\mathrm{E}}_{\mathrm{V}}{}_m^{}$ were sensitive to GPP and NEE. ${\alpha }_q$ and $V_{\mathrm{m}}^{25}$ were sensitive to the GPP but not to the NEE, and ${ER}_0$ and $Q_{10}$ were only sensitive to the NEE. The results indicated that the sensitive parameters and sensitive degree varied for different model outputs.

3.3. Temporal Characteristics of Parameter SA

The total sensitivity index variation of the primary model parameters with the DOY for two outputs, namely GPP and NEE during the entire growing season was shown in Figure 3 and Figure 4.

As for GPP as shown in Figure 3, the TSI of $E_{\mathrm{V}}{}_m^{}$ and $D_0$ obviously changed over time. TSI of $E_V{}_m^{}$ decreased from DOY10 to DOY150, and was less than 0.01 from DOY160 to DOY240, and then increased from DOY260 to DOY340. SIs of $D_0$ was less than 0.01 before DOY100, and increased to the maximum at DOY220, and then declined. The SIs of the main parameters ${\alpha }_{\mathrm{q}}$ and $f_{ci}$ vary weakly with DOY, and TSI of $V_m^{25}$ did not change with DOY, So TSI of $V_m^{25}$ was not shown in Figure 3. As for NEE (Figure 4), $E{R}_0$and $Q_{10}$ played key roles, and the variation in them correlated with time. TSI of ${\mathit{ER}}_0$ and $Q_{10}$ were very small between DOY110 and DOY210, and increased quickly to 0.423 and 0.463, respectively.

The analysis of temporal change of SIs showed that the parameters were not invariable and changed prominently in different growth periods.

The study found that the variation in sensitivity had a trade-off effect. The ${\alpha }_q$, $f_{ci}$, and $E_V{}_m^{}$ declined with time, whereas ${\mathit{ER}}_0$ and $Q_{10}$ increased and became more sensitive from DOY110 to DOY210. After DOY210, ${\mathit{ER}}_0$ and $Q_{10}$ declined, and $f_{ci}$, ${\mathrm{E}}_{\mathrm{V}}{}_m^{}$ and ${\alpha }_q$ increased.

3.4. Parameter Uncertainty Analysis based on Sensitivity Analysis

To analyze the sources of uncertainty in the parameters, two experiments were designed for comparison. Based on the sensitivity analysis, four parameters $f_{ci}$, $V_m^{25}$, ${\alpha }_q$, and ${\mathrm{E}}_{\mathrm{V}}{}_m^{}$ had the highest sensitivity when the model’s output was the GPP. We designed two parameter estimation experiments to analyze uncertainty in the model based on the sensitivity analysis. Experiment I estimated four sensitive parameters and experiment II estimated 16 based on MCMC. The results are shown in Figure 5.

The results of GPP optimization are shown in Figure 5. The results obtained using the four most sensitive parameters through MCMC were almost identical to those obtained using 16 parameters. The RMSE values were 1.657 and 1.496, respectively, and the bias values of the simulated GPP were −0.03 and −0.055, respectively.

4. Discussion

It is worthy to further discuss the experimental results and related studies of sensitive parameters, the variability of parameter SIs due to different model outputs (GPP and NEE) and the relationship between parameter sensitivity and environmental factors at different growth stages.

4.1. Comparing to Previous Studies of Sensitive Parameters

In the flux-based ecological model, the sensitive parameters were $V_m^{25}$, ${\alpha }_q,\mathrm{and}{f}_C{}_i$ for GPP, and ${\mathit{ER}}_0$, $Q_{10}$, and $f_C{}_i$ for NEE. The variation of $V_m^{25}$ and SIs of $V_m^{25}$ were the main factor causing the uncertainty of GPP simulation. The results of related studies showed that varying with species, seasonal, nitrogen (N), and light intensity [58,59]. Walker found that global GPP varied between 108.1~128.2 PgC yr⁻¹ with nine implementations of the global distribution of $V_m^{25}$ [59], and greater effort was needed to pay close attention to the sensitivity of the variation in $V_m^{25}$ in the growing season and Fluxnet-stations. So, future research of the comparison of the SIs of parameters, especially $V_m^{25}$, would be in demand. $f_C{}_i$ was also important, and the most sensitive parameter overall for the GPP and NEE output, and a recent study indicates that $f_C{}_i$ was not a constant and varies along with environmental gradients [61], which was coincided with the parameter SA experiment. According to the above, the SA and other studies need to discover some parameters vary with the environment, the initial setting of the parameters as a constant was thus a crucial source of uncertainty [1,62].

4.2. Analysis of Variation of Parameters SIs

The variation of parameter SI was influenced by parameter range. The SA in 3.1 results showed that smaller parameters were more susceptible to ranges. The narrow range will limit the sensitivity index of the small-value parameters. Therefore, the value range of principle is very important for accurate SI. For this model, the recommended parameter ranges are given in Table 1. Similar to the previous study [19], the SIs of ±30% perturbation range was applicable when the specific range cannot be given.

The above results showed that the most sensitive parameters of GPP were different from that of NEE, and the reason was needed to discuss. GPP and NEE were the main products of photosynthesis and respiration processes, and their physiological processes are expressed by formulas. According to formulas (15) and (16), NEE was a direct result of ER, while ER was directly calculated with ER₀ and Q₁₀. Thus, the two parameters were more sensitive to NEE but insensitive to GPP.

4.3. Environmental Factors Analyses of Temporal Variance of Parameter SIs

The numerical experiment results showed that the SIs of key parameters have temporal variations on the GPP and NEE. More importantly, the mild trade-off effect was found in the parameter SIs during the growth season. An earlier study also found that sensitivities of several soil respiration-associated parameters have strong seasonal variations [29]. Therefore, temporal characteristics and the trade-off effect were essential features of parameters’ SIs.

As for GPP (Figure 3a), the SIs of $E_{\mathrm{V}}{}_m^{}$ and $D_0$ changed over the growing season. SIs of $E_V{}_m^{}$ increased from DOY10 to DOY150, weakened from DOY160 to DOY240, and increased from DOY160 to DOY240.

The study found that the variation in sensitivity had a trade-off effect (Figure 4). The SIs of ${\alpha }_q$, $f_{ci}$, and $E_V{}_m^{}$ declined with DOY, whereas SIs of ${\mathit{ER}}_0$ and $Q_{10}$ increased and became more sensitive from DOY110 to DOY210. After DOY210, SIs of ${\mathit{ER}}_0$ and $Q_{10}$ declined, and SIs of $f_{ci}$, ${\mathrm{E}}_{\mathrm{V}}{}_m^{}$, and ${\alpha }_q$ increased.

The dynamic evolution of the parameter SI indicated that the sensitive parameters from DOY10 to DOY340 for different outputs (GPP and NEE) were fluctuant. To further explore the causes of this phenomenon, the correlation between average air temperature(T_a), minimum air temperature (T_min), maximum air temperature (T_max), precipitation (P), PAR and showed SI of ${\mathsf{\alpha }}_{\mathrm{q}}$, D₀, $f_{ci}$ and ${\mathrm{E}}_{\mathrm{V}}{}_m^{}$ was analyzed. The result illustrates that only T_a was the more important environmental factor that influences the temporal variation of parameter, as shown in Figure 6 and Figure 7.

Parameters SI of $f_{ci}$ and ${\mathit{ER}}_0$ did not change significantly with an increase in T_a, this was due to the two parameters that were not closely related to temperature. So, the temperature change is not the cause of variation of these two parameters in the whole growth season. Unfortunately, what causes changes of SI of f_ci for NEE during the growing season Figure 7b was unclear and remained to be discussed in future studies.

Parameters SI of ${\mathrm{E}}_{\mathrm{V}}{}_m^{}$ decreased with an increase in T_a, and the trend of $Q_{10}$ was identical to that of T_a. SA experiment found that the effect of $Q_{10}$ for the model output was different. However, $Q_{10}$ was regarded as a constant in previous studies and had caught uncertainty in CO₂ simulation [63,64]. Therefore, it was worth further investigate that the parameters were taken as constants and the uncertainty analysis resulting from parameters on the model output in the existing model.

4.4. SA Improves the Effect of Parameter Optimization

The SA showed that each parameter has a different degree of influence on the result and varied with time. The results can help us better understand the effects of parameter optimization and reduce uncertainty in the model. Therefore, the perspective of parameter SA needs to be emphasized in uncertainty analysis, parameter optimization, and model structure.

Yuan design a two-step Bayesian inversion method to improve parameter estimates when model parameters converge undesirably [65]. In this method, whole parameters were optimized, and according to the convergence effect of the optimization, split into two or more steps. However, the SA recognized still more sensitive parameters and improved the two-step method. This study can help enhance the comprehension of the structure of ecological models. It can also be used to reformulate the controls of the GPP in next-generation ecological models.

5. Conclusions

A comprehensive SA was conducted for an ecosystem model, FBEM, based on a globally quantitative SA algorithm, EFAST. The sensitive parameters were distinguished from all parameters by quantifying the parameter’s SIs and ranking their importance of the FBEM. The following findings were obtained by various numerical experiments. First, the effects of the range of parameter variation on the SIs were significant. It is found in the experiment that the SIs of parameters were changed with a different value range, therefore, the suggested range of parameter sensitivity analysis in FBEM was given. Second, Sensitive parameters of the FBEM model for GPP and NEE were identified, the most sensitive parameters were ${\alpha }_q$, $f_C{}_i\mathrm{and}$ $V_m^{25}$ for GPP, while ${\mathit{ER}}_0$, $Q_{10}$ and $f_C{}_i$ were most sensitive for NEE. Second, an uncertainty analysis of the results of comparative experiments on parameter optimization was used to understand the performance of the model and the dynamics of the parameters. The results of this SA and uncertainty analysis are discussed in a more general framework. Furthermore, the temporal characteristics of the parameters’ SI on GPP and NEE were subsequently described: the sensitivity index clearly changed in different growth stages of the plants. SI of $E_V{}_m^{}$ gradually decreased while SI of D0 and ${\alpha }_q$ parameters were changing with a different value range, kept increasing at the beginning of the growth season with the model output GPP. SI of $f_{ci}$, Q₁₀ and $E_V{}_m^{}$ were variation in the growth season for NEE. The reasons for SI variation were discussed, temporal characteristic of ${\alpha }_q$, Q₁₀ and $E_V{}_m^{}$ were found to be mainly explained by the change of T_a. Moreover, more importantly, the mild trade-off effect was observed in the variation of the parameter SIs during the growing season in the EBEM. This study provided an improved understanding of the uncertainty of the ecological model caused by the parameters’ sensitivity and it also insights into potential approaches for the improvement of GPP simulation.