Dynamic Model for Caragana korshinskii Shrub Aboveground Biomass Based on Theoretical and Allometric Growth Equations

: As one of the ways to achieve carbon neutralization, shrub biomass plays an important role for natural resource management decision making in arid regions. To investigate biomass dynamic variations of Caragana korshinskii , a typical shrub found in the arid desert area of Ningxia, northwest China, we combined a nonlinear simultaneous (NLS) equation system with theoretical growth (TG) and allometric growth (AG) equations. On the basis of a large biomass survey dataset and analytical data of shrub stems, four methods (NOLS, NSUR, 2SLS, and 3SLS) of the NLS equations system were combined with the TG and AG equations. A model was subsequently established to predict the AGB growth of C. korshinskii . The absolute mean residual ( AMR ), root mean system error ( RMSE ), and adjusted determination coefﬁcient (adj- R 2 ) were used to evaluate the performance of the equations. Results revealed that the NSUR method of the NLS equations had better performance than other methods and the independent equations for BD and H growth and AGB. Additionally, the NSUR method exhibited extremely signiﬁcant differences ( p < 0.0001) when compared with the equations without heteroscedasticity on the basis of the likelihood ratio (LR) test, which used the power function (PF) as the variance function. The NSUR method of the NLS equations was an efﬁcient method for predicting the dynamic growth of AGB by combining the TG and AG equations and could estimate the carbon storage for shrubs accurately, which was important for stand productivity and carbon sequestration capacity.


Introduction
Biomass is investigated when studying the carbon storage of forest ecosystems, as it is an important indicator of stand productivity and carbon sequestration capacity and plays an indicative role in forest quality assessment [1,2]. Shrubs are indispensable and important species in the forest ecosystem and play essential roles in ecological processes. Aboveground biomass (AGB) is expressed on an area basis and is central to many ecological processes and services provided by shrublands [3]. It is the main forest vegetation type found in the arid aeolian sand regions (AASRs) of China. Shrubs are widely used for ecological protection and restoration measures and reconstruction projects due to their well-developed root systems and remarkable abilities in wind and sand fixation, soil and water conservation, and drought resistance. Thus, studying shrub biomass is of great importance for ecosystem restoration, regeneration, and security in the AASRs [4].
In recent years, scholars have constructed many biomass models [5,6]. The first study on biomass was conducted in 1876 when Ebermeryer in Germany investigated the dry matter productivity of the Bavarian Forest by measuring the amount of leaf litter and wood weight of certain tree species. The earliest measurements and studies on shrub biomass have been traced back to the early 1960s [7,8]. Lufafa et al. [9] studied the shrub

Study Area
The study area is located in the AASRs, Yanchi County (37 • 4 -38 • 10' N, 106 • 30 -107 • 47 E), Ningxia Province, northwestern China. It is a typical zone consisting of agricultural to pastoral areas with a temperate continental climate [32]. The annual average temperature is 7.7 • C, the mean annual precipitation is 297 mm, and the annual evaporation is 2100 mm [33,34]. The soil consists mainly of sierozem and aeolian sand. Due to sparse precipitation, large evaporation, and a serious shortage of water resources in the county, the tree species mainly include grassland, sandy, and halophyte vegetation, including C. korshinskii, Nitraria sibirica, and Hedysarum scoparium.

Data Collection
C. korshinskii is an important shrub species for soil and water conservation and sandfixation afforestation in northwest and northeast China. It is a perennial legume shrub with strong drought resistance and has a planting area that exceeds 40% [35,36]. A total of 39 temporary sample plots of C. korshinskii were selected in the study area ( Figure 1). Sample plots were square and varied in size, ranging from 100 to 400 m 2 . The BD (mm) at 10 cm aboveground, H (m), crown width (m), and total number of stems (N) of all shrubs in the sample plot were measured. Then, two or three individual shrubs from each sample plot were selected for analysis on the basis of the average BD, H, crown width, and number of stems. A total of 87 individual shrubs were collected. The aboveground parts were collected, and dry samples of the stems and leaves were transported to the laboratory at a constant 85 • C and constant weight. We calculated the moisture content of each component according to the relationship between the fresh and dry mass of each sample, then calculated the dry matter mass of each part of the sample and added each part to obtain the AGB of C. korshinskii (Table 1). Three to five stems with the largest BD and H from individual shrubs were selected, segmented into 20 cm lengths, and cut into disks, which were scanned using WinDENDRO software to count the number of annual rings and measure the diameter and H in increments.

Study Area
The study area is located in the AASRs, Yanchi County (37°4'-38°10´ N, 106°30´-107°47´ E), Ningxia Province, northwestern China. It is a typical zone consisting of agricultural to pastoral areas with a temperate continental climate [32]. The annual average temperature is 7.7 °C, the mean annual precipitation is 297 mm, and the annual evaporation is 2100 mm [33,34]. The soil consists mainly of sierozem and aeolian sand. Due to sparse precipitation, large evaporation, and a serious shortage of water resources in the county, the tree species mainly include grassland, sandy, and halophyte vegetation, including C. korshinskii, Nitraria sibirica, and Hedysarum scoparium.

Data Collection
C. korshinskii is an important shrub species for soil and water conservation and sandfixation afforestation in northwest and northeast China. It is a perennial legume shrub with strong drought resistance and has a planting area that exceeds 40% [35,36]. A total of 39 temporary sample plots of C. korshinskii were selected in the study area ( Figure 1). Sample plots were square and varied in size, ranging from 100 to 400 m 2 . The BD (mm) at 10 cm aboveground, H (m), crown width (m), and total number of stems (N) of all shrubs in the sample plot were measured. Then, two or three individual shrubs from each sample plot were selected for analysis on the basis of the average BD, H, crown width, and number of stems. A total of 87 individual shrubs were collected. The aboveground parts were collected, and dry samples of the stems and leaves were transported to the laboratory at a constant 85 °C and constant weight. We calculated the moisture content of each component according to the relationship between the fresh and dry mass of each sample, then calculated the dry matter mass of each part of the sample and added each part to obtain the AGB of C. korshinskii (Table 1). Three to five stems with the largest BD and H from individual shrubs were selected, segmented into 20 cm lengths, and cut into disks, which were scanned using WinDENDRO software to count the number of annual rings and measure the diameter and H in increments.

Test of Data Normality
The Shapiro-Wilk test was used to test the normality of BD growth, H growth, and AGB data ( Table 2). Because the data showed non-normal distribution (p < 0.05), the gnls function in R was applied to establish the generalized nonlinear (GNL) models of BD growth, H growth, and AGB. As a model that describes changes in biometric variables (BD and H) with age, the TG equation reflects the basic law of shrub stand growth with certain assumptions based on biological characteristics. Differential or calculus equations of the total growth curve of the shrub were established, and the initial or boundary conditions were substituted to obtain the special solution of the differential equations. Currently, the most widely used TG equations are the five aforementioned functions (Table 3).

Equations Formulas
Gompertz Note: A is the shrub age; I is the increments of BD and H at A year; a 1 , b 1 , and c 1 are model parameters for BD growth; and a 2 , b 2 , and c 2 are model parameters for H growth.

Allometric Growth Equation
The AG equation is used to estimate shrub biomass with easily measurable factors, such as BD and/or H, which saves manpower and material resources and is less destructive when compared to the direct harvesting method [37]. Three AG equations of AGB were established in the study (Table 4).

Equations Independent Variables Model Parameters
Note: W is the AGB; n is the number of stems, and BD m , and H m are the means of BD and H for an individual shrub, respectively; and a 3 and b 3 are model parameters.

Nonlinear Simultaneous Equations
The TG equations of BD and H with age were established using the five functions ( Table 3). The AG equations of AGB were constructed using BD and H as the independent variables and the power function (PF) ( Table 4). Then, the best equations for the TG and AG equations were selected after evaluation. Finally, the nonlinear ordinary least square (NOLS), nonlinear seemingly uncorrelated regression (NSUR), two-stage least squares (2SLS), and three-stage least squares (3SLS) methods were used to build a dynamic system of the NLS equations for AGB by combining the best TG and AG equations ( Table 5). The formula of the NLS equations in this study is as follows: where f are the best performed equations for BD growth, H growth, and AGB models according to the evaluation indices. Table 5. The four methods of the NLS system.

Method
Instruments Objective Function Covariance of θ Note: γ is a column vector of the residuals for each equation, S is the variance-covariance matrix between the equations ( σ ij = e i e j / (T − k i ) · T − k j ), X is the matrix of the partial derivatives with respect to the coefficients, V is the matrix of the instrument variables Z(Z Z) −1 Z, Z is the matrix of the instrument variables, and I m is the n×n identity matrix. The NSUR and 3SLS methods require two solutions. The first solution for NSUR is an NOLS solution that obtains the variance-covariance matrix and fits all the equations simultaneously. The second solution for 3SLS uses the variance-covariance matrix from the 2SLS solution and fits all the equations simultaneously.
NOLS is widely used to estimate parameters in a single equation and can thus lead to biased and inconsistent estimations by neglecting correlations among different equations when estimating parameters for each equation [38]. NSUR is an efficient parameter estimation technique when the error components of a system of unrelated equations are correlated [39]. The parameter variances estimated from large samples using NSUR may be less than the parameter variances obtained by NOLS. Meanwhile, 2SLS and 3SLS are adopted for exact identification and overidentified structural equation models. However, 2SLS neglects correlations, and there is some information loss among different equations for one-by-one estimations [40]. As a system method of estimation, 3SLS is effective for the simultaneous estimation of all parameters of different equations, and it considers the correlation of random errors. A system of NLS equations can be written as follows: where ε t is the residuals of y observations and the function evaluated by the coefficient estimates.

Heteroscedasticity Test and Elimination
The existence of heteroscedasticity in models makes parameter estimation by ordinary regression biased, leading to increased error and parameter variation coefficients, and thereby affecting the reliability of the models. Alternatively, residual graph analysis is an intuitive and convenient analysis method. When the regression model has heteroscedasticity, the distribution of the points on the residual graph shows a certain trend. In this study, when there was heteroscedasticity, the variance functions, including the PF, exponential function (EF), and constant plus PF (CPF), were used to deal with variance heterogeneity. To compare the models with and without variance functions, the Akaike information criterion (AIC), Bayesian information criterion (BIC), and negative two times the logarithm of the likelihood (−2LL) were used; the smaller the values, the better the model. The formulas of the criterions are as follows: where L is the likelihood of the equations, r is the numbers of model parameters, and n i is the number of samplers.

Evaluation of the Models
The TG and AG equations were fit by nonlinear regression. Then, a set of NLS equations combining the TG and AG equations were estimated to obtain the AGB estimation model of C. korshinskii. Finally, the variance functions were used to eliminate heteroscedasticity in the fitting. Leave-one-out cross-validation was performed to evaluate two statistical criteria. Three indices, the AMR, root mean square error (RMSE), and adjusted coefficient of determination (adj-R 2 ), were calculated as model accuracy metrics. The formulas of the fit statistics are as follows: where y ij is BD growth, H growth, or AGB prediction, and y is the average of the observations.

Theoretical Growth Equations of Basal Diameter and Height
Five TG equations were constructed to fit the BD and H growth of C. korshinskii ( Table 6). The results showed that the five equations had good performances when the adj-R 2 was >0.8. However, Gompertz growth equations showed better fitting accuracy than others for both BD and H growth. The AMR, RMSE, and adj-R 2 were 1.4767, 1.9767, and 0.8636 for the BD growth model, respectively. The AMR, RMSE, and adj-R 2 were 0.2378, 0.2939, and 0.8248 for the H growth model, respectively. Therefore, the Gompertz equation was selected as the optimal TG equation and built in the follow-up study. On the basis of easily measurable factors, the total BD (nBD m ), total shrub H (nH m ), and total BD squared multiplied by H ((nBD m ) 2 H m ) were used to establish the AG equations for the AGB of C. korshinskii ( Table 7 was selected as the best model for the AG equation. Table 7. Summary of evaluations from fitting AGB of C. korshinskii shrubs using GNL regression.

Nonlinear Simultaneous Equations of Dynamic Models
In this study, the best performed TG equations of BD and H and the best performed AG equation of AGB were combined using the NOLS, NSUR, 2SLS, and 3SLS methods to construct the NLS equations. Table 8 shows the parameters and evaluation indices after fitting each model. Results revealed that the same estimation parameters and evaluations were obtained by using GNL regression and the NOLS method, which was the same estimation when using the 2SLS and 3SLS methods. However, the NSUR and 2SLS or 3SLS methods showed different outcomes from NOLS. Namely, the NSUR method had better performance than the others. For BD growth of C. korshinskii shrubs, the AMR and RMSE of BD NSUR were 1.3481 and 1.8046, which decreased by 9.54% and 9.53% when compared to Equation (4), respectively. The adj-R 2 of BD NSUR was 0.9055 and increased by 4.63% when compared to Equation (4). For H growth of C. korshinskii shrubs, the AMR and RMSE of H NSUR were 0.2171 and 0.2683, which decreased by 9.53% and 9.54% when compared to Equation (5), respectively. The adj-R 2 of H NSUR was 0.8644 and increased by 4.58% when compared to Equation (5). For the AGB of C. korshinskii shrubs, the AMR and RMSE of AGB NSUR were 1.0514 and 1.3584, which decreased by 4.49% and 4.48% when compared to Equation (6), respectively. The adj-R 2 of AGB NSUR was 0.8246 and increased by 2.82% when compared to Equation (6). Therefore, the NSUR method of the NLS equation system was selected to establish the NLS equations for BD growth, H growth, and AGB.  Figure 2 shows the relationships between the observed and fitted values of BD growth, H growth, and AGB models of C. korshinskii shrubs according to NLS. The rates of R 2 for the three models were 98.66%, 98.65, and 96.14%.

Models with Heteroscedasticity
Biomass and sample plot data often have measurement errors, as well as individual differences among samples, which models sometimes neglect, thereby leading to heteroscedasticity of the estimated results. In this study, the PF, EF, and CPF were used to eliminate heteroscedasticity. Then, the AIC, BIC, and −2LL were used to compare model performances. The value of −2LL was produced by the log of the LR, which provides a ratio of the probability of correct predictions and incorrect predictions. Results revealed that there were significant differences (p < 0.0001) detected between the models with and without heteroscedasticity for BD growth, H growth, and AGB ( Table 9). The models that used PF as the variance function had the smallest AIC, BIC, and −2LL for BD growth, H growth, and AGB.

Models with Heteroscedasticity
Biomass and sample plot data often have measurement errors, as well as individual differences among samples, which models sometimes neglect, thereby leading to heteroscedasticity of the estimated results. In this study, the PF, EF, and CPF were used to eliminate heteroscedasticity. Then, the AIC, BIC, and −2LL were used to compare model performances. The value of −2LL was produced by the log of the LR, which provides a ratio of the probability of correct predictions and incorrect predictions. Results revealed that there were significant differences (p < 0.0001) detected between the models with and without heteroscedasticity for BD growth, H growth, and AGB ( Table 9). The models that used PF as the variance function had the smallest AIC, BIC, and −2LL for BD growth, H growth, and AGB. Note: a LR was calculated with respect to the models without heteroscedasticity and the models using PF as the variance function for the AGB of C. korshinskii shrubs. b LR was calculated with respect to the models without heteroscedasticity and the models using EF as the variance function for the AGB of C. korshinskii shrubs. c LR was calculated with respect to the models without heteroscedasticity and the models using CPF as the variance function for the AGB of C. korshinskii shrubs. Figure 3 shows the standardized residuals of the BD growth, H growth, and AGB models, which had poor fit for all shrub components with obvious trends (trumpet residuals). However, the plots of the models with heteroscedasticity using PF as the variance function had better fit and no obvious trends. Table 10 shows the evaluations of the NLS equation of BD and H growth and AGB of C. korshinskii shrubs. The Gompertz equation indicated there was a relationship between the AGB and age and numbers of stems using GNL regression. Results revealed that the AMR and RMSE of the NLS equation were 1.0514 and 1.3584, which decreased by 22.31% and 19.02%, respectively, when compared to the Gompertz growth equation for AGB. The adj-R 2 was 0.8246 from the NLS equations, which increased by 10.64%. Additionally, we found that the NLS equation greatly improved the model accuracy and reduced forecast estimation errors.  Figure 4 shows the predicted BD growth, H growth, and AGB of C. korshinskii shrubs from the NLS equations and predicted AGB versus age and number of stems from the Gompertz equations using GNL regression. There were significant differences between the AGB estimates from the NLS equations and the Gompertz growth equation using GNL regression. Figure 3 shows the standardized residuals of the BD growth, H growth, and A models, which had poor fit for all shrub components with obvious trends (trumpet re uals). However, the plots of the models with heteroscedasticity using PF as the vari function had better fit and no obvious trends.  Table 10 shows the evaluations of the NLS equation of BD and H growth and A of C. korshinskii shrubs. The Gompertz equation indicated there was a relationship tween the AGB and age and numbers of stems using GNL regression. Results reve The NLS equations system provided a feasible research method for estimating Caragana shrub layer biomass; scientifically supported regional ecosystem protection, restoration, and reconstruction; and accurately evaluated ecosystem service functions in the AASRs.

Models Comparisons
Gompertz equations using GNL regression. There were significant differences between the AGB estimates from the NLS equations and the Gompertz growth equation using GNL regression. The NLS equations system provided a feasible research method for estimating Caragana shrub layer biomass; scientifically supported regional ecosystem protection, restoration, and reconstruction; and accurately evaluated ecosystem service functions in the AASRs.

Discussion
As important components of forest resources, shrubs play unique roles in mitigating global climate change [41]. Monitoring shrub biomass and carbon storage has increasingly

Discussion
As important components of forest resources, shrubs play unique roles in mitigating global climate change [41]. Monitoring shrub biomass and carbon storage has increasingly attracted global attention [3,42]. There are many direct harvesting, model estimation, and remote sensing interpretation methods for obtaining shrub biomass. Although the accuracy of precise harvesting is very high, it requires considerable manpower and material resources and is destructive to the forest. Additionally, it greatly increases the vulnerability of the desert ecosystem and does not contribute to its protection. Advanced methods, such as remote sensing technology and geographic information systems, have irreplaceable advantages in the estimation of forest biomass and net growth on a large and even global scale. However, due to the limitations of remote sensing data in space, spectrum, and radiation resolution, the accuracy of estimating ground forest biomass using remote sensing data is unreliable [43]. Therefore, due to the establishment of shrub biomass prediction models using easily measurable factors, these models have become effective methods for measuring shrub biomass in the AASRs.
In this study, we applied the GNL regression approach [44] to develop TG equations of BD and H on the basis of the Gompertz, Logistic, Mitscherlich, Richards, and Korf equations. Then, the optimal equations for BD and H growth were selected on the basis of the AMR, RMSE, and adj-R 2 values. The Gompertz equation showed the best performance for both BD and H growth of C. korshinskii shrubs (Table 6) [20,45]. Moreover, the nBD m , nH m , and (nBD m ) 2 H m were used as the independent variables to build the AGB AG equation. The equation using (nBD m ) 2 H m as the independent variable performed better than the others (Table 7) [10,11,46]. Finally, the NOLS, NSUR, 2SLS, and 3SLS methods of the NLS equations were constructed and combined BD and H growth from the TG equation and AGB from the AG equation (Table 8).
In this study, we found identical fit and evaluation when using the NOLS and GNL regression approach. It also showed the same prediction for the methods of 2SLS and 3SLS. Tang et al. [47] demonstrated that, in some situations, the NOLS is more efficient than the GNL regressions. Despite having a slight lower efficiency, a main advantage of 2SLS is that this can be estimated consistently with the existence of heteroscedasticity. What is more, there were theoretical advantages for both 2SLS and 3SLS when compared to the NSUR method. However, they required large sample datasets and overidentified estimation, which was difficult to realize in the forest biomass equations [48]. The duration of the equation fittings for the NLS equations methods was ordered as follows: 3SLS (16.87 s) > 2SLS (10.92 s) > NSUR (3.21 s) > NOLS (1.26 s). Therefore, the NSUR method of the NLS equations was determined to have the best performance [49,50]. The NLS equations system could significantly improve the model accuracy and reduce the estimation error of the predictions [23]. The estimation error also met the accuracy requirements on a regional scale; reduced the error between biomass and carbon storage estimation; and could serve as a reference for forest management, the rational utilization of resources, and mitigation of climate change.
As an important component of shrub biomass, AGB accounts for 60-70% of the biomass of an individual shrub [51], while the underground biomass of C. korshinskii shrubs accounts for 30-40% of the total shrub biomass [52]. However, the variation in underground biomass growth for C. korshinskii shrubs remains unclear. Future studies should focus on identifying variations in shrub biomass to assist with the improvement of shrub biomass estimates when modeling C. korshinskii growth.

Conclusions
In this study, the NLS equation systems of BD growth, H growth, and AGB models were used to quantify C. korshinskii shrub AGB growth on the basis of the theorical growth and allometric growth equations in northwestern China. It was found that the Gompertz equation had better performance than the other equations using GNL regression for the BD and H growth models. For the AGB models, using (nBD m ) 2 H m as the independent variable resulted in better evaluations. The NSUR method of the NLS equations had better performance than the other three methods and the independent equations, which thereby provides an efficient method for predicting AGB growth of C. korshinskii shrubs. Moreover, the equations with heteroscedasticity accounted for multiple sources of heteroscedasticity found in the data, thus making the NLS equations approach an attractive option for C. korshinskii shrub biomass growth estimation. The NLS equations fitted in this study are appropriate for modeling new shrubs that have ecological characteristics and growth features similar to C. korshinskii shrubs in the AASRs.
Author Contributions: Conceptualization, H.X. and X.J.; methodology, H.X.; software, X.J.; validation, H.X. and X.J.; formal analysis, X.J.; investigation, B.W. and X.W.; writing-original draft preparation, X.J.; writing-review and editing, H.X.; visualization, X.J.; supervision, H.X. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: Data available on request due to privacy restrictions. The data presented in this study are available on request from the corresponding author. The data are not publicly available due to be collected by privacy.