In this study, an empirical allometric equation for estimating forest biomass based on LAI products and ground-based observations of forest biomass was constructed and used to simulate the forest biomass in the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) from 2000 to 2022. Based on the simulated forest biomass, we analyzed the spatial and temporal changes in the forest biomass in the GBA and the drivers of these changes. The specific research goals of this study include: (1) constructing a forest biomass remote sensing estimation model applicable to the GBA, (2) revealing the spatial and temporal characteristics of forest biomass in the GBA over the past 23 years, (3) quantifying the effects of LUC and FBD changes on the total forest biomass in the GBA, and (4) investigating the key environmental controls of the drivers spatiotemporal FBD variations.

The Guangdong-Hong Kong-Macao (GBA) in the south of China (21° 25'-24° 30' N, 111° 12'- 115° 35' E) is known as one of the four global major bay areas (Fig. 1) and one of the most economically advanced (contributes 12% of China's total GDP) regions in China. The GBA is a highly urbanized yet ecologically diverse region, characterized by coastal and estuarine ecosystems, urban ecosystem, croplands and forest. To ensure long-term resilience, the GBA has committed to sustainable and green development, striving to balance rapid economic growth with robust ecological protection. In 2020, the forest area in the GBA was 26,000 km, accounting for 47% of the GBA's total area (5.59 × 10 km). The GBA has a subtropical climate with an average annual temperature of 22.3 °C and an average annual total rainfall of 1832 mm yr. The entire area is dominated by evergreen broad-leaved forests. The elevation is high in the northwest of the GBA and low in the southeast, with mountains distributed in the north and plains mainly in the estuary of the Pearl River. Over the past decades, the forest ecosystem in the GBA has experienced notable climate change and human perturbations. Accurately quantifying the dynamics of GBA's forest biomass is essential for revealing its impact on ecological and socio-economic significance.

The observed specific leaf area (SLA) data, comprising 4131 records within the GBA, were obtained from the TRY database (Table 1). The inventory data on the carbon density at 42 sites (kg C m) for both AGB and BGB were obtained from Xu et al. (2018). The inventory data of the total forest stock volume and forest area data in each district/county of the GBA in 2021 were obtained from the Guangdong Provincial Bureau of Statistics (https://gdzd.stats.gov.cn/, in Chinese). Shenzhen has not publicly released its inventory data on the forest stock volume and forest area. The inventory data on the forest biomass in each district of Shenzhen in 2013 were obtained from Wen et al. (Table 1).

The LAI images were obtained from the MOD15A2H.v061 product (https://lpdaac.usgs.gov/products/mod15a2hv061/). The original LAI data have a temporal resolution of 8 days and a spatial resolution of 500 m. These LAI data have good accuracy across diverse biomes, with a root-mean-square error (RMSE) ranging from 0.66 to 0.69. We calculated the annual mean LAI images over the GBA using the original 8-day LAI dataset using the Google Earth Engine. The land use data over the GBA at a spatial resolution of 30 m from 2000 to 2022 were extracted from the land use dataset produced by Yang et al.. The precipitation (Prec) and temperature (Temp) data over the GBA at a spatial resolution of 0.01 were extracted from the ChinaMet dataset produced by Hu et al., Zhang et al., and Hu et al. have indicated that ChinaMet dataset has a good accuracy of Prec in China, with a RMSE ranging from 4 to 12. Data of Atmospheric CO concentrations was obtained from the National Oceanic and Atmospheric Administration (NOAA, https://gml.noaa.gov/ccgg/trends/). The digital elevation model (DEM) data over the GBA at a spatial resolution of 30 m was extracted from Shuttle Radar Topography Mission (SRTM). The aspect and slope were derived from DEM using the 3D Analyst tool in ArcGIS Pro. The soil properties, including pH, silt, sand, and clay fractions were obtained from the GSDE dataset produced by Shanguan et al..

The leaf biomass density (LBD) can be calculated from the LAI and SLA using the following equation:

where LAI is the leaf area index, and SLA is the specific leaf area. In this study, the average value (= 163.55 cm·g) of 4131 SLA measurements in the GBA was used in Eq. 1 to calculate the from the LAI (Fig. 2). To estimate the uncertainties of the simulated LBD, we also estimated the LBD based on an SLA equaling the mean plus standard deviation (SD) of the 4131 SLA measurements (= 93.05 cm·g) and another SLA equaling the mean minus SD of the 4131 SLA measurements (= 234.05 cm·g; Fig. 2).

Many observations have shown that plant growth generally follows an allometric growth relationship. This provides the possibility and theoretical basis for estimating the AGB and BGB based on leaf biomass. Here we first extracted the LAI data at 42 sites distributed across the GBA (Fig. 1) from the MODIS LAI data. Then we calculated the LBD at each of the 42 forest sites from the extracted LAI using Eq. 1. Finally, we trained two empirical allometric functions for calculating the AGB density (AGBD) and BGB density (BGBD) from the LBD respectively using ordinary least squares (OLS) regression in Origin2023 (Fig. 3a, b). Note that, we have conducted a leave-one-out cross-validation to evaluate the accuracy and reliability of the fitted empirical allometric functions, and to access the uncertainties in the optimized parameter values of the fitted empirical functions (Fig. 3c, d). Through leave-one-out cross-validation, we obtained 42 parameters values of the empirical allometric function for AGBD (Eq. 2) and BGBD (Eq. 3), respectively. The means of the 42 fitted parameter values of Eqs. 2 and 3 were finally used to investigate the spatiotemporal variations of forest biomass in GBA. The maximum and minimum of the 42 fitted parameter values were used to quantify the potential uncertainties in our simulated forest biomass.

We simulated the forest AGB and BGB in each year from 2000 to 2022 across the GBA based on the MODIS LAI and forest distribution data obtained from using Eqs. 1-3 (Table 1). The simulated forest biomass based on the MODIS LAI product was independently validated using the inventory data (Table 1) on the forest biomass in each district and county of the 11 cities in the GBA. The Guangdong Provincial Bureau of Statistics provided the inventory data on the total forest standing stock (m) and total forest area (hm) for each district and county of the 11 cities in the GBA. We first calculated the area-averaged density of the forest standing stock (DSD, m·hm) in each city in the GBA with these data. Previous studies indicated that there is a strong linear relationship between the area-averaged biomass density (DB, t·hm) and DSD, although the specific linear regression function describing this relationship depends on the specific forest type. In this study, we collected five empirical linear regression functions from previous studies to calculate the DB based on the inventory-based DSD. The five regression functions for calculating the DB can be described as follows:

where i is the index of each of the five functions; α and β are the slope and interception of each linear regression function. DB is the Inventory-based biomass calculated by function i. The specific values of α and β in each of the five regression functions can be found in Table 2.

The LAI-based total forest biomass calculated based on Eqs. 1-3 in each district and county of the 11 cities in the GBA was evaluated against the average result of the inventory-based biomass obtained from the Guangdong Provincial Bureau of Statistics (Table 1). The evaluation result shows that the LAI-based forest biomass can be compared well with the inventory-based biomass (Fig. 4). The regression slope was close to 1.0 (0.99), and the determination coefficient (R) of the regression function equaled 0.96.

We further investigated the influences of the FBD (t·hm) and LUC on the forest biomass changes across the GBA from 2000 to 2022. Firstly, the actual total forest biomass in the GBA for a specific year was calculated as follows:

where i is a specific year during 2000-2022; Biomass is the total forest biomass in year i (t); A is the forest area in year i (hm); and FBD is the FBD in year i (t·hm).

We launched a new simulation to investigate the separate effect of the LUC on the forest biomass dynamics in the GBA, in which the forest area A changed year by year, while the FBD was fixed to the state in the year 2000:

We further investigated the key environmental controls of spatial and interannual variations in FBD by partial correlation analysis. Partial correlation analysis is a statistical technique that measures the degree of association between two variables while controlling the influence of one or more additional variables. This method is particularly valuable in ecological studies when environmental factors are interrelated. To investigate the key environmental drivers of the interannual variation in FBD, we calculated the partial correlation coefficient between FBD in each city of GBA and the corresponding mean annual atmospheric CO concentration, temperature, solar shortwave radiation downward, and annual total precipitation. To explore the key environmental controls of the spatial variation in FBD, we calculated the partial correlation coefficient between mean FBD in each 0.25° × 0.25° km grid cell across the GBA and climatic (mean annual temperature and annual total precipitation), soil (pH and texture), and topographic (elevation, slope, and aspect).