Carbon stock and sink economic values of forest ecosystem in the forest industry region of Heilongjiang Province, China

The values of forest carbon stock (CSV) and carbon sink (COV) are important topics in the global carbon cycle. We quantitatively analyzed the factors affecting changes in both for forest ecosystem in 2000−2015. With multiple linear stepwise regression analysis, we obtained the factors that had a significant impact on changes of CSV and COV, and then the impacts of these variables on CSV and COV were used for further quantitative analysis using the vector autoregressive model. Our results indicated that both stand age and afforestation area positively affect CSV and COV; however, the forest enterprise gross output value negatively affects CSV. Stand age has the largest long-term cumulative impact on CSV and COV, reaching 40.4% and 9.8%, respectively. The impact of enterprise gross output value and afforestation area on CSV and COV is the smallest, reaching 4.0% and 0.3%, respectively.


Introduction
Climate change stimulated by emissions of CO 2 has become a significant issue worldwide. As the largest carbon pool in terrestrial ecosystems, forests play an important role in maintaining global carbon balances and in balancing global warming (Dixon et al. 1994;Wingfield et al. 2015;Cui et al. 2018). Trees absorb CO 2 and release oxygen through photosynthesis in the growth process, and the absorbed CO 2 is fixed in the plant and soil in the form of biomass, making forests the most important carbon sink or carbon pool of terrestrial ecosystems. Globally, forest vegetation and soils contain approximately 1146 MT carbon, accounting for about 45% of the total earth's carbon (Bonan 2008). Carbon is also significantly exchanged between forests and the atmosphere. Carbon storage (i.e., the stock of carbon stored in plant biomass and soil) and sequestration (i.e., the removal of atmospheric carbon per unit of time by plants and soils) are recognized ecosystem services that contribute to climate change mitigation (Pan et al. 2011;Alongi 2012;Caldeira 2012). Therefore, quantifying economic values of total carbon stocks and annual carbon sinks of forest ecosystem are required for carbon trading and forest carbon projects.
Total forest carbon stock mainly includes vegetation, soil and litter carbon stocks, an important parameter for assessing carbon sequestration capacity and carbon budgets (Ravindranath and Ostwald 2009). The value of forest carbon stock (CSV) is to express the forest fixed amount of CO 2 in the monetary terms which helps people to understand the utility of the forest carbon sink function ). In addition, the value of forest carbon sinks (COV) is often used to assess the value of a forest's ecological services. As an important ecological function, carbon fixation and oxygen release by forest vegetation play an important regulatory role in the material cycle and energy flow of ecosystems by providing and increasing O 2 concentrations while fixing and reducing atmospheric CO 2 , and maintaining the balance between atmospheric CO 2 and O 2 (Deng et al. 2011). Therefore, estimating COV in forest vegetation has attracted considerable attention (Zhang et al. 2018). However, there is a common challenge to understanding the factors affecting CSV and COV.
Several studies have analyzed the factors affecting the economic value of forest carbon stock, CO 2 fixation and O 2 release. Adams et al. (1999) indicated that the most important factor affecting the value of forest carbon sinks in the short-term is forest management; meanwhile, forest area is the decisive factor affecting long-term forest carbon sinks. Xu et al. (2010) found that timber production was the most important factor affecting the value of forest carbon sinks in China. However, despite the need for such research, there is a lack of emphasis on identifying the environment, policy, and economic factors affecting CSV and COV.
The forest industry region of Heilongjiang Province (HLJFI) is located in Northeast China ( Fig. 1), with a total forest area of over 10.0 million ha, which is about 31.0% of the total forest area of China with a forest cover of 83.9%. We quantified the economic values of total forest carbon stock, carbon fixation, and oxygen release during 2000−2015 by using national forest inventory and net primary productivity (NPP) data. In addition, we econometrically assessed the potential effects and importance of environment, policy, and economic factors on the CSV and COV in HLJFI.

Forest inventory data
We estimated forest biomass, carbon stock and carbon density from 2000 to 2010 with three consecutive forest inventories (the sixth to the eighth inventories 2000, 2005, and 2010, respectively) in HLJFI. Meanwhile, the forest biomass, carbon stock and carbon density from 2010 to 2015 was predicted by using the forest inventory data from 2000 to 2010, due to the lack of available forest inventory data in 2015. The data set was measured with 5-year intervals and periodically to provide statistics on stand age, height (H), diameter at breast height (DBH), and forest type.

Climate data
Annual average temperature and precipitation data were collected from the China Meteorological Data Network, including 26 stations in HLJFI. We calculated the average values of these stations as the climate factors.

Economic and policy data
The economic and management data were developed by the Heilongjiang forest industry. These data sets record economic factors, including forestry investment, number of employees, enterprise gross output value, and management factors such as forest tending area, afforestation area (AA), and reforestation area from 2000 to 2015.

NPP data
The annual NPP from 2000 to 2015 in HLJFI was produced by Wang et al. (2017) using the Integrated Terrestrial Ecosystem Carbon (InTEC) model. The InTEC model has been validated and applied widely (Ju et al. 2007;Zhang et al. 2012Zhang et al. , 2015.

Forest fire disturbance
The annual forest fire disturbance area in HLJFI from 2000 to 2015 was collected from Heilongjiang Provincial Forestry Department, the data include various information such as geographic coordinates of the fire location, fire-discovery time, fire-extinguishing time, burning area, and fire cause.

Calculating the carbon stock from 2000 to 2010
Total forest carbon stock mainly includes vegetation, soil and understory carbon stocks. For vegetation carbon stock, the volume-derived biomass method is effective to estimate forest carbon storage; therefore, we first calculated the vegetation biomass (includes above-ground and under-ground) using the relationship between stand biomass and stocking volumes for the main tree species in Heilongjiang Province, which was produced by Dong et al. (2015). A ratio of 0.5 was then used to convert biomass to carbon density. We used vegetation carbon storage for soil and understory carbon stock combined with conversion coefficients for calculation, as Eqs. 1 and 2.
where C si represents the soil carbon density of ith forest type, C ui represents the understory carbon density of ith forest type, C vi is the vegetation carbon density of ith forest type, α is the carbon conversion coefficient of vegetation to understory vegetation, and β is the carbon conversion coefficient of vegetation to the soil. α and β are 0.195 and 1.244, respectively (IPCC, 1996).
The total carbon stock of the forest ecosystem in 2000, 2005 and 2010 is calculated as the sum of vegetation carbon stock, soil carbon stock, and understory carbon stock (Eq. 3).
where C i represents the total carbon stock of ith forest type, and A i is the area of ith forest type.

Predicting the carbon stock in 2015
Due to the lack of available forest inventory data, an effective method by Zhang et al. (2020) was used to predict the total forest carbon stock in 2015 by combining an expanded forest area and the average annual growth per unit area during 2000−2010 (Eqs. 4−6). Finally, assume that the increase rate of forest carbon stocks is constant during each 5-year interval period of forest inventory, and calculate annual total forest carbon stocks from 2000 to 2015. (1) where C 2015 represents predicted carbon stock in 2015 (Tg); C R is the predicted carbon stock due to an average annual increase in 2015 (Tg); C S is the predicted carbon stock due to increased forest area (Tg);

Calculating CSV and COV
In this study, the carbon taxation method is solely employed to estimate annual total carbon stock price, CO 2 fixation, and O 2 release. We adopted the Swedish carbon tax rate (150 USD t −1 ) as the unit price of forest carbon sequestration from the atmosphere, and converted it into yuan RMB according to the current exchange rate. The price of forest oxygen release is 1,000 yuan t −1 (Nation Forestry and Grassland Administration 2020). The annual CSV can be calculated as the unit price of forest carbon sequestration multiplied by the annual forest carbon stock from 2000 to 2015.
We defined the annual COV as the sum of the price of forest carbon sequestration from the atmosphere and oxygen released into the atmosphere in unit time (Eq. 7) (Nation Forestry and Grassland Administration 2020).
where V C and V O represent the annual value of forest carbon fixation and oxygen release, respectively; P C and P O represent the unit price of forest carbon fixation and oxygen release, respectively; M C is the carbon content of CO 2 (27.27%); F is the annual soil carbon fixation (F/NPP = 0.02/0.49) (Woodburry et al. 2007); and S is the area of the ith forest type.

Evaluating influencing factors
We used multiple linear stepwise regression analysis (MLR) combined with 11 variables to obtain the factors that have a significant impact on the changes to CSV and COV; the variables include stand age, annual average temperature and precipitation, annual forest fire disturbance area and annual economic and policy factors. The impact of these factors on CSV and COV was then used for further quantitative analysis.
We employed the vector autoregressive model (VAR) to measure contribution rate of variables filtered by MLR on CSV and COV. The VAR is an economic model used to capture the linear interdependencies among multiple time series, and it can predict more accurately than traditional structural models (Sims 1980). In the VAR model, all variables enter the model in the same way; each variable has an equation that explains its evolution according to its past values and error term of all the variables in the model, as it describes the dynamic evolution of several variables from their shared history. VAR models have been widely used in determining and analyzing the relationship between multiple time series economic indicators (Amuka et al. 2016;Uwamariya and Gasana 2018). The structure of the VAR model is as follows: where Y t is a h vector of endogenous variables, X t is a p vector of exogenous variables, k is the lag orders, A 1 ,…, A k are matrices of coefficients to be estimated, and ε t is a k vector of innovations.

CSV and COV during 2000−2015
The total forest carbon stock increased from 715.5 Tg in 2000 to 1041.0 Tg in 2010, and we predicted that the total forest carbon stock would continue to increase to 1240.9 Tg in 2015.
The CSV increased from 858.6 billion yuan RMB in 2000 to 1489.0 billion in 2015, an average annual growth of nearly 5.0%. The COV increased from 51.5 billion yuan RMB in 2000 to 60.7 billion in 2015, an average annual rate of growth of nearly 1.0% (Fig. 2).
The results of MLR indicated that the variables that have a significant impact on the CSV are stand age (AGE), enterprise gross output value (EGOV), and afforestation area (AA). The variables that have a significant impact on the COV are stand age and the afforestation area. Both models are significant (p < 0.05), and the R 2 of the CSV model and the COV model are 97.0% and 83.0%, respectively.

VAR model result
We used CSV and COV as response variables and the factors fitted by the MLR method as independent variables to fit two VAR models separately. To eliminate heteroscedasticity, we chose the method of taking the natural logarithm to each variable that is commonly used in the time series in the normalization method.

Variable stability test
We first tested the stability of the original variable series and found that the original series was not stable. Therefore, we made a further difference series test. Consequently, by repeatedly testing the difference series of these dependent variables and independent variables, all were stable for the first-order difference series (Table 1).

Model stability test
The reliability of the VAR model result depends on how the variables in the model can be identified over time. The stability condition is satisfied in VAR only if all the characteristic roots lie within the unit circle. The result of the VAR stability test of the CSV-VAR model and COV-VAR model are presented in Fig. 3a, b, showing that all the values of the two models lie within the unit circle, thereby satisfying the VAR stability condition.

VAR lag order selection criteria
Before the co-integration test, it is crucial to determine the lag period because the optimal lag of the co-integration test is generally the optimal lag of VAR minus 1. Therefore, first we need to determine the optimal lag order of the VAR model. Tables 1 and 2 show the five evaluation statistical indicators of the 0−2 order VAR model: the values of LR, FPE, AIC, SC, and HQ, and the optimal lag period given by each evaluation indicator is marked with "*". If the test results of these statistical evaluation indicators are inconsistent, the lag order of VAR is determined based on the majority principle. It can be seen that all the indicators select the optimal lag order 2 in Tables 2 and 3. Consequently, the lag order of the CSV-VAR model and COV-VAR model are all defined as 2.

Johansen co-integration test
When conducting time series analysis, it is traditionally required that the series should be stable, that is, there is no random trend or definite trend, otherwise it will produce a false regression problem. However, in reality, the time series  is usually non-stationary and can be differentiated to make it stable (such as the first-order difference processing of each data before), but this will lose the total long-term information, and this information is necessary for the analysis of the problem. Therefore, co-integration was used to solve this problem, and finally to determine the long-term stable relationship between various variables. Tables 4 and 5 show the trace statistics results of the two models, respectively. There is only one co-integration equation for the CSV-VAR model and COV-VAR model at the 0.05 significance level.

Impulse response function test
The impulse response is an important aspect of the dynamic characteristics of the VAR model system. It describes the trajectory of each variable's change or impact on itself and on all other variables, and helps us to know the pattern of response to shocks, that is, whether the response is in an upward or downward movement. The shock from the exogenous variable can transmit to the endogenous; therefore, impulse response provides information on the period-toperiod response of the endogenous variable to the innovation/shocks in the exogenous.
For CSV, Fig. 4 indicates that, except for the slight negative impulse in the seventh and eighth period, the stand age factor had a positive impulse in other periods, and reached its maximum in the third period. However, the effect of stand age on carbon prices began to gradually decrease after the fourth period. The EGOV was always a negative impulse on CSV (Fig. 4), and after the ninth period, the impact of EGOV on CSV is almost 0. The impact of AA on CSV showed a positive trend over the whole period, and reached its maximum in the fourth period. Figure 5 shows that the impact of stand age on COV had an obvious positive trend before the fourth period, but turned into a negative trend during the fourth and sixth periods. After that, the impact of stand age on COV was close to the zero axis. Figure 5 also indicates that AA had a negative impact on COV during the first two periods. In the subsequent period, the impact of AA on COV was almost positive.   From the impulse response function perspective, each influencing factor, positively or negatively, impacts CSV or COV. To further quantitatively analyze the degree of contribution of each influencing factor in the long-term CSV or COV change process, the variance decomposition needs to be conducted.
The changes of CSV affected by itself showed a downward trend, from 100% at the beginning to 49.5% in the last period. In contrast, the impact of other variables on the CSV showed a gradual upward trend. The influence of stand age, EGOV and AA on CSV reached 40.4%, 9.7% and 0.3% in the tenth period, respectively (Table 6). Table 7 shows that the change of COV was gradually decreasing due to its own impact, from 100% in the first period to 86.2% in the tenth period. The influence of stand age on the COV showed an overall upward trend; in the long run, the cumulative contribution of age to COV changes reached 9.7%. Similarly, the impact of AA on the COV also showed an upward trend, but the contribution was relatively smaller, reaching 4.0% in the tenth period.

Discussion
Numerous studies have shown that stand age is one of the most important factors that affects carbon accumulation in forest ecosystem. Our impulse response function test result showed that stand age had a positive impulse on CSV. This is because, with the increase of stand age, biomass carbon stock continues to accumulate. In our study, soil carbon storage and understory carbon storage are directly proportional to vegetation carbon storage; hence, total forest ecosystem carbon stock will increase with stand age, which affects the CSV.
Similarly, stand age significantly affects changes in net primary productivity (NPP). Wang et al. (2018a) reported that NPP increased rapidly during the early stages of stand growth, reached a peak when mature and then declined after maturity. During a young forest's growth, NPP increases rapidly with age caused mainly by increases in biomass growth due to faster cell division. For old forest stands, where photosynthetic products were mainly consumed by foliage and fine root turnover, NPP begins to decrease. Our study used NPP and the carbon tax method to calculate COV, and the change curve of the impact of stand age on COV is consistent with the change curve of NPP-age; therefore, forest age is an important influencing factor of COV change.
The impact of AA on CSV is generally positive. The effect of the first two periods is weak but increases rapidly from the third period and reaches a peak in the fourth period. Our result results align well with Wu et al. (2015), who also indicated a lag in the impact of AA on carbon stock; the initial effect is relatively small. However, it also continues to increase to a maximum as AA and forest biomass continue to increase.
The increase in AA has a slightly positive impact on COV in the first four periods, then the impulse curve fluctuates up and down around the 0 axis. Huang et al. (2015) indicated that during the first four years of afforestation, NPP increased rapidly and then the growth rate decreased. Afforestation can sometimes negatively impact NPP because afforestation reduces biodiversity and soil fertility, further affects the local water cycle, and increases the risk of forest fires, pests, and other ecological problems (Thuille and Schulze 2006).
Our results indicate that EGOV is an important factor affecting the carbon value of forests and always has a negative impact. This is agreement with Xu and Jia (2012) which showed, as EGOV increases, forest carbon stock will decrease, which in turn affects CSV. Wang et al. (2018b) also considered that the EGOV by excessive consumption of forest resources is not conducive to forest carbon sinks.
Due to data limitations, our research only analyzed the impact of environmental, economic, and policy factors on the CSV and COV over 15 years. Despite the limitation, our study still highlights the main factors and their influence on CSV and COV. Our results show that stand age and afforestation will positively impact CSC and COV, which indicates that large-area afforestation with young trees is of considerable significance for increasing CSV and COV. Meanwhile, our results also indicate that a reduction in forest resource consumption is contributes to the carbon sink of forest ecosystem.

Conclusion
We used multiple linear regression to analyze the factors that significantly impact the values of forest carbon stock and carbon sink in HLJFI from 2000 to 2015, and adopted the VAR model to further quantitatively analyze the impact of these factors on the values of forest carbon stock and carbon sink changes. Among environmental, economic, and policy factors, stand age and afforestation area have significant effects on the values of forest carbon stock and carbon sink, while enterprise gross output value has a significant impact on the values of forest carbon stock. For the values of forest carbon stock changes, the increase in stand age and afforestation area have a promotional effect on the values of forest carbon stock; however, the enterprise gross output value will negatively affect the values of forest carbon stock in the long-term. For the values of forest carbon sink changes, the impact of stand age and afforestation area is relatively weak, and younger and more productive forest stands are conducive to the values of forest carbon sink growth.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.