A global gridded municipal water withdrawal estimation method using aggregated data and artificial neural network

Municipal water use is defined as the water used for domestic, household purposes or public services (Ritchie & Roser 2017). Municipal water withdrawal (MWW), meaning the amount of water withdrawn from surface or ground water sources for municipal use, currently accounts for 12% of the total water withdrawal globally, while this percentage varies dramatically across regions and countries, ranging from 0.45% in Somalia and 100% in Monaco (FAO 2021).

It is expected that domestic water use will increase significantly over the first half of this century (Wada et al. 2016b; Boretti & Rosa 2019). These trends have brought heavey burdens to water resources management, especially in highly urbanized areas (Wang et al. 2021c). Therefore, it is very important to estimate MWW with high spatial resolution to assess the pattern and distribution. Moreover, understanding the spatial patterns of MWW is an important step to improve water resource utilization efficiency, reduce water shortage, and development of mitigation adaptation strategy of water resources (Vandecasteele et al. 2014; Bierkens et al. 2015). Furthermore, the development of high-resolution global hydrological models requires reliable estimation of water withdrawals for validation.

However, there is a lack of high-resolution MWW data. Currently, the spatial granularity of MWW data is too coarse, mostly reported at the region (e.g., country, state, or basin) level. An important reason is that it is not easy to obtain high quality MWW data in many regions, especially for developing countries, because it is costly to measure and collect data of the water use of a particularly large number of domestic water users while the domestic water and industrial water are mixed in the water supply network. In addition, the comparability of MWW data from different countries and regions is always in doubt because of the inconsistency of statistical caliber and statistical method (Ellingson et al. 2019).

Consequently, there is a huge gap between the importance and availability of high-resolution MWW data. Fine-resolution estimation and mapping of MWW is necessary and valuable. In general, the process of redistribute coarse region level (e.g., administrative unit) census or survey data to a finer scale (e.g., pixel level) is called spatial disaggregation (spatialization or downscaling).

Currently, multiple methods have been used to produce gridded MWW. Most studies use population or population density as a proxy to disaggregate MWW. For instance, Guo et al. (2013) established the relationship model between per capita domestic water and per capita GDP, and in combination with the population distribution map, estimated the 4-kilometer grid distribution of domestic water in parts of Northwest China. Huang et al. (2018) used the population density maps as the proxy for disaggregating domestic water withdrawal from administrative division level to grid level (0.5° resolution). Europe public water withdrawal at 5 × 5 km resolution for 2006 was estimated by Vandecasteele et al. (2014) using population and tourism density as the proxy. Moreover, gridded MWW estimation and projection modules have been incorporated in many global hydrological models, such as the PCR-GLOBWB (Wada & Bierkens 2014), WaterGAP (Flörke et al. 2013), and H08 (Hanasaki et al. 2013). Wada et al. (2016a) estimated domestic water use by multiplying the number of persons in a grid cell with the country-specific per capita domestic water extraction on a global scale with a high-resolution of 0.1°. The methods applied in hydrological models are mostly driven by population numbers, GDP per capita, and per capita water use intensities (Wada et al. 2016b). A high-resolution water demand method for households and industries was constructed with a spatial resolution of 30 arc-seconds by Lips (2020) following the estimation approach of Wada et al. (2016b).

In summary, the current studies or models have four main limitations. Firstly, the resolution is still too coarse, mostly in 0.5° spatial resolution (about 50-km in the equator) (e.g., Flörke et al. 2013; Hanasaki et al. 2013; Wada & Bierkens 2014; Huang et al. 2018). Secondly, some studies created high-resolution MWW data but most at local or regional scale (e.g., Guo et al. 2013; Vandecasteele et al. 2014; Lips 2020). Thirdly, most of the studies establish the relationship between MWW and the most influencing factors in simple empirical ways relying on some empirical parameters or formulas. Lastly, existing studies mostly convert their spatial aggregation problem into a standard supervised regression problem by first aggregating the covariates within each census level.

For standard supervised regression, each training example has an individual output in training set. In the disaggregation problem, the training set consists of subsets of examples and the individual real target of each training example is unknown. Alternatively, for each subset of examples, an aggregated target is known. This framework is called aggregated output learning problem (Musicant et al. 2007). Some approaches are proposed to handle this kind of problem, such as aggregated linear model (Hernández-González et al. 2019), linear exponential model (Derval et al. 2020), aggregated Gaussian processes (Zhu et al. 2022a), and machine learning models (Musicant et al. 2007; Zhang et al. 2021). As far as we know, there is no such kind of approaches applied in MWW disaggregation.

There are many factors that can influence MWW, such as climate and environmental conditions, socioeconomic factors, technology development, and policy interventions (Wang et al. 2021c). With the development of satellite remote sensing technology, massive geospatial data, which can be used to derive covariates in MWW estimation, are available at high spatial-temporal resolutions.

This study aims to develop a framework for estimating MWW at global scale with a high resolution of 0.1°. Different models, including traditional regression model and machine learning models, are constructed and evaluated to estimate global gridded MWW using aggregated MWW data and gridded covariates. There are three new aspects in this study. Firstly, aggregated output learning models are constructed and verified. Secondly, to methodologically compare different models in estimating gridded MWW with high spatial resolution at global scale. Thirdly, to provide a more accurate gridded MWW dataset which is conducive to the global water resources management, especially to the water resources management in developing countries.

In this study, the country level and sub-national level for MWW data of the year 2015 or around 2015 were collected from various datasets (Table 1): (I) Country level MWW data from FAO AUASTAT database (FAO 2021); (II) United States state level and county level data from USGS (Dieter et al. 2018). Since many counties are very small, we aggregated counties into 279 zones based on the principle of spatial proximity; (III) China provincial level, prefecture level and basin level data are collected from the Water Resources Bulletins issued by departments at different levels; (IV) Brazil municipality level data are acquired from Brazil National Water Agency; (V) Russian district MWW data are estimated from Shiklomanov et al. (2011).

Domestic water use is influenced by both natural conditions such as climate and the availability of water and social conditions such as the income level and the habit of water use. The choice of covariates should reflect the natural and social factors that affect domestic water use. Climate factor, topographic factors, and NDVI can reflect the natural conditions, while population density, night light, fossil fuel CO₂ emissions, NO₂ density and HDI can reflect socioeconomic conditions. Our research goal is to produce gridded MWW estimation, thus we only consider gridded dataset (e.g. remote sensing and reanalysis data). Considering the availability of data, a total of 13 geospatial covariates are used, among which three are related to climate, four to topography, and six others. The detailed geospatial covariates and their respective data sources are listed in Table 2.

Climate-related covariates include precipitation, mean temperature, and net solar radiation for the year 2015. We derive annual precipitation from the Multi-Source Weighted-Ensemble Precipitation (MSWEP) dataset (version 2; 0.1° spatial resolution) (Beck et al. 2017). MSWEP ensembles a variety of gauge-, satellite-, and reanalysis-based precipitation datasets to achieve better accuracy. For mean temperature, we employ the CRU TS v4.04 dataset (Harris et al. 2020). In addition, net solar radiation was derived from the Surface Irradiances of Edition 4.0 Clouds and the Earth's Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) data product (Kato et al. 2018).

For topographic variables, elevation, slope, and standard deviations of elevation and slope are extracted from the Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010) (Danielson & Gesch 2011), which is an enhanced elevation model covering the whole global released by the U.S. Geological Survey and the National Geospatial-Intelligence Agency. GMTED2010 has three different spatial resolutions of 30, 15, and 7.5 arc-seconds, and the 30 arc-seconds one was used in this study. Firstly, slope data is produced from the elevation raster. For each 0.1 × 0.1° grid cell, means and standard deviations of elevation and slope were subsequently calculated from the elevation and slope data.

Other geospatial covariates include the Normalized Difference Vegetation Index (NDVI), population density, nighttime light index, CO₂ emissions, NO₂ vertical column density, and the Human Development Index (HDI). Mean NDVI for the year 2015 is derived from GIMMS3g dataset (Fensholt & Proud 2012). WorldPop developed peer-reviewed open and high-resolution population distributions data. In this study, population density is derived from WorldPop dataset (Lloyd et al. 2019). Nighttime light can reflect the extent and intensity of human activities and provide the potential to assess socioeconomic development. VIIRS Stray Light Corrected Nighttime Day/Night Band Composites Version 1 (Mills et al. 2013) is employed to acquire nighttime light index information. Satellite-based fossil fuel CO₂ emissions and NO₂ concentrations can reflect the anthropogenic influence on atmosphere, which may help us to estimate MWW. In this study, ODIAC Fossil Fuel Emission Dataset 2019 (Oda 2015) and Merged TM4NO2A version 2.3 (Georgoulias et al. 2019) are used to derive fossil fuel CO₂ emissions and NO₂ concentrations. HDI is often used to describe the development status of an area, and are thus included in this study. HDI for the year 2015 is retrieved from Kummu et al. (2018).

Due to the different resolutions of these gridded datasets, resampling or aggregation procedures are performed prior to the model construction and application. The gridded dataset with a lower resolution is resampled to 0.1° cell size using the bilinear method, while an aggregation procedure is utilized to generate a reduced-resolution for dataset with higher resolution than 0.1°. ArcGIS software is used to perform these spatial analysis procedures.

Data pre-processing is a crucial step that helps enhance the quality of data which directly affects the ability of our model to learn. In this study, the features are first transformed using inverse hyperbolic sine transformation and then standardized. The pre-processing of the MWW and its covariates is performed using ArcGIS software complemented with Python scripts.

The task of this study is to use information of region level aggregated observation of MWW and the spatially gridded covariates to estimate the gridded MWW.

Suppose we have k regions that the aggregated output is known, the region set is represented by ⁠. For region ⁠, the gridded inputs are represented by ⁠, where ⁠, n_r is the number of pixels in region r, and the corresponding gridded outputs are represented by ⁠.

We have data of a set of region-level MWW ⁠, where each ⁠. We denote region-level average covariates as ⁠.

Artificial neural network (ANN) is one of the most common supervised machine learning models. The excellent performance and adaptability of machine learning has demonstrated its potential in many fields (Rasouli et al. 2012; Povak et al. 2014; He et al. 2016; Pelletier et al. 2016; Lamorski et al. 2017; Yan et al. 2019). Compared with traditional approaches (including traditional statistical models and empirical models), machine learning approaches have proven to be more effective, precise, and flexible. Machine learning uses one or more algorithms to explore the relationship between responses and their related predictors, with no need to consider the explicit mathematical form of the model. Complex nonlinear relationships can be easily handled with machine learning, which may facilitate the discovery of the underlying mechanisms (Zhu et al. 2022b). Machine learning models rely on ancillary and remotely sensed data are applied in aggregated data spatialization (Musicant et al. 2007; Gervasoni et al. 2018; Qiu et al. 2019; Derval et al. 2020; Šimbera 2020). ANN is an artificial neural network with multiple hidden layers between the input and output layers (LeCun et al. 2015). Conventional ANNs are typically feedforward networks in which information moves in only a forward direction from the input layer through hidden layers and to the output layer. Recently, ANN has been widely used in various fields because of its high accuracy and ability to model complex and non-linear relationships.

The second and third types are both pixel level learning models by using aggregated output information. The general way to formulate a model can break into two parts: the relationship between gridded inputs and unknown gridded outputs, and the relationship between gridded outputs and aggregated outputs. The first part postulates a relationship (f) between the gridded inputs and outputs, even though the gridded outputs are not directly observed. The second part converts the unobserved gridded outputs to the region-level outputs using sum or weighted sum method.

In order to find θ, the loss function is globally minimized by optimization algorithms. In this study, optimization is done using the AdamW algorithm (Loshchilov & Hutter 2019). Adam, for Adaptive moment estimation, is a variation of the standard Stochastic Gradient Descent algorithm; it uses first-order gradients and estimations of the first and second moments of these gradients to regularize learning steps.

For model building, we used MWW data of 821 regions, including 279 aggregated zones in the United States, 355 prefectural level divisions in China, seven districts in Russia and other 180 countries (Figure 3). To evaluate the accuracy of the different models, 2 × 5-fold cross validation procedures were performed in this study. The procedure randomly divided the 821 regions into five folds, and each fold was then used once as the validation regions, with the rest as the training regions. This procedure repeats 2 times. The results can be used as a measurement of model generalization ability. After the models had been trained, global gridded MWW were generated using global geospatial covariates as the input to the models.

Three statistical methods comprising coefficient of determination (R²), normalized mean absolute error (NMAE), and normalized root mean square error (NRMSE) are used to evaluate model performance. Compared to absolute performance metrics (e.g., MAE and RMSE), relative model performance metrics (e.g., NMAE and NRMSE) are more intuitive and comparable regarding datasets with different scales. There are many normalization methods in common use. One can normalize by the mean, the range, the standard deviation, or the interquartile range of observations. In this study, we simply use the mean of the observations as the normalization method. Furthermore, Taylor diagram is also used to evaluate the model performance, which shows standard deviation, correlation coefficient, and root mean square error (RMSE) at the same time.

In order to evaluate the different model, the performance metrics of R², NMAE, NRMSE are calculated during the 2 × 5-fold cross validation procedure for the training regions which include 279 merged regions in the United States, 355 prefectures and cities in China, seven districts in Russia and other 180 countries from FAO AQUASTAT database. The accuracy statistics of 2 × 5-fold cross-validation for different models are given in Table 3. Note that values in bold represent the best results in terms of specific performance metric. Within the six models studied, the two better models are Lin_M with the highest R² (0.94) and the lowest NRMSE (90.51%), and NN_M with the second highest R² (0.92) and the lowest NMAE (31.09%). It can be judged from the performance indicators that model Lin_M and model NN_M have good ability to simulate the distribution of MWW globally. As previously stated we use MAE as the loss function, according to this criterion, NN_M achieved the best performance, followed by Lin_M and LinExp_M. The performance of pixel-level indirect models (Lin_M, LinExp_M, and NN_M) are found to be more optimal than other models.

US county and state level MWWs estimated from our models are compared with that from USGS. In these comparisons, 3223 counties and 53 states are taken into consideration. The results are shown in Table 4 and Figure 5(a). For both levels, the NN_M models show the best accuracy with higher R² and lower NMAE and NRMSE. For comparisons by county, the R², NMAE and NRMSE for NN_M model are 0.76, 45.37%, and 172.77%, respectively. While the state level comparison is much better, the R², NMAE and NRMSE for NN_M model are 0.99, 10.39%, and 14.82%, respectively.

According to the water resources zoning of China, there are 10 first-level water resources zones and 82 second-level water resources zones. The MWW data was acquired from Water Resources Bulletin of different basins. Due to some MWW data of the second-level water resources zones not able to be accessed, we used data of the first-level water resources zones instead. The results of the above models by China water resources zone are shown in Table 5 and Figure 5(c). The NN_M produces better estimates than other models, which achieved the highest R² (0.98) and the lowest NMAE (10.18%) and NRMSE (16.01%).

For provincial level comparison of China, we achieved similar results as other areas. The NN_M model shows the best accuracy with a lower NMAE and NRMSE. The R², NMAE and NRMSE for NN_M model are 0.98, 6.64%, and 10.57%, respectively.

Brazilian municipalities are administrative subdivisions of the Brazilian states. Brazil currently has 27 states and 5570 municipalities. The model performance by Brazil municipality is somewhat different from previous cases (Table 6 and Figure 5(e)). Here it is shown that Lin_M is the best model, although NN_M is still relatively good and acceptable.

Overall, it is renowned that the NN_M model performed in a superior manner at most regions. Lin_M also performs well and is stable. In addition, the pixel-level indirect models (Lin_M, LinExp_M, and NN_M) are found to be more optimal in estimating MWW in terms of the results presented in Figure 5 and Tables 4,⁵–6.

Through the above verification and performance comparison, we conclude that NN_M is the best model and Lin_M is the second best model. Thus, we present the estimation results of global gridded MWW only considering these two models.

Global total MWW for the year 2015 from Lin_M and NN_M are estimated at 449.81 and 476.96 km³ respectively (Table 7), which are very close to estimations from other studies. According to Wada et al. (2011) and Gleick (2012), world total MWW for the year 2000 was 453.15 km³. Global total domestic water use is estimated for the year 2010 to be approximately 450 km³ by Wada et al. (2016a), and 422 km³ by Huang et al. (2018). World Resources Institute (Gassert et al. 2014) estimated world total MWW as 534.59 km³ for the year 2010. Table 7 shows the world and continental total MWW from our models.

The excellent performance of the NN_M in different regions and scales demonstrates the predictive accuracy of the ANN model. Generally, machine learning models have better performance than traditional regression models (Povak et al. 2014; Lamorski et al. 2017; Yan et al. 2019). The results from the cross-validation show that the performance of the NN_M reaches satisfactory levels, indicating that the employed technique is suitable for estimating MWW.

Limitations of the machine learning models are also identified. The behavior of the machine learning models may be less intuitive to interpret than traditional regression models because its algorithm cannot be fully described mechanistically. To address this problem, interpretation methods such as feature importance quantification and partial dependence analysis can provide a deeper understanding of the models. Recently, a newly interpretation approach, Shapley Additive exPlanations (SHAP) proposed by Lundberg & Lee (2017), has shown promising performance in terms of its interpretability (Matin & Pradhan 2021; Wang et al. 2021a, 2021b). SHAP can perform local and global interpretability simultaneously, and it has a solid theoretical foundation compared with other methods. Unfortunately, this new approach currently does not support our aggregated output learning problem. Overfitting is another machine learning common problem; use grid search or random search method, sometimes even by trial and error, to find the best model parameters.

Neural networks are often considered as ‘black box’ as they are difficult or impossible to interpret. Here we use the permutation feature importance (PFI) to compare the relative importance of input variables. PFI is a global interpretation method which describes the average behavior of a machine learning model. PFI measures the increase in the prediction error of a model after shuffling a feature (variable) (Molnar 2018).

In this study, six different models, region-level linear regression model (RegLin), pixel-level direct linear-exponential model (LinExp), pixel-level direct ANN model (NN), pixel-level indirect linear model (Lin_M), pixel-level indirect linear-exponential model (LinExp_M), and pixel-level indirect ANN model (NN_M), are constructed and evaluated in estimating gridded MWW.

• (1)

  The pixel-level indirect models (Lin_M, LinExp_M, and NN_M) show better performance than the other two types of models.

• (2)

  NN_M and Lin_M can provide good estimation of global gridded MWW, with verification NMAE of 31.09%, 31.34% respectively by global country/region. NN_M is more accurate in MWW estimation, whereas Lin_M is more interpretable.

• (3)

  Region-level uncertainty analysis indicates the robustness of our model at the region scale. For major regions, CV of cross-validation is less than 0.2.

• (4)

  However, there is relatively large uncertainty of estimation results in the pixel level, shown by the high CV of gridded MWW for NN_M model 20 replicates test. Although it does not influence much the suitability of the methods and the reliability of the estimation results in regions with relatively high municipal water use, there is need for further study as to how this happens and how to solve the problem.

• (5)

  Future research should draw more attention on time series of MWW due to its much more application scenarios in hydrological modelling and water resources management.

This research was carried out with support from the National Key Research and Development Program of China (2021YFE0103900) and the National Natural Science Foundation of China (41901047).

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.