Combined impacts of climate change and human activities on blue and green water resources in a high-intensity development watershed

Sustainable management of blue and green water resources is vital for the stability and sustainability of watershed ecosystems. Although there has been extensive attention paid to blue water (BW), which is closely related to human beings, the relevance of green water (GW) to ecosystem security is typically disregarded in water resource evaluations. Specifically, comprehensive studies are scarce on the detection and attribution of variations of blue and green water in the Dongjiang River basin (DRB), an important source of regional water supply in the Guangdong–Hong Kong–Macao Greater Bay Area (GBA) of China. Here we assess the variations of BW and GW scarcity and quantify the impacts of climate change and land use change on BW and GW in DRB using the multi-water-flux calibrated Soil and Water Assessment Tool (SWAT). Results show that BW and green water storage (GWS) in DRB increased slowly at rates of 0.14 and 0.015 mm a⁻¹, respectively, while green water flow (GWF) decreased significantly at a rate of −0.21 mm a⁻¹. The degree of BW and GW scarcity in DRB is low, and the per capita water resources in more than 80 % of DRB exceed 1700 m³ per capita per year. Attribution results show that 88.0 %, 88.5 %, and 39.4 % of changes in BW, GWF, and GWS result from climate change. Both climate change and land use change have decreased BW, while climate change (land use change) has decreased (increased) GWF in DRB. These findings can guide the optimization of the allocation of blue and green water resources between upper and lower reach areas in DRB and further improve the understanding of blue and green water evolution patterns in humid regions.

1 Introduction

Land use and land cover change (LUCC) and climate variability may alter hydrological processes in watersheds (Berezovskaya et al., 2004; Chagas et al., 2022; Konapala et al., 2020; Tan et al., 2022a), which successively affect variations of regional water resources (Hoek van Dijke et al., 2022; Pokhrel et al., 2021; Stocker et al., 2023; Suzuki et al., 2021), potentially leading to ecosystem degradation and severe water shortage crises (Aghakhani Afshar et al., 2018; Zuo et al., 2015). With the development of society and the economy, there is an increasing need for water resources to accommodate human water utilization, encompassing agricultural, domestic, and industrial water usage. Water scarcity and spatiotemporal mismatch between regional water supply and demand in certain regions are becoming increasingly severe, significantly affecting sustainable development in these regions (Cook et al., 2014). Quantifying water resources in a changing environment is crucial for guiding efficient and sustainable water use.

Previous studies on water resource assessment explored the effects of climate change and anthropogenic factors on available water resources, including streamflow (Ahiablame et al., 2017; Tan et al., 2023), baseflow (Ficklin et al., 2016; Tan et al., 2020), lake water (Acero Triana and Ajami, 2022; Tao et al., 2020), and groundwater (Han et al., 2020). Falkenmark and Rockström (2006) introduced a novel perspective on water resource assessment by categorizing water resources as blue water (BW) or green water (GW). BW is the total of deep aquifer recharge and river streamflow, such as water in lakes and rivers. Water users such as industries, agriculture, and municipal users can directly utilize BW. By contrast, GW is the portion of precipitation that is not drained to the river for streamflow generation. GW is temporarily retained in the soil before eventually being released back into the air by evapotranspiration (ET). GW encompasses both green water flow (GWF) and green water storage (GWS) (Veettil and Mishra, 2018; Zang and Liu, 2013). Traditional water resource assessments concentrate on available water resources and only consider BW while neglecting GW (Dai et al., 2022), although GW is also essential. GW supplies about 80 % of the total water resources, sustaining crop growth and the sustainable development of forest and grassland ecosystems in arid regions or during dry seasons (Li et al., 2018; Schuol et al., 2008). Green water scarcity can lead to ecosystem degradation and intensify competition between human needs and ecosystems for water resources (Falkenmark et al., 2003; Veettil and Mishra, 2018). Compared to traditional streamflow assessment methods, water resource scarcity assessment methods based on the framework of BW and GW are more appropriate for maintaining sustainable water resource management (Cooper et al., 2022; Liu et al., 2017). Recently, some studies characterized water scarcity by assessing variations of BW and GW. For example, Veettil and Mishra (2020) assessed blue water scarcity and green water scarcity to show the water security status of counties in the United States. Hoekstra et al. (2012) used the concept of the BW footprint to study water scarcity issues. Schyns et al. (2019) used the GW footprint to investigate green water scarcity and found that the increasingly severe shortage of GW poses a significant threat to natural ecosystems.

The impacts of climate change and human activities on the hydrological cycle processes in watersheds have attracted widespread attention (Ahiablame et al., 2017; Chouchane et al., 2020; Cooper et al., 2022; Tan et al., 2022b; Veettil and Mishra, 2016). Changes in land use alter the underlying surface conditions. For example, afforestation or deforestation may exacerbate or alleviate global or regional climate change and thus affect hydrological cycle processes (Bai et al., 2020; Lian et al., 2020; Qiu et al., 2023). Changes in land use often lead to alterations in land–atmosphere interactions, and vegetation cover changes are essential for regulating climate systems and land ecosystems (Foley et al., 2005; Huang et al., 2020). Large-scale greening could modify geophysical interactions between the atmosphere and the ground, impacting larger or local regional hydrological cycles. Land degradation (Walters and Babbar-Sebens, 2016), deforestation (Lee et al., 2011), and urbanization (Mohan and Kandya, 2015; Zhang et al., 2018) also have far-reaching effects on the climate and hydrological cycle.

Climate change is also crucial to the variations in BW and GW resources. Precipitation is the source of BW and GW, and factors such as temperature, solar radiation, and potential evapotranspiration significantly influence the changes in BW and GW in watersheds, especially in GWF (Pandey et al., 2019; Schewe et al., 2014). For a single watershed, BW depends directly on precipitation and evapotranspiration (GWF) (Shen et al., 2017; Vano et al., 2012). Furthermore, precipitation intensity can have a significant impact on the redistribution of precipitation, BW, and GW by altering infiltration and runoff generation processes (Eekhout et al., 2018; Nearing et al., 2005). Therefore, it is crucial to quantify the effects of climate change and LUCC on BW and GW resources in a watershed for effective water resource planning and management.

Water resource management is the primary issue to be addressed for water security. Hydrological models are important tools for meeting various needs in water resource management. Hydrological model simulation is an effective method for evaluating changes in blue and green water resources. As a widely used semi-distributed parametric hydrological model, the SWAT (Soil and Water Assessment Tool) model is increasingly being used in water resource management at the watershed scale. Based on the SWAT model, researchers simulated the spatiotemporal changes in blue and green water resources in Iran (Ahiablame et al., 2017), the Yangtze River basin (Nie et al., 2023), the Poyang Lake basin (Liu et al., 2023), and India (Sharma et al., 2023). Some studies have also used model simulations to analyze the effects of climate change and human activities on water resource changes in the Meki River basin (Hordofa et al., 2023), China (Liu et al., 2022), and Ningxia (Ahiablame et al., 2017). However, most of the hydrological models used in the study were calibrated and validated using only observed streamflow data without checking the accuracy of other simulated water variables, which can lead to uncertainties in modeling soil moisture and evapotranspiration (Nie et al., 2023).

The Dongjiang River basin (DRB) is a crucial water source region for core cities in the Greater Bay Area (GBA), such as Shenzhen, Hong Kong, and Huizhou. Given the significant BW demand from agriculture, domestic utilization, and industry as well as the GW demand from over 18 000 km² of forested land, the water resource stress in DRB is extremely high, although DRB is located in wet South China (Liu et al., 2018). The growing mismatch between increasing water demand and decreasing water supply, along with seasonal and pollution-induced water scarcity issues, is becoming increasingly prominent (Yang et al., 2018). However, the majority of current studies on water resources of DRB focus on changes and scarcity of surface water and groundwater (BW) while overlooking the critical role and spatiotemporal variations of GW (Huang et al., 2022; Jiang et al., 2023; Wu et al., 2021). With the high-intensity urbanization and climate change in DRB, changes in BW and GW resources in DRB remain unknown.

This study aims to analyze the influence of climate change and LUCC on BW and GW in DRB. The objectives of this study are (a) to build the SWAT model for DRB hydrological simulation, (b) to quantitatively evaluate the spatial and temporal variation of BW and GW in DRB, (c) to assess the status of water scarcity in DRB using the framework of BW and GW resources, and (d) to estimate the effects of climate change and LUCC on BW and GW in DRB.

2 Materials and methods

2.1 Study area

The Dongjiang River is an important tributary of the Pearl River, located between longitude 113°25^′–115°52^′ E and latitude 22°26^′–25°12^′ N. It originates in Xunwu County, Jiangxi Province, flows through Jiangxi and Guangdong provinces, and goes across major cities, including Longchuan, Heyuan, Dongguan, and Shenzhen. The trunk stream of the Dongjiang River has a total length of 562 km. DRB covers a watershed area of 3.5×10⁴ km². DRB is in the subtropical monsoon climate zone with adequate precipitation and high temperatures. The average annual precipitation ranges from 1500 to 2400 mm, and the average temperature of the basin is 21°C (Wu et al., 2019a). The altitude of the basin decreases from the northeast to the southwest. Regions of the upper reaches of DRB are dominated by mountains and hills, those of the middle reaches of DRB are dominated by hills and plains, and those of the lower reaches of DRB are dominated by plains.

Previous hydrological simulation studies of DRB mainly used the Boluo hydrometric station as the outlet of the watershed (He et al., 2013; Wu et al., 2019a), so this study only analyzes the area of DRB where water flows to the Boluo station (Fig. 1). The Boluo hydrometric station is the main control station in the lower reaches of the Dongjiang River. The Boluo hydrometric station occupies a drainage area of 25 325 km², which is 71.7 % of the total area of DRB. Since the 1950s, more than 896 reservoirs, ponds, dams, and other water conservancy facilities have been constructed in DRB. Among them, the Baipenzhu Reservoir, Fengshuiba Reservoir, and Xinfengjiang Reservoir are the three largest reservoirs in the basin, with a cumulative storage capacity of 17 048×10⁶ m³. The Dongjiang–Shenzhen Water Supply Project constructed in 1964 diverts water from the Dongjiang River to Shenzhen and Hong Kong to provide freshwater resources for municipal use. Over 70 % of Hong Kong's freshwater supply comes from the Dongjiang River. Therefore, it is crucial to comprehend the shifts in water resources within DRB for projecting future available water resources for the development of GBA.

Figure 1 Location and characteristics of the study area: (a) location of the watershed, spatial distribution of the hydrometeorological stations, and digital elevation model (Farr et al., 2007). (b) Land use map (Xu et al., 2018).

2.2 Methodology

2.2.1 SWAT model

The SWAT model was adopted to simulate hydrological processes and estimate the amount of BW and GW for DRB (Arnold et al., 1998; Neitsch et al., 2002). The SWAT model is widely applied to simulate streamflow and surface runoff (Arshad et al., 2022; Martínez-Salvador and Conesa-García, 2020; Nie et al., 2023). The SWAT model is also widely utilized for exploring changes in BW and GW (Dai et al., 2022; Liang et al., 2018; Schuol et al., 2008).

In SWAT modeling, DRB was divided into 63 subbasins (Fig. S1 in the Supplement), and each subbasin was then categorized into hydrological response units (HRUs) depending on land use, soils, and slope. The Soil Conservation Service (SCS) curve number method was used for flow partitioning according to land use, soil type and antecedent soil moisture. The Penman–Monteith method was used to calculate potential evapotranspiration, which comprehensively considered various climatic factors such as solar radiation, air temperature, wind speed, and relative humidity (Arnold et al., 1998; Neitsch et al., 2002).

2.2.2 Model calibration and validation

To reduce the influence of hydraulic engineering, the SWAT model was calibrated and validated by utilizing monthly restored natural streamflow at the Boluo and Heyuan hydrometric stations. The optimum model parameters are shown in Table 1. All the selected parameters are automatically calibrated with 500 simulations via SWAT-CUP. The warmup period for model simulations is the first 2 years of the simulation period. Reconstructed natural streamflow in 1970–1979 was used to calibrate the model, and monthly time series of reconstructed natural streamflow, ET from Global Land Evaporation Amsterdam Model (GLEAM) and soil moisture data from ERA5 during 1980–1989 were used to validate the model. The calibration period for this study was 1970–1979, and the validation period was 1980–1989. Three metrics, i.e., the determination coefficient (R²), the percentage bias (PBIAS), and Nash–Sutcliffe efficiency (NSE), were applied to evaluate the simulation performance of the SWAT model:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^n{\left(Q_{\mathrm{nat}}-Q_{\mathrm{sim}}\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^n{\left(Q_{\mathrm{nat}}-Q_{\mathrm{ave}}\right)}^{\mathrm{2}}}},\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{PBIAS}={\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{sim}}}-Q_{\mathrm{ave}}}{Q_{\mathrm{ave}}}}\times \mathrm{100},\end{array}$$

$$\begin{array}{}\text{(3)}& R^{\mathrm{2}}={\left[{\displaystyle \frac{\sum _{i=\mathrm{1}}^n\left(Q_{\mathrm{nat}}-Q_{\mathrm{ave}}\right)\left(Q_{\mathrm{sim}}-\stackrel{\mathrm{‾}}{Q_{\mathrm{sim}}}\right)}{\sqrt{\sum _{i=\mathrm{1}}^n{\left(Q_{\mathrm{nat}}-Q_{\mathrm{ave}}\right)}^{\mathrm{2}}\sum _{i=\mathrm{1}}^n\left(Q_{\mathrm{sim}}-\stackrel{\mathrm{‾}}{Q_{\mathrm{sim}}}\right)}}}\right]}^{\mathrm{2}},\end{array}$$

where Q_nat, Q_ave, Q_sim, and $\stackrel{\mathrm{‾}}{Q_{\mathrm{sim}}}$ are monthly natural streamflow, mean monthly natural streamflow, simulated streamflow, and mean monthly simulated streamflow, respectively. n is the total number of time steps.

Table 1 Range of the main parameters and their optimal values obtained from the model calibration.

Note: the symbols V and R denote a replacement and a relative change to the default parameter value.

This study reconstructed the natural monthly streamflow series of the basin by combining the inflow and outflow of the three major reservoirs (Xinfengjiang Reservoir, Fengshuba Reservoir, and Baipenzhu Reservoir) in DRB, based on the watershed water balance (Tu et al., 2018):

$$\begin{array}{}\text{(4)}& Q_{\mathrm{nat}}=Q_{\mathrm{o}}+\mathrm{\Delta }Q=Q_{\mathrm{o}}+Q_{\mathrm{in}}-Q_{\mathrm{out}},\end{array}$$

where ΔQ is the total reduced water volume and Q_o, Q_in, and Q_out are the observed streamflow, reservoir inflow, and reservoir outflow, respectively.

2.3 Calculation of blue and green water and water security indicators

2.3.1 Calculation of blue and green water

BW is calculated from the sum of water yield (SWAT output WYLD) and groundwater storage. The former refers to the amount of water that leaves the HRU and enters the channel. The latter represents the net amount of water recharged to aquifers (SWAT output GW_RCHG) and the amount of aquifer water discharged to the main channel (SWAT output GW_W) during a time step (Hordofa et al., 2023). GW can be divided into two components, i.e., GWF, which is the actual evapotranspiration (SWAT output ET) from the HRU, and GWS, which is the soil water moisture (SWAT output SW) (Nie et al., 2023; Veettil and Mishra, 2018). The calculation of the Green Water Index (GWI) involves dividing the quantity of GW by the sum of BW and GW (Ding et al., 2024).

2.3.2 Blue and green water scarcity

Blue water scarcity (BWSC) is determined by the quotient of BW withdrawal and availability. The estimation of BW withdrawals (BWWs) in this study involved the multiplication of the aggregate population in each subbasin by the combined water consumption per person (Liang et al., 2020). The population of each subbasin was extracted from the population raster data. Blue water availability (BWA) represents the quantity of water that can be utilized without negatively impacting the river ecosystems. Exhaustive exploitation of BW in rivers may adversely impact river ecosystems. Previous studies have generally used environmental flow requirements (EFRs) as a suitable metric for sustaining robust ecosystems (Honrado et al., 2013). According to the studies of Richter (2010) and Richter et al. (2012), extracting more than 20 % of the water from rivers may result in ecological degradation. Therefore, 20 % of streamflow can be deemed BW and used for water supply (Veettil and Mishra, 2016). The calculations of EFR, BWA, and BWSC are as follows:

$$\begin{array}{}\text{(5)}& {\mathrm{EFR}}_{(a,t)}=\mathrm{0.8}\times Q_{{\mathrm{mean}}_{(a,t)}},\end{array}$$

where EFR_(a,t) is the EFR for subbasin “a” during time “t”; Q_mean is the long-term monthly average streamflow.

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{\mathrm{BWA}}_{(a,t)}=Q_{(a,t)}-{\mathrm{EFQ}}_{(a,t)}\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{BWSC}=\mathrm{BWW}/\mathrm{BWA}\end{array}$$

Green water scarcity (GWSC) is defined as the ratio between green water footprint (GWFO) and green water availability (GWA). GWFO denotes the actual evapotranspiration from the watershed. GWA is the soil moisture that is available for evapotranspiration and vegetation transpiration and is equal to the initial soil moisture (Liang et al., 2020). GWSC can be formulated as

$$\begin{array}{}\text{(8)}& {\mathrm{GWSC}}_{(a,t)}={\mathrm{GWFO}}_{(a,t)}/{\mathrm{GWA}}_{(a,t)},\end{array}$$

where GWSC is the green water scarcity, GWFO_(x,t) is the actual evapotranspiration, and GWA_(a,t) is the initial soil moisture.

Based on the blue water scarcity and green water scarcity, the water scarcity of a region is categorized as mild scarcity, moderate scarcity, severe scarcity, and extreme scarcity, with thresholds set at 100 %, 150 %, and 200 %, respectively.

2.3.3 Regional water stress

The Falkenmark index (FLK) (Falkenmark et al., 1989) is a widely used measure of water stress, defined as the proportion of BWA to the overall population. The Falkenmark index is classified as no stress, stress, scarcity, and absolute scarcity based on per capita water use. Absolute scarcity is regarded as occurring in areas where the indicator threshold is lower than 500 m³ per capita per year, and no stress is thought to occur in areas where the threshold is higher than 1700 m³ per capita per year.

2.4 Calculation of relative contributions

2.4.1 Scenario design and simulation

Three scenarios were constructed to assess the impacts of climate change and LUCC on BW and GW by changing climate conditions (land use) while keeping land use (climate conditions) for the three scenario simulations (Table 2). The land use map was fixed when simulating the influences of climate change on blue and green water (S2–S1), while climate conditions were fixed when simulating the influences of LUCC on blue and green water (S3–S2). The climate conditions and the land use were altered when assessing the joint influences of climate change and LUCC on blue and green water (S3–S1).

Table 2 Scenario settings for the simulation of effects of climate change and LUCC on blue and green water.

2.4.2 Relative contribution rate calculation

The influences of climate change and LUCC on the changes in blue and green water in different periods are evaluated utilizing the relative contribution (RC) in this study (Li et al., 2021).

The climate change contribution to BW and GW change is estimated by

$$\begin{array}{}\text{(9)}& {\mathrm{RC}}_{\mathrm{C}}={\displaystyle \frac{\left|X_{\mathrm{2}}-X_{\mathrm{1}}\right|}{\left|X_{\mathrm{2}}-X_{\mathrm{1}}\right|+\left|X_{\mathrm{3}}-X_{\mathrm{2}}\right|}}\times \mathrm{100}\mathit{\%},\end{array}$$

where X₁, X₂, and X₃ are the amounts of water including BW, GWF, and GWS, respectively, for scenarios S1, S2, and S3.

The contributions of LUCC to changes in BW and GW are estimated by Eq. (11).

$$\begin{array}{}\text{(10)}& {\mathrm{RC}}_{\mathrm{L}}={\displaystyle \frac{\left|X_{\mathrm{3}}-X_{\mathrm{2}}\right|}{\left|X_{\mathrm{3}}-X_{\mathrm{2}}\right|+\left|X_{\mathrm{2}}-X_{\mathrm{1}}\right|}}\times \mathrm{100}\mathit{\%}\end{array}$$

2.5 Data

The dataset used in this study consists of three parts: (1) hydrometeorological data; (2) geospatial data encompassing DEMs, soil types, and land uses; and (3) socioeconomic data encompassing per capita water consumption and population data.

Observed monthly streamflow data of the two hydrological stations in the study were collected for the years 1970–2000 from the Boluo station and Heyuan station, and the observed streamflow time series of these two hydrological stations have no missing data. Monthly inflow and outflow data of the three major reservoirs in DRB were also collected. All the hydrological data were obtained from the Guangdong Provincial Hydrological Bureau. Meteorological daily precipitation, temperature, and other data for 1968–2017 from 21 meteorological stations in the watershed were obtained from the National Meteorological Information Center of the China Meteorological Administration. Monthly actual ET data for SWAT model validation were obtained from the Amsterdam Evapotranspiration Model dataset with a spatial resolution of 0.25°×0.25° (Martens et al., 2017). Monthly soil moisture data for SWAT model validation were obtained from the European Centre for Medium-Range Weather Forecasts ERA5-Land dataset with a spatial resolution of 0.1°×0.1° (Hersbach et al., 2023). The actual evapotranspiration and soil moisture of the watershed equal the average of all grids included in DRB.

The 90 m resolution DEM data and 30 m resolution land use data at 10-year intervals (i.e., 1980, 1990, 2000, 2010, and 2015) are obtained from the Resource and Environment Science and Data Platform of the Chinese Academy of Sciences (Xu et al., 2018). Soil data are obtained from the 1 km resolution Harmonized World Soil Database from the Food and Agriculture Organization of the United Nations (Fischer et al., 2008).

The annual per capita integrated water consumption data of DRB from 2000 to 2017 were acquired from the Water Resources Bulletin of Guangdong Province. The population data in 2000, 2005, 2010, and 2015 were obtained from the 1×1 km spatial raster data of the Resource and Environment Science and Data Center of the Chinese Academy of Sciences (Xu, 2017b).

3 Results

3.1 Model performance

3.1.1 Streamflow reconstruction

The difference between the monthly average observed streamflow and the monthly average natural streamflow is small (Fig. 2). The monthly average measured streamflow and natural streamflow at the Heyuan station are 492.1 and 507.9 m³ s⁻¹, respectively, while the monthly average measured streamflow and natural streamflow at the Boluo station are 768.4 and 796.7 m³ s⁻¹, respectively. The difference between the measured streamflow and the natural streamflow mainly occurs in November, December, January, and February (when the measured streamflow is greater than the natural streamflow) and in May, June, and July (when the measured streamflow is less than the natural streamflow) (Fig. 2a and c).

Figure 2 Observed streamflow and natural streamflow processes at the Heyuan and Boluo stations from 1970 to 2000. (a) Annual distribution of streamflow at the Heyuan station, (b) streamflow process at the Heyuan station, (c) annual distribution of streamflow at the Boluo station, and (d) streamflow process at the Boluo station.

3.1.2 Model calibration and verification

The SWAT model shows sufficient accuracies in simulating streamflow, actual evapotranspiration, and soil moisture changes in DRB and can better simulate both seasonal and interannual changes in streamflow. During the calibration period, both stations achieved R² above 0.9, a NSE exceeding 0.8, and a PBIAS of less than 14 % (Fig. 3). Both stations had a simulated streamflow R² of greater than 0.8 during the validation period. The NSE values for streamflow simulation at the Heyuan station and Boluo station of the validation were 0.81 and 0.74, respectively. The model performs well in simulating the ET and soil moisture. Since the GLEAM ET data and ERA5 soil moisture data are raster data with a spatial resolution of 0.25×0.25°, considering the influence of data accuracy on the results, this study uses the watershed scale to validate the simulation results of ET and soil moisture. In the validation period, the R² and NSE for the simulation of evapotranspiration were 0.92 and 0.8, respectively (Fig. S2), while the R² and the NSE for the soil moisture simulation were both greater than 0.6. These validation results show that the model can be used to simulate hydrological regimes in DRB.

Figure 3 Simulated and observed monthly streamflow at the (a) Heyuan and (b) Boluo gauge stations during the calibration and validation periods.

3.2 LUCC and climate variability in DRB

The LUCC in DRB is mainly the decrease in cultivated land and the increase in urban land. The land use in DRB primarily consisted of forest land (18 875–18 833 km²), which forms more than 70 % of DRB. From 1980 to 2015, the urban land and water areas showed increases of 469.4 km² (137 %) and 17.4 km² (2.8 %), while the grassland, cultivated land, and forest land showed decreases of 41.3 km² (4.3 %), 487.5 km² (10.8 %), and 42.1 km² (0.2 %), respectively (Table 3).

Table 3 Land use transfer matrix in DRB from 1980 to 2015.

DRB exhibited significant regional differences in multiyear average precipitation, temperature, and potential evapotranspiration. The precipitation exhibited an increasing trend from the center to the south and north of DRB. The temperature and potential evapotranspiration showed an overall distribution pattern of greater values in the south and minor values in the north of DRB (Fig. 4). The multiyear average precipitation for the entire DRB was 1790.1 mm, with annual precipitation ranging from 1236.2 to 2567.5 mm. The regions with the highest multiyear average annual precipitation are located in the southeast of DRB, where annual precipitation exceeds 2200 mm, while the regions with the lowest precipitation are in the northeast of the watershed. The average annual temperature in DRB ranged from 19.5 to 21.3 °C, and the average annual potential evapotranspiration ranged from 1101.5 to 1320.6 mm. The south of DRB is predominantly urban and characterized by the urban heat island effect, while the north of DRB is mountainous with higher elevations, leading to the spatial distribution of temperatures.

Figure 4 Spatial distribution of the annual mean (a) precipitation, (b) temperature, and (c) potential evapotranspiration in DRB from 1960 to 2017. Each triangle represents the Mann–Kendall test result at a meteorological station.

The average temperature and potential evapotranspiration at DRB meteorological stations exhibited significant variations, while precipitation showed a relatively minor trend (Fig. 4). Overall, basin-averaged precipitation and potential evapotranspiration showed a nonsignificant decreasing trend, while temperatures showed a significant increasing trend. There was no significant change trend of precipitation for all stations in DRB (Fig. 4a). Twenty out of 21 meteorological stations in the region showed statistically significant increasing trends in average temperature, indicating a warming trend (Fig. 4b). Nine stations showed a significant decreasing trend in potential evapotranspiration that was primarily located in northern DRB (Fig. 4c).

The mean precipitation, temperature, and potential evapotranspiration of DRB can be obtained from the precipitation, temperature, and potential evapotranspiration of stations using the Tyson polygon method. The interannual variation of annual precipitation in DRB showed an insignificant decreasing trend (−0.51 mm a⁻¹). The annual mean temperature showed a significant increasing trend (0.024 °C a⁻¹). The annual potential evapotranspiration showed a significant decreasing trend (−0.38 mm a⁻¹) (Fig. S3).

3.3 Blue and green water resources

The average annual BW and GW were 1240.8 and 840.7 mm, respectively. The DRB water resources were dominated by BW, representing 60.1 % of the total water resources, and BW was 1.48 times higher than that of GW resources. The average GWF and GWS were 689.3 and 151.4 mm, respectively.

The annual BW resources in the subbasins of DRB ranged from 893.7 to 1990 mm, showing an increasing trend from the center to the south and north of DRB, aligning with the spatial distribution of precipitation (Fig. 5a). The regions with abundant BW resources are located in the central and southeastern parts of DRB (>1300 mm), and the BW in the upper reaches is comparatively low (<1100 mm). Differences in the spatial distribution of BW are primarily caused by differences in the spatial distribution of precipitation. Overall, the GWF and GWS are more evenly distributed in the subbasins than BW. The annual GWF in the subbasins of DRB ranged from 573.6 to 923.6 mm. The subbasins with higher GWF are primarily located in the Xinfengjiang Reservoir area in the middle reaches (>700 mm), while the low GWF subbasins are located in the southwest of DRB (<600 mm) (Fig. 5b). The land use in the subbasins where the Xinfengjiang Reservoir is located is primarily water areas, with a higher water evaporation rate than other regions, resulting in a greater GWF in this area than in the other regions. The annual GWS in the subbasins of DRB ranged from 126 to 190.6 mm. The subbasins with higher GWS are mainly located in the lower part of DRB (>150 mm) (Fig. 5c). The distribution pattern of GWS resources has a close relationship with the soil type of the watershed. The upper reaches and the northwestern part of the watershed are mostly red soil, while the middle and lower reaches are dominated by reddish soil. Reddish soil has a smaller water storage capacity than red soil, loses water faster, and has weaker water conservation and water supply performance than red soil. This is the primary factor in the north–south discrepancies in the number of GWS resources in DRB. In addition, the southern region mostly has large- and medium-sized cities. As urban construction land expands, the land use type in the region gradually changes to urban land and industrial land and the solidification of road surfaces reduces the area of bare soil in the region, resulting in a decrease in GWS resources. The annual GWI (Fig. 5d) showed a spatial pattern opposite to BW, decreasing from 0.45 in the upper reaches to 0.3 in the lower reaches. The highest GWI is found in the upper reaches, which is due to the relatively low rainfall in the upper reaches and the lush vegetation with significant plant interception and transpiration, resulting in a higher proportion of total evapotranspiration than in the middle and lower reaches. The central part of the basin has the highest precipitation, leading to a lower GWI. The southern part of the watershed has the highest temperature, and the evapotranspiration is high. Meanwhile, the lower reaches have a large proportion of agricultural and urban land, and crop irrigation can increase evapotranspiration.

Figure 5 Spatial distribution of the mean (a) BW, (b) GWF, (c) GWS, and (d) GWI in DRB over 1970–2017.

In DRB, there was no significant increasing trend in either BW or GWS, while GWF exhibited a significant decreasing trend. The annual trend rate of BW in DRB was 0.14 mm a⁻¹, with an annual fluctuation range of 713.6–2017.5 mm during 1970–2017. The minimum BW occurred in 1991, while the maximum BW was recorded in 2016 (Fig. 6a). The GWF in DRB from 1970 to 2017 exhibited a significant decreasing trend (−0.57 mm a⁻¹) (Fig. 6b). The minimum GWF occurred in 2005 (603.1 mm), while the maximum GWF was recorded in 1974 (721.3 mm). In contrast, the GWS in DRB from 1970 to 2017 has been slowly increasing at a rate of 0.015 mm a⁻¹ (Fig. 6c). The annual fluctuation in GWS was smaller than BW and GWF. The GWI in DRB from 1970 to 2017 showed no significant decreasing trend at a rate of −0.0003 % a⁻¹ (p>0.05) (Fig. 6d), implying that the redistribution of precipitation in DRB might change slowly.

Figure 6 Interannual variation of (a) BW, (b) GWF, (c) GWS, and (d) GWI in DRB during 1970–2017.

3.4 Blue and green water scarcity

The average blue water scarcity level in DRB was low (22.4 %) during 1970–2017. The blue water scarcity levels in various subbasins ranged from 0.1 % to 206 %. The multiyear average blue water scarcities, except for one subbasin in the southwest, were all low (<100 %) (Fig. 7a). This indicates that blue water scarcity is not common in DRB at the annual scale. Regions with relatively high blue water scarcity (>20 %) are mostly situated in the upper reaches of various tributaries within the watershed, where river streamflow is relatively low. The area with the highest blue water scarcity (206 %) is located in the 63rd subbasin of Shenzhen and Huizhou, reaching a moderate level of blue water scarcity. This region has a large population, with a much higher blue water demand than other areas. Additionally, this subbasin is situated in the upper reaches of the primary tributary of DRB, resulting in a limited supply of BW resources. Although the northern parts of subbasins 55, 56, and 61 have large populations, these subbasins are situated downstream of the main Dongjiang River with a higher streamflow, leading to lower BWSC levels. The average GWSC in the entire basin from 1970 to 2017 was low (41.4 %). The blue water scarcity levels in the various subbasins ranged from 31 % to 104 %. The vegetation cover in DRB is high, and DRB thus has relatively high rates of vegetation transpiration and interception evaporation. The basin experiences a GWSC of nearly 50 %, indicating a potential occurrence of GWSC. The areas with higher GWSC are primarily situated in the middle reaches for DRB (Fig. 7b), where water surface evaporation is high, resulting in their GWSC exceeding 100 %. The evaporated water in these areas originates from the reservoirs, not the soil, leading to an overestimation of the GWSC in these subbasins.

Figure 7 Spatial distribution of the mean (a) BWSC, (b) GWSC, and (c) FLK index in DRB over 1970–2017.

Furthermore, the FLK index was used to quantify population-driven water resource scarcity. F1–F4 represent absolute scarcity, scarcity, stress, and no stress, respectively. The results showed that most regions in DRB have no water scarcity pressure (Fig. 7c). However, the 63rd subbasin experienced absolute water scarcity and the 52nd subbasin experienced water scarcity. There were six lower-reach subbasins and four upper-reach subbasins facing water stress. DRB receives ample precipitation, resulting in a relatively large river flow and generally leading to a higher FLK index. As a result, the basin faces lower water resource pressure.

Figure 8 Hierarchical clustering tree of (a) BWSC and (b) GWSC.

This study also further identified hotspots of BWSC and GWSC in DRB using hierarchical clustering of BWSC and GWSC in each subbasin. Figure 8 shows the clustering tree results for BWSC and GWSC. When the standardized distance was set to 500, all the subbasins could be divided into four categories according to blue water scarcity: (1) the first category consisted of 27 subbasins, such as 32, 56, and 28, where the blue water scarcity level was lowest (<20 %). (2) The second category comprised subbasin 63, which has the most severe blue water scarcity (206 %). (3) The third category comprised seven subbasins, such as 52, 58, and 60, all located in the lower reaches, with relatively high blue water scarcity levels (40 %–100 %). These subbasins are mostly located in the tributaries of the lower reaches, with a relatively large population and smaller river streamflow compared to the mainstem of the Dongjiang River. (4) The fourth category consisted of 28 subbasins, such as 59, 62, and 8, with blue water scarcity levels ranging from 20 % to 40 %. Similarly, hierarchical clustering was conducted for GWSC. When the standardized distance was set to 500, GWSC in the subbasins could be divided into three categories: (1) the first category consisted of 56 subbasins, such as 37, 56, and 29, with relatively low GWSC levels, all below 50 %, indicating low GWSC. (2) The second category consisted of subbasins 32 and 33, where the predominant land use type was water areas, leading to higher GWSC due to high water surface evaporation. (3) The third category consisted of subbasins 47, 31, 54, 30, and 36, where the water area proportion in these subbasins was larger than in the others, leading to significant influences from water surface evaporation. Figure S4 shows the annual variation of blue water scarcity and green water scarcity in the basin. Except for some subbasins, the blue and green water scarcity in most of them is less than 50 %. The degree of green water scarcity is higher than that of blue water scarcity in most of the subbasins. Only subbasin 63 downstream experienced a severe blue water scarcity.

Figure 9 Interannual variation of (a) BWSC, (b) GWSC, and (c) the FLK index in DRB during 1970–2017.

The interannual variations in BWSC and GWSC in DRB showed distinct regional differences. The BWSC in the basin slowly increased at a rate of 0.3 % a⁻¹ (Fig. 9a). The BWSC in the lower reaches slowly increased at a rate of 1.1 % a⁻¹, while the BWSC in the upper and middle reaches slowly decreased at −0.47 % a⁻¹ and −0.1 % a⁻¹, respectively. The GWSC in the upper, middle, and lower reaches of DRB showed a decreasing trend, with the basin-scale GWSC decreasing significantly at a rate of −0.04 % a⁻¹ (Fig. 9b). Despite the acceleration of urbanization and a significant increase in population in the middle and lower reaches of the watershed, blue water availability and the amount of obtainable BW have been increasing. Additionally, the annual per capita water consumption in the basin has decreased from 481.0 m³ in 2000 to 245.0 m³ in 2020. As a result, the rate of increase in BWSC in the watershed has been relatively small. In contrast, the GWF in DRB demonstrated a significant decreasing trend, and the GWS increased slowly. Therefore, the GWSC in DRB demonstrated a significant decreasing trend. Meanwhile, the FLK index of the watershed showed a significant decreasing trend (−285.3 m³ yr⁻¹), which means that the per capita water resources in the watershed have significantly decreased (Fig. 9c). This is due to the rapid population growth in the watershed and the slow increase in available water resources.

3.5 Impacts of LUCC and climate change on blue and green water

To examine the impacts of climate change and LUCC on BW and GW change, this study set three climate conditions and land use scenarios to explore this effect by comparing the scenarios (Table 3). The combined impacts of climate change and LUCC on BW and GWS in DRB were superimposed, and the combined effect on GWF was a negatively synergistic effect. Figure 10 shows the variations in BW and GW under the impacts of climate change (S2–S1) and LUCC (S3–S2), their combined effects (S3–S1), and the relative contributions of climate change and LUCC to the BW and GW changes in DRB during 1970–2017. Under the joint influences of climate change and LUCC, BW decreased by 4.5 mm a⁻¹. Of this decrease, climate change resulted in a loss in BW of 3.9 mm a⁻¹, contributing 88.0 %, while LUCC led to a loss in BW of 0.5 mm a⁻¹, contributing 12.0 % (Fig. 10a). The effect of climate change on BW variation is much greater than that of LUCC at the basin scale. Under the combined influences of climate change and LUCC, GWF decreased by 17.0 mm a⁻¹. Of this decrease, climate change accounted for a decrease in GWF of 19.5 mm a⁻¹, contributing 88.5 % to the decrease, while LUCC led to an increase in GWF of 2.5 mm a⁻¹, contributing 11.5 % (Fig. 10b). Overall, the influence of climate change on GWF changes in the watershed is significantly more pronounced than that of LUCC. Under the joint influences of climate change and LUCC, GWS increased by 0.7 mm a⁻¹. Of this increase, climate change contributed to an increase in GWS of 0.3 mm a⁻¹, accounting for 39.4 %, while LUCC contributed to an increase in GWS of 0.4 mm a⁻¹, accounting for 60.6 % (Fig. 10c). DRB is situated in a humid region with high GWS, resulting in small fluctuations of GWS in response to precipitation changes. The fluctuations of GWS are primarily influenced by soil properties and land use. In general, the effect of climate change on the GWS change in DRB is smaller than the effect of LUCC.

Figure 10 Effects and relative contribution of climate change and LUCC to the changes in (a) BW, (b) GWF, and (c) GWS in DRB during 1970–2017.

Under the coupled influences of climate change and LUCC, the BW and GW resources in DRB have changed. However, there were differences in the joint impacts of climate change and LUCC on BW and GW. Both climate change and LUCC have led to the decrease in BW in the watershed, and the combined effect of climate change and LUCC on BW equals the sum of their individual effects. Climate change, such as a decrease in potential evapotranspiration, has resulted in a decrease in GWF in DRB, while LUCC has led to an increase in GWF. Therefore, the joint impacts of climate change and LUCC on GWF were partially offset, resulting in the joint impacts of climate change and LUCC on GWF being less than the sum of their absolute individual effects. Both climate change and LUCC have led to an increase in GWS in DRB, and the joint impacts of climate change and LUCC on GWS equal the sum of their individual effects.

4 Discussion

This study used the SWAT model to simulate changes in BW and GW resources in DRB over the past 5 decades and their responses to climate change and LUCC. It also assessed water resource security in the basin. The findings revealed that GWF exhibited a decreasing trend and that BW and GWS exhibited an increasing trend. Liu et al. (2010) similarly found an increasing trend in annual surface runoff in DRB. Potential evapotranspiration in DRB showed a decreasing trend, which may be the main cause of the significant decrease in GWF in the basin (Fig. S3), and similar conclusions are drawn in He et al. (2013).

We show that water resources in DRB are dominated by BW with a mean annual GWI of 0.4, which is the same as what many studies show in humid areas (Nie et al., 2023). Although the GWI in humid areas is much smaller than that in arid areas, the ratio of GW in DRB still reaches 40 %, so it is imperative to incorporate GW into the water resource assessment system. The GWI in the upper and middle reaches of DRB exceeded 0.4, while that in the lower reaches was only about 0.3. These results mean that, to ensure the appropriate utilization of water resources, effective water management in the upper and middle reaches of DRB should consider GW planning, while water management in the lower reaches should mainly consider BW. The assessment results of BWSC and GWSC in DRB similarly illustrate this issue. The GWSC in the upper and middle reaches was bigger than that in the lower reaches of DRB, while the BWSC in the lower reaches of DRB was bigger than in the upper and middle reaches (Fig. 9).

There are robust correlations between BW and precipitation, GWF, and potential evapotranspiration in DRB. Climate change plays a dominant role in variations of BW and GWF. BW is more sensitive to precipitation and potential evapotranspiration. GWF shows sensitivity to changes in potential evapotranspiration, and GWS is influenced by both precipitation and potential evapotranspiration (Ahiablame et al., 2017; He et al., 2015). Of course, some studies in arid regions show that GWF is mainly affected by precipitation (Ahiablame et al., 2017), which may be linked to the hydrothermal conditions of the basin. There is sufficient precipitation in DRB, where the GWF changes are mainly energy-limited and the effect of precipitation on GWF is smaller.

Although BW and GW are mainly affected by climate change, the influences of LUCC on them cannot be ignored. The reaction of water resources to LUCC is exceedingly intricate and involves various hydrological processes, including runoff yield, infiltration, and groundwater (Cuo, 2016; Zhang and Shangguan, 2016). As there is a strong compensatory effect of diverse land use in the hydrological system, particularly in expansive watersheds, this could create a strong resistance to GW and BW conversion (Lin et al., 2015). A decrease in forest land or an increase in cultivated and urban land could lead to a rise in BW and a decline in GW in the watershed. Veettil and Mishra (2018) demonstrated that there is a 10 % rise in forest land cover and a 1.4 % drop in BW, indicating a negative elasticity between the two. However, the effect of urban land on streamflow in different periods showed the opposite effect. On the one hand, the increase in urban land results in increases in impermeable area and thus surface runoff in the basin, but at the same time the increase in urban land may also reduce groundwater discharge to streamflow. At the same time, LUCC often results in changes in vegetation. Vegetation variations affect the water cycle by altering canopy interception (Shao et al., 2018; Wu et al., 2019b), transpiration (Chen et al., 2023), and canopy evaporation and ameliorating soil structure (Qiu et al., 2022). Therefore, increasing vegetation often increases infiltration and soil moisture and reduces surface runoff.

There are several limitations and uncertainties in this study. (1) Since the quantity of the BW and GW is derived from the output results of the model simulations, including water yield, ET, soil moisture, and groundwater, the precision of the outcomes depends largely on the precision of the model simulations. Given the absence of observed evapotranspiration and soil moisture data for DRB, this study calibrated and validated the SWAT model using only monthly streamflow, which may weaken these results to some extent. To enhance the credibility of the model, this study also utilized widely used actual evapotranspiration data (GLEAM) and soil moisture (ERA5-Land) during model validation at a basin scale. The findings indicated that the simulation performance is relatively good and meets the accuracy requirements for simulation. (2) Climate change, LUCC, and large reservoir operation are the primary factors influencing the changes in hydrological conditions in DRB. The contributions of reservoir regulation, LUCC, water resource utilization, and climate change to the distribution of intra-annual flow are 33.5 %, −9 %, 4.5 %, and 1 %, respectively, during 1956–2009 (Tu et al., 2015). The operation of reservoirs, including large reservoirs like the Xinfengjiang Reservoir, is one important reason for hydrological changes in DRB (Lin et al., 2014; Zhang et al., 2015). The reservoir module was not established when constructing the SWAT model in this study. To obtain natural BW and GW volumes in the watershed and mitigate the impact of hydraulic engineering, reconstructed natural streamflow based on observed flow was utilized for model calibration and validation. However, hydraulic engineering significantly influences the annual allocation of BW. The flow restoration considered the impacts of the three major reservoirs on the Dongjiang River and did not consider the impacts of other minor hydraulic projects and human water consumption. (3) The calculations of BWSC and the FLK index both include environmental flows. This study represented the proportion of environmental flow in streamflow as 80 %. Some studies have suggested that assuming environmental flow to be 80 % of the total water resources in a basin may overestimate water scarcity (Liu et al., 2017; Richter et al., 2012). Therefore, we varied the proportion of environmental flow and assessed the degree of BWSC using 60 % and 70 % proportions. Results show that only the 63rd subbasin changed from severe BWSC to moderate to high BWSC, while other subbasins remained with low BWSC. Therefore, the threshold for environmental flow has a minor impact in this paper. The assessment of BWSC and per capita water resources did not take into account the water demand of cities such as Shenzhen and Hong Kong, although the water supply for these cities primarily comes from the Dongjiang River through the Dongjiang–Shenzhen Water Supply Project. (4) The hydrological modeling approach utilized in this study is a frequently used method for quantitative analysis of attribution. Nevertheless, it implies independence between climate change and LUCC and does not adequately distinguish between the impacts of these two components. Such restrictions are diffusely recognized as existing (Dey and Mishra, 2017). Despite this recognized limitation, hydrological modeling methods have been widely used in numerous similar studies, yielding credible results (Li et al., 2021; Nie et al., 2023).

5 Conclusion

This study analyzed the spatiotemporal evolution of BW and GW, assessed the water security, and evaluated the effects of climate change and LUCC on BW and GW in DRB using the SWAT model. The conclusions can be outlined as follows:

1. During 1970–2017, grassland, cultivated land, and forest land in DRB decreased by 4.3 %, 10.8 %, and 0.2 %, respectively, while urban land and water areas increased by 137 % and 2.8 %, respectively. The annual precipitation and potential evapotranspiration showed a nonsignificant decreasing trend, while the annual average temperature showed a significantly increasing trend.

2. The annual BW, GWF, and green storage in DRB from 1970 to 2017 were 1240.8, 840.7, and 151.4 mm, respectively. BW (0.14 mm a⁻¹) and GWS (0.015 mm a⁻¹) in DRB showed no significant increasing trend, and GWF (−0.57 mm a⁻¹) showed a significant decreasing trend.

3. The levels of the annual BWSC and GWSC in DRB were low, and per capita water resources exceeded 1700 m³ per capita per year. BWSC displayed a nonsignificant increasing trend, while GWSC and the FLK index displayed a significant decreasing trend, especially in the lower reaches.

4. Climate change was the major driving factor of changes in BW and GWF, and LUCC was the major driving factor of GWS change. Climate change contributed 88.0 %, 88.5 %, and 39.4 % of the changes in BW, GWF, and GWS in DRB. Both climate change and LUCC decrease (increase) BW (GWS), while climate change (LUCC) decreases (increases) GWF in DRB.