**Carbon Mass Change and Its Drivers in a Boreal Coniferous Forest in the Qilian Mountains, China from 1964 to 2013**

#### **Shu Fang 1,2, Zhibin He 1,\*, Jun Du 1, Longfei Chen 1, Pengfei Lin 1,2 and Minmin Zhao 1,2**


Received: 16 November 2017; Accepted: 22 January 2018; Published: 25 January 2018

**Abstract:** Carbon storage of mountain forests is vulnerable to climate change but the changes in carbon flux through time are poorly understood. Moreover, the relative contributions to carbon flux of drivers such as climate and atmospheric CO2 still have significant uncertainties. We used the dynamic model LPJ-GUESS with climate data from twelve meteorological stations in the Qilian Mountains, China to simulate changes in carbon mass of a montane boreal forest, and the influence of temperature, precipitation, and CO2 concentration from 1964 to 2013 on carbon flux. The results showed that the carbon mass has increased 1.202 kg/m<sup>2</sup> from 1964 to 2013, and net primary productivity (NPP) ranged from 0.997 to 1.122 kg/m2/year. We concluded that the highest carbon mass proportion for this montane boreal forest was at altitudes 2700–3100 m (proportion of ecosystem carbon was between 93–97%), with maximum carbon density observed at 2700–2900 m. In the last 50 years, the increase in precipitation and in CO2 concentration is expected to increase carbon mass and NPP of *Picea crassifolia* Kom. (Pinaceae) (Qinghai spruce). The effect of temperature on NPP was positive but that on carbon mass was not clear. The increase in CO2 concentration over the past 50 years was a major contributor to the increase in carbon storage, and drought was the foremost limiting factor in carbon storage capacity of this montane boreal forest. *Picea crassifolia* forest was vulnerable to climate change. Further studies need to focus on the impact of extreme weather, especially drought, on carbon storage in *Picea crassifolia* forests.

**Keywords:** carbon mass; NPP; *Picea crassifolia*; climate change

#### **1. Introduction**

Climate change influences the carbon and water cycles in mountain areas [1]. Forests in mountain areas provide important ecological and socio-economic services; however, carbon storage is vulnerable to the effects of climate change and may decrease its ecological service function over time [2,3]. Research on the response of mountain ecosystems, especially of alpine forests and tree species, to climate change, is lacking [1,2]. Additionally, the impacts of climate change on montane forests are examined distinctively with dendroclimatic techniques [2–4]. High mountain boreal forests, which have numerous organic pools stored aboveground and in the permafrost, play an important role in regional carbon budgets and are exposed to rapid climate change; however, research on carbon pools at high elevations and controls on carbon flux through time lag behind [5,6] due primarily to insufficient data [7]. Ecosystem models are useful tools for describing and quantifying water, matter, and momentum fluxes between the biosphere and the atmosphere [8]. Using biomass models to

obtain mass data can facilitate simulation of forest carbon dynamics in areas where such data are lacking, while analyses of the montane boreal forest biome will improve the accuracy of global carbon models [3].

Climate change is primarily related to changes since the 1950s in solar radiation, temperature, and precipitation that have not been observed before [9]. Climate change and related disturbances may substantially impact the species ranges, population sizes, and extinction risks in mountain forests, and even shift the carbon sink activity of a forest to a net carbon source [10–12]. In a simulation of a species distribution in Mediterranean mountains, the cold-adapted species decreased in significance, and the species range and communities changed under global warming [13]. Extreme climates may change plant intrinsic characteristics and increase the possibility of mountain ecosystem shifts from forest to shrubland or grassland [14]. These climate change-induced species conversions will result in a long-term loss of carbon stock in mountain forest ecosystems [15]. One stark example of this was observed with mountain pine beetle in western North America; this insect outbreak, caused by climate change, compromised the ability of forests to take up and store carbon [12].

The response of mountain forests to climate change depends heavily on relative changes in temperature, precipitation, CO2, and their interaction effect [10]; the relative contributions of individual historical drivers can be assessed by high-resolution climate data and modeling [16,17]. Sensitivity analyses of climate factors focused on tropical and temperate areas. However, analyses of the driving factors of carbon storage changes in recent years at high latitudes and in high elevation areas are still lacking [18]. The sensitivity to climate change of carbon storage in mountain forests is based on initial climatic variability; thus, increased mean annual temperature at high elevations results in sensitivity of tree growth and changes in forest productivity [3]. Examined by three ecological models (LPJ-GUESS, ForClim, LandClim), carbon storage across all elevations revealed a lower sensitivity than other ecosystem services to a 2 ◦C warming [19]. However, under the influence of global warming, mountain forests at high altitudes exhibited positive growth, whereas drought stress led to a negative effect on carbon mass at lower altitudes [4]. The impacts of climate change on mountain forests were also closely related to the frequency and intensity of extreme weather, especially precipitation-related events [5]. Extreme warming events may influence mountain boreal plant activity, plant litter production, soil moisture, and insect life cycles, while extreme weather events such as precipitation extremes and severe storms may lead to damage and may influence the recovery capacity of mountain boreal forest ecosystems [6]. Uncertainties exist as to the main climatic controls on carbon changes in long-term simulations of carbon mass in mountain forests [9,10]. Investigating carbon flux sensitivity to climate change will reduce the uncertainty of the past factors of climate change and help guide forest management. For a better understanding of the effects of climate change on boreal forests, it is necessary to quantify the driving parameters of carbon changes in boreal forests during the last 50 years.

The Qilian Mountains, located in the northwestern part of China at the north-eastern edge of the Tibetan Plateau, are typical arid mountains, and the source of key inland rivers; these mountains are sensitive to climate change [20]. The average temperature rise is 80% higher than the national rate (0.14 ◦C/10a), and extreme climate and precipitation events occurred more frequently here in the past 50 years than in other parts of China [21]. Under climate change, environmental degradation, soil erosion, and loss of biodiversity have increased, and vegetation cover and growth in the Qilian Mountains have changed significantly in recent years [22]. Montane Boreal forest is the dominant forest type, and it plays important roles in the hydrology and biogeochemistry of the area due to high altitudes (2000–5500 m) and minimal human-induced ecological deterioration [23]. Between 1957 and 2007, temperature increased by 0.29 ◦C/10a and precipitation by 5.5 mm/10a, resulting in a shift of the tree line of the montane boreal forest in this area to higher elevations [24]. Simple regression models of phenological dynamics in this montane boreal forest revealed that the length of the growing season can be expected to increase during the next two decades [25]. This prompted new research efforts to determine the relationships between montane boreal forest carbon dynamics and climate

change. Ecosystem stability of the montane boreal forest in the Qilian Mountains supports downstream economy and production because the forest serves as the most important water conservation entity in the area. Better understanding of the role of climate and atmospheric CO2 in determining carbon flux in this montane boreal forest, and the relationship between the NPP and biomass is critical for guiding ecology policy to optimize carbon storage.

To address this knowledge need, we used climate and CO2 concentration data to reconstruct carbon storage in living biomass for the montane boreal forest in the Qilian Mountains, China from 1964–2013. With this, we examined the carbon stock in this forest type at different elevations and its carbon flows during the last 50 years. Subsequently, we used the LPJ-GUESS model to explore the sensitivity of forest carbon stocks in this region to temperature, precipitation, CO2, and extreme conditions to determine the main past factor controlling the change in carbon flux in the last 50 years.

#### **2. Materials and Models**

#### *2.1. Area Description and Data Collection*

The Qilian Mountains are one of the major arid mountain ranges in northwestern China, located in the northeastern Tibetan Plateau, between 93◦33 36–103◦54 00 E and 35◦50 24–40◦01 12 N. The distribution of vegetation is strongly controlled by altitude and can be separated into five types (from low to high altitude): desert steppe, forest steppe, sub-alpine shrubby meadow, alpine cold desert, and ice/snow zone. In the forest steppe, the montane boreal forest is located at high altitudes where the level of human impact is low [26]. The dominant boreal species, *Picea crassifolia*, grows on shady and partly-shady north slopes at altitudes ranging from 2600 to 3400 m, and accounts for as much as 76% of the total forested area in the forest steppe zone [27]. In total, twelve sampling sites were selected across an elevation gradient from 2300–3500 m at intervals of 200 m to analyze the carbon change of montane boreal forest by altitude.

We selected twelve meteorological stations located within the tree-growing areas to extract the climate parameters (Figure 1). In an earlier study of the spatial distribution of *Picea crassifolia* biomass and carbon storage in the Qilian Mountains, the average carbon content was calculated as 0.52 [28]; this was similar to the results based on a synthesis of carbon contents for temperate and boreal conifer wood (*n* = 36) of 50.8 ± 0.6% [29]. In this study, we used a carbon content of 0.52 for *Picea crassifolia*.

**Figure 1.** Study area and locations of meteorological stations used in the study in the Qilian Mountains.

The climate data, including mean monthly temperature, monthly precipitation, mean monthly daily proportion of sunshine hours, and monthly rainy days, were obtained from the climate data repository of Chinese meteorological stations (http://data.cma.cn/site/index.html), and of Qinghai Province (Table 1). Atmospheric CO2 concentrations for 1964–2013 were derived from the atmospheric carbon dioxide mixing ratios from NOAA (National Oceanic and Atmospheric Administration) ESRL (Earth System Research Laboratory) Carbon Cycle Cooperative Global Air Sampling Network (http://www.esrl.noaa.gov/gmd/ccgg/mbl/data.php).


**Table 1.** Site information.

#### *2.2. The LPJ-GUESS Model*

LPJ-GUESS is a validated vegetation carbon dynamics model [30]. The basic use of the model is to simulate carbon balance of different tree species and to predict changes in carbon pools and fluxes under global climate change [31,32]. A LPJ-GUESS simulation of the spatial and temporal patterns of carbon fluxes associated with regrowth after agricultural abandonment indicated that semi-arid regions, where carbon balance was strongly associated with both precipitation and temperature, were key to the understanding and predicting of the global carbon cycle [33]. Further, changes in the spatial patterns of wildfires, estimated by LPJ-GUESS, were critical to the proposal for a more reasonable climate policy [34]. The model can simulate and predict the responses of plant water fluxes to elevated CO2 at leaf and stand scales [35]. In addition to natural vegetation, agricultural crop production and crop responses to climate change were analyzed in an effort to increase food security [36]. In China, the model has been applied in subtropical, temperate, and mountain zones [37–39], and to the entire country [40]. Carbon simulation by LPJ-GUESS using an individual climate factor demonstrated carbon release by management, and CO2 as the most important driver for carbon change in Lady Park Wood in the last 20 years [16].

Plant functional types (PFTs), which can simplify species diversity of vegetation, represent groups of species with similar functional traits. Among forest carbon dynamics simulation models, both, the individual-based model GUESS (General Ecosystem Simulator) [30], and the area-based dynamic global vegetation model LPJ (Lund-Potsdam-Jena) [41], can be successful in predicting the carbon flux of PFTs. Subsequently, a new vegetation dynamics model, LPJ-GUESS, was developed by combining the two other models [30].

Two model frameworks exist, the cohort and population. We employed the "cohort mode" in this study. In the "cohort mode", individual trees are distinguished, but are identical within each cohort (age class) [30]. Population processes and disturbances are modeled stochastically, and stand characteristics are averaged over 100 patches of 0.1 ha, representing "random samples" of the simulated stand. The model is driven by short-wave radiation (photosynthetically active light), temperature, precipitation, CO2 concentration of the air, and soil conditions. The CO2 influences assimilation rate, and soil conditions modify plant water uptake; soil data are already provided in the model. Species information includes the physiological characteristics of PFTs, such as prescribed allometric relationships. Based on physiology, morphology, phenology, and on the response to disturbance and to bio-climatic limiting factors, the model defines ten PFTs, of which 8 are woody, and two are herbaceous [42]. We used the full set of PFTs and added a new PFT for *Picea crassifolia* to reproduce the current stage of the vegetation.

#### *2.3. Model Forcing and Simulation Protocol*

The simulation normally follows two or three phases and begins with "bare ground", which means that the modelled area is bare, with no vegetation present. The first phase of the simulation is known as the "spin-up"; in this phase, input data are normally based on the first few years of available (historical) data. During the subsequent, "historical" phase, the model uses "observed" climate and CO2 data as input. We chose two phases, "spin-up" and "historical" to run the model for the last fifty years. Climate data for an initial 300-year "spin-up" phase were not available; such "spin-up" equilibrates the initial vegetation and carbon pools with climate at the beginning of the study period. The model was first "spun-up" for 300 years, recycling the observed time series from 1964–1993 [30,43]. The "historical" period then ran from 1964 to 2013; the simulation time interval for the LPJ-GUESS model is five years.

In a study to investigate the ability of LPJ-GUESS to reproduce features of real vegetation, the model was demonstrated to not require site-specific calibration and could be used to simulate vegetation dynamics on a regional basis or under past or future climates and atmospheric CO2 levels for reparameterization becauase plant growth is modeled mechanistically [31]. Many of the parameters of species are decided by its PFTs, such as whether it is needle-leaved or broad-leaved; or whether it is boreal or temperate species. The species-specific parameters such as tree longevity are defined in the article. We referred to the relevant parameters for the boreal forest to define the parameters of *Picea crassifolia* [31,44]. We specified the mean length of the life of foliage at 11.8 years, as determined in a field investigation [45], and tree longevity at 250 years, as described in the literature [46]. Data used for biomass calculations for the Qilian Mountains were obtained from the literature [28,47–49] and were used to verify the accuracy of those that were calculated.

#### *2.4. Sensitivity Analysis of Temperature, CO2, and Extreme Whether*

Sensitivity analysis is an important aspect of evaluating the resilience of ecosystems to climate change and can be combined with models to measure the effect of changing input parameters under some climate scenarios and extreme weather conditions [50,51]. Based on the simulation of LPJ-GUESS, further research was conducted to determine the relative effect of climate variables and CO2 concentration on the carbon flux of montane boreal forests in Qilian Mountains for the period 1964–2013.

The temperature increased during the simulation period, but this trend could not be fit by a linear model. Further, the trend in precipitation was not clear. Thus, we first examined the climate effect by simulating the temperature, precipitation, and CO2 independently of each other. Using these climate simulations as a baseline, the uncertainties of these parameters were then compared. The influence of the temperature was observed for an increase or decrease of 1 and 2 ◦C; precipitation was changed multiplicatively, because it is a zero-based variable, by increasing or decreasing by 10% and 20%. The influence of CO2 was examined by removing its trend and increasing the value by 50 ppm and 100 ppm.

Because the climate variables are strongly related to each other, a comparison of the effects of climate change on carbon flux based on changing a single parameter in the data was incomplete [6,16]. Hence, to preserve the relationship of temperature, precipitation, and radiation, we represented extreme weather as follows to determine the impact of climate factors on carbon mass. Local weather

in the study area was ranked in terms of temperature or precipitation levels (annual mean temperature from April to October, and annual mean precipitation from May to September), and the top five warmest, coldest, wettest, and driest years; these were then cycled through the model repeatedly to simulate an extreme climate.

#### **3. Results**

#### *3.1. Changes in Carbon Mass with Altitude*

The monthly mean temperature and precipitation in the simulated period followed normal distribution. Warm periods were concentrated in April to October, and high precipitation was concentrated in May to September; almost all of the stations exhibited the same trends. Mean annual temperature at the twelve stations increased by 0.29–3.69 ◦C between 1964 and 2013, with an average increase of 0.73 ◦C/10a. Differences in annual precipitation among altitudes were not apparent.

Carbon stocks of the montane boreal forest were concentrated at 2700–3100 m and exhibited a single peak of 11.787 kg/m<sup>2</sup> at 2700–2900 m in 2013, with the average carbon mass at 2700–3100 m of 10.503 kg/m2. The calculation of the biomass of *Picea crassifolia* was different for the different measurement methods and research areas, but the results of our simulation were within the range of other studies (Table 2). The proportion of montane boreal forest was >28% at all altitudes. Specifically, montane boreal forest accounted for 30% at 2300 to 2700 m, increasing to 93% at 2700–2900 m and 97% at 2900–3100 m, declining to 37%, and again reaching 49% at 3300–3500 m (Figure 2). The carbon stored in the ecosystem at 2300–2700 m was 7.068–7.591 kg·C/m2, declining to about 0.242–0.645 kg·C/m<sup>2</sup> at 3100–3500 m.


**Table 2.** Carbon mass calculations for *Picea crassifolia* in the literature.

**Figure 2.** Carbon mass proportion at different elevations for twelve meteorological stations. Data are the means of simulated values for each range in elevation.

We investigated the NEE (Net ecosystem exchange) and carbon mass change in the mountain boreal forest ecosystem (at 2700–3100 m where the boreal forest accounted for the largest carbon mass) (Figure 3). In the model, a positive value of NEE indicated that the ecosystem was a carbon source, while a negative value indicated a carbon sink. Boreal forest was a carbon sink before the 1980s; sink strength at elevation 2700–2900 m declined, until the sink activity became a carbon source after 1980, and it became a relatively strong source in the 1990s. However, the strength of the carbon source weakened and the forest ecosystem at 2900–3100 m became a carbon sink again after the year 2000. We analysed the change in carbon mass during 1964–2013 and found the there was an upward trend for this montane boreal forest. The biomass of *Picea crassifolia* forest below 3100 m decreased in the study period, with the greatest decline of 2.595 kg/m<sup>2</sup> at 2700–3100 m, and least of 0.125 kg/m<sup>2</sup> at 2300–2500 m. The carbon mass increased by 0.035 kg/m<sup>2</sup> at 3100–3300 m, and by 0.256 kg/m2 at 3300–3500 m. As the NEE indicated, the decreasing trend of carbon mass for montane boreal forest at 2700–2900 m stabilized in the 2000s and even exhibited a slight increase; carbon mass at 2900–3100 m was stable before the 1990s, and increased after that.

**Figure 3.** The carbon flows change in montane boreal forest ecosystem from 1964–2013. kg/m2/a: kg/m2/year.

#### *3.2. The Sensitivity Analysis of Carbon Mass of Montane Boreal Forest to Climate Variables and CO2*

The simulation showed that the carbon carrying capacity of *Picea crassifolia* increased in the study period, given no changes to any variable. Mean carbon mass increased by 1.202 kg/m<sup>2</sup> from 1964–2013, and the 5-year (simulation time interval for the LPJ-GUESS model) mean NPP increased by 0.023 kg/m2/year to between 0.997 and 1.122 kg/m2/year in the simulation period. Subsequently, each of the climatic variables was changed one at a time.

Following a change in temperature, the trend in NPP did not change, but the NPP value was higher than before (Figure 4). The effect of lower temperature on NPP was greater than that of higher temperature. When the temperature increased by 1 ◦C, 5-year mean NPP increased by 0.090 kg/m2/year; when the temperature increased by 2 ◦C, mean NPP increased by 0.148 kg/m2/year per five years. When the temperature decreased by 1 ◦C, mean NPP decreased 0.113 kg/m2/year per five years; when the temperature decreased by 2 ◦C, mean NPP decreased by 0.292 kg/m2/year per five years. However, the change in carbon stocks was not consistent with NPP when the temperature changed; increased or decreased temperature appeared to suppress carbon stocks of the montane boreal forest. Mean carbon stock decreased between 1964–2013 by 0.768 kg/m<sup>2</sup> per five years; the increments decreased from 12.37% to 6.32% as the temperature increased by 1 ◦C. As the temperature

increased by 2 ◦C, mean carbon stock decreased by 0.893 kg/m<sup>2</sup> per five years, and carbon mass was reduced by 8.24%. With the temperature decrease of 1 ◦C, carbon stock increased by 0.026 kg/m2 per five years and the increments increased to 13.28%. With the temperature decrease of 2 ◦C, mean carbon stock decreased by 1.307 kg/m2 per five years, and carbon mass decreased by 10.11% during the last 50 years.

**Figure 4.** LPJ-GUESS output showing the effects on stored carbon (top panel) and NPP (bottom panel) of changing the temperature only. Temperature was modified by adding or subtracting 1 (or 2) ◦C to each daily climate value. Lines represent the means of the simulations. Temp: the abbreviation of temperature.

The trend in NPP changed over the course of simulation following changes in precipitation (Figure 5); the more precipitation, the higher the NPP, and the less precipitation, the lower the NPP. The impact of reduced precipitation was greater than that of increased precipitation. A 20% increase in precipitation resulted in a mean 5-year NPP increase of 0.024 kg/m2/year, while a 20% decrease reduced NPP by 0.080 kg/m2/year per five years.

The changes in carbon stock were consistent with NPP. Mean carbon stocks increased 0.242 kg/m<sup>2</sup> per five years as precipitation increased by 10%, and 1.165 kg/m2 per five years as the precipitation increased by 20%. As precipitation decreased by 10 and 20%, mean carbon stocks decreased 1.517 and 2.969 kg/m2 per five years, respectively. The increase in carbon mass resulting from increased precipitation was 24.39%, and 24.93%, respectively, (basic value 12.37%), and carbon mass decreased 19.17%, and 34.01%, respectively, if precipitation decreased.

Mean NPP decreased when the CO2 concentration remained unchanged (de-trended) and increased when the CO2 concentration increased (Figure 6). Carbon mass decreased when the CO2 concentration remained unchanged, but it increased as the CO2 concentration levels rose. Overall, the CO2 concentration had an increasingly positive effect as atmospheric levels rose. Carbon mass decreased by 0.240 kg/m<sup>2</sup> for the simulation years if CO2 concentrations did not increase

after 1964. Mean carbon storage increased by 0.746 and 0.866 kg/m2 per five years for CO2 concentration increases of 50 and 100 ppm. De-trended CO2 supressed carbon mass by about 2.45%, and increased CO2 promoted carbon storage from 12.37% to 18.98%, and to 24.74% in response to 50 and 100 ppm, respectively.

**Figure 5.** LPJ-GUESS output showing the effects on stored carbon (top panel) and NPP (bottom panel) of changing precipitation only. Precipitation was changed multiplicatively because it is a zero-based variable. Lines represent the means of the simulations. Prec: the abbreviation of precipitation.

**Figure 6.** LPJ-GUESS output showing the effect on stored carbon (top panel) and NPP (bottom panel) of changing CO2 concertration only. CO2 concertration was de-trended by using 1964 levels throughout and was modified by adding values to each year value. Lines represent the means of simulations.

Indeed, in the simulation period, the NPP and carbon mass of the montane boreal forest increased. Combined with the analysis of single-factor changes, carbon stocks decreased as temperature rose but increased as CO2 concentration increased. The actual temperature increased in the study area from 1964 to 2013, precipitation fluctuated, and CO2 concentration increased, and the positive effect of CO2 concentration contributed more to the carbon stocks than the negative effect of temperature.

#### *3.3. The Sensitivity Analysis of Carbon Mass of Montane Boreal Forest to Extreme Weather*

We analyzed the extremes in weather events from 1964–2013 in stations we chose by ranking the mean value of heat and water in concentration months (annual mean temperature from April to October, and annual mean precipitation from May to September) of each year.

Under warm and wet conditions, NPP increased, and under warm conditions alone, NPP increased more (Figure 7). Under cold and dry conditions, NPP decreased, and it decreased further under dry conditions alone. However, the living biomass was greatly suppressed under dry conditions alone and decreased by 42% compared to the no-change simulation; mean decrease in NPP was 1.816 kg/m<sup>2</sup> per five years. Under the other three conditions, biomass production performed better than under the actual climate. Biomass production increased (10% at the end of period) most notably under wet conditions, with a mean increase of 1.085 kg/m<sup>2</sup> per five years.

**Figure 7.** LPJ-GUESS output showing the effects of extreme weather conditions on stored carbon (top panel) and NPP (bottom panel). Local weather data were used to determine the most extreme high /low temperature and high/low rainfall from a 5-year period.

#### **4. Discussion**

#### *4.1. Model Applicability and Carbon Mass at Different Altitudes*

Our LPJ-GUESS model simulation showed that the main distribution of *Picea crassifolia* in the Qilian Mountains was at 2700 to 3100 m, and the average carbon storage was 10.503 kg/m2, with a maximum of 11.787 kg/m2. The biomass pools at elevations 2300–2700 m contained more carbon than those at 3100–3500 m for the different dominant species.

The vertical distribution pattern of *Picea crassifolia* influenced the distribution of carbon density and NPP. Altitude was a dominant factor influencing the soil organic carbon concentration and community pattern of *Picea crassifolia* [52]; further, altitude exerted a strong influence on the growth of *Picea crassifolia* by affecting the microclimate, including air temperature and humidity, and soil moisture [24]. Trees were low and forest density was often low at low altitudes, increasing gradually at mid-altitudes, and exhibiting a pattern of isolated trees or scattered patches at the upper growth limit [53]. The most suitable conditions for the growth of *Picea crassifolia* forest were observed between 2800 and 2900 m [49], with grassland between 3000 to 3700 m [54]. Density and basal area of *Picea crassifolia* were also higher between 2650 to 3100 m than at other altitudes, and beyond 3100 m, the density decreased with an increase in altitude [55].

Evergreen conifers in cold and high-altitude zones of the montane boreal forest have lower carbon biomass due to low photosynthesis and respiration rates than that in warmer habitats [56]. Carbon mass in the biomass pools examined in this study, i.e., 5.462 kg·C/m2, was similar to the estimate for northern montane boreal forests ranging from 4.2 to 5.3 kg·C/m2 [47]. The focus of current research in the Qilian Mountains is on the main tree species, *Picea crassifolia*, but results are not uniform. Initial surveys at the Sidalong forest farm on the Qilian Mountain showed that *Picea crassifolia* biomass was about 24.298 kg/m<sup>2</sup> (Chang et al., 1995). Recently, the aboveground carbon stocks of *Picea crassifolia* at the north-eastern edge of Qilian Mountains were calculated as 4.3 kg/C·m2 [47]. The estimate of the above-ground biomass in southern Qilian Mountains was between 0.1885 and 22.065 kg/m2 [57]. Our simulation of *Picea crassifolia* NPP was between 0.52–0.58 kg·C/m2/year, which is in the range of 0.3–0.7 kg·C/m2/year for boreal ecosystem productivity examined by remote sensing [58], and higher than the mean NPP (0.38 kg·C/m2/year) of Qilian Mountains [59]. Hydrothermal conditions drove carbon density differences in boreal forests; also, carbon density first increased, and then decreased with stand age, with the highest value at age 183 years [28]. Thus, carbon density calculations were different due to vertical distribution of vegetation patterns and the age of the forest in this study.

The carbon mass and NEE change in *Picea crassifolia* forest ecosystem we calculated revealed that this mountain forest tended to becoming a carbon source from a carbon sink in the last 50 years; this trend was slow and even reversed in the 2000s. The carbon mass of *Picea crassifolia* decreased at low elevations and increased at high elevations, showing an upward trend in vegetation distribution, but the highest carbon mass proportion region for the forest did not change in the study period. Similar to other boreal forests, the mountain boreal forest also tended to become a carbon source under climate change [60–62], but this trend in the Qilian Mountains slowed in the 2000s.

The LPJ-GUESS model has been used in a large number of studies, and its performance has been evaluated several times. The model has also been used in China. Combined LPJ-GUESS and High Accuracy Surface Modeling (HASM) allowed for an economical estimation of forest biomass and the research used climate data from 735 meteorological stations in Chinese mainland from 1950–2010 and compared the results with the Seventh National Forest Resources Inventory data in China [40]. The results presented in this article indicated that the LPJ-GUESS model was suitable for these 735 meteorological stations; stations used in our study were contained in those stations. A carbon balance calculation by LPJ-GUESS of three deciduous forests in a mountainous area near Beijing showed that the model can be applied to a warm temperate forest in China [37]. Also, the forest production and carbon dynamics study of Masson Pine Forest in the Jigongshan region demonstrated that the model can simulate growth dynamics of subtropical forests [39]. Finally, the simulation

of the carbon cycle of a *Larix chinensis* forest at Taibai Mountains, China, indicated the model was suitable for analyzing vegetation characteristics in mountainous areas [38]. Similar to previous biomass calculations (Table 2), we found that the LPJ-GUESS model was suitable for simulating carbon change of the montane boreal forest in northwestern China. The carbon density calculated by the model is based only on foliage, fine roots, and sapwood and heartwood; thus, the results of the model can be expected to be lower than results of some field measurements.

#### *4.2. The Effects of The Climate Variables on Carbon Cycling in Montane Boreal Forest*

Carbon density and NPP in this montane boreal forest ecosystem were concentrated at 2700–3100 m. Mean carbon mass increased by 13.37%, with a 5-year mean of 1.202 kg/m2 from 1964 to 2013; mean NPP increased by 0.023 kg/m2/year per five years.

In this study, NPP clearly increased as temperature increased and under warm conditions. The effect of temperature on carbon storage was complex in this montane boreal forest. Increasing temperature limited the storage of carbon, while decreasing temperature promoted its storage. Under warm conditions, carbon storage increased, and if temperature decreased considerably, carbon mass was reduced. The temperature sensitivity of biomass allocation may be an important but not obvious regulator of the carbon cycle in the boreal forest [19,56]. With an increase in mean annual temperature, and only a modest expected increase in precipitation, a shift in boreal forest may be observed to a woodland/shrubland ecosystem type, which is more suited for such an environment [63]. Further, the overall warming trend may lead to earlier seasonal plant growth; however, incomplete development of green tissues exposed to colder environments may negatively influence tree growth and carbon storage [6]. To the contrary, some researchers concluded that increased temperature and longer growing seasons will increase the NPP and enhance carbon uptake of boreal forests; this response, however, was weak in interannual variations [64]. Boreal forests store vast amounts of carbon in soil, but the temperature effects on soils are challenging to measure. The temperature relationship with above-ground carbon density was affected by baseline temperature such that, at mean annual temperature <8 ◦C, the relationship was positive, and negative in regions with a mean annual temperature >10 ◦C for mature boreal forests [65].

Precipitation is positively related to the biomass and NPP of the montane boreal forest. Precipitation has a much greater impact on montane boreal forests than temperature; especially in reduced precipitation and in dry environments, both biomass and NPP decreased significantly and hindered productivity and biomass storage. Drought conditions likely affected the radial growth of *Picea crassifolia* forest in Qilian Mountains in the last half-century [66]. Drought was the major driver of the release of total original carbon, which reduced soil respiration and NPP in a large boreal watershed [67]. Extreme drought decreases carbon assimilation and reduces the carbon sink strength of forests [6]. Under extremely dry conditions, boreal forests may degrade to low-productivity open woodlands [63].

In this study, the effect of CO2 concentration on the biomass and NPP of montane boreal forest was also consistent, and the parameters were positively correlated. We showed that, if the CO2 concentration remained at the level present 50 years ago, the ability of the montane boreal forest to produce and store carbon would decrease. However, the capacity of forests to store carbon improved with increasing CO2 concentration, and the higher the CO2 amounts, the greater the carbon storage capacity. Although not as pronounced as that in tropical and temperate forests, the increase in CO2 concentration led to an increase in NPP and carbon stocks in boreal forests both at local and global scales [62,68]. A physiologically-based forest model showed that the total NPP and the carbon storage in biomass improved under elevated CO2 concentration due to an increase in net photosynthetic rates at leaf-level and smaller shifts in carbon residence time [56,67]. Using combined satellite and ground observations for 1950–2011, Denos and others (2013) demonstrated that the sequestration of atmospheric CO2 increased due to higher net CO2 uptake associated with spring and fall growth extensions in northern ecosystems. Increasing atmospheric CO2 is expected to increase boreal forest carbon in western China [69]. A simulation by the BIOME-BGC model in *Picea crassifolia* forest showed that the effect of CO2 concentration on NPP was more significant than that of climate change for future climate [70].

#### **5. Conclusions**

The carbon mass of montane boreal forest simulated by LPJ-GUESS was comparable, though lower than that reported in other studies. Carbon storage was concentrated at the altitude of 2700–3100 m in the Qilian Mountains. Within this altitudinal range, carbon storage of *Picea crassifolia* increased by 1.202 kg/m2 per five years, and NPP was between 0.52–0.58 kg·C/m2/year during the last 50 years. *Picea crassifolia* exhibited a trend toward climbing higher in elevation and becoming a carbon source in the last 50 years.

The boreal forest is significantly affected by climate change and has a slow recovery process. Under steady precipitation, carbon storage increased with increasing atmospheric CO2 concentration despite the negative effects of warming on montane boreal forest in the Qilian Mountains in the last 50 years. A lack of water is the greatest threat for carbon storage in *Picea crassifolia*. Under the changing CO2 concentration and precipitation conditions, carbon mass and NPP were positively correlated. The relationship between climate factors, CO2 concentration, extreme conditions, and carbon storage was largely impacted by stand age, geographical location, altitude, and the environment; therefore, further research needs to examine the drivers of interannual variability in the carbon cycle at different scales and under different conditions of climate change.

**Acknowledgments:** We are very grateful to Kathryn Piatek for her comments and editorial assistance. This work was supported by the National Natural Science Foundation of China (No. 41621001, 41522102 and 41601051) and the Foundation for Excellent Youth Scholars of "Northwest Institute of Eco-Environment and Resources", Chinese Academy of Sciences.

**Author Contributions:** Zhibin He and Shu Fang conceived and designed the experiments; Shu Fang performed the experiments, analyzed the data and wrote the paper; Jun Du, Longfei Chen, Pengfei Lin, and Minmin Zhao optimized the experiment and modified the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

#### **References**


© 2018 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/).

### *Article* **Forest Floor and Mineral Soil Respiration Rates in a Northern Minnesota Red Pine Chronosequence**

#### **Matthew Powers 1,\*, Randall Kolka 2, John Bradford 3, Brian Palik <sup>2</sup> and Martin Jurgensen <sup>4</sup>**


Received: 31 October 2017; Accepted: 25 December 2017; Published: 29 December 2017

**Abstract:** We measured total soil CO2 efflux (RS) and efflux from the forest floor layers (RFF) in red pine (*Pinus resinosa* Ait.) stands of different ages to examine relationships between stand age and belowground C cycling. Soil temperature and RS were often lower in a 31-year-old stand (Y31) than in 9-year-old (Y9), 61-year-old (Y61), or 123-year-old (Y123) stands. This pattern was most apparent during warm summer months, but there were no consistent differences in RFF among different-aged stands. RFF represented an average of 4–13% of total soil respiration, and forest floor removal increased moisture content in the mineral soil. We found no evidence of an age effect on the temperature sensitivity of RS, but respiration rates in Y61 and Y123 were less sensitive to low soil moisture than RS in Y9 and Y31. Our results suggest that soil respiration's sensitivity to soil moisture may change more over the course of stand development than its sensitivity to soil temperature in red pine, and that management activities that alter landscape-scale age distributions in red pine forests could have significant impacts on rates of soil CO2 efflux from this forest type.

**Keywords:** carbon cycling; *Pinus resinosa*; soil respiration; stand age

#### **1. Introduction**

Soil respiration represents about 70% of ecosystem respiration in temperate forests [1,2], and includes a combination of respiration from plant roots, mycorrhizae, and microorganisms in the leaf litter, humus, and mineral soil. Variables such as soil temperature [3,4], soil moisture [2,5–7], litter quality and quantity [4,8], and local stand structure [3,9] exert strong controls over soil respiration in forests. Soil temperature and soil moisture may also covary seasonally, or show varying relationships across sites, leading to confounded effects on soil respiration [10,11]. Changes in these variables that occur over the course of stand development could also lead to changes in soil respiration as forests age. Soil temperatures and moisture availability, for instance, may vary between young, open-canopied regenerating stands and stands with a dense, closed canopy representative of the stem exclusion phase of development. As land management agencies begin to incorporate C storage and sequestration into their management goals there is an increasing need to understand how developmental changes influence belowground C cycling.

Soil C cycling can display a large amount of within-system and among-system variability as forests age and develop. Estimates of soil respiration derived from relationships between net primary production and net ecosystem production suggest a characteristic age-related trend in which

the heterotrophic components of soil respiration (micro- and macro-fauna) in temperate forests are expected to decline with age, and the autotrophic portion of soil respiration (plant roots and associated mycorrhizae) are expected to peak in middle-aged stands as net primary production (NPP) and, therefore, substrate supplies peak [12]. Direct measurements of soil respiration in stands of different ages, however, have shown that total soil CO2 efflux (RS) can (1) initially increase with age, peaking in young-intermediate-aged, closed canopy stands before declining with age in mature to older stands [4,13,14]; (2) decrease with age [15]; (3) increase steadily from young to old-growth stands [6,16]; or (4) show opposing age-related trends across the geographic range of the same forest type [17]. This variability in age-related soil respiration trends among different systems underscores the need for ecosystem-specific studies.

While respiration from roots and associated mycorrhizae make up the largest proportion of soil respiration [18], estimates of respiration associated with decomposition of the forest floor layers range from 5–48% of total soil respiration in temperate forests [4,18–22]. There are no clearly established patterns of respiration from the forest floor layers (RFF) across stands of different ages, but RS generally has positive correlations with forest floor mass or thickness [3,15,23]. This suggests that RFF could increase with stand age in forests that are characterized by slow litter decay rates and a resulting increase in forest floor thickness over time. Forests dominated by evergreen conifers, for instance, often accumulate forest floor mass and thickness as they age [24,25], which could lead to increased RFF in older stands. The accumulation of forest floor layers during stand development can also exert a large indirect influence on mineral soil respiration by altering soil temperature and moisture or providing leachate of labile substrate into the mineral soil [26,27].

Understanding how stand age influences C cycling is particularly important because forest management activities have direct impacts on landscape-scale age distributions. Different rotation lengths, for instance, can result in dramatically different abundances of young vs. old stands. This would have significant impacts on soil C losses for systems that show a high degree of variability in soil respiration among age classes. The collection of harvest residues and residual wood as a feedstock for biofuel production can also reduce forest floor and mineral soil C and nutrient stocks [28–30], and litter removal clearly impacts both soil respiration and other soil processes [20,31].

We conducted an experiment to characterize age-related differences in total soil respiration (RS) and respiration from the forest floor layers (RFF) in red pine (*Pinus resinosa* Ait.) stands aged 9–123 years. Productivity in managed red pine systems peaks between 130 and 140 years across a range of basal areas [32], so we expected RS to increase steadily with age across our chronosequence due to increases in autotrophic respiration linked to increasing NPP [33]. We also predicted that RFF would increase with age across the chronosequence because forest floor C (and thus, forest floor mass) increases steadily with age up to at least 150 years in red pine forests [25], and forest floor mass is positively correlated with RS [3,15,23].

#### **2. Materials and Methods**

#### *2.1. Study Area*

Our study sites included four red pine stands on the Chippewa National Forest in northern Minnesota, USA. The study sites included an open, 9-year-old stand of red pine saplings that were planted following the clearcutting of the previous red pine stand (Y9), a 31-year-old plantation in the early stages of stem exclusion (Y31), a 61-year-old plantation with a well-developed shrub layer and some tree recruitment in the understory (Y61), and a 123 year-old naturally-regenerated stand with a well-develop shrub and understory tree layer (Y123). Three of the stands were located in close proximity (within 1000 m of one another), and the fourth (Y61) was located 5 km to the west.

Forest floor and soil characteristics, live tree density, and basal area varied somewhat by stand age (Table 1). Red pine represented 94–98% of the basal area of stems ≥2.5 cm in diameter in Y31, Y61, and Y123, with eastern white pine (*Pinus strobus* L.), paper birch (*Betula papyrifera* Marsh.), balsam fir

(*Abies balsamea* L.), and northern red oak (*Quercus rubra* L.) each representing anywhere from 0–3.4% of the remaining basal area within these stands. Y9 had very little basal area of stems ≥2.5 cm in diameter. Hardwoods represented < 2% of basal area in any stand. Beaked hazel (*Corylus cornuta* Marsch.), balsam fir, juneberries (*Amelanchier* spp.), and blueberries (*Vaccinium* spp.) were the dominant understory shrubs and saplings, with a much more developed shrub, sapling, and herbaceous layer in the young, open-canopied Y9 than the other stands.

Study sites were selected to minimize variability in edaphic and physiographic characteristics that could contribute to differences in the soil environment or local microclimate. All stands were located on coarse-textured, excessively drained outwash sands classified as Menahga mixed, frigid Typic Udipsamments to minimize variability in physical soil characteristics that could contribute to differences in soil respiration [34]. Topographic variation across the study sites was very limited with less than 10 m variability in elevation among sites, slopes less than 10%, and plots within each site located on predominantly south or southeast facing aspects. The study area has an average annual temperature of 3.9 ◦C and average annual precipitation of 700 mm.

**Table 1.** Forest floor, soil (0–30 cm), and overstory structure characteristics of 9, 31, 61, and 123 year-old red pine stands used to study age-related changes in soil respiration, forest floor respiration, and litter decomposition.


FF = forest floor; TPH = trees ha<sup>−</sup>1; pH A = soil pH from the surface to the top of the A horizon; pH A-30 = soil pH from the top of the A horizon to a depth of 30 cm; all pH values measured in a CaCl2 solution. Forest floor and soil %C and %N measured on a LECO Total Elemental Analyzer. Note that %C and %N data were not ash-corrected. Forest floor mass is an estimate based on a companion study in the same stands [35] and mass from the actual collars in this study.

#### *2.2. Field and Laboratory Methods*

We selected a chronosequence of 4 red pine stands at 9, 31, 61, and 123 years of age. We placed 25 cm (inside diameter) PVC (polyvinyl chloride) collars into the soil in November of 2008 after leaf fall of deciduous trees and shrubs (approximately 6 months prior to our first sampling). Nine collars were placed at three randomly located sample plots in each stand, for a total of 27 collars per stand. The nine collars in each plot were placed in groups of three at a distance of four meters from plot center along transects centered at azimuths of 60, 180, and 300 degrees at each plot.

Since the presence of leaf litter can exert a large indirect influence on mineral soil respiration by altering soil temperature and moisture or providing leachate of labile substrate into the mineral soil [26,27], we used forest floor manipulation treatments to examine direct and indirect impacts of forest floor removal on RS and estimates of RFF. In this context, "direct" implies the additive contribution of CO2 efflux from litter decomposition to total RS, while "indirect" implies the effects that modification of the soil environment due to the presence of an intact litter layer could have on both autotrophic and heterotrophic contributions to RS from the mineral soil. In April of 2010, we removed the forest floor (down to mineral soil) from one randomly selected collar at each transect (the "no litter" treatment) to evaluate how complete forest floor removal impacts soil respiration. Since complete forest floor removal could create indirect impacts on soil respiration by altering the soil environment, we also removed the forest floor from a second PVC collar at each transect, and replaced it with a removable, 20-mesh (841 micron) bag (filled with litter collected from an area near each transect of collars equal in size to the cross-sectional area of the collars themselves. In this treatment (the "removable bag" treatment) the mesh bags were removed immediately prior to taking soil respiration measurements, with the intention of maintaining mineral soil conditions similar to those under intact

litter while removing the direct contribution of CO2 efflux from the forest floor from estimates of soil respiration. The third collar in each transect was left as an untreated control.

Total soil CO2 efflux (RS) was measured on a monthly basis on all 27 collars in each stand from May through October in 2009 (pre-manipulation; i.e., prior to litter removal and removable bag installation) and from April through October in 2010 (post-manipulation; i.e., after litter removal and removable bag installation) using a custom-built, closed-chamber system with an LI-820 gas analyzer (LI-Cor, Inc., Lincoln, NE, USA) connected to a datalogger programmed to record the CO2 concentrations in the system at 2 s intervals. Equipment malfunctions prevented us from completing the August 2009 measurement.

Since our closed-chamber system calculates CO2 flux rates by first scrubbing CO2 out of the chamber and then measuring the rate of increase as chamber CO2 rises back to a pre-set upper boundary, measurement times on individual collars varied based on instantaneous flux rates. RS measurements were conducted over a 2–3 day period each month, with all collars in each stand being sampled on a single day in each month. Individual stands were sampled at alternating the times of day (0930–1200 or 1200–1500) from one month to the next in an effort to average out variability associated with taking measurements at different locations at different times of day.

Temperature readings from the top 10 cm of soil suggest soils were frozen from mid-late November until mid-April during both years of measurement, so we assume that soil respiration was minimal during the November to April period. All fresh litter inputs were removed from the no litter treatments on a bi-weekly to monthly basis throughout our 2010 measurements, and respiration measurements in the removable bag treatment were made five to ten minutes after removing the bags from the collars to allow the CO2 environment within the collar to equilibrate after removing the bag. We measured mineral soil temperature at a depth of 10 cm next to each collar, and measured volumetric soil moisture content in the top 12 cm of the mineral soilusing a HydroSense Soil Water Measurement System consisting of a CS620 sensor attached to a CD620 display (Campbell Scientific, Inc., Logan, UT, USA). Soil moisture measurements were taken at four points around each collar during each soil respiration measurement. After litter treatment in 2010, temperature and moisture measurements for the no litter and removable bag treatments were made within similar treatment areas next to each treatment collar to capture treatment effects without disturbing the soil environment within the collars used for respiration measurements.

We estimated RFF based on the differences of RS measurements between the control and treatment collars as Equations (1) and (2):

$$R\_{\text{FF NO LITTER}} = R\_{\text{S CONTROL}} - R\_{\text{S NO LITTER}} \tag{1}$$

and

$$\text{R\\_F\\_REMOVALE\\_BAG} = \text{R\\_CONTROL} - \text{R\\_REMOABLE\\_BAG} \tag{2}$$

where RFF is respiration from the forest floor layers, RS CONTROL is the measured respiration from the control collar in a transect, RS CONTROL is the measured respiration from the "no litter" treatment collar in that transect (i.e., collars that received a complete litter removal treatment), and RS REMOVABLE BAG is the measured respiration from the "removable bag" collar in that same transect (i.e., collars in which the litter was removed from the collar and replaced with a mesh bag filled with litter that could be left in place between measurements and removed immediately prior to measurements).

#### *2.3. Analytical Methods*

We analyzed litter treatment and stand age effects on RS, RFF, soil temperature, and soil moisture using linear mixed models. Models included litter treatment, stand age, month of measurement, and the interactions among these three factors as fixed effects along with random effects to account for the spatial nesting of litter treatments within plots and the temporal nesting of measurements within plots. We used averages of litter treatments or litter types from the three transects at each plot as the

dependent variable in each analysis. Respiration data, temperature data, and moisture data from 2009 and 2010 were analyzed separately because the measurements spanned different months in the two years. Bonferroni-adjusted t-tests were used to evaluate significant differences between treatment levels. Model assumptions were evaluated with residual plots.

We analyzed relationships between RS and soil temperature based on the exponential growth function (Equation (3)):

$$\text{Rs} = \beta\_1 \exp(\beta\_2 \times T) \tag{3}$$

where RS is soil respiration and T is temperature. To include nested random effects (which do not have built-in support for mixed-effects nonlinear models in SAS version 9.2 (SAS Institute, Cary, NC, USA)) to account for repeated measurements on collars representing different treatments within plots, we log-transformed the exponential growth function, and analyzed the soil respiration—temperature relationship using linear mixed models of the form (Equation (4)):

$$\text{LN}(\text{R}\_{\text{S}}) = \text{LN}(\beta\_1) + (\beta\_2) \times \text{T} \tag{4}$$

where LN refers to the natural logarithm, RS is soil respiration, and T is soil temperature. To test for potential stand age and litter treatment effects on the relationship between soil respiration and temperature, we fit models that allowed β<sup>1</sup> and β<sup>2</sup> to vary both individually and in tandem across the levels of both stand age and litter treatment. We used AIC scores to choose the model with the best fit to our data, but also considered the *p*-values of individual model terms to judge the significance of age and litter treatment effects. All models relating RS to temperature included random effects to account for the use of repeated measurements and the nesting of treatments within plots in each stand. We estimated the sensitivity of soil respiration to temperature using Q10 values. We calculated Q10 as Equation (5):

$$\mathbf{Q}\_{10} = (\mathbf{R}\_2/\mathbf{R}\_1)^{(10/(\mathbf{T}\_2 - \mathbf{T}\_1))} \tag{5}$$

where R1 and R2 are modeled soil respiration rates at temperatures T1 and T2, respectively.

We also used linear mixed models to analyze relationships between soil respiration and soil moisture. All soil moisture models followed the general form (Equation (6)):

$$\mathcal{R} = \beta\_1 + \beta\_2 \times \text{SM} \tag{6}$$

where R is either an absolute, or one of two normalized measure of soil respiration and SM is volumetric soil moisture. Thus, we ran three sets of models relating soil respiration to soil moisture including (1) one set including our direct, field-based RS estimates as the dependent variable; (2) a second set using a normalized expression of soil respiration (we refer to this as RSN) calculated as the ratio of observed RS to the value predicted by our best-fitting temperature model based on soil temperatures at the time of field RS sampling; and (3) a third set expressing soil respiration as the difference between our observed value from field measurements of RS, and the value predicted by our best-fitting temperature model based on soil temperatures at the time of field RS sampling (we refer to this variable as RSDIF).

Like our respiration-temperature models, we fit respiration-moisture models that allowed β<sup>1</sup> and β<sup>2</sup> to vary by both individually and in combination across the levels of stand age and litter type, then used AIC scores to choose the best model. We used residual plots to analyze model assumptions regarding the normality and homogeneity of error variances for all of our regression models, and applied transformations when necessary. We used Bonferroni-adjusted *t*-tests for all multiple comparisons. We calculated a conditional R2-type goodness of fit value for the respiration—temperature and respiration—moisture models as Equation (7):

$$\mathbf{R}^2 \llcorner = 1 - \left[ \Sigma \left( \mathbf{y}\_{\text{i}\overline{\mathfrak{y}}} - \mathbf{\hat{y}}\_{\text{i}\overline{\mathfrak{y}}} \right)^2 / \Sigma \left( \mathbf{y}\_{\text{i}\overline{\mathfrak{y}}} - \overline{\mathfrak{y}} \right)^2 \right] \tag{7}$$

where yij and yˆij are the observed and predicted respiration rates for each individual subject and y is the overall mean of the observed respiration rates [36]. All statistical analyses were performed using SAS version 9.2 (SAS Institute, Cary, NC, USA) at α = 0.05 significance level.

#### **3. Results**

#### *3.1. Total Soil Respiration*

Litter treatment had significant impacts on RS (*p* = 0.010 for post-manipulation measurements), and stand age (*p* < 0.001 and *p* = 0.001), month of measurement (*p* < 0.001 for both years), and the stand age *x* month of measurement interaction (*p* < 0.001 for both years) were also significant in the pre-manipulation and post-manipulation measurement years. There were no significant interactions involving litter treatment in either year of measurement. There were no significant differences in RS among collars assigned different litter treatments prior to litter treatment (*p* = 0.285), but RS was significantly lower for the no litter treatment than for the control treatment after litter treatment (Table 2).

**Table 2.** Effects of forest floor removal treatments on total soil respiration (RS), forest floor respiration (RFF) (calculated after litter treatment in 2010 only), and soil moisture (SM) before (2009) and after (2010) experimental litter removal treatment in red pine stands. Different letters indicate significant differences among treatments within a column. Respiration rates reported in μmol CO2 m−<sup>2</sup> s<sup>−</sup>1.


Averaged across months, RS was highest in Y9 and lowest in Y31 in 2009, and generally lowest in Y31, but similar among other stands in 2010 (Table 3). However, the age effect was variable across the individual months of measurement (Figure 1). From May through September of 2009 (Figure 1e) and June through August of 2010 (Figure 1f), RS was generally lowest in Y31. There were no differences in RS among stands in October of 2009, or in April, May, or October of 2010, when soil temperatures were at their lowest for each year. Y61 had higher RS than any other stand in September of 2010, but was generally similar to the Y9 and Y123 stands in other months. Monthly variability in RS generally paralleled soil temperature trends in both years. In 2009, RS monthly patterns of Y61 and Y123 closely followed the temperature trend, although RS in Y9 and Y31 began declining in July while soil temperatures remained high. In 2010, RS patterns closely followed the monthly temperature trends, showing a steady increase through August followed by a decline in September and October.

**Table 3.** Mean values of total soil respiration (RS), volumetric soil moisture (SM), soil temperature (T), and forest floor respiration (RFF) before experimental litter removal (2009) and after litter removal (2010) in red pine stands aged 9 (Y9), 31 (Y31), 61 (Y61), and 123 (Y123) years. Values represent annualized averages across multiple months of measurement in each year. RFF was calculated only after litter removal (i.e., only in 2010). Different letters indicate significant differences among treatments within a column. Respiration rates reported in μmol CO2 m−<sup>2</sup> s<sup>−</sup>1.


**Figure 1.** Soil temperatures (a, b), soil moisture (c, d), and soil respiration (RS, e, f) from two years of measurement in 9, 31, 61, and 123 year-old red pine stands in northern Minnesota, USA. Error bars indicate standard error. An asterisk indicates a month with significant differences among stand ages.

#### *3.2. Forest Floor Respiration*

Our estimates of RFF were more than five times higher for the no litter treatment than for the removable bag treatment (*p* = 0.044, Table 2). RFF was not significantly affected by stand age (Table 3) or month of measurement, and there were no significant interactions between stand age, month of measurement, and litter treatment type (i.e., RFF LITTER vs. RFF REMOVABLE BAG). RFF represented 13.1% of RS in the no litter treatment compared to 3.9% of RS in the removable bag treatment (*p* = 0.018, Table 2). This percentage varied by month of measurement (*p* < 0.001), but the month effect was not consistent across stands of different ages (*p* = 0.027). RFF generally represented the smallest percentage of RS in April when soil temperatures were lowest across all stands, but differences among other months were highly variable across stands (Figure 2). There were no significant interactions involving litter treatment, so monthly values displayed in Figure 2 have been expressed as the mean of RFF estimates calculated using the removable bag treatment and the mineral soil treatment.

**Figure 2.** CO2 efflux from forest floor layers (RFF, top panel) and the percentage of total soil respiration (RS) represented by the forest floor (bottom panel) in 9, 31, 61, and 123 year-old red pine stands in northern Minnesota, USA. Error bars indicate standard error. An asterisk indicates a month with significant differences among stand ages. RFF estimates displayed here reflect the mean of values calculated from two separate treatments (complete litter removal and temporary removal of litter in a mesh bag).

#### *3.3. Soil Temperature and Moisture*

Litter treatment did not have a significant effect on soil temperature in either pre- or post-treatment measurements (*p* = 0.999 and 0.372, respectively), but stand age (*p* < 0.001 for both years) and month of measurement (*p* < 0.001 for both years) each did, and the stand age effect varied across months (*p* < 0.001 for both years). There were no significant interactions involving litter treatment in either year of measurement. Averaged across months, soil temperatures were generally highest in Y9, lowest in Y31 and Y123, and intermediate to high in Y61 (Table 3), but these patterns changed in October, when soil temperatures were lowest in Y9 and highest in Y61 (Figure 1). The temporal trends of soil temperature were somewhat different from 2009 to 2010. In 2009, soil temperatures increased from May to June, were similar in June, July, and August, and then declined rapidly (Figure 1a). In 2010, soil temperatures rose steadily from April through August then declined in September and October (Figure 1b). Soil temperatures in the spring of 2009 were somewhat higher than those during the spring of 2010 but the reverse was true during late summer and autumn.

Litter treatment had a significant impact on soil moisture in 2010 (*p* < 0.001), and stand age, month of measurement, and the stand age *x* month of measurement interaction were each significant in both the pre-manipulation and post-manipulation years of measurement (*p* < 0.001 for all). There were no significant interactions involving litter treatment in either year of measurement. Although there were no differences between collars assigned different litter treatments prior to manipulation in 2009 (*p* = 0.174), soil moisture was higher in the no litter treatment than in the removable bag and control treatments after litter treatment in 2010 (Table 2). In 2009, soil moisture was typically highest in Y9, and similar among other treatments (Table 3, Figure 1c). In 2010, soil moisture was similar among stands during April and October, highest in Y9 during May and August, and generally low in Y61 during July and September (Figure 1d).

#### *3.4. Relationships between Soil Respiration, Temperature, and Moisture*

There was a significant, exponential relationship between RS and temperature (*p* < 0.001, R2 <sup>C</sup> = 0.787), but stand age and litter treatment did not have significant effects on the RS—temperature relationship (Figure 3). There was also a significant, logarithmic relationship between RS and soil moisture (*p* < 0.001), but model predictions were not well correlated with field observations (R2 <sup>C</sup> = 0.083). Additionally, the shape of the RS—moisture relationship varied with stand age (*p* < 0.001 for both soil moisture and the soil moisture *x* stand age interaction). Our model predicted positive relationships between RS and moisture in Y9 and Y31, but negative relationships between these two variables in Y61 and Y123 stands (Figure 4).

**Figure 3.** Relationships between soil respiration (RS) and soil temperature in red pine stands in northern Minnesota, USA.

**Figure 4.** Relationships between soil respiration (RS) and soil moisture in 9 (Y9), 31 (Y31), 61 (Y61), and 123 (Y123) year-old red pine stands in northern Minnesota, USA.

RSN was not significantly correlated with soil moisture (*p* = 0.407), and age, litter treatment, and the interactions among these three variables did not have significant effects on RSN in any of the models we tested (Figure 5). Further, no RSN models including any combination of soil moisture, age, litter treatment and their interactions performed better than the null model for RSN that included only the intercept and mixed-effects terms accounting for the spatial nesting of collars within plots and the repeated measurements of individual collars and plots over time. RSDIF did have a significant, negative relationship with soil moisture (*p* = 0.0137), but stand age, litter treatment, and their interactions were never significant, and no models including stand age or litter treatment as fixed effects improved on the fit of the general RSDIF to soil moisture model (Figure 5). Further, the correlation between RSDIF and soil moisture was extremely weak (R<sup>2</sup> <sup>C</sup> = 0.009).

**Figure 5.** Relationships between two normalized estimates of soil respiration (RS) and soil moisture in red pine stands in northern Minnesota, USA. RS Observed (top panel) indicates direct, field-based measurements of soil CO2 efflux, and RS modeled (bottom panel) indicates estimates of soil CO2 efflux derived from our best-fitting regression model relating RS to temperature. RSDIF was calculated as the differences between RS Observed and RS Modeled.

#### **4. Discussion**

Our results suggest that stand age may have significant impacts on soil carbon cycling in red pine ecosystems, but stand age had different effects on RS and RFF. Our young, closed-canopy stand in the early stages of stem exclusion (Y31) showed a distinct seasonal pattern of RS compared to a young stand regenerating after a clearcut, a mature stand, and an old stand. In contrast, stand age had little impact on RFF or RFF's contribution to RS. Our results also suggest that the forest floor influences soil respiration both directly, through the contribution of CO2 efflux associated with the decomposition of litter material, and indirectly, by altering the physical environment within the upper mineral soil.

#### *4.1. Stand Age Effects on Soil Respiration*

Our results do not support the increasing trend of RS that was predicted based on linkages between productivity and respiration [33]. Although theory predicts a possible decline in RS from middle-aged to older forests [12], several studies have reported positive relationships between stand age and soil respiration rates spanning various stages of forest development [6,13,16,17]. We found no evidence of either of these patterns, but we did find that soil respiration rates were often lower in Y31 than in the other stands during the middle of the growing season, particularly in July and August when soil temperatures were highest. Y31 often had lower soil temperatures than the other stands, particularly Y9 and Y61. As other studies have indicated [10,19,37], we found that RS increased exponentially with increasing soil temperature, so low soil temperatures are a likely cause of the low soil respiration rates in Y31. In our study system, Y31 is representative of plantations in the stem exclusion stage of development. Dense canopies that severely limit understory light availability are common characteristics of this phase [38], which could explain the low soil temperatures we observed in Y31. Although we did not directly measure either understory light availability or overstory leaf area index in these stands, we anecdotally observed that Y31 had very limited understory forb and shrub development in comparison with the other sites, which could be indicative of lower light levels on the forest floor, and thus, lower soil temperatures.

Although the RS-temperature relationship appeared relatively stable across our stands, we did find evidence that the RS-soil moisture relationship varied among different-aged stands. Some studies have reported that soil respiration's sensitivity to temperature can change with stand age or successional status [13,16,19], but few have reported different age-related trends in the RS-moisture relationship. Like Irvine and Law [6], we found that soil respiration rates were generally higher in our older stands than in our younger stands during periods of low soil moisture. As soil moisture increased, however, soil respiration rates increased in Y9 and Y31 stands, but declined in Y61 and Y123. While the cause of declining RS as soil moisture increased in Y61 and Y123 stands is unclear, the effect (i.e., decreasing differences in RS between young and older stands as soil moisture increases) is consistent with Irvine and Law's [6] proposal that the larger trees in older stands may avoid drought impacts on metabolic processes by accessing water in deeper soil horizons. The possibility of reduced drought impacts on metabolic processes in older stands is also consistent with findings from foliar gas exchange studies that suggest larger trees in older stands show less physiological response to drought than smaller trees in younger stands [39–41]. Thus, the age-related trends we observed for the RS-moisture relationship are likely driven by more pronounced physiological responses (including respiratory functions associated with growth and maintenance) to drought in younger stands. However, we caution that relationships between RS and soil moisture in our dataset were extremely weak, so soil moisture availability was not likely a major driver of RS variability in our system. Soil moisture typically exerts strong controls on RS in forests only when soils are very dry or very wet [10,42], so it is possible that the 5–15% range of volumetric soil moisture content observed in this study simply did not produce conditions that would limit either decomposer activity or root respiration.

Although we found only limited evidence of age-related variability in the RS-moisture relationship, seasonal patterns of RS appeared to closely track soil temperature trends, regardless of stand age. The plateaued pattern of RS characteristic of most stands in 2009 and the late summer peak in RS observed in all stands in 2010 paralleled the monthly soil temperature variability for those two years, and the high RS values observed in Y61 during September 2010 occurred during a period with higher soil temperatures in that stand than in any other stand. In contrast the seasonal dips in soil moisture observed in September of 2009 and June of 2010 did not appear to have a large influence on RS and the correlation coefficient for our RS-moisture relationship was an order of magnitude lower than that of our RS-temperature relationship. Annualized averages of RS across months of measurement in both years also closely tracked stand-level average soil temperatures, and the generally low soil temperatures in the dense-canopied Y31 appear to be a primary driver of this stand's low RS, particularly when compared with the young, open-canopied Y9, which had considerably higher average soil temperatures.

Additionally, our two normalized estimates of RS (i.e., RSN and RSDIF), which expressed observed RS relative to predictions from our RS-temperature model, showed little to no relationship with soil moisture. This suggests that some of the apparent soil moisture impacts on RS across stand of different ages that we observed were likely driven by covarying temperature trends as was observed in other studies e.g., [10,11]. Collectively, these results underscore the importance of soil temperature as a driver of RS variability e.g., [3,4], and suggest that variations in stand structure that influence soil temperatures (such as the transition from an open-canopy during cohort establishment or stand initiation phase to a dense, closed-canopy during stem exclusion) could potentially exert significant controls on soil carbon cycling.

Although we did not directly test the relationships, variability in soil temperature paralleled differences in temperatures across the two years, while soil moisture generally tracked precipitation levels in the wetter 2010 sampling period, but not as closely in the drier summer of 2009 (Table 4). That both soil temperatures and soil respiration followed the summer plateau pattern of air temperatures in 2009 and the gradual rise in air temperatures through August of 2010, regardless of stand age, speaks to the strong controls that climate exerts on soil processes.


**Table 4.** Mean monthly temperature (T) and precipitation (P) data for our study area during the months of RS measurement in 2009 and 2010 [43].

#### *4.2. Forest Floor Contributions to Soil Respiration*

Our results indicate that forest floor removal reduced RS by 4–13% over the course of our April through October sampling (when averaged across RFF estimates from the no litter and removable bag treatment), but do not suggest the anticipated age-related effects on either RFF or the forest floor's contributions to RS. While several studies have found that respiration from organic horizons represents a larger percentage of RS than our estimates [18–20,22,44], our results fall within the range reported by others [4,22,24], and support the general consensus across these studies that root respiration and decomposition of organic compounds within the mineral soil represent the largest contributions to RS.

The reasons for the relatively low forest floor contributions to RS found in our study are not clear, although substrate quality could be a factor. Like the mineral soil, respiration rates in the forest floor are negatively correlated with the C:N ratio of the substrate [4], and three of our four stands had relatively high forest floor C:N ratios. Forest floor C:N ratios ranged from 39–45 in our Y31, Y61, and Y123 stands. In contrast, forest floor C:N ratios range from 31–32 in young and old-growth Appalachian hardwood forests [4], 34–35 for Douglas-fir (*Pseudotsuga menziesii* Franco) stands in the western United States [24], 29–32 in young to middle-aged Chinese hardwood forests [45], and 21–27 for various components of the forest floor in a tropical montane cloud forest [41]. While the forest floor C:N ratio of our youngest stand was somewhat lower than the other three stands, this stand had a much smaller forest floor layer than the older stands. Forest floor thickness is positively correlated with respiration rates [3,15], so the thinner forest floor may have offset the effects of higher quality litter in Y9, which could explain why RFF estimates were similar across stands despite the difference in forest floor C:N ratios.

The result that RS in the removable bag treatment was similar to RS for the control treatments suggests that at least some portion of the difference between control and no litter treatments was due to the influence that forest floor layers exert on the mineral soil environment [31]. Soil temperatures, however, were similar among the three litter treatments. Soil moisture was higher in the no litter treatment, but soil moisture generally has a positive relationship with RS [5,7], and had little to no

relationship with temperature-normalized estimates of RS in our study, so increased soil moisture is not a likely explanation for lower respiration rates in the no litter treatment. If anything, increased soil moisture in the no litter treatment may have compensated in part for reductions in RS associated with the absence of a decomposing forest floor. If the difference in soil respiration rates between no litter and control collars was due entirely to the absence of respiration contributions from the forest floor, however, we would expect a similar difference to be apparent in the removable bag treatment. The apparent importance of indirect forest floor impacts on RS could suggest that RFF and mineral soil respiration estimates from studies that do not account for the forest floor's environmental influences may contain considerable error, at least in systems with relatively thick forest floor layers, such as the red pine forest studied here. The lack of soil temperature differences between the no litter and removable bag treatments, and limited impact of soil moisture on RS observed at our study sites, however suggest that indirect controls of the soil respiration other than soil microclimate drove the differences in RFF estimates between the no litter and removable bag treatments that we observed.

Negative estimates of RFF during some individual measurement periods were an artifact of our RFF estimation technique, which relied on the difference in soil CO2 efflux measurements between control collars and adjacent collars where litter was removed. Although efforts were made to place adjacent collars in similar environments and adjacent collars generally had similar soil environments (as measured by soil temperature and moisture), differences in autotrophic contributions to soil CO2 efflux associated with variable live root biomass beneath collars may have contributed error to our RFF estimates. However estimates of the contribution of RFF to total RS in a study that used similar methods of estimating RFF while also controlling for root respiration via trenching were broadly similar to our estimates [18], Additionally, litter contributions to total RS have been shown to be very low during the pre-growth and pre-dormancy period [22], which suggests that even small differences in autotrophic respiration could have contributed to the generally negative estimates of RFF for our earliest and latest measurement periods.

#### **5. Conclusions**

Our results have some important implications about C cycling in red pine ecosystems. First, soil respiration responses to changes in temperature appear to be relatively constant across stand ages in red pine. Differences in soil temperature also appear to contribute to low soil respiration rates in dense, young, closed-canopy red pine plantations when compared to open-canopied regenerating stands and older stands that have transitioned out of the stem exclusion phase of development. This suggests that changes in rotation lengths for red pine management could impact soil C dynamics by modifying the proportional representation of age classes across the landscape and consequently, shifting the relative abundance of stand structures with relatively low soil respiration rates (e.g., our young, closed canopy plantation) vs. those with higher, or more seasonally-variable soil respiration rates (e.g., our young, open-canopied stand and our more mature stands).

**Acknowledgments:** We would like to thank Doug Kastendick, Nate Aspelin, Eric Doty, and Ryan Kolka for their assistance in locating field sites, installing respiration collars and litter bags, collecting soil respiration data, and processing litter samples in the laboratory. Funding for this project was provided by the United States Department of Agriculture Forest Service Northern Research Station and the Agenda 2020 Program. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the United States Government.

**Author Contributions:** R.K., J.B., B.P. and M.J. conceived and designed the experiments; M.P. performed the experiments and analyzed the data; M.P. wrote the paper with significant feedback from R.K., J.B., B.P. and M.J.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Law, B.E.; Ryan, M.G.; Anthoni, P.M. Seasonal and annual respiration of a ponderosa pine ecosystem. *Glob. Chang. Biol.* **1999**, *5*, 169–182. [CrossRef]


© 2017 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/).

## *Article* **Carbon Sequestration in Protected Areas: A Case Study of an** *Abies religiosa* **(H.B.K.) Schlecht. et Cham Forest**

**Pablo I. Fragoso-López 1, Rodrigo Rodríguez-Laguna 2, Elena M. Otazo-Sánchez 1, César A. González-Ramírez 1, José René Valdéz-Lazalde 3, Hermann J. Cortés-Blobaum <sup>1</sup> and Ramón Razo-Zárate 2,\***


Received: 23 August 2017; Accepted: 28 October 2017; Published: 12 November 2017

**Abstract:** The effects of global climate change have highlighted forest ecosystems as a key element in reducing the amount of atmospheric carbon through photosynthesis. The objective of this study was to estimate the amount of carbon content and its percentage capture in a protected *Abies religiosa* forest in which the study area was zoned with satellite image analysis. Dendrometric and epidometric variables were used to determine the volume and increase of aerial biomass, and stored carbon and its capture rate using equations. The results indicate that this forest contains an average of 105.72 MgC ha<sup>−</sup>1, with an estimated sequestration rate of 1.03 MgC ha−<sup>1</sup> yr<sup>−</sup>1. The results show that carbon capture increasing depends on the increase in volume. Therefore, in order to achieve the maximum yield in a forest, it is necessary to implement sustainable forest management that favors the sustained use of soil productivity.

**Keywords:** climate change; protected forest; carbon sequestration; *Abies religiosa*

#### **1. Introduction**

The negative impacts of anthropogenic emissions of greenhouse gases such as irreversible damage to ecosystems, increased pressure on water resources, alterations in food production, and damage to human health, among others, have been reported in different studies [1–10]. The need to stabilize the carbon content of the atmosphere has been manifested in a series of international and local agreements and policies, such as the Kyoto Protocol and the Treaty of Paris. The purpose of these agreements and policies is to reduce emissions of greenhouse gases (GHG), with mechanisms to optimize carbon sinks.

Currently, forests store about 800 gigatons of carbon (GtC) [11] and it is estimated that by 2050 they could sequester up to an additional 87 GtC [12,13]. It was estimated that in the period between 2000 and 2007, the carbon sequestration rate of the world's forests averaged 4.1 GtCyr−<sup>1</sup> [14], corresponding to approximately 30% of fossil fuel emissions in 2010 [15].

Globally, protected forests have been proposed as a potentially cost-effective strategy to counter deforestation and degradation [16–18], favoring carbon permanence in the forest. Countries with the greatest threats from their forests due to degradation and devastation have increased their percentage of protected areas (Figure 1) in an attempt to conserve the environmental services of their forests [19,20]. Out of 3984 million hectares of forests in the world, 13.25% have a protected area status [21], and this percentage is mainly because of many of these protected sites partially fulfilling their conservation objectives [22], primarily derived from the budgetary constraints in which most of these areas operate. Financial resources managers for protected areas are increasingly emphasizing cost-effective aspects such as ecosystem services, including carbon sequestration [23].

**Figure 1.** Percentage of protected area per country in 2015 in relation to its total forest area. Own elaboration with information of [24].

Currently, decision-makers can use a large number of methods to assess protected area management and prioritize investments, and there are over 70 methods that have been developed to provide standards for obtaining indicators, analyses, and interpretations [25,26]; however, none of these provide an estimation of carbon sequestration rates.

The present study aimed to estimate the carbon content in the above-ground biomass (baseline) of a forest of *Abies religiosa* (H.B.K.) Schlecht. et Cham., which is part of El Chico National Park, Hidalgo, Mexico. In addition, it attempted to ascertain the carbon capture rate from the annual volumetric increase of the species in the area. This information is necessary for designing the strategies that allow the environmental objectives that contribute to the reduction of the negative effects in climate change to be fulfilled.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study was carried out in the federal zone of the protected natural area called "El Chico National Park" located in the western area of the mountain range of Pachuca, in the state of Hidalgo, Mexico. Geographically, it is located between the extreme coordinates 20◦11 57 and 20◦12 02 north latitude and 98◦43 08"and 98◦43 06" west longitude.

The area is owned by the Mexican nation and comprises 1,833,000 hectares [27]. The climate is C (m) (w) b (i ) gw, which is described as: a temperate-sub-humid climate with a fresh and long summer; the average annual temperature varies between 12 and 18 ◦C. The rainfall regime is in summer, and the percentage of winter rain in relation to the annual total is less than 5% [28].

The soil type is constituted by associations that group the following soil units: Humic cambisol-Androsol ocrico-Litosol, which has a volcanic origin association typical of mountainous zones, and humic Andosol-Humic cambisol, corresponding to forest soils associated with *Abies* and *Quercus* forests [29].

The predominant type of vegetation is *Abies religiosa* forest, which covers 67% of the park surface. The shrub stratum is dominated by the species *Archibaccharis hieracioides* Blake, *Baccharis conferta* H.B.K., *Eupatorium hidalgense* Rob., *Fuchsia thymifolia* H.B.K., *Ribes affine* H.B.K., *Salvia elegans* Vahl., *Senecio angulifolius* D.C., and *Stevia monardifolia* H.B.K. [30].

#### *2.2. Description of the Species*

The *Abies* forest is considered a component of the Leopold boreal forest [31] because of its similarity in terms flora, fauna, physiognomy, and ecological conditions to large forest masses covering northern parts of North America and Eurasia, and is also known as Taiga. This type of forest usually develops in an interval of altitude between 2400 and 3600 m, and its humidity requirement is high, registering precipitations superior to 1000 mm annually. They are dense forests, of heights that can reach up to 50 m, with canopies of a triangular contour; the high density conditions reduce the amount of light reaching the interior, which limits the development of shrubs and herbaceous species [32].

#### *2.3. Forest Zoning*

This process was performed to identify the areas where the species of interest is dominant. The stands were also identified where the species is mixed with other genera; however, these were excluded from the study. This concept is also known as stratification and consists of the division of the forest area into portions or spatial units called stands, with sections that possess similar physical and biological characteristics [33,34].

The identification of the stands was made through a RapidEye-4 satellite (provided by the Space Agency of the Mexican Government) image analysis dated 25 February 2015, with a spatial resolution of5m[35]. Initially, the image was corrected atmospherically and radiometrically [36]. Subsequently, it was estimated at the pixel level "Red Edge Normalized Difference Vegetation Index" (*RedEdgeNDVI*) [37–41].

$$RedEdge\_{NDVI} = \frac{R710 - R705}{R710 + R705} \tag{1}$$

where R710 is the band 4 RedEdge or near red and R705 corresponds to the band 3 RedEdge. The NDVI has been used in the elaboration of the process of forest logging, and it is mentioned that it has shown acceptable results in the identification of vegetal associations.

The RedEdge NDVI results were analyzed along with altitude, latitude, and exposure values to determine the final logging through an overlay positioning process.

The obtained image analysis was validated with the records obtained through field trips. The existing vegetation types and the boundaries between them were determined by georeferenced points.

#### *2.4. Field Information*

#### Sampling Design and Characteristics of Sites

A random sampling design was used to define the location of 33 circular sampling sites of 1000 m2, in which a total of 682 trees were measured. The shape and size of the sites were suitable for the purpose pursued since they have shown good results for the calculation of volumetric stocks or biomass content [42,43].

At each sampling site, the diameter information at the breast height (dap) of all trees with a diameter equal to or greater than 7.5 cm and the total height of each tree [44] was recorded. For the determination of the number of annual growth rings in 2.5 cm, known as passage time, an average tree per diameter category was selected at each sampling site by means of a 250-mm Haglöf® Presser (Haglöf, Langsele, Sweden) drill and the results were calculated as the annual volumetric increase per hectare [45], which was subsequently used for estimating carbon sequestration.

#### *2.5. Information Processing*

#### 2.5.1. Calculation of Volumetric Stocks

With the records of each stand, we proceeded to estimate the value of the variables basal area and timber volume [36,46]. The basal area for each tree was obtained by the following equation:

$$B\_A = \frac{\pi}{4} \ast D\_{\text{bh}}^2 \tag{2}$$

where *BA* is the basal area (m2) and *Dbh* corresponds to the diameter at breast height (m), which were grouped in diametric categories of 5 cm.

The volume per hectare per stand was determined through a procedure known in the forest as the average hectare [47] by employing the following equation:

$$\nabla s = \sum\_{i=1}^{DC} \left( B\_{Ai} \ast \overline{H} \ast M\_{Ci} \ast E\_{Fi} \ast N\_{Ti} \right) \tag{3}$$

where *Vs* corresponds to the total volume (trunk + branches + leaves + roots) of the woodland per hectare (m3); *DC* is the number of diametric categories in the stand; *i* is the diameter category; *BA* is the basal area of the average tree by diameter category (m2); *H* corresponds to the average height of the trees in the stands for each diameter category (meters); *MC* is the morphic coefficient for gender Abies; *EF* is the expansion factor of stem to total volume (branches, leaves, and roots) [48]; and *NTi* is the number of trees per hectare per diameter category.

#### 2.5.2. Volumetric Increment

The calculation of the increment of the forest mass for each of the stands was obtained by means of the "Klepac Fast" method, which relates the number of trees per diameter category per hectare, the volume of the tree type for each diameter category, and the step time for each case [45,49]. This method considers a series of calculations between these three variables to obtain the percentage of the increase by the following equation:

$$P = \frac{1000}{Dbh} \ast \frac{1}{T} \tag{4}$$

where *P* corresponds to the percentage of increase, *Dbh* is the diameter at breast height, and *T* is the step time. It is worth mentioning that this parameter was used to obtain the carbon sequestration rate.

#### 2.5.3. Aerial Biomass Content

This information was obtained by the equation developed by Avendaño et al. [50], for the calculation of biomass for *Abies religiosa* based on the diameter at breast height.

$$B = \begin{array}{c} 0.0713 \ D\_{b\text{h}}^{2.5104} \end{array} \tag{5}$$

where *B* is the total biomass of the tree (Mg) and *Dbh* is the diameter at breast height (m).

#### 2.5.4. Carbon Content

The estimated carbon content for each tree was calculated using the equation developed for this species, which is based on the diameter at breast height [50].

$$A\_{\rm CC} = 0.0332 \, D\_{bh}^{2.5104} \tag{6}$$

where *ACC* is the carbon content for *Abies religiosa* (Mg) and *Dbh* is the diameter at breast height (m).

The carbon content for each diameter category was obtained by the following equation:

$$D\mathbf{c}\_{\mathfrak{C}} = \mathcal{A}\_{\mathfrak{C}\mathfrak{C}} \ast \mathcal{N}\_{\mathbf{T}} \tag{7}$$

where *Dccc* is the carbon content by diameter category (Mg), *Acc* is the carbon content per tree (Mg), and *NT* is the number of trees of the corresponding diametric category.

The carbon content per hectare was obtained by adding the carbon content from all the diametric categories.

2.5.5. Carbon Sequestration Rate

After calculating the carbon content per hectare and after having determined the volumetric increments, the carbon capture rate for each stand was determined, considering that the increment parameter refers to the volume increase per unit time. Once the amount of biomass that this forest can generate in a certain time period was obtained, the carbon stored was calculated using the following equation:

$$
\overline{\mathbb{C}}\_{SR} = \left( \overline{\mathbb{C}}\_{\mathbb{C}} \, / \overline{\mathbb{V}}\, \right) \* \mathbb{C}\_{AI} \tag{8}
$$

where *CSR* is the rate of carbon sequestration per hectare per year (Mg), *CC* is the carbon content (MgC ha−1), *Vs* corresponds to the volumetric stocks (m3 ha−1), and *CAI* is the current annual increase (m3).

#### **3. Results**

#### *3.1. Forest Zoning and Field Information*

Of the total area studied, eight stands were identified, with *Abies religiosa* covering an area of 1229.65 hectares. The rest of the area (603.35 hectares) corresponds to rock formations and other types of vegetation such as *Quercus*, *Juniperus*, and grassland. The study was addressed to *Abies religiosa* since it is the conservation species of interest in this protected area (Figure 2).

**Figure 2.** Forest zoning obtained through the analysis of satellite images in combination with altitude information, orientation, and slope.

The Table 1 show the values of the NDVI together with the land orientation, as well as the slope range that exists in each one of the stands, for which the criteria for the realization of the forest zoning were the values of NDVI, orientation, and species composition.



\* Refers to other uses or plant formations, "-" refers to not obtained.

#### *3.2. Volumetric Stocks*

The total volume is 757,681.69 m3, and the average volumetric stock per hectare is 616.18 m3. Table 2 details the *Abies religiosa* volumes for each stand.


**Table 2.** Volume in m3 per stand.

Stands 1, 2, 3, 4, and 7 have a density of more than 200 trees per hectare, a basal area greater than 30 m2, and a volume that exceeds 670 m<sup>3</sup> ha−1. Stands 5, 6, and 8 have densities of less than 150 trees per hectare and volumes below 470 m3 ha<sup>−</sup>1, which is due to the presence of disturbances such as forest fires and pests that have affected the density per unit area.

#### *3.3. Volumetric Increment*

Table 3 presents the results of the calculation of volumetric increase, for which we considered the records of 165 trees. It can be observed that stands 4, 5, and 6 present smaller times of passage in relation to the rest of the stands, probably due to the age and density of the forest mass. The highest productivity is found in stands 1, 4, and 5, where the area of these stands also influences the result.


**Table 3.** Result of calculating volumetric increments.

#### *3.4. Carbon Content*

The estimate of the total carbon content is 130,004.02 Mg, with an average per hectare of 105.72 MgC. The results for each stand are shown in Table 4.


**Table 4.** Carbon content per hectare and per stand in megagrams.

#### *3.5. Carbon Sequestration*

The result for the carbon capture rate in the *Abies religiosa* aerial biomass is 1267.66 MgC yr−1, considering the studied area, with an average of 1.03 MgC yr−<sup>1</sup> per hectare.

Stands 1 and 4 have higher volumetric increments. Consequently, the sequestered carbon is higher in relation to the rest of the stands (Figure 3), and this can be attributed to the north orientation of the surface where the insolation is less and more moisture is available. Stand 6 is located in a transition zone with the presence of other tree genres that were not counted but that compete with the *Abies*, decreasing its density, due to the fact that the volume and increase of biomass is reduced. Stand 8 presents a situation similar to the previous stand, only in this case the mixture of genera is due to anthropic activities (reforestation).

**Figure 3.** Annual carbon sequestration per hectare in megagrams.

Carbon sequestration in stand 1 is the one with the highest catch mainly due to the surface area and annual increment, followed by stand 4. The strata with the least carbon sequestration are the ones with the lowest surface area and the lowest density of trees (Figure 4).

**Figure 4.** Annual carbon capture above-ground per stand.

#### **4. Discussion**

#### *4.1. Carbon Content*

The area of study corresponds to a protected area decreed in 1898, in which the conservation policies that prevent the modification of the forest structure have caused the longevity of this mass, placing it in the last stage of development (old fustal), and this is derived from the age and diameters present. Diameters greater than 1.10 m were recorded, while the smallest diameter registered was 7.5 cm. On average, the number of trees per hectare is 193, with a volume of 616.18 m3. Table 5 presents comparative information of estimated carbon in protected areas with the presence of *Abies religiosa*.

**Table 5.** Comparison of studies in protected forests where they have calculated the carbon content in aerial biomass in *Abies religiosa*.


There are several hypotheses about the carbon content in forest ecosystems. This storage is dependent on the amount of existing biomass, which depends on the age, diameter, and height of the trees. When woodland density is affected or altered for some cause such as: pests, diseases, forest fires, clandestination, or harvesting, the carbon content per unit area is also altered [56]. The silviculously managed forests have young stands as a result of the application of silvicultural treatments [57], where the stands with a higher density and basal area are those containing more volume of above-ground biomass, and consequently, more carbon content.

The records consulted about carbon contents in protected temperate forest are ≥115 MgC ha−<sup>1</sup> [51–55], whereas this study estimated more than 105.72 MgC ha<sup>−</sup>1. These differences can be attributed to several factors such as the quality of the site studied, and the density and age of the trees, among others. Unlike the preserved forests, disturbed forests contain less carbon, due to the affected forest mass. Aguirre et al. [57] mentions that for a managed forest of *Pinus patula* with an orientation that is contemporary, the content is 63.98 MgC ha<sup>−</sup>1.

In the case of protected forests or silvically managed forests, the ecosystem service of carbon storage is fulfilled, whereas in the case of harvested forests, the carbon content is lower and mainly depends on applied silvicultural treatments. For the production of coetaneous masses, the content is variable and dependent on the age, diameter, height, and density of the trees in a given stand. The present study considers there to be 48% more carbon content in El Chico than the data provided by Aguirre et al. [57]. Masera et al. [58] mentions that managed forests with a temperate climate contain 118 MgC ha−<sup>1</sup> or 10% more than the estimated data in this study. The differences between Aguirre and Masera information can be mainly attributed to the applied silvicultural system (intensive or conservative).

#### *4.2. Carbon Sequestration*

Forest ecosystems sequester carbon and are considered as an option for the mitigation of the effects of an increasing atmospheric CO2 load [59–62]. The potential for carbon of any forest species depends on the maximum amount of biomass it can produce per unit of time. In species of accelerated growth, this parameter is relatively fast reached, whereas in species with a slow growth, the period of time required to reach the maximum biomass content is longer, and consequently, carbon sequestration is higher [63].

The volumetric increases represent the parameter which is able to determine the rate of carbon sequestration, and for the case presented in this study, the *Abies religiosa* forest captures 1.03 MgC ha−<sup>1</sup> yr−1, which is equivalent to 3.78 MgCO2 ha−<sup>1</sup> [64], with an annual average increase of 2.92%. Because it is located within a protected area, the forest mass has not been altered, which has caused the increment curve to decline. As the age of the mass increases, the increments decrease [45], and as a consequence, the potential for carbon capture is also affected.

Compared with protected forests, sustainably harvested forest areas capture more carbon [65], which is mainly due to the management of the age factor within the masses. Similarly, the use of harvested biomass for the production of long-lasting products retains carbon for long periods of time [66,67]. There are records of carbon capture in managed forests that exceed the results obtained in this study. Liu et al. [68] estimated net biomass productivity for forests in the Appalachian region where the forest harvest exists, and reported data ranging from 1.8 to 6.2 MgC ha<sup>−</sup>1. Zhang et al. [69] estimated a value of 2.4 MgC ha−<sup>1</sup> captured for a Massonian *Pinus* forest in China. For the specific case of *Abies religiosa* with the information of Manzanilla et al. [70], we can estimate 3.1 MgC ha−1. The differences between these investigations and the estimated rate in this study can be attributed to factors such as: forest management, location of the studied area, and species, among others. Navarro et al. [65] mentioned that these results cannot be cofrontable since they correspond to different ecosystems, each with their own particularities.

Unlike carbon content, the rate of sequestration is influenced by forest management techniques; in protected forests, the amount of stored carbon is higher, but the catch rate is reduced. With sustainable forest management practices, this difference can be balanced.

Forest management should be included in protected forests, and it is desirable to consider aspects that relate the conservation of species and their habitats with the carbon storage outside forest areas. The extraction of biomass for the elaboration of long-lasting products such as furniture or infrastructure is a way to reduce the risk of leakage within this type of ecosystem, reducing the carbon content inside a warehouse, so that emissions risks are reduced and capture is encouraged by stimulating increases in forest mass in order reach the maximum biomass production potential.

#### **5. Conclusions**

The methodology used for the evaluation of a protected forest of *Abies religiosa* as aerial carbon storage and its capture capacity is adequate and reproducible in areas with similar conditions where it is not possible to use destructive methods. The use of high resolution satellite images combined with the analysis of physical aspects of the terrain allowed for detailed zoning directly related to the density of the forest mass, and its amount of biomass and carbon.

The results presented denote that the amount of carbon stored above-ground is directly related to the density and degree of disturbance. In this case, because it is a protected forest where the silvicultural activities are restricted, the productivity of the mass is lower, and consequently, the rate of carbon sequestration decreases.

The quantified carbon additionality is a parameter that depends on several factors, for which the age of the forest mass is considered one of the most important, and in the young and vigorous forests, the rate of sequestration is higher than in forests of advanced stages of development. In this study, *Abies religiosa* caught 1.03 MgC ha−<sup>1</sup> yr<sup>−</sup>1, with an annual average increase of 2.92%.

The implementation of silvicultural techniques governed by forest management with correct principles, foundations, and objectives greatly contributes to the reduction of atmospheric carbon. Within protected forests, the application of forest management techniques makes it possible to obtain sustainable forests and maximize the potential of forest soils to increase this important ecosystem service.

**Acknowledgments:** The authors thank the Mexican Government for providing the RapidEye satellite image used for stratification, as well as for the personnel who manage the El Chico National Park for the facilities granted to carry out this study.

**Author Contributions:** P.I.F.-L. and H.J.C.-B. drafted the document and elaborated, and designed the experiment; C.A.G.-R. contributed to the experimental design planning, data analysis, revision, and translation of the manuscript. E.M.O.-S. contributed to the revision of the document, and support and advice in the development of research. J.R.V.-L., R.R.-L., and R.R.-Z. were involved in document revision, data analysis, field phase support, and cartographic phase review.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2017 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/).

## *Article* **Nonlinear Variations of Net Primary Productivity and Its Relationship with Climate and Vegetation Phenology, China**

#### **Jian Yang 1,2, Xin Chang Zhang 1,\*, Zhao Hui Luo <sup>2</sup> and Xi Jun Yu <sup>2</sup>**


Received: 8 August 2017; Accepted: 19 September 2017; Published: 25 September 2017

**Abstract:** Net primary productivity (NPP) is an important component of the terrestrial carbon cycle. In this study, NPP was estimated based on two models and Moderate Resolution Imaging Spaectroradiometer (MODIS) data. The spatiotemporal patterns of NPP and the correlations with climate factors and vegetation phenology were then analyzed. Our results showed that NPP derived from MODIS performed well in China. Spatially, NPP decreased from the southeast toward the northwest. Temporally, NPP showed a nonlinear increasing trend at a national scale, but the magnitude became slow after 2004. At a regional scale, NPP in Northern China and the Tibetan Plateau showed a nonlinear increasing trend, while the NPP decreased in most areas of Southern China. The decreases in NPP were more than offset by the increases. At the biome level, all vegetation types displayed an increasing trend, except for shrub and evergreen broad forests (EBF). Moreover, a turning point year occurred for all vegetation types, except for EBF. Generally, climatic factors and Length of Season were all positively correlated with the NPP, while the relationships were much more diverse at a regional level. The direct effect of solar radiation on the NPP was larger (0.31) than precipitation (0.25) and temperature (0.07). Our results indicated that China could mitigate climate warming at a regional and/or global scale to some extent during the time period of 2001–2014.

**Keywords:** net primary production; spatiotemporal patterns; climate change; phenology; China

#### **1. Introduction**

Net primary productivity (NPP) is the net amount of carbon accumulated by plants in a given period, and has been regarded as one of the main components of the carbon cycle [1]. Due to a variety of direct and/or indirect anthropogenic activities (e.g., land clearing and conversion) and nature disturbances (e.g., fire, pests) as well as global and regional climate change, forested ecosystems have undergone substantial changes in cover and have increasingly shown declines in health over recent decades [2,3]. As a sensitive indicator of forest cover, function and health, NPP loss may affect the composition of the atmosphere, fresh water availability, biodiversity [4,5], and the ecological adjusting mechanism of energy supply and distribution [6].

Quantifying the inter-annual variability in NPP and the interactions with climate factors would help us understand the terrestrial carbon dynamics and underlying mechanisms in responses to climate change [7]. Recently, many studies have been undertaken on the spatiotemporal variation of NPP and its relationship with climate factors on global and regional scales using a range of approaches from observational [2,8,9] to a suite of remote sensing based methods [10,11]. However, the NPP results were diverse in trends and magnitudes, even in the same region among models [10,12] due to different inputs. Moderate Resolution Imaging Spaectroradiometer (MODIS) annual NPP products have been widely used [13,14]; however, uncertainties from inputs and the algorithm resulted in biases and restrained data use regionally or globally to some extent. For instance, 70–80% accuracy of MODIS land cover products (MOD12Q1) are an assimilated meteorological dataset, not observed data with coarse spatial resolution, the cloud-contaminated MODIS FPAR/LAI (MOD15A2), and weaknesses in the MOD17 algorithm [15,16]. Therefore, it is essential to compare MODIS derived NPP to other models and/or field observed NPP, especially for China, which encompasses a wide range of ecosystems and climates. In addition, the responses of different ecosystems to different magnitudes of climate change are still far from clear, especially in China [17]. Moreover, few studies [12] have investigated the effect of solar radiation as an integrated surrogate for the effects of both day length and sunlight intensity [18] on vegetation NPP, especially when temperature and precipitation are considered simultaneously. In addition, previous studies have shown that vegetation phenology, an important factor that affects plant productivity, has changed dramatically due to climate changes and anthropogenic interference [19], but only a few studies have explored the effect of phenophase variation on NPP [20]. At present, due to the diversity of the trends and magnitudes of vegetation phenophases at different scales, its effects on NPP are still unclear, especially in China, where the vegetation phenophases are diverse at both the regional and biome levels [19]. Furthermore, a linear regression method has been applied by most studies to analyze NPP trend [10,21] despite the trend always showing a non-linear trend, or has one or more turning points within the time period [22,23]. These limitations and/or gaps have impeded our understanding of the dynamic relationship and consequently researchers may have underestimated future changes in plant productivity and/or the carbon cycle throughout China.

China encompasses a variety of ecosystems and climates. The regional climate ranging from tropical to cold-temperate, and from humid in the south to extremely dry in the northwest [24]. Land cover types are diverse, including a broad range of tropical, temperate and boreal forests, grassland, cropland and desert [25]. In recent decades, China has experienced dramatic changes in climate such as remarkably strong El Niño events [26], the freezing low temperatures in early 2008 [24], and frequent occurrences of severe droughts [27]. Meanwhile, land use and land cover changes have occurred at unprecedented rates due to quick economic development, dramatic urbanization [28], and implementation of several large scale forest plantation programs [29]. These changes have resulted in large variations in China's terrestrial ecosystem productivity and have definitely adjusted the terrestrial carbon cycle in China [30]. However, whether the temporal trend of NPP continuously increased or decreased during the study period is still unclear given the possible changes of climate derivers and anthropogenic activities. Additionally, the correlations between NPP and climate derivers as well as vegetation phenophases in recent decade still remain unclear [31].

Due to climatic variability, topographic complexity, natural ecosystem diversity, and intensive human disturbance, China is becoming one of the most critical and sensitive regions in the global carbon cycle for determining the carbon budget at regional and global scales. Furthermore, it provides a good opportunity to identify the effects of climate change on NPP to forecast the potential biosphere feedback to nature in the climate system.

Therefore, two simple models were applied to estimate NPP and the results accompanied by MODIS derived NPP were all compared with the field observed NPP. Then, the results that showed less biases with the field observed NPP were selected to analyze the vegetation NPP dynamics and its relationship with both climate and vegetation phenology in China. More specifically, the authors aim was to (1) explore which model's result was more accurate, and the quantity of uncertainty of MODIS derived NPP in China; (2) understand the spatial pattern of NPP, and investigate whether the temporal trend of NPP was continuously increasing or decreasing in China for the period 2001–2014 given climate change and anthropogenic activities; and (3) estimate the effects of climatic driving factors and vegetation phenology changes on vegetation NPP.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study focused only on Mainland China. The climate in China is extremely diverse, ranging from tropical regions in the south to subarctic in the north [11]. Therefore, the whole study area was divided into three sub-regions (Figure S1) according to the climatic regionalization of China [32]: (1) Northern China, with mean annual temperature ranges from −4 ◦C to 14 ◦C and total annual precipitation ranges from 200 mm in the northwest to 1000 mm in the southeast; (2) Southern China, where mean annual temperature ranges from 14 ◦C to 22 ◦C and total annual precipitation ranges from 1000 to 2000 mm [33]; and (3) the Tibetan Plateau, which has been called the Third Pole of the World, with an average altitude close to 4000 m above sea level, a mean annual temperature ranging from −5 ◦C to 12 ◦C and precipitation ranging from >800 mm to <200 mm [34].

#### *2.2. Dataset and Data Processing*

#### 2.2.1. Annual Net Primary Productivity (NPP)

NPP derived from MODIS from 2001–2014 in China were used (MOD17A3, 1 km), which have been further used as an important data source for plant productivity monitoring and assessment [35]. The MODIS annual NPP algorithm relies on the summation of the daily estimation of Gross Primary Productivity (GPP) computed globally minus growth and maintenance respiration [14]. To estimate annual NPP, first, the daily estimates of maintenance respiration for leaves and fine roots are subtracted from the daily GPP values. These daily reduced GPP estimates are then summed for each year, and estimates for annual maintenance respiration of living wood tissue and annual total growth respiration are subsequently subtracted, resulting in annual NPP estimates. Full details of the algorithm in MODIS derived NPP can be found in Reference [14]. However, due to the deficiencies and uncertainty of MODIS derived NPP as above-mentioned, the Miami model and Thornthwaite Memorial model were applied in this study, as they are simple to operate with only a few parameters that were easy to obtain, and have a higher spatial and temporal resolution of climate data in China. Additionally, they have previously been applied in China [11,13]. The Miami model is calculated as follows (Equations (1)–(3)):

$$NPP\_{T,P} = \min\{f\_1(T), f\_2(P)\}\tag{1}$$

$$f\_1(T) = \frac{k\_1}{1 + e^{k\_2 - k\_3 \times T}}\tag{2}$$

$$\,\_1f\_1(P)\,\,=k\_4\,\times\,\left(1-e^{-k\_5\times P}\right)\tag{3}$$

where *K*1, *K*<sup>2</sup> and *K*<sup>3</sup> in Equation (2) are the temperature response parameters with values of 3000, 1.315 and 0.119, respectively and *T* is the annual average temperature (◦C); *K*<sup>4</sup> and *K*<sup>5</sup> in Equation (3) are the precipitation response parameters with values of 3000 and 0.000664, respectively, and *P* is the annual average precipitation (mm). The parameters of *Ki* (*i* = 1, 2, ... , 5) were calculated using the least squares method based on the field testing of NPP, and the relative temperature and precipitation collected at 50 locations scattered across five continents [36].

The Thornthwaite Memorial model is expressed as follows (Equations (4)–(6)):

$$NPP \;= \; 3000 \times \left(1 - e^{-0.0009695 \times |v - 20|}\right) \tag{4}$$

$$v = \frac{1.05 \times R}{\sqrt{1 + \left(1 + 1.05 \frac{R}{L}\right)^2}}\tag{5}$$

$$L = \text{3000} + 25 \times T + 25T^3 \tag{6}$$

where *v* represents the mean annual actual evapotranspiration (mm); *L* represents the mean annual evapotranspiration (mm); and *T* and *R* represent the mean annual temperature and mean annual precipitation, respectively.

Due to the different inputs and parameters in the three models, all the estimated NPP results were compared with the field observed NPP (Table S1), which was calculated from flux tower data based on the eddy covariance method, a micrometeorological method with high-frequency data collection. This method provides direct measures of net carbon fluxes between vegetated canopies and the atmosphere over short and long timescales with minimal disturbance to the underlying vegetation. The observation data were collected from ChinaFLUX sites [37] directly and indirectly from the literature. High-precision CO2 concentration measurements and surface carbon flux measures, made at eddy covariance measurement sites, can be used to improve and validate the algorithms being used by remote sensing and ecosystem models [38]. The observations almost covered the typical ecosystem types in China. Therefore, observations from carbon flux data were applied as a reference to compare the results from the three models. Sites only having one year, the corresponding yearly data for observations, and three simulated values were used directly; and for sites longer than one year, the corresponding multi-year averaged data were applied. Next, two statistical indicators, namely root mean square error (RMSE), and mean absolute error (MAE), were used to evaluate the performance of each method above-mentioned. The best estimated result was then considered in the following analysis of the NPP relationship with climatic factors and vegetation phenology. The two statistical indicators were calculated as follows (Equations (7) and (8)):

$$\text{RMSE} = \sqrt{\frac{\sum\_{i=1}^{n} (NPP\_i - NPP^\*)}{n}} \tag{7}$$

$$\text{MAE} = \frac{1}{n} \sum\_{i=1}^{n} |NPP\_i - NPP^\*| \tag{8}$$

where *NPPi* is the estimated NPP, *NPP*<sup>∗</sup> is the field observed NPP; and *n* is the sample size.

#### 2.2.2. Climate Data and Phenology Extraction

Monthly meteorological data from 2001–2014, including temperature, precipitation, and solar radiation, were acquired from the China Meteorological Data Sharing Service System (downloaded from [39]) All meteorological data used in this study were verified by China's Meteorological Information Center (located in Beijing, China) [40], thus false or missing data from some of the stations were eliminated [11]. There are 653 meteorological stations recording temperature and precipitation data (301 in Northern China, 258 in Southern China, and 94 in the Tibetan Plateau), and 99 stations recording solar radiation (50 in Northern China, 38 in Southern China, and 11 in the Tibetan Plateau; Figure S2). The Kriging method was used for the spatial interpolation of climate data across the study area [41].

Vegetation phenology metrics (Start of Season (SOS) and End of Season (EOS) for each year) were extracted using TIMESAT software (which is widely used for simulating vegetation phenology [42,43]) by applying the Savitzky–Golay (S–G) method to generate smooth time-series MODIS EVI data. We adopted an adaptation strength of 2.0, no spike filtering, seasonal parameter of 0.5, S–G window size of 2, and amplitude season start and end of 20% to calculate the phenology parameters. Length of Season (LOS) was calculated as the difference between the SOS and EOS values.

#### 2.2.3. Land Cover Data

A 1-km spatial resolution land cover product, MOD12Q1, was applied in this study to analyze NPP and its variation across vegetation types. The MOD12Q1 land cover dataset which was extensively applied to monitor land use and land cover change [44], is mainly based on the International Geosphere–Biosphere program (IGBP) classification system that obtains a classification algorithm of

the decision tree and artificial neural network [45]. The MOD12Q1 land cover classes were further reclassified into seven major land cover classes (Figure S1) in the present study: evergreen needle-leaf forest (ENF), evergreen broadleaf forest (EBF), deciduous needle-leaf forest (DNF), deciduous broadleaf forest (DBF), farmland grassland and meadow (GM), and shrubs. However, the data only reflected the land cover classifications and could not consider changes in land cover that occurred over our study period. It was estimated that vegetation changes might be relatively stable over a short time period of approximately 10 years at a regional or global scale [42].

#### *2.3. Data Processing*

#### 2.3.1. Trend Analysis and Turning Point Year Detection

The trends of NPP during 2001–2014 were calculated at both pixel level and regional level. Due to autocorrelation among the inter-annual time series data, a robust non-parametric Mann–Kendall (M–K) trend analysis [46] was applied. This method did not require the independence and normality of the time series data [47], which has been widely used in trend analysis [19]. Previous studies have reported that the M–K test statistic *Z* was approximately normally distributed when the sample size was *n* ≥ 8. A positive or a negative *Z* value indicated an increasing or a decreasing trend respectively, which were all monotonic [46]. The formulas for the M–K method are described in detail in Reference [46]. In addition, trends of the NPP were tested at a significance level of *α* = 0.05.

The Theil–Sen median slope estimator was applied to estimate the rate of change of NPP, which was more appropriate for assessing the rate of change in short or noisy time series [48]. The Theil–Sen median slope was computed as (Equation (9)):

$$\beta\_i = \text{Median}(\frac{\mathbf{x}\_j - \mathbf{x}\_k}{j - k}) \text{ for } i = 1, \dots, N \tag{9}$$

where *xj* and *xk* are the data values at times *j* and *k* (*j>k*), respectively.

To detect the timing and magnitude of NPP changes, a linear regression model and a piecewise linear regression model (Equation (10)) were used. This method has been widely used [41,49] as it can detect potential turning points (TPs) in a trend of time-series data. In addition, TP was limited to the years 2004 to 2010 to avoid obtaining a too-short segment before or after the TP [50]. The maximum number of TPs was specified as one given that too many TPs would make the result more complex and thus create more uncertainty in the understanding of NPP trends [23]. The models used were:

$$Y = \begin{cases} \begin{array}{l} \beta\_1 \times t + \beta\_0 + \varepsilon & (t \le a) \\ \beta\_1 \times t + \beta\_2 \times (t - \alpha) + \beta\_0 + \varepsilon & (t > a) \end{array} \tag{10}$$

where *Y* is the NPP; *t* is the year; *α* is the turning point of the NPP time series; *β*<sup>0</sup> is the intercept; *β*<sup>1</sup> is the magnitude of the NPP trend before the TP; (*β*<sup>1</sup> + *β*2) is the magnitude of the NPP trend after the TP; *α* is the year of the TP; and *ε* is the residual random error.

The two models were then fitted to the NPP time-series data using the least-squares method, and the Akaike Information Criterion (AIC) [51] was used to determine whether the TPs were significant as it provides a means for model selection [41,51]. The AIC values of these two models were calculated as (Equation (11)):

$$AIC = n \times \log\left(\frac{RSS}{n}\right) + 2k + \frac{2k(k+1)}{n-k-1} \tag{11}$$

where *RSS* is the residual sum of squares for the estimated model; *k* is the number of parameters; and *n* is the sample size.

Finally, ΔAIC was defined as the difference of the AIC2 of the piecewise linear regression model (model 2) and the AIC1 of linear regression model (model 1); if the ΔAIC was less than −2, then model 2 was significantly preferred [8].

#### 2.3.2. Correlations Relating Climatic Factors, LOS and NPP

Partial correlation analysis was used to explore the relationships between NPP and climatic factors (mean annual temperature, annual cumulative precipitation, and annual cumulative solar radiation from 2001–2014) as well as LOS, which excluded the effects of other variables. The significance of the correlation coefficients was tested at a significance level of 0.05.

To understand the direct and indirect effects of climatic factors, LOS on NPP variation (excluding multi-collinearity between climatic factors and LOS), a structural equation model (SEM) with standardized data was applied to examine the influences of climate and LOS change on NPP variation. The SEM was initiated by including all possible relationships. The least significant relationship was then removed stepwise until all relationships were significant and the fit of the model did not increase further [52]. SEM was conducted by the "sem" R package [53], which was then visualized by the R package "semPlot".

#### **3. Results**

#### *3.1. Accuracy Assessment of NPP Estimation*

The accuracy of three modeled NPP results (MODIS derived, the Miami model, and the Thornthwaite Memorial model) were assessed by a comparison with the in situ observations. Validation analysis indicated that the modeled NPP results were all very significantly (Figure 1; *p* < 0.01) correlated with the in situ NPP measurements. The *R*<sup>2</sup> value for the MODIS derived NPP was 0.39, which was higher than that of Miami model (*R*<sup>2</sup> = 0.27) and the Thornthwaite Memorial model (*R*<sup>2</sup> = 0.28). Furthermore, the RMSE and MAE deceased to 0.19 kg C m−<sup>2</sup> a−<sup>1</sup> and 0.14 kg C m−<sup>2</sup> a<sup>−</sup>1, respectively for the MODIS NPP data compared to 0.29 kg C m−<sup>2</sup> a−<sup>1</sup> and 0.23 kg C m−<sup>2</sup> a−<sup>1</sup> for the Miami model, as well as 0.31 kg C m−<sup>2</sup> a−<sup>1</sup> and 0.26 kg C m−<sup>2</sup> a−<sup>1</sup> for the Thornthwaite Memorial model, respectively. The results demonstrated that the NPP derived from MODIS had an improved performance of NPP simulation by increasing *R*<sup>2</sup> by 46.18% and 40.46% when compared with the Miami model and Thornthwaite Memorial models, respectively. In addition, the NPP derived from MODIS decreased the relative RMSE and MAE by 34.01% and 38.34%, respectively, when compared with the Miami model, and by 39.71% and 45.88%, respectively, when compared with the Thornthwaite Memorial model.

#### *3.2. Spatial Pattern of NPP*

The spatial pattern of mean annual NPP derived from MODIS for the period 2001–2014 is shown in Figure 2a. The spatial pattern of annual NPP was uneven, which showed gradients decreasing from the south to the north and from the east to the west. The highest value occurred in Southern China, with values generally higher than 0.6 kg C m−<sup>2</sup> in most of that area. In contrast, annual NPP was usually lower than 0.2 kg C m−<sup>2</sup> in the Tibetan Plateau. For the remaining regions, annual NPP ranged between 0.2 and 0.6 kg C m−2. The standard deviation (Figure 2b) shows a similar spatial pattern with mean annual NPP. The highest value (more than 0.05 kg C m<sup>−</sup>2) was mostly distributed in Southern China. The values in Northern China ranged from 0.02–0.05 kg C m<sup>−</sup>2, with the exception of the northern and eastern parts of Northeast China, which was more than 0.05 kg C m<sup>−</sup>2. The lowest standard deviation value was found in the Tibetan Plateau, with values mostly less than 0.02 kg C m<sup>−</sup>2.

**Figure 1.** Comparison between the in situ Net Primary Productivity (NPP) site observations and NPP derived from MODIS (**a**); Miami model (**b**); and Thornthwaite Memorial model (**c**).

**Figure 2.** Spatial distribution of mean Net Primary Productivity (NPP) (**a**) and standard deviation (**b**) in China between 2001–2014.

*Forests* **2017**, *8*, 361

At the biome level, NPP differs in terms of vegetation types (Figure 3). Generally, forest ecosystems have a higher annual NPP than grassland (*p* < 0.001). Among the forest ecosystems, evergreen forests show a higher NPP than deciduous forests (*p* < 0.001). The highest NPP was found in EBF, with an average annual value of 0.975 ± 0.362 kg C m−2, which was almost twice that of the DBF (0.542 ± 0.301 kg C m<sup>−</sup>2) and DNF (0.329 ± 0.103 kg C m−2). ENF had the second largest NPP (0.605 ± 0.277 kg C m<sup>−</sup>2), followed by Farmland (0.422 ± 0.159 kg C m−2) and Shrub (0.317 ± 0.273 kg C m<sup>−</sup>2), respectively. Grassland had the lowest annual NPP, with an average value of 0.162 ± 0.196 kg C m<sup>−</sup>2.

**Figure 3.** Mean annual Net Primary Productivity (NPP) for each vegetation type over 2001–2014. The error bar in each column indicates the standard deviation. Different lowercase letters indicate significant differences in NPP among vegetation types (one-way ANOVA and *t*-test, *p* < 0.001).

#### *3.3. Temporal Trend of NPP*

Figure 4 details the temporal trends in annual NPP across China. Over a 14-year period (2001–2014), 17.69% of total pixels displayed either significantly increased or decreased trends (Figure 4b; *p* < 0.05) in the whole study region. Pixels with a positive (i.e., increased NPP) trend in NPP accounted for 78.94% of the total pixels, and 15.18% of the total pixels exhibited a significantly positive trend (Figure 4b,c; *p* < 0.05). The increasing rate of NPP in China mainly occurred between 0 and 0.4 × <sup>10</sup>−<sup>2</sup> kg C m−<sup>2</sup> <sup>a</sup>−<sup>1</sup> (Figure 4c). At a regional scale, 79.15% of the total pixels showed a positive trend in Northern China, and the increasing rate was mostly distributed in 0–0.6 × <sup>10</sup>−<sup>2</sup> kg C m−<sup>2</sup> <sup>a</sup>−<sup>1</sup> (Figure 4d). Similarly, 88.12% of the total pixels displayed an increasing trend in the Tibetan Plateau, where an increasing magnitude mostly occurred at 0–0.4 × <sup>10</sup>−<sup>2</sup> kg C m−<sup>2</sup> <sup>a</sup>−<sup>1</sup> (Figure 4f). In contrast, negative trends occurred in 57.64% of the total pixels in Southern China, with an absolute magnitude mainly ranging from 0.2 × <sup>10</sup>−<sup>2</sup> kg C m−<sup>2</sup> <sup>a</sup>−<sup>1</sup> to 0.8 × <sup>10</sup>−<sup>2</sup> kg C m−<sup>2</sup> <sup>a</sup>−<sup>1</sup> (Figure 4e).

**Figure 4.** Trends in annual Net Primary Productivity (NPP) within China between 2001 and 2014 (**a**); and pixels with significant (*p* < 0.05) and very significant (*p* < 0.01) trends are shown in (**b**). The count distributions of NPP trends in China (**c**); Northern China (**d**); Southern China (**e**) and the Tibetan Plateau (**f**) are also shown. A positive trend indicated an increase in NPP and vice versa; DS and DVS represent decreased significantly and very significantly, respectively; and IS and IVS represent increased significantly and very significantly, respectively. The red lines in figures (**c**), (**d**), (**e**) and (**f**) represent boundary for negative slop and positive slop.

#### *3.4. Turning Point Year of NPP*

#### 3.4.1. Turning Point Year of NPP at National and Regional Scale

To clarify the characteristics of the NPP trends, the averaged trends of the NPP and the turning point year at national and regional levels were calculated over the study period (Figure 5). At the pixel scale, 87.26% of the total pixels showed abrupt changes in inter-annual NPP variation (ΔAIC < −2; Figure S3b), most of which appeared at year 2004 or year 2010 (which may also have been caused by criterion we set up in the calculation), accounting for 29.36% and 31.78%, respectively (Figure S3a). At a national scale, the linear regression model suggested an increasing trend of 0.0007 kg C m−<sup>2</sup> a−<sup>1</sup> (*p* = 0.3118) from 2001 to 2014 (Figure 5a). However, two distinct periods were clearly identified by the piecewise linear regression model: a trend change from a marginal significant increasing trend of 0.0104 kg C m−<sup>2</sup> a−<sup>1</sup> (*p* = 0.0746) before 2004 to a non-significant increasing trend of 0.0012 kg C m−<sup>2</sup> a−<sup>1</sup> (*p* = 0.2775) after 2004 (Figure 5a). The information criterion of the piecewise regression model was less than that of the linear regression model (ΔAIC = −3.1797). At a regional scale, a marginal significant increasing trend of 0.0025 kg C m−<sup>2</sup> a-1 (*p* = 0.0509) in Northern China was found, and the turning point year occurred in 2010 (ΔAIC < −2). However, the increasing trends were all non-significant before and after 2010, with a rate of 0.0006 and 0.0068 kg C m−<sup>2</sup> a<sup>−</sup>1, respectively (Figure 5b). In Southern China, a rate of −0.0007 kg C m−<sup>2</sup> <sup>a</sup>−<sup>1</sup> (*<sup>p</sup>* = 0.6559) was found for annual NPP variation, and a turning point year occurred in 2010 (ΔAIC < −2). However, the variation trend was opposite (Figure 5c). Although a significant increasing trend of 0.0016 kg C m−<sup>2</sup> a−<sup>1</sup> (*p* < 0.01) was found in the Tibetan Plateau over the past 14 years (Figure 5d), the increasing trend slowed down after 2006 (ΔAIC < −2).

**Figure 5.** Inter-annual variations of Net Primary Productivity (NPP) in China at national level (**a**) and regional level (Northern China (**b**); Southern China (**c**); and the Tibetan Plateau (**d**)). Trends estimated by the least-squares linear regression are shown. The blue and pink lines indicate the linear fit before and after the turning year, respectively. The red line indicates linear fit during the period 2001–2014. TP represents the turning point year.

#### 3.4.2. Turning Point Year of NPP at Biome Scale

At the biome level, the variation of NPP was diverse (Figure 6). All vegetation types had two clearly distinct periods over the study period (ΔAIC < −2), with the exception of the EBF (ΔAIC = −1.9314) which only showed a non-significant decreasing trend of −0.0001 kg C m−<sup>2</sup> <sup>a</sup>−<sup>1</sup> (*p* = 0.9314) during 2001–2014. Farmland and grassland all increased, with a rate of 0.0011 kg C m−<sup>2</sup> a−<sup>1</sup> (*p* = 0.3734) and 0.0016 kg C m−<sup>2</sup> a−<sup>1</sup> (*p* < 0.01), respectively. In addition, the increasing trend before 2004 was larger than that after 2004 for both farmland and grassland. An increasing trend of NPP was found for DNF, ENF, and DBF, but the opposite trend was found for shrub, despite a turning point year occurring in 2010 (all ΔAIC < −2 and may also have been caused by the criterion we set up in the calculation) for all above vegetation types. In addition, a decreasing NPP occurred for ENF, DBF and shrub during 2001–2010, while an increasing NPP occurred for DNF during the same time period. After 2010, the trend of these four vegetation types was all positive, despite the existence of different magnitudes.

**Figure 6.** Inter-annual variations of Net Primary Productivity (NPP) in China at the biome level (Farmland (**a**), Grassland (**b**), Deciduous needle-leaf forest (**c**), Evergreen needle-leaf forest (**d**), Deciduous broadleaf forest (**e**), Shrub (**f**), and Evergreen broadleaf forest (**g**)). Trends estimated by the least-squares linear regression are shown. The blue and pink lines indicate the linear fit before and after the turning year, respectively. The red line indicates linear fit during the period 2001–2014. TP represents the turning point year.

#### *3.5. Relationship Relating Climatic Factors, LOS and NPP*

#### 3.5.1. Correlation Analysis

Climatic factors were generally positively correlated with NPP on a national scale (Figure 7). Positive partial correlations were observed for 62.66%, 70.79% and 57.60% of the total pixels with temperature, precipitation, and solar radiation, respectively. The positive correlation coefficients for the climatic variables ranged from 0.1–0.7. About 56.23% of the total pixels showed positive partial correlations between NPP and LOS, and the positive partial correlation coefficient was mostly distributed between 0.1 and 0.5.

At a regional level, 62.97% of the total pixels in Northern China displayed a negative relationship between NPP and temperature, while precipitation, solar radiation and LOS all generally showed a positive relationship with the NPP (Table 1). Over 57% of the total pixels in Southern China indicated a positive relationship between the NPP and climatic factors as well as LOS. However, the relationship in the Tibetan Plateau was diverse. Specifically, over 83% of the total pixels for temperature and about 58% of the total pixels for precipitation displayed a positive relationship with the NPP, while in general, an opposite relationship occurred between the NPP and solar radiation. In addition, the relationship between the NPP and LOS in the Tibetan Plateau was ambiguous.

**Table 1.** Partial correlation between the Net Primary Productivity (NPP) and climatic factors as well as Length of Season (LOS).


+ Represent positive relationship; − Represent negative relationship.

**Figure 7.** The partial correlation coefficients and corresponding frequency distribution between the NPP and temperature (**a**,**b**), precipitation (**c**,**d**), solar radiation (**e**,**f**) and LOS (**g**,**h**) in China from 2001 to 2014. The red lines in figures (b), (d), (f) and (h) represent boundary for negative relationship and positive relationship.

#### 3.5.2. Structural Equation Models

Our structural equation model (SEM) did not detect a relationship between temperature, precipitation and solar radiation, which indicates that there was a strong relationship among the climatic factors, LOS and NPP. In the SEM (χ<sup>2</sup> = 1.0138, d.f. = 5, *<sup>p</sup>* = 0.9614, AIC = 21.0138, BIC = −1815), climatic factors explained a total of 23% of the variation in LOS, and temperature and solar radiation had a positive effect on LOS, while precipitation had the opposite effect (Figure 8). Moreover, the effect of climatic factors and LOS on NPP was 0.63, and explained a total of 56% of the variation in NPP. Among the climatic factors, the direct effect of solar radiation on NPP was the largest (0.31), followed by precipitation (0.25), and temperature (0.07).

**Figure 8.** Structural equation model relating climatic factors, Length of Season (LOS) and Net Primary Productivity (NPP) in China. Single headed arrows indicate directional relationships, while double headed arrows indicate covariances. Numbers in brackets are *R*<sup>2</sup> values.

#### **4. Discussion**

#### *4.1. Uncertainties in NPP Estimates*

Net primary productivity (NPP), an indicator of the accumulation of atmospheric CO2 in terrestrial ecosystems, plays a crucial role in global change [54]. The accurate estimation of NPP is a crucial step for reliably quantifying carbon fluxes between the atmosphere and terrestrial ecosystems [31], especially for regions with different topography and climatic conditions. Although deficiency and biases exist in the MODIS derived NPP, the validation illustrated that the model applied by MODIS was superior to the Miami model and Thornthwaite Memorial model, with an increased R2 value and a decreased RMSE and MAE. The RMSE and MAE between the MODIS derived NPP and observations were only 0.19 kg C m−<sup>2</sup> and 0.14 kg C m<sup>−</sup>2, respectively. This indicated that the MODIS derived NPP may be suitable for analysis with relatively low regional biases in China for the period of 2001–2014 to some extent despite coarse global climate data applied in MOD17. Additionally, it illustrated that the parameters input and algorithms may affect the results. For instance, the Miami model is one of the first global empirical models that only utilizes temperature or precipitation in the model [55], while the Thornthwaite Memorial model determines NPP for a particular location as the actual evapotranspiration functions [36]. The interactions among climatic factors and other influential factors such as vegetation types are not included in these two models [11]. Furthermore, there is no mechanism to account for changing vegetation density in the two models [56]. However, the model applied by MODIS incorporates biogeochemical principles in a mechanistic modeling environment and the vegetation feedback to climate conditions through changes in Leaf Area Index and absorbed radiation, and has the advantage of providing spatially continuous estimates with a consistent methodology, which is important for any large-scale studies [35]. It should be noted that NPP is not easy to measure and only 26 field sites were included in this study, which may cause biases due to the larger study area. In addition, heterogeneity, stand density, and pixel resolution may also

affect validation results. Therefore, increased observation sites and an extended observation period has become essential for NPP estimation. Additionally, forest inventory data is the only data source that provides large-scale consistent productivity assessments [57]. Therefore, the application of forest inventory data rather than numerous field observed data in this study would be more preferable to validate the suitability of MODIS NPP in China, and more efforts should be made in the future. What is more, temporal analysis of match between models and observations based on a long time period were also a good choice for NPP estimates and comparison.

#### *4.2. Spatiotemporal Variation of NPP*

In general, all three NPP estimates showed similar spatial patterns, with mean annual NPP higher in the southeast and gradually decreasing towards the northwest (Figure 3 and Figure S2). This finding was consistent with previous studies [10–12,58]. This may have been due to the climatic gradient (e.g., temperature and precipitation) from the southeast to the northwest in China, which is more favorable for vegetation growth in southern China [11,59]. The lowest NPP value in northwest China and the Tibetan Plateau, which was consistent with References [10,11], may have resulted from low temperature and the absence of precipitation. In addition, vegetation type was also a key factor that affected NPP spatial distribution. In Southern China, the dominant vegetation type are evergreen forests (Figure S1) which have a higher productivity than that of the grassland and desert vegetation (Figure 4) distributed mostly in Northern China and the Tibetan Plateau. Productivity in the southeast Tibetan Plateau was higher than that of other areas in the Tibetan Plateau as demonstrated by our conclusion above.

The total amount of NPP increased at a rate of 0.009 Pg C a−<sup>1</sup> in China from 2001–2014, which was similar with Reference [10], who found a rate of 0.008 Pg C a−<sup>1</sup> in China from 1999–2010. The decreases in NPP over Southern China were more than offset by increases in NPP over Northern China and the Tibetan Plateau. The overall increased NPP could help to mitigate climate warming at regional scale and/or global scale to some extent. In terms of trend variation, a continually increasing trend of NPP in China from 2001–2014 was found, despite the occurrence of a turning point year which was similar with the previous finding in Reference [10]. This may be due to the fact that over the past three decades, the central government of China has decided to combat severe environmental degradation, including declining vegetation cover and expanding desertification [29]. To realize this goal, the Three North Shelter Forest System project, the Beijing-Tianjin Sand Source Control Program, and the Grain for Green Project [29] were implemented. Alternatively, as the turning point is related to ENSO events/cycles, it is necessary to validate this conclusion in the future.

At a regional scale, NPP trends in China showed a prominent geographical heterogeneity. The trends before and after the turning point year in Northern China and the Tibetan Plateau were all positive. In contrast, the trend in Southern China was the opposite, which may have been due to ecorestoration, forest/grassland protection, and reforestation in Northern China and the Tibetan Plateau after 2000 [29]. However, in Southern China, freezing low temperatures in early 2008, the severe drought in 2009 [31], and decreased solar radiation [12] may have all contributed to the decreasing NPP. Due to the physiological and/or local environmental conditions, the trend and turning point year among vegetation types were diverse.

#### *4.3. Diverse NPP Correlations with Climate Drivers and LOS*

Vegetation NPP is influenced by a variety of factors, and the primary one is climate-related. Temperature, precipitation, and solar radiation are the main climate drivers for vegetation growth [59], but the relationships vary with spatial scale. Nevertheless, limited efforts have been made to investigate the relative roles of different climate variables in the NPP [10]. Generally, temperature has been found to be positively correlated with the NPP at a national level, which is consistent with previous studies [10,60]. This response may in part be related to the increased activity of photosynthetic enzymes [61] and improved capacity for photosynthesis and growth [62]. However, the relationship

between temperature and NPP in Northern China was negative. This may have resulted from higher temperatures leading to increased water scarcity caused by accelerated evaporation, which would work against vegetation growth, especially for arid and semi-arid areas in Northern China [49]. Additionally, the positive relationship between precipitation and NPP in Northern China illustrated this conclusion. Precipitation showed a positive correlation with NPP at both the national scale and regional scale, which may have occurred as more precipitation increases soil moisture, satisfying the water requirements for vegetation growth and productivity increase, especially in late spring and summer, which usually have higher temperatures and more evapotranspiration [63].

Little is known about the light-related physiological mechanisms that regulate NPP due to the difficulty in distinguishing the effects of day length (i.e., photoperiod) and light intensity [18], especially when studies are extended to consider their correlations with temperature and precipitation. Therefore, in this study, solar radiation was considered as an integrated surrogate for both day length and sunlight intensity [18] to investigate the effects of solar radiation on the NPP. A positive correlation was observed between the NPP and solar radiation at both the national and regional scales, except for a negative relationship which occurred in the Tibetan Plateau. Generally, increased solar radiation provides sufficient materials and energy for vegetation photosynthesis and solar radiation is typically accompanied by warmer temperatures and sufficient sunlight intensity, both of which can enhance photosynthetic capacity and promote vegetation NPP. The negative correlations between solar radiation and NPP in the Tibetan Plateau may be partly related to larger areas of melting seasonal snow and permafrost soils caused by abundant solar radiation, leading to higher soil moisture content and creating an anaerobic soil environment within the plant root zone, thus limiting vegetation growth [64]. Alternatively, increased solar radiation may enhance surface soil evaporation and limit water availability for herbaceous plants that have shallow root systems.

Vegetation phenology has long been regarded as an important factor that affects vegetation productivity. Consistent with previous studies in References [20,65], our results also indicated a positive relationship between the NPP and LOS at both national and regional levels. This suggests that an extension of LOS is one of the most important factors that affect plant productivity [66].

#### *4.4. Shortcomings and Uncertainties*

In this study, we used data from only 99 solar radiation stations, as meteorological stations in the Tibetan Plateau are scarce. The sparse distribution of climate data limited the detail and accuracy of the relationship between climatic factors and plant productivity. Therefore, the results are likely to be subject to some ambiguity.

Vegetation NPP is influenced by a variety of factors. In this study, only three climatic factors and LOS were considered. Although these factors have explained a total of 56% of the variation of NPP (Figure 8), other climate-related factors, such as sunshine duration, soil temperature, precipitation characteristics (e.g., effective precipitation, precipitation intensity), CO2 concentrations, N enrichment and deposition (e.g., policy and planning changes), should be considered in future studies [67]. Moreover, other phonological metrics (e.g., Start of Season, End of Season), species competition, and disturbances such as anthropogenic activities (irrigation, fertilization [68,69], harvest, land use/land cover change), wildfires, plant diseases, pests, floods, and droughts as well as the time-lag effect of the above–mentioned variables can vary by region [50], and should also be considered in future studies.

The results presented here indicate the complexity of vegetation NPP dynamics and the response or reaction strategy of vegetation to climate change during the study period. However, the trend for NPP before and after turning points were mostly non-significant despite all turning points being significant (except evergreen broadleaf forest) based on ΔAIC. This may be caused by a short time period, or a non-linear variation on NPP during the study period. Therefore, the conclusions mentioned above should be made with more caution, and more effort should be made to explore the NPP trends (e.g., non-linear trend) in the future to clarify the direction of dynamic NPP trends [70]. In addition, it is difficult to estimate future NPP dynamics using mathematical models (e.g., linear regression analysis) due to the nonlinear characteristics of NPP dynamics in a long time series, the uncertainty about future climate change, and the time-lag effect of NPP responses to climate change. To date, the Hurst exponent, which has been widely used in hydrology, climatology, economics, geology, and geochemistry, could be considered to predict the NPP variation. The results presented here (Figure S4) indicate that NPP in China will continue to increase based on the Hurst exponent. Unfortunately, this exponent characterizes trends based only on past and present environmental conditions, without considering future environmental change, especially in developing areas, or giving a time period for dynamic NPP trends. These limitations restrict our ability to predict and anticipate future NPP variation, the carbon cycle and its relationship with and response to climate change and human activity, and hence more effort must be made to predict future NPP variation.

#### **5. Conclusions**

In this study, the spatiotemporal patterns of NPP and its correlation with climatic factors as well as vegetation phenology during 2001–2014 in China were investigated, and the main conclusions can be summarized as follows:


**Supplementary Materials:** The following are available online at www.mdpi.com/1999-4907/8/10/361/s1, Table S1: Sites information of the field Net primary productivity (NPP) data used in this study; Figure S1: Vegetation types (DNF, deciduous needle-leaf forest; ENF, evergreen needle-leaf forest; EBF, evergreen broadleaf forest; DBF, deciduous broadleaf forest; and GM, grassland and meadow) as well as three sub-region divisions (**a**), and provinces distribution (**b**) in China; Figure S2: Distribution of meteorological stations across China, red points represent stations for temperature and precipitation, and black triangles represent stations for solar radiation; Figure S3: Spatial distribution of Turing Point (TP) years (**a**) and corresponding ΔAIC values (**b**) of the NPP in China for the period 2001–2014; Figure S4: Spatial distribution of the Hurst exponent in China from 2001 to 2014.

**Acknowledgments:** This study is supported by the National Natural Science Foundation of China (Grant No. 41431178), the Natural Science Foundation of Guangdong Province, China (Grant No. 2016A030311016), the Innovation Project of Guangdong Province Water Resources Department (Grant No. 2015-02) and the Central Fund Supporting Nonprofit Scientific Institutes for Basic Research and Development (PM-zx021-201407-007).

**Author Contributions:** Xinchang Zhang and Xijun Yu outlined the research topic, assisted with manuscript writing, and coordinated the revision activities. Jian Yang and Zhaohui Luo performed data collection, data analysis, the interpretation of results, manuscript writing, and coordinated the revision activities.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Haberl, H.; Erb, K.H.; Krausmann, F.; Gaube, V.; Bondeau, A.; Plutzar, C.; Gingrich, S.; Lucht, W.; Fischer-Kowalski, M. Quantifying and mapping the human appropriation of net primary production in earth's terrestrial ecosystems. *Proc. Natl. Acad. Sci. USA* **2007**, *104*, 12942–12947. [CrossRef] [PubMed]


© 2017 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/).

*Article*
