Improving the Simulation Accuracy of the Net Ecosystem Productivity of Subtropical Forests in China: Sensitivity Analysis and Parameter Calibration Based on the BIOME-BGC Model

Abstract

Subtropical forests have strong carbon sequestration potential; however, the spatiotemporal patterns of their carbon sink are unclear. The BIOME-BGC model is a powerful tool for forest carbon sink estimation while the numerous parameters, as well as the localization, limit their application. This study takes three typical subtropical forests (evergreen broadleaf forest, EBF; evergreen needleleaf forest, ENF; and bamboo forest, BF) in China as examples, assesses the sensitivity of 43 ecophysiological parameters in the BIOME-BGC model both by the Morris method and the extended Fourier amplitude sensitivity test (EFAST), and then evaluates the net ecosystem productivity (NEP) estimation accuracy based on the dataset of the fiveFi long-term carbon flux sites of those three typical forests from 2000 to 2015. The results showed that (1) both sensitivity analysis methods can effectively screen out important parameters affecting NEP simulation while the Morris method is more computationally efficient and the EFAST is better in the quantitative evaluation of sensitivity. (2) The highly sensitive parameters obtained using the two methods are basically the same; however, their importance varies across sites and vegetation types, e.g., the most sensitive parameters are k for the EBF and ENF and R_act25 for the BF, respectively. (3) The optimized parameters successfully improved the NEP simulation accuracy in subtropical forests, with average correlation coefficients increased by 25.19% and normalized root mean square error reduced by 21.74% compared with those simulated by original parameters. This study provides a theoretical basis for the optimization of process model parameters and important technical support for accurate NEP simulations of subtropical forest ecosystems.

1. Introduction

Forests are the largest carbon reservoirs in terrestrial ecosystems and help maintain global carbon balance and mitigate global climate change [1,2]. Net ecosystem productivity (NEP), an indicator of carbon sinks in terrestrial ecosystems, is the difference between the net primary productivity (NPP) of vegetation and heterotrophic respiration [3] and can quantitatively describe the net carbon exchange between terrestrial ecosystems and the atmosphere.

Models coupled with remote sensing capabilities merge the inherent advantages of remote sensing, such as real-time, dynamic, broad spatially synchronous monitoring, and multiresolution attributes with the quantitative simulation capabilities of models that address diverse aspects of the ecosystem’s carbon cycle. This fusion facilitates expansive simulations of the carbon cycle processes, establishing it as a paramount methodology for extensive forest carbon sink monitoring [4]. The process model is developed based on vegetation physiological ecology and ecosystem processes and functions [5]. As a paradigmatic ecosystem process model, BIOME-BGC simulates physiological and ecological processes, such as photosynthesis, respiration, and decomposition, across a range of scales [6]. Its intrinsic mechanistic rigor, high degree of integration, and excellent extrapolation capabilities [7] render it an invaluable asset for terrestrial ecosystem studies, with broad applications both within and beyond national borders [8,9,10]. However, many of the model’s input parameters cannot be derived from real measurements and errors in parameter estimation lead to large uncertainties in the simulation results. Therefore, the targeted optimization and calibration of key parameters can significantly reduce the workload of model calibration and improve model accuracy and predictive capability [11].

Global sensitivity analysis can simultaneously test the influence of multiple parameter changes and interactions between parameters on the model output and is necessary for model optimization and parameter calibration. Qualitative global sensitivity analyses provide a qualitative measure of the effect of model parameters on simulation outcomes and include the Fourier amplitude sensitivity test (FAST) [12] and the Morris method [13]. Quantitative global sensitivity analysis calculates the contribution rate of model parameters to the uncertainty of simulation outcomes and mainly includes the Sobol method [14] and the extended Fourier amplitude sensitivity test (EFAST) [15]. The Morris method can screen out sensitive parameters under a small sample size and has good applicability to models with many analysis parameters and large computational loads. It is a qualitative global sensitivity analysis method widely used at present. Wang et al. [16] employed the Morris method to identify the most critical parameters in a C-RIVE model under various nutrient conditions. DeJongea et al. [17] demonstrated that the computationally less costly Morris method can effectively screen for sensitive parameters of the CERES-maize crop model under different irrigation conditions. The EFAST combines the capabilities of the Sobol method, which calculates the parameter interactions, and the efficient sampling of the FAST. The EFAST analyzes not only the impact of each parameter’s variation on model output but also the influence of interactions between multiple parameters on model outcomes [18], finding widespread application in numerous nonlinear models. Li et al. [19] used the EFAST to assess the contributions of various input parameters in the DSSAT-CERES model to wheat yield and quality. Gilardelli et al. [20] employed the EFAST to evaluate the impact on the yield of parameter value variations in crop models simulating different crops under climate change. The assessment of sensitive parameters affecting the model using suitable methods is a prerequisite for model calibration and optimization, as well as the basis for localized application of the model.

The NEP of subtropical forest ecosystems in the East Asian monsoon region, which plays an important role in the global carbon cycle and carbon sink function, is 0.72 Pg C·a⁻¹ [21]. China is an important distribution area for subtropical forests and evergreen broad-leaved forests and evergreen coniferous forests are widely distributed in subtropical China [21]. The bamboo forest is known as the second-largest forest in the world and its carbon stock accounts for ~0.94% of the global forest ecosystem carbon stock [22]. China is known as the “world’s bamboo kingdom”, with abundant bamboo resources widely distributed in subtropical areas, such as Zhejiang, Fujian, and Jiangxi.

Simulations of the carbon and water cycles of various ecosystems found great uncertainty while monitoring subtropical forest carbon sinks, mainly for the following reasons. (1) When the BIOME-BGC model has been used to simulate the carbon cycle, the bamboo forest has not been treated as a separate forest, introducing uncertainty into the simulation. (2) Many studies have reported sensitivity analyses of BIOME-BGC parameters. Yan et al. [23] applied the EFAST to analyze and obtain sensitive parameters affecting the gross primary productivity (GPP) of the Changbaishan Forest. Tatarinov et al. [24] analyzed the effects of site and physiological–ecological parameters on the NPP of temperate forests in Europe. However, previous studies have predominantly focused on discerning the effects of parameter modifications on outputs, such as GPP and NPP. Notably, there is a paucity of evidence of the implications of model parameter fluctuations for NEP outputs, which amplifies the uncertainties in the simulation outcomes. (3) The calibration parameters of the BIOME-BGC model vary for different vegetation types under different environmental and site conditions. Raj et al. [25] investigated the sensitivity of the parameters in the BIOME-BGC model to the carbon fluxes of Douglas fir in the Speulderbos forest in the Netherlands and found that FLNR and FRC:LC (New fine root C: leaf C) were the key factors. Liu et al. [26] found that k was the most sensitive to estimate carbon fluxes in rubber forest ecosystems on Hainan Island. And FLRN was a low-sensitivity parameter but FRC:LC was not found to be the key control parameter. Kkumar and Raghubanshi [27] showed that FLNR and FRC:LC contributed less to the modeling of carbon fluxes in dry tropical forests. However, these studies suggest that the results of sensitivity analyses in other regions cannot be directly applied to simulate the carbon cycle in subtropical forests because of the climate, soil, and vegetation types.

In this study, we investigated the sensitivity of the physiological ecological parameters in the BIOME-BGC model to simulate the NEP of three forest ecosystems at five sites in subtropical China using the Morris method and EFAST. We aimed to screen key control parameters, improve the accuracy of model calibration, and provide scientific support for accurately assessing the carbon sink capacity of subtropical forests.

2. Materials and Methods

2.1. Study Area

Subtropical China (22–34° N, 98–122° E) is located south of the Qinling and Huaihe rivers, north of the Leizhou Peninsula, and east of the Hengduan Mountains, accounting for approximately 1/4 of the country’s national land area (Figure 1). Subtropical China belongs to the humid monsoon zone on the east coast, with the average annual temperature ranging from −1 to 24 °C and the average annual precipitation ranging from 450 to 2125 mm. The terrain is low in the west and high in the east and the soil type is reddish-yellow loam. Typical subtropical forests are evergreen broadleaf forests (EBFs), evergreen needleleaf forests (ENFs), and bamboo forests (BFs). Therefore, five flux observation sites were selected in this study for these three forests to participate in the parameter sensitivity analysis and parameter optimization, including Dinghushan (DHS), Tianmushan (TMS), Qianyezhou (QYZ), Anji (AJ), and Taihuyuan (THY). The information in

Table 1. Information of five flux sites used in this study; the three typical forests include the evergreen broadleaf forest (EBF), evergreen needleleaf forest (ENF), and bamboo forest (BF).

Site Name	Lat (°N)	Lon (°E)	Plant Functional Type
Dinghushan (DHS)	23.17	112.53	EBF
Tianmushan (TMS)	30.32	119.48	EBF
Qianyanzhou (QYZ)	26.74	115.06	ENF
Anji (AJ)	30.28	119.40	BF
Taihuyuan (THY)	30.26	119.59	BF

.

2.2. Subtropical Forest Flux Data and Meteorological Factors

Because the flux observation station may have problems, such as data loss caused by equipment failure, in order to ensure the accuracy and continuity of the observation data, we selected the data from 2000 to 2015, which are more complete than other years. Therefore, daily-scale NEP and meteorological data were collected from five observation stations from 2000 to 2015 and the meteorological factors were daily maximum temperature (°C), daily minimum temperature (°C), daily precipitation (mm), daily relative humidity (%), and total solar radiation (WJ-m⁻²). Data from the DHS, TMS, and QYZ sites were obtained from ChinaFlux (http://chinaflux.org, accessed on 10 September 2023). Data from the AJ and THY sites were obtained from the Carbon Flux Observing System, including the open-circuit vorticity correlation system and the conventional meteorological observation system. A detailed account of the instrumentation and the associated apparatus was provided by Mao et al. [28]. The observations from the carbon flux tower were sampled at a frequency of 10 Hz and processed into daily-scale NEP as follows. The flux data were corrected using the method described by Papale et al. [29] and processed using EdiRe software. Data interpolation of missing data was performed using the nightly data-based approach [30,31], whereby if the proportion of missing data to the total daily data volume was >20%, then the data for that day were discarded; if it was <20%, then the daily average value of NEP was multiplied by 24 to obtain the daily value of NEP. Daily meteorological data are half-hourly statistically recorded by the conventional meteorological observation system into daily-scale data.

2.3. Soil and Topographic Data

The soil data included the percentage of soil silt content and available soil water-holding capacity, which were calculated based on the soil texture obtained from the Harmonized World Soil Database (HWSD1.2) for the subtropics. Subtropical soil texture data were obtained from the Chinese soil dataset with a spatial resolution of 0.05° provided by the HWSD1.2. The topographic data included the Chinese subtropical digital elevation model and slope and direction data. Digital elevation model data were obtained from The ASTER Global Digital Elevation Model version 3 [32] and the slope and slope direction were calculated using ArcMap10.2. Finally, corresponding soil and terrain data were extracted based on the latitude and longitude of the flux observation sites.

2.4. BIOME-BGC Model Description

BIOME-BGC is a physiological–ecological process model that evolved from the forest-BGC model [33] (Figure 2). The model applies mechanistic modeling of photosynthetic enzymatic reactions [34], using a two-leaf model for the canopy and simulating GPP through the Farquhar Photosynthesis and Leuning Conductance models. Autotrophic respiration is divided into two components: maintenance and growth. Maintenance respiration was calculated using a model that relates plant tissue nitrogen content to the Q10 relationship. Growth respiration was obtained by multiplying the photosynthetic fixed organic matter of each plant tissue by a fixed coefficient. NEP is the difference between NPP and heterotrophic respiration; NPP is obtained by subtracting autotrophic respiration from GPP. Heterotrophic respiration originates mainly from the decomposition of soil organic matter and plant litter. The input data driving this model consist of three parts: (1) initialization files, such as elevation data, soil data, and atmospheric CO₂ concentration; (2) meteorological data at a daily scale, including minimum and maximum temperatures, precipitation, solar radiation, and saturated water vapor pressure; and (3) physiological and ecological parameters.

The BIOME-BGC model adopts a constant rate of growth and dieback for evergreen forests, does not consider the effects of bamboo forest size or anthropogenic disturbances on carbon cycling, and is not suitable for simulating carbon cycling in bamboo forest ecosystems [35]. Mao et al. [36] improved the BIOME-BGC model by introducing a forest age factor, bamboo whip carbon pool, bamboo forest climatic processes, the carbon allocation factor, and four bamboo forest management measures (fertilizer application, bamboo shoot digging, bamboo tip hooking, and bamboo harvesting).

The model was run in two steps, firstly a spin-up simulation, starting with very low initial values of carbon and nitrogen and running until the ecosystem reached a steady state, judged by an annual change in the soil carbon pool of ≤0.0005 kg·m⁻² [37]. The spin-up results were used as the initial values for carbon, nitrogen, and water reservoirs in the normal simulation. In this study, the original BIOME-BGC model was used for the EBF and ENF carbon cycle simulations and the improved BIOME-BGC model was used for the BF carbon cycle simulation.

Figure 2. Structure of the BIOME-BGC model [38].

Figure 2. Structure of the BIOME-BGC model [38].

2.5. Sensitivity Analysis

Sensitivity analyses of the physiological and ecological parameters of BIOME-BGC were conducted using the Morris method and EFAST to qualitatively and quantitatively analyze the effects of the input parameters on NEP.

2.5.1. Morris Method

The Morris method calculates the sensitivity of each parameter by differentiation, assuming that the system model is $Y=$ f(x₁, x₂,…, x_k) and k is the number of dimensions of the model parameters. Morris’s method first maps the range of values of each parameter to the interval [0, 1] and divides it into p levels, which constitute the sample space of the k dimensions and p levels. The elemental effects of each parameter are defined as follows:

$$S_i=\frac{[y(x_1,\cdot \cdot \cdot ,x_{i-1},x_{i+\Delta },x_{i+1},\cdot \cdot \cdot ,x_k)-y(x)]}{\Delta }$$

(1)

where $\Delta$ is a value between $1/p-1$ and $1-1/(p-1)$.

Using this method to calculate the elemental effects for each variable one by one requires a total of $k+1$ operations on the model; however, a single calculation does not show the whole of the sample space so it is usually necessary to generate different samples to calculate the S_i several times and to find the mean µ_i and the standard deviation σ_i the formula is as follows:

$${\mu }_i=\frac{{\displaystyle \sum _{j=1}^r{S}_{ij}}}{r}$$

(2)

$${\sigma }_i=\sqrt{\frac{{\displaystyle \sum _{j=1}^r(S_{ij}}-{\mu }_{\mathrm{i}}{)}^2}{r}}$$

(3)

where r is the number of repetitions.

To eliminate the positive and negative due to non-monotonic functional relationships of µ_i values canceling each other out, ${\mu }_i^{\ast }$ is usually used for evaluation [39]; see Equation (4):

$${\mu }_i^{*}=\frac{{\displaystyle \sum _{j=1}^r|S_{ij}}|}{r}$$

(4)

where ${\mu }_i^{\ast }$ reflects the sensitivity of x_i and the larger the ${\mu }_i^{\ast }$, the greater its influence on the model output results; σ_i reflects the interaction between x_i and other parameters.

2.5.2. EFAST

The EFAST quantifies the contribution of parameters to the variance of the model results. The model was assumed to be $Y=f\left(x\right)=f$(x_i,j,k,…) and the input parameter was (x_i,j,k,…) By estimating the sensitivity of the variance contribution of each input parameter change to the output $Y$, the total variance of the model output V_Y can be decomposed as follows [40]:

$$V_Y={\displaystyle \sum _i{V}_i}+{\displaystyle \sum _{i\ne j}{V}_{ij}}+\cdot \cdot \cdot +{\displaystyle \sum _{i\ne j\ne \cdot \cdot \cdot \cdot \cdot \cdot \ne k}{V}_{ij\cdot \cdot \cdot k}}$$

(5)

where V_Y is the total variance of the model results $Y$; V_i is the variance of x_i; and $V_{ij}~V_{ij\dots k}$ is the variance of the parameter interactions.

$$V_i=V[E(Y/x_i)]$$

(6)

where $E(Y/x_i)$ is the conditional expectation of $Y$ over $x_i$ and $V\left[E\right(Y/x_i\left)\right]$ is the variance of the conditional expectation of $Y$ over x_i.

$$V_{ij}=V[E(Y/x_i,x_j)]-V_i-V_j$$

(7)

where $E(Y/x_i,x_j)$ is the conditional expectation of $Y$ over $x_i$ and $x_j$ and $V\left[E\right(Y/x_i,x_j\left)\right]$ is the variance of $Y$ over $x_i$ and $x_j$.

The variance of the conditional expectation is called the main effect, which reflects the significance of $x_i$ to the variance of the model results $Y$. Therefore, the sensitivity index can be obtained by calculating the ratio of the variance of each parameter and parameter interaction to the total variance. The first-order sensitivity index S_i can reflect the direct contribution of $x_i$ to the total variance of the model results and can be defined as:

$$S_i=\frac{V_i}{V_Y}$$

(8)

Similarly, the second-order sensitivity index S_ij for $x_i$ is defined as:

$$S_{ij}=\frac{V_{ij}}{V_Y}$$

(9)

The total sensitivity index S_iT is the sum of the sensitivity indices of each order, reflecting the direct contribution of $x_i$ and the indirect contribution to the total variance of the model output through the interaction coupling between the parameters:

$$S_{iT}=S_i+S_{ij}+\cdot \cdot \cdot \cdot \cdot \cdot +S_{ij\cdot \cdot \cdot \cdot \cdot \cdot k}$$

(10)

The greater the sensitivity index of $x_i$; the greater the contribution of the parameter to the model results $Y$, both directly and indirectly through the interaction coupling between parameters; and the greater the explanatory power for the variance of changes in the model results.

2.5.3. Parameter Value Establishment

This study focused on the physiological and ecological parameters of this model. First, the parameters involved in the sensitivity analysis were selected and their baseline values were determined. The fire mortality rate was set to zero and the parameters were used to characterize markers, such as vegetation type or photosynthetic characteristics; parameters with fixed combination ratios were not included. Finally, 38, 38, and 43 parameters were selected for the EBF, ENF, and BF to participate in the sensitivity analysis, which were categorized in this study (

Table 2. Photosynthesis-related parameters involved in sensitive analyses of the BIOME-BGC model.

Symbol	Description	Unit
Photosynthesis biophysics parameters
FLNR	Fraction of leaf N in Rubisco	kgNRub·kgNleaf−1
Ko25	Michaelis constant of oxidation reaction at 25 °C	-
Ract25	Rubisco activity at 25 °C	μmol·mg·Rubisco−1·min−1
Kc25	Michaelis constant of carboxylation reaction at 25 °C	-
Q10kc	The Q10 temperature coefficient of kc	-
Q10ko	The Q10 temperature coefficient of ko	-
Q10Ract	The Q10 temperature coefficient of Rubisco	-
Allocation of carbon parameters
FRC:LC	New fine root C: leaf C	kgC·kgC−1
LWC:TWC	New live wood C: total wood C	kgC·kgC−1
SC:LC	New stem C: new leaf C	kgC·kgC−1
CRS:SC	New coarse root C: stem C	kgC·kgC−1
CGP	Current growth proportion	-
Canopy structure biophysics parameters
Wint	Water interception coefficient	LAI−1·d−1
k	Light extinction coefficient	-
SLA	Average specific leaf area	m·kgC−1
SLAshd:sun	Ratio of shade SLA: sunlit SLA	-
LAIall:proj	Ratio of all sides to projected leaf area	-
Stomatal conductance biophysics parameters
Gsmax	Maximum stomatal conductance	m·s−1
Gbl	Boundary layer conductance	m·s−1
LWPf	Leaf water potential: completion of gs reduction	MPa
LWPi	Leaf water potential: start of gs reduction	MPa
VPDf	Vapor pressure deficit: completion of gs reduction	Pa
VPDi	Vapor pressure deficit: start of gs reduction	Pa
Gcut	Cuticular conductance	m·s−1

and

Table 3. Respiratory-related and other physio-ecological parameters involved in sensitivity analysis of the BIOME-BGC model.

Symbol	Description	Unit
Heterotrophic respiration biophysics parameters
C:Nlitter	C: N of falling leaf litter	kgC·kgN−1
LFG	Litterfall period	-
LWT	Annual live wood turnover fraction	a−1
Maintenance respiration biophysics parameters
MRpern	Maintenance respiration in kg C/day per kg of tissue N	kgC·kgN−1·d−1
C:Nleaf	C: N of leaves	kgC·kgN−1
C:Nfr	C: N of fine roots	kgC·kgN−1
C:Ndw	C: N of dead wood	kgC·kgN−1
C:Nlw	C: N of live wood	kgC·kgN−1
Vegetation chemical parameters
Tt	Transfer growth period	-
Llab	Leaf litter labile proportion	-
Llig	Leaf litter lignin proportion	-
FRlab	Fine root labile proportion	-
FRlig	Fine root lignin proportion	-
DWlig	Dead wood lignin proportion	-
Management measure parameters
Pdsw	Excavation percentage of winter bamboo shoots	%
Pobtr_total	The ratio of hook tip carbon storage to total leaf carbon storage	%
SCRages	Selective cutting ratio of each age	-
Fer	Apply fertilizer	kgN·hm−2
Pdss	Percentage of bamboo shoots harvested	%

). The baseline values of the parameters were obtained using three methods: (1) some of the parameters were obtained by actual measurements in the sample plots, (2) other plant physiological parameters were obtained by conducting literature searches for each parameter and calculating the mean values, and (3) the parameter values of the EBF and ENF provided by White et al. [41] were used for parameters that could not be determined from the literature. Secondly, the value range of each parameter was determined because BIOME-BGC did not give the specific value range of each parameter so we floated 20% on the basis of the benchmark value of each parameter [41], setting the value range of the parameters as (x − Δx, x + Δx), where x is the benchmark value of the parameter Δx = 0.2x using a uniform distribution; the parameters are independent of each other, as shown in Supplementary Materials (Tables S1 and S2) [41,42,43,44,45,46,47,48,49,50,51,52].

2.5.4. Sensitivity Analysis Steps

The sensitivity analysis was conducted with SimLab (Version 2.2) (Figure 3). The value range, distribution form, and sampling time were set for each parameter. Using the Monte Carlo method to generate parameter samples by random sampling of the parameters in Table 2, the distribution form of each parameter was assumed to be uniform [53]. The sampling number of the Morris method was $r\times (K+1)$ ($K$ is the number of parameters; $r$ is the number of sampling repetitions, where r_min = 4 and r = 10) and the total number of sampling times for the five observation points was 2050. The number of samples for the EFAST was $r\times K$ ($K$ is the number of parameters; n is the number of sampling repetitions, where the method specifies that the results are valid when n ≥ 65; in this study, n = 150) so the total number of samples for the five observation points was 30,000.

For exporting and reformatting the sampling results, combined with meteorological, soil, and topographic data, different parameter combinations were used to drive the BIOME-BGC model to simulate mean NEP values from 2000 to 2015. Sensitivity analyses were conducted using the Morris method and EFAST to determine the sensitivity of the different parameters. For the results of the Morris method, the ${\mu }^{\ast }$ values of all parameters were accumulated to obtain the average value ${\mu }_{avg}$; when ${\mu }^{\ast }$ >${\mu }_{avg}$, the parameter is sensitive. For the results of the EFAST, the sensitivity index was divided into three levels, >0.2 for the high-sensitivity parameter, >0.1 and ≤0.2 for the medium-sensitivity parameter, and the rest of them were insensitive parameters [24].

2.6. Model Accuracy Validation

In this paper, the low-sensitivity parameters are assigned fixed values and the screened high-sensitivity parameters are randomly sampled in the interval of plus or minus 100% based on the baseline values and then crossed with each other to obtain different parameter combinations. The root mean square error (RMSE) of the NEP results simulated with each set of parameters is compared with the observed values and the parameter combination with the smallest RMSE is selected as the input parameters of the final BIOME-BGC model.

The simulation results were validated using daily flux observations. The model simulation accuracy was verified using the correlation coefficient (R) between the simulated and measured values and the NRMSE (%) using the following equations:

$$R=\frac{{\displaystyle \sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y})}}{\sqrt{{{\displaystyle \sum _{i=1}^N(x_i-\overline{x})}}^2}\sqrt{{\displaystyle \sum _{i=1}^N{(y_i-\overline{y})}^2}}}$$

(11)

$$NRMSE=\frac{RMSE}{M_{\max }-M_{\min }}\times 100\%$$

(12)

where N is the number of observed data, $x_i$ is the simulated value, $\overline{x}$ is the mean value of the simulated value, $y_i$ is the measured value, $\overline{y}$ is the mean value of the measured value, and M_max and M_min are the maximum and minimum values of the measured dataset, respectively. Regaring R, the closer the value is to 1, the better the fit between the simulated and measured results and the smaller the value of the NRMSE, the smaller the deviation of the simulated value from the observed value, and the better the simulation result of the model.

3. Results

3.1. Uncertainty Analysis

The uncertainty of the input data was calculated using the EFAST, as shown in Figure 4. The EBF (DHS site) had the highest NEP values, with a distribution ranging from 392 to 852 gC∙m⁻²∙a⁻¹. The ENF (QYZ site) had a relatively scattered distribution of NEP, with a range of 150 to 690 gC∙m⁻²∙a⁻¹. The NEP distribution of the BF (THY site) was relatively concentrated.

Table 4. Statistical information of the annual average NEP for three forest types at five sites.

Statistics	DHS-EBF	TMS-EBF	QYZ-ENF	AJ-BF	THY-BF
Mean	551.90	612.79	449.89	558.19	529.23
Standard deviation	76.47	89.15	104.46	79.47	67.65
Coefficient of variation (%)	13.86	14.55	23.22	14.24	12.78

shows that there was no significant difference in the coefficient of variations (CV) between the two EBF sites, with a 13.86% CV for DHS and a 14.55% CV for TMS. The uncertainty in the ENF (CV, 23.22%) was significantly higher than that of the EBF and BF (CV = 14.24% and CV = 12.78%).

3.2. Morris Sensitivity Analysis

Figure 5 shows the sensitivity indices of the physiological and ecological parameters at each of the five sites (${\mathsf{\mu }}^{\ast },\mathsf{\sigma }$). For the EBF, k was the most influential parameter at both sites. In addition to k, the order of sensitive parameters at DHS was SLA, G_smax, SC:LC, Q_10Ract, and G_bl; for TMS, the sensitive parameters were in the order SLA, G_smax, G_bl, FLNR, and R_act25. The order of the top three parameters was the same for both sites; however, the order of the less-sensitive parameters differed. The sensitive parameters affecting the ENF at QYZ were k, MR_pern, K_c25, Q_10Rac, G_bl, FLNR, and SC:LC. Additionally, k determines the amount of light energy that can be sequestered by the canopy, which directly affects the total photosynthetic primary productivity of the stand and shows high sensitivity in both the EBF and ENF. For the BF, AJ and THY were similar in terms of topography, soil, climate, and other conditions and their sensitivity parameters and rankings were also highly similar. The sensitivity parameters and rankings of AJ were R_act25 > K_o25 > Q_10Rac > FLNR > FRC:LC > K_c25 > SCR_ages. Compared with AJ, THY had the fourth highest rank of K_o25 and the remaining sensitivities were consistent.

3.3. EFAST Sensitivity Analysis

The physiological and ecological parameters of S_i and S_iT for each of the five observation sites are shown in Figure 6.

For the EBF, the highly sensitive parameters affecting the NEP simulation results at DHS were k, SLA, G_smax, and SC:LC and the moderately sensitive parameters were G_bl, Q_10Ract, and FRC:LC. Among them, the S_iT of k was much larger than that of S_i; the S_iT of G_smax, G_bl, and Q_10Ract was >0.1 and that of S_i was <0.1, suggesting that these four parameters mainly influence the NEP simulated values through interaction with other parameters. The S_iT of FRC:LC was <0; whereas, that of S_i was >0, indicating that FRC:LC directly affects the NEP simulation results. Additionally, k was shown to be a highly sensitive parameter that affected the NEP simulation results at TMS and the moderately sensitive parameters were FLNR, G_bl, G_smax, R_act25, C:N_leaf, and SLA. In addition, k, SLA, FLNR, and G_bl were common sensitive parameters at DHS and TMS. The EBF physiological and ecological parameters of different sites showed similar combinations of sensitive parameters.

The high-sensitivity parameters affecting the NEP simulation results at QYZ were k and MR_pern and the medium-sensitivity parameters were Q_10Ract, G_bl, FLNR, SC:LC, and K_c25, where the S_iT of SC:LC was >0.1 and that of S_i was <0.1. This suggests that SC:LC does not directly affect the results but rather varies synergistically with the other parameters, affecting the simulated NEP.

For the BF, the high-sensitivity parameters affecting the NEP simulation results at AJ were K_o25, SCR_ages, R_act25, and FRC:LC and the medium-sensitivity parameters included Fer, K_c25, Q_10Ract, FLNR, and k. The high-sensitivity parameters affecting the NEP simulation results at THY were R_act25, Q_10Ract, and SCR_ages and the medium-sensitivity parameters included FRC:LC, K_o25, FLNR, and Fer. The differences in the S_i and the S_iT of the sensitive parameters at both sites were small and they mainly affected the NEP of the BF through direct effects. Among these, SCR_ages and Fer are bamboo-forest-specific parameters that significantly affected the NEP simulation results for the BF.

3.4. Comparative Analysis of Morris and EFAST Results

As shown in Figure 7, the combination of the two sensitivity analysis methods led to ten key sensitive parameters affecting NEP: k, SC:LC, SLA, FRC:LC, G_smax, MR_pern, K_o25, R_act25, Q_10Ract, and SCR_ages (for the BF only). The most sensitive parameter of the EBF and ENF calculated using the Morris method was k. The key control parameter of the BF (THY) was R_act25, which is consistent with the results of the EFAST. At DHS, the EFAST added FRC:LC over the Morris method. At TMS, C:N_leaf was added; at QYZ, the types of sensitive parameters were the same but the ordering was different. At AJ and THY, the EFAST added the bamboo-forest-specific parameter Fer over the Morris method. The types of sensitive parameters screened by the two methods that have a large impact on the NEP were highly similar; however, the order of importance of sensitive parameters was different. The EFAST generally calculates more sensitive parameters than the Morris method, partly because of the number of samples and partly due to the effects caused by the coupling of the parameters. The key control parameters vary from site to site and the key control parameters of each site show the same sensitivity at other sites with different degrees of sensitivity. Additionally, k is the common sensitive parameter in QYZ, DHS, TMS, and AJ and showed high sensitivity in the EBF (QYZ) and ENF (DHS and TMS), with a S_iT of 0.4 or more, the first in the Morris method. In addition to k, MR_pern was the most sensitive parameter for the ENF (QYZ) (S_iT = 0.34 and ${\mu }^{\ast }$= 0.54), SLA was the most sensitive parameter for the EBF (DHS) (S_iT = 0.39 and ${\mu }^{\ast }$ ranged from 0.65 to 0.75), and K_o25 and R_act25 were the most influential parameters for the BF (AJ and THY) (S_iT = 0.48, S_iT = 0.56 and ${\mu }^{\ast }$ ranged from 0.5 to 1).

3.5. BIOME-BGC Model Validation

The carbon flux values simulated by the parameter localized model are closer to the flux observations than the original model. The original model had unrealistic peaks while the accuracy was significantly improved by the parameter-localized model (Figure 8). The R of NEP simulated using the original model ranged from 0.45 to 0.66 and NRMSE ranged from 15.88 to 28.43% for the five sites; the R of NEP simulated using the parameter-optimized model ranged from 0.61 to 0.75 and NRMSE ranged from 14.83 to 21.09%. Compared with the original model, the average R of the localized model improved by 25.19% and the average NRMSE reduced by 21.74%. The localized model can better simulate the changing characteristics of carbon fluxes in subtropical forests.

4. Discussion

4.1. Uncertainty Analysis

Uncertainty in the BIOME-BGC model arises from many sources, including the model structure, input variables, and model parameters [54]. Uncertainty in the model structure may be due to an insufficient understanding of ecosystem processes and the simplification of ecological process mechanisms by modelers [55,56]. The effects of forest fires and anthropogenic activities [8], for example, on carbon and nitrogen water cycles were not considered in this simulation, which may generate uncertainty. Uncertainty in the input variables may be caused by measurement errors, insufficient data collection, or statistics [57,58]. There are two main reasons for uncertainty in the model parameters. First, some input parameters cannot be directly observed, especially those related to carbon and nitrogen ratios and partitioning. Second, obtaining a range of parameter variations is difficult (e.g., Q_10kc, C:N_leaf, MR_pern). For some parameters that were difficult to obtain, we reviewed the literature and directly used default model values, which may lead to uncertainty. The range of parameter values has a significant effect on the results of sensitivity analyses. A smaller baseline value is more likely to be affected by the parameter range and smaller parameter ranges limit their sensitivity indexes [59]. The wide range of variability in the NEP for the ENF indicates that when generating multiple sets of sample parameters, some combinations contained parameters that were outside the normal range of physiological parameters of the ENF (Figure 4 and Table 4), representing forests experiencing unsustainable or disturbed growth during the years simulated by the model, which led to lower NEP values in the simulations. However, sensitivity analyses are not affected by these conditions owing to the presence of available and useful spin-up phases [53]. Therefore, we controlled the parameter fluctuation amplitude to 20%, which ensured consistency in the parameter range to a certain extent.

When comparing the results of the variance-based method (e.g., EFAST) with those of the Morris method, the different forms of distribution of the input parameters can lead to uncertainties in the sensitivity indices [14,60]. When the input parameters of the models are independent of each other and obey a uniform distribution, the sensitivity indices show a high correlation. However, when the input parameters are arbitrarily distributed, the results will greatly differ between the methods and ${\mu }_i^{\ast }$ will not represent S_iT well [61,62]. In this study, the distribution form of the input parameters was set to be uniform, and the results showed that the types of sensitive parameters screened by the two methods were basically similar. However, there were differences in the ordering of the sensitive parameters, which was also reflected in the following studies. Ciffroy and Benedetti [63] analyzed the parameters of a geochemical morphology model using the two methods and the types of sensitive parameters calculated by the two methods were similar but the ordering was slightly different. Peng et al. [18] used the Morris method and EFAST to analyze the global sensitivity of the parameters of the DCSEM model and the sensitive parameters were similar. Xue et al. [64] applied the two methods to a process model (CROBAS model) and the sensitivity parameters were similar but with slight differences in importance. The differences between the results are, firstly, due to the different sample sizes used for the two methods (the sample size and computational volume of the EFAST are larger). Secondly, the results of the sensitivity analysis of the output variables are affected by the strength of the interactions between the parameters while the Morris method cannot judge the interaction between the parameters based on the magnitude of ${\mu }_i^{\ast }$, thus underestimating the sensitivity of some parameters. Both methods can effectively screen out the important parameters affecting the output of the model. The Morris method has a significant advantage when the sample size is small, parameters are numerous, and the time required for the model operation is long. The EFAST can further quantitatively analyze the contribution of each input parameter and the interactions between the parameters to the simulation results [18]. Therefore, for parameter sensitivity analyses of complex process models, the Morris method can be used for qualitative studies and the EFAST for quantitative evaluations.

4.2. Sensitive Parameters

Sensitivity analyses based on both the Morris method and EFAST showed that k, SC:LC, SLA, FRC:LC, G_smax, MR_pern, K_o25, R_act25, Q_10Ract, and SCR_ages (for the BF only) significantly affected the simulated NEP. However, there may be differences in parameter sensitivity owing to forest type, growing environment, geographic location, and other factors. Liu et al. [26] analyzed the parameter sensitivity of rubber trees in Hainan Province, China, and found that the key sensitive parameter affecting the carbon flux output of the model was k, followed by SLA and C:N_leaf, with FLNR, FRC:LC, and SC:LC also showing sensitivity. Kumar and Raghubanshi [27] reported k, SLA, SC:LC, and G_smax as highly sensitive parameters affecting carbon fluxes in dry tropical forests, with k being the most influential. Raj et al. [25] found that Douglas fir carbon fluxes in the Speulderbos forest in the Netherlands were consistently the most sensitive to FLNR, with the second most sensitive parameter being FRC:LC, followed by SC:LC and SLA. Yan et al. [23] simulated forest carbon fluxes in the Changbaishan region and showed that FLNR was the most sensitive; SLA and G_smax also showed high sensitivity but k did not. Ritika et al. [65] reported sensitivity analyses for Indian forests and yielded similar parameters, with changes in SLA values having the greatest impact on carbon flux simulations. In general, the main physiological and ecological parameters affecting ecosystem carbon fluxes in the BIOME-BGC model were those affecting photosynthesis and respiration; however, the maximum sensitivity parameters of their respective responses differed because of the environment, vegetation type, and other factors. This is also true for the results of the present study, which illustrates the reasonableness and reliability of the results.

For the study area and forest types used in this study, k showed high sensitivity in both EBF and ENF. Additionally, k reflects the extent to which the light passing through the vegetation layer is subjected to processes, such as absorption, scattering, and reflection. When k decreases, more light passes through the vegetation layer to reach the ground, thus affecting the intensity of photosynthesis and evapotranspiration in plants and contributing to an increase in carbon fixation and net carbon balance. Conversely, when k increases, the absorption, scattering, and reflection of light by vegetation increases, allowing less light to pass through the vegetation layer to reach the ground, thus reducing carbon fixation and the net carbon balance. In addition to k, MR_pern was the most sensitive parameter in the ENF; MR_pern represents the coefficient of the linear relationship between tissue N and maintenance respiration and is directly related to maintenance respiration, which, in turn, affects NEP [66]. SLA was the most influential parameter in the EBF and mainly affects plant photosynthesis by influencing the area and intensity of light received by the leaves. Changes in SC:LC values affect the calculation of leaf carbon content, subsequently affecting plant growth and photosynthesis. G_smax is an important parameter for plants to regulate gas exchange and photosynthesis. When G_smax decreases, the concentration of CO₂ in the leaf decreases, which, in turn, reduces the rate of photosynthesis [67]. Q_10Ract, K_o25, and R_act25 affect the photosynthetic efficiency of leaves by controlling the levels of photosynthesis-related enzymes that influence the output of carbon flux [68]. FRC:LC is also an important parameter and when FRC:LC increases (more carbon in the leaves is transferred to the fine roots), leaf growth is inhibited, which, in turn, affects the accumulation of carbon in the plants [69]. SCR_ages, a bamboo-forest-specific parameter, significantly reduces the aboveground carbon stocks in bamboo forests and has a non-negligible effect on stand productivity [70,71]. Atmospheric carbon enters the biosphere through photosynthesis and returns to the atmosphere through respiration and these two physiological processes are crucial to ecosystem carbon cycles [42]. Therefore, it is reasonable to assume that the parameters affecting photosynthesis and respiration are the most sensitive to the carbon flux output of the BIOME-BGC model and the parameters affecting photosynthesis and respiration should be considered first for correction when the model is extrapolated to other regions. The results of this study clarified the key sensitive parameters of carbon fluxes in different regions and vegetation types and provided data to support the accurate simulation of carbon fluxes and the correction of model parameters at the regional scale.

5. Conclusions

The key sensitive parameters affecting NEP in subtropical forest ecosystems were explored and analyzed in the BIOME-BGC model using the Morris method and EFAST and the localized BIOME-BGC model was validated using flux observation data. The following conclusions were drawn:

• Both the Morris method and EFAST can effectively screen out the important parameters affecting the output of the model. The Morris method has a significant advantage when the sample size is small, the parameters are numerous, and the computational effort is high; however, it is only a qualitative method of parameter sensitivity analysis. The EFAST method allows further quantitative analysis of the contribution of each input parameter and the interaction between parameters to the simulation results; however, its computational efficiency is lower than that of the Morris method. For parameter sensitivity analyses of complex process models, the Morris method is used for qualitative studies and the EFAST method is used for quantitative studies;
• The parameters k, SC:LC, SLA, FRC:LC, G_smax, MR_pern, K_o25, R_act25, Q_10Ract, and SCR_ages (only for the BF) significantly affected the simulated subtropical forest NEP. Priority should be given to these parameters in model parameter optimization and correction to reduce computation and improve model accuracy;
• Compared with the flux observation data, the parameter-optimized BIOME-BGC model significantly improved the simulation ability of the original model for the NEP of subtropical forest ecosystems in China; the average R of the NEP increased by 25.19% and the average NRMSE reduced by 21.74%.

This study deepens the understanding of the process model and provides a reliable solution for accurate spatial and temporal simulations of the subtropical forest carbon cycle.