Spatial pattern and driving factors of biomass carbon density for natural and planted coniferous forests in mountainous terrain, eastern Loess Plateau of China

Abstract

Background

Understanding the spatial pattern and driving factors of forest carbon density in mountainous terrain is of great importance for monitoring forest carbon in support of sustainable forest management for mitigating climate change.

Methods

We collected the forest inventory data in 2015 in Shanxi Province, eastern Loess Plateau of China, to explore the spatial pattern and driving factors of biomass carbon density (BCD) for natural and planted coniferous forests using Anselin Local Moran’s I, Local Getis-Ord G* and semivariogram analyses, and multi-group structural equation modeling, respectively.

Results

The result of spatial pattern of BCDs for natural forests showed that the BCD was generally higher in the north but lower in the south of Shanxi. The spatial pattern for planted forests was substantially different from that for natural forests. The results of multi-group SEM suggested that elevation (or temperature as the alternative factor of elevation) and stand age were important driving factors of BCD for these two forest types. Compared with other factors, the effects of latitude and elevation on BCD showed much greater difference between these two forest types. The difference in indirect effect of latitude (mainly through affecting elevation and stand age) between natural and planted forests was to some extent a reflection of the difference between the spatial patterns of BCDs for natural and planted forests in Shanxi.

Conclusions

The natural coniferous forests had a higher biomass carbon density, a stronger spatial dependency of biomass carbon density relative to planted coniferous forests in Shanxi. Elevation was the most important driving factor, and the effect on biomass carbon density was stronger for natural than planted coniferous forests. Besides, latitude presented only indirect effect on it for the two forest types.

Introduction

The role of forests as terrestrial sinks of atmospheric carbon dioxide has received increasing attention since the late 1990s (FAO 2010, 2015). Assessment of carbon density and storage of forests has been a focus in recent years (Fang 2001; Pan et al. 2011). Mountain forests account for 28% of the total global forest area (Price et al. 2000), and are considered to be important terrestrial sinks of carbon dioxide (Chang et al. 2015; Zald et al. 2016). Understanding the spatial patterns and driving factors of forest carbon density in mountainous terrain is of great importance for monitoring forest carbon in support of sustainable forest management for mitigating climate change (Zald et al. 2016).

The spatial pattern information of forest carbon density can play an important role in evaluating carbon sequestration potentials and forest management (Fu et al. 2015; Lin et al. 2017; Verkerk et al. 2019). Ren et al. (2013) examined the spatial pattern of carbon density in forest ecosystems in Guangdong Province, and found that the spatial distribution of forest carbon density was uneven, and it was a reflection of the difference in forest management and economic and social development. Fu et al. (2014) observed a clear spatial pattern of forest litter carbon density in Zhejiang Province using Anselin Local Moran’s I and geostatistical interpolation. A similar result was found for forest biomass carbon density in the same province (Fu et al. 2015). Using the same methods, Dai et al. (2018) further found that the forest carbon density decreased from southwest to northeast in Zhejiang Province, roughly in line with the topographic feature across this province. Using Anselin Local Moran’s I, Local Getis-Ord G* and geostatistical interpolation, Lin et al. (2017) analyzed the spatial variability of forest carbon density in Jiangle County, Fujian Province. However, little attempt has been made to compare the spatial patterns of carbon densities between natural and planted forests (Wang et al. 2018a).

Forest carbon density and their spatial patterns in mountainous terrain may be the result of complex interaction among various driving factors. Temperature and precipitation are two critical climatic factors driving biomass carbon density. For example, structural equation modeling (SEM) revealed that forest biomass carbon density were primarily controlled by climate factors, and increased with rising temperature and precipitation along the north-south transect of eastern China (Wen and He 2016); while multi-group SEM showed that mean annual temperature had negative indirect effect on carbon density in old-growth tropical forests across the globe (Durán et al. 2015). Meanwhile, previous studies have revealed that stand age and canopy coverage are important driving factors of forest carbon density. For instance, Xu et al. (2010) reported that forest age had significant effect on biomass carbon density for most forest types in China. The average above-ground biomass carbon density of mature forests increased with stand age for the forests younger than 450 years (Liu et al. 2014); while the biomass carbon density of ponderosa pine forest increases rapidly until 150–200 years (Law et al. 2003). The SEM indicated that canopy coverage and stand age were the most important driving factors of biomass carbon density in subtropical forests (Xu et al. 2018a, 2018b). Moreover, the effect of elevation on forest carbon density has attracted great attention, especially in mountainous areas. For example, Gairola et al. (2011) found a significant positive relationship of forest biomass carbon density with elevation in moist temperate valley slopes of the Garhwal Himalaya (India). Liu and Nan (2018) reported a positive relationship between carbon density and elevation for three natural coniferous forests in the Guandi Mountain, Shanxi Province. In addition, as another important geographical factor, latitude may have indirect effect on forest biomass carbon density by affecting plant growth and productivity through temperature and precipitation effects. It was found that the forest biomass carbon density decreased with increasing latitude along the north-south transect of eastern China (Wen and He 2016); while the biomass carbon density of Moso bamboo forests linearly increased with latitude (Xu et al. 2018a, 2018b).

Shanxi Province is a mountainous terrain in north China, characterized by six mountain ranges, namely the Hengshan and Wutai Mountains in the north, the Zhongtiao Mountains in the south, the Lüliang Mountains in the west, the Taihang Mountains in the east, and the Taiyue Mountains at the center (Fig. 1a). Coniferous forests are widely distributed in the province, prominently in the high mountains. In the past few decades, China has launched a series of large-scale afforestation and forest protection programs, including the Three-North Shelter Forest Program since 1978, the Program for Conversion of Cropland into Forests since 1999, and the Natural Forest Protection Program since 2000. With the implementation of these national projects, the areas of natural and planted coniferous forests (NCFs and PCFs, respectively), and their contributions to regional carbon stocks may have changed substantially in this region.

Fig. 1

Maps showing the general topography in Shanxi Province (a) and the distribution of the sample plots for the investigation of biomass carbon density of natural and planted coniferous forests (NCFs and PCFs) in the region (b). The number in parentheses denotes the number of the sample plots

Previous studies have examined the carbon density and carbon stocks of Pinus tabulaeformis forest in Shanxi on a local (Sun 2011; Cheng et al. 2012; Chi et al. 2014) or provincial scale (Wang et al. 2014a). Recently, Liu and Nan (2018) examined the carbon stocks of three natural coniferous forests (Larix principis-rupprechtii forest, Picea meyerii forest and Pinus tabulaeformis forest) along an altitudinal gradient from 1200 to 2700 m in the Guandi Mountain. Wang et al. (2018a) analyzed the spatial patterns of carbon densities of natural and planted forests in the Lüliang Mountains. An integrated study of carbon density of coniferous forests is urgently needed for a proper evaluation of the current capacity of carbon sequestration of natural versus planted coniferous forests and for a better understanding of how various factors affect the carbon density of natural versus planted coniferous forests in the mountainous region.

Within this context, the main objectives of this study were to quantify the biomass carbon densities and stocks of natural and planted coniferous forests in Shanxi Province, China, and to explore the spatial patterns and driving factors of biomass carbon density of natural and planted coniferous forests in the mountainous terrain. We first used the national forest inventory data in 2015 to quantify the biomass carbon densities based on the method of biomass expansion factor. We then characterized the spatial patterns of biomass carbon density of natural and planted coniferous forests using Anselin Local Moran’s I, Local Getis-Ord G* and semivariogram analyses. Finally, we used multi-group structural equation modeling to evaluate relative importance of driving factors selected in explaining variations of biomass carbon density, to estimate the direct and indirect effects of each driving factor, and to reveal the similarities and differences between the models for natural and planted coniferous forests.

Materials and methods

Study region

Shanxi Province is located in the middle reach of the Yellow River, and the east of the Loess Plateau in north China, at 34°34′–40°43′ N and 110°14′–114°33′ E, with a total area of 156,271 km² (Fig. 1a). As a mountainous plateau, the mountains and hills make up 80.3% of the province’s land area (Fan and Wang 2011). Climatically, Shanxi belongs to the temperate monsoon climate zone (Guo et al. 2015) and is influenced by the summer and winter monsoons, with a hot and wet summer and a dry cold winter. The total area of forests is about 32,100 km², accounting for 20.05% of the province’s land area. Coniferous forests account for 43% of all the forests in the province (Editorial Board of Forests in Shanxi 1992; Ma 2001).

Data collection and preprocessing

Forest inventory data

The data from a total of 508 field sample plots were used in this study. Of all the plots from the 9th national forest inventory data in 2015, 289 and 219 are natural and planted coniferous forests plots, respectively (Table 1, Fig. 1b). The permanent plots (each with an area of 667 m²) were established systematically based on a 4 km × 4 km grid (Xiao 2005). The data collected in each plot included forest origin, forest type, forest stand factors (stand age, age group, diameter at breath height of 1.3 m (DBH), tree height and forest coverage), and geographic location (latitude, longitude and elevation). For trees with DBH ≥ 5 cm, their DBH were recorded. A total of 12 dominant coniferous tree species were identified over the study region. According to the dominant species, all the natural (planted) coniferous forests were classified into five categories: Pinus, Platycladus, Larix, Picea and other forests (Fig. 1b).

Table 1 Descriptive statistics of stand age, stand coverage, elevation, latitude, annual mean temperature and annual precipitation for natural and planted coniferous forests in the study region

For each forest plot (each forest stand), the carbon density (CD) (Mg∙ha^− 1) of living trees was estimated using the biomass expansion factor (BEF) method (Wang et al. 2018a):

$$ CD=\left(V\times WD\times BEF\right)\times \left(1+R\right)\times CF $$

(1)

where V is the volume of biomass (m³∙ha^− 1), i.e., the sum of the above-ground trunk volumes of individual living trees; WD is the wood density (Mg∙m^− 3), BEF is the biomass expansion factor, R is the root-shoot ratio, and CF represents carbon content in oven-dried biomass, i.e., the ratio of carbon (Mg) to oven-dried biomass (Mg). The V for individual trees was calculated using the diameter at breath height of 1.3 m and the height of each tree based on the formula for each tree species (Wang et al. 2018b).

In this study, the biomass carbon density (BCD, Mg∙ha^− 1) represented only the carbon density of living trees, and did not include that from shrubs, herbs, litter, and soil. The mean BCD of each forest type was estimated by dividing the sum of biomass carbon densities of the forest plots by the number of the forest plots for the forest type. The total biomass carbon stock (Mg) for each forest type was estimated by multiplying the mean biomass carbon density by the area of the forest type, which was approximated by multiplying the representative area (4 km × 4 km) of a sample plot by the number of the forest plots.

Climate data

The data of annual mean temperature (TEMP) and mean annual precipitation (PRCP) for the period 1981–2010 at the 61 meteorological stations were obtained from the meteorological database of the scientific data platform of Shanxi Province. The TEMP (PRCP) for each sample plot was computed by interpolation of the measurements (1981–2010 means) of the 61 stations using regression-kriging interpolation method (Lamsal et al. 2011; Plouffe et al. 2015), which combined regression of the climatic variable on topography variables (slope degree, slope aspect and slope position) with kriging of the regression residuals (Eldeiry and Garcia 2010). Similar as in the previous study (Wang et al. 2018a), the TEMP of each plot from the interpolation was further corrected according to the difference between the plot elevation and the interpolated elevation using the temperature lapse rate of 4.89 °C∙km^− 1 as the correction factor (Wang et al. 2014b).

Data analysis

Statistical analyses

The dataset of mean BCDs for each forest type met the normal distribution and passed the homogeneity of variances test, and then One-way ANOVA with Duncan’s multiple tests was used to examine the difference in mean BCD between/among forest types.

Spatial pattern analyses

Anselin Local Moran’s I, Local Getis-Ord G* and semivariogram analyses were used to investigate the spatial patterns of forest BCDs from different aspects. Anselin Local Moran’s I can be used to identify spatial clusters and outliers (Fu et al. 2015). The significant positive values of the index show clusters that represent positive spatial autocorrelation; the significant negative values of this index show outliers that represent negative spatial autocorrelation; and the “not significant” values represent no spatial autocorrelation. The clusters consist of high-high clusters (high values in a high value neighborhood) and low-low clusters (low values in a low value neighborhood); and the outliers include high-low outlier (a high value surrounded by low values) and low-high outlier (a low value surrounded by high values) (Fu et al. 2015; Lin et al. 2017; Dai et al. 2018). The statistic of Anselin Local Moran’s I can be expressed as:

$$ {I}_i=\frac{x_i-\overline{X}}{S_i^2}\sum \limits_{j=1,j\ne i}^n{w}_{i,j}\left({x}_j-\overline{X}\right) $$

(2)

$$ {S}_i^2=\frac{\sum \limits_{j=1,j\ne i}^n{\left({x}_j-\overline{X}\right)}^2}{n-1} $$

(3)

where x_i and x_j are the values of attribute for observations i and j, respectively, $\overline{X}$ is the mean of the corresponding values of the attribute, w_{i, j} is the spatial weight between observations i and j, which can be defined as the fixed distance, which was obtained based on the largest global Moran’s I value, indicating the strongest spatial autocorrelation of carbon density (Fu et al. 2014). The w_{i, j} is given the same weight within the distance, while those outside the distance band are given the weight of 0. The symbol of n is the total number of observations.

Meantime, Local Getis-Ord G* can be used to detect spatial hot spots and cold spots that represent clusters of similar values that are significantly larger and smaller than the mean, respectively. There is a Z score for each target value in the test. For statistically significant positive Z scores, the larger the Z score is, the more intense the clustering of high values (hot spot). For statistically significant negative Z scores, the smaller the Z score is, the more intense the clustering of low values (cold spot). The “not significant” positive or negative Z scores represent no spatial autocorrelation (Manepalli et al. 2011). Getis-Ord Gi* can be described as:

$$ {G}_i^{\ast }=\frac{\sum \limits_{j=1}^n{w}_{i,j}{x}_j-\overline{X}\sum \limits_{j=1}^n{w}_{i,j}}{\sqrt{\frac{\sum \limits_{j=1}^n{x}_j^2}{n}-{\left(\overline{X}\right)}^2}\sqrt{\frac{n\sum \limits_{j=1}^n{w}_{i,j}^2-{\left(\sum \limits_{j=1}^n{w}_{i,j}\right)}^2}{n-1}}} $$

(4)

where x_i, x_j, $\overline{X}$, w_{i, j}, and n are the same as in Eq. 2.

Anselin Local Moran’s I can be used to study a certain spatial regularity, while Local Getis-Ord Gi* enables differentiation between clusters of similar values that are high or low relative to the mean, thus aiding in the detection of unusual events (Lin et al. 2017). The two indicators of exploring spatial autocorrelation were carried out using the modules of Cluster and Outlier Analysis (Anselin Local Moran’s I) and Hot Spot Analysis (Getis-Ord Gi*) of Spatial Statistics Toolbox within software ARCGIS 10.2.

Furthermore, semivariogram analysis was used to explore the spatial heterogeneity of forest biomass carbon density. Semivariogram analysis is a useful tool for analysis and interpretation of spatial data in ecology (Robertson 1987; Rossi et al. 1992) and has been used in quantification of spatial heterogeneity of forest carbon density (Wang et al. 2018a). We performed this analysis based on the standard principles of semivariogram model fitting in geostatistics (Wang 1999) using GS+ 9.0 Geostatistics Software (Gamma Design Software, LLC). This analysis was not only performed for the data of BCD but also for the data from the cluster and outlier analysis, and hot spot analysis to examine further the spatial heterogeneity of BCD of natural versus planted forests over the study region.

Multi-group SEM

Multi-group SEM was used to compare the direct and indirect effects of the main driving factors on BCD between NCFs and PCFs. Before constructing the multi-group SEM model, we tested the collinearity of the driving factors, and eliminated annual mean temperature (TEMP) from further analysis (see Annex 1 of Additional file 1).

We began with the construction of a priori model for the multi-group SEM procedure. The priori model was constructed to include all hypothesized causal relationships between BCD and potential driving factors based on the known effects of driving factors and the relationships among driving factors of forest growth and biomass accumulation (Bollen 1989; Shipley 2004; Xu et al. 2018a, 2018b). As shown in Annex 2a of Additional file 1, five driving factors were included in the priori model: stand age, stand coverage, mean annual precipitation, latitude, and elevation.

The second step was to find a good-fitting structural equation model (Grace et al. 2010) to the data of all coniferous forests based on the priori model. We compared the actual and estimated covariance matrices to determine whether the model was acceptable using the chi-square (χ²) test, comparative fit index (CFI) and root square mean error of approximation (RMSEA). If the model is not acceptable, the model is further revised by removing or adding a connection between variables according to the modification indices and corresponding parameter changes expected. The process is performed continuously until finding a good-fitting model (see Annex 2b of Additional file 1).

The final step was to fit a multi-group SEM model for NCFs and PCFs. (1) The most constrained multi-group SEM was first performed, for which individual path coefficients were forced to be equal between groups (Shipley 2004). If the multi-group model was rejected (P < 0.05), indicating that at least one path is not equal between groups (Shipley 2004), and further analysis will be performed. Otherwise, if it was acceptable (P > 0.05), indicating that there is no any difference between groups; and no further analysis will be carried out. (2) By taking a full free multi-group model as a starting point, fitting a series of nested models to find which paths were equal or which paths differ across groups. For the nested models, only one path of the model could be constrained to be equal between groups at a time (Table 2). The difference in the maximum likelihood χ² statistics between models before and after adding an equal constrained path was used to test if this change is acceptable. (3) Multi-group SEM was fitted by constraining those path coefficients acceptable across groups. The SEM and multi-group SEM were implemented using the R version 3.4.4 and lavaan package version 0.6–3 (Rosseel 2012).

Table 2 Multi-group comparison of path coefficients between natural and planted coniferous forests. A non-significant value, highlighted in bold, indicates that the path contribution to the model is equal between forest types

Results

Biomass carbon density and storage

The mean BCD of NCFs was significantly larger than that of PCFs (26.46 Mg∙ha^− 1 versus 18.69 Mg∙ha^− 1), with an average BCD of 23.07 Mg∙ha^− 1 for all coniferous forests (Table 3). A significant difference was detected among the average BCDs of the four main types of forests (Table 3). For the four main types of forests: Pinus, Platycladus, Larix and Picea forests, the mean BCD for natural forests was significantly larger than that for planted forests for the first three types of forests, but not for the last one. The average BCD of natural and planted forests was the largest for Picea forests, followed by Larix and Pinus forests, and it was the smallest for Platycladus forests.

Table 3 Mean biomass carbon density (Mean BCD) and biomass carbon storage (BCS) of natural, planted and all coniferous forests in the study region

The BCD varied with forests age. For all the natural or planted coniferous forests, the BCD increased from young to near-mature forests, while it slightly decreased for mature forests (there was no over-mature forests) relative to near-mature forests (Table 4). Although the average BCD of NCFs appeared higher than that of PCFs for every age group, a significant difference was only detected for young and half-mature age groups (Table 4).

Table 4 Mean biomass carbon density (Mg∙ha⁻¹) of individual age groups of natural, and planted coniferous forests in the study region

It was also observed that the total amount of biomass carbon storage (BCS) for all the coniferous forests was 18.48 × 10⁶ Mg (Table 4). The BCS of NCFs was nearly two times that of PCFs (12.04 × 10⁶ Mg versus 6.44 × 10⁶ Mg), and accounted for 65.1% and 34.9% of the total BCS, respectively. For the four main types of forests, the Pinus forests had the largest contribution (65.9%) to the total BCS, followed by Larix and Picea forests (17.4% and 12.6%, respectively), and Platycladus forests had the smallest contribution (3.9%). Of the total BCS, half-mature forests accounted for 54.6%, and near-mature, young, and mature forests 28.6%, 10.5% and 6.3%, respectively (Annex 4 of Additional file 1).

Spatial autocorrelation and heterogeneity of biomass carbon density

Anselin Local Moran’s I analysis showed that the number of the high-high clusters of BCD for NCFs was more than twice that for PCFs (Fig. 2). The high-high cluster areas for the NCFs were identified in the Wutai, Guancen and Guandi Mountains. In contrast, the high-high cluster areas for the PCFs were observed in the Hengshan, Wutai, Guancen and Luya Mountains. The low-low clusters for NCFs were mainly distributed in the southernmost part of Lüliang Mountains, the southernmost part of Taihang Mountains, the south of Taiyue Mountains and the northeast of Zhongtiao Mountains, while the only one low-low cluster for PCFs was observed in the middle section of Lüliang Mountains. Besides, some high-low outliers were found in the Taiyue, and Zhongtiao Mountains, and the southernmost part of the Taihang Mountains for NCFs, but only two high-low outliers were observed in Guandi Mountain for PCFs. As a whole, the total number of clusters of BCD for NCFs was substantially larger than that for PCFs (81 and 14, respectively), accounting for 28.0% and 6.4% of the forest stands for them, respectively. Noticeably, the high-high and low-low cluster areas of BCD for NCFs were distributed in the north and south mountainous areas of Shanxi, respectively; while both high-high and low-low clusters of BCD for PCFs were distributed in the north mountainous areas of Shanxi.

Fig. 2

Maps showing the spatial clusters and outliers of biomass carbon density of natural and planted coniferous forests (NCFs and PCFs) in Shanxi Province. The number in parentheses is the number of sample plots (forest stands)

Moreover, the Getis-ord G* analysis showed that for NCFs, all the hot-spots of BCD were distributed in the north of 37°N in Shanxi (Fig. 3), including some areas of Luya, Guandi, Wutai and Taihang Mountains, while all the cold-spots of BCD were located in the south of 37°N in the province, including some areas of the southernmost part of Lüliang Mountains, the southernmost part of Taihang Mountains, the south of Taiyue Mountains and the northeast part of Zhongtiao Mountains. In contrast, for PCFs, there were a small number of hot- and cold-spots of BCD that were distributed in the north of 37°N in Shanxi. The hot-spots occurred in some areas of the Wutai, Hengshan and Luya Mountains, and the cold-spots in the middle section of the Taihang Mountains. This suggests that relative to the regional average, the BCD was generally higher in the north but lower in the south for NCFs, while the BCD was only higher or lower in some mountainous areas in the north of 37°N for PCFs over the province. The hot-spot and cold-spot patterns for NCFs and PCFs were in line with the high-high and low-low clusters revealed by Anselin Local Moran’s I analysis for them, respectively.

Fig. 3

Maps showing the spatial cold and hot spots of biomass carbon density of natural and planted coniferous forests (NCFs and PCFs) in Shanxi Province. The number in parentheses is the number of sample plots (forest stands)

Furthermore, semivariogram analysis showed a clear difference in the spatial heterogeneity of biomass carbon density between natural and planted forests (Table 5). The smaller value of nugget-to-sill for natural than for planted forests (0.692 versus 0.887) indicated that the spatial dependence of biomass carbon density was substantially higher for the former than the latter across the entire region. Noticeably, this analysis also showed a much higher spatial dependence of BCD for natural than for planted forests based on the data from the cluster and outlier analysis (sill = 0.5540 and 0.1592, respectively) and the hot spot analysis (sill = 0.6580 and 0.0427, respectively).

Table 5 Estimated parameters of semivariogram based on the data of biomass carbon density and from cluster and outlier analysis for natural and planted coniferous forests

Influences of various factors on biomass carbon density

The multi-group SEM model accounted for 56.0% and 57.0% of the variation in BCD for NCFs and PCFs, respectively (Fig. 4). ELE had the greater total effect on BCD than other driving factors for both NCFs and PCFs in Shanxi, with the total effect stronger for NCFs (0.409) than PCFs (0.339). The effects of ELE on BCD included the direct effect and indirect effect by affecting AGE for both NCFs and PCFs, with direct effect stronger than indirect effect for NCFs (0.348 versus 0.062) and PCFs (0.264 versus 0.076). Meanwhile, AGE, COV and PRCP also showed positive direct effects on BCD, and AGE and PRCP had indirect positive effects by affecting COV. The total effects of the three driving factors decreased in the order: AGE > COV > PRCP for both NCFs and PCFs. Moreover, LAT presented only indirect effects on BCD for NCFs and PCFs (Annex 5 of Additional file 1), with the total effect stronger for NCFs than PCFs. The result of multi-group analysis generally indicated that ELE and AGE were important driving factors of BCD for both NCFs and PCFs. Compared with other factors, the effects of both latitude and elevation on BCD showed much greater difference between NCFs than PCFs.

Fig. 4

Multi-group structural equation model for natural and planted coniferous forests. The factors influencing the biomass carbon density (BCD) of natural and planted coniferous forests (NCFs and PCFs) include stand age (AGE), stand coverage (COV), elevation (ELE), latitude (LAT) and annual precipitation (PRCP). Solid arrows indicate significant paths (P < 0.05), while dashed arrows are non-significant paths. The numbers along the arrows are standardized regression weights indicating the effect size (black, significant; grey, non-significant). The thickness of the arrow is proportion to the size of the effect. R² denotes the proportion of variation explained. Dark green arrows represent equal effect between NCFs (panel a) and PCFs (panel b); while dark red arrows represent unequal effect between NCFs (panel a) and PCFs (panel b)

Discussion

Increasing biomass carbon density and storage

The current study revealed that the BCD was significantly greater for natural than for planted coniferous forests in Shanxi. The higher BCD for the former than the latter can be primarily attributed to the distinct difference in stand age. The stand age for NCFs was nearly twice that for PCFs (at an average age of 54.1 and 28.5 years, respectively). For the largely same reason, the BCD was higher for each type of the natural coniferous forests than for its planted counterparts, because the stand age for the former was at least 1.7 times that for the latter. For example, the huge difference in BCD between natural and planted Platycladus forests (11.74 Mg∙ha^− 1 versus 1.95 Mg∙ha^− 1) was obviously related to the larger stand age for natural than for planted Platycladus forests (at an average age of 47.9 and 12.9 years, respectively).

Li and Lei (2010) have reported that natural and planted forests accounted for 83.05% and 16.95% of the total carbon stock across China, respectively. In our study, the representative area of a sample plot is 15.8 km², and the total area of natural and planted coniferous forests is 45.50 and 34.46 ha, respectively. The natural and planted coniferous forests accounted for 65.1% and 34.9% of the total carbon stock, respectively (Table 3). The large contribution of planted coniferous forests to the carbon stock was chiefly due to vast changes in land use in past decades. The implementation of the Three-North Shelter Forest Program and the Program for Conversion of Cropland into Forests has led to a huge increase in planted coniferous forests. Moreover, because majority of the natural coniferous forests and most of the planted coniferous forests (67.8% and 83.1%, respectively) were at young and half-mature ages under the current inventory, it is expected that the total carbon stock of coniferous forests will increase in future, and that planted coniferous forests will contribute more to future carbon stock in Shanxi.

Spatial patterns of biomass carbon density

The spatial pattern of BCD for NCFs was clearly different from that for PCFs. Although the spatial variation of BCD was related to natural and anthropogenic factors for both NCFs and PCFs, the spatial pattern of BCD for the former was likely shaped to a greater extent by natural influences, while that for the latter was probably formed to a greater extent due to anthropogenic influences. For NCFs, the spatial pattern of BCD was first closely related to the spatial distribution patterns of forest types across the region, and the associated difference in BCD among forest types. As shown in Fig. 1b, most of the natural Picea and Larix forests are distributed in the north Shanxi, while majority of the natural Pinus and Platycladus forests are located in the south Shanxi. More important, the BCDs for the former were much higher, while the BCDs for the latter were much smaller compared with the average BCD for NCFs as a whole (Tables 3 and 4). Meanwhile, the establishment of natural reserves for the protection of natural Picea and Larix forests might have enhanced the formation of the high-high clusters of BCD or hot spots of BCD in the core areas of natural reserves, where the forests had a continuous distribution and were well protected. On the other hand, the spatial distribution patterns of planted coniferous forests were to great extent different from those of their natural counterparts (Fig. 1b) as a result of reforestation and particularly afforestation. For instance, the planted Pinus forests, which accounts for 67.6% of all the planted forests, were distributed across the entire region, at lower elevations generally. Because most of the planted coniferous forests were still at the young and half-mature stages as above mentioned, the BCDs for planted Pinus and Platycladus forests were generally lower than the average BCD for PCFs as a whole (Tables 3 and 4), and only some of planted Picea and Larix forests had a higher than average BCD. Therefore, only a small number of the high-high clusters of BCD or hot spots of BCD were identified for planted coniferous forests in the north mountainous areas of Shanxi.

Effect of elevation versus temperature on biomass carbon density

We found that ELE was the most important driving factor of BCD among the five factors tested, and effect of ELE was greater for NCFs than for PCFs. Elevation gradient is associated with changes in temperature and precipitation, and forest type (Sanaei et al. 2018). ELE can affect forest canopy, stem density and stand basal area, and therefore affect aboveground biomass (Xu et al. 2018a, 2018b) through regulating moisture and soil water availability (Fisk et al. 1998). Based on the SEM, Xu et al. (Xu et al. 2018a, 2018b) reported that ELE was the most important abiotic driving factor of vegetation carbon stocks, while canopy density and forest age were the most crucial driving factors of vegetation carbon stocks in Zhejiang Province, China. The unusually large effect of ELE on BCD for NCFs in Shanxi was probably due to that the range of elevation (483 to 2560 m) over this mountainous region is much larger than that in Zhejiang Province (2 to 1530 m), and noticeably, the Picea and Larix forests are generally distributed at a higher elevation than Pinus and Platycladus forests across the entire region (Table 1); and the BCD for the former was generally greater than for the latter, as indicated by the much higher mean BCDs for the former, while much lower mean BCDs for the latter compared with the overall average BCD of natural coniferous forests (Table 4). In addition, the positive correlation of ELE with AGE and COV also enhanced the total positive effect of ELE on BCD. This was probably due to that the forests at higher elevations have usually experienced little disturbance, and can grow steadily and present higher BCD until the mature stage.

However, it was worthy note that we eliminated annual mean temperature (TEMP) in the multi-group SEM because of the collinear relationship between ELE and TEMP, and the stronger correlation of biomass carbon density with ELE than with TEMP (Annex 1 of Additional file 1). To examine further what the result could be if TEMP rather than ELE was selected, we performed the multi-group SEM for NCFs and PCFs again (Fig. 5 and Annexes 3, 6 and 7 of Additional file 1). It was interesting that similar result was found from the analysis except for that the effect of elevation was replaced by the effect of TEMP, which showed a significant direct negative effect on BCD for NCFs and PCFs (Fig. 5). The total effects of the driving factors were − 0.367 (TEMP), 0.244 (LAT), 0.200 (AGE), 0.177 (COV) and 0.110 (PRCP) for NCFs; and − 0.367 (TEMP), 0.248 (AGE), 0.177 (COV), 0.110 (PRCP) and 0.038 (LAT) for PCFs (Annex 7 of Additional file 1). This result suggested that TEMP could be to great extent taken as an alternative factor to ELE in the multi-group SEM for NCFs and PCFs. The negative effect of TEMP on BCD is consistent with the latest study in the Lüliang Mountains (Wang et al. 2018a).

Fig. 5

The multi-group structure equation model for natural and planted coniferous forests when annual mean temperature (TEMP) was taken an alternative factor to elevation (ELE). See Fig. 4. for details on the description of the model

Effects of other factors on biomass carbon density

AGE and COV were other important factors in modulating the BCD in NCFs and PCFs. The strong positive effect of AGE on BCD could be attributed to cumulative tree growth over time (Ali et al. 2016). Meanwhile, as hypothesized, we found that AGE was significantly positively related to COV, which had a strong direct effect on BCD. Our findings were in agreement with Xu et al. (Xu et al. 2018a, 2018b) who found that canopy density and forest age were the dominant determinants of vegetation carbon density in subtropical forest ecosystems. The pronounced effect of AGE on BCD was a good reflection of the process that the BCD increased with stand age from young to near-mature ages (Annex 4 of Additional file 1). Although the BCD declined slightly at the mature stage, the forests at this stage only accounted for a small proportion of the total forests.

Multi-group SEM revealed that LAT affected BCD indirectly and the total positive effect for NCFs was much greater than that for PCFs. Similarly, Guo and Ren (2014) reported that tree biomass had stronger latitudinal trend in natural forests than in planted forests across China. It was interesting that further analysis showed much stronger correlation, for NCFs than for PCFs, of LAT with the corresponding BCDs derived from Hot spot analysis (r = 0.893, P < 0.001; r = 0.361, P < 0.001) and Clusters and Outliers analysis (r = 0.543, P < 0.001; r = 0.215, P = 0.001). This implied that the relationship between LAT and BCD detected by multi-group SEM well reflected the spatial pattern of BCDs for NCFs and PCFs.

Different from the previous studies in which the relationships of latitude and other factors with vegetation carbon density were analyzed separately (Guo and Ren 2014; Wen and He 2016), our study incorporated LAT into SEM to assess the indirect effects of LAT via other driving factors. In this way, we found that the final total indirect effect of LAT on BCD was the combined effects of LAT on ELE, AGE and PRCP (Fig. 4). Therefore, the introduction of LAT into SEM could help us to understand the mechanism of effect of LAT on BCD. Noticeably, we found that the difference in the total effect of LAT between NCFs and PCFs was mainly due to the different effect of LAT on ELE and AGE (Fig. 4). The stronger positive effect of LAT on BCD for NCFs than PCFs was due to the larger positive effect of LAT on ELE for the former than the latter and no effect of LAT on AGE for NCFs but a large negative effect on AGE for PCFs. The larger difference in the effect of LAT on BCD between NCFs and PCFs was to some extent a mirror of the difference in the spatial patterns of BCDs between these two types of forests in the study region.

Conclusions

The present study found that the biomass carbon density and carbon stock were significantly greater for natural than planted coniferous forests in Shanxi. However, the planted coniferous forests will contribute more to future total carbon stock. Spatially, the higher biomass carbon density for the two forest types was clustered in the north mountainous areas, and the clustered number was larger for natural than planted coniferous forests. Meanwhile, the biomass carbon density of only natural coniferous forests was detected the clustered area in the south mountainous areas. Obviously, a much higher spatial dependence of biomass carbon density for natural than for planted forests. The effect of each driving factor on biomass carbon density was not always consistent between forest types. Elevation was the most important driving factor, and the effect on biomass carbon density was stronger for natural than planted coniferous forests. Latitude presented only indirect effects on biomass carbon density, and the difference between the two forest types was to some extent a mirror of the difference in the spatial patterns of BCDs between them in the study region.