Spatiotemporal Evolution of the Carbon Fluxes from Bamboo Forests and their Response to Climate Change Based on a BEPS Model in China

Abstract

Carbon flux is the main basis for judging the carbon source/sink of forest ecosystems. Bamboo forests have gained much attention because of their high carbon sequestration capacity. In this study, we used a boreal ecosystem productivity simulator (BEPS) model to simulate the gross primary productivity (GPP) and net primary productivity (NPP) of bamboo forests in China during 2001–2018, and then explored the spatiotemporal evolution of the carbon fluxes and their response to climatic factors. The results showed that: (1) The simulated and observed GPP values exhibited a good correlation with the determination coefficient (R²), root mean square error (RMSE), and absolute bias (aBIAS) of 0.58, 1.43 g C m⁻² day⁻¹, and 1.21 g C m⁻² day⁻¹, respectively. (2) During 2001–2018, GPP and NPP showed fluctuating increasing trends with growth rates of 5.20 g C m⁻² yr⁻¹ and 3.88 g C m⁻² yr⁻¹, respectively. The spatial distribution characteristics of GPP and NPP were stronger in the south and east than in the north and west. Additionally, the trend slope results showed that GPP and NPP mainly increased, and approximately 30% of the area showed a significant increasing trend. (3) Our study showed that more than half of the area exhibited the fact that the influence of the average annual precipitation had positive effects on GPP and NPP, while the average annual minimum and maximum temperatures had negative effects on GPP and NPP. On a monthly scale, our study also demonstrated that the influence of precipitation on GPP and NPP was higher than that of the influence of temperature on them.

1. Introduction

Dynamic change in the carbon cycles of terrestrial ecosystems is a core component of climate change and regional sustainable development [1]; it plays an important role in the global carbon balance. Because of the impacts of various environmental and biological factors (such as climate change, vegetation distribution, and land-use change), the carbon cycles of terrestrial ecosystems show significant spatial heterogeneity [2]. Forest ecosystems are an important component of terrestrial ecosystems and play an important role in improving and maintaining the ecological environment, in addition to regulating the global carbon balance [3,4]. Therefore, it is essential to quantify carbon fluxes in forest ecosystems and explore their response to environmental factors in the carbon cycles of terrestrial ecosystems. Carbon flux [5] is the basis for determining the carbon source/sink of an ecosystem and plays an important role in the global carbon cycle and carbon balance. Ecosystem productivity is an important indicator for quantitatively describing the carbon sequestration capacity of an ecosystem, which mainly includes gross primary production (GPP) and net primary production (NPP) [6]. GPP refers to the amount of organic carbon fixed by photosynthesis per unit time and unit area of green plants [7,8]. It reflects the carbon sequestration ability of vegetation, and is the largest carbon flux in terrestrial ecosystems [9]. NPP is the organic matter or energy remaining for vegetation growth after deducting the organic matter consumed by vegetation autotrophic respiration (RA) on the basis of GPP [10]. It can directly reflect the production capacity and ecological environment quality of surface vegetation in the natural environment, and is an important indicator for evaluating the carbon sink of ecosystems, in addition to regulating ecological processes [11,12].

The methods to obtain carbon fluxes of forest ecosystems mainly include sample site inventory, eddy covariance technology, and model simulation. Sample site inventory can estimate carbon fluxes more accurately, but the estimation of the carbon fluxes of forest ecosystems requires long-term field measurements, which consume substantial amounts of time and labor [13]; therefore, it limits the estimation of carbon fluxes in forest ecosystems. Eddy covariance technology has the advantages of being a long-term, continuous, and non-destructive method [14], and has been widely used to estimate the carbon fluxes of forest ecosystems [15]. However, the number of flux observation sites is limited and the area of effective observation is very small. There are uncertainties in expanding it from the site to the regional scale, and it is often hindered by topography and climate conditions; therefore, eddy covariance technology has limitations in studying carbon fluxes in forest ecosystems at the regional scale.

Model simulation is an important method of evaluating carbon fluxes in forest ecosystems [16,17]. Remote sensing for earth observation technology has the characteristics of real-time, dynamic, and large-area synchronous monitoring, in addition to rich information [13]. It readily records the dynamic changes in environmental conditions, vegetation distribution patterns and activities, and land use in the form of electromagnetic information. This provides the necessary parameters of vegetation (such as NDVI and LAI) and environmental variables for the carbon flux model, and becomes a powerful method to study the distribution, seasonal change, and interannual change in carbon fluxes [18]. Therefore, the application of remote sensing data in the model estimation helps achieve cross-scale simulation of the carbon cycle process and reflect the spatial distribution and dynamic changes in the carbon budget at the regional and global scales. It increases the reliability and operability of vegetation carbon flux estimation and has become an important research topic [18,19,20,21,22]. Ecological process models simulate the effects of biological vegetation processes such as canopy photosynthesis, absorption, transpiration, and changes in soil moisture content on carbon fluxes, and have become an important method for carbon flux simulation. Common ecological process models include the Biome-BGC model [23], the BEPS model [24], and the InTEC model [25]. In recent years, scholars have studied the carbon fluxes of forest ecosystems in different regions using different ecological process models combined with remote sensing data. For example, Du et al. [26] used an improved Biome-BGC model with remote sensing data to simulate the above-ground carbon storage of bamboo forests in Zhejiang Province from 2003 to 2014, and analyzed its spatiotemporal patterns and influencing factors. Zhang et al. [27] used remote sensing data and the BEPS model to study the spatiotemporal distribution characteristics of GPP and NPP in terrestrial ecosystems in East Asia. Zheng et al. [28] used the InTEC model to simulate the NEP of the forests in Zhejiang Province during 1985–2015, and analyzed the response of climatic factors such as temperature, precipitation, relative humidity, and radiation.

The BEPS model is an ecological process model based on the FOREST-BGC model [29]. It integrates multi-source data as model inputs and is a good choice for simulating terrestrial ecosystem productivity with higher accuracy on larger spatial scales. The BEPS model successfully solved the problem of spatiotemporal scale conversion by using remote sensing data. Additionally, it solved the overestimation problem of the FOREST-BGC model by introducing the clumping index and advanced canopy transmission model [30]. Compared with other ecological process models, it may have the most potential to adequately address the spatiotemporal dynamics of carbon fluxes because of its strong theoretical basis and practical applicability [31]. Previously, it was used to simulate the productivity of the boreal forest ecosystem in Canada [32]. Presently, the model has been frequently modified and improved, and has been widely used to simulate the carbon fluxes of different regional terrestrial ecosystems at various spatiotemporal scales [33,34,35,36,37,38].

Bamboo belongs to a family of perennial graminaceous plants. There are approximately 150 genera and 1225 species of bamboo forests in the world, and the total area of bamboo forests worldwide accounts for more than 30 million ha [39], making it “the second largest forest in the world”. China is located in the center of bamboo distribution in the world. It has the richest bamboo resources in the world in terms of the number of species (more than 500 varieties of 39 species) and area [40]. According to the ninth National Forest Resources Inventory (2014–2018), China’s bamboo forest area is 6,411,600 ha [41], accounting for approximately 20% of the world’s bamboo forest area. Compared with the eighth National Forest Resources Inventory (2009–2013), their area of bamboo forest has increased by more than 400,000 ha. It is known as the “Bamboo Kingdom” [40,42]. Bamboo forests have a great carbon sequestration capacity and differ from other forests in mitigating climate change, and their impact on global climate change has become an important concern [43,44]. Several scholars have explored carbon cycling in bamboo forests, and have synthesized information concerning primary production [45,46], carbon stocks [47,48], and biomass [49,50].

Although relevant studies on bamboo forests’ carbon cycles have been conducted, the characteristics of bamboo forests’ carbon dynamics and their response to changing environmental conditions are still poorly understood [47,51,52]. Previous studies mainly focused on estimating carbon fluxes at the scale of sites, regions, and provinces, while relatively few studies have done so at the national scale. In addition, some studies lack simulations of physiological and ecological processes, leading to massive errors in the estimated results [53]. Therefore, the study of carbon fluxes from bamboo forests in China is essential for the study of the carbon cycles of forest ecosystems under the global climate background. The objectives of this study include (1) driving the BEPS model to simulate the carbon fluxes of bamboo forests in China from 2001 to 2018; (2) exploring the spatiotemporal evolution of bamboo forests’ carbon fluxes in China and the driving influence of climate change on carbon fluxes of bamboo forests in China.

2. Materials and Methods

2.1. Study Area

China has a vast territory and diverse climate types (Figure 1). The country has a north–south temperature gradient and an east–west precipitation gradient driven by the summer monsoon [54]. Bamboo forests are a unique and important forest type in subtropical regions of China, and are widely distributed across Zhejiang, Fujian, Jiangxi, Hunan, Sichuan, Anhui, Hubei, Guangdong, Guangxi, and other provinces.

2.2. Flux Measurement Sites

The flux observation sites are located in Zhejiang Province, which are the Anji Moso bamboo flux measurement site (30.46° N, 119.66° E) and the Lin’an Lei bamboo flux measurement site (30.30° N, 119.58° E) (Figure 1a). The height of the observation tower at Anji was 40 m, and the vegetation type around the flux tower for 1 km × 1 km was dominated by 1–4-year-old Moso bamboo forests. The height of the observation tower at Lin’an was 20 m, and the vegetation types around the flux tower were mainly 2–3-year-old Lei bamboo forests. The carbon flux data were continuously measured by an eddy covariance system of flux measurement sites. The system consists of an open-path infrared CO₂/H₂O gas analyzer (Li-7500, LiCor Biosciences Inc., Lincoln, NE, USA) and a three-dimensional sonic anemometer (CAST3, Campbell Scientific Inc., Logan, UT, USA). According to the principle of the eddy covariance system, 30-min carbon flux data were calculated online and stored.

2.3. Data Acquisition and Processing

The required BEPS model input data included bamboo forest information in China, MODIS leaf area index (MODIS LAI) (

Table 1. MODIS data and descriptions.

MODIS	Abbreviation	Time	Spatial Resolution	Time Resolution	To Use
MOD13A2	NDVI	2018	1000 m	16 days	Extract the bamboo forest
MOD09A1	REF	2018	500 m	8 days	Extract the bamboo forest
MOD15A2	LAI	2001–2018	1000 m	8 days	Model input

), the available soil water-holding capacity (AWC), daily meteorological data, and the biological parameters of bamboo forests. All data were reprojected to the WGS84 coordinate system with a spatial resolution of 1 km.

2.3.1. MODIS Data and Preprocessing

MODIS is a new generation of optical and infrared remote sensing instruments that “integrate image and spectrum” in the current world. It is widely used in the carbon cycles of terrestrial ecosystems because of its high time and spectral resolutions. This study uses MODIS normalized difference vegetation index (MODIS NDVI) (MOD13A2), MODIS land surface reflectance (MODIS REF) (MOD09A1), and MODIS LAI (MOD15A2) from NASA (https://ladsweb.modaps.eosdis.nasa.gov, accessed on 13 May 2020) to extract information on bamboo forests in China and simulate the carbon fluxes of bamboo forests in China. The MODIS data are shown in Table 1.

The MODIS Reprojection Tool (MRT) was used to preprocess MODIS data, such as mosaicking, format conversion, reprojection, and resampling. MOD09A1 was reprojected to the WGS84 coordinate system, and the spatial resolution was resampled to 1 km using the nearest neighborhood method. After resampling, these data were clipped to the boundaries of China.

2.3.2. Bamboo Forest Distribution Data of China

The distribution information of Chinese bamboo forests in 2003, 2008, 2014, and 2018 was extracted. The information on Chinese bamboo forests from 2003, 2008, and 2014 has been extracted in our previous study [42]. On this basis, we extracted information on Chinese bamboo forests from 2018. The flow chart of bamboo extraction is shown in Figure 2.

The main process is as follows: First, a total of 23 multi-temporal MODIS NDVI images are available. In order to further improve the MODIS NDVI data quality, these 23 images were composited into 12 multi-temporal images by selecting a maximum of two corresponding pixels of two neighboring MODIS NDVI images as the value of a new pixel (NDVI_max12) [55]. Then, a minimum noise fraction (MNF) transform [56] was employed to convert the NDVI_max12 data to obtain the principal component variables of NDVI _max12 data (NDVI _max12 MNF), and the first six bands with a cumulative contribution rate greater than 90% (NDVI_max12 MNF_1-6) were retained for classification. Second, according to the image texture and spectral information features, the five types of samples (forest, farmland, water, bare land, and residential land) were selected by visual interpretation [57], and then the study area was classified by the maximum likelihood classification (MLC). On this basis, the forest information in China was extracted by masking. Third, using the forest information in China to extract the normalized difference vegetation of forests (NDVI_{forest_12}) and the land surface reflectance of forests (REF_{forest_7}), MNF was then performed on them to obtain the principal component variables of NDVI_{forest_12} data (NDVI_{forest_12} MNF) and REF_{forest_7} data (REF_{forest_7} MNF). We retained the bands with a cumulative contribution rate greater than 85%, that is, the first nine bands of NDVI_{forest_12} MNF data (NDVI_{forest_12} MNF_1-9) and the first five bands of REF_{forest_7} MNF data (REF_{forest_7} MNF_1-5). On this basis, according to the training samples of bamboo forests, broad-leaved forests, and coniferous forests, the corresponding attribute values were extracted as the characteristic variables to construct a decision tree model, and the information on Chinese bamboo forests was extracted by using the constructed decision tree model (Figure 3). In this study, 85 bamboo forest survey samples from Zhejiang Province in 2019 and 440 bamboo forest samples from China selected from Landsat 8 images in 2018 were used as bamboo forest verification samples for point-by-point verification. Finally, the least-squares mixed-pixel decomposition method [58] was used to obtain the abundance information of bamboo forests in China. The results were presented in Figure 1a.

It was verified that the accuracy of bamboo forest extraction was 76.54–81.56%, and that the extracted area was close to the inventory area of forest resources (

Table 2. Extraction accuracy evaluation and the comparison of estimated and inventory bamboo forest area of China.

Year	Classification Accuracy Evaluation				Bamboo Forest Area (104 ha)
	Bamboo Forest Samples	Correctly	Incorrectly	User’s Accuracy (%)	Estimate	Inventory
2003	387 [42]	309	78	79.84	486.56	495.32 [42]
2008	414 [42]	328	86	79.23	545.14	548.73 [42]
2014	536 [42]	435	101	81.16	639.22	610.65 [42]
2018	525	402	123	76.54	669.83	656.08 [41,42]

Note: the results of the ninth National Forest Resources Survey do not have data from Taiwan, so the bamboo forest area of Taiwan is based on the results of the eighth National Forest Resources Survey.

), which laid the foundation for the simulation of GPP and NPP for bamboo forests in China. Only the information of bamboo forests in 2003, 2008, 2014, and 2018 was extracted; therefore, the bamboo forest information of the unclassified year was replaced by a similar year from which the bamboo forest information was extracted.

2.3.3. MODIS LAI Data

Leaf area index (LAI) is an important input parameter for simulating the carbon cycles of forest ecosystems, and is closely related to the photosynthesis, steaming, water utilization, and productivity formation of vegetation [59]. Remote sensing technology is an important method for obtaining a large-scale LAI. However, MODIS LAI data are susceptible to the influences of factors such as the atmosphere, which leads to an irregular reduction in data. To reduce data noise and improve data quality, the locally adjusted cubic-spline capping (LACC) [60] algorithm was used to smooth the clipped MODIS LAI data. Then, the smoothed MODIS LAI data were assimilated by the particle filter (PF) algorithm [61]. The assimilated MODIS LAI data were shown in Figure 1b.

2.3.4. Soil Data

The soil texture data map was provided by the Chinese Academy of Sciences (http://www.soil.csdb.cn, accessed on 11 December 2020). AWC is an important factor in terms of plant growth, affecting stomatal conductance and photosynthesis [62]. In this study, based on the empirical relationship, an AWC map with a 1 km resolution was obtained from a soil data thematic map. The spatial distribution of AWC was shown in Figure 1c.

2.3.5. Meteorological Data

Meteorological data from 2001 to 2018 were obtained from the National Meteorological Information Center of the China Meteorological Administration (http://data.cma.cn, accessed on 7 August 2020), and mainly included minimum temperature (T_min), maximum temperature (T_max), precipitation (Pre) solar radiation, and relative humidity. These meteorological factors are the main environmental factors in the carbon–water cycles [63]. The inverse distance weighting method was used to interpolate the data of 824 meteorological sites observed in the study area (Figure 1a) into spatial data with a 1 km resolution to obtain the grid cells of daily scale meteorological data of the study area. They are shown in Figure 1d–f. Among them, the temperature was corrected by the digital elevation method, and it was assumed that the temperature decreased by 6.5 °C for each per-kilometer increase in altitude. Solar radiation was simulated based on the measurements of sunshine duration at each site, following the methods of Ju et al. [64]. Monthly and annual meteorological data based on interpolated daily scale meteorological data were obtained.

2.3.6. Biological Parameters

The major biological parameters of bamboo forests used in the BEPS model are shown in

Table 3. Major biological parameters used as inputs into the BEPS model for simulating the CO₂ fluxes of bamboo forests.

Symbol	Unit	Description	Value	Reference
Ω	-	Clumping index	0.5	Measurement
Sarea		Specific leaf area	27	Measurement
Vm,25	umol m−2s−1	Maximum carboxylation rate at 25 °C	50	Iteration
Q10,leaf	-	Q10 for leaf	1.4	Iteration
Q10,stem	-	Q10 for stem	1.3	Iteration
Q10,root	-	Q10 for root	1.2	Iteration
Mleaf	kg C m−2	Average carbon storage of leaf	0.15	[66]
Mstem	kg C m−2	Average carbon storage of stem	1.76	[66]
Mroot	kg C m−2	Average carbon storage of root	1.15	[66]

. The clumping index (Ω) and specific leaf area (S_area) came from the measured data of the flux observation station. The maximum carboxylation rate at 25 °C (V_m) and the Q₁₀ for leaves, stems, and roots were calculated based on an iteration method. The initial value of the four parameters was established according to Chen et al. [65], and the iteration range for each parameter was set as ± 100%. The iteration step was defined as 1 for V_m and 0.1 for the other three parameters. The average carbon storage of leaves, stems, and roots was calculated using the methods of Zhou and Jiang [66]. Bamboo forests are a special type of forest. The photosynthesis capacity of bamboo forests is similar to C₃ trees [67]. Therefore, for constant parameter values, we referred to Feng et al. [62] to simulate the carbon cycle of bamboo forests.

2.4. BEPS Model Simulation and Evaluation

2.4.1. BEPS Model Description

The BEPS model is mainly composed of four parts: energy transmission, carbon cycle, water cycle, and physiological regulation sub-models [68]. It combines ecology, plant physiology, meteorology, and other disciplines to simulate the relationship between the photosynthesis, respiration, carbon distribution, water balance, and energy balance of vegetation [63], which demonstrates the combination of remote sensing data and ecological process models. The main feature of this model is that the instantaneous Farquhar photosynthetic model at the leaf scale is converted into the daily total photosynthetic model through the integration of stomatal conductance to realize the time scale expansion. Then, according to the principle of light transmission in the canopy, the vegetation canopy leaves were divided into shaded and sunlit leaves to simulate the radiation budget of the corresponding leaves. This helps achieve the expansion from the leaf scale to the canopy space scale. Detailed descriptions of the BEPS model can be found in Liu et al. [69] and Chen et al. [65]. The main simulation process of the model is as follows:

(1) The LAI_sunlit and LAI_shade are calculated as follows:

where LAI_canopy is the total LAI of the canopy; LAI_sunlit and LAI_shade are the canopy LAIs of sunlit and shaded leaves, respectively; θ is the daily mean solar zenith angle; and Ω is the clumping index:

$${\mathrm{LAI}}_{\mathrm{sunlit}}=2\cos \mathsf{\theta }[1-\exp (-0.5\mathsf{\Omega }\mathrm{LAI}/\cos \mathsf{\theta }\left)\right]$$

(1)

$${\mathrm{LAI}}_{\mathrm{shade}}{=\mathrm{LAI}}_{\mathrm{canopy}}{-\mathrm{LAI}}_{\mathrm{sunlit}}$$

(2)

(2) The photosynthesis rate is calculated as follows:

$$\mathrm{A}=\min ({\mathrm{W}}_{\mathrm{c}},{\mathrm{W}}_{\mathrm{j}})-{\mathrm{R}}_{\mathrm{d}}$$

(3)

$${\mathrm{W}}_{\mathrm{c}}={\mathrm{V}}_{\mathrm{m}}\frac{{\mathrm{C}}_{\mathrm{i}}-\Gamma }{{\mathrm{C}}_{\mathrm{i}}+{\mathrm{K}}_{\mathrm{c}}(1+{\mathrm{O}}_2/{\mathrm{K}}_{\mathrm{o}})}$$

(4)

$${\mathrm{W}}_{\mathrm{j}}=\mathrm{J}\frac{{\mathrm{C}}_{\mathrm{i}}-\Gamma }{4({\mathrm{C}}_{\mathrm{i}}+2\Gamma )}$$

(5)

$${\mathrm{R}}_{\mathrm{d}}=0.015{\mathrm{V}}_{\mathrm{m}}$$

(6)

where A is the net photosynthesis rate; W_c and W_j are the Rubisco-limited and RuBP-limited gross photosynthesis rates, respectively; R_d is the daytime leaf dark respiration; V_m is the maximum carboxylation rate at 25 °C; C_i and O₂ are the intercellular CO₂ and oxygen concentrations in the atmosphere, respectively; Γ is the CO₂ compensation point, without dark respiration; K_c and K_o are the Michaelis–Menten constants for CO₂ and O₂, respectively; and J is the electron transmission rate.

(3) The total canopy photosynthesis rate is evaluated as follows:

$${\mathrm{A}}_{\mathrm{canopy}}={\mathrm{A}}_{\mathrm{sunlit}}{\mathrm{LAI}}_{\mathrm{sunlit}} +{\mathrm{A}}_{\mathrm{shade}}{\mathrm{LAI}}_{\mathrm{shade}}$$

(7)

where A_canopy is the total photosynthesis rate of the canopy; A_sunlit and A_shade are the photosynthesis rates of sunlit and shaded leaves, respectively; and LAI_sunlit and LAI_shade are the LAIs of sunlit and shaded leaves, respectively.

(4) The GPP and NPP values are determined as follows:

$$\mathrm{GPP}={\mathrm{A}}_{\mathrm{canopy}} \times {\mathrm{L}}_{\mathrm{day}} \times {\mathrm{F}}_{\mathrm{GPP}}$$

(8)

$$\mathrm{NPP}=\mathrm{GPP}-{\mathrm{R}}_{\mathrm{a}}$$

(9)

$${\mathrm{R}}_{\mathrm{a}}={\mathrm{R}}_{\mathrm{m}}+{\mathrm{R}}_{\mathrm{g}}={\mathrm{R}}_{\mathrm{m},\mathrm{i}}+{\mathrm{R}}_{\mathrm{g},\mathrm{i}}$$

(10)

where GPP is gross primary productivity; NPP is net primary productivity; R_a is the autotrophic respiration of the vegetation; L_day is the length of the day; F_GPP is a scale factor for converting photosynthesis into GPP; R_m and R_g are the maintenance breathing rate and growth respiration rates, respectively; i is the different parts of vegetation (i = 1, 2, and 3 for leaves, stems, and roots, respectively); and R_m,i and R_g,i are the maintenance and growth respiration rates of different parts, respectively.

2.4.2. Evaluation of Simulation Results

In this study, the results of the BEPS model simulation were evaluated with precision using the determination coefficient (R²), root mean square error (RMSE), and absolute bias (aBIAS). The formulas for the calculation are as follows [13]:

$${\mathrm{R}}^2=1-\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{m}}_{\mathrm{i}}-{\mathrm{o}}_{\mathrm{i}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{o}}_{\mathrm{i}}-\overline{{\mathrm{o}}_{\mathrm{i}}}\right)}^2}$$

(11)

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{m}}_{\mathrm{i}}-{\mathrm{o}}_{\mathrm{i}})}^2}$$

(12)

$$\mathrm{aBIAS}=\frac{1}{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{m}}_{\mathrm{i}}-{\mathrm{o}}_{\mathrm{i}}\right|$$

(13)

where m_i is the simulated value; o_i is the observed value; and $\overline{{\mathrm{o}}_{\mathrm{i}}}$ is the average value of the observed value. Generally, the larger the R² value, the smaller the RMSE and aBIAS values, the higher the accuracy, and vice versa.

2.5. Spatiotemporal Evolution Analysis of Carbon Fluxes

2.5.1. Variation Coefficient of Carbon Fluxes

The variation coefficient (CV) is the ratio of the standard deviation to the average, which reflects the stability of a set of data. The higher the value of the CV, the more unstable the data, that is, the greater the fluctuation, and vice versa. To analyze the spatial fluctuations in carbon fluxes of bamboo forests during 2001–2018, the CVs of the GPP and NPP of each pixel were calculated as follows [70]:

$$\mathrm{CV}=\sqrt{\frac{\frac{1}{\mathrm{n}-1}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{P}}_{\mathrm{i}}-\overline{\mathrm{P}})}^2}{\overline{\mathrm{P}}}}$$

(14)

In Equation (14), CV is the variation coefficient; n = 18, and is the number of monitoring years; P_i is the value of each pixel of the GPP or NPP image in the i-th year (where i = 1, 2,…, n); and $\overline{\mathrm{P}}$ is the average value of each pixel of GPP or NPP. According to the calculation results, by performing the Jenks natural breaks classifications in ArcGIS software [71] the results of the CV were divided into five levels: low fluctuation (CV <= 0.1246), lower fluctuation (0.1246 < CV ≤ 0.2342), medium fluctuation (0.2342 < CV ≤ 0.4132), higher fluctuation (0.4132 < CV ≤ 0.7364), and high fluctuation (CV > 0.7364).

2.5.2. Trend Slope of Carbon Fluxes

To quantitatively study the trends of carbon fluxes of bamboo forests in China from 2001 to 2018 a linear regression analysis was used to calculate the trends of GPP and NPP of each pixel, as follows [72,73]:

$$\mathrm{slope}=\frac{\mathrm{n} \times {{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left(\mathrm{i}\times {\mathrm{P}}_{\mathrm{i}}\right) - {{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{i} \times {{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{P}}_{\mathrm{i}}}{\mathrm{n} \times {{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{i}}^2 - {\left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{i}\right)}^2}$$

(15)

In Equation (15), slope is the trend slope; n = 18, and is the number of monitoring years; and P_i is the GPP or NPP of bamboo forests in the i-th year, (i = 1, 2,…, n). The value of the trend slope indicates the rate of increase or decrease. When slope > 0, the GPP and NPP increase, and when slope < 0, the GPP and NPP decrease.

To analyze whether the variation trend of the GPP and NPP was significant, the F-test was used to test the significance of the variation trend of GPP and NPP. The variation trend was divided into five levels: significantly reduced (slope < 0, p < 0.01), reduced (slope < 0, 0.01 < p < 0.05), basically stable (p > 0.05), increased (slope > 0, 0.01 < p < 0.05), and significantly increased (slope > 0, p < 0.01).

2.6. Analysis of Spatiotemporal Responses of Carbon Fluxes to Climate Change

2.6.1. Partial Correlation Analysis of Carbon Fluxes to Climate Change

A correlation analysis reveals the closeness of the relationship between the study variables. Partial correlation analysis refers to the calculation of the correlation between two variables without considering the influence of other variables [74]. Partial correlation analysis can better reflect the impact of a single climate factor on carbon fluxes. Therefore, this study uses a pixel-based partial correlation analysis to calculate the partial correlation coefficients (PPCs) of GPP and NPP with climatic factors, and analyze the response between carbon fluxes and climatic factors. To determine the PCCs, we first calculated the correlation coefficient using the following formula [75]:

$${\mathrm{R}}_{\mathrm{xy}}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left[\right({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\left)\right({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\left)\right]}{\sqrt{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}})}^2{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}})}^2}}$$

(16)

In Equation (16), R_xy is the correlation coefficient between variables x and y; n is the number of study years; xi and yi represent the values of variables x and y in the i-th year, respectively, (i = 1, 2,…, n); and $\overline{\mathrm{x}}$ and $\overline{\mathrm{y}}$ represent the mean value of variables x and y, respectively. The range of R_xy is [–1,1]; when R_xy > 0 the two variables are positively correlated, and when R_xy < 0 the two variables are negatively correlated. The larger the absolute value of R_xy the higher the correlation, and vice versa. Generally, 0.3 and 0.6 are the distinction points of the absolute value of the correlation coefficient, as the weak correlation (0 < |R_xy| ≤ 0.3), low correlation (0.3 < |R_xy| ≤ 0.6), and significant correlation (0.6 < |R_xy| ≤ 1).

Based on the evaluated correlation coefficient, the PPC was calculated as follows [72]:

$${\mathrm{R}}_{\mathrm{ab},\mathrm{cd}}=\frac{{\mathrm{R}}_{\mathrm{ab},\mathrm{d}}-{\mathrm{R}}_{\mathrm{ac},\mathrm{d}} \times {\mathrm{R}}_{\mathrm{bc},\mathrm{d}}}{\sqrt{(1-{\mathrm{R}}_{\mathrm{ac},\mathrm{d}}^2) \times (1-{\mathrm{R}}_{\mathrm{bc},\mathrm{d}}^2)}}$$

(17)

In Equation (17), R_ab,cd represents the PCC between variables a and b when variables c and d are fixed; R_ab,d, R_ac,d, and R_bc,d represent the PCC between variables a and b, variables a and c, and variables b and c, respectively, when the variable d is fixed. The higher the PCC, the greater the influence of the variable on GPP and NPP. A t-test was used to test the significance of the PPC.

2.6.2. Path Analysis of Climate Change to Carbon Fluxes

In order to analyze the direct and indirect effects of climate factors (temperature and precipitation) on carbon fluxes of bamboo forests, path analysis [76] was used to calculate the direct and indirect path coefficients of temperature and precipitation on the carbon fluxes. The formulas for the calculation are as follows:

$${\mathrm{P}}_{\mathrm{i}\to \mathrm{y}} =\frac{{\mathrm{b}}_{\mathrm{i}}{\mathrm{S}}_{\mathrm{i}}}{{\mathrm{S}}_{\mathrm{y}}}$$

(18)

$${\mathrm{P}}_{\mathrm{j}\to \mathrm{i}\to \mathrm{y}} ={\mathrm{r}}_{\mathrm{ij}}{\mathrm{P}}_{\mathrm{i}\to \mathrm{y}}$$

(19)

where ${\mathrm{P}}_{\mathrm{i}\to \mathrm{y}}$ is the direct path coefficient, b_i is the regression coefficient, S_i is the standard deviation of variable i, S_y is the standard deviation of variable y, ${\mathrm{P}}_{\mathrm{j}\to \mathrm{i}\to \mathrm{y}}$ is the indirect path coefficient of variable j acting on variable y through variable i, and r_ij is the correlation coefficient between variable i and variable j.

3. Results

3.1. BEPS Model Validation

In this study, the observed carbon flux data from the Anji site and the Lin’an site during 2011–2014 were used to validate the BEPS model. The daily scale carbon flux data were obtained by accumulating the observed 30-min carbon flux data. The evaluation results are shown in Figure 4, where the R², RMSE, and aBIAS were 0.58, 1.43 g C m⁻² day⁻¹, and 1.21 g C m⁻² day⁻¹, respectively. There was a good correlation between the simulated and observed values of GPP. Therefore, the BEPS model could be considered to be suitable to simulate the productivity of bamboo forests in China.

3.2. Spatiotemporal Evolution of Carbon Fluxes from Bamboo Forests in China

3.2.1. Temporal Evolution Trend

The variation trends of the monthly and annual average GPP and NPP during 2001–2018 are shown in Figure 5. The GPP and NPP exhibited similar temporal variation characteristics. At the monthly scale, the average values of the GPP and NPP of bamboo forests showed unimodal changes. At the annual scale, the average values of GPP and NPP were 904.02 g C m⁻² yr⁻¹ and 716.88 g C m⁻² yr⁻¹, respectively, and the ranges in variation were 764.42–994.61 g C m⁻² yr⁻¹ and 600.03–788.25 g C m⁻² yr⁻¹, respectively. The annual average values of GPP and NPP were the lowest in 2003 and the highest in 2007. During the statistical period, the overall variation trends of GPP and NPP were similar, showing an increasing trend, and the increasing trend was not significant (p > 0.05); the growth rates were 5.20 g C m⁻² yr⁻¹ and 3.88 g C m⁻² yr⁻¹, respectively.

3.2.2. Spatial Distribution Characteristics

The spatial distribution of the mean GPP and NPP values of bamboo forests in China is shown in Figure 6. From Figure 6 we can see that the mean GPP and NPP values had strong spatial heterogeneity. On the whole, the GPP and NPP present a distribution characteristic of being more in the south and east, and less in the north and west. In addition, the spatial distributions of GPP and NPP were compared and it was found that the high GPP and NPP values of bamboo forests were mainly concentrated in northwestern Zhejiang, central Fujian, western Jiangxi, and so on, and the proportion of high-value distribution was gradually increasing. The low GPP and NPP values were mainly distributed in Guizhou, Shanxi, Yunnan, and other regions where the distribution of bamboo forests is relatively scattered.

3.2.3. Analysis of the Fluctuation in Carbon Fluxes

To analyze the fluctuations in carbon fluxes from bamboo forests in China from 2001 to 2018, we calculated the CVs of GPP and NPP and prepared a spatial distribution diagram of the fluctuations based on the CV classification results, as shown in Figure 7. The zones with low fluctuations in GPP and NPP accounted for the largest proportion, at 43.71% and 42.37%, respectively. These were followed by zones with lower fluctuation, at 36.37% and 37.77%, respectively. The areas of higher fluctuation and high fluctuation were considerably small and scattered, among which the areas of higher fluctuation accounted for 3.78% and 3.81%, respectively, and those of high fluctuation accounted for 1.36% for both. In addition, by comparing Figure 7a,b, it was found that the spatial fluctuations in GPP and NPP exhibited evident consistency, where GPP and NPP had low fluctuations, and vice versa.

3.2.4. Analysis of the Trend Slope of Carbon Fluxes

The spatial distribution of the trend slope and the significance of GPP and NPP from 2001 to 2018 are shown in Figure 8.

Figure 8a,c show the spatial distribution of the trend slope and the significance of GPP, respectively. From Figure 8a,c, it can be seen that GPP exhibits an increasing trend (slope_gpp > 0) at 57.58% and a significant increasing trend (slope_gpp > 0, p < 0.01) at 30.32%, mainly distributed in northwestern Zhejiang, western Jiangxi, central Fujian, southwest Anhui, and central Sichuan. GPP exhibits a decreasing trend (slope_gpp < 0) at 42.42% and a significant decreasing trend (slope_gpp < 0, p < 0.01) at 20.53%, mainly distributed in southwestern Zhejiang, eastern Jiangxi, eastern Anhui, and western Guangdong.

Figure 8b,d show the spatial distribution of the trend slope and the significance of NPP, respectively. As shown in Figure 8b,d, NPP shows an increasing trend (slope_npp > 0) at 57.56% and a significant increasing trend (slope_npp > 0, p < 0.01) at 30.32%. NPP shows a decreasing trend (slope_npp < 0) at 42.44% and a significant decreasing trend (slope_npp < 0, p < 0.01) at 20.54%. By comparing Figure 8c,d, it can be seen that the regions with significantly increased and decreased NPP are consistent with the regions that had significantly increased and decreased GPP.

In summary, the spatial distribution of the trend slope of GPP and NPP was similar, the spatial distribution range of the increasing trend was larger than the spatial distribution range of the decreasing trend, and the areas of approximately 30% showed a significant increasing trend, indicating that the carbon fluxes of bamboo forests in China had been gradually increasing over the past 20 years.

3.3. Analysis of Climate Drivers of Carbon Fluxes of Spatiotemporal Evolution

3.3.1. Partial Correlation between Carbon Fluxes and Climate Factors

Climatic factors are important environmental factors that affect the growth of bamboo forests. To quantitatively analyze the influence of climatic factors on carbon fluxes of bamboo forests, the PPC of GPP and NPP with the Pre, T_min, and T_max of bamboo forests in China from 2001 to 2018 were calculated. The results are presented in Figure 9.

Figure 9a,d show the spatial distribution of the PCC of GPP and NPP with Pre, respectively. The proportions of the study area with positive correlations of GPP and NPP with Pre were 52.32% and 54.76%, respectively, mainly distributed in central Zhejiang, northwestern Jiangxi, Chongqing, and Sichuan. The proportions with negative correlations were 47.68% and 45.24%, respectively, mainly distributed in southeast Anhui, northwestern Zhejiang, and eastern Guangxi. Overall, Pre was mainly positively correlated with GPP and NPP, that is, the amount of precipitation considerably promoted the growth of bamboo forests. The proportions with significant (p < 0.05) correlations of the PCC of GPP and NPP with Pre were only 5.74% and 5.67%, respectively.

Figure 9b,e show the spatial distribution of the PPC of GPP and NPP with T_min, respectively. The areas where the GPP and NPP were positively correlated with T_min were 44.15% and 43.68%, respectively, and were mainly distributed in northwestern Zhejiang, central Hunan, and Hubei. Meanwhile, in 55.85% and 56.32% of the areas the GPP and NPP, respectively, showed a negative correlation with T_min, and were mainly distributed in central Zhejiang, Guangdong, and Guangxi. By comparing Figure 9a,d, in addition to Figure 9b,e, it was observed that the PPC of GPP and NPP with Pre and T_min had opposite spatial distribution patterns. Where GPP and NPP were positively correlated with Pre the correlation was negative with T_min, and vice versa. The proportions with significant (p < 0.05) correlations of the PCC of GPP and NPP with T_min were only 6.13% and 6.10%, respectively.

Figure 9c,f show the spatial distribution of the PCC of GPP and NPP, respectively, with T_max. GPP and NPP were positively correlated with T_max, accounting for 49.18% and 48.72%, respectively, mainly distributed in Guangdong, Guizhou, and western Jiangxi, and they negatively correlated with T_max, accounting for 50.82% and 51.28%, respectively, mainly distributed in Guangxi, Anhui, Yunnan, and western Hunan. Overall, there was a nonsignificant negative correlation of GPP and NPP with T_max, which indicated that high temperature somewhat affected the growth of bamboo forests. The proportions with significant (p < 0.05) correlations of the PCC of GPP and NPP with T_max were only 5.13% and 5.20%, respectively.

In summary, a certain correlation existed for the GPP and NPP of bamboo forests in China with precipitation and temperature, and, overall, they were positively correlated with Pre, negatively correlated with T_min, and had an insignificant negative correlation with T_max. In addition, there were evident spatial differences in the correlation of GPP and NPP with climatic factors, and the PPC with Pre and T_min exhibited complementary characteristics.

3.3.2. The Impact of Climate Factors on Carbon Fluxes on a Monthly Scale

The variations in GPP and NPP with temperature and precipitation on a monthly scale are shown in Figure 10. The values of GPP and NPP exhibited different characteristics owing to the influence of hydrothermal conditions. From February to July, with the temperature and precipitation gradually increasing, bamboo forests entered the growing season; therefore, the values of GPP and NPP showed a rapid increase trend. After August, the decrease in GPP and NPP was caused by the gradual decrease in temperature and precipitation, in addition to the fall of bamboo leaves. In December, January, and February the temperature and precipitation are lower, and the values of GPP and NPP were also decreased to the smallest values of the year. In summary, the values of GPP and NPP are closely related to temperature and precipitation, and good hydrothermal conditions are conducive to the growth of bamboo forests.

To further analyze the impact of temperature and precipitation on carbon fluxes of bamboo forests, we conducted a path analysis of the impact of temperature and precipitation on GPP and NPP on a monthly scale. The results are shown in

Table 4. Path analysis of temperature and precipitation to GPP and NPP.

	Climate Factors	Correlation Coefficient	Direct Path Coefficient	Indirect Path Coefficient		Partial Correlation Coefficient
				→ Temperature	→ Precipitation
GPP	Temperature	0.649 **	0.387 **	-	0.26	0.399
	Precipitation	0.659 **	0.415 **	0.24	-	0.423
NPP	Temperature	0.562 **	0.301 **	-	0.15	0.293
	Precipitation	0.602 **	0.412 **	0.11	-	0.386

Note: **, p < 0.01.

. It can be seen that temperature and precipitation have a significant correlation with GPP and NPP. According to the correlation coefficient and partial correlation coefficient, the influence of precipitation on GPP and NPP is higher than that of the influence of temperature on them. In addition, according to the direct path coefficient and indirect path coefficient, the direct influence of precipitation on GPP and NPP is higher than that of the direct influence of temperature on GPP and NPP.

4. Discussion

The simulated value of the carbon fluxes of bamboo forests had a good correlation with the observed value of the flux observation station (Figure 4), and the R² was 0.58. Other than that, in order to further prove the reliability of this study we compared the simulated NPP with a related study (see

Table 5. Comparison of the simulated NPP results in this study with the simulated results of other studies.

Site	Model	Mean NPP (g C m−2 y−1)	Reference
China	BEPS	716.88	This study
Tianmu Mountain, Zhejiang	CASA	740	[77]
Anji, Zhejiang	Triplex-Flux	835.58	[45]
Fujian	BEPS	788.6	[78]

). It can be seen from Table 5 that our simulated mean value of NPP was slightly lower than that of related studies. Due to the fact that structures, mechanisms, and input parameters varied for different models, there are variances in the simulation results of different models. Additionally, there may be differences due to different study areas and periods. Of course, this study also has some shortcomings, the following aspects of which can be analyzed. Firstly, the simulated results of the BEPS model largely depend on the quality of the input data; deficiencies in the input data will affect the accuracy of the simulation results. The resolution of the data in this study is low, so there may be limitations in simulating the carbon flux of bamboo forests in China. Secondly, in this study, the bamboo forest abundance data were used to drive the BEPS model, which solved the influence of mixed pixels on the carbon flux simulation to some extent. However, the phenomenon of “different objects with same spectrums” in remote sensing images will affect the result of bamboo forest extraction. Thirdly, we only used the observed data of two carbon flux observation stations to verify the simulated results of bamboo forest carbon fluxes in China. Therefore, there are limitations on the spatial scale. Finally, the carbon fluxes of bamboo forests were not only affected by climate factors but also by human activities and geographic factors (such as slope, aspect, and elevation). This study only considered the impact of climate factors, so there may still be a certain gap between the simulated results and the real situation.

Due to the rapid growth of bamboo forests and their high ecological, economic, and social value, some areas promoted the reclamation of wasteland and the plantation of bamboo forests [79], which increased the total area of bamboo forests in China. Therefore, the GPP and NPP of bamboo forests also increased. As shown in Figure 5, certain fluctuations occurred in the annual average GPP and NPP values of bamboo forests, which might be related to climate change. For example, in 2003 there was less precipitation and large-scale drought occurred in the summer (Figure 10), which was not conducive to the growth of bamboo forests, leading to low GPP and NPP values in that year. In 2009 and 2010 the values of GPP and NPP were low, which may be related to the natural large-scale low-temperature, snow, and ice disaster in South China in 2008 [80].

This study found that the spatial distribution range of carbon fluxes of bamboo forests was increasing larger than that which was decreasing. The areas with increases were mainly distributed in northwestern Zhejiang, western Jiangxi, central Fujian, and other regions. The reason for the increase in carbon fluxes may be that under the combined influence of favorable factors (such as a warm climate, abundant precipitation, and sufficient radiation) the growth ability of bamboo forests is relatively strong. In addition, these regions have significantly developed bamboo industries and advanced bamboo forest management techniques; therefore, bamboo forests in these regions have increased rapidly with higher productivity. The areas with decreased GPP and NPP were mainly distributed in southwest and central Zhejiang, eastern Jiangxi, northeastern Fujian, and eastern Guangdong. On the one hand, bamboo forests might be reduced due to urban expansion in some areas. On the other hand, because bamboo forests mostly have a scattered distribution, when they are distributed across a small area the difficulty of bamboo forest extraction would be increased, affecting the simulation results of carbon fluxes from bamboo forests.

Climate change has an important impact on vegetation growth. An evident coupling relationship was observed between vegetation and climatic factors [81]. At present, many scholars have analyzed the effects of climatic factors on the carbon fluxes of different vegetation from different spatiotemporal scales, and have found that there is a correlation between carbon fluxes and climatic factors [82,83]. Bamboo forests have a warm and humid climate and are very sensitive to hydrothermal changes. Related scholars have conducted studies on the impact of climatic factors on the carbon fluxes of bamboo forests. For example, Li et al. [4] analyzed the relationship between the carbon fluxes and climatic factors (temperature and precipitation) of bamboo forests in Zhejiang Province from 2011 to 2015, and found that lower precipitation and higher temperatures may have a negative impact on the carbon fluxes from bamboo forests. Chen et al. [84] used eddy correlation technology to continuously observe the carbon fluxes of bamboo forests in Anji, and found that high temperature and drought caused a significant decrease in the carbon fluxes of bamboo forests. These results are consistent with the results of this study on the driving influence of climatic factors and the carbon fluxes of bamboo forests.

5. Conclusions

This study utilized remote sensing data to drive the BEPS model to simulate the carbon fluxes from bamboo forests in China during 2001–2018, and analyzed the spatiotemporal evolution pattern of carbon fluxes and the response of climatic factors to these changes. Our study showed that the simulated values had a good correlation with the observed values, and the R², RMSE, and aBIAS were 0.58, 1.43 g C m⁻² day⁻¹, and 1.21 g C m⁻² day⁻¹, respectively. It provided a feasible way for the study of bamboo forest carbon cycles on a large spatial scale. In addition, our study also suggested that climate change was a driver that affected the spatiotemporal dynamic evolution of carbon fluxes in bamboo forests, and its driving effect exhibited evident spatial variations. This provided a theoretical basis of bamboo forests to cope with climate change.

However, this study still has some limitations. For example, (1) the low resolution of the data limited the simulation of bamboo forest carbon flux; (2) fewer flux observation sites may lead to certain deficiencies in verifying the model simulation results; and (3) we only considered the impact of climatic factors (temperature and precipitation) on the carbon fluxes of bamboo forests. In the future, these limits can be further improved to better simulate the carbon fluxes of bamboo forests in China.