Abstract
Urbanization has a strong signal on the hydrologic cycle, leading to reduced infiltration and faster and larger runoff. However, the detailed analysis of the contribution of urbanization to different quantiles of discharge is still lacking; particularly, less is known in watersheds that have been experiencing such large and rapid urbanization as those in China. Here, we focus on the Wenyu watershed, a rapidly urbanizing basin located in the Beijing metropolitan area. Using a statistical attribution framework, we examine the hydrological response to the increasing urbanization across a wide range of discharge quantiles, from low to high flows; moreover, we perform analyses at the seasonal scale to capture differences in the physical processes at play during the year. Results indicate that the selected GAMLSS (generalized additive models for location, scale, and shape) models can capture well the variability in streamflow in this highly urbanized basin, with the average Spearman correlation coefficients between observations and the median of the fitted models of 0.84, 0.79, and 0.81 in spring, summer, and winter, respectively. Overall, urbanization played a different role in the different seasons and discharge quantiles. More specifically, we find its strongest impact to be in winter and spring, and for low and median quantiles. The role of precipitation is the strongest in summer, and it increases as we move toward the upper tail of the discharge distribution, especially above the 55th percentile, for which precipitation is selected as the only important predictor. Recycled water, on the other hand, tends to play a more dominant role in winter and spring.
HIGHLIGHTS
Urbanization has a strong signal on the hydrologic cycle; however, less is known in watersheds that have been experiencing such large and rapid urbanization. In this research, we focus on a fast urbanizing basin, and the hydrological responses to the increasing urbanization have been examined.
The quantification of the impacts of drivers for different discharge quantiles and seasons has been examined.
Urbanization plays a decreasing role towards the upper tail of the distribution.
Graphical Abstract
INTRODUCTION
Urban watershed hydrology has attracted significant attention, especially over the past few decades (e.g., Smith et al. 2013; Pumo et al. 2017; Arnone et al. 2018; Hopkins et al. 2019; Kelleher & Mcphillips 2019; Shen et al. 2019; Blum et al. 2020; Dudley et al. 2020) because of its significant societal and economic implications. In fact, 2009 was the first year that more than half of the world's population was living in urban and suburban areas (United Nations 2010); moreover, the proportion of the urban population is projected to exceed 80% by 2030 (Lee & Heaney 2003). Land-use modifications and water utilization and reuse, accompanied by the concentration of population and industries, can have a significant impact on the regional hydrological cycles (e.g., Defries & Eshleman 2004). Therefore, it is important that we include the role of urbanization in describing the changes in streamflow across a wide range of discharges (i.e., from low to high flows).
Previous studies examined the different impacts of urbanization (e.g., Damodaram et al. 2010; Jacobson 2011; Qin et al. 2013; Pumo et al. 2017; Arnone et al. 2018; Hou et al. 2019), finding that it leads to increased urban flood risk and decreased baseflow and groundwater recharge in peri-urban catchments (e.g., Hamel et al. 2013; Miller et al. 2014; Prosdocimi et al. 2015). Sheng & Wilson (2009) examined the effects of urbanization on streamflow in the Los Angeles area and found that 90% of the storm rainfall became runoff in urban watersheds, while the forested watersheds retained 25% of the rainfall. Not only does urbanization amplify the adverse effects of flooding by increasing runoff from impervious surfaces and blocking water flow, but it also increases flood risk through the occupation of floodplains (e.g., Yin et al. 2015).
There are two main approaches to evaluating the role of urbanization in the hydrological cycle. One method is based on the use of hydrological and hydraulic models. Pumo et al. (2017) used a real-time integrated basin simulator based on a triangular irregular network to compare and contrast the role of urbanization and climate change on the hydrology of two watersheds, highlighting the interaction between the two drivers (see also Arnone et al. (2018) for results related to extreme discharge events). Suriya & Mudgal (2012) used a hydrological modeling system (HEC-HMS) to assess the impact of different land uses on the flooded areas and the water depth for a given rainfall amount. They found that the increased impervious areas led to higher flood peaks and volumes. Seidl et al. (2020) used the Storm Water Management Model (SWMM) to evaluate the runoff behavior of four urban catchments in Brazil with different land uses, finding much higher runoffs in the most urbanized watershed. Wang et al. (2020) used the soil and water assessment tool (SWAT) to assess the impact of land-use changes and climate change on extreme hydrological events in the Yangtze River delta, China. They found that land-use changes are expected to decrease low flows and increase high flows. Zope et al. (2015) used HEC-HMS to evaluate the impact of land use on the flood hydrograph; they found marginal increases in the peak runoff and volumes for increasing urbanization. Zhang et al. (2020) developed a hydraulic model to evaluate flood inundation to account for the rapid urbanization in the Luojiang river basin, representing China's southeastern coastal watersheds. One potential limitation of these models is represented by the difficulties in their calibration and the large requirements in terms of input data to characterize all the different processes at play in great detail.
Another approach to evaluate the role of urbanization is through statistical models. Even though this approach is not able to describe all the detailed processes that influence urban hydrology, these models are relatively straightforward to develop and have shown good performance (e.g., Villarini et al. 2009b; Prosdocimi et al. 2015; Slater & Villarini 2017; Zhang et al. 2018; Blum et al. 2020). Although previous work used the statistical models to broadly assess the impact of urbanization and climate change on the urban hydrological cycle (e.g., Yang et al. 2013), most studies have focused on their impacts on extreme events (e.g., extreme precipitation and extreme floods) in urbanized watersheds. A comprehensive study on the hydrological response of urbanized watersheds is still lacking, especially for basins that have experienced extremely rapid changes.
Here, we focus on the Wenyu river basin, which is known as the ‘mother river’ of Beijing. Its urbanized area has increased from about 4% in 1985 to approximately 52% in 2016, and it hosts the largest portion of human activities in Beijing. Beijing, as the political, cultural, and economic center of China, is one of the most urbanized and economically developed cities in China. Urban expansion is known to affect the hydrological cycle system in urban regions, and the policy associated with ‘Sponge City’ was proposed in 2013 in China for flood control and waterlogging relief in urban areas. Therefore, the Wenyu river basin represents an ideal case study to examine the hydrological response to such massive urbanization. In particular, we will develop a statistical framework to attribute the observed changes in discharge properties to drivers associated with climate and land use/land cover. Rather than focusing only on extremes, we will perform analyses across a wide range of quantiles of the discharge distribution (i.e., from low to high flows) and at the seasonal level. The objective of this study is to better understand the relationship between urbanization and streamflow in this highly urbanized watershed, which can provide the scientific basis for flood management and China's sponge city construction and contribute to a more sustainable and integrated water management.
DESCRIPTION OF THE STUDY AREA AND DATA
Study area
The Wenyu river basin is located in the central area of Beijing, China (Figure 1), with a total length of 47.5 km and a drainage area of 2,478 km2 (Zhang et al. 2011). The climatic conditions of this river basin belong to the continental monsoon climate, with hot and rainy summers and cold and dry winters. The average annual precipitation is 619 mm, with most of the precipitation mainly concentrated between June and September, and accounting for about 84% of the entire year (Liao et al. 2018; Yang et al. 2018). The Wenyu river basin hosts the most intensive human activities in Beijing and has been experiencing extreme urbanization over the past three decades. With the rapid development of Beijing since the 1980s, its land use has undergone a significant change, with the impervious area increasing from about 4% in 1985 to approximately 42% in 2016 (Figure 2).
Boundary, river network, and distribution of hydrological and meteorological stations in the Wenyu river basin.
Boundary, river network, and distribution of hydrological and meteorological stations in the Wenyu river basin.
Spatial distribution of impervious areas (in red) in the Wenyu river basin from 1985 to 2016. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/wcc.2022.041.
Spatial distribution of impervious areas (in red) in the Wenyu river basin from 1985 to 2016. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/wcc.2022.041.
Data description
We used the daily streamflow data at the Beiguan dam (the location of this station is given in Figure 1), covering the period from 1985 to 2016. The streamflow data were provided by the Beijing Hydrology Center, and their quality was checked to ensure there were no missing streamflow data. The whole year was divided into four seasons: spring (March to May), summer (June to August), fall (September to November) and winter (December to February).
The weighted value of the average daily precipitation at the 10 stations (Figure 1) was calculated as the basin average precipitation and then accumulated to obtain the seasonal precipitation. Figure 3 summarizes the inter-annual variability in precipitation with Sen's slope methods (Sen 1968), highlighting decreasing trends in summer and increasing trends in fall, which are significant at the 0.01 level. The seasonal precipitation in spring has a slightly downward trend, whereas that in winter has a slightly upward trend during the study period (neither of them statistically significant).
Accumulated seasonal precipitation for (a) spring, (b) summer, (c) fall, and (d) winter. The red lines represent the trend computed based on Sen's slope estimator, with S representing the value of Sen's slope. The solid (dashed) lines represent results that are (not) significant at the 0.01 level (also indicated by ‘**’). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/wcc.2022.041.
Accumulated seasonal precipitation for (a) spring, (b) summer, (c) fall, and (d) winter. The red lines represent the trend computed based on Sen's slope estimator, with S representing the value of Sen's slope. The solid (dashed) lines represent results that are (not) significant at the 0.01 level (also indicated by ‘**’). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/wcc.2022.041.
We used the essential urban land-use categories in the China (EULUC-China; available from http://data.ess.tsinghua.edu.cn; Gong et al. 2020) dataset to capture the changes of the artificial impervious areas in the Wenyu river basin. This dataset was derived from the Google Earth Engine platform at a spatial resolution of 30 m.
Temperature data were derived from the China Meteorological Forcing Dataset (CMFD), which was developed by the Institute of Tibetan Plateau Research of the Chinese Academy of Sciences (http://westdc.westgis.ac.cn/data/7a35329c-c53f-4267-aa07-e0037d913a21). The CMFD is a gridded dataset covering the period from January 1979 to the present, with temporal and spatial resolutions of 3 h and 0.1°. The CMFD dataset, which is derived from the fusion of different sources (e.g., remote sensing, reanalysis, and observation), includes seven variables (i.e., precipitation rate, near-surface temperature, near-surface pressure, near-surface air-specific humidity, near-surface wind speed, downward short-wave, and long-wave radiation) (Ren et al. 2018). Previous studies showed that the CMFD dataset is one of the most widely used climate datasets for China because of its long record and good quality (He et al. 2020). We used this dataset to compute an average seasonal temperature for the basin.
We also used the annual recycled water amount data for the basin covering the period 1985–2016; this quantity includes the amount of treated domestic and industrial sewage. We used the annual groundwater level in Beijing to reflect the groundwater situation in the Wenyu river basin, which was measured at the end of every year. The recycled water amount and the groundwater-level data are provided by the Beijing Hydrology Center.
We considered precipitation, impervious areas, temperature, and groundwater levels as the explanatory predictors in summer and fall. The recycled water amount was used as a covariate in the modeling of the winter and spring discharge because there is limited winter and spring rainfall in Beijing, and also because this is one of the largest rivers in the area to receive the recycled water from the sewage treatment plants of Beijing. Furthermore, the relationship between precipitation and discharge was severely affected by the water resource management in winter and spring to meet the storage requirement in this highly urbanized basin. An impervious area is not included as a potential predictor in winter and spring because of its strong relation to recycled water amounts. We have normalized all predictors, so that they spanned the same range.
METHODOLOGY
Development of GAMLSS models
We use the generalized additive models for location, scale, and shape (GAMLSS; Rigby & Stasinopoulos 2005) to describe the relationship between predictors and predictand. We selected these models because of their flexibility in terms of distributions and dependence of the parameters on the covariates. These models have been used extensively in the hydrometeorological literature (e.g., Villarini et al.2009a, 2009b; Machado et al. 2015; Zhang et al. 2016; Li et al. 2017; Slater & Villarini 2017; Li et al. 2018; Zou et al. 2018). Here, we focus on two distributions (gamma and lognormal) because they are bounded at zero and can cover a broad range of shapes, from highly skewed to almost symmetric. All analyses are performed with the gamlss package in R (Rigby & Stasinopoulos 2005; Stasinopoulos & Rigby 2007). For each season, our response variables are represented by the quantile of the discharge distribution ranging from Q0.05 to Q0.95 (with an interval of 0.05), similar to Slater & Villarini (2017). The predictors are the seasonal precipitation, antecedent wetness (i.e., it is the precipitation for the season prior to the one of interest and represents a simple proxy for antecedent conditions), seasonal temperature, annual percentage of impervious area, annual recycled water amount, and annual groundwater level. Let us consider Y as the Q0.5 (i.e., median of the seasonal discharge distribution in a given season). We consider a total of 31 models in summer and fall (Supplementary Table 1), and 23 models in spring and winter (Supplementary Table 2), in which the parameters μ and σ vary as a function of the different predictors (with the exception of model 1, in which the parameters are constant).
Selection of the best models
To select a model among those in Supplementary Tables 1 and 2, we used the Schwarz–Bayesian Criterion (SBC; Schwarz 1978). Then, we evaluated the quality of these selected models by examining the normality of the residuals and the visual examination of the model's fit. In particular, we examined the worm plots (Buuren & Fredriks 2001) and computed the first four moments of the residuals and their Filliben correlation coefficient. We also computed the Pearson and Spearman correlation coefficient between the median of the fitted models and the observations, and performed leave-one-out validation to assess the robustness of our results.
Impact of the urbanization
RESULTS
For each season and discharge quantile, we fitted the models in Supplementary Tables 1 and 2 for both the gamma and lognormal distributions. In the main text, we will focus on the results for the Q0.05, Q0.5, and Q0.95, and the results for the other quantiles are provided in the Supplementary Material. As an example, Figure 4 shows the results for each of the four seasons and selected discharge quantiles. Overall, the low quantiles show an increasing trend in each season (Figure 4, left column), especially since the beginning of the 21st century. The median of the discharge distribution in spring and winter shows a large increasing trend, with a much more muted trend in the other two seasons (Figure 4, middle column). The results for the 95th percentile do not show trends over the study period, with the exception of the spring for which we see an increasing trend.
Modeling of the Q0.05 (left column), Q0.50 (middle column), and Q0.95 (right column) in each of the four seasons (by row). The red circles are the observations. The black line represents the median of the fitted model, while the dark (light) gray areas represent the areas between the 5th and the 95th (25th and 75th) percentiles. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/wcc.2022.041.
Modeling of the Q0.05 (left column), Q0.50 (middle column), and Q0.95 (right column) in each of the four seasons (by row). The red circles are the observations. The black line represents the median of the fitted model, while the dark (light) gray areas represent the areas between the 5th and the 95th (25th and 75th) percentiles. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/wcc.2022.041.
We found that the gamma distribution was mostly selected in spring and summer, while the lognormal distribution was in the fall; during winter, the gamma distribution was selected for the low quantiles, while the lognormal was for the high ones. Overall, these models can well describe the year-to-year variability in discharge, as well as the overall tendencies. Their goodness of fit is further supported by the other metrics we considered. Supplementary Table 3 shows that the residuals have a mean close to 0, a variance close to 1, a coefficient of skewness close to 0, and a coefficient of kurtosis close to 3; moreover, the Filliben correlation coefficients are all close to 1. These quantitative results are further supported by the worm plots, which are shown in Figure 5 for the Q0.05, Q0.5, and Q0.95.
Worm plots for the modeling of the Q0.05 (left column), Q0.50 (middle column), and Q0.95 (right column) in each of the four seasons (by row).
Worm plots for the modeling of the Q0.05 (left column), Q0.50 (middle column), and Q0.95 (right column) in each of the four seasons (by row).
We computed the Spearman and Pearson correlation coefficients between observations and the median of the fitted models as a measure of accuracy (Figure 6(a) and 6(b)). In general, the models tend to perform better in spring, summer, and winter, with the average Spearman correlation coefficients of 0.84, 0.79, and 0.81, respectively (the corresponding results based on the Pearson correlation coefficients are 0.81, 0.79, and 0.76). The performance decreases in the fall, where the average Spearman and Pearson correlation coefficients are 0.40 and 0.41, respectively. To examine the robustness of our results, we used leave-one-out cross-validation. As shown in Figure 6(c) and 6(d), the very good performance of our models is retained in all the seasons, with the exception of the fall.
Summary of the correlation coefficient between observations and median of the fitted distribution for different discharge quantiles and seasons. The left (right) panels show the results based on the Spearman (Pearson) correlation coefficient. The results in the top row are based on the entire period, while those in the bottom row are based on leave-one-out cross-validation.
Summary of the correlation coefficient between observations and median of the fitted distribution for different discharge quantiles and seasons. The left (right) panels show the results based on the Spearman (Pearson) correlation coefficient. The results in the top row are based on the entire period, while those in the bottom row are based on leave-one-out cross-validation.
Up to this point, we have focused our attention on the performance of our fitted models. Given that these models can describe the observations well, we can focus on examining the role of the different predictors. Figures 7 and 8 summarize the dependence of the location and scale parameters μ and σ on the different predictors for the different seasons and quantiles. The intercept of μ increases for increasing discharge quantiles, capturing the increasing discharge magnitude. On the other hand, the intercept of σ does not present consistent tendencies across seasons and quantiles (Figure 8, left column), which is consistent with previous studies (e.g., Villarini & Strong 2014; Li et al. 2017).
Values of the estimated intercept (left panels) and coefficients for the different covariates (right panels) for the location parameter μ in each season.
Values of the estimated intercept (left panels) and coefficients for the different covariates (right panels) for the location parameter μ in each season.
Values of the estimated intercept (left panels) and coefficients for the different covariates (right panels) for the scale parameter σ in each season.
Values of the estimated intercept (left panels) and coefficients for the different covariates (right panels) for the scale parameter σ in each season.
Figure 7 (right column) shows the dependence of the location parameter μ on the different predictors. For spring, the coefficients of recycled water amount, which decreases for increasing discharge quantiles, are higher than those of antecedent wetness. This indicated that the impact of recycled water amount on lower quantiles is larger than that on high quantiles, and it tends to be larger than the role of antecedent wetness. For summer, the values of the coefficients of precipitation showed increasing trends as we move from low to high flows, while the values of the coefficient for the impervious area tend to follow the opposite trajectory. This indicates that the impact of precipitation tends to increase as we focus on higher and higher discharge quantiles, while the opposite is true for impervious areas; this is particularly true above the 55th percentile, for which impervious area is no longer selected as an important predictor. For fall, precipitation, impervious area, and antecedent wetness are the three selected covariates. The impervious area is selected only up to the median, further suggesting that the role of urbanization is mostly felt in the low to median flow, rather than on the upper tail. As we focus on the higher quantiles, precipitation and antecedent wetness tend to play a much more significant role. For winter, similar to the results for spring, we find that recycled water amount is the most important predictor, with overall larger values for lower quantiles.
Figure 8 (right column) shows the dependence of the scale parameter σ on the different predictors. In the spring, the recycled water amount has a negative coefficient for all quantiles, and antecedent wetness is selected only for discharge quantiles from 0.30 to 0.45. Precipitation, impervious area, and temperature have negative values in the summer, while temperature is selected only in the middle of the discharge distribution; the parameter σ becomes constant when we move beyond the 80th percentile. In the fall, precipitation and antecedent wetness are the covariates selected most consistently, especially up to the median. Finally, in the winter, there is an overall weak dependence of the σ parameters on covariates, with the exception of the higher quantiles, for which recycled water amount and temperature are selected.
Figure 9 summarizes the cumulative contribution of urbanization for different discharge quantile, season, and year. There are some interesting patterns that we can pick out. Across all seasons, as we get closer to the present time, the cumulative impact increases (except 2004 and 2005 in winter and spring because the recycled water amount in these 2 years is relatively low), reaching values up to 60–80% increases.
The cumulative contribution of urbanization for different quantiles of discharges in each season for (a) spring, (b) winter, (c) summer, and (d) fall.
The cumulative contribution of urbanization for different quantiles of discharges in each season for (a) spring, (b) winter, (c) summer, and (d) fall.
In general, the impact of urbanization was different for different discharge quantiles. It is the main cause of the increasing discharge in spring and winter. The amount of urban living sewage has increased significantly in the past decades because of the rapid increase of urban population in this metropolitan river basin, leading to increasing discharge. For summer and fall, the impervious area plays the largest role in lower discharge quantiles (below the 30th percentile in summer and below the 50th percentile in fall). Extreme floods usually occurred in the summer in Beijing, and precipitation was the major driver, with no significant role played by urbanization.
DISCUSSIONS AND CONCLUSIONS
In this study, we developed a statistical framework to attribute the changes in a large set of discharge quantiles (from low to high flows). We focused on the Wenyu river basin, within the Beijing metropolitan area, and performed analyses at the seasonal scale from 1985 to 2016. Precipitation, antecedent wetness, temperature, impervious areas, recycled water amount, and groundwater levels were considered as potential predictors to identify the main drivers for different discharge quantiles in this highly urbanized watershed. Our conclusions can be summarized as follows.
The GAMLSS model can capture the variability of streamflow in the urban watershed. The discharge data are better described in spring, summer, and winter, with an average Spearman correlation coefficient between observations and the median of the fitted model of 0.84, 0.79, and 0.81, respectively; the performance decreases in the fall, with the average Spearman correlation coefficient of 0.40. Summer is the core of the flood season in Beijing, and the model captures it well. In spring and winter, the recycled water amount represents the key driver of discharge, in part due to the low precipitation during those seasons. The performance of the statistical models decreased in the fall, consistent with the prior work (Slater & Villarini 2017). We evaluated the robustness of our results using leave-one-out cross-validation, with the results of these analyses supporting the goodness of fit of our models.
In terms of predictors, the recycled water amount was the main factor for streamflow in spring and winter, and the antecedent wetness has an impact on the higher quantiles in spring. In summer, both precipitation and impervious area are the main driving factors for discharge quantiles below the median, and only precipitation is a relevant factor for the upper half of the summer discharge distribution. For fall, precipitation, impervious area, and antecedent wetness are the three driving factors, with the impervious area (precipitation and antecedent wetness) being more relevant for the lower (upper) half of the distribution. The impact of urbanization was different for different discharge quantiles. Overall, we found the strongest impact of urbanization to be in winter and spring, and for low and medium quantiles. The role of precipitation is the strongest in the summer, and it increases as we move toward the upper tail of the discharge distribution. Recycled water, on the other hand, tends to play a more dominant role in winter and spring.
This study has proposed and developed a statistical framework for the attribution of discharge (from low to high flow) in a basin that has been experiencing extreme urbanization over the past few decades. In this study, the comprehensive potential factors have been considered, and the impact of climate change and urbanization on hydrology in the very complex urban watershed have been assessed. The results showed that climate change and urbanization have important effects on urban watershed hydrology, and it is consistent with the previous studies (Damodaram et al. 2010; Jacobson 2011; Qin et al. 2013). The advantage of this study lies in a more detailed evaluation of the impact of climate change and urbanization, and our findings suggest that the biggest impact of impervious areas is not on the extreme floods, but rather on the lower part of the distribution. Therefore, it is important to allow water to infiltrate into the ground. Alternative approaches beyond the reduction of impervious areas should be considered to effectively mitigate the impacts of floods on the urban areas. This is an important issue to deal with, given that a large number of people live in urban areas, increasing the societal and economic impacts associated with flooding. These conclusions have important significance for the Sponge city's construction and urban water management in China.
The limitation of this study is caused not only by the fact that there are parameter uncertainties and limitations of the statistical models, but also that only one watershed has been considered due to the shortage of data. The representativeness of the Wenyu river basin is still limited, although this river basin has experienced a very obvious process of urbanization.
AUTHOR CONTRIBUTIONS
M.F.R. and G.V. designed the experiments and analyzed the data; M.F.R. performed the experiments; Z.X.X. and B.P. worked on the paper and interpreted the results; Y.C.W. and L.G.D. provided the observation data.
FUNDING
This work was financially supported by the National Key R&D Program of China (2017YFC1502701) and the National Natural Science Foundation of China (51879008). The authors wish to thank the Beijing Hydrology Center for providing the observation data.
CONFLICT OF INTEREST STATEMENT
The authors declare no conflicts of interest.
DATA AVAILABILITY STATEMENT
All relevant data are available from an online repository or repositories. The essential urban land use categories in China were downloaded from: http://data.ess.tsinghua.edu.cn. The temperature data were downloaded from: http://westdc.westgis.ac.cn/data/7a35329c-c53f-4267-aa07-e0037d913a21).