Forest is one of the most widely distributed terrestrial ecosystem, which has a direct effect on absorption of CO₂ and other greenhouse gases and acts an essential indicator in global carbon cycle and climate regulation [1]. Forest biomass refers to the total amount of organic matter on a unit area at a certain moment, which accounts for about 85%-90% of total biomass on earth and thus is the most essential index indicating the productivity of the terrestrial ecosystem [2,3]. Forest biomass is not only an important indicator of forest carbon sink capacity, but also a major parameter of forest carbon budget [4- 6]. Although forest biomass carbon pool plays an indispensable role in the regional carbon budget, the understanding of the change of forest biomass carbon stock and its response to the climate change is still poor due to limited data and spatial heterogeneity [1,7,8]. Therefore, estimating forest biomass as well as its change is a substantial part of the study on carbon cycle and climate change.

Climate change, especially high temperature, drought, flood, fire and other extreme disasters can disturb the growth of forest [9-13]. The suitability of meteorological factors to the forest can be summarized as optimum, the upper limit and the lower limit. When the meteorological factors varying within the optimum range, the vegetation would grow to the best condition; whereas inhibited vegetation or even death of vegetation could occur upon meteorological factors near or beyond the upper or lower limit [14]. Due to the global climate change, extreme drought occurred more frequently in recent years, which directly influenced the productivity and composition of the vegetation. Studies have verified that the same vegetation has distinct response to different drought frequency and intensity [15], in addition, the response of vegetation which is different types but nearby to distinct frequency and intensity drought are also different [16]. Therefore, there is a clear need to study productivity change of forest vegetation in response to specific drought.

Severe drought occurred in five provinces in the southwestern China, which endured for a long time and deeply impacted a considerably large area. In particular, the extreme drought occurred in 2010 was a disaster to the agriculture, industry, city development and other economic and social activities in the area, resulting in death of forest vegetation and thus a serious blow of local ecological environment. In the past, many studies analyzed the impact of drought on vegetation in the region. Yang et al. [17] provided the most seriously affected area using the distribution of the site monthly precipitation anomaly in the past 50 years and the percentage of abnormal precipitation. This article comprehensively explored the reason for the drought and the change of the corresponding affected area over time. Zhao et al. [15] divided drought into different stages using Palmer drought severity index (PDSI) in the southwestern area from 2009-2010 and further studied the response of forest, grassland and savanna to different levels of drought. Severe drought can have different impacts on vegetation during different stage of drought. Zhang et al. [18] verified that the drought which occurred in the spring would make the enhanced vegetation index (EVI) and gross primary productivity (GPP) reduced, indicating that spring drought have an important impact on vegetation productivity and terrestrial carbon cycle. Although there are a handful studies on the impact of extreme drought on vegetation growth in the southwestern China, there are fewer studies focusing on the impact of drought especially years continuous drought on forest biomass carbon sink. Biomass is a major feature of terrestrial carbon cycle, therefore, evaluating the impact of extreme climate change on forest BCD is an important and urgent task. This study first calculated the actual biomass carbon density using the method of biomass conversion factors (BEF) and accompanied with forest volume and areas based on 3rd to 8th forest resources statistics of China. China forest BCD with high accuracy was then estimated using Inventory-satellite-based method in conjunction with remote sensing data NDVI and elevation. At last, we used PDSI as the drought index to analyze the correlation between BCD anomaly and PDSI from 2009 to 2013 and verified the impact of drought disturbance on forest biomass.

Study area

The study area is located in the southwestern China, including five provinces: Yunnan, Guizhou, Chongqing, Guangxi and Sichuan. Since 2009, the precipitation in this region was continuously below normal level, and a rare continuous drought occurred from 2009 to 2012, resulting in death of a large area of forest, which was explored by quite a few studies [15,17,18]. Figure 1 shows the distribution of different types of vegetation in this region. The red square area are our study area which stand for the most seriously affected part, and they are calculated using the monthly precipitation anomaly spatial distribution data in the past 50 years and the percentage of anomaly precipitation [17].

Figure 1: The location of the five provinces in southwestern China and the spatial distribution of land cover. The pixels are generated using unchanged pixels in land cover data.

NDVI data

The NDVI, which is defined as a normalized ratio of the near infrared and red bands, is widely used as a proxy of canopy greenness. The Global Inventory Monitoring and Modeling Systems (GIMMS3g) NDVI product was used to calculate NDVI at 0.083° spatial resolution and 15 day temporal resolution [19-21]. The method of Maximum Value Composition (MVC) was selected in the process of data synthesis, the data were corrected for the effect of atmospheric gases, thin cirrus clouds and aerosols. As NDVI can effectively reflect the biological and physicochemical characteristics of vegetation, such as biomass, coverage and chlorophyll content, it can be used to monitor the Earth’s terrestrial photosynthetic vegetation activity in support of vegetation change, climatic impact and biophysical interpretations [22].

Forest inventory data

The forest volume data were collected from forest inventory data for 1984-2013 in China. Sample areas in each province are periodically checked in China’s national forest inventory method. These forest resource data can be used to investigate the macroscopic status and dynamics of forest resources, which is important for the development and adjustment of forestry policy, planning and programs. Eight national forest inventories were completed from 1950 until 2013. The Ministry of Forestry completed the first national forest statistics in the 1960s. The second national forest inventory was conducted in 1977-1981, and the third and the fourth inventories were conducted in 1984-1988 and 1989-1993. The forest inventory then was conducted every five years (1994-1998, 1999-2003, 2004-2008, and 2009-2013). These inventory data provide a scientific and reliable basis for national dynamic monitoring and management of forest resources.

This research used 3^rd to 8^th forest resources statistics of China (1984-2013). The data of land vegetation included stand forest, economic forest and bamboo, among which stand forest could provide the area of the dominant species and the volume data while the economic and bamboo area data were estimated from the data of different provinces. It is worth notice that the definition of stand forest was changed in 1994: the canopy density for stand forest was adjusted from greater than 0.3 to no less than 0.2. In order to make consistent comparisons, we constructed a linear relationship between area and total carbon stock based on both standards [2,23]. To achieve more accurate conversions, we used power function conversion with canopy density equal to or greater than 0.2 as the uniform standard (as shown below) in the following analysis.

Area0.2=1.29^*Area_0.3 ^0.995 (R²=0.996, N=30) (1)

Carbon0.2=1.147^*Carbon_0.3^0.995 (R²=0.996, N=30) (2)

where Area_0.2 and Area_0.3 indicated the provincial stand forest area (in 10⁴ ha) with canopy density greater than 0.3 and 0.2, respectively; Carbon_0.2 and Carbon_0.3 indicated the biomass carbon stock (T_g C) under the two standards.

Forest type data

There are multiple ways to classify the land cover and MODIS classification data is the most important one [19]: the vegetation map of the People’s Republic of China with a scale of 1:1000000 [24]. In MODIS database, MCD12Q1 data is updated annually with 500 m spatial resolution. The data recorded from 2001 to 2011 reveal the dynamic change of vegetation, and therefore this study used MCD12Q1 data.

In order to reduce the error in the analysis of impacts of climate change on BCD due to the change of land classification, we selected the invariant pixels in the MODIS database corresponding to the vegetation types from 2001 to 2011. In the study area, the invariant pixels can be divided into four categories: mixed forest is major vegetation which accounts for 69.7%, savanna accounts for 7.9%, grassland accounts for 9.6% and cropland accounts for 9.0%. This study focused on exploring the impact of drought disturbance on forest biomass, so the main vegetation we focused on was mixed forest.

PDSI data

The quantification of drought is a difficult task, as we usually identify a drought by its effects on different systems (agriculture, water resources, ecosystem), but there is not a unique physical variable we can measure to quantify drought intensity [25]. The analysis of the spatial and temporal features of drought is the prerequisite to evaluate the response of forest ecological system to drought. However, the spacetime features of drought depend on specific assessment indicators, and the climate and vegetation factors are widely different in different assessment indicators. So the rationality of the assessment index is the prerequisite for an accurate assessment of drought and its impact. In the assessment of drought and its impact, the drought intensity, duration, frequency and area determine the effects of drought and are the core for regional drought analysis [26]. Drought is usually caused by a variety of complex physical mechanism and the indices for assessing drought are usually constructed in specific region and for a certain period of time with temporal and spatial characteristics [27].

In this article, we used the monthly scale with 0.5D spatial resolution as the PDSI to measure the drought level. The Palmer Drought Severity Index (PDSI) can be used to monitor drought conditions, calculating with precipitation, temperature, and soil moisture data [28,29]. According to the principle of soil water balance, droughts are defined to be continuous abnormal water deficit or water deficit due to the fact that actual water supply continues to be less than the local suitable water supply.

Forest BCD estimation

In order to analyze the relationship between meteorological factors and forest BCD, we need to first obtain the BCD of all forest pixels in China. First we counted the area and volume of every province from forest resources statistics of China. The ways to calculate the biomass of the three forest types are also different. The biomass density of economic forest is 23.7 Mg/ha (1 Mg=1000 kg or 10⁶ g). Bamboo can be divided into moso and other bamboo; biomass density of moso is 81.9 Mg/ha while biomass density of other bamboo is 53.1 Mg/ha [30]. For calculation of stand forest BCD, we used biomass expansion factor (BEF). The results showed that biomass can be calculated from forest volume data using BEF as shown below [31]:

B=aV=b (3)

Where B is the total biomass (Mg ha^-1); V stand for forest volume (m³ ha^-1); a and b are conversion coefficients of specific types of forest. In order to achieve a more accurate calculation of forest biomass and improve the potential deficiencies including the short-term statistics, without considering the influence of forest age and insufficient data [31,32], we updated the conversion coefficient (Table 1). In this case, the forest volume data can be more accurately converted to biomass using 0.5 as the carbon conversion coefficient of biomass. Previous studies have shown that multiple linear regression model of each province can be constructed in terms of the corresponding mean BCD, NDVI and the latitude and longitude coordinates. Moreover, elevation also has a direct impact on the vegetation and thus BCD. We therefore estimated the BCD of each pixel position using the multiple linear regression relationship between BCD and these variables.

Table 1: Biomass conversion coefficients [31].

Impact of drought

In order to evaluate the impact of the drought occurred from 2009 to 2013 on forest carbon stock and to detect the unusual change of BCD during this period, we analyzed the correlation between BCD indicating vegetation carbon absorption and PDSI indicating the drought level.

We have found that the 2009-2010 drought could be divided into different stages: first stage, middle stage and last stage and the middle stage (Dec 2009-Sep 2010) was the most serious one [15]. Therefore, in the analysis of 2009-2013 drought, we estimated the drought level from the PDSI averaged over the PDSI detected from Dec of previous year to Sep of the following year. The 2010 drought severely affected the vegetation in the five provinces in the southwestern China. However, the impact of drought was recovered in the following several years and the response of vegetation to drought also recovered somehow.

Z value is often used for anomaly detection. In order to show the unusual change of forest vegetation due to the drought, we calculated the Z value of BCD from 2009 to 2013 as shown below:

(4)

Where X is the NDVI, in our article that is BCD, Xn is the average NDVI, and σ (X) stands for the standard deviation. We used the mean value of BCD of 2001-2008 as the vegetation carbon absorption in normal years to calculate the anomaly from 2009 to 2013.

Biomass estimation

In this study, forest inventory data and corresponding remote sensing data were used to develop the forest biomass estimation model. First the BCD of each province was calculated by BEF and area conversion method. Then the average values of maximum NDVI for each province were calculated to establish the relationship between the average BCD and the NDVI data of each province. Previous studies showed that spatial location had an important impact on biomass when estimating remote BCD [5]. Location information was calculated as the average altitude and longitude of the province and was tested in the regression models. Elevation is another most important parameters impacting vegetation. The effect of elevation on biomass estimation was also tested by introducing the variable of elevation. The optimal regression equation was identified according to a stepwise regression method through the least squares regression analysis.

BCD= 7 4 6 . 0 6 3+ 8 6 . 8 3 7 ^*NDVI - 1 7 . 3 2 6 ^* l n (NDVI ) - 14.65^*lon+0.067^*lon²-0.347^*lat+0.0086^*dem

(R²=0.86, P<0.001) (5)

Forest inventory data in 1994-2013 were used to evaluate BCD estimation (Figure 2), and the R² is 0.86. Evaluated results show that the model is reliable in estimation of the China’s BCD.

Figure 2: Comparison between observed and estimated results of BCD. Model evaluation using forest inventory data in 1994-2013. In this figure, correlation coefficient between the two estimates was 0.86, with a good coincidence between the two estimates for most provinces.

Elevation data was mainly from the global climate data network (http://www.worldclim.org/), and the resolution was 500 m. The study area was a typical plateau mountainous region, in which the average elevation of the mixed forest was 1896 m.

BCD changes in the study area

To evaluate the spatial changes of forest biomass carbon density in our study area, we calculated mean of BCD from 2009 to 2013 (Figure 3). It can be seen that the biomass carbon density in our study area shows significant spatial differences. The average BCD in study area was 39.02 Mg C/ha. Figure 3 illustrates the distribution of BCD, indicating a highly spatial heterogeneity which reflects differences in ages and climatic conditions. The higher BCD occupied the northwestern of region, mean BCD is 50-70 Mg C/ha, while BCD relatively lower in the eastern coastal and central areas.

Figure 3: Spatial distribution of mean BCD from 2009-2013. The higher carbon density appeared in the northwestern while lower in the eastern region.

The average altitude of the study area is 1896 m, and we make two groups corresponding to lower than 1896 m and higher than 1896 m, then statisticed the frequency distribution of the mean BCD in 2000-2013 (Figure 4). From the distribution frequency of vegetation, the result of the independent sample T test show that the frequency distribution of BCD has no significant difference with altitude. The average value of BCD of all pixels below the mean altitude is 36.18 Mg C/ha, while the average value of BCD higher than mean altitude is 39.43 Mg C/ha.

Figure 4: BCD distribution frequency of vegetation at different altitude.

The effect of drought on BCD

In order to analyze the impact of drought of different levels on vegetation BCD, we calculated the mean PDSI from 2009 to 2013 (Figure 5) as well as the spatial distribution and frequency variation of BCD anomaly (Figure 6). In the evealing fact that remarkable PDSI box chart, the points from top to bottom represented the maximum value, 99% site, 75% site, mean value, 25% site, 1% site and the minimum value, respectively. The statistical analysis showed that, compared to 2009, the mean value of all PDSI in 2010 decreased by 1.292. The mean value of 2011 and 2012 PDSI pixels were still lower than that of 2009, which declined by 0.502 and 0.339, respectively, but increased compared to 2010. The mean value returned to the 2009 level until 2013. PDSI data clearly suggested that the drought occurred in 2010 was the most serious one and less severe drought still occurred in 2011 and 2012.

Figure 5: The figure shows the distribution of all invariant PDSI from 2009 to 2013. It is clear that PDSI mean value greater than 0 in 2009, which was at normal level.

Figure 6: BCD anomaly spatial distribution (left) and frequency distribution (right). BCD anomaly spatial distribution was calculated using Z value method and the data from previous Dec to Oct of current year.

The distribution and frequency of BCD Anomaly showed that the distribution of BCD in different years were significantly different. In 2009, the frequency of BCD Anomaly mainly distributed from -1 to 1 and the mean value was -0.129; the spatial distribution showed slight anomaly in the center of the study area. However, compared to 2009, BCD significantly declined in the southwestern part of the study area in 2010. The BCD Anomaly frequency mainly distributed from -4 to 1 and the mean value was -1.004. In 2011, the BCD in northwestern part recovered but the pixels indicating BCD anomaly in the southwestern part of the study area increased remarkably. The BCD Anomaly mainly distributed from -4 to 0 and the mean value was -0.959. In 2012, the overall situation recovered. In 2013, more than 72% of the pixels were greater than 0, the mean value of BCD Anomaly reached 0.371 and the vegetation growth completely returned to normal.

We implemented a correlation between PDSI and BCD anomaly in 2010 and the result was shown below. It is clear from the Figure 7 that BCD anomaly decreased with the decrease of PDSI and R2 reached 0.37 (P<0.001). This indicated that the response of forest was strong upon strong drought disturbance, which led to weakened vegetation carbon sequestration capacity and decline of BCD anomaly.

Figure 7: Linear relationship between PDSI and BCD Anomaly. The figure showed that BCD Anomaly declined with the decline of PDSI. This suggested that BCD Anomaly were more seriously affected by more severe drought.

In this work, we combined PDSI and BCD Anomaly to analyze the change of BCD Anomaly upon drought disturbance. First of all, in the biomass calculation of five provinces using BEF method, we adopted the most comprehensive conversion coefficients, which not only covers the longest time span, but also takes into account the age and segments of stand. Multiple nonlinear regression was then performed using remote sensing and survey data to improve the spatial resolution in calculation of biomass. At last, the linear analysis was implemented between PDSI and BCD Anomaly to detect the relationship between drought and BCD.

Biological carbon sequestration is an important biological process, in which vegetation absorb CO2 through photosynthesis and transform it into organic matter. Biomass is used as the index indicating the quality of organic matter per unit area of vegetation. Biological carbon sequestration can be influenced by internal factors such as age of stand forest, diversity of tree species as well as external disturbances. The driving force for the change of forest carbon stock is the dynamic changes of biomass due to the forest vegetation activity. Due to distinct external environmental factors, the carbon stock and carbon density can also be greatly different. Therefore, the study on the dynamic changes of forest carbon stock upon different levels of drought disturbance is important not only for evaluating the role of forest in regional-scale carbon cycle, but also for restoration, protection and management of forest.

Usually, drought, especially extreme drought will affect the growth of vegetation. Sergio M. Vicente-Serrano[31] also proves that the response of vegetation to drought depends on characteristic drought time-scales for different biomes. In the study of Zhang[14], climate and forest disturbances would impact the forest biomass to a certain extent. The same with these studies, our work provided valuable information and evidence for the relationships between forest biomass and drought disturbances.

In our estimation, we not only used the new conversion coefficients to calculate the stand biomass with thorough consideration of the age of stand and diversity of tree species, and also increased the biomass of bamboo and economic forest area based on the detailed table of bamboo and economic forest area, which achieved a more complete biomass data set of these provinces.

The NDVI value used in the remote sensing estimation model can be used to represent the best vegetation growth status, but it cannot represent the organic matter in the biological carbon cycle. NDVI in the remote sensing estimation model cannot be used in the assessment of vegetation biomass change upon drought disturbance. Moreover, NDVI only accounts for 12.8% of the biomass while the elevation accounts for 56%, geographical location accounts for 17.5%, and therefore NDVI alone cannot be used to represent BCD. Elevation is one of the important parameters affecting the growth of vegetation, which is also involved in our study. Compared to reference [5], in which the interpretation of all factors to BCD was 64%, the interpretation in this study is improved by 22% and the accuracy of BCD estimation is remarkably improved.

Figure 2 showed that remote sensing estimation is smaller with respect to the forest statistical data, which is due to the fact that the national statistical data of forest resources is mainly obtained by manual operation and thus a certain error is inevitable. Take, for example, Tibet, due to the large area and wide distribution of forests, human error may be inevitable in the survey.

In our remote sensing estimation model, we used the annual maximum value of NDVI to estimate the forest BCD, but there may be still certain uncertainty. When the age of tree is small, NDVI can somehow represent the changes in forest biomass. But due to the saturation limit of NDVI [39], it may not increase upon the growth of vegetation, but the biomass will further increase due to photosynthesis. This is one of the limitations in this study.

In this article, we used remote sensing method to estimate the biomass in five provinces in the southwestern China from 2009 to 2013, and then calculated the variation of biomass carbon density. The results showed that the accuracy can reach R2=0.86 using BEF to calculate the forest biomass followed by remote sensing estimation. After 2000, BCD of the forest vegetation in the study area increased with the average annual increase rate of 0.102 Mg C /ha, but it remarkably reduced from 2009 to 2012 and especially plunged from 2009 to 2010.This suggests that the forest carbon sequestration capability weakened and forest carbon stock declined during this time, which almost returned to normal in 2013. PDSI is an important index indicating drought level. PDSI in the study area declined from 2009 - 2010, was lower than normal level in 2011 - 2012 and recovered in 2013. BCD also showed anomaly during this time, and in particular BCD anomaly was serious in 2010. Correlation analysis between BCD anomaly and PDSI showed that BCD anomaly decreased upon decrease of PDSI, indicating a more serious drought and BCD anomaly that lead to greater impact on BCD and thus forest carbon cycle.

This work was supported by the High Technology Research and Development Program of China (Grant No. 2013AA122801) and the Fundamental Research Funds for the Central Universities (Grant No. 2014KJJCB17).

The analysis was performed by XZ and HW. All authors contributed with ideas, writing and discussions.