Relationship between Net Primary Productivity and Forest Stand Age under Different Site Conditions and Its Implications for Regional Carbon Cycle Study

Abstract

Net primary productivity (NPP) is a key component in the terrestrial ecosystem carbon cycle, and it varies according to stand age and site class index (SCI) for different forest types. Here we report an improved method for describing the relationships between NPP and stand age at various SCI values for the main forest types and groups in Heilongjiang Province, China, using existing yield tables, biomass equations, and forest inventory data. We calculated NPP as the sum of four components: Annual accumulation of live biomass, annual mortality of biomass, foliage turnover, and fine root turnover in soil. We also consider the NPP of understory vegetation or moss. These NPP-age relationships under different site conditions indicate that the NPP values of broadleaved and coniferous, as well as broadleaved mixed forests increase rapidly and reach a maximum when in young forests. However, for coniferous forest types, the maximum NPP generally occurs in mature forests. In addition, a higher SCI leads to a higher NPP value. Finally, we input these NPP-age relationships at various SCI values into the Integrated Terrestrial Ecosystem Carbon (InTEC) model to modify NPP modeling to estimate NPP in Heilongjiang Province in China from 2001 to 2010. All of the results showed that the methods reported in this study provide a reliable approach for estimating regional forest carbon budgets.

1. Introduction

Net primary production (NPP) not only reflects biological properties and stand structures for a forest ecosystem, but also describes the influence of environmental factors on forest growth [1]. Growth is substantially affected by CO₂ fertilization, climate variability, nitrogen (N) deposition [2,3], and disturbance factors, such as wildfires, harvesting, and insects [4,5,6]. Disturbances are major determinants of forest carbon stocks and uptake. They generally reduce land carbon stocks but also initiate a regrowth legacy that contributes substantially to the contemporary rate of carbon stock increase in forestlands [7].

The interannual variability of NPP is not only affected by climatic variability, tree species composition, and forest succession [8], but is also closely related to forest age. Many previous studies have demonstrated that stand age impacts on NPP can be substantial at both local and national scales [9,10,11]. For example, NPP rapidly increases in young forest stands, reaches a maximum when mature, and gradually declines in older stands [12,13,14]. However, several recent studies have provided evidence that not all forest types exhibit this decline [15,16]. Additionally, NPP is not only affected by forest age, but also by factors related to environmental site conditions. Chen et al. [17] showed that NPP increased faster and the value of NPP was higher at more productive sites. Hence, both quantitative analyses of the successional change of NPP along with age structure and an understanding of the relationships between NPP, age, and site conditions, are critically important for improving forest carbon cycle estimation.

Traditional forest inventories have provided a large amount of ground data on tree growth, including diameter at breast height (DBH) and tree height (H). Fang et al. [18] used forest inventory data from 1949 to 1998 to estimate forest biomass and carbon (C) storage in China. Similarly, stand yield tables were developed and designed based on forest inventories, and contain the most readily available information regarding the effects of stand age on tree growth. Tree growth can be combined with stand yield tables to estimate the relationships between NPP and age for various site conditions. Due to averaging and smoothing procedures which are generally applied in yield table development, NPP-age relationships are averaged over long-term mean environmental conditions. Therefore, this method can produce mean NPP-age relationships [17].

Biogeochemical process models are effective tools used to help understand the response of forest ecosystems to various environmental drivers. The Integrated Terrestrial Ecosystem Carbon (InTEC) model is a process-based biogeochemical model produced by Chen et al. [19] which integrates the effects of both disturbance and non-disturbance factors in long-term C budget simulations [20,21]. The model has been validated and used to simulate the historical carbon budget of forest environments. For instance, it was used to simulate the historical change of C dynamics and to analyze spatiotemporal carbon sink distributions in Canada [22]. Then, the model was calibrated further for China [3,23]. Recently, the InTEC model has been used to estimate the C balance and distribution in North American forest ecosystems [5,6].

Studies have used a variety of methods to develop the NPP-age relationships for different forest types and regions. Chen et al. [17] used yield tables to quantify mean age-NPP relationships in black spruce stands in Ontario, Canada. Wang et al. [24] developed functions describing the relationship between national mean NPP and stand age for five typical forest ecosystems in China in 2001, using stand age information derived from forest inventory data and NPP simulated by BEPS (Boreal Ecosystem Productivity Simulator). Zaehle et al. [25] analyzed the first-order effects of age on aboveground productivity of forests in Europe using a global terrestrial biosphere model (the Lund-Potsdam-Jena model, LPJ). He et al. [26] developed NPP-age relationships for major forest biomes in the U.S. by combining data from the Forest Inventory and Analysis (FIA), remote sensing, and species-specific traits.

Until now, none of the existing regional C cycle modeling has incorporated systematic NPP-age relationships at various site class index (SCI) values in Heilongjiang Province, China. Heilongjiang Province (121°11′ E to 135°05′ E and 43°25′ N to 53°33′ N) is located in Northeast China (Figure 1). The province has a total forest area of over 20 million ha, which is about 45.8% of the total land area of the province, and the forest coverage rate is 45.7%. The total volume of living trees is 1.76 billion m³ and the average volume per hectare is 78.6 m³. In addition, the forest fire hazard in Heilongjiang Province is the highest of any area in China [27]. Therefore, the objectives of this study were to develop systematic NPP-age relationships at various SCI values for major forest types and groups in Heilongjiang Province, China, by combining data from forest inventories and stand yield tables, and to examine whether the new NPP-age relationships at various SCI values could estimate NPP through C cycle modeling.

2. Data Set

2.1. Forest Inventory Data and Yield Tables

The Forest Inventory data set was derived from the eighth national forest inventory data recorded from 2009 to 2014 in Heilongjiang Province [28]. The data set was measured periodically to provide statistics on stand age, H, DBH, and forest type.

Yield tables reflect the growth processes of forest stands with various stand structures and forest treatments. The yield tables of Heilongjiang Province, China were developed by the Forest Survey Scheme Designing Institute in 2010, and they provided information on the forest, such as forest type, mean age, DBH, H, stand density (S), volume (V) for different forest types and groups, and SCI. Based on this information, the standing live biomass can be calculated. A five-year age interval is used in the yield tables.

2.2. Input Data in the InTEC Model

The Integrated Terrestrial Ecosystem Carbon (InTEC) model is a process-based biogeochemical model driven by monthly climate data, vegetation parameters, and forest disturbance information. The model can estimate the effects of changing climate, atmospheric composition, disturbance, and forest recovery on the long-term C and N cycles in forest ecosystems.

To drive the InTEC model, a series of datasets were created in this study. All spatial data are listed in

Table 1. Descriptions of input data used in the InTEC model.

Data Set	Description	Resolution
Climate data	Monthly mean temperature, relative humidity, total precipitation and solar irradiance during 1901–2010.	1 km
CO2 concentrations	Measured data between 1958 and 2010 and Canadian Global Coupled Model (CGCM2) before 1958 [34]	1 km
Nitrogen deposition	Chemical transport model [22]	1 km
Forest stand age map	Developed from forest inventory data in 2010	1 km
Forest type map	Developed from CAS In 2006	1 km
Soil data	Wilt point and field capacity of soil water, soil depth, fraction of clay, silt, and sand	1 km
Leaf area index	Developed from SPOT4 VEGETATION data in 2003	1 km
Reference NPP	Developed from BEPS model in 2003	1 km
SCI map	Developed from forest inventory data in 2010	1 km

and employed in the UTM projection WGS-84 coordinate system and interpolated to a one kilometer resolution.

The monthly mean temperature, relative humidity, and total precipitation for Heilongjiang in China from 1901 to 2010 were obtained from 0.5° global data set interpolated by the UK Climate Research Unit [29]. The dataset was produced by measurements taken at available stations of the National Meteorological Administration in China. Monthly solar irradiance data for the period prior to 1948 was produced by the Bristow-Campbell model [30] derived from the historical temperature, humidity, and precipitation data. For the period from 1948 to 2010, monthly solar irradiance data from the T62 Gaussian reanalysis data of the U.S. National Center for Atmospheric Research (NCAR) were used [31].

The annual atmospheric CO₂ concentrations from 1958 to 2010 were taken from the dataset obtained at the Mauna Loa Observatory (20° N, 156° W) [32]. The pre-1958 CO₂ concentrations were estimated based on the Canadian Global Coupled Model (CGCM2) [33].

Spatial nitrogen (N) deposition data in 1993 were obtained from a dataset that was simulated based on a chemical transport model (TM3). This dataset included the predicted value of N deposition in 1860, 1993, and 2050, and the spatial resolution was 5° in longitude and 3.75° in latitude. Spatial N deposition from 1901 to 2010 (except 1993) was calculated based on historical greenhouse gas emissions and the average N deposition data in 1993 [21].

A map of forest stand age in 2010 (Figure 2) was produced from the eighth national forest inventory data recorded from 2009 to 2014 in Heilongjiang Province [28]. Stand age for each pixel was estimated by interpolation (kriging) [34].

The forest type map of Heilongjiang at a one km resolution (Figure 3) was obtained from the Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences, published in 2007 [35]. The scale of the map is 1:1,000,000.

The physical properties of the soil, including the field capacity of soil water, soil depth, wilting point, and the fractions of clay, silt, and sand, were included. Field capacity and wilting point were derived from the International Geosphere-Biosphere Programme, Global Gridded Surfaces of Selected Soil Characteristics [36]. Soil depth was derived from the global soil texture dataset from Oak Ridge National Laboratory Distributed Active Archive Center, Tennessee, U.S. [37]. The fractions of clay, silt, and sand were obtained from the Harmonized World Soil Database (HWSD) constructed by the Food and Agriculture Organization of the United Nations (FAO) and the International Institute for Applied Systems Analysis (IIASA) [38].

A maximum leaf area index (LAI) map was produced in 2003 by Deng et al. [39] using SPOT-VEGETATION data based on an algorithm that considers the bidirectional distribution reflectance function (BRDF). Detailed descriptions of the 2003 LAI map and data processing methods are reported elsewhere [39].

NPP in 2003 was selected as the reference year. Wang et al. [40] produced a NPP distribution map in 2003 using the BEPS model, and the reference year NPP data were used to calibrate the initial value of NPP in InTEC.

Because SCI was defined as the average height of a forest stand at a predetermined age, the national forest inventory data can offer both age and height information for producing SCI distribution maps in Heilongjiang Province (Figure 4). The minimum SCI is eight, the maximum is 22, and the average is 14. We assume that SCI is unchanged during the model simulation period.

3. Methods

3.1. SCI and NPP Calculation

SCI is defined as the mean tree height in a specific base stand age (20–50 years old) of a site. It can be calculated using the following two equations which were published in yield tables: Equation (1) is for Pinus koraiensis and Pinus sylvestris forests, and Equation (2) is used for other forest types and groups:

$$SCI=H\exp \left[-k\left(1/t_I-1/t\right)\right]$$

(1)

$$SCI=H\frac{{\left[1-\exp (-k{t}_I)\right]}^c}{{\left[1-\exp (-kt)\right]}^c}$$

(2)

where H represents the mean tree height for a specific forest stand; t_I is the base stand age; t is the current mean age of forest; and k and c are regression coefficients. The t_I, k, and c values are listed in

Table 2. Coefficient estimates for Equations (1) and (2) for nine forest types and three forest type groups. t_I is the base stand age in years, and k and c are regression coefficients.

Forest Type and Group	${\mathit{t}}_{\mathit{I}}$	$\mathit{k}$	$\mathit{c}$
Pinus koraiensis	30	21.468	-
Larix gmelinii	30	0.023	0.837
Pinus sylvestris	30	17.454	-
Quercus mongolica	50	0.012	0.579
Populus davidiana	30	0.032	0.805
Tilia amurensis	50	0.029	0.838
Betula platyphylla	40	0.009	0.663
Betula davurica	40	0.021	0.775
Populus nigra	20	0.041	0.620
Coniferous forest	80	0.004	0.625
Broadleaved forest	40	0.011	0.541
Mixed forest	50	0.007	0.579

and were published in yield tables.

We estimated the NPP of a forest over a year as the sum of four components. The first was change in biomass (∆B_tree), including stems (∆B_s), branches (∆B_b), foliage (∆B_f), and roots (∆B_r). The second component was mortality, which includes standing dead trees and downed dead wood (M_tree). The third component was tissue turnover (L_tree), including the turnover of foliage (L_f) and fine roots (L_fr) in soil. The last component was the NPP of the understory vegetation and moss (NPP_und). Hence, NPP over a year could be calculated using the following equation [17]:

$$NPP=\Delta B_{tree}+M_{tree}+L_{tree}+NP{P}_{und}$$

(3)

3.2. The First Component of NPP—Change in Biomass

The biomass of several forest components can be calculated using the stand information in yield tables and biomass equations. In the single-tree additive system of biomass equations for the main tree species in Heilongjiang Province, that was produced by Dong et al. [41,42], the biomass equations contained only DBH and were constructed based on allometric equations.

For three mixed forest type groups, a stand biomass equation produced by Dong et al. [42] was used to estimate the biomass of the main forest type groups in Heilongjiang Province. These biomass equations contained both tree basal area (G) and H.

Growth equations include the Logistic equation, Mitscherlich equation, Gompertz equation, Richards equation [43], and Korf equation [44]. In this study, we found Korf equation fitted the biomass data well for broadleaved mixed forests and coniferous and broadleaved mixed forests with an R² larger than 0.997. Meanwhile, Richards equation produced better model fitting for the other forest types (R² = 0.995–0.999). The Richards and Korf equations were computed by Equations (4) and (5), respectively. Hence, the change in biomass at age i (∆B(i)) could be calculated by the first-order derivative of the Equations (4) and (5).

$$B=A\times {(1-e^{-bt})}^c$$

(4)

$$B=A\times e^{{-bt}^{-c}}$$

(5)

where B is biomass in t/hm², including B_t, B_s, B_b, B_f, and B_r; t is stand age in years; and a, b, and c are coefficients of tree growth equations that are listed in Appendix A

Table A1. The Regression Coefficients for Equation (6) for nine pure forest types.

Forest Type	SCI	A	b	c	R2	RMSE	NO.
Betula platyphylla	SCI = 10	5236.44	0.02	1.38	99.99	12.36	13
	SCI = 12	6201.10	0.02	1.40	99.99	16.83	13
	SCI = 14	7222.93	0.02	1.43	99.99	16.24	13
	SCI = 16	8366.93	0.02	1.41	99.99	16.89	13
	SCI = 18	8900.00	0.03	1.49	99.97	36.34	13
Tilia amurensis	SCI = 10	5682.84	0.02	1.18	99.96	25.28	15
	SCI = 12	6950.75	0.02	1.17	99.95	38.33	15
	SCI = 14	8238.55	0.02	1.17	99.96	41.40	15
	SCI = 16	9410.55	0.02	1.21	99.97	40.57	15
	SCI = 18	10,702.49	0.02	1.22	99.97	48.42	15
Betula davurica	SCI = 8	2999.99	0.01	1.07	99.81	21.74	13
	SCI = 10	4200.00	0.01	1.10	99.78	36.52	13
	SCI = 12	5699.99	0.01	1.09	99.79	48.89	13
	SCI = 14	7499.97	0.02	1.10	99.79	66.78	13
	SCI = 16	9499.95	0.02	1.07	99.85	71.88	13
Pinus koraiensis	SCI = 6	6035.83	0.09	5.20	99.92	53.30	21
	SCI = 8	6398.17	0.11	5.89	99.94	50.06	21
	SCI = 10	6687.73	0.11	6.33	99.93	58.34	21
	SCI = 12	6957.79	0.12	6.55	99.88	78.06	21
Larix gmelinii	SCI = 8	6287.62	0.05	2.28	99.96	28.06	17
	SCI = 10	7905.53	0.05	2.33	99.96	38.48	17
	SCI = 12	9504.50	0.05	2.30	99.96	46.80	17
	SCI = 14	10,808.28	0.06	2.39	99.97	46.98	17
	SCI = 16	12,222.37	0.06	2.37	99.95	68.27	17
	SCI = 18	13,096.63	0.06	2.51	99.97	53.06	17
	SCI = 20	14,095.76	0.07	2.56	99.97	68.49	17
Populus davidiana	SCI = 8	4999.99	0.02	1.22	99.79	49.68	15
	SCI = 10	6799.99	0.02	1.20	99.84	57.37	15
	SCI = 12	8999.98	0.02	1.17	99.89	61.34	15
	SCI = 14	10,599.97	0.02	1.20	99.83	94.82	15
	SCI = 16	12,999.98	0.02	1.18	99.87	98.02	15
	SCI = 18	15,099.99	0.02	1.20	99.86	119.74	15
	SCI = 20	17,299.98	0.02	1.20	99.83	155.35	15
Populus nigra	SCI = 6	2300.00	0.05	1.34	99.94	11.00	13
	SCI = 8	2700.00	0.07	1.55	99.86	24.85	13
	SCI = 10	3500.00	0.06	1.47	99.89	27.75	13
	SCI = 12	3900.00	0.08	1.56	99.86	37.17	13
	SCI = 14	4700.00	0.07	1.41	99.96	22.53	13
	SCI = 16	5300.00	0.07	1.36	99.98	19.64	13
Quercus mongolica	SCI = 8	9300.00	0.02	1.18	99.52	147.15	19
	SCI = 10	11,300.00	0.02	1.16	99.61	157.85	19
	SCI = 12	12,999.99	0.02	1.13	99.71	158.48	19
	SCI = 14	15,699.99	0.01	1.08	99.78	155.56	19
	SCI = 16	16,999.93	0.01	1.09	99.76	181.02	19
	SCI = 18	18,599.98	0.01	1.09	99.77	197.61	19
Pinus sylvestris	SCI = 6	3143.28	0.06	3.16	99.98	11.90	21
	SCI = 8	4861.32	0.06	3.32	99.97	21.60	21
	SCI = 10	6792.03	0.06	3.42	99.97	29.06	21
	SCI = 12	8998.50	0.06	3.50	99.98	38.02	21
	SCI = 14	11,622.27	0.06	3.41	99.97	52.19	21
	SCI = 16	13,999.99	0.06	3.60	99.98	59.10	21
	SCI = 18	16,500.00	0.06	3.81	99.96	93.98	21

and

Table A2. The Regression Coefficients for Equations (6) and (7) in three main forest type groups.

Forest Type Group	SCI	A	b	c	R2	RMSE	NO.
Coniferous forest	SCI = 12	8499.97	0.01	1.19	99.98	22.00	29
	SCI = 14	9399.97	0.01	1.22	99.96	38.23	29
	SCI = 16	10,999.99	0.01	1.20	99.97	40.89	29
	SCI = 18	12,999.99	0.01	1.16	99.98	31.24	29
	SCI = 20	13,999.99	0.01	1.18	99.97	48.28	29
	SCI = 22	15,999.88	0.01	1.14	99.99	36.86	29
	SCI = 24	16,999.96	0.01	1.15	99.98	54.04	29
Broadleaved forest	SCI = 8	9800.00	7.90	0.48	99.78	60.84	17
	SCI = 10	12,000.00	8.08	0.49	99.84	63.76	17
	SCI = 12	15,500.00	7.77	0.47	99.92	52.61	17
	SCI = 14	16,500.00	8.10	0.50	99.92	64.67	17
	SCI = 16	19,000.00	8.00	0.49	99.94	61.86	17
	SCI = 18	21,998.88	7.74	0.48	99.97	53.23	17
	SCI = 20	23,999.99	7.74	0.48	99.97	55.17	17
Mixed forest	SCI = 8	8000.00	12.16	0.63	99.69	70.56	23
	SCI = 10	10,000.00	12.12	0.63	99.79	73.14	23
	SCI = 12	11,500.00	12.58	0.65	99.82	81.08	23
	SCI = 14	13,900.00	11.70	0.63	99.90	70.77	23
	SCI = 16	15,000.00	12.37	0.66	99.89	85.60	23
	SCI = 18	16,000.00	13.09	0.70	99.87	103.89	23
	SCI = 20	17,000.00	13.78	0.73	99.85	122.83	23

.

3.3. The Second Component of NPP—Mortality

Mortality represents a substantial component of C cycling. The yield tables do not include information about standing dead trees and the downed dead wood. In this study we used the average mortality rate of biomass ($\overline{MR}\%$) for different age groups for nine forest types and three forest type groups in Heilongjiang Province that were published by Jia [45]. We used $\overline{MR}\%$ instead of a mortality equation due to insufficient mortality information in the yield tables. Average mortality rates for the main forest types and groups of Heilongjiang Province are only a crude approximation of the actual mortality. We limited mortality to a constant value after a forest was considered mature. Past research had indicated that there was no relationship between mortality and SCI [46]. The average mortality rates for the 12 forest types and groups are listed in

Table 3. The average mortality rate in percent volume change for nine forest types and three forest type groups.

Forest Type and Group	Age Group	$\overline{\mathit{M}\mathit{R}}\mathbf{\%}$	Forest Type and Group	Age Group	$\overline{\mathit{M}\mathit{R}}\mathbf{\%}$
Pinus koraiensis	Young	0.51	Betula platyphylla	Young	0.75
	Middle aged	0.25		Middle aged	1.19
	Near mature	0.50		Near mature	0.99
	Mature	0.30		Mature	1.29
Larix gmelinii	Young	0.71	Betula davurica	Young	0.33
	Middle aged	0.70		Middle aged	0.64
	Near mature	0.83		Near mature	0.40
	Mature	0.57		Mature	0.35
Pinus sylvestris	Young	0.11	Populus nigra	Young	1.08
	Middle aged	0.33		Middle aged	1.06
	Near mature	0.24		Near mature	0.66
	Mature	0.64		Mature	0.62
Quercus mongolica	Young	1.45	Coniferous forest	Young	1.12
	Middle aged	0.88		Middle aged	0.84
	Near mature	0.52		Near mature	0.56
	Mature	0.98		Mature	0.78
Populus davidiana	Young	1.08	Broadleaved forest	Young	0.54
	Middle aged	1.06		Middle aged	1.09
	Near mature	0.66		Near mature	0.95
	Mature	0.62		Mature	1.06
Tilia amurensis	Young	0.35	Mixed forest	Young	0.72
	Middle aged	0.39		Middle aged	1.08
	Near mature	0.46		Near mature	0.86
	Mature	0.37		Mature	1.93

.

3.4. The Third Component of NPP—Tissue Turnover

We calculated the foliage turnover as the product of foliage biomass and a species-dependent turnover rate as [26]:

$$L_f=B_f\times T_f\times C_f$$

(6)

where B_f is the biomass of foliage and T_f is the foliage turnover rate that differs by forest types. For broadleaved forest (DBF), T_f was 1. C_f is the C content in foliage, (C_f = 0.5) [47].

Thornton et al. [48] found a relationship between the turnover of foliage and the turnover of fine roots, so L_fr could be calculated using the following equation:

$$L_{fr}=L_f\times e$$

(7)

where e is the new fine root carbon to new leaf carbon allocation. T_f and e are obtained from White [49] and are listed in

Table 4. New fine root carbon to new leaf carbon allocation values and foliage turnover rates for 12 forest types and groups.

Forest Type and Group	Tf	e
Pinus koraiensis	0.26 ± 0.15	1.4 ± 1.5
Larix gmelinii	1	1.4 ± 1.5
Pinus sylvestris	0.383 ± 0.02	1.33 ± 0.79
Quercus mongolica	1	1.2 ± 1.2
Populus davidiana	1	1.2 ± 1.2
Tilia amurensis	1	1.2 ± 1.2
Betula platyphylla	1	1.26
Betula davurica	1	1.26
Populus nigra	1	1.2 ± 1.2
Coniferous forest	0.26 ± 0.15	1.4 ± 1.5
Broadleaved forest	1	1.2 ± 1.2
Mixed forest	1	1.3 ± 0.1

.

3.5. The Fourth Component of NPP—The NPP of the Understory

The NPP of understory vegetation could not be calculated from yield tables and biomass equations, so we selected some sample plots on Maoer mountain in Heilongjiang Province to fit the change in NPP_und with age. Fifty plots for three forest type groups were selected based on stand age, from young to mature. The plot size was 20 m × 30 m, further divided into four smaller plots that were located in each of the four quadrants of the plot. The area of these small plots was 1 m × 1 m. The tissues of annual herbs, new twigs, new foliage, and all fine roots with a diameter smaller than 2 mm were collected from the small plots. We first measured the fresh weight of these tissues then samples were dried at 70 °C and weighed to the nearest 0.1 g. NPP_und was calculated as the sum of the dry biomass of herbs, new twigs, new foliage components, and new fine roots. We used nonlinear regression (SPSS 2017) to estimate NPP_und as a function of stand age as [17]:

$$NP{P}_{und}=A\times e^{-k\times age}$$

(8)

where A and k are empirical coefficients. This equation was applied for all ages in the range from 10 to 60. For all other ages, the maximum value of the equation for the age <10 and the minimum value for the age >60 were assigned due to the lack of data. The estimates of coefficients are listed in

Table 5. Coefficient estimates and goodness-of-fit statistics of the NPP_und equations.

Forest Type and Group	A	k	R2	RMSE	NO.
Coniferous forest	82.372	0.033	0.61	16.1	15
Broadleaved forest	92.729	0.034	0.65	12.3	25
Mixed forest	65.823	0.017	0.42	13.4	10

.

3.6. NPP-Age Relationships

We used nonlinear regression (SPSS 2017) to estimate NPP as a function of stand age and SCI as Equation (9). The equation was obtained by [26] that could fit the change in NPP at different ages and SCI, for each forest type and group:

$$NPP(age,SCI)=a\left[1+\frac{bSCI{(\raisebox{1ex}{$age$}\!\left/ \!\raisebox{-1ex}{$c$}\right.)}^d-1}{e^{(\raisebox{1ex}{$age$}\!\left/ \!\raisebox{-1ex}{$c$}\right.)}}\right]$$

(9)

where a, b, c, and d are coefficients depending on site conditions.

The R² and RMSE of the relationship between NPP and age and SCI were listed in

Table 6. The four regression coefficients for Equation (9) for the nine forest types and three forest type groups.

Forest Type and Group	a	b	c	d	R2	RMSE	NO.
Pinus koraiensis	100	0.7021	8.0106	2.0028	88.82	39.58	88
Larix gmelinii	100	1.1143	15.3176	1.2058	97.91	28.06	126
Pinus sylvestris	34.8865	1.6934	11.9259	1.8703	98.62	15.9	154
Quercus mongolica	41.038	0.928	21.2796	0.4599	97.22	11.93	120
Populus davidiana	35.3185	0.7109	46.1129	0.0698	98.15	11.19	112
Tilia amurensis	120.7272	0.2344	47.6497	0.0283	98.19	11.61	80
Betula platyphylla	54.6661	0.6702	27.1883	0.4812	99.12	6.76	70
Betula davurica	23.1453	0.7132	34.692	0.0621	98.11	5.9	70
Populus nigra	132.9282	0.3969	9.6522	0.4974	99.45	6.58	84
Coniferous forest	91.8927	0.2705	30.1884	0.5506	98.36	7.26	210
Broadleaved forest	57.6312	0.6305	15.9528	0.4594	98.46	9.01	147
Mixed forest	74.5847	0.502	28.5017	0.3161	99.18	7.35	168

, where the R² and the RMSE illustrated that the developed method could fit NPP-relationships at various SCI values well.

4. Results

4.1. Results for NPP-Age Relationships

The increase in biomass reached a peak in 10 to 20 years for broadleaved forests and in 20 to 30 years for coniferous forests (Figure 5 and Appendix A Figure A1). Afterward, the biomass increase began to decrease and was almost unchanged after the forest reached maturity; meanwhile, the foliage and fine root turnover became the main components of NPP at older ages. The mortality varied greatly with stand age; it generally reached a maximum after maturity due to competition from neighboring trees. The mortality of broadleaved forests was higher than that of coniferous forests. For the three mixed forest type groups, the maximum biomass increase in coniferous mixed forests was the least, and the age at which biomass growth reached a peak was the oldest. Broadleaved mixed forests had the largest biomass increase, with the peak growth age occurring at approximately 17 years, while in coniferous mixed forests and coniferous and broadleaved mixed forests it was approximately 20 to 30 years.

Different SCI values had different biomass increase curves, and from Figure 5 we found the sites that were the most productive due to having the better site conditions, shown by having higher maximum biomass increases and faster rates of biomass growth. The maximum biomass increase at maximum-SCI was up to two to three times greater than the minimum-SCI observed. The other components of NPP at the better sites were also higher than in the lower site conditions.

We divided the 12 forest types and groups into three main forest type groups: broadleaved forest (BF) (including BP, BD, TA, QM, PD, PN, and BMF), coniferous forest (CF) (including PK, LG, PS, and CMF) and coniferous and broadleaved mixed forest (MF). We compared all of the NPP for different age classes (young, mature, and old stands) in two categories of SCI values for the 12 forest types and groups, as shown in Figure 6 and Appendix A Figure A2. The NPP-age relationships at various SCI values for different forest types and groups presented similar trends. NPP increased with stand age initially, reached a peak, then declined due to slower biomass increases, and finally stayed relatively stable as determined by foliage and fine-root turnover in short time periods. However, for different forest types and groups, the age classes at which NPP reached a maximum were different. For BF and MF, the NPP increased rapidly and reached a maximum when young; thereafter, NPP began to decline gradually when mature and old. For CF, the maximum NPP occurred when a mature stand generally, and the range of NPP decline was smaller than BF and remained at a value of approximately 70% of the maximum in old age. Within the same forest type, higher SCI values led to a higher NPP and the extent of variation was greater than in the less-productive sites, where the trend was similar to the biomass increase. In reality, the effects of stand break-up or species succession on NPP could be significant, but they were not included in our analysis due to lack of data for the seed regeneration of trees.

4.2. Results for Combined NPP-Age Relationships into the InTEC Model

The average NPP of CF, BF, and MF were 330 g C·m⁻²·y⁻¹, 296 g C·m⁻²·y⁻¹, and 305 g C·m⁻²·y⁻¹, respectively (Figure 7). From Figure 8, we found that the age of the majority of CF were between 40 and 60 years old, which was approaching the mature stand stage of CF. At this age, CF had greater carbon fixation capability than the other two forest type groups. However, 60% BF were between 20 and 60 years old, and at this time the NPP of BF began to decline. The NPP of MF was almost the same as CF, due to MF was younger than the other two categories, and this was during the period of the highest carbon fixation.

For the purpose of analyzing the influence of different SCI on NPP, we calculated the average NPP from 2000 to 2010 using SCI values of 10–12, 12–14, and 14–16 for the three main forest type groups (Figure 9), because we found that the distribution of SCI were concentrated mainly between 10 and 16. The average NPP of all the three forest type groups presented a similar fluctuating trend at various values of SCI. The larger values of SCI led to a higher NPP. For CF, the variation of NPP along SCI was not significant. For BF, the increase in NPP along SCI was more obvious than CF, where there was a 24% increase from minimum to maximum SCI. MF had the largest NPP increase, with an average NPP increase of 46% from minimum to maximum SCI.

We first compared the InTEC results of NPP in 2010 against the NPP before the NPP-age relationships were modified by the new NPP-age relationships at various SCI values in InTEC (Figure 10). After this modification with the new NPP-age curves, the average value of NPP for Heilongjiang Province was 389 g C·m⁻²·y⁻¹. Nearly the entire north region of Heilongjiang had an NPP value above 400 g C·m⁻²·y⁻¹. From Figure 3, we could find the main forest type was Larix gmelinii forest (in green), covering 78% in the north region of Heilongjiang, which had a stronger carbon sequestration capability than other forest types. The average NPP of the middle and south regions of Heilongjiang was divided into two levels and showed an increasing trend from north to south, which agreed with the trend of SCI (Figure 4). Before modification, the average value of NPP for Heilongjiang Province was 389 g C·m⁻²·y⁻¹, showing only a small difference in the provincial average. However, there was a significant difference in the spatial distribution of NPP. Before this modification with the new NPP-age curves, we observed an obvious overestimate of NPP in the north and south region of Heilongjiang compared with that after modification. This is because the original InTEC model included only one curve (SCI = 14), however, we could find that (Figure 4) the real SCI in the north and south region of Heilongjiang should be lower than 14. This showed that the NPP-age relationships with different site conditions had significant differences and the NPP-age relationships which developed in this study make the spatial distribution of NPP more reliable for carbon cycle modeling.

Using the national forest inventory data, the NPP of every plot could be estimated using the same method as with yield tables. We compared the NPP results to measurements derived by [50] using InTEC, before and after the SCI modification of the NPP-age relationship (Figure 11). The R² value of the linear regression increased from 40% to 53%, and RMSE reduced from 77 to 69 g C·m⁻²·y⁻¹ after modification of the NPP-age relationships in InTEC with those derived from this study. These results showed the suitability of the method and that the new NPP-age relationships could improve the accuracy of NPP estimation using InTEC.

4.3. Uncertainties

The model for total living biomass had the highest accuracy (

Table 7. The Mean absolute percent errors (MAE%) of NPP and several components for the three main forest type groups.

Forest Type Group	NPP and Components	Mean Absolute Percent Errors (MAE%)	No.
Coniferous forest	Total biomass	14	173
	Foliage biomass	14	173
	Mortality	43	345
	NPPund	49	15
	NPP	73	31
Broadleaved forest	Total biomass	15	2522
	Foliage biomass	43	2522
	Mortality	67	497
	NPPund	41	25
	NPP	46	93
Mixed forest	Total biomass	12	519
	Foliage biomass	10	519
	Mortality	61	276
	NPPund	66	10
	NPP	58	54

), and the mean absolute percent errors (MAE%) of the total estimated total biomass ranged from 12% to 15%, with a mean of about 14%. Broadleaved forest had the lowest precision among the three type groups. Estimated foliage biomass precision was obviously lower than total biomass; the MAE% ranged from 10% to 43%, with a mean of about 22%. Due to the yield tables not including information about mortality, hence we estimated mortality with the average mortality rate of biomass. We used the mortality data which was reported by Jia to estimate the error of mortality in our study. The MAE% for mortality was nearly 57% (ranging from 43% to 67%). We could not collect the measured data about tissue turnover (foliage and fine roots), however, based on a similar study published by Chen [17], the error for foliage turnover could be as high as 70% and 33% for the fine-root turnover. Due to the varying complexity of understory for different forest type groups, the estimate for NPP_und also contained a high level of error, with a mean of about 52% (ranging from 41% to 66%).

The calculated NPP in this study was validated using NPP database created by Yu based on field measurements, which has been widely used to validate forest NPP [40,50]. The MAE% of our estimation of stand NPP ranged from 46% to 73%, with a mean of about 59%. The error was lower in the early stages of development because total NPP was dominated by the growth of the biomass components that could be estimated with comparatively low error. As stands developed, however, NPP became increasingly dominated by the turnover of roots and litterfall, both of which had a high level of associated error.

In NPP-age relationships at various SCI values, other errors may arise from errors inherent in the yield tables, as the age and SCI in the yield tables were also affected by sample sizes, site locations, environmental variability, and stand density. The parameters of T_f and e could not be obtained for the study region of China, so we used values from White et al. [49], which also brought associated error that could not be ignored.

5. Discussion

Our results showed that NPP varied with stand age and SCI, presenting a generally consistent trend of NPP-age relationships at various SCI values for different forest types and groups across the Heilongjiang Province, China. The NPP-age relationship curves at various SCI values were determined by components including the increase in total biomass, mortality, the fine root and foliage turnover rates, and the NPP of understory vegetation and moss. This method was similar in principle to that of Gower et al. [51] and Bernier et al. [52], who estimated stand-level NPP in a specific year using measured annual radial increased, biomass allometric equations, and measured foliage and root turnover rates. In the growth of young forests, NPP increasing rapidly with age was caused mainly by increases in biomass growth. Zha et al. [14] reported that aboveground NPP (ANPP) increased with age, stabilizing after 25 years. These results aligned well with those of our study. For old stand forests, where the products of photosynthesis were mainly consumed by foliage and fine root turnover, NPP began to decrease. Ryan et al. [13] reported that the range of NPP decreased by one half to one third from the peak growth period. Our results showed that it remained steady at 30% to 60% of the maximum NPP.

NPP-age relationships varied with site quality; higher productivity stands reached maximum NPP earlier than less productive stands. In addition, the decrease in NPP with age was also more dramatic for more productive sites. Our results were similar to those reported by Chen et al. [17].

We compared our NPP-age relationships with similar studies for China and the U.S. [24,26]. Our NPP-age curves were similar to the curves derived for China and U.S. forests except for DBF. Wang et al. [24] reported that the age at which the NPP reached peak value was always older, especially for DBF, which did not show a decline until the age of 120. However, our NPP-age relationships for DBF showed the decline of NPP at approximately 20 years of stand age. Three possible reasons included: first, there were only 37 samples from Heilongjiang Province in Wang’s study; second, Heilongjiang Province has fertile soil and earlier forest maturity than other regions; and lastly, the NPP-age relationships derived from BEPS modeling used in Wang et al. [24] may cause more error than that from yield tables.

Two studies had found old forests to be as productive as young forest stands [15,16]. In our study, the predicted NPP of old forests declined by an average value of 33% from its peak for all of forest types and groups studied. This may be caused by many other site factors, including climate, soil, and drainage. These factors caused the decline of biomass increase in old stands.

We implemented these NPP-age relationships into the InTEC model to estimate the NPP in Heilongjiang from 2000 to 2010. Our results showed that the average NPP of coniferous forest type groups was higher than that for the other two forest type groups. This aligned with the research of Mao et al. [53], which simulated NPP based on MODIS and AVHRR NDVI data sets. Our results also indicated the younger stands and higher SCI values may lead to a higher NPP, which agrees with the findings of Chen et al. [17]. However, the average NPP of the three main forest type groups was slightly lower than found in the research of Mao et al. [53]. MODIS NPP had been proven to oversimulate because the product was strongly driven by climate data [54].

6. Conclusions

In this study, we developed a specific method for estimating the relationships between NPP and age, at different SCI values for 12 forest types and groups, based on data from yield tables, biomass equations, and forest inventories. We derived the stand-level NPP through quantifying the four components: total biomass annual increment, tree mortality, turnover of foliage and fine roots, and the NPP of understory vegetation and moss.

Two similarities in the temporal patterns of the NPP-age relationships at various SCI values occurred between various forest types and groups: first, NPP increased rapidly during early stages of stand growth, reached a peak when young or mature, and then declined after maturity; secondly, high productivity sites had higher maximum NPP than poor productivity sites.

We implemented these NPP-age relationships at various SCI values into the InTEC model to estimate the NPP for three major forest biomes: coniferous forest, broadleaved forest, and mixed forest in Heilongjiang Province, China, from 2001 to 2010. During this period, the average NPP of coniferous forest was larger than broadleaved forest and mixed forest due to a greater carbon fixation capability and being productive at a mature age.

In this study, we also compared the results from InTEC modeling after modification with the new NPP-age relationships at various SCI values derived from this study. The model results showed that the new NPP-age relationships could better fit the NPP changes with age and SCI variation than the old ones used in the model.

However, there were still some weaknesses in our study: First, because mortality information was not included in the yield tables, we used the average mortality rate of biomass as the first order approximation; and second, the parameters of foliage and fine root turnover rates were obtained from White et al. [49] for North American forests, because of the lack of local measurements for these parameters. Our future studies could be directed toward acquiring this information in our region.

Although large uncertainties still exist in the development of these NPP-age relationships, this was nevertheless the first study in China to develop localized NPP-age relationships at various SCI values. These relationships have many potential uses for the analysis of forest management because they provide a new, independent, and comprehensive source of information on forest growth. This study should have implications for other forest regions of the world where SCI varies considerably.