Open AccessArticle

Parameter Optimization of the 3PG Model Based on Sensitivity Analysis and a Bayesian Method

Chenjian Liu

^1,2,3,

Xiaoman Zheng

and

Yin Ren

^1,*

Key Laboratory of Urban Environment and Health, Fujian Key Laboratory of Watershed Ecology, Institute of Urban Environment, Chinese Academy of Sciences, Xiamen 361021, China

College of Resources and Environment, Fujian Agriculture and Forestry University, Fuzhou 350002, China

School of Public Administration, Fujian Agriculture and Forestry University, Fuzhou 350002, China

Author to whom correspondence should be addressed.

Forests 2020, 11(12), 1369; https://doi.org/10.3390/f11121369

Submission received: 30 October 2020 / Revised: 17 December 2020 / Accepted: 17 December 2020 / Published: 21 December 2020

(This article belongs to the Section Forest Inventory, Modeling and Remote Sensing)

Download

Browse Figures

Figure 1
Map of the study area. "> Figure 2
Flowchart of the methodology followed in this study. 3PG: Physiological Principles in Predicting Growth model; MCMC: Markov chain Monte Carlo method; stemNo: stand stocking; WF: foliage biomass; WR: root biomass; WS: stem biomass including branches and bark; LAI: canopy leaf area index (LAI); standvol: stand volume excluding branches and bark; MAI: mean annual volume increment; avDBH: stand-based mean diameter at breast height (DBH). "> Figure 3
The effects of the adjustment of model parameters on their sensitivity to stand stocking (stemNo) and diameter at breast height (DBH). "> Figure 4
Sensitivity grades of model parameters to stand stocking and DBH. "> Figure 5
The posterior distributions of seven parameters. Note: The x-axis represents the value of the parameter and the y-axis represents the probability density. "> Figure 6
Comparison of stand biomass before and after adjustment. Note: the red dots represent the observational data of Cunninghamia lanceolata trees from the Shunchang Forest Farm, the thin line represents the predicted stand biomass value (according to the observational data, the 3PG model can automatically generate a prediction curve), and the thick line represents the simulated stand biomass value. ">

Versions Notes

Abstract

Sensitivity analysis and parameter optimization of stand models can improve their efficiency and accuracy, and increase their applicability. In this study, the sensitivity analysis, screening, and optimization of 63 model parameters of the Physiological Principles in Predicting Growth (3PG) model were performed by combining a sensitivity analysis method and the Markov chain Monte Carlo (MCMC) method of Bayesian posterior estimation theory. Additionally, a nine-year observational dataset of Chinese fir trees felled in the Shunchang Forest Farm, Nanping, was used to analyze, screen, and optimize the 63 model parameters of the 3PG model. The results showed the following: (1) The parameters that are most sensitive to stand stocking and diameter at breast height (DBH) are nWs(power in stem mass vs. diameter relationship), aWs(constant in stem mass vs. diameter relationship), alphaCx(maximum canopy quantum efficiency), k(extinction coefficient for PAR absorption by canopy), pRx(maximum fraction of NPP to roots), pRn(minimum fraction of NPP to roots), and CoeffCond(defines stomatal response to VPD); (2) MCMC can be used to optimize the parameters of the 3PG model, in which the posterior probability distributions of nWs, aWs, alphaCx, pRx, pRn, and CoeffCond conform to approximately normal or skewed distributions, and the peak value is prominent; and (3) compared with the accuracy before sensitivity analysis and a Bayesian method, the biomass simulation accuracy of the stand model was increased by 13.92%, and all indicators show that the accuracy of the improved model is superior. This method can be used to calibrate the parameters and analyze the uncertainty of multi-parameter complex stand growth models, which are important for the improvement of parameter estimation and simulation accuracy.

Keywords:

Chinese fir; Markov chain Monte Carlo (MCMC); parameter estimation; stand models

1. Introduction

Experimental observations and model simulations are common methods for estimating the pattern and variability of the global carbon cycle in order to understand the key processes and control mechanisms of this cycle. To reduce the uncertainties of the parameters of ecological models and improve the ability of such models to simulate and predict ecosystem processes and changes, researchers have successively carried out a series of parameter estimation studies for terrestrial ecosystem models focusing on sample plots at regional and global scales [1,2,3,4,5,6].

As an important part of the carbon pool, forests play a pivotal role in the carbon cycle of terrestrial ecosystems and the global carbon cycle. Determining how to obtain the key parameters of forest growth and estimate forest biomass have become hot research topics. Achieving accurate model simulations depends on the accurate acquisition of many model parameters, weather-driven data, and site parameters. Among these, model parameters are the main source of uncertainty in the simulation results. Therefore, the accurate acquisition of model parameters is a prerequisite for model application and the improvement of model prediction accuracy [7,8].

It is a parameter estimation problem to obtain the parameters of a model from observed values. For linear equations and simple nonlinear equations, the least-squares method can be used to estimate the parameters. The Physiological Principles in Predicting Growth (3PG) model is based on physiological processes and simulates the gradual decline of solar radiation, carbon balance, water balance, and many other physiological processes, and involves many equations [9]. The most commonly used parameter optimization methods include Markov chain Monte Carlo (MCMC), the annealing method (AM), the genetic algorithm (GA), and particle swarm optimization (PSO).

The MCMC is a kind of Monte Carlo method which is performed by computer under the framework of Bayesian theory. In 1953, Metropolis et al. [10] considered the common Boltzmann distribution sampling problem in physics and proposed the MCMC method (also known as the Metropolis method) for the first time. To improve the sampling efficiency of the MCMC method, Chib et al. created the Metropolis–Hastings (M-H) algorithm by modifying the acceptance rate in sampling based on the Metropolis algorithm (the acceptance rate on both sides of the detailed and stable condition equation was enlarged in the same proportion), thereby increasing the skip acceptance rate in sampling [11,12]. To deal with the high-dimensional distribution of parameters, Smith et al. proposed a sampling method with a higher acceptance rate using the Gibbs random area; the sampling efficiency of this method under high-dimensional conditions was significantly improved [13].

The MCMC parameter estimation method is widely used in ecological research [14,15,16]. For parameter estimation using nonlinear optimization methods (e.g., the 3PG model), the choice of optimization parameters and the amount of calculation are very important. MCMC can combine prior knowledge of the parameters and observational data to obtain the posterior distribution of the model parameters; the posterior values of the parameters can then be used as the parameter calibration result, and the optimized model can be compared with the original model.

This study selects the widely used 3PG model as the research object. Taking a fir forest as an example, firstly, the parameters of the model are screened via sensitivity analysis. Then, the MCMC parameter optimization method is used to optimize the sensitive model parameters. The results can be used to develop a universal method for the optimization of forestry model parameters and provide guidance for the parameter correction and application of forestry models.

2. Materials and Methods

2.1. Study Site

This paper takes Shunchang County in Nanping as the research area. Shunchang County is located between 26°39’ and 27°12’ N latitude and between 117°30’ and 118°14’1 E longitude. It is situated in the northwest of Fujian Province and covers an area of 1985 square kilometers. The climate is a mid-subtropical maritime monsoon climate and is affected by the continental climate. The average annual temperature is 18.5 °C and the annual average precipitation is 1756 mm. The frost-free period is 305 days. The annual average sunshine is 1740.7 h. The county’s forest coverage rate is 75.6% and the main forest vegetation types are fir, Masson pine, and broad-leaved forest (Figure 1).

The input data of the 3PG model include model parameters, site parameters, stand parameters, meteorological data, and observational data [17]. In this study, we cut down Cunninghamia lanceolata trees from 0 to 9 years old in the Shunchang Forest Farm in Nanping City to obtain some model parameters, stand parameters, and observational data. Some model parameters were obtained from the forest resources inventory data of Nanjing County for 2003, while others came from literature reference values and default values. The meteorological data were provided by the Weather Bureau of Fujian Province, including data from 1994 to 2002. The monthly maximum temperature, minimum temperature, average rainfall, and average frost days were interpolated for the study area using the ANUSPLIN (Hutchinson) [18]. The monthly average meteorological data are listed in Table 1.

Figure 2 shows a flowchart of the methodology followed in this study. The 3PG stand model was selected as the ecological prediction model, and the input parameters of the model were obtained by observation of Cunninghamia lanceolata at the sample-plot scale. Based on parameter calibration and model localization, the 3PG model can accurately simulate the vegetation growth process, biomass, diameter at breast height (DBH), and the leaf area index (LAI). Moreover, the MCMC method based on Bayesian theory can obtain the posterior distribution of the model parameters and can thus be used to estimate the parameters; the posterior distribution of parameters can quantitatively express the uncertainty of model parameters under the existing observation conditions. Through the optimization of model parameters, the simulation of future changes in biomass of Chinese fir forest was improved.

2.2. Methods

2.2.1. 3PG Model

The 3PG model was developed by Landsberg and Waring in 1997 and is based on the photosynthetic physiological processes of plants [19]. The model takes the stand-scale as the spatial scale and the month as the timescale and considers the complete carbon balance in the actual environment as well as the climatic conditions, site conditions, management measures, and tree physiological characteristics [20]. The model is mainly composed of three modules. The first module is the carbon fixation module. This module mainly uses the Beer–Lambert extinction formula to simulate the absorption of solar radiation by trees, and then uses the canopy quantum efficiency conversion to estimate primary production. This process accounts for the forest age, temperature, soil moisture content, frost days, and soil type. The effects of soil fertility and other factors on carbon fixation were simulated dynamically by a series of functional relationships. The second module is the biomass allocation model. This module mainly simulates the distribution of fixed carbon among leaves, trunks, and roots. By accounting for the regeneration of the root system, the leaf litterfall rate, and natural stand-thinning, this module mainly calculates the biomass allocation and the changes of various tree organs using an allometric growth equation and the 3/2 self-thinning rule. The third module is the water balance module, which mainly consists of dynamic models related to soil water, such as models of rainfall, evaporation, and artificial irrigation, according to the Penman–Monteith equation. This module considers the influence of tree age, solar radiation, the difference in water vapor pressure, and canopy quantum efficiency on water balance. Detailed descriptions of the parameters of the 3PG model are given in Table 2.

2.2.2. Sensitivity Analysis

Sensitivity analysis is an important step in the application of the 3PG model. The optimization of non-sensitive parameters will not improve the accuracy of the model and will increase the amount of calculation. The sensitivity analysis of stand model parameters can adequately distinguish and define the sensitivity and importance of model parameters and provide a basis for the selection of sensitive parameters. Before parameter optimization, sensitivity analysis of the parameters was carried out, the parameters which have the largest influence on the model simulation were selected, and then the parameters were optimized. This can not only improve the efficiency of the optimization algorithm but also reduce the calculation time.

There is a built-in sensitivity analysis table in the 3PG model, which can be used for the sensitivity analysis of site factors and model parameters. To test the influence of changes to parameter values on the output of the 3PG model, parameter sensitivity analysis was performed. The initial condition was that the other operating parameters of the model remain unchanged. Only by changing the values of the parameters, the parameters that are sensitive to the biomass and growth of the Chinese fir plantation as simulated by the 3PG model were explored, and the selected sensitive parameters were optimized by MCMC.

The model output sensitive to the parameter is as follows:

λ_{1} (X_{1}, p) = \frac{p}{X} \frac{\partial X}{\partial p}

(1)

Where λ₁ is the relative sensitivity, X is the model output value, and p is the model parameter. If the parameter is not sensitive to the model output value, λ₁ is 0. When p increases, X also increases, and λ₁ is positive; otherwise, λ₁ is negative.

The approximate value of λ₁ that can be obtained by the finite difference method is expressed as follows:

λ_{1} (X, p) = \frac{p_{0}}{X_{0}} \frac{X_{+} - X_{-}}{2 δ_{p}}

(2)

Where δ_p indicates the size of the change range of the parameter p, and X₀ = X(P₀).X_-=X(p₀-δ_p),X₊=X(p₀+δ_p)

2.2.3. MCMC

Bayesian theory can combine the prior knowledge of model parameters and the corresponding observations of the model output to realize the posterior estimation of model parameters. The MCMC method involves constructing a Markov chain with the parameter posterior distribution as a stationary distribution under the framework of Bayesian theory, so as to obtain the posterior samples of the parameters and infer the numerical characteristics of the parameters based on these samples. Bayesian theory is expressed in the following formula:

p (θ / y) = \frac{f (y / θ) g (θ)}{\int f (y / θ) g (θ) d (θ)}

(3)

where θ and y represent the parameters and simulated output values of the 3PG model (e.g., biomass and diameter at breast height), respectively; p(θ/y) is the posterior probability density function of the parameters; and f(y/θ) is the observed data. The conditional probability density under the prior parameter value is also called the likelihood function. g(θ) is the prior distribution of the parameter. Formula (3) can be changed to the following:

p (θ / x, y) = \frac{f (y / θ, x), g (θ / x)}{\int f (y / θ, x) g (θ / x) d (θ)} \approx f (y / θ, x) g (θ)

(4)

Commonly used MCMC sampling methods include the Metropolis algorithm, the Metropolis–Hastings (M-H) algorithm, the Gibbs sampling algorithm, and the adaptive Metropolis algorithm. This study used the M-H sampling method. The steps for this method are as follows: (1) Randomly generate initial estimates of parameters within a range of values; (2) generate new parameter values based on the posterior distribution of parameters assumed in advance; (3) calculate the acceptance probability; (4) generate a uniformly distributed random number U in the interval [0,1]; (5) if p ≥ U, the model accepts the new parameter value; otherwise, it will reject the new parameter value; and (6) repeat steps (2)–(5) until enough samples are obtained, and finally obtain the posterior estimation of the parameters.

The initial values of parameters are based on the observed values, and the parameters are roughly adjusted manually by a trial-and-error method to make the trends of stand stocking and diameter at break height approximately the same. The prior distribution can be divided into two categories: conjugate prior and non-informative prior. Non-informative prior refers to a kind of priori constructed when the value range of parameters and their status in the model are known but nothing is known about other parameters. In this paper, we chose the non-informative prior which accords with the uniform distribution.

The purpose of using the MCMC method is to obtain a posterior sample of model parameters by combining the measured data with the prior knowledge of model parameters. We input the initial parameter values and posterior parameter values into the 3PG model, respectively, comparing the simulated values with the measured values. The root-mean-square error (RMSE) between simulation value and the corresponding observed value was taken as the accuracy evaluation index. The formula is as follows:

{RMSE}_{stem} = \sqrt{\sum_{i = 1}^{i = n} {({Stem}_{i} - {Stem}_{i}^{obs})}^{2} / n}

(5)

where Stem represents stem biomass including branches and bark, Stem_i represents the i-th Stem simulation value, Stem^obs represents the i-th corresponding observed value, and n represents the number of Stem observations.

3. Results

3.1. 3PG Model

By changing the values of the parameters (taking the default value of ±30% for the upper and lower bounds) and observing the sensitivity of the results, the sensitive parameters were selected (Figure 3). The results showed that the sensitive parameters were nWs (power in the stem mass vs. diameter relationship), aWs (constant in the stem mass vs. diameter relationship), alphaCx (maximum canopy quantum efficiency), k (extinction coefficient for the absorption of PAR by canopy), pRx (maximum fraction of NPP to roots), pRn (minimum fraction of net primary productivity (NPP) to roots), and CoeffCond (defines stomatal response to vapor pressure deficit (VPD)).

Through the finite difference method of Formula (2), the sensitivities of the above screened parameters to stand stocking and DBH were obtained, and the sensitivity levels were divided and ranked according to the distribution of values [19] (Figure 4 and Table 3).

The error lines above the points indicate that the model parameters have a positive impact on the stand stocking and DBH, which indicates that the parameter values increase, so stand stocking and DBH increase, while the error lines below the points indicate that the model parameters have a negative impact on the stand stocking and DBH, which indicates that the parameter values increase, and stand stocking and DBH decrease. Ranking scheme: ≤0.075: Grade 0; 0.075–0.25: Grade 1; 0.25–0.4: Grade 2; ≥0.4: Grade 3.

The results showed that alphaCx, CoeffCond, and pRn corresponded to stemNo sensitivity Grade 3, pRx to stemNo sensitivity Grade 2, k to stemNo sensitivity Grade 1, nWs and aWs to stemNo sensitivity grade 0, nWs and AlphaCx to DBH sensitivity Grade 3, CoeffCond to DBH sensitivity Grade 2, and aWs, k, pRx, and pRn to DBH sensitivity Grade 1.

Finally, seven parameters were selected as the parameters to be calibrated, namely nWs, aWs, alphaCx, k, pRx, pRn, and CoeffCond.

3.2. MCMC Result

According to the default value of the 3PG model and analytical tree data, the initial values and ranges of the seven parameters to be optimized were established (Table 4).

The seven most sensitive parameters in the biomass simulation using the 3PG model were calibrated via the MCMC method, and the range of each parameter was expanded by a certain number of times. After 50,000 iterations of parameter expansion by 300%, it was found that the values of the final parameters are relatively stable—that is, further expanding the range or increasing the number of iterations does not lead to a large change in the parameter values. Compared with the prior distribution of each parameter, the posterior distribution differed greatly (Figure 5). The posterior probability distributions of nWs, aWs, alphaCx, pRx, pRn, and CoeffCond approximately conform to a normal or skewed distribution, and the peak value is prominent, which indicates that the uncertainty of these parameters is relatively small, and k’s peak value is not prominent, showing an irregular distribution.

For parameters with a prominent, near-normal, or skewed distribution, the peak value is taken as the result (Table 4). We took the following values as input values: nWs = 2.76, aWs = 0.42, alphaCx = 0.16, pRx = 0.56, pRn = 0.24, and CoeffCond = 0.11. For the one parameter that was not prominent, namely k, the default value was taken as the model input value, that is, k = 0.5.

3.3. Comparative Analysis

Subsequently, the initial values and posterior values of parameters were input into the model for running and were compared with the observed values. The results are shown in Figure 6. In the 3PG simulation, the average measured stand biomass over the nine study years was 164.9 t/ha, and the simulated values before and after adjustment were 138.0 t/ha and 160.1 t/ha, respectively. Compared with the initial value, the accuracy of the stem simulation was improved by 13.92% (posterior value), and the deviation decreased from 16.36% to 2.44%. This shows that the posterior value calibrated by MCMC can achieve a better fit with the observational data.

The results show that the RMSEs((Root Mean Square Error) of the stem values with the initial value (default parameter) and posterior value in the model simulation are 1.24 and 0.98, respectively; the RMSEs of the height values are 0.34 and 0.32, respectively; and the RMSEs of the DBH values are 0.71 and 0.69, respectively. All indicators show that the proposed model has a higher accuracy than the initial 3PG model (Table 5).

4. Discussion

4.1. Parameter Sensitivity

In this study, the sensitivity of 63 model parameters to stand stocking (stemNo) and DBH was analyzed by sensitivity analysis in the 3PG model. According to the sensitivity results for stemNo and DBM, the sensitive parameters were power in the stem mass versus diameter relationship (nWs), the constant in the stem mass versus diameter relationship (aWs), the maximum fraction of NPP to roots (pRx), and the minimum fraction of NPP to roots (pRn), which are all related to the allocation relationship and proportion of biomass, as well as the maximum canopy quantum efficiency (alphaCx), which defines stomatal response to VPD (CoeffCond), and the extinction coefficient for PAR absorption by canopy (k).

These sensitivities can be attributed to the fact that the canopy quantum efficiency and extinction coefficient affect the interception of light energy, photosynthetic production, and respiration, which is the main reason for limiting the activity of photosynthesis [20,21]; they also affect the power value of the relationship between dry biomass and DBH, the constant value of the relationship between dry biomass and DBH, the maximum net primary production allocated to the roots, and the minimum net primary production allocated to the roots. The results showed that the biomass and DBH were influenced by the distribution and proportion of biomass. Leaf stomata were the main outlet for water vapor from the plant body to the atmosphere, and the response of stomata to VPD affects transpiration, photosynthesis, and respiration; therefore, stand stocking and DBH were more sensitive to the above parameters.

L-J Esprey et al. [20,22] found that the parameters Y, alphaCx, MaxCond, rhoMax, fracBB0, fracBB1, fCalpha700, Topt, and k were sensitive to stand stocking and DBH. Meanwhile, Xiaodong Song et al. [23,24] found that nWS, aWS, alphaCx, MaxCond, k, tWaterMax, FR, pRx, pRn, CoeffCond, tRho, fracBB1, and fracBB2 were sensitive parameters. However, in the present study, Y, rhoMax, fracBB0, fracBB1, and Topt were not found to be sensitive parameters. The main reason for this is that the species, environment, and method in this work are different from those in the aforementioned two studies. Although local sensitivity analysis is simple and easy to conduct, it cannot test the influence of the interaction between model parameters on the output of the model. Therefore, global sensitivity analysis should be used for analysis in subsequent studies.

4.2. Parameter Optimization using MCMC

According to the results of the sensitivity screening, seven parameters in the 3PG model which have a large influence on stand stocking and DBH were selected for optimization. Compared with the prior distribution of the parameters, the posterior distribution of the parameters was greatly different. The more the scope of the posterior distribution was narrowed, the smaller the uncertainty of the optimized parameters, which indicates the validity of the calibration-parameter selection. The results show that the parameters tend to be stable and the model can simulate the observed values adequately. The results show that the MCMC method can be used to obtain the optimal parameters stably, and additionally verify the feasibility of this method for parameter adjustment in the 3PG model.

For example, for the parameter k, the posterior distribution is flat and irregular, and the peak value is not prominent. This shows that different parameter combinations can obtain the same model output value; this phenomenon may be due to the redundancy and correlation of model parameters, model structure error, and input-output error. In addition to the parameters that are sensitive to stand stocking and DBH, the 3PG model can also screen the parameters that are sensitive to other parameters, such as leaf biomass, root biomass, stand volume, and leaf area index. The methods that were used for sensitivity analysis in the present study can also be used for the selection and optimization of these other sensitive parameters [20].

4.3. Model Uncertainty

The uncertainty of model parameters is mainly due to the fact that some parameters in the model are difficult to obtain directly. The goal of traditional parameter estimation methods is to find a set of optimal parameters in some specific model structures. In the model calibration, most researchers adjust some specific parameters and choose the parameters with the smallest error as the calibration parameters according to the error between the model simulation and the measured value. However, this method is subjective [25,26]. Meanwhile, some scholars construct the objective function of the difference between the simulated value and the observed value to minimize the difference between the model simulation and the measured value, so as to obtain the estimated value of the model parameters. As it is impossible to obtain accurate model parameters in complex process-based models, it is not feasible to obtain only one parameter estimate using an optimization algorithm [27]. The MCMC method based on Bayesian theory can obtain the posterior distribution of model parameters and accordingly has been widely applied [28,29,30]. By minimizing the objective function, the best fitting degree between the model output and the actual observation data can be achieved, which can effectively reduce the uncertainty of model parameters, improve the accuracy of model simulation, and enhance the practical application value of ecological models.

In addition to the model parameters, the structure of the model itself and the meteorological driving data are the main sources of uncertainty in the 3PG model. The uncertainty caused by the model structure is mainly due to the fact that it is difficult to quantitatively and accurately describe the processes of photosynthesis, biomass allocation, and water balance, and that the effects of factors such as plant respiration, extreme weather, and disasters are not considered in the model, which also affects the simulation of stand growth. Meteorological data are important driving data for ecological models. In this paper, to obtain a spatiotemporally continuous meteorological driving dataset for the study area, an interpolation method was used. Due to the uneven spatial distribution of discontinuous macro-phenomena such as temperature and rainfall, the uncertainty caused by the interpolation method will also affect the stand prediction, which is one of the bottlenecks restricting the practical application of the model for stand management. The atmospheric system is highly nonlinear and chaotic, and therefore uncertainty is inevitable in weather forecasting, and the same is true for the 3PG model, which is driven by weather [31,32,33].

The uncertainty of simulation values increases with increasing simulation time. To overcome this problem, a sequential parameter estimation method or data assimilation method can be used to dynamically optimize the model parameters or state variables according to new observational data that are obtained in the future so as to reduce error propagation and more accurately predict the carbon-cycle processes of ecosystems under future climatic conditions. Additionally, researchers should further study how to quantify and reduce the uncertainty caused by errors in the model input data and by the model structure itself.

5. Conclusions

In this study, the sensitivity of 63 parameters in the 3PG model to stand stocking and DBH was analyzed based on the analytical data of trees from 0 to 9 years old, meteorological data from 1994–2003, and forest inventory data from the Shunchang Forest Farm in Nanping. The seven parameters that had the greatest influence on biomass and DBH were selected, and these parameters were then optimized using the Markov chain Monte Carlo (MCMC) method. The conclusions are as follows:

(1): Among the 63 parameters of the 3PG model, the parameters that are most sensitive to stand stocking and DBH were nWs, aWs, alphaCx, k, pRx, pRn, and CoeffCond.
(2): The parameters that have the greatest influence on stand stocking are alphaCx, CoeffCond, pRn, pRx, k, nWs, and aW, and the parameters with the greatest influence on DBH are nWs, alphaCx, CoeffCond, aWs, pRn, pRx, and k.
(3): The posterior probability distributions of nWs, aWs, alphaCx, pRx, pRn, and MaxCond have an approximately normal or skewed distribution with a prominent peak value; however, the peak value of k is not prominent, showing an irregular distribution.
(4): Compared with the simulation results using the default parameters, the RMSEs of the stem values with the initial value (default parameter) and posterior value in the model simulation are 1.24 and 0.98, respectively; the RMSEs of the height values are 0.34 and 0.32, respectively; and the RMSEs of the DBH values are 0.71 and 0.69, respectively.

Author Contributions

C.L. designed the study, finished model running, drew the figures and tables, and wrote the main body of the manuscript; X.Z. collected the data, provided the methodology, and analyzed the data; Y.R. edited the draft, acquired the funding, and revised the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the National Science Foundation of China (31670645, 31972951, 41801182, and 41807502), the National Social Science Fund (17ZDA058), the National Key Research Program of China (2016YFC0502704), Fujian Provincial Department of S&T Project (2018T3018), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA23020502), the Ningbo Municipal Department of Science and Technology (2019C10056), the Xiamen Municipal Department of Science and Technology (3502Z20130037 and 3502Z20142016), the Key Laboratory of Urban Environment and Health of CAS (KLUEH-C-201701), the Key Program of the Chinese Academy of Sciences (KFZDSW-324), and the Youth Innovation Promotion Association CAS (2014267).

Acknowledgments

We acknowledge Xiaodong Song from Zhejiang University for their contribution in making the simulation possible. The Forestry Bureau of Nanping City provided the forest inventory data. The Weather Bureau of Fujian Province provided meteorological data. We thank Zhengwei Li for his advice and assistance in data processing. Thanks to Zhifeng Wu for his spiritual support. Thanks to Kuo Liao provided resources. Especially, the authors would also like to thank the anonymous reviewers for their very important suggestions, which were particularly helpful for improving this paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

Luo, Y.; White, L.W.; Canadell, J.G.; DeLucia, E.H.; Ellsworth, D.S.; Finzi, A.; Lichter, J.; Schlesinger, W.H. Sustainability of terrestrial carbon sequestration: A case study in Duke Forest with inversion approach. Glob. Biogeochem. Cycles 2003, 17. [Google Scholar] [CrossRef]
Raupach, M.R.; Rayner, P.J.; Barrett, D.J.; DeFries, R.S.; Heimann, M.; Ojima, D.S.; Quegan, S.; Schmullius, C.C. Model-data synthesis in terrestrial carbon observation: Methods, data requirements and data uncertainty specifications. Glob. Chang. Biol. 2005, 11, 378–397. [Google Scholar]
Damian, J.B. Steady state turnover time of carbon in the Australian terrestrial biosphere. Glob. Biogeochem. Cycles 2002, 16, 55. [Google Scholar]
Wang, Y.P.; Barrett, D.J. Estimating regional terrestrial carbon fluxes for the Australian continent using a multiple- constraint approach II. The Atmospheric constraint. Tellus B Chem. Phys. Meteorol. 2003, 55, 270–289. [Google Scholar]
Rayner, P.J.; Scholze, M.; Knorr, W.; Kaminski, T.; Giering, R.; Widmann, H. Two decades of terrestrial carbon fluxes from a carbon cycle data assimilation system (CCDAS). Glob. Biogeochem. Cycles 2005, 19, GB2026. [Google Scholar] [CrossRef] [Green Version]
Williams, M.; Schwarz, P.A.; Law, B.E.; Irvine, J.; Kurpius, M.R. An improved analysis of forest carbon dynamics using data assimilation. Glob. Chang. Biol. 2005, 11, 89–105. [Google Scholar]
Luo, Y.; Weng, E.; Wu, X.; Gao, C.; Zhou, X.; Zhang, L. Parameter identifiability, constraint, and equifinality in data assimilation with ecosystem models. Ecol. Appl. 2009, 19, 571–574. [Google Scholar]
Lin, J.C.; Pejam, M.R.; Chan, E.; Wofsy, S.C.; Gottlieb, E.W.; Margolis, H.A.; McCaughey, J.H. Attributing uncertainties in simulated biospheric carbon fluxes to different error sources. Glob. Biogeochem. Cycles 2011, 25. [Google Scholar] [CrossRef]
Almeida, A.C.; Landsberg, J.J.; Sands, P.J. Parameterisation of 3-PG model for fast-growing Eucalyptus grandis plantations. For. Ecol. Manag. 2004, 193, 179–195. [Google Scholar]
Metropolis, N.; Rosenbluth, A.W.; Rosenbluth, M.N.; Teller, A.H.; Teller, E. Equation of state calculations by fast computing machines. J. Chem. Phys. 1953, 21, 1087–1092. [Google Scholar] [CrossRef] [Green Version]
Chib, S.; Greenberg, E. Understanding the Metropolis-Hastings algorithm. Am. Stat. 1995, 49, 327–335. [Google Scholar]
Hastings W, K. Monte Carlo sampling methods using Markov chain and their applications. Biometrika 1970, 57, 97–109. [Google Scholar]
Geman, S.; Geman, D. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE T. Pattern. Anal. 1987, 6, 721–741. [Google Scholar]
Zobitz, J.M.; Desai, A.R.; Moore, D.J.P.; Chadwick, M.A. A primer for data assimilation with ecological models using Markov Chain Monte Carlo (MCMC). Oecologia 2011, 167, 599–611. [Google Scholar] [PubMed]
Ren, X.; Honglin, H.; Moore, D.J.P.; Zhang, L.; Liu, M.; Li, F.; Yu, G.; Wang, H. Uncertainty analysis of modeled carbon and water fluxes in a subtropical coniferous plantation. J. Geophys. Res. Biogeosci. 2013, 118, 1674–1688. [Google Scholar]
Ricciuto, D.M.; King, A.W.; Dragoni, D.; Post, W.M. Parameter and prediction uncertainty in an optimized terrestrial carbon cycle model: Effects of constraining variables and data record length. J. Geophys. Res. 2011, 116. [Google Scholar] [CrossRef] [Green Version]
Landsberg, J.-J.; Waring, R.-H. A generalised model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning. For. Ecol. Manag. 1997, 95, 209–228. [Google Scholar]
Hutchinson, M.F. The application of thin plate smoothing splines to continent-wide data assimilation. Data Assim. Syst. 1991, 27, 104–113. [Google Scholar]
Coops, N.-C.; Waring, R.-H. Assessing forest growth across southwestern Oregon under a range of current and future global change scenarios using a process model, 3-PG. Glob. Chang. Biol. 2001, 7, 15–29. [Google Scholar]
Esprey, L.-J.; Sands, P.-J.; Smith, C.-W. Understanding 3-PG using a sensitivity analysis. For. Ecol. Manag. 2004, 193, 235–250. [Google Scholar]
Battaglia, M.; Sands, P. Application of sensitivity analysis to a model of Eucalyptus globulus plantation productivity. Ecol. Model. 1998, 111, 237–259. [Google Scholar] [CrossRef]
Battagalia, M.; Sands, P.J. Process-based forest productivity models and their application in forest management. For. Ecol. Manag. 1998, 102, 13–32. [Google Scholar] [CrossRef]
Song, X.; Bryan, B.A.; Almeida, A.C.; Paul, K.I.; Zhao, G.; Ren, Y. Time-dependent sensitivity of a process-based ecological model. Ecol. Model. 2013, 265, 114–123. [Google Scholar] [CrossRef]
Song, X.; Bryan, B.A.; Paul, K.I.; Zhao, G. Variance-based sensitivity analysis of a forest growth model. Ecol. Model. 2012, 247, 135–143. [Google Scholar] [CrossRef]
Mertens, J.; Madsen, H.; Feyen, L.; Jacques, D.; Feyen, J. Including prior information in the estimation of effective soil parameters in unsaturated zone modelling. J. Hydrol. 2004, 294, 251–269. [Google Scholar] [CrossRef]
Seidel, S.J.; Palosuo, T.; Thorburn, P.; Wallach, D. Towards improved calibration of crop models: Where are we now and where should we go. Eur. J. Agron. 2018, 94, 25–35. [Google Scholar] [CrossRef]
Gauch, H.G.; Hwang, J.G.; Fick, G.W. Model evaluation by comparison of model-based predictions and measured values. Agron. J. 2003, 95, 1442–1446. [Google Scholar] [CrossRef] [Green Version]
Irmak, A.; Jones, J.W.; Batchelor, W.D.; Paz, J.O. Estimating spatially variable soil properties for application of crop models in precision. Agriculture 2001, 44, 1343. [Google Scholar] [CrossRef]
Romanowicz, R.J.; Beven, K.J. Comments on generalised likelihood uncertainty estimation. Reliab. Eng. Syst. Saf. 2006, 91, 1315–1321. [Google Scholar] [CrossRef]
Moulin, S.; Bondeau, A.; Delecolle, R. Combining agricultural crop models and satellite observations: From field to regional scales. Int. J. Remote Sens. 1998, 19, 1021–1036. [Google Scholar] [CrossRef]
Iizumi, T.; Yokozawa, M.; Nishimori, M. Parameter estimation and uncertainty analysis of a large-scale crop model for paddy rice: Application of a Bayesian approach. Agric. For. Meteorol. 2009, 149, 333–348. [Google Scholar] [CrossRef]
Marin, F.; Jones, J.W.; Boote, K.J. A stochastic method for crop models: Including uncertainty in a sugarcane model. Agron. J. 2017, 109, 483–495. [Google Scholar] [CrossRef]
Kalnay, E.; Kanamitsu, M.; Kistler, R.; Collins, W.; Deaven, D.; Gandin, L.; Iredell, M.; Saha, S.; White, G.; Woollen, J.; et al. The NCEP/NCAR 40-year reanalysis project Bull. Am. Meteor. Soc 1996, 77, 437–471. [Google Scholar] [CrossRef] [2.0.CO;2" target='_blank'>Green Version]

Figure 1. Map of the study area.

Figure 2. Flowchart of the methodology followed in this study. 3PG: Physiological Principles in Predicting Growth model; MCMC: Markov chain Monte Carlo method; stemNo: stand stocking; WF: foliage biomass; WR: root biomass; WS: stem biomass including branches and bark; LAI: canopy leaf area index (LAI); standvol: stand volume excluding branches and bark; MAI: mean annual volume increment; avDBH: stand-based mean diameter at breast height (DBH).

Figure 3. The effects of the adjustment of model parameters on their sensitivity to stand stocking (stemNo) and diameter at breast height (DBH).

Figure 4. Sensitivity grades of model parameters to stand stocking and DBH.

Figure 5. The posterior distributions of seven parameters. Note: The x-axis represents the value of the parameter and the y-axis represents the probability density.

Figure 6. Comparison of stand biomass before and after adjustment. Note: the red dots represent the observational data of Cunninghamia lanceolata trees from the Shunchang Forest Farm, the thin line represents the predicted stand biomass value (according to the observational data, the 3PG model can automatically generate a prediction curve), and the thick line represents the simulated stand biomass value.

Table 1. Monthly average meteorological data.

Month	Maximum Temperature/°C	Minimum Temperature/°C	Precipitation/mm	Solar Radiation (MJm²d⁻¹)	Frost Days/d
January	22.45	–1.88	36.29	25.45	0
February	23.94	2.16	89.38	23.27	0
March	24.53	4.02	915.97	19.85	0
April	25.51	10.69	146.07	14.36	0
May	26.96	13.95	148.8	10.08	0.67
June	27.6	18.39	224.55	7.91	2.25
July	28.58	21.43	187.58	8.46	3.75
August	27.97	21.99	295.55	11.82	2.5
September	27.18	17.36	210.4	15.62	0.33
October	25.9	9.75	33.53	20.44	0
November	23.79	3.83	21.32	23.61	0
December	21.43	–0.44	35.42	25.64	0

*: multiplication sign.

Table 2. Detailed descriptions of the parameters of the Physiological Principles in Predicting Growth (3PG) model.

Parameter Name	Definition	Unit	Symbol	Value
pFS2	Ratio of foliage:stem partitioning at B = 2 cm	-	p2	1
pFS20	Ratio of foliage:stem partitioning at B = 20 cm	-	p20	0.15
aWS	Constant in stem mass vs. diameter relationship	-	aS	0.095
nWS	Power in stem mass vs. diameter relationship	-	nS	2.4
pRx	Maximum fraction of NPP to roots	-	ηRx	0.8
pRn	Minimum fraction of NPP to roots	-	ηRn	0.25
gammaF0	Litterfall rate at t = 0	month⁻¹	γF0	0.027
gammaF1	Litterfall rate for mature stands	month⁻¹	γF1	0.001
tgammaF	Age at which litterfall rate has median value	month⁻¹	Ft	12
Rttover	Average monthly root turnover rate	month⁻¹	γR	0.015
Tmin	Minimum temperature for growth	°C	Tmin	8.5
Topt	Optimum temperature for growth	°C	Topt	16
Tmax	Maximum temperature for growth	°C	Tmax	40
kF	Number of production days lost for each frost day	days	kF	0
m0	Value of m when FR = 0	-	m0	0
fN0	Value of fN when FR = 0	-	fN0	1
fNn	Power of (1–FR) in fN	-	nfN	0
CoeffCond	Defines stomatal response to VPD	mbar	kD	0.05
fCalpha700	Quantum efficiency at 700 ppm	-	fC700	0.7
fCg700	Canopy conductivity at 700 ppm	-	fCg700	0.7
SWconst	Moisture ratio deficit which gives f_θ= 0.5	-	c	0.7
SWpower	Power of moisture ratio deficit in f_θ	-	n	9
MaxAge	Maximum stand age used to compute relative age	yr	tx	50
nAge	Power of relative age in f_age	-	nage	4
rAge	Relative age to give f_age = 0.5	-	rage	0.95
MinCond	Minimum canopy conductance	m s⁻¹	gCn	0
MaxCond	Maximum canopy conductance	m s⁻¹	gCx	0.02
LAIgcx	Canopy LAI for maximum canopy conductance	m² m⁻²	LCx	3.33
BLcond	Canopy boundary layer conductance	m s⁻¹	gB	0.2
gammaN0	Seedling mortality rate (t = 0)	yr⁻¹	N0	0.03
gammaNx	Mortality rate for older stands (large t)	yr⁻¹	N1	0.001
tgammaN	Age at which γN = 1/2(γN0+γN1)	yr	tN	12
ngammaN	Shape of mortality response	-	nN	0.015
wSx1000	Maximum stem mass per tree at 1000 trees/ha	kg/tree	wSx1000	300
thinPower	Power in self-thinning law	-	nN	3/2
mF	Fractions of mean foliage, root, and stem biomass pools per tree on each dying tree	-	mF	0
mR		-	mR	0.2
mS		-	mS	0.2
SLA0	Specific leaf area at a stand age of 0	m² kg⁻¹	σ0	11
SLA1	Specific leaf area for mature stands	m² kg⁻¹	σ1	4
tSLA	Age at which specific leaf area = (SLA₀+SLA₁)/2	yr	t	2.5
MaxIntcptn	Maximum fraction of rainfall intercepted by canopy	-	iRx	0.15
LAImax-Intcptn	LAI for maximum rainfall interception	m² m⁻²	Lix	0
k	Extinction coefficient for PAR absorption by canopy	-	k	0.5
fullCanAge	Age at full canopy cover	yr	tc	0
alphaCx	Maximum canopy quantum efficiency	-	Cx	0.06
Y	Ratio NPP/GPP	-	Y	0.47
fracBB0	Branch and bark fraction at a stand age of 0	-	p	0.75
fracBB1	Branch and bark fraction for mature stands	-	p	0.15
tBB	Age at which p_BB = 1/2(P_BB0+ P_BB1)	yr	t_BB	2
rhoMax	Minimum basic density for young trees	t m⁻³	ρ0	0.5
rhoMax	Maximum basic density for older trees	t m⁻³	ρ1	0.5
tRho	Age at which p = 1/2 density of old and young trees	yr	t_ρ	4
aH	Constant in stem diameter to height relationship	yr	aH	0
nHB	Power of DBH in stem height relationship	-	nHB	0
nHN	Power of stocking in stem height relationship	-	nHN	0
aV	Constant in stem diameter to volume relationship	-	aV	0
nVB	Power of DBH in stem volume relationship	-	nVB	0
nVN	Power of stocking in stem volume relationship	-	nVN	0
Qa	Intercept in net radiation vs. solar radiation relationship	W m⁻²	Qa	−90
Qb	Slope of net radiation vs. solar radiation relationship	-	Qb	0.8
gDM_mol	Molecular weight of dry matter	g mol⁻¹		24
molPAR_MJ	Conversion of solar radiation to PAR	mol MJ⁻¹		2.3

Note: DBH: diameter at breast height. LAI: leaf area index. NPP: net primary productivity. GPP: gross primary productivity. VPD: vapor pressure deficit. PAR: photosynthetically active radiation.

Table 3. Results of the sensitivity analysis for stand stocking (stemNo) and DBH predicted by the 3PG model.

Parameter	Sensitivity Value (λ₁)		Grade (λ₁)
Parameter	StemNo	DBH	StemNo	DBH
nWS	–0.06	–3.10	0	3
aWS	–0.02	–0.21	0	1
alphaCx	1.17	0.48	3	3
k	0.21	0.09	1	1
pRx	–0.35	–0.13	2	1
pRn	–0.50	–0.21	3	1
CoeffCond	–0.73	–0.31	3	2

Table 4. Initial values, ranges, and posterior distributions of the seven parameters to be optimized.

Parameter	Unit	Initial Value	Range	Mode	Mean + SD
nWS	-	2.4	[0,8]	2.76	3.08 ± 0.17
aWS	-	0.095	[0,2]	0.42	0.58 ± 0.02
alphaCx	-	0.06	[0,0.6]	0.16	0.26 ± 0.1
k	-	0.5	[0,1.5]	1.04	0.80 ± 0.24
pRx	-	0.6	[0,2]	0.56	0.71 ± 0.09
pRn	-	0.25	[0,1.6]	0.24	0.59 ± 0.1
CoeffCond	1/mbar	0.05	[0,0.65]	0.11	0.24 ± 0.01

Table 5. Comparison of parameters before and after optimization.

	Height		DBH		Stand Biomass
	Before Optimization	After Optimization	Before Optimization	After Optimization	Before Optimization	After Optimization
RMSE	0.34	0.32	0.71	0.69	1.24	0.98

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Liu, C.; Zheng, X.; Ren, Y. Parameter Optimization of the 3PG Model Based on Sensitivity Analysis and a Bayesian Method. Forests 2020, 11, 1369. https://doi.org/10.3390/f11121369

AMA Style

Liu C, Zheng X, Ren Y. Parameter Optimization of the 3PG Model Based on Sensitivity Analysis and a Bayesian Method. Forests. 2020; 11(12):1369. https://doi.org/10.3390/f11121369

Chicago/Turabian Style

Liu, Chenjian, Xiaoman Zheng, and Yin Ren. 2020. "Parameter Optimization of the 3PG Model Based on Sensitivity Analysis and a Bayesian Method" Forests 11, no. 12: 1369. https://doi.org/10.3390/f11121369

APA Style

Liu, C., Zheng, X., & Ren, Y. (2020). Parameter Optimization of the 3PG Model Based on Sensitivity Analysis and a Bayesian Method. Forests, 11(12), 1369. https://doi.org/10.3390/f11121369

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu