The GBA is located in South China (21°30′-23°40′N, 112°12′-113°48′E) and covers an area of 56,000 km² . The GBA includes two special administrative regions [i.e., Hong Kong (HK) and Macao (MC)] and nine cities in Guangdong Province [i.e., Guangzhou (GZ), Dongguan (DG), Huizhou (HZ), Shenzhen (SZ), Foshan (FS), Zhongshan (ZS), Zhuhai (ZH), Jiangmen (JM), and Zhaoqing (ZQ)]. The GBA has a subtropical climate with a mean annual temperature of 22.3°C and a mean annual rainfall of 1832 mm. The center of the GBA is the alluvial plain of the Pearl River Delta, surrounded by hills in the west, north, and east (Li and Wen, 2018; Lyu et al., 2018; Yang et al., 2020). The GBA is a highly integrated urbanized region with the economy aggregate almost equal to that of the New York Bay Area (Wang et al., 2021); however, the development levels among cities are still uneven. Among the 11 cities, HK, SZ, and GZ are far ahead of other cities in terms of annual gross domestic product (GDP). In 2019, the GBA had a total population of over 72 million and a GDP of 11.59 trillion Yuan (Liu F et al., 2020). The rapid urbanization and high-intensity land development in the GBA has put tremendous pressure on regional hydrological and ecological systems, making this a key policy concern that required attention.

This study will select the Soil and Water Assessment Tool (SWAT) as the hydrological model, and this model needs the topographic, soil properties, land use, and meteorological data as the input data. The topographic data are from the Shuttle Radar Topography Mission at 90 m × 90 m resolution (https://www2.jpl.nasa.gov/srtm/). The soil data (scale = 1:1,000,000) are from the Harmonized World Soil Database, and there are totally 49 soil types (http://www.fao.org/home/en/). The land use data during 1992–2018 are from the European Space Agency, with a resolution of 30 m (https://www.esa.int/), and the main land use types included are Urban aeras, Cropland, Herbaceous, Tree or shrub cover, Bare areas, and Water bodies. The daily precipitation, maximum and minimum temperature, mean wind speed, solar radiation and relative humidity data during 1979–2018 are from the China Meteorological Assimilation Datasets for the SWAT model (CMADS, http://www.cmads.org/), with the resolution of 1/3°; its format is designed to match the input data format of the SWAT model, and thus, the CMADS data can be directly used without format conversion (Liu et al., 2018; Meng et al., 2019; Song et al., 2020; Wang N et al., 2020; Zhang L et al., 2020; Liu et al., 2021). The daily SF data during 2000–2013 at the Boluo hydrological station (the outlet of the East River Basin, see Figure 1) are provided by the Pearl River Water Resources Commission (Shi et al., 2018; Liu S et al., 2020).

With the SWAT model, we divided the East River Basin into 25 sub-basins and the GBA into 290 sub-basins. Combined with 49 soil types and 6 land use types, taking HRU as the basic simulation unit, the hydrological elements of each sub-basin were modeled and aggregated, and they are converged to the outlet of the whole river basin. Sensitivity analysis and automatic calibration were performed by using the SUFI-2 in the SWAT-CUP, which indicated that 10 parameters listed in Table 2 were the most sensitive.

In this study, we divided the SF data into the calibration period of 2000–2006 and the validation period of 2007–2013 to ensure that the SF data used for calibration and validation contain a certain range of hydrological events (i.e., high, normal, and low flow years) in each period (Wang et al., 2010). Figure 2A shows the comparison of the observed and simulated SF data at the Boluo hydrological station during 2000–2006 (calibration period) and 2007–2013 (validation period). The results show that NSE, R² and PBIAS were 0.91, 0.92, 24.1% for the calibration period, while 0.89, 0.9 and 29.6% for the validation period. In general, R² > 0.85, NSE >0.90 and PBAIS <30% indicate that the SWAT model is suitable for the study area. Meanwhile, the results show that P-factor > 0.7, R-factor < 1.5, which proves that the simulation performance is good. Figures 2B,C show the linear fittings of the simulated SF and observed SF during the calibration and validation periods, respectively. The results show that the R² values were greater than 0.9, indicating that the simulation results of this study were reliable (Meng et al., 2019; Guo et al., 2020; Tan et al., 2020; Wang Q et al., 2020).

Figure 3 shows the land use/cover maps in 1992 and 2018 (Figure 3A) and the area changes of land use/cover from 1992 to 2018 (Figure 3B). Compared with 1992, the area of urban areas in the GBA had a dramatic increase of 256.74% in 2018, while the areas of all other land use/cover types decreased (i.e., −13.21%, −67.84%, −2.06%, −78.49%, and −7.83% for cropland, herbaceous, tree or shrub cover, bare areas, and water bodies, respectively. It is worth noting that the main types of land use/cover in the GBA are cropland and tree or shrub cover, which account for about 70% of the GBA and are mainly distributed in ZQ, JM and HZ. In contrast, urban areas are mainly concentrated in the Pearl River Delta Plain (Figure 3A), and its area increased rapidly from 1992 to 2018 due to the rapid socioeconomic development (Figure 3B). The cities with the top five urbanization expansion areas are GZ, SZ, DG, FS, and ZS.

In addition, the conversion matrix of LUCC indicates that most herbaceous and bare areas have changed to urban areas, and only a small proportion of cropland and water bodies have changed to urban areas (Table 3). As a result, the area of urban areas has presented a substantial growth. It is worth mentioning that only the conversion rates from other land use/cover types to urban areas were relatively larger, and the interconversions of different land use/cover types would have different effects on hydrological variables, leading to impacts on regional hydrological cycle in the GBA (Tu et al., 2015).

To isolate the individual impacts of climate change and LUCC on regional hydrological processes, this study firstly determined the change points in the time series of climate variables, and then determined the baseline period and the interference period. Generally, the contribution of potential evapotranspiration is much smaller than that of precipitation, especially on rainy days; moreover, potential evapotranspiration can be closely related to temperature. Therefore, we only analyzed the precipitation and temperature changes but did not analyze the potential evaporation change in this study. In addition, Xue et al. (2021) demonstrated that radiation could be negligible when simulating streamflow with the SWAT model. Therefore, we did not consider the radiation changes in this study. Figures 4A,B show the temporal variations of precipitation and temperature in the GBA during 1979–2018 by using Mann-Kendall method, while Figures 4C,D show the results of change point test for precipitation and temperature by using sliding t-test method. Change points around 2017 and 1997 were found in the time series of precipitation (Figure 4C) and temperature (Figure 4D) from 1979 to 2018. Combining the change point tests for precipitation and temperature, the middle year 1997 was regarded as the change point from 1979 to 2018. Therefore, in this study, years before 1997 were classified as the baseline period and years after 1997 were classified as the interference period. Then, the temporal trends of precipitation and temperature during these two periods could be obtained, as shown in Figures 4A,B, respectively. Different trends could be observed for these two climatic variables since 1979, and the annual changing rates of precipitation and temperature were not quite large. Precipitation showed an increasing trend in the baseline period but a decreasing trend in the interference period, with the changing rates of 57.28 mm/10a and −39.44 mm/10a, respectively. Temperature showed increasing trends in both baseline and interference periods, and the increasing trend in the interference period (i.e., 0.33°C/10a) was slightly greater than that in the baseline period (i.e., 0.22°C/10a). In addition, the climatic variables in the GBA differed in spatial distribution characteristics during the 40-year period. Figure 5 shows the spatial distribution characteristics of precipitation and temperature in the GBA during 1979–2018. Compared with the baseline period, decreasing precipitation could be identified in the northwestern, northeastern, and central-eastern parts of the GBA, while increasing temperature could be identified in the northeastern part of the GBA.

The SWAT model was used to simulate the changes of three hydrological variables (i.e., SF, SM, and ET) in the GBA during the study period. Since the SF differences in Humen and Modaomen between 1992 and 2018 were relatively large, this study will take these two estuaries as examples to analyze the variation characteristics of hydrological variables. During 1979–1997, the SF and SM in Humen increased with the changing rates of 41.97 m³/10a and 6.70 mm/10a, respectively; the SF and SM in Modaomen increased with the changing rates of 40.35 m³/10a and 9.23 mm/10a, respectively. In contrast, during 1998–2018, the SF and SM of Humen and Modaomen showed decreasing trends, with the changing rates of −47.60 m³/10a and −3.10 mm/10a for the SF and SM in Humen, and −45.18 m³/10a and −1.47 mm/10a for the SF and SM in Modaomen. However, the ET trends were significantly different from those of SF and SM. Both Humen and Modaomen showed decreasing trends from 1979 to 1997 with the changing rates of −11.17 mm/10a and −16.18 mm/10a, respectively. During 1998–2018, Humen showed an increasing trend (15.57 mm/10a), while Modaomen showed a decreasing trend (−2.29 mm/10a). Generally, the changes of the SF, SM, and ET in the GBA might be related to both climate change and LUCC. The decrease of precipitation and the increase of temperature in the GBA could lead to the increase of the ET and then affect the SF. Meanwhile, LUCC would affect not only the ET, but also the SF and SM. Different land use types have different ET rates due to different vegetation cover, leaf area index, root depth and albedo. Moreover, soil properties are one of the most important factors affecting infiltration capacity, and LUCC can indirectly affect infiltration and then the SM through changing the soil properties (Wan and Yang, 2007).

Based on the differences of hydrological variables between S1 and other three scenarios, the contributions of climate change and LUCC to hydrological cycle can be evaluated. Figure 6 shows the spatial distributions of these three hydrological variables (i.e., SF, SM, and ET) under four scenarios. These changes were not only the result of climate change, but also the direct human activities including LUCC.

Known from Figure 7, during 1979–2018, the SF in these eight estuaries all showed decreasing trends due to climate change and LUCC. Except for Hongqimen, which had the contribution of climate change to the SF slightly larger than that of LUCC, the contributions of LUCC to the SF in the other seven estuaries were all greater than those of climate change to the SF. This may be related to the changes of climatic variables in the interference period. Compared with the baseline period, precipitation in the interference period decreased by 39.44 mm/10a, and temperature in the interference period increased by 0.33°C/10a, which could synthetically lead to the decrease in the SF. However, the variation amplitudes of climatic variables were still smaller than that of LUCC, indicating that other land use/cover types were mostly converted to urban areas and the area of urban areas increased by 256.74% from 1992 to 2018 (Figure 3B). Therefore, the overall contribution of climate change to the SF was a little smaller than that of LUCC (Figure 7). For the two estuaries with the relatively large SF differences (i.e., Humen and Modaomen), the impacts of climate change and LUCC on the SF were similar. The relative contributions of climate change and LUCC to the decrease in the SF were −47% and −48.6% in Humen, and −47% and −48.3% in Maodaomen. It is worth noting that the sum of the contributions of climate change and LUCC to the SF was less than 100%, which may be due to the interaction between climate change and LUCC (Jiang et al., 2007; Li et al., 2009).

Due to the absence of the observed ET and SM data, this study only quantitatively analyzed the simulated ET and SM. We calculated the relative contributions of climate change and LUCC to the simulated ET and SM by using Eqs 4, 5, and the relevant results are listed in Table 4. Known from Figure 6, the changes in the SM and ET mainly occurred in the eastern and western parts of the GBA. Limited by the length of the article, this study only showed the detailed analyses of Humen and Modaomen, and other estuaries basically had similar characteristics.

Under the impact of climate change, ET showed increasing trends in both Humen and Modaomen, with the percentage changes of 1.7 and 1.0%, respectively. In contrast, SM showed opposite trends in both Humen and Modaomen, with the percentage changes of −1.5% and −3.8%, respectively. As shown in Table 4, climate change is the main driver of SM changes, with precipitation as the main source of SM (Li, 2014; Luo et al., 2021). Generally, SM increases rapidly after the occurrence of precipitation and then gradually decreases due to plant transpiration and ET (Holsten et al., 2009). Temperature can indirectly affect SM through controlling ET. In the context of global warming, the increase in ET caused by higher temperature is universally accepted (Holsten et al., 2009; Wang et al., 2017). In general, SM can be regarded as the balance between precipitation and ET (He et al., 2012; O’Gorman, 2012). Under the impact of LUCC, ET showed decreasing trends in both Humen and Modaomen, with the percentage changes of −3.6% and −1.5%, respectively. The variation amplitudes were larger than those under the impact of climate change. SM also showed increasing trends in both Humen and Modaomen, with the percentage changes of 10.5 and 2.3%, respectively. Under the combined impacts of climate change and LUCC, ET showed decreasing trends in both Humen and Modaomen, with the percentage changes of −2.0% and −0.7%, respectively, while SM showed increasing trends in both Humen and Modaomen, with the percentage changes of 8.3 and 3.6%, respectively. The changes in ET and SM were closely related to climate change and LUCC in the study area. Influenced by urbanization, part of precipitation is retained by impervious surfaces and directly returned to the atmosphere as ET; part of precipitation infiltrates into the underground, which is absorbed by plant roots and then transpires to the atmosphere. On the one hand, urbanization may lead to a decrease in the amount of ET and infiltration by reducing permeable surfaces, plant transpiration, and canopy retention (Chen et al., 2011; Zhang and Liang, 2018). On the other hand, conversion from herbaceous to cropland may result in an increase in the SM for vegetation growth (Atchley and Maxwell, 2011; Chang et al., 2014).

The SWAT model was to quantify the impacts of climate change and LUCC on regional hydrological processes under four scenarios during 1979–2018. The simulation results were calibrated and validated based on the observed SF data of the East River Basin, which were consistent with those in the previous results (Wu and Chen, 2013). As mentioned before, the overlap area between the GBA and the East River Basin accounts for 20.74% of the GBA and 38% of the East River Basin, which may affect the representativeness of the results to some extent. Moreover, although the SWAT model could well simulate the three hydrologic variables in the GBA, there were still some uncertainties in the model simulations. First, the SWAT-CUP used the SUFI-2 algorithm to evaluate the inherent uncertainties in the SWAT model, which might be due to the uncertainties in model calibration, forcing data, model structure, and model parameter estimation (Rouholahnejad et al., 2012; Ashraf Vaghefi et al., 2014; Abbaspour et al., 2017). Due to the complexity of hydrological processes, errors in hydrological data and imperfections in the model structure can bring uncertainties in the parameter calibration of the SWAT model, which may further affect the accuracy of the simulation results. Second, human activities such as reservoir operation and crop management were not considered in this study, resulting in the differences between the simulated and observed SF (Baker and Miller, 2013; Guo et al., 2016; Wu et al., 2019). Third, the simulation results were only based on two land use/cover maps in 1992 and 2018, which might cause underestimations or overestimations of LUCC in certain years and then the uncertainties in ET and SM (Li et al., 2009). Fourth, this study used the sliding t-test approach to determine the baseline period and the interference period; however, different assumptions and methods may lead to different change points, further affecting the follow-up analysis. Finally, our model was calibrated and validated based on only SF data, and thus, ET and SM simulations need further validation if the observed ET and SM data were available.

The results in this study indicated that the contribution of LUCC to regional hydrological process was greater than that of climate change. Unprecedented urban expansion has occurred in the GBA during the past decades, which may surely increase the area of impervious surface in this region. Land use/cover types in the GBA have mainly changed from cropland, herbaceous, tree or shrub cover, bare areas, water bodies to urban areas, and such results are consistent with those in the previous studies (Yang et al., 2020; Zhang J et al., 2020; Wang et al., 2021).

In addition, the conversions among different land use/cover types can affect the hydrological variables. In general, tree or shrub cover can produce less SF and SM than other land use/cover types, but more ET (Li et al., 2009). Since the biomass per unit area of tree or shrub cover and herbaceous is much lower than that of cropland, the water storage capacity, transpiration capacity, and thermal capacity of tree or shrub cover and herbaceous are all lower than that of cropland (Li Z et al., 2020). Therefore, the conversions from tree or shrub cover and herbaceous to cropland led to the changes in the SF, ET, and SM in this study. Normally, urbanization may lead to an increase in the SF. Although most of other land use/cover types in the GBA were converted to urban areas, the proportion of urbanization area in the GBA was still relatively small (i.e., about 13%). Therefore, the increase in the SF caused by urbanization may be smaller than the decrease in the SF caused by the transformation processes of other land use/cover types, leading to the overall decreasing SF in the GBA. Moreover, with the development of urbanization, a large number of natural surfaces have been transformed into impervious surfaces, which may largely affect the ecological environment of cities and generate a series of ecological and environmental problems, such as urban heat island effect and urban flooding disasters (Shi et al., 2020). Therefore, the proportion of impervious surfaces should be controlled, and green space should be increased in the future development of cities in the GBA.

This study evaluated the temporal and spatial variation characteristics of climate variables from 1979 to 2018, showing that precipitation has decreased and temperature has increased during the past 40 years. The hydrological simulation results showed the decrease in the SF, increase in the ET, and decrease in the SM due to climate change. Such results were consistent with those in the previous studies (Ji and Duan, 2019; Liu Z et al., 2020). Precipitation is an important component of hydrological cycle, which can supply water and control the ET to some extent. During precipitation, the water vapor content near the surface is close to the saturated water vapor, so the evaporation from the soil surface can be ignored. Moreover, according to the Clausius-Clapeyron equation, saturated water vapor increases along with increasing temperature, which increases the water vapor pressure difference between the surface and air, and further enhances the ET (Guo et al., 2016; Li Q et al., 2017; Shrestha et al., 2019; Yang et al., 2017). As a result, the decrease in precipitation and the increase in temperature directly led to the decrease in the SF and SM but increase in the ET (Figure 7 and Table 4). Previous studies have shown that the impacts of climate change on the regional SF mainly originates from the changes in precipitation, while precipitation is an important source of the ET and SM. With the gradually worsening trend of global warming, the effect of high temperature on the regional hydrological cycle aggravates rapidly (Wang et al., 2017; Liu et al., 2019). Such changes to the regional hydrological cycle caused by climate change have caused a redistribution of water resources in time and space, which has important implications for the rational utilization and planning management of water resources.

This study also quantified the combined impacts of climate change and LUCC on regional hydrological cycle. The results indicated that the impact of LUCC on the SF was greater than that of climate change in the GBA. Such results were inconsistent with those in the previous studies on adjacent regions, which may be due to the different change points, lengths of data, and hydrological models (Wu et al., 2017; Mendoza-Ponce et al., 2018; Wang and Stephenson, 2018). For example, reported that the earlier the change point, the smaller the relative contribution of climate change and the larger the relative contribution of human activities.

However, there is an interaction between climate change and LUCC. Climatic variables affect LUCC, and in turn LUCC affects regional climate (Li G et al., 2017). The coupled climatic variables, land use/cover types, and hydrological models are complex and complicated by the interactions and feedbacks between subsystems through physical processes that separate the individual impacts of climate change and LUCC on hydrological processes from their combined impacts. In this study, urbanization-induced LUCC has direct and indirect effects on regional climate. For other land use/cover types, in the vertical direction, tree or shrub cover can reshape the regional surface roughness, increase atmospheric dynamic turbulence, transport water vapor to a higher position, and make the atmospheric water saturation in the clouds higher. However, the reduction of tree or shrub cover will reduce the atmospheric water vapor and absorbed heat, which may further hinder cloud condensation in wet and cool atmosphere (Li Z et al., 2020).

In addition, the spatial distributions of hydrological variables under climate change and LUCC showed that, unlike the SM, there was spatial heterogeneity in the ET changes caused by climate change and LUCC in the GBA (Figure 6). There was a strong regional pattern, e.g., climate change increases the ET in northeastern, northwestern, and southwestern regions, while LUCC decreases the ET in northeastern and southwestern regions.

We need to fully aware of the following limitations, which are mainly related to three aspects. First, our assessment is only on an annual scale, which does not reflect peak SF within a year, and runoff on an interannual scale is averaged to include runoff during both abundant and dry periods. Second, due to limitations of the SWAT model and available data, other human activities (e.g., reservoir operations) were not considered in this study. Reservoirs, such as the Xinfengjiang Reservoir in the East River Basin, would most likely affect the streamflow simulation. For example, storage change in reservoirs may increase the streamflow in dry season but decrease the streamflow in rainy season. Finally, this study ignores the interactions between climate change and LUCC by separating their individual impacts from their combined impacts. Therefore, the impacts of other human activities and the interactions between climate change and LUCC remain unresolved issues. Nevertheless, aware of the above limitations, other findings can be explored for quantifying the impact of climate change and LUCC on regional hydrological processes in the GBA based on this study, which has important implications for sustainable water resources management in developed regions.