Global Sensitivity Analysis of the LPJ Model for Larix olgensis Henry Forests NPP in Jilin Province, China

Abstract

Parameter sensitivity analysis can determine the influence of the input parameters on the model output. Identification and calibration of critical parameters are the crucial points of the process model optimization. Based on the Extended Fourier Amplitude Sensitivity Test (EFAST) and the Morris method, this paper analyzes and compares the parameter sensitivity of the annual mean net primary productivity (NPP) of Larix olgensis Henry forests in Jilin Province simulated by the Lund–Potsdam–Jena dynamic global vegetation model (LPJ model) in 2009–2014 and 2000–2019, and deeply examines the sensitivity and influence of the two methods to each parameter and their respective influence on the model’s output. Moreover, it optimizes some selected parameters and re-simulates the NPP of Larix olgensis forests in Jilin Province from 2010 to 2019. The conclusions are the following: (1) For the LPJ model, the sensitive and non-influential parameters could be identified, which could guide the optimization order of the model and was valuable for model area applications. (2) The results of the two methods were similar but not identical. The sensitivity parameters were significantly correlated (p < 0.05); parameter k_rp was the most sensitive parameter, followed by parameters α_m, α_a and g_m. These sensitive parameters were mainly found in the photosynthesis, water balance, and allometric growth modules. (3) The EFAST method had a higher precision than the Morris method, which could calculate quantitatively the contribution rate of each parameter to the variances of the model results; however, the Morris method involved fewer model running times and higher efficiency. (4) The mean relative error (MRE) and mean absolute error (MAE) of the simulated value of LPJ model after parameter optimization decreases. The optimized annual mean value of NPP from 2010 to 2019 was 580 g C m⁻² a⁻¹, with a mean annual growth rate of 2.13%, exhibiting a fluctuating growth trend. The MAE of the simulated value of LPJ model after parameter optimization decreases.

1. Introduction

Net primary productivity (NPP), as a critical variable characterizing the terrestrial ecological process, can reflect the plant’s production capacity in the natural state and play an essential function in the global carbon cycle [1,2,3,4]. There is a lack of comprehensive and systematic research on the dynamic simulation of NPP of Larix olgensis Henry forests in Jilin Province. Based on the 3-PG model, Xie et al. studied the dynamic changes of NPP with forest age in 30 fixed plots in Jilin Province during a rotation period [5]. Li et al. simulated the NPP of Larix olgensis forests in Jilin Province from 2000 to 2019 using the LPJ-DGVM model (the Lund–Potsdam–Jena dynamic global vegetation model, hereinafter referred to as the “LPJ model”) and proposed that the model parameters would bring great error to the simulation process [6]. The LPJ model, as an ecological process model, has been extensively used in NPP simulation in China and abroad [7,8]. Its advantage is that it can simulate the mechanism behind vegetation life and biogeochemical cycle processes, which has become an effective tool to simulate large-scale (regional to global) vegetation geographical distribution, productivity, and carbon balance, and thus predict the potential impact of climate change on the terrestrial ecosystem [9,10,11,12]. However, the LPJ model process is complex and contains several vegetation physiological and ecological parameters. So, the parameter values of the model cannot meet the heterogeneous conditions, such as climate and vegetation types, in various study areas. Some parameters must be further tested and rectified for extrapolation or localization. The calibration of each model parameter entails a huge workload, and the improvement effect of model accuracy is poor. Therefore, it is necessary to screen the sensitivity of the model’s critical parameters of the model to improve its efficiency before calibrating its parameters, which has also long been the research focus in the identifying and optimizing the critical parameters of the process model [13,14,15,16,17].

Sensitivity analysis of model parameters is an effective approach to solve this problem [18,19]. Sensitivity analysis explains the influence of an output variable on the change of a parameter by changing the value of the input variable (parameter), calculating the influence degree and the sensitivity degree of the input parameter on the output variable, and determining which parameters have a significant impact on the model output variable [20]. Generally, only the parameters that significantly impact the output variables are calibrated in model calibration. The parameters with low sensitivity are set as fixed values to improve the model calibration accuracy and reduce the data-processing workload [20,21,22]. Sensitivity analysis methods include the local and global sensitivity analysis methods. Local sensitivity analysis evaluates the local response of the model outputs; it is easy to operate and is extensively used in the sensitivity analysis of the model parameters, but ignores the coupling effect between the parameters. The analysis results are one-sided, such as OAT (one at a time), suitable for linear and near-linear models [22,23,24]. Meanwhile, global sensitivity analysis is an improvement of the local sensitivity analysis method; it considers the response in the whole parameter space and tests the influence of multiple parameter changes on the simulation results, which is more suitable for LPJ with a highly nonlinear model. At present, the commonly used global sensitivity analysis methods include the Sobel method [25], the Morris method [13,26], the Fourier Amplitude Sensitivity Test (FAST), and the Extended Fourier Amplitude Sensitivity Test (EFAST) [27,28,29,30]. The EFAST method and the Morris method are more widely used. The EFAST combines the advantages of the Sobel and the FAST to calculate the interaction of various parameters accurately. The Morris method calculates the interaction based on the idea of differentiation, and the advantage is that the number of calculations is small.

This study aims to research the impact of the parameters of the LPJ model on the simulation of Larix olgensis forests NPP in Jilin Province through the EFAST and Morris method, which are extensively used to measure the global sensitivity of the model parameters, to provide technical support for improving the calibration efficiency of the parameters and the model’s localized application.

2. Material and Methods

2.1. Site Description

The LPJ model is applied in Jilin Province (121.63–131.32° E, 40.87–46.30° N, Figure 1), which is located in the middle of Northeast China, in the north temperate zone, with a diverse climate and complex landform. The monthly mean temperature is −15–22 °C and the annual precipitation is 400–800 mm. From east to west, there is an obvious climate change trend from humid to semi-humid and then to semi-arid. Taking the Dahei Mountain as the boundary, the province is divided into western plains and eastern mountains, with the topography being high in the southeast and low in the northwest. Forests are mainly distributed in the middle and eastern mountainous areas and low mountains and hilly areas. The tree species are mainly coniferous trees such as Larix olgensis, Pinus koraiensis Sieb. et Zucc., Picea asperata Mast., Pinus sylvestris var. mongolica Litv., and broad-leaved trees such as Quercus mongolica Fisch. ex Ledeb., Betula platyphylla Suk., and Populus spp. The shrubs mainly include Salix rorida Lacksch., Spiraea salicifolia L., Lonicera japonica Thunb., Corylus heterophylla Fisch. ex Trautv., etc., and the herbs are mainly Gramineae and Compositae.

2.2. Climate Data

The climate data were collected from the China Meteorological data network (http://data.cma.cn (accessed on 1 December 2019)), with the monthly data of 41 meteorological stations in Jilin Province and its surrounding area within 100 km from 2000 to 2019 (Figure 2). It includes monthly mean temperature, monthly precipitation, monthly wet days, and monthly sunshine percentage. The Kriging interpolation method was used to batch interpolate discrete station data into continuous grid data.

2.3. Description of the LPJ Model and Input Parameter Setting

2.3.1. Description of LPJ-DGVM

The LPJ model was developed based on the BIOME series models to simulate regional and even global terrestrial ecosystems. Beginning from the principle of vegetation dynamics, the model that simulates the process is based on plant photosynthetic biochemical reactions, canopy energy balance, allometric growth distribution of biomass, and soil water balance. Then, the fast reactions processes of ecosystem plants, such as photosynthesis, respiration, and evaporation, were simulated grid-by-grid and year-by-year in the same model framework. In the meantime, the slow reactions processes, including resource competition, tissue turnover, allometric distribution, and litter decomposition, were simulated. In addition, the natural disturbance factors and death effects of the ecosystem were considered to establish the simulated population and individual deaths. These processes are primarily affected by the simulated area’s environmental conditions, species composition, and species attributes [9,10]. The LPJ model can simulate eight plant functional types (PFT) and two herb vegetation functional types according to different phenological, physiological, and morphological characteristics. To realize the scale conversion from vegetation individual to population in the simulation unit, the model takes the mean individual of PFT as the calculation unit, and each simulation unit is composed of one or more PFTs. According to the scale conversion rules, the carbon stocks of woody plants are allocated to four vegetation tissue pools (leaves, sapwood, heartwood, and fine roots), aboveground or underground litter pools, and two soil carbon pools. Vegetation photosynthesis is calculated on the canopy scale, and the NPP is calculated according to organic matter accumulated by vegetation photosynthesis (gross primary productivity, GPP) minus its respiratory consumption (R_m) (Equation (1)) [10].

$$\mathrm{NPP}=\mathrm{GPP}-R_m$$

(1)

The model is driven by climate data (e.g., monthly mean temperature, monthly precipitation, monthly wet days, and monthly sunshine percentage), and must input soil attribute data and atmospheric CO₂ concentration data in the simulation area (see Section 2.2 for the acquisition and processing of meteorological data). The soil attribute data were derived from the International Food and Agriculture Organization (FAO) soil dataset. Considering a short simulation time interval, the global mean atmospheric CO₂ concentration, as observed by the Mauna Loa station in Hawaii since 1951, was used as the historical mean concentration of the study area. The most sensitive climatic factor for vegetation distribution is the monthly mean temperature of the coldest month. The species appear only when the actual monthly mean temperature is within the growth temperature range and reaches the lower limit of the accumulated temperature required for growth. Based on the secondary unit (vegetation type) of the classification system in the Vegetation Map of the People’s Republic of China [31], it is defined as a northern deciduous, coniferous forest in the original LPJ model, which is inconsistent with the actual distribution. Hence, in this study, Larix olgensis forests were divided into temperate and cold temperate deciduous coniferous forest, and its bioclimatic limiting parameters were modified. Moreover, the allometric growth parameter k_latosa [32,33] was modified based on [34]; referring to [35], the empirical parameter g_m [36] in the water demand function was modified; referring to [37], the asymptotic maximum mortality rate (k_mort1) and the coefficient (k_mort2) in the death function [38] were modified.

The LPJ model begins from “bare land” (no vegetation cover). The vegetation cover and soil carbon pool reach an approximate equilibrium after 1000 years of “spin up”. This study uses 20 years of climate data from 2000 to 2019 to run this process. Considering that the Larix olgensis forests in the study area are natural secondary forests and plantations, the model is “zeroed” in 1001. Hence, the biomass of four vegetation tissue pools (leaves, sapwood, heartwood, and fine roots) was zeroed, the litter carbon pool and soil carbon pool were simulated as usual, and the NPP of the Larix olgensis forests in Jilin Province from 2000 to 2019 were extracted for analysis.

2.3.2. Input Parameters Setting

Due to several physiological and ecological parameters in the LPJ model, the sensitivity analysis of all parameters is inefficient and unnecessary. Therefore, it is necessary to classify and screen model parameters in advance [17]: (1) some marker parameters need not be optimized, such as vegetation and leaf types; (2) deterministic parameters need not be optimized, such as the leaf-to-root ratio under non-water-stressed conditions (1 by default) and the water scalar value at which leaves shed due to a drought-deciduous PFT (0 by default); and (3) the indeterminable parameters at a specific range, such as fine root carbon-to-nitrogen ratio and fine root turnover, need not be optimized and must use default values. After the screening, 27 adjustable vegetation physical and chemical parameters that do not change with vegetation types were selected from eight aspects: allometry, photosynthesis, water balance, respiration, mortality, soil and litter decomposition, establishment, and parameters varying with Larix olgensis. The 13 main parameters changing with vegetation types were selected for sensitivity analysis, and the change with vegetation mainly refers to the physical and chemical parameters of the Larix olgensis forests different from the other vegetation types in this paper.

Table 1. Definition and intervals of the LPJ model parameters used for the sensitivity analysis.

No	Parameter	Default Value	Lower Value	Upper Value	Description	Reference
1	kallom1-A	100	90	110	parameter in Equation (4)	[37]
2	kallom2-A	40	36	44	parameter in Equation (3)	[39]
3	kallom3-A	0.50	0.45	0.55	parameter in Equation (3)	[39]
4	klatosa-A	4000	3600	4400	parameter (cross-sectional area) in Equation (1)	[34]
5	krp-A	1.60	1.44	1.76	parameter in allometric Equation (4)	[40]
6	θ-P	0.70	0.63	0.77	colimitation (shape) parameter	[41]
7	λmax,C4-P	0.40	0.36	0.44	optimal ratio of intercellular to ambient CO2 concentration in C4 plants	[41]
8	λmax,C3-P	0.80	0.72	0.88	optimal ratio of intercellular to ambient CO2 concentration in C3 plants	[42]
9	α C3-P	0.080	0.072	0.088	intrinsic quantum efficiency of CO2 uptake in C3 plants	[43]
10	α C4-P	0.0530	0.0477	0.0583	intrinsic quantum efficiency of CO2 uptake in C4 plants	[41]
11	a C3-P	0.0150	0.0135	0.0165	leaf respiration as a fraction of Robisco capacity for C3 plants	[43]
12	a C4-P	0.020	0.018	0.022	leaf respiration as a fraction of vmax for C4 plants	[41]
13	αa-P	0.450	0.405	0.495	fraction of PAR assimilated at ecosystem level relative to leaf level	[44]
14	q10ko-P	1.20	1.08	1.32	q10 for temperature-sensitive parameter ko	[42]
15	q10kc-P	2.10	1.89	2.31	q10 for temperature-sensitive parameter kc	[42]
16	q10tau-P	0.570	0.513	0.627	q10 for temperature-sensitive parameter tau	[42]
17	αm-W	1.40	1.26	1.54	Priestley-Taylor coefficient (Demand)	[36]
18	gm-W	3.26	2.934	3.586	empirical parameter in demand function	[35]
19	rgrowth-R	0.25	0.225	0.275	Proportion of growth respiration per unit NPP	[45]
20	kmort1-M	0.050	0.045	0.055	asymptotic maximum mortality rate (/year)	[37]
21	kmort2-M	0.50	0.45	0.55	mortality equation	[37]
22	fair-D	0.70	0.63	0.77	fraction of litter decomposition going directly into the atmosphere	[46]
23	finter-D	0.980	0.882	1.078	fraction of litter entering fast soil decomposition pool	[47]
24	klitter10-D	0.350	0.315	0.385	litter decomposition rate at 10 deg C (/year)	[47]
25	ksoil_fast10-D	0.030	0.027	0.033	fast pool decomposition rate at 10 deg C (/year)	[47]
26	ksoil_slow10-D	0.0010	0.0009	0.0011	slow pool decomposition rate at 10 deg C (/year)	[47]
27	estmax-E	0.120	0.108	0.132	maximum sapling establishment rate	[37]
28	para1-L	0.70	0.63	0.77	fraction of roots in upper soil layer	[10]
29	para4-L	0.30	0.27	0.33	canopy conductance component (gmin, mm/s) not associated with photosynthesis	[42]
30	para7-L	100	90	110	maximum foliar N content (mg/g)	[42]
31	para8-L	0.120	0.108	0.132	fire resistance index	[10]
32	para24-L	−4.0	−3.6	−4.4	low temperature limit for CO2 uptake CO2	[10]
33	para25-L	15.0	13.5	16.5	lower range of temperature optimum for photosynthesis	[10]
34	para26-L	30.0	27.0	33.0	upper range of temperature optimum for photosynthesis	[10]
35	para27-L	42.0	37.8	46.2	high temperature limit for CO2 uptake	[10]
36	para28-L	−13.0	−11.7	−14.3	minimum coldest monthly mean temperature	[10]
37	para29-L	3.0	2.7	3.3	maximum coldest monthly mean temperature	[10]
38	para30-L	900	810	990	minimum growing degree days (at above 5 °C)	[10]
39	para31-L	1000	900	1100	upper limit of temperature of the warmest month	[10]
40	para36-L	0.060	0.054	0.066	interception storage parameter, unitless	[10]

Note: -A represents the allometry module, -P represents the photosynthesis module, -W represents the water balance module, -R represents the respiration module, -D represents the decomposition module, -E represents the establishment module, -M represents the mortality module, -L represents the parameters varying with the Larix olgensis module. The equations marked in the table are from the literature [10].

provides the specific parameters, meanings, and value range required for model simulation. The upper and lower limits of the parameters were set to 10% of the model’s default value; moreover, the parameter distribution was a uniform distribution. The selection and value of the parameters refer to the existing literature.

2.4. Sensitivity Analysis Methods and Scheme

2.4.1. EFAST Method

The EFAST is a sensitivity analysis method based on variance decomposition proposed by Saltelli et al., combining the advantages of FAST and Sobol’s method [14]. It analyzes the direct and indirect impacts of the change in each parameter on the model output and tests the impact of the simultaneous change of multiple parameters on the model results. The frequency curve of the Fourier series is obtained by Fourier transform for sampling, while the variance of the model results caused by the single parameter and parameter interaction is obtained. The main formulae are as follows:

Assuming that there is a model $y=f(x_1,x_2,\cdots ,x_m)$, the dimension is reduced to $y=f(s)$ through an appropriate conversion function, and the Fourier transform can be performed on f(s):

$$x_i=0.5+\frac{\arcsin (\sin (w_{i}s+{\phi }_i))}{\pi }$$

(2)

$$y=f(s)={\displaystyle \sum _{p=-\infty }^{+\infty }[A_p\cos (ps)+B_p\sin (ps)]}$$

(3)

$$A_p=\frac{1}{2\pi }{\displaystyle {\int }_{-\pi }^{\pi }f(s)\cos (ps)ds}$$

(4)

$$B_p=\frac{1}{2\pi }{\displaystyle {\int }_{-\pi }^{\pi }f(s)\sin (ps)ds}$$

(5)

where s is scalar, $s\in [-\pi ,\pi ]$, and ds is the differential of s; $w_i$ is the integer frequency corresponding to $x_i$; ${\phi }_i$ is the random initial phase of the parameter, ${\phi }_i\in [0,2\pi ]$; p is the Fourier transform parameter, $p\in Z$; A_p and B_p are Fourier amplitudes.

The spectrum of the Fourier series is

$${\mathsf{\Lambda }}_p=A_p{}^2+B_p{}^2$$

(6)

The model output variance $V_i$ caused by the change of the parameter $x_i$ is

$$V_i=V[E(y|x_i)]={\displaystyle \sum _p{\mathsf{\Lambda }}_p{w}_i}=2{\displaystyle \sum _{p=1}^{+\infty }{\mathsf{\Lambda }}_p}{w}_i$$

(7)

where $E(y|x_i)$ is the expectation of y corresponding to the parameter $x_i$. The total variance V(y) of the model is

$$V(y)={\displaystyle \sum _{p{\in }^Z}{\mathsf{\Lambda }}_p}=2{\displaystyle \sum _{p=1}^{+\infty }{\mathsf{\Lambda }}_p}$$

(8)

The scalar s is sampled at $[-\pi ,\pi ]$, and the sampled value is converted into the value of each parameter through conversion to obtain the Fourier amplitudes A_p and B_p approximately:

$$A_p\approx \frac{1}{N_s}{\displaystyle \sum _{k=1}^{N_s}f(s_k)\cos (p^{\prime }{s}_k)}$$

(9)

$$B_p\approx \frac{1}{N_s}{\displaystyle \sum _{k=1}^{N_s}f(s_k)\sin (p^{\prime }{s}_k)}$$

(10)

where $N_s$ is the sample size, $s_k$ is the kth sampling value of s, and $p^{\prime }$ is the variation parameter, $p^{\prime }\in \{-\frac{N_s-1}{2},\cdots ,-1,0,1,\cdots ,\frac{N_s-1}{2}\}$. The total variance $V(y)$ of all parameters of the model can be decomposed into

$$V(y)={\displaystyle \sum _{i=1}^k{V}_i}+{\displaystyle \sum _{1\le i\le j\le k}{V}_{ij}+\cdots +V_{1\cdots k}}$$

(11)

where y is the output result (i.e., NPP), k is the number of parameters, V_i is the variance of the ith parameter, V_ij is the variance of the output result caused by the interaction of the ith and jth parameters, and V_1⋯k is the variance of the output result caused by the interaction of all parameters.

Then, the first-order sensitivity index (S_i) and the total order sensitivity index (ST_i) were calculated according to the variance of each parameter. The formulae are the following:

$$S_i=\frac{V_i}{V(y)}$$

(12)

$$S{T}_i=\frac{V(y)-V_{-i}}{V(y)}$$

(13)

where $V_{-i}$ is the variance of all parameters except the ith parameter.

The greater the sensitivity index of the parameter, the more significant the direct or indirect impact (or contribution) of the parameters on the model results, and vice versa. Generally, parameters with a first-order sensitivity index greater than 0.05 and a global sensitivity index greater than 0.10 can be deemed sensitive parameters to the model output results [16,26,48].

2.4.2. Morris Method

The Morris method or the element effects (EE) method was proposed by Morris [13] in 1991, and then improved by Campolongo by optimizing the sampling strategy, whose calculation cost is significantly reduced [49]. This method improves the shortcomings of the local sensitivity analysis method. It originated from the experimental design of one change. Only one parameter is changed as the different situations reset for each parameter, all situations of the parameters are sampled, and the sensitivity of the parameters is calculated by differential methods. It is suitable for model sensitivity analysis with several input factors or high calculation costs. It can render better consideration to efficiency and accuracy and distinguish sensitive parameters and their direct interaction under fewer sampling conditions. The formula is as follows:

$$d(x_1,\cdots ,x_k,\mathsf{\Delta })=[y(x_1,\cdots ,x_{i-1},x_i+\mathsf{\Delta },x_{i+1},\cdots ,x_k)-y(x_1,\cdots ,x_k)]/\mathsf{\Delta }$$

(14)

where d represents the unit effect—that is, the change rate of the output results of two groups of samples with one different parameter; y(x) is the model output, x₁,⋯, x_k is the k parameter variables, $\mathsf{\Delta }$ is the parameter change value with the value of 1/(p − 1), and p is the horizontal quantity, equivalent to the quartile.

Then, we calculate the mean value of the unit effect (d) of each parameter (μ*) and standard deviation (σ). A more significant μ* indicates that the overall sensitivity of the input parameters is significant, and a more significant σ implies that the interaction or nonlinear effect of the input parameters is strong. Considerably, this parameter is sensitive to the model output only when μ* > μ_mean (μ_mean is the mean of each parameter μ*).

2.4.3. Sensitivity Analysis Scheme

The sensitivity analysis of the model parameters is realized through the sensitivity and uncertainty analysis software SIMLAB (version 2.2), and the program is written in the Python environment to realize the batch operation of the LPJ model parameters. The specific scheme is as follows. Firstly, the SIMLAB software sets the number of parameters, parameter distribution, global sensitivity analysis method, and sampling times. A total of 3880 (n·k) samples are taken by EFAST, where k is the number of parameters and n is the number of repetitions of the EFAST method. It is specified that n ≥ 65 and n in this study is 97. A total of 410 ((k + 1) × r) samples are taken by the Morris method, and k is the number of parameters and r is the number of repetitions of the Morris method, where r_min = 4, r = 10 in this study). Secondly, the automatic call of the LPJ model parameters and the output of simulation results are realized in the Python environment. Finally, the global sensitivity analysis of the output results is conducted using SIMLAB software. In addition, during the LPJ model, to improve the batch-processing efficiency, 13 parameters with varying physical and chemical parameters of Larix olgensis forests are read separately in a conditional cycle, and other types are simulated by default. Meanwhile, this paper applies two stages, 2009–2014 and 2000–2019, for sensitivity analysis to avoid the contingency of the output results of a particular stage in parameter sensitivity analysis. It uses the Spearman rank correlation method to analyze the NPP parameters of the sensitivity analysis output results of the EFAST and Morris method.

2.5. Remote Sensing Observation Data of NPP

NPP remote sensing data were downloaded from the MOD17A3H dataset from 2000 to 2019 (https://e4ftl01.cr.usgs.gov/ (accessed on 1 February 2020)) with a spatial resolution of 500 m × 500 m, and the unit is kg C/m² (hereinafter uniformly converted to g C/m² and called the MODIS NPP). The data were mainly used to verify the simulation results of the LPJ model.

2.6. Measured Data and Accuracy Evaluation Index

The measured data mainly include the measured data of the fixed and temporary sample plots of Larix olgensis forests. The fixed sample plot data were acquired from the eighth and ninth national forest resources continuous inventory data of Jilin Province (hereinafter referred to as “FCI” data), and the survey times were 2009 and 2014, distributed in 213 sample plots of young forest (66), middle-aged forest (63), near mature forest (39), and mature forest (45). According to the biomass equation of the main tree species in Northeast China studied by Dong et al. [50] and Wang et al. [51], the sample wood biomass in the sample plot was calculated, summarized, and divided by the sample plot area to get the biomass per unit area of the sampled plot. The annual mean value of NPP was calculated by multiplying the difference in biomass per unit area before and after the sample plot by carbon content rate (uniformly 0.5 [52]), and then divided by interval time (from 2009 to 2014). These data mainly aim to verify the model optimization results. The temporary sample plot data were from 880 Larix olgensis forests in 19 sample plots obtained from the field survey in 2019. The equations for diameter at breast height (D)-tree height, and D-crown area of the Larix olgensis forests were established, primarily used to modify the sensitive parameters screened by the LPJ model.

The mean relative error (MRE) and the mean absolute error (MAE) verify the accuracy of the annual mean NPP simulated by the LPJ model (Equations (15) and (16)). The smaller the MRE and MAE, the higher the accuracy.

$$\mathrm{MRE}=\frac{1}{n}{\displaystyle \sum _{i=1}^n\frac{\left|{\mathrm{NPP}}_{\mathrm{obs}}-{\mathrm{NPP}}_{\mathrm{LPJ}}\right|}{{\mathrm{NPP}}_{\mathrm{obs}}}}$$

(15)

$$\mathrm{MAE}=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|{\mathrm{NPP}}_{\mathrm{obs}}-{\mathrm{NPP}}_{\mathrm{LPJ}}\right|}$$

(16)

where n is the total number of years, i is the ith year, ${\mathrm{NPP}}_{\mathrm{LPJ}}$ is the LPJ model-simulated NPP, and ${\mathrm{NPP}}_{\mathrm{obs}}$ is the measured value (2009–2014 signifies the measured NPP, while 2010–2019 provides the MODIS NPP).

3. Results

3.1. LPJ Model Output NPP Based on the EFAST and Morris Methods

Based on the EFAST and Morris method, NPP was simulated after parameter sampling in 2009–2014 and 2000–2019, obtaining the simulation results of 3880 and 410 groups of parameter values, respectively (Figure 3). The simulation results of both methods fluctuated around their mean values in each period. The annual mean values of the EFAST method were 598 g C m⁻² a⁻¹ and 546 g C m⁻² a⁻¹, while the annual mean values of the Morris method were 617 g C m⁻² a⁻¹ and 566 g C m⁻² a⁻¹. Thus, the simulation value based on the Morris method was slightly higher than that based on the EFAST method. The simulated annual mean value of the two methods in 2000–2019 was greater than that of the MODIS NPP data (572 g C m⁻² a⁻¹), while the mean in 2009–2014 was less than that of the MODIS NPP data (573 g C m⁻² a⁻¹). The simulated values of the two methods in 2000–2019 were slightly larger than those in 2009–2014, primarily because the simulated values in 2012 were significantly lower than the MODIS NPP. It is probably caused by the lower temperature than the other years and a continuous cloudy climate.

3.2. Importance Ranking of the LPJ Model to NPP Sensitivity Based on the EFAST Method

Based on the EFAST method, sensitivity analysis was conducted between the mean NPP and the input parameters of the model in 2009–2014 and 2000–2019 (Figure 4), depicting that the parameter sensitivity was very similar in both stages. The global sensitivity index and the first-order sensitivity index of parameter k_rp in the two periods were the largest, 0.709 and 0.666 in 2009–2014 and 0.711 and 0.665 in 2000–2019. It could explain more than 70% and 60% of the variation variance, which is significantly higher than the sensitivity of the other parameters. Meanwhile, it depicted that the direct contribution rate of this parameter on the model results and the coupling effect with the other parameters were the largest. The sensitivity index of α_m in 2009–2014 and 2000–2019 was second only to k_rp, while the global sensitivity index and the first-order sensitivity index were greater than 0.1. The global sensitivity indices of α_a and g_m from 2009 to 2014 were 0.104 and 0.100, while the first-order sensitivity indices were 0.085 and 0.078.

Meanwhile, the global sensitivity indices of the parameters from 2000 to 2019 were 0.121 and 0.100, and the first-order sensitivity indices were 0.105 and 0.076. These four parameters could be identified as the sensitivity parameters, exhibiting a strong coupling effect. The first-order sensitivity index was greater than 0.05, while the global sensitivity index was greater than 0.10. Additionally, the sensitivity ranking of NPP simulated by the LPJ model parameters based on EFAST is as shown in

Table 2. Sensitivity analysis rankings of the LPJ model parameters for the NPP output response based on the EFAST and Morris method.

NO	Parameters	2009–2014		2000–2019
		EFAST Rank	Morris Rank	EFAST Rank	Morris Rank
1	kallom1-A	7	7	7	7
2	kallom2-A	9	9	10	9
3	kallom3-A	5	6	5	6
4	klatosa-A	8	12	8	12
5	krp-A	1	1	1	1
6	θ-P	12	23	12	13
7	λmax,C4-P	35	31	34	37
8	λmax,C3-P	32	25	33	26
9	α C3-P	6	3	6	2
10	α C4-P	22	36	24	40
11	a C3-P	30	13	25	15
12	a C4-P	37	20	38	27
13	αa-P	3	4	3	3
14	αm-W	2	2	2	4
15	gm-W	4	5	4	5
16	kmort1-M	29	14	26	14
17	kmort2-M	31	37	30	38
18	fair-D	27	21	29	17
19	finter-D	13	15	13	16
20	estmax-E	23	11	23	11
21	rgrowth-R	15	8	11	8
22	q10ko-P	16	32	17	32
23	q10kc-P	20	24	18	20
24	q10tau-P	17	10	9	10
25	klitter10-D	11	38	22	36
26	ksoil_fast10-D	19	33	19	31
27	ksoil_slow10-D	21	39	20	35
28	para1-L	40	17	40	28
29	para4-L	33	40	35	34
30	para7-L	26	19	28	29
31	para8-L	39	22	39	23
32	para24-L	10	28	15	24
33	para25-L	25	18	27	18
34	para26-L	28	35	31	39
35	para27-L	18	27	14	33
36	para28-L	38	16	37	19
37	para29-L	36	26	36	22
38	para30-L	34	34	32	30
39	para31-L	24	30	21	21
40	para36-L	14	29	16	25

; k_allom3, α_C3, k_allom1, k_latosa, and k_allom2 ranked high in sensitivity, indicating that they also largely contributed to the model results. The sensitivity indices of the other parameters were small and had little difference, implying that the coupled correlation between these parameters was small and relatively independent of each other, and that the default values could be applicable.

To conclude, the parameters that significantly impacted the NPP model output by EFAST generally existed in the photosynthesis module, water balance module, and allometric growth module. Among them, the parameters k_rp, k_allom1, k_allom3, and k_allom2 were mainly related to the allometric growth of individual vegetation, acting on the D-tree height equation and the D-crown equation, respectively, which directly reflect the growth characteristics and were the basis of NPP. Parameter α_C3 was the inherent quantum rate of CO₂ absorption by C₃ plants, acting on photosynthesis. It was one of the factors controlling the initial slope of the photosynthesis equation and represented the maximum efficiency of vegetation absorption and the conversion of incident light energy. Parameter α_a denoted the proportion of photosynthetic effective radiation assimilated and absorbed at the ecosystem level related to leaves, mainly applied to calculate the absorption of photosynthetic effective radiation by vegetation. Those were the two parameters with strong sensitivity in the photosynthesis module that acted on an essential link in photosynthesis and regulated the efficiency of vegetation in converting photosynthetic effective radiation into NPP. The parameters α_m and g_m were essential parameters controlling water balance in the water demand equation, and were directly involved in vegetation evapotranspiration, related to the daily net photosynthetic value, and were sensitive to NPP.

3.3. Importance Ranking of the LPJ Model to NPP Sensitivity Based on the Morris Method

Sensitivity analysis was conducted between the output mean NPP and the input parameters of the model in 2009–2014 and 2000–2019 based on the Morris method (Figure 5). The value of the sensitivity index of parameter k_rp in the two periods (μ* and σ) was large and not shown. The parameter sensitivity of this method was also highly similar in 2009–2014 and 2000–2019, which could exclude the varying effects of the output results of the LPJ model in a particular period on parameter sensitivity. Consistent with the sensitivity analysis results based on EFAST, the most sensitive parameter was k_rp, whose μ* and σ were 148.2 and 135.3 in 2009–2014 and were 177.2 and 145.6 in 2000–2019, respectively, far exceeding the mean (i.e., μ_mean, 16.4 and 17.5, respectively), which was significantly higher than the sensitivity of the other parameters. Additionally, σ was significantly greater than σ_mean (the mean values of σ were 10.7 and 11.0, respectively), indicating that the parameter k_rp had the most remarkable overall sensitivity and interaction to the model results. The overall sensitivity of α_m, α_C3, and α_a was second only to k_rp, whose μ* was more than 55, and σ was greater than the σ_mean, signifying that the three parameters also interacted with the other parameters. In contrast, α_a with a smaller σ had less interaction with the other parameters. The μ* > μ_mean of g_m, k_allom3 plus k_allom1, with larger σ, could be considered sensitive parameters with a strong coupling effect. In addition, the sensitivity ranking of the k_allom2 and r_growth parameters was also higher (Table 2), which impacted the model output results. Based on the Morris method, the parameter sensitivity, which had a significant impact on the NPP of the output results of the model in the two periods, was the same, and it also existed in the photosynthesis, water balance, and allometric growth modules, which were the same as the EFAST results.

3.4. Comparison of Analysis Results Based on the EFAST and Morris Methods

The sensitivity analysis results of the LPJ model parameters based on the EFAST method and the Morris method were the same. For NPP of Larix olgensis forests in 2009–2014 and 2000–2019, the parameter with the most extensive sensitivity index was k_rp. The parameters α_m, α_a, and g_m were likewise sensitive parameters, showing that the two methods could screen the sensitive parameters relative to the output results. Thus, it was preferred to calibrate the methods in future parameter localization research. Accordingly, the sensitive parameters obtained by the Morris method were slightly more than the EFAST method, mainly including α_C3, k_allom3, k_allom1, and k_allom2, which were likewise very high in sensitivity rankings based on the EFAST method. Among the sensitivity parameters based on the Morris method, α_C3 was one of the sensitivity parameters with a high σ and μ*, representing the proportion of leaf respiration in the carboxylation capacity of C₃ vegetation and acting on the photosynthesis module, and consequently having a particular impact on NPP.

The Spearman rank correlation coefficient applied to the correlation coefficient between the ordered variables and was highly suitable for correlation analysis of the sensitivity rankings of the Morris method and EFAST. The correlation analysis of the NPP parameters of the sensitivity analysis output results of the EFAST (ST_i) and the Morris method (μ*) using the Spearman rank correlation method indicated that the rank correlations of the two methods were very significant or significant (

Table 3. Correlation coefficients between the Morris method (μ*) and EFAST (ST_i).

Parameters of Analysis	2009–2014	2000–2019
All parameters	0.501 **	0.674 **
Both sensitively	0.883 **	0.738 *

Note: ** indicates the correlation is significant at the level of 0.01 and * is at the 0.05 level.

). The correlation coefficients of all parameters were 0.501 and 0.674, respectively, and the correlation analysis of the parameters sensitive to both methods showed a very significant correlation (0.833) in 2009–2014 and a significant correlation (0.738) in 2000–2019.

Based on the advantages of EFAST, the first-order and interactive sensitivity index contribution ratio of NPP in 2009–2014 and 2000–2019 to the parameters jointly sensitive to the two methods were examined (Figure 6). The results showed that the contributions of the screened sensitivity parameters to NPP in the two periods were the same, mainly coming from the direct contribution rates of each parameter, which were between 78–93% and 78–94%, respectively. Conversely, the interaction was small, and the greater the sensitivity index, the greater the direct contribution rate to NPP. The sums of the first-order sensitivity index contribution rates of the four common sensitivity indices in the two periods were consistent, at 89%. Thus, the four parameters could explain 89% of the NPP variation variance, reflecting its significance to model optimization.

3.5. Model Optimization and Accuracy Verification

Since EFAST could quantitatively calculate the contribution rate of each parameter to the variance of the model results, we selected a group of parameters close to the annual MODIS NPP from 2010 to 2019 originating from the 3880 combinations generated in the sensitivity analysis. The evaluation function was the cumulative mean error (CME, Equation (17)) between the simulated value and the observed value. Utilizing the measured data of the temporary sample plots, k_rp, k_allom3, k_allom1, and k_allom2, among the more sensitive parameters were optimized by the power function (y = ax^b), or Equations (3) and (4) in Sitch et al. [10]. The results were k_rp = 1.31, k_allom1 = 82.35, k_allom3 = 0.48, and k_allom2 = 42.879 (Figure 7). The LPJ model was rerun to compute the annual mean NPP of Larix olgensis forests in Jilin Province from 2009 to 2014 and from 2010 to 2019 with the other parameters using the default values. It was compared with the measured data of the fixed sample plots and the MODIS NPP in the same stage (

Table 4. Table of precision evaluation for annual mean NPP.

Periods	Methods	Before Optimization		After Optimization
		MRE (%)	MAE (g C m−2 a−1)	MRE (%)	MAE (g C m−2 a−1)
2009–2014	LPJ–Measured	50.8	193.2	40.2	153.0
2010–2019	LPJ–MODIS NPP	11.4	65.5	9.8	56.7

Note: MRE—mean relative error; MAE—mean absolute error.

).

$$\mathrm{CME}=({\displaystyle \sum _{i=1}^n\left|{\mathrm{NPP}}_{\mathrm{LPJ}}{-\mathrm{NPP}}_{\mathrm{obs}}\right|})/n$$

(17)

The NPP simulation of the LPJ model before and after optimization were compared and analyzed combined with the MODIS NPP and the measured data (Table 4). The results indicated that the MRE and MAE of the LPJ-simulated and MODIS NPP significantly decreased during 2009–2014 and 2010–2019. Compared with the measured NPP, the optimized simulated MRE and MAE decreased by 10.6% and 40.2 g C m⁻² a⁻¹. The MRE and MAE between the optimized simulation values and the MODIS NPP further decreased, indicating that the optimization effect of the above parameters was relatively good. Similarly, it can reasonably be extended to other areas or other species for simulation according to local climate and species characteristics. In addition, compared to the measured values of NPP, the simulated values of the LPJ model were close to the MODIS NPP. The MRE and MAE were smaller than the measured data, determining that the measured values of NPP from 2009 to 2014 were different from the simulated values of the LPJ model; this was because the measured values did not include the productivity of shrubs and grasses, and that human activities, such as new planting and logging, were not considered in the simulation. Based on the “FCI” data, there was abnormal cutting of sample trees in some sample plots, most of which were Larix olgensis forests, even reaching more than 50% of the total number of standing trees in the sample plot. This factor was not considered in the model.

On the other hand, the measured value was the mean value of various places from 2009 to 2014, not the measured value of that year. Following parameter optimization, the mean annual NPP value of Larix olgensis forests in Jilin Province from 2010 to 2019, as simulated by the LPJ model, was 580 g C m⁻² a⁻¹, with an annual mean growth rate of 2.13%, slightly higher than that of the MODIS NPP (2.02%). Still, there was no significant difference (p > 0.05), exhibiting a fluctuating growth trend (Figure 8).

4. Discussion

4.1. Comparison between the EFAST and Morris Method

Based on the EFAST and Morris method, the parameter sensitivity analysis of the annual mean NPP of Larix olgensis forests in Jilin Province simulated by the LPJ model in 2009–2014 and 2010–2019 showed that the sensitivity parameters screened by the two methods were highly correlated. The parameter k_rp was the most sensitive parameter, followed by parameters α_m, α_a, g_m, and α_C3. The parameters k_allom3, k_allom1, and k_allom2 had a higher sensitivity ranking, and were mainly distributed in the photosynthesis, water balance, and allometric growth modules. They were indispensable parameters for controlling the positive succession of the entire ecosystem, as NPP could best reflect the production capacity of the vegetation ecosystem in a natural state. Thus, it was reasonable to screen these sensitive parameters by both methods. Pappas’s study of LPJ-GUESS, based on the Morris method, indicated that the photosynthetic parameters were susceptible, which impacted the parameters of forest structure (growth, establishment, and mortality) [27]. This was also consistent with the results analysis of sensitivity when simulating carbon and water flux in China with LPJ-DGVM model based on the EFAST method, whose results showed that α_C3 was the most important parameter of the carbon flux module for the output results GPP, NPP, and R_h, representing the utilization efficiency of light [25]. In this paper, the sensitivity ranking of α_C3 was also very high. In addition, Ma showed that α_a and g_m were also sensitive in the study about the LPJ model [53]. Zhang et al. also showed that the parameters of photosynthesis module and water balance module needed to be optimized first when using a simulated annealing algorithm to optimize the parameters of the Biome-BGC model [17].

Overall, the methods used in this study had roughly the same results in the NPP parameter sensitivity analysis, which could screen the parameters that had a significant impact on the target variables, and the regularity was relatively similar. Li et al. [26] used the EFAST and Morris method to analyze the sensitivity of CROPGRO-Tomato model under different irrigation treatments, which showed the most sensitive parameters both could be screened, and the results of both methods were highly correlated for the output variables. Wu et al. [48] and Song et al. [54] analyzed the sensitivity of the CROPGRO-Cotton model and CERES-Wheat model, respectively, based on both methods, which also showed that there was a significant correlation between the parameter sensitivity results obtained by both methods. However, both methods had their strengths and weaknesses. The Morris method required fewer running times and a higher model efficiency, significantly reducing the working time and more operability in model localization. The most significant feature of EFAST was to sample the parameters through the search function with a high accuracy. It could quantitatively analyze the first-order and global sensitivity indices of each parameter and calculate quantitatively the contribution rate of each parameter to the variance of the model results. Likewise, it required a large number of samples. Vazquez-Cruz et al. [28] explored the sensitivity analysis of the TOMGRO model based on the EFAST and Sobel methods, which also showed the superiority of the EFAST method.

In this paper, a global sensitivity analysis of the NPP parameters was carried out for the output results of the LPJ model. The results of the two sensitivity analysis methods show that the parameters of the photosynthesis module have a great impact on the output of NPP. Referring to the research of Wu et al. [48], it was found that the number of sensitive parameters of the two methods for different output results are different. In the sensitivity analysis of the CROPGRO model, it was found that for the maximum leaf area of the output results, the number of sensitive parameters screened by the Morris method is more than that by the EFAST method, but on the contrary for aboveground dry matter quality. In view of the different levels of correlation between the two methods with different output results, the choice of sensitivity analysis method of the model output results should be fully considered in the future optimization of the LPJ model. This study is based on the fact that the sampling times of the EFAST is nearly 10 times that of the Morris method. Especially for the complexity of the model process, the more times the sensitivity analysis is, the lower the operation efficiency of the model is.

Therefore, if the number of model parameters was enormous and only the overall sensitivity analysis results were required, the Morris method could fully meet the requirements. If a more detailed analysis was necessary, it was recommended to select the sensitivity analysis results of the EFAST with more scientific sampling and more operations as the research analysis [26]. The most sensitive parameters screened by the two methods were the same, and the specific one depended on the situation.

4.2. Uncertainty of Sensitivity Analysis Method

Various physiological and ecological parameters have caused uncertainty that the cumulative distribution function could quantify in the model simulation. Taking the sensitivity analysis based on the EFAST as an example, the annual mean NPP in 2009–2014 and 2000–2019 showed positive skewness (Figure 9); the peaks were all left, with a distribution range of 399.2–819.6 g C m⁻² a⁻¹ and 439.6–877.2 g C m⁻² a⁻¹, respectively, and the highest frequencies were 520–540 g C m⁻² a⁻¹ and 580–600 g C m⁻² a⁻¹, with frequencies of 562 and 518, respectively. The median was within this interval (

Table 5. Uncertainty analysis statistics for the different periods as determined from the EFAST global sensitivity analysis.

Period of NPP	The MODIS NPP	Medium Value	95% Confidence Interval Lower Limit	95% Confidence Interval Upper Limit	Standard Deviation
2009–2014	573	536	418	675	66
2000–2019	580	587	463	734	69

), and the 95% confidence intervals of the NPP in the two periods were 418–675 g C m⁻² a⁻¹ and 463–734 g C m⁻² a⁻¹, respectively. The simulation results were mainly concentrated on 480–560 g C m⁻² a⁻¹ and 540–620 g C m⁻² a⁻¹, which were close to the MODIS NPP in the same period but were relatively large output values in both periods.

Referring to the settings of other models, such as the CERES model [54] and CROPGRO model [22,26,48], this study set the upper and lower limits of the LPJ model parameters to 10% of the model’s default values. The sensitivity of the selected parameters changing with PFT type was weak; in particular, several climate parameters did not show the difference, such as the temperature limit of CO₂ absorption, optimal temperature, minimum temperature, and maximum temperature limit of photosynthesis. The most likely reason for this was that the setting range was too small, resulting in the sensitivity of the parameters not being reflected. However, Zhang et al. set the default value to float up and down by 20–50% when optimizing some physiological and ecological parameters of the BIOME model, after which the R² of the linear fitting equation between the simulated value of the optimized model and the flux observation value was improved [17]. Maybe need to adjust the parameters’ range in sensitivity analysis in next study.

Meanwhile, Shin et al. pointed out that the parameter value range and distribution type directly affect the sensitivity ranking of the parameters [55]. Vanuytrecht et al. indicated that it is necessary to explore the value and distribution patterns of some parameters without specific physical or biological significance from their possible minima to maxima during correction when they analyzed the parameter sensitivity of the AquaCrop model [56]. The distribution patterns of the parameters involved in sensitivity analysis will also affect the analysis results. This study assumes a uniform distribution and due to the lack of knowledge, some biological significance may be ignored. Based only on the previous results, there may be limited parameter range and distribution settings. Therefore, we can expand the value ranges of some parameters, set varying parameter change ranges, and analyze the impact of parameter change ranges on the output variables sensitivity, to determine the best floating ratio in the next step.

4.3. Uncertainty of the LPJ Model

The vegetation distribution simulated by the optimized LPJ model in Jilin Province was consistent with the actual distribution. The Larix olgensis forests were mainly distributed in the eastern mountainous areas and the central low mountain and hilly area of Jilin Province, with a similarity to the measured value of NPP, indicating that the model simulation effect was good. It was close to and had no significant difference from the MODIS NPP, but slightly lower than the simulation value (726 g C m⁻² a⁻¹) of the net primary productivity of the Changbai Larix olgensis forests in Northeast China by He [57] using the Biome-BGC model from 1998 to 2010. At the same time, the research indicated that the simulation value of the Biome-BGC model was higher than the measured NPP of its sample plot, and the research scope of the model’s different parameter settings and data sources further affected the simulation results. Additionally, compared to [58], using the LPJ model to examine the net primary productivity of natural vegetation in China, the NPP in Changbai Mountain was 600–700 g C m⁻² a⁻¹, which was proximate to the results of this study.

Although the LPJ model, after parameter optimization, estimated the NPP of Larix olgensis forests in Jilin Province, which was similar to the estimation and remote sensing inversion results of the other models, this model’s simulation results were reasonable in the productivity of Larix olgensis forests in Jilin Province. There were still differences and uncertainties in the simulation results between different models due to data limitations.

The possible uncertainties regarding the results of this study are as follows. Firstly, there may be some errors in the screening of the measured data used for verification. In addition, the biomass calculation method used here only calculate the aboveground biomass, and there is a lack of litter biomass. According to Wang et al., the litter biomass of Larix olgensis accounted for 4.255% [59]. So, the NPP simulation from the plot was underestimated. Secondly, some physiological and ecological parameters in the optimized LPJ model were the selected optimal combination based on referring to the settings in previous studies; this has had an impact on the simulation results. Thirdly, the model did not consider the impact of human activities on NPP simulation, including the positive role of ecological construction in natural forest protection projects on NPP growth, natural disasters, and the decline caused by abnormal logging. Fourthly, the resolution of the LPJ model simulation results was 0.1° × 0.1°, while the resolution of the MODIS NPP was 500 m × 500 m and the sample plot size was 0.06 hm², implying that the conversion between different scale data brought some differences to the simulation results. Fifthly, although the model introduced the fire interference module, it could not fully explain the mechanism of triggering a fire and proposed various theoretical assumptions, which was a problem for most models. Additionally, the model version used in this paper had not considered the nitrogen cycle process; it was well known that the significance of the nitrogen cycle to forest growth would affect the estimation accuracy of NPP. Sixthly, the model set 10 vegetation function types, ignored other more detailed types such as farmland, and assumed that the same vegetation function type used the same set of physiological and ecological parameters. It did not consider the differences among specific tree species and changes over time, ignoring vegetation functional types’ influence on the geographical environment, site conditions, age, climate, and other factors. Undoubtedly, this was also a difficulty in model optimization. The improvement of the model is to better simulate the dynamics of NPP or biomass change in the future. It is more meaningful to carry out sensitivity analysis according to different climatic region conditions. We find that the sensitivity parameters of the LPJ model simulation results under different climate stresses are different in the research of Pappas et al. [27], so different climate changes can be added to future research, such as temperature, precipitation, CO₂ concentration, and so on.

The LPJ model is an ecological process model with several physiological and ecological parameters. Although the model results are affected by many aspects, it is found that the MRE or MAE were markedly reduced compared with the observation data or the measured data in this study; this shows that the optimization direction of this paper was feasible. We especially improved the model performance by optimizing the parameters of the tree species. Next, when estimating the NPP of vegetation on a smaller scale, it is vital to consider the landform, the climate attributes of the study area, and the improved model of the physiological and ecological attributes of the study species. Then, one should examine the applicability and optimization scheme of the model under extreme climate conditions. In the future, we can determine the law of parameter differences in different vegetation types to enable a model with better local applicability and a higher accuracy, which is also one of the current research directions.

5. Conclusions

The results of the two methods were similar but not identical. The sensitivity parameters were significantly correlated (p < 0.05), and parameter k_rp was the most sensitive parameter, followed by the parameters α_m, α_a, and g_m. The sensitive parameters screened by the Morris method were more than by the EFAST, which were α_C3, k_allom3, k_allom1, and k_allom2, respectively, with a very high sensitivity ranking in the EFAST. These sensitive parameters were mainly found in the photosynthesis, water balance, and allometric growth modules. The EFAST method had a high precision, which could quantitatively analyze the first-order and global sensitivity indices of each parameter, and quantitatively calculated the contribution rate of each parameter to the variance of the model results. Still, it required a large number of samples. The Morris method necessitated fewer model running times and higher efficiency. Overall, the most sensitive parameters screened by the two methods were the same, and the specific one depended on the situation.

The MRE and MAE between the Larix olgensis forests NPP simulated by the LPJ model with optimized parameters and the validation data markedly decreased in Jilin Province from 2010 to 2019. There is no doubt that this model optimization is effective. Meanwhile, the optimized annual mean value of NPP from 2010 to 2019 was 580 g C m⁻² a⁻¹, with an annual mean growth rate of 2.13%, exhibiting a fluctuating growth trend, slightly higher than that of the MODIS NPP (2.02%).