Vegetation Dynamics under Rapid Urbanization in the Guangdong–Hong Kong–Macao Greater Bay Area Urban Agglomeration during the Past Two Decades

: Detection of long-term vegetation dynamics is important for identifying vegetation improvement and degradation, especially for rapidly urbanizing regions with intensive land cover conversions. The Guangdong–Hong Kong–Macao Greater Bay Area (GBA) urban agglomeration has experienced rapid urbanization during the past decades with profound impacts on vegetation, so there is an urgent need to evaluate vegetation dynamics across land use/cover change (LUCC). Based on the normalized difference vegetation index (NDVI) during 2001–2020, we used coefﬁcient of variation, Theil–Sen median trend analysis, and Hurst exponent to analyze the spatiotemporal change and future consistency of vegetation growth among the main LUCC in the GBA. Results demonstrated that low NDVI values with high ﬂuctuations were mainly distributed in the central urban areas, whereas high NDVI values with low ﬂuctuations were primarily located in the peripheral hilly mountains. The area-averaged NDVI showed an overall increasing trend at a rate of 0.0030 year − 1 , and areas with vegetation improvement (82.99%) were more than four times those with vegetation degradation (17.01%). The persistent forest and grassland and the regions converted from built-up to vegetation displayed the most obvious greening; NDVI in over 90% of these areas showed an increasing trend. In contrast, vegetation browning occurred in more than 60% of the regions converted from vegetation to built-up. Future vegetation change in most areas (91.37%) will continue the existing trends, and 80.06% of the GBA was predicted to develop in a benign direction, compared to 19.94% in a malignant direction. Our results contribute to in-depth understanding of vegetation dynamics during rapid urbanization in the GBA, which is crucial for vegetation conservation and land-use optimization.


Introduction
Vegetation, an important medium of material circulation and energy flow in terrestrial ecosystems, links the atmosphere, hydrosphere, and pedosphere. Vegetation change greatly affects ecosystem functions and ecosystem services via influencing carbon and nitrogen cycle, water cycle, and climate regulation [1][2][3][4][5][6]. Therefore, in-depth understanding of vegetation dynamics is essential to maintain ecosystem functions and ecosystem services, which is of critical significance for improving environmental quality and promoting human welfare.
Long-term records of remotely sensed spectral vegetation indices derived from sensors such as Advanced Very High Resolution Radiometer (AVHRR), Systeme Probatoire The GBA (21 • (Figure 1). Subtropical monsoon climate dominates the study area with the mean annual temperature, precipitation, and sunshine duration ranges of 20.74-22.67 • C, 1349-2753 mm, and 1484-2294 h, respectively, over the past 60 years. Rapid urbanization has been happening in the GBA for the last several decades and more than 12% of the national GDP was created with less than 1% of the national territory and 5% of the national population in 2018. During the urbanization process, land use and cover has changed dramatically ( Figure 1), and the extensive occupation of ecological land by built-up land tends to impose a negative effect on vegetation and ecosystem services [36].
Remote Sens. 2022, 14, x FOR PEER REVIEW 3 of 18 coefficient of variation, linear regression analysis, the Theil-Sen median trend analysis with a Mann-Kendall significance test, and Hurst exponent based on R/S analysis were used to (1) investigate spatiotemporal fluctuations, trends, and sustainability of vegetation and (2) compare the differences in vegetation dynamics over main LUCC types. The following hypotheses were put forward for this study: (1) vegetation dynamics in the GBA vary across LUCC and (2) urbanization has large and complex impacts on vegetation dynamics. Results of this study are expected to gain insight into the spatiotemporal changes in vegetation in the GBA over the past 20 years, which is beneficial for evaluation of vegetation improvement and degradation under rapid urbanization, thereby providing guidance for future vegetation conservation and restoration.

Study Area
The GBA (21°34′-24°24′N, 111°21′-115°25′E) is located at the Pearl River Delta region in South China with an area of ~56,000 km 2 . The GBA is one of the largest urban agglomerations in the world and consists of 9 mainland cities and 2 special administrative regions: Guangzhou (GZ), Shenzhen (SZ), Zhuhai (ZH), Foshan (FS), Huizhou (HZ), Dongguan (DG), Zhongshan (ZS), Jiangmen (JM), Zhaoqing (ZQ), Hong Kong (HK), and Macao (MC) (Figure 1). Subtropical monsoon climate dominates the study area with the mean annual temperature, precipitation, and sunshine duration ranges of 20.74-22.67 °C, 1349-2753 mm, and 1484-2294 h, respectively, over the past 60 years. Rapid urbanization has been happening in the GBA for the last several decades and more than 12% of the national GDP was created with less than 1% of the national territory and 5% of the national population in 2018. During the urbanization process, land use and cover has changed dramatically (Figure 1), and the extensive occupation of ecological land by built-up land tends to impose a negative effect on vegetation and ecosystem services [36].

NDVI Data
The 16-day MODIS NDVI data (MOD13Q1) at 250 m spatial resolution for 2001-2020 were acquired from the Land Processes Distributed Active Archive Center (LP DAAC) within the U.S. National Aeronautics and Space Administration (NASA) Earth Observing System Data and Information System (EOSDIS) (https://lpdaac.usgs.gov, (accessed on 16 August 2021)). In order to reduce contamination caused by clouds and haze, an adaptive

NDVI Data
The 16-day MODIS NDVI data (MOD13Q1) at 250 m spatial resolution for 2001-2020 were acquired from the Land Processes Distributed Active Archive Center (LP DAAC) within the U.S. National Aeronautics and Space Administration (NASA) Earth Observing System Data and Information System (EOSDIS) (https://lpdaac.usgs.gov, (accessed on 16 August 2021)). In order to reduce contamination caused by clouds and haze, an adaptive Savitzky-Golay filter method was applied to smooth the annual NDVI cycle and reconstruct a high-quality NDVI time-series via the TIMESAT software package [37]. The maximum value composite (MVC) method was then used to yield the annual maximum NDVI dataset for vegetation dynamic analyses. The land use/land cover maps for 2000 and 2020 were obtained from the Resource and Environment Science and Data Center, Chinese Academy of Sciences (http://www.resdc.cn, (accessed on 6 December 2021)). Derived from Landsat images at a pixel resolution of 30 m, land use types were classified into forest, grassland, cropland, waterbody, built-up, and unutilized land. The spatial resolution of land use data was resampled to match NDVI data on ArcGIS 10.6 for subsequent analyses.

Night Light Data
Night light (NL), a useful and comprehensive proxy to monitor energy consumption, and socioeconomic and demographic activities, is widely used to investigate long-term urbanization processes. The night light dataset over the period from 2001 to 2020 was provided by National Tibetan Plateau Data Center (http://data.tpdc.ac.cn, (accessed on 29 March 2022)). A Night-Time Light convolutional Long Short-Term Memory (NTLSTM) network was applied to produce the annual Prolonged Artificial Nighttime-light DAtaset (PANDA), and the high confidential level of the data quality was guaranteed with the Root Mean Squared-Error (RMSE) at 0.73, the coefficient of determination (R 2 ) at 0.95, and the linear slope at 0.99 between modeled and original images [38]. Compared with previous night light products, the PANDA correlated better with socioeconomic indicators such as Gross Domestic Product (GDP), population, and built-up areas. The spatial resolution of PANDA was approximately 1 km, which was resampled to match NDVI data on ArcGIS 10.6 for subsequent analyses.

Coefficient of Variation
The coefficient of variation refers to the ratio of the standard deviation to the mean NDVI time series, which is usually used to characterize vegetation fluctuations and stability. The coefficient of variation in this study was calculated as follows: where CV is the coefficient of variation of NDVI time series; σ NDVI and µ NDVI are the standard deviation and the mean value of NDVI during the study period, respectively; n is the number of years in the study period; NDVI i is the NDVI value of year i.

Trend Analysis
At the regional level, the area-averaged NDVI of each year for the whole study area was calculated and the linear regression analysis was performed to detect the overall vegetation change trend with the slope coefficient, which can be calculated as below: where slope is the inter-annual trend of NDVI; n is the number of years in the study period; and NDVI t is the averaged NDVI value of year t. At the pixel level, the Theil-Sen median trend analysis with a Mann-Kendall significance test was employed to reveal the vegetation change trend of each pixel. Based on non-parametric techniques, this method does not require the samples to follow certain distributions and it is insensitive to outliers, thereby showing robustness in the linear trend analysis of long time series data [39][40][41]. First, the Theil-Sen estimator was used to estimate the slope of the Theil-Sen median: where S NDVI is the median from the set of slopes of the NDVI time series; NDVI i and NDVI j are the NDVI values of year i and year j, respectively. S NDVI > 0 denotes an increasing trend of NDVI, and S NDVI > 0 indicates a decreasing trend of NDVI during the period of 2001-2020. Second, the Mann-Kendall test was used to check the significance of the trend: where n is the length of years in the study area; sgn is a sign function. Z value, calculated using S statistics, was used to test the significance of the NDVI trend. For a given significance level α, |Z| > Z 1−α/2 means that the change trend of the time series is significant. This study selected the significance level α = 0.05, then the trend of vegetation dynamics could be classified into four grades: (1) significant decrease (S NDVI < 0 and |Z| > 1.96), (2) slight decrease (S NDVI < 0 and |Z| ≤ 1.96), (3) slight increase (S NDVI > 0 and |Z| ≤ 1.96), and (4) significant increase (S NDVI > 0 and |Z| > 1.96).

Hurst Exponent
The Hurst exponent based on R/S analysis is generally used to analyze the long-term dependence and predict the future consistency of time series, which has been widely applied in hydrology, climatology, economics, and vegetation dynamics [7,[42][43][44]. The calculation of the Hurst exponent normally includes six steps: (1) Assume there is a time series {ξ(t)}, t = 1, 2, . . . , n.
(2) Define the mean sequence of the time series.
where c is a proportional parameter and H value is the Hurst exponent. The Hurst exponent analysis was adopted in this study to investigate the consistency of NDVI time series. The consistency of the NDVI trend was estimated by the H value ranging from 0 to 1. If 0 < H < 0.5, anti-consistency was predicted for the NDVI time series, indicating a contrary vegetation trend in the future with the H value closer to 0 for more anti-consistency. If H = 0.5, the NDVI time series was stochastic without consistency. If 0.5 < H < 1, consistency was predicted for the NDVI time series, indicating a continuous vegetation trend in the future with an H value closer to 1 for more consistency.

Vegetation Dynamics of LUCC Types
A land use/land cover transfer matrix was used to detect the land transformation information between 2000 and 2020 [45]. During the past 20 years, most of the land use/land cover in the GBA remained unchanged and the area with land use/cover conversions accounted for 7.45% of the whole area of the GBA ( Figure 1). Among the changed land, 85.68% was from vegetated land (cropland, forest, and grassland) to built-up land (Table 1). In addition, there was also a small amount of land transformed from built-up to vegetation. Since the conversion between vegetation and non-vegetation had great impacts on NDVI, this study focused on the vegetation dynamics for the main LUCC types: forest persistence (FP), grassland persistence (GP), cropland persistence (CP), built-up persistence (BP), vegetation to built-up (VB), and built-up to vegetation (BV). The fluctuations, trends, and future consistency of NDVI time series were then analyzed and compared for these LUCC types.

Correlation Analysis
In order to disentangle the relationships between vegetation dynamics and socioeconomic development, Pearson correlation analysis was performed via software Origin 2018 and the correlation coefficient (R) between the annual NDVI and the annual NL for different LUCC types was calculated based on the linear regression as follows: where R is the correlation coefficient between NDVI and NL, NL i and NDVI i are the areaaveraged NL and NDVI values of year i, NL and NDV I are the average NL and NDVI for the whole study period, and n is the number of years.

Spatial Heterogeneity and Fluctuations in Vegetation Dynamics
The NDVI averaged over the period 2001-2020 displayed a distinct spatial pattern: low values mainly distributed in the central region and high values in the peripheral region ( Figure 2a). The NDVI was classified into five grades by an equal-interval classification method, and the area percentage was 0.53%, 10.09%, 15.26%, 39.92%, and 34.20% for each NDVI grade from low to high values. Regions with low NDVI values (NDVI < 0.4) were primarily distributed in the highly urbanized Pearl River Estuary including the southern GZ, FS, ZS, ZH, DG, SZ, and MC, while regions with high NDVI values (NDVI ≥ 0.6) were widely located in ZQ, HZ, JM, HK, and the northern GZ. Area statistics indicated that NDVI differentiated greatly across LUCC types ( where R is the correlation coefficient between NDVI and NL, NLi and NDVIi are the areaaveraged NL and NDVI values of year i, and are the average NL and NDVI for the whole study period, and n is the number of years.

Spatial Heterogeneity and Fluctuations in Vegetation Dynamics
The NDVI averaged over the period 2001-2020 displayed a distinct spatial pattern: low values mainly distributed in the central region and high values in the peripheral region ( Figure 2a). The NDVI was classified into five grades by an equal-interval classification method, and the area percentage was 0.53%, 10.09%, 15.26%, 39.92%, and 34.20% for each NDVI grade from low to high values. Regions with low NDVI values (NDVI < 0.4) were primarily distributed in the highly urbanized Pearl River Estuary including the southern GZ, FS, ZS, ZH, DG, SZ, and MC, while regions with high NDVI values (NDVI ≥ 0.6) were widely located in ZQ, HZ, JM, HK, and the northern GZ. Area statistics indicated that NDVI differentiated greatly across LUCC types ( Table 2). High NDVI values (NDVI ≥ 0.6) accounted for 97.47%, 91.23%, and 71.82% of FP, GP, and CP regions, respectively. Overall, NDVI for most vegetated areas was higher than 0.6, whereas the NDVI decreased significantly when vegetation changed to built-up with 88.18% of the VB region having NDVI lower than 0.6. In comparison, nearly 60% of the BP region had low NDVI values (NDVI < 0.4), but this percentage sharply decreased if built-up converted to vegetation.    (Table 2). Forest showed the highest stability with nearly half of the FP region having the lowest fluctuations, followed by GP and CP regions. Most of the BP region experienced moderate to high fluctuations. In addition, changed regions between vegetation and built-up had relatively higher vegetation fluctuations than the unchanged regions, especially for the VB region with areas undergoing the highest fluctuations accounting for 36.62%.

Trends in Vegetation Dynamics
At the regional scale, linear regression analysis showed that the area-averaged NDVI generally presented a significant increasing trend (p < 0.001, n = 20) at a rate of 0.0030 year −1 during 2001-2020 (Figure 3a). The minimum value of the area-averaged NDVI was 0.6550 and the maximum value was 0.7117, which appeared in 2004 and 2015, respectively. At the pixel scale, the spatial pattern of NDVI trend over the past 20 years illustrated that NDVI increased in 82.99% of the study area, against 17.01% for the NDVI decreased area (Figure 3b). NDVI significantly increased over 60% of the study area, which was mainly distributed in ZC, HZ, JM, SZ, HK, and northern GZ. In contrast, areas with a significant decreasing trend were chiefly in FS, DG, SZ, ZH, and southern GZ, only accounting for 6.35% of the study area. Additionally, areas without vegetation change hardly existed (less than 0.01%).
Vegetation trends among main LUCC types demonstrated great differences (Figure 3c). Forest and grassland exhibited the largest area of significant increasing NDVI, covering 76.64% and 73.11% of FP and GP regions, much higher than other LUCC types. Areas with increasing NDVI could reach to approximately 70% for CP and BP regions, in which the significant increasing areas accounted for 45.38% and 38.09%, respectively. Instead, more than 60% (30.94% for significant decrease and 30.24% for slight decrease) of the VB region experienced an NDVI decrease, while an NDVI increase was observed over 96% (68.97% for significant increase and 27.59% for slight increase) of the BV region, even higher than FP and GP regions.

Future Changes in Vegetation Dynamics
The Hurst exponent analysis indicated that the future dynamics of NDVI in most of the study area would continue to keep the same change trend as the past 20 years (Table 3, Figure 4a). Areas with a consistent NDVI trend (0.5 < H < 1) occupied 91.37% of the region, compared to areas with the anti-consistent trend (0 < H < 0.5) accounting for merely 8.63%. On the basis of the trend analysis and the Hurst exponent, the future development of NDVI was summarized as eight states, which could be concluded as two main directions: benign and malignant ( Table 3). The consistent significant increasing trend was the most important state of future NDVI change, accounting for 59.28% of the study area, followed by the consistent slight increasing trend (17.93%). Together with the small area of anti-consistent decreasing trends, the percentage of areas where NDVI was predicted to develop in the benign direction as a whole reached 80.06%, over four times that in the malignant direction (19.94%). NDVI in most of the areas with malignant development presented a continuous decreasing trend, including the slight decrease (7.91%) and the significant decrease (6.25%).
Moreover, areas that were predicted to convert from an increasing trend to a decreasing trend covered 5.78% of the study area.

Future Changes in Vegetation Dynamics
The Hurst exponent analysis indicated that the future dynamics of NDVI in most of the study area would continue to keep the same change trend as the past 20 years (Table  3, Figure 4a). Areas with a consistent NDVI trend (0.5 < H < 1) occupied 91.37% of the region, compared to areas with the anti-consistent trend (0 < H < 0.5) accounting for merely 8.63%. On the basis of the trend analysis and the Hurst exponent, the future development  to develop in the benign direction as a whole reached 80.06%, over four times that in the malignant direction (19.94%). NDVI in most of the areas with malignant development presented a continuous decreasing trend, including the slight decrease (7.91%) and the significant decrease (6.25%). Moreover, areas that were predicted to convert from an increasing trend to a decreasing trend covered 5.78% of the study area.  Future changes in NDVI among LUCC varied dramatically (Figure 4b). Except for the VB region, NDVI in most of the other regions was predicted to sustain a continuously increasing trend. Areas with a future persistent significant increasing trend accounted for 75.62% and 72.39% of FP and GP regions, respectively, followed closely by the BV region with an area proportion of 68.97%, much higher than that of other LUCC types. In addition, NDVI in 44.88% of the CP region and 37.84% of the BP region would continue to increase dramatically, while this percentage was only 14.66% for the VB region. On the Future changes in NDVI among LUCC varied dramatically (Figure 4b). Except for the VB region, NDVI in most of the other regions was predicted to sustain a continuously increasing trend. Areas with a future persistent significant increasing trend accounted for 75.62% and 72.39% of FP and GP regions, respectively, followed closely by the BV region with an area proportion of 68.97%, much higher than that of other LUCC types. In addition, NDVI in 44.88% of the CP region and 37.84% of the BP region would continue to increase dramatically, while this percentage was only 14.66% for the VB region. On the contrary, the area proportion of NDVI with a consistent significant decreasing trend in the VB region was highest at 30.58%, and this percentage was 25.39% for the area with a consistent slight decreasing trend. Collectively, NDVI in approximately 90% of the FP and GP regions and 70% of the CP, BP, and BV regions was predicted to develop in a benign direction, while NDVI in more than half (58.15%) of the VB region would develop in a malignant direction in the future.

Spatiotemporal Variations in Vegetation Dynamics
Vegetation dynamics over the past 20 years (2001-2020) in the GBA generally displayed a distinct spatial pattern: low coverage, high fluctuations, and decreasing trends in the center, and high coverage, low fluctuations, and increasing trends in the peripheral regions (Figures 2 and 3b), in agreement with previous studies [20,46]. This pattern coin-cided with the topography from alluvial plains in the central estuary area to mountains around its periphery. Firstly, topography is one of the most important determinants of vegetation distribution due to its significant effect on hydrological processes [47], local hydrothermal conditions [7], and soil quality [48]. Thus, forests with high NDVI values are usually distributed in high-elevation mountainous zones. Secondly, human activities are closely related to topography [20]. As negative relationships are commonly observed between elevation/slope and human activities, the land use in the GBA has formed a spatial pattern of "built-up-cropland-forest" from low to high elevations and small to large slopes (Figure 1). Urban expansion and infrastructure construction will inevitably cause vegetation degradation, so high fluctuations and decreasing trends of vegetation change mainly happen at lower elevations and smaller slopes. Thus, it can be concluded that urbanization plays a pivotal role in forming the spatiotemporal patterns of vegetation coverage, fluctuations, and trends in the GBA.
This study found that vegetation at the regional scale exhibited a significant linear greening trend, in line with a previous study [20]. However, the increasing rate of NDVI in our results (0.0030 year −1 ) was smaller than that in the previous work (0.0033 year −1 ). This is possibly due to the difference in the study period. Our analyses were conducted over a longer time series (2001-2020) than the previous study (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016), and the browning trend after 2015 would reduce the overall rate of vegetation greening. In addition, our findings also indicated that the majority (91.37%) of the study area would retain the same vegetation trend in the future, no matter the increasing trend or the decreasing trend (Table 3, Figure 4a). This high proportion of future consistency implied a relatively stable state of long-term dependence of the NDVI time series.

Vegetation Trends among LUCC Types
Rapid economic development has been lasting for several decades in the GBA urban agglomeration, and human activities were widely regarded as one of the main causes of vegetation change. The correlation coefficient (R) between NDVI and NL was used to explore the influences of socioeconomic development on vegetation growth during the past 20 years ( Figure 5). For the whole GBA, NDVI was significantly correlated with NL (p < 0.01), and the positive slope indicated that socioeconomic development was generally beneficial for vegetation greening (Figure 5a), consistent with previous studies highlighting a positive correlation between vegetation and economic growth (GDP, population) [20]. However, the impacts of socioeconomic development on vegetation varied for different LUCC types. Except for the negative correlation for the VB region and no obvious correlation for the BP region between NDVI and NL, significantly positive relationships were observed for other LUCC types ( Figure 5).
Urbanization was widely regarded as the principal cause of vegetation browning in the core areas of the GBA [26,49,50], which was also concluded by our results that the decreasing trend of NDVI appeared in 61.18% of the areas converted from various types of vegetation to built-up lands (i.e., the VB region) throughout 2001-2020 ( Figure 3). Results of correlation analysis between NDVI and NL also confirmed that socioeconomic activities inhibited vegetation growth in the newly developed built-up areas (Figure 5f). However, it was worth noting that vegetation in most (67.73%) of the unchanged built-up areas (i.e., the BP region) presented a greening trend, contradicting the general knowledge that economic development and population growth result in vegetation degradation [51,52]. One possible explanation is the re-feeding of economy to ecology. Taking the most developed GZ and SZ as an example, vegetation of old urban districts in both cities generally showed a significant increasing trend although the NDVI values in these areas were relatively low, and the greening trend was predicted to continue in the future (Figure 4). These results agree with other studies that also observed urban greening in most built-up areas in GZ and SZ recently [53,54]. Nevertheless, not all cities experienced greening in the BP regions, and some cities even presented vegetation degradation in most of the BP region such as ZQ. This may be related to the developmental status of cities that urban environments in cities with higher levels of economy are more likely to have positive effects on indirect enhancement of vegetation growth [55]. The lack of an obvious correlation between NDVI and NL for the BP region also confirmed this complicated vegetation change with socioeconomic growth (Figure 5e). Thus, those highly developed cities (i.e., GZ and SZ) were more likely to implement more investment in urban greening and green development to meet the growing demand for urban human settlement environments [56,57]. This was also supported by a prior work that highlighted a "win-win" relationship between economic growth and vegetation improvement in all municipalities of the Pearl River Delta, especially in recent years [20]. Another reason is the remarkable urban heat island effect in the Pearl River Delta region, which was previously deduced to enhance vegetation productivity [28]. In addition, the fertilization effect of carbon dioxide and nitrogen deposition also contributes to vegetation greening [58,59], particularly in urban areas with intensive transport emissions. Recent studies concluded the general negative direct impacts (land cover transformation) and widespread positive indirect impacts (changes in environments) of urbanization on vegetation growth [27,55], which can also be supported by our findings that vegetation greening occurred in most unchanged urban areas (i.e., the BP region) and vegetation browning happened in most newly developed urban areas (i.e., the VB region). The phenomenon of vegetation improvement in urban cores and vegetation degradation in newly developed peripheries was also observed in other regions [60,61]. As such, vegetation responses to urbanization are highly complex.
The persistent vegetated areas over the past 20 years in the GBA manifested obvious vegetation greening and high future consistency (Figures 3 and 4). For the FP and GP regions, vegetation improved in over 90% of their ranges, and most of them were statistically significant. This was believed to be highly related to the ecological restoration projects such as the Grain to Green Program and the Natural Forest Conservation Program [62,63]. Artificial afforestation and forest facilitation via these projects and other ecological policies were beneficial for vegetation improvement. The significantly positive relationships between NDVI and NL for the FP and GP regions also revealed that more attention would be paid and more investment would be implemented to these ecological projects to facilitate better vegetation restoration and environment improvement (Figure 5b,c). For example, free grain and cash payments were provided for farmers who returned their farmlands on steep topographic slopes or with low yields to forest and grassland via policies such as the transregional compensation for ecological conservation [64,65]. Furthermore, these ecological projects, implemented since 1998, have formed extensive young forests, and forest regrowth and succession might consistently contribute to vegetation greening in the future [6,66]. For the CP region, most areas also showed a significant or slight increasing trend, accounting for 45.38% and 27.40%, respectively. Agricultural intensifications such as multiple cropping, fertilizer usage, effective irrigation, and agricultural mechanization were regarded as the main reason for cropland greening [6,57,67]. Moreover, cropland reclamation also contributed considerably to vegetation enhancement. In conclusion, vegetation greening was commonly observed in the unchanged vegetated areas (i.e., the FP, GP, and CP regions) during the urbanization process (Figure 5b-d), and the area percentage of significant increase was ranked as FP > GP > CP (Figure 3). This was comparable to a previous study that found that the contribution rate to vegetation productivity improvement was highest in the unchanged forest, followed by unchanged grassland and farmland in the six provinces along the southeast coast of China (including the GBA) during 1982-2015 [30]. Environmental efficiency of urban agglomeration development could be evaluated by the slope of linear regression [51], and the steeper slope between NDVI and NL for the FP region indicated a lower cost of forest per unit of socioeconomic growth than the GP and CP regions. Nevertheless, vegetation browning was still detected in some regions and more attention should be paid to these regions to prevent further degradation. To sum up, vegetation dynamics in the GBA vary greatly among different LUCC types during urbanization due to diverse natural and human-induced causes.

Limitations and Uncertainties
Some limitations and uncertainties remain in this study. First, the sensitivity to the presence of short-range dependence may affect the reliability of the Hurst exponent due to the short time of the NDVI observation [16,43]. Second, the vegetation change is more likely to be complex and nonlinear under the influence of climate change and anthropogenic disturbance [68,69], which may have great impacts on the accuracy of the predictions for future consistency of vegetation trends. Moreover, the duration of the consistency of future vegetation trends cannot be estimated by the Hurst analysis [8]. Therefore, a longer period of vegetation dynamics should be deeply studied and the abrupt change in vegetation trends should be necessarily considered in the next work. Additionally, vegetation dynamics are sophisticated and influenced by multiple factors including climate change and human activities, and usually perform distinctively for different seasons [10,50,70], so future efforts are needed to reveal the underlying mechanisms with a finer temporal resolution.

Conclusions
This study evaluated the spatiotemporal heterogeneity, fluctuations, trends, and future consistency of vegetation dynamics by quantifying the coefficient of variation, Theil-Sen median trend analysis, and Hurst exponent of NDVI time series over the period 2001-2020 in the GBA, a rapidly developing urban agglomeration in South China. Results demonstrated that vegetation exhibited a significant greening trend at the regional scale during the past 20 years. Vegetation improvement generally appeared in the peripheral mountainous regions with relatively high NDVI values and low fluctuations, while vegetation degradation mainly happened in the central estuary area with relatively low NDVI values and high fluctuations. Future vegetation change in most areas will continue the existing trends, and 80.06% of the study area was predicted to develop in a benign direction, against 19.94% in a malignant direction. Vegetation dynamics among LUCC types varied remarkedly. The persistent forest and grassland contributed the most to vegetation greening with NDVI in over 90% of areas showing an increasing trend, as well as the regions converted from built-up to vegetation. Most of the persistent cropland and built-up land also presented vegetation enhancement but with a lower proportion at approximately 70%. Nevertheless, more than 60% of the regions converted from vegetation to built-up experienced vegetation degradation. Correlation analysis revealed that vegetation growth was positively correlated with socioeconomic development except for the persistent built-up areas and newly developed built-up areas. Multiple anthropogenic and natural factors were regarded as main reasons for these vegetation changes.