Detecting Water Diversion Fingerprints in the Danjiangkou Reservoir from Satellite Gravimetry and Altimetry Data

The Danjiangkou Reservoir (DJKR) is the freshwater source for the Middle Route of the South-to-North Water Diversion Project in China, and its water level and storage changes are important for water resource management. To maximize the potential capacity of the Gravity Recovery and Climate Experiment (GRACE) mission, an improved Lagrange multiplier method (ILMM) is first proposed to detect terrestrial water storage anomalies (TWSA) in the small-scale basin (DJKR). Moreover, for the first time, water diversion fingerprints are proposed to analyze the spatiotemporal pattern of the TWSA in the DJKR. The results indicate that the increased water level and storage signals due to the DJKR impoundment in 2014 can be effectively detected by using the ILMM, and they agree well with the results from altimetry and in situ data. Additionally, the water diversion fingerprints due to the DJKR impoundment are inferred, and describe the progression of spatiotemporal variability in water storage. The results show that water storage decreased in the upper Hanjiang River and increased in the DJKR as well as to the east of it during the period 2013–2015. Our research provides a scientific decision-making basis for monitoring the water resources of the DJKR and managing the South-to-North Water Diversion Project.


Introduction
Surface water that is stored in lakes, reservoirs, and rivers plays a large role in agricultural irrigation, aquaculture, hydroelectric power, disaster prevention and mitigation, human life, and industrial activities [1]. However, surface water resources are unevenly distributed in China, with abundant freshwater in the south and water scarcity in the north. To address this issue, the Chinese government decided to implement the South-to-North Water Diversion Project (SNWDP) in 2002, which includes three routes: East, Middle, and West. In particular, the Middle Route of the SNWDP (MRSNWDP) was devised to divert freshwater from the Danjiangkou Reservoir (DJKR) to Henan and Hebei provinces and the Beijing and Tianjin municipalities in China [2]. The Hanjiang River is the longest tributary of the Yangtze River [3], and approximately 70% of the freshwater is diverted to the DJKR. The DJKR is The main contributions and novelties of this research are as follows: (1) the relationship between changes in water level from in situ observations and the areas from Landsat images is determined in the DJKR; (2) a new method of inferred leakage correlation and ILMM is proposed to detect the TWSA in the DJKR; (3) the signal from the DJKR impoundment in 2014 is captured by using GRACE and multi-source altimetry missions for the first time; and (4) the water diversion fingerprints due to DJKR impoundment are derived and verified for the first time.

Study Area
The DJKR (32 • 36 -33 • 48 N, 110 • 59 -111 • 49 E) is located upstream of the Hanjiang River at the junction between the Hanjiang River and Danjiang River and forms a "V" shape. It includes the Hanjiang and Danjiang Reservoir areas, and distributes water to Hubei and Henan provinces in China (Figure 1). With a total drainage area of approximately 95,000 km 2 and an average water storage volume of approximately 39.48 billion cubic meters, the DJKR is the largest artificial freshwater lake in Asia [23]. In addition, the DJKR is a top-grade multi-purpose reservoir with flood control, electricity generation, navigation, and agricultural irrigation functions. In this study, data from altimetry missions were used to monitor the water level change in the DJKR (the red rectangular region in Figure 1), and GRACE mission data were used to detect the TWSA in the total drainage basin. The main contributions and novelties of this research are as follows: (1) the relationship between changes in water level from in situ observations and the areas from Landsat images is determined in the DJKR; (2) a new method of inferred leakage correlation and ILMM is proposed to detect the TWSA in the DJKR; (3) the signal from the DJKR impoundment in 2014 is captured by using GRACE and multi-source altimetry missions for the first time; and (4) the water diversion fingerprints due to DJKR impoundment are derived and verified for the first time.

Study Area
The DJKR (32°36′-33°48′ N, 110°59′-111°49′ E) is located upstream of the Hanjiang River at the junction between the Hanjiang River and Danjiang River and forms a "V" shape. It includes the Hanjiang and Danjiang Reservoir areas, and distributes water to Hubei and Henan provinces in China (Figure 1). With a total drainage area of approximately 95,000 km 2 and an average water storage volume of approximately 39.48 billion cubic meters, the DJKR is the largest artificial freshwater lake in Asia [23]. In addition, the DJKR is a top-grade multi-purpose reservoir with flood control, electricity generation, navigation, and agricultural irrigation functions. In this study, data from altimetry missions were used to monitor the water level change in the DJKR (the red rectangular region in Figure 1), and GRACE mission data were used to detect the TWSA in the total drainage basin.   The location of the DJKR has a subtropical monsoon climate, which has remarkable transitional climatic characteristics. The annual mean temperature and precipitation are 14.4-15.7 • C and 800-1000 mm, respectively, and 80% of the rainfall occurs between May and October [24]. The upper Hanjiang River is the source of 70% of the runoff of the DJKR, and the remaining 30% comes from the Danjiang River, which is the longest tributary of the Hanjiang River. The DJKR has a complex topography, with mountains and hills forming approximately 97% of the area.
To implement the MRSNWDP, the height of the Danjiangkou Dam (DJKD) was increased from 162.0 m to 176.7 m, and its normal water level and storage capacity were increased from 157.0 m to 170.0 m and from 1.7 × 10 10 m 3 to 2.9 × 10 10 m 3 , respectively [25]. The DJKR was impounded in October 2013, and began diverting water to Beijing, Tianjin, Henan, and Hebei in December 2014.

Data
The area changes in the DJKR were determined by using Landsat 7 images, and the TWSA were detected by using GRACE data. The Global Land Data Assimilation System (GLDAS) [26], Climate Prediction Center (CPC) [27], WaterGAP Global Hydrology Model (WGHM) [28], and Community Land Model (CLM4.0) [29] were used to determine the scale factors [30]. Altimetry data were used to determine the water-level change in the DJKR, and in situ water-level data were used to establish the relationship with the area of the DJKR and verify the water-level changes identified from altimetry data and the TWSA detected from GRACE data. In addition, precipitation, evapotranspiration, GLDAS, and CLM4.0 were used to support the GRACE-TWSA-based detection of water diversion fingerprints in the DJKR.

Satellite Images
Data from the Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), Operational Land Imager (OLI), and Moderate Resolution Imaging Spectroradiometer (MODIS) were successfully applied to estimate the area of the lakes, reservoirs, and rivers [31]. Although MODIS data a have high temporal resolution, their spatial resolution is approximately 250-1000 m, and the narrowest width of the DJKR is approximately 300 m. The spatial and temporal resolutions of Landsat multispectral images (TM, ETM+, and OLI) are 30 m and 16 days, respectively. Landsat spectral bands 1-7 were used. Here, Landsat 7 multispectral images (path/row 138/037, TM and OLI) captured during the period of 2011-2014 were collected from the United States (U.S.) Geological Survey (USGS) server [32].

GRACE Data
The ILMM requires an error in the satellite measurement, which is the formal error (Wahr et al., 2006) of the GRACE-variable gravity field model. It was not provided by the Center for Space Research (CSR) [33,34]. To compare with the results of the MASCON from the Jet Propulsion Laboratory (JPL), GRACE data were obtained from the JPL Release 06 from April 2002 to March 2016. The C 20 term of the GRACE time-variable gravity field model was determined from satellite laser ranging observational data [35]; the degree-one harmonic coefficients (Earth's geocenter) were estimated from Swenson et al. [36], and correction for the glacial isostatic adjustment (GIA) was done following A. Geruo et al. [37].
TWSA were also obtained from the GRGS (Groupe de Recherche de Géodésie Spatiale, etc. [38]) with a degree and order of 80, which applied DDK5 filtering [39] to the data from the CSR, GeoForschungsZentrum Postdam (GFZ), and JPL. The method used to estimate the TWSA from the above products is described in Wahr et al. [40]. The results from the JPL fan-shaped (FAN) filter [41], JPL mass concentration (MASCON) [42], and CSR MASCON [43] were also compared.

Laser Altimetry Data
The Ice, Cloud, and land Elevation Satellite (ICESat) mission, which is part of the National Aeronautics and Space Administration (NASA) Earth Observing System (EOS), was launched in January 2003 from Vandenberg Air Force Base; it ended in August 2010. The sole instrument on board ICESat was the Geoscience Laser Altimeter System (GLAS), which is a laser altimeter that provides a high-precision dataset. As the ICESat orbited, the GLAS produced a series of laser spots of approximately 70 m in diameter that were separated by nearly 170 m along the spacecraft ground track [44]. From 2003 to 2009, ICESat/GLAS observed 19 campaigns with a ground track cycle of 8 or 91 days. Here, the data from GLA14 (Version 34) [45] were used.

Radar Altimetry Data
In this study, we used the 20-Hz Environmental Satellite (Envisat) Geophysical Data Record (GDR) data product in version 2.1 [46] and the 40 Hz Satellite with ARgos and ALtiKa (SARAL) GDR data product in the Calibration/Validation (Cal/Val) phase. Since the quality of the Envisat data might not be good at the mission-commissioning phase, this study applied data that were obtained between January 2003 and September 2010 with a repeat cycle of 35 days. The SARAL reached the historical orbit of Envisat at the end of October 2013 [47] and transitioned to a drifting orbit in July 2016. Thus, for the Ka band, we chose altimetry data that were collected between November 2013 and March 2016 with a repeat cycle of 35 days.

Model Data
(1) GLDAS: The Global Land Data Assimilation System (GLDAS) [26] includes four monthly land surface models: the CLM (the Community Land Model), Mosaic (MOS), VIC (the Variable Infiltration Capacity), and NOAH (National Centers for Environmental Prediction/Oregon State University/Air Force/Hydrologic Research Lab Model), which do not include surface water storage in lakes, reservoirs, and river channels. The average of the four monthly land surface models with soil moisture states in a 1 • × 1 • grid was used to estimate the water storage in the top 2 m of the soil layer. GLDAS include soil moisture, snow, and vegetation canopy storage, excluding surface water storage in lakes, reservoirs, and river channels. Therefore, GLDAS was subtracted by the GRACE data to infer the total change in surface water storage.
(2) CPC: The Climate Prediction Center (CPC) [27] with soil moisture states in 0.5 • × 0.5 • with monthly time steps were used to estimate the water storage.
(4) CLM4.0: The National Center for Atmospheric Research (NCAR) Community Land Model (CLM4.0) [29] simulates the partitioning of mass and energy from the atmosphere, the redistribution of mass and energy within the land surface, and the export of freshwater and heat to the oceans. The spatial and temporal resolutions of CLM4.0 are 0.9 • × 1.25 • and monthly, respectively. Components of terrestrial water storage output by the CLM4.0 include soil moisture, snow, vegetation canopy storage, channel storage in rivers, and climate-driven change; human activities are excluded. GRACE-TWSA is the total water storage change. Therefore, human-induced TWSA can be determined by GRACE subtracted by CLM4.0.

Precipitation, Evapotranspiration, and Water-Level Data
The monthly precipitation data are from the Tropical Rainfall Measuring Mission (TRMM), which has spatial and temporal resolutions of 0.25 • × 0.25 • and monthly, respectively [48], as well as the Danjiangkou weather station (Danjiangkou, China) (China Meteorological Data Service Center [49]) and the Xiaotao hydrological station (Xiantao, China) (Chang Jiang Water Resources Commission of the Ministry of Water Resources [50]). Here, evapotranspiration data are from different products, i.e., monthly evapotranspiration products from GLDAS, Moderate Resolution Imaging Spectroradiometer (MODIS), and the Danjiangkou weather station. The GLDAS data integrate satellite data and Land Surface Model (LSM) data to generate a global distribution of land surface states (e.g., evapotranspiration). Evapotranspiration data from MODIS with spatial and temporal resolutions of 1 • × 1 • and monthly are based on the Penman-Monteith method, in which remote sensing data and meteorological observations are combined [51]. The in situ water level data at the DJKD observation station are from the Hydrology and Water Resources Survey Bureau of China [52].

Deriving the DJKR Area Change from Landsat 7
The steps for capturing the area of the DJKR from Landsat 7 images are as follows: (1) preprocess the original images via radiation correction, geometric correction, and atmospheric correction [3], transform them to the Universal Transverse Mercator (UTM) projection, and obtain the reflectance data; (2) extract the water body information by using the modified normalized difference water index (MNDWI) [53]; (3) determine the threshold segmentation by using the Otsu method [54]; and (4) remove the water data that are outside of the reservoir, count the total number of pixels, multiply by the area of the pixels, and obtain the area of the DJKR as the final output.
Changes in water volume (i.e., changes in mass) cannot be estimated by only using water-level data. Thus, the relationship between the change in area and water level should be obtained. The most widely used method is to combine the area retrieved from satellite images and the water level from multi-source altimetry missions or from in situ observations.

Estimating the Water-Level Change in the DJKR from Altimetry Data
In Frappart et al.'s study for European ENVIronmental SATellite (ENVISAT) validation over the Amazon basin [55], the elevation change retrieved by the ice-1 algorithm [56] agreed best with in situ gauge data from the inland water level. Here, the altimetry data (ENVISAT and SARAL) were retrieved by the ice-1 algorithm [57] to perform further analysis. To retrieve the elevation change in the water level in the DJKR, geophysical range corrections, including (European Centre for Medium-Range Weather Forecasts) modeled wet and dry troposphere delays, the ionosphere delay, solid Earth tides, and pole tides were applied [56]. Due to the potential contamination by radar echoes from land reflectors and the locking loss of onboard trackers near the coasts, only data near the center of the lake were used. The Envisat and SARAL tracks are shown in Figure 1.
The ICESat laser altimeter measured elevations in 18 cycles from 2003 to 2009, with each cycle lasting 12-15 days. Previous studies, such as that by Nicolas et al. [58], have shown that the mean accuracy of the ICESat-derived elevations over flat deserts is ±15 cm. Therefore, ICESat has sufficient accuracy to determine changes in the lake and reservoir water levels. ICESat observations were edited using the following empirical procedure: (1) remove the ICESat elevations that exceed the minimum and maximum Shuttle Radar Terrestrial Mission (SRTM) elevations within the study region; (2) remove the elevations at two consecutive along-track footprints that have a difference exceeding 15 m (the rationale for this approach is that allowing a maximum lake/reservoir slope of 5 • for two neighboring along-track footprints spaced at approximately 172 m results in a maximum elevation difference of approximately 15 m); and (3) determine a smoothed value h n at a given point using n neighboring points (n = 50 in this study). Along a sufficiently long ground track, the differences between the original and smoothed heights were computed to determine the standard deviation (δ n ).
Then, if the original height (h i ) fit the condition h i − h n > 3δ n , then h i was removed. The final ICESat tracks are shown in Figure 1. Moreover, the biases between the different altimetry missions (ENVISAT, SARAL, and ICESat) were corrected by eliminating systematic errors, and the surface water-level results from different altimetry missions in different places were the regional average values in the DJKR water bodies.

Improved Lagrange Multiplier Method to Infer TWSA in the DJKR
Due to the flight characteristic factors of the GRACE task and the limits of the background and current mathematical models, the GRACE-TWSA must be filtered to eliminate the north-south striping error [40]: where λ is the geocenter longitude; θ is the geocenter latitude; a is the semi-major axis of the reference ellipsoid; l, m are the degree and order; L is the highest degree; ρ ave is the average Earth density (5517 kg/m 3 ); ρ w is the water density (1000 kg/m 3 ); k l is the Love number [59]; P lm is the completely normalized Legendre association function [60]; ∆C lm and ∆S lm are the residual spherical harmonic coefficient after subtracting the long-term mean of Earth's gravity field, which is about the normal Earth gravity field; and W l is the smoothing coefficients of the Gaussian filter [40,61].
Smoothing methods can also cause signal attenuation and leakage. Therefore, the scale factor (S), leakage (l m ), and bias (b m ) were obtained on the basis of the hydrological models from Landerer and Swenson [30], Longuevergne et al. [62], and Klees et al. [63]: where θ is the colatitude, λ is the geocentric longitude, dΩ = sin θdθdλ, M(θ, λ) is the value from the models, and R(θ, λ) is a function of the study region: A bar over M(θ, λ) and R(θ, λ) indicates their first filter.
To determine the leakage, the outside region of the study area is defined as: The above methods are highly dependent on the hydrological model, which can cause scaling factors to differ considerably, especially in semi-arid, arid, and over-irrigated areas [64]. GRACE-TWSA can be considered a hydrological signal, but this signal differs from hydrological models. If the two types of data are fitted, signal aliasing can occur. Therefore, a scale factor that is independent of the hydrological model was proposed by Dutt Vishwakarma et al. [65], so it is called the DV method (DVM) herein: The TWSA (Equation (1)) in the spatial domain can be expressed as: The first filter is: The abbreviation expression is: where l is the leakage and is calculated by: With f = sF: Thus, the leakage from Equation (11) can be transformed by using Equation (7): Here, we used f = sF to derive a new expression of leakage l (Equation (13)). Therefore, from Equations (7), (10), and (13), we can determine the surface mass change expression, which does not need the hydrological model to restore the signal and correct the leakage: According to Equations (7), (9), (13), and (14), the key to improving the accuracy of mass change is the determination of the optimal smooth kernel for the study area.
To estimate the regional surface mass changes using GRACE data, Swenson et al. [66] proposed the Lagrange multiplier method (LMM), which reduced the impact of GRACE observation errors.
The key step in this method is the determination of the Lagrange parameter λ, which is obtained by the minimum satellite measurement error and signal leakage error, and only the formal error is needed.
The LMM is a spatial smoothing method [66] that requires the restoration and correction of the signal. Here, in the improved LMM (ILMM), the signals were restored and corrected after applying the LMM. The process of the ILMM comprises the following steps: (1) the TWSA are first determined by using the LMM; (2) the restored signal and corrected signal are determined by using Equations (7) and (13); and (3) the TWSA results are finally determined from Equation (14). More details on the Lagrange multiplier method are in Appendix A, Swenson et al. [66], and Chao et al. [14].
From the above descriptions, the surface mass change can be inferred by different methods.

Inferring Total Surface Water Storage and Human-Induced Surface Water Storage Anomaly
The total TWSA from GRACE can be disaggregated into soil moisture anomaly (SMA), surface water storage anomaly (SWSA), and groundwater storage anomaly (GWSA) [6]: Here, our study area is the DJKR; the total surface water storage anomaly (TSWSA) is defined as the sum of the SWSA and GWSA, and the human-induced surface water storage anomaly (HSWSA) can be represented by the difference between the TSWSA and climate-driven surface water storage anomaly (CSWSA). Arranging Equation (18) [67,68]: SMA can be obtained from GLDAS or CLM4.0, and CSWSA can be obtained from CLM4.0 [68,69]. Here, TSWSA is determined by combining GRACE and GLDAS, and CSWSA is determined by combining GRACE and CLM4.0, which is the same as approach as that in the studies of Voss et al. [67] and Joodaki et al. [68].

The Area Change in the DJKR from Landsat
The surface area of the DJKR expansion or reduction responds to a water level increase or decrease, respectively, which is proven in Figure 2a. The surface area and water-level change can be affected by climate variability and human activities. The natural (climate-driven) variability is mainly from climate-derived precipitation and evapotranspiration changes, which result in seasonal changes in surface area and water level. Human-induced changes include irrigation, drinking water, and power generation, which lead to the consumption of water resources, and artificial impoundment increases water resources. From Figure 2a, the surface area and water level during the period of 2012-2013 both declined without a seasonal effect, possibly as a result of human activities (such as power generation). The surface area and water level both increased in 2014 because of artificial impoundment. Figure 2b shows a good linear correlation between the changes in the DJKR area and the water level, with an R 2 value of 0.99. Once an optimal relationship between the changes in area and water level is determined, it can be used to fill in and predict missing data on the area or water level, thus generating a longer and more efficient time series that can be used to verify the results from the GRACE data. This approach has useful applications in water resource management, flood control, and environmental monitoring of the MRSNWDP.

Results of the Scale Factor and Leakage
The SFM depends on the hydrological model and the smoothing method, but the DVM is independent of the hydrological model and only depends on the optimal average kernel. Here, the scale factors were obtained by different hydrological models and smoothing methods from the SFM and different smoothing methods from the DVM (Table 1). As shown in Table 1, the scale factors are greater when the smoothed radius increases, which shows that more signals are attenuated for larger smoothed radii. The scale factors from different hydrological models are basically consistent with each other. The scale factors from the WGHM and CLM4.0 are larger than those from the CPC and GLDAS. The scale factor from the smoothed LMM is smaller than that from the 300-km Gauss filter by using the SFM; thus, the LMM can retain a more effective signal than that from the 300-km Gauss filter.
The scale factor results from the DVM show that as the smoothed radius increases, the scale factors also quickly increase. According to Equation (7), the scale factors largely depend on the study area. When the area is smaller, more and quicker signals are attenuated. The scale factor determined by using the DVM with a smoothed 200-km Gauss and the scale factor from the LMM are consistent with each other and closer to the results of the SFM than different hydrological models. This result shows that the key to scale factors from the DVM is to establish an approach that obtains the optimal regional average kernel, and the optimal kernel based on the LMM method can effectively preserve the signal and improve the spatial resolution.
The red line in Figure 3 represents the results of the DJKR leakage obtained from the LMM and Equation (13), and the other line represents data from the SFM (Equation (3)) with the smoothed  The SFM depends on the hydrological model and the smoothing method, but the DVM is independent of the hydrological model and only depends on the optimal average kernel. Here, the scale factors were obtained by different hydrological models and smoothing methods from the SFM and different smoothing methods from the DVM (Table 1). As shown in Table 1, the scale factors are greater when the smoothed radius increases, which shows that more signals are attenuated for larger smoothed radii. The scale factors from different hydrological models are basically consistent with each other. The scale factors from the WGHM and CLM4.0 are larger than those from the CPC and GLDAS. The scale factor from the smoothed LMM is smaller than that from the 300-km Gauss filter by using the SFM; thus, the LMM can retain a more effective signal than that from the 300-km Gauss filter.
The scale factor results from the DVM show that as the smoothed radius increases, the scale factors also quickly increase. According to Equation (7), the scale factors largely depend on the study area. When the area is smaller, more and quicker signals are attenuated. The scale factor determined by using the DVM with a smoothed 200-km Gauss and the scale factor from the LMM are consistent with each other and closer to the results of the SFM than different hydrological models. This result shows that the key to scale factors from the DVM is to establish an approach that obtains the optimal regional average kernel, and the optimal kernel based on the LMM method can effectively preserve the signal and improve the spatial resolution.
The red line in Figure 3 represents the results of the DJKR leakage obtained from the LMM and Equation (13), and the other line represents data from the SFM (Equation (3)) with the smoothed 300-km Gauss and WGHM model. As shown in Figure 3, because the area of the DJKR is small, the leakage error is large. Therefore, correcting the signal leakage is crucial for determining the TWSA in small-scale basins. 300-km Gauss and WGHM model. As shown in Figure 3, because the area of the DJKR is small, the leakage error is large. Therefore, correcting the signal leakage is crucial for determining the TWSA in small-scale basins.    Figure 4 shows the time series of the TWSA in the DJKR according to different data sources and different post-processing methods. The results in Figure 4b are from different filtering methods (such as the FAN and DDK5 filters) by applying different data (such as CSR, JPL, and GFZ) and different MASCON products. As shown in Figure 4, the amplitude is slightly different, but the periodic signals are basically consistent with each other. The periodic signals of the TWSA from the ILMM are more obvious than those from other methods. The difference in the TWSA amplitude between January 2014 and November 2014 according to the ILMM is about 30 cm, which is the largest of the results of different methods. 300-km Gauss and WGHM model. As shown in Figure 3, because the area of the DJKR is small, the leakage error is large. Therefore, correcting the signal leakage is crucial for determining the TWSA in small-scale basins.  Figure 4 shows the time series of the TWSA in the DJKR according to different data sources and different post-processing methods. The results in Figure 4b are from different filtering methods (such as the FAN and DDK5 filters) by applying different data (such as CSR, JPL, and GFZ) and different MASCON products. As shown in Figure 4, the amplitude is slightly different, but the periodic signals are basically consistent with each other. The periodic signals of the TWSA from the ILMM are more obvious than those from other methods. The difference in the TWSA amplitude between January 2014 and November 2014 according to the ILMM is about 30 cm, which is the largest of the results of different methods.

DJKR Water Level and Storage Changes from GRACE, Altimetry, and Hydroclimatic Data
As shown in Figure 5a, the periodic signals between water levels from the altimetry mission (ICESat, ENVISAT, and SARAL/Altika), TWSA from in situ observations of changes in the water level and area, and GRACE-TWSA agree well with each other. Figure 5a reveals that the 2006-2009 water-level data from ICESat and ENVISAT are in good agreement. The objective of the ICESat mission was to monitor mass changes in the ice sheets at the poles, which led to sparse laser footprints at middle and low latitudes. The signals of the water level from the ICESat mission were not well captured, especially between 2007-2008 (Figure 5a). Between 2012-2014, the water level declined as a result of human activities, such as power generation. This can be immediately captured by in situ data. Human activities cannot be reflected in real time because of the limitation of GRACE spatiotemporal resolution. Therefore, there are some peak TWSA values and human-induced TWSA, but no signals are detected in the in situ water level. The amplitudes between GRACE-TWSA and the water storage based on water level change and area data are slightly different.

DJKR Water Level and Storage Changes from GRACE, Altimetry, and Hydroclimatic Data
As shown in Figure 5a, the periodic signals between water levels from the altimetry mission (ICESat, ENVISAT, and SARAL/Altika), TWSA from in situ observations of changes in the water level and area, and GRACE-TWSA agree well with each other. Figure 5a reveals that the 2006-2009 water-level data from ICESat and ENVISAT are in good agreement. The objective of the ICESat mission was to monitor mass changes in the ice sheets at the poles, which led to sparse laser footprints at middle and low latitudes. The signals of the water level from the ICESat mission were not well captured, especially between 2007-2008 (Figure 5a). Between 2012-2014, the water level declined as a result of human activities, such as power generation. This can be immediately captured by in situ data. Human activities cannot be reflected in real time because of the limitation of GRACE spatiotemporal resolution. Therefore, there are some peak TWSA values and human-induced TWSA, but no signals are detected in the in situ water level. The amplitudes between GRACE-TWSA and the water storage based on water level change and area data are slightly different.  The signals for the annual, semi-annual, and trend were removed by the least-squares fitting method, and the residual is called the non-seasonal signal [15], which is shown in Figure 5b. Figure 5b shows the non-seasonal signal of GRACE-TWSA, human-induced TWSA, the SARAL/Altika altimetry, and in situ (DJKD) water-level changes. From Figure 5a, the human-induced TWSA are weaker than those from GRACE-TWSA before 2011, but they were the same after 2011. The non-seasonal signal of GRACE-TWSA and human-induced TWSA from Figure 5b show the same results, indicating that the TWSA in the DJKR were influenced mostly by human activities after 2011 [50]. Figure 5 also shows a steep water level increase due to the impoundment in 2014, which is detected in the GRACE data and Altika altimetry data. Moreover, they agree well with the in situ data. According to the above results, the DJKR impoundment signal in 2014 was captured by GRACE, altimetry missions, and in situ data, indicating that GRACE and altimetry can be effectively used to monitor the surface water change in a small reservoir. Figure 6a shows the change in precipitation in the Hanjiang River basin (HRB) from the TRMM, the Xiantao hydrological control station, and the Danjiangkou city weather station. The Xiantao hydrological control station exhibits hydrological information for the entire HRB. The results show that the HRB precipitation data from the TRMM and the Xiantao hydrological control station are basically the same, and the correlation coefficients between them are over 0.8. They agree well with the precipitation data from the Danjiangkou city weather station, except for amplitude differences in the summer; thus, the data from the TRMM are sufficient to quantify the HRB precipitation. From the strong relationship between the monthly precipitations data from the TRMM and the Danjiangkou city weather station (e.g., high correlation coefficients), the rainfall situation for the entire HRB and DJKR are basically the same. Figure 6b shows the evapotranspiration in the Hanjiang River basin and the DJKR from the four GLDAS models, MODIS, and the Danjiangkou city weather station, and it indicates that evapotranspiration has a large uncertainty. As shown in Figure 6, the water in the HRB was sufficient. Precipitation and evapotranspiration were steady, without heavy rainfall and/or low evaporation in 2014.
The signals for the annual, semi-annual, and trend were removed by the least-squares fitting method, and the residual is called the non-seasonal signal [15], which is shown in Figure 5b. Figure  5b shows the non-seasonal signal of GRACE-TWSA, human-induced TWSA, the SARAL/Altika altimetry, and in situ (DJKD) water-level changes. From Figure 5a, the human-induced TWSA are weaker than those from GRACE-TWSA before 2011, but they were the same after 2011. The non-seasonal signal of GRACE-TWSA and human-induced TWSA from Figure 5b show the same results, indicating that the TWSA in the DJKR were influenced mostly by human activities after 2011 [50]. Figure 5 also shows a steep water level increase due to the impoundment in 2014, which is detected in the GRACE data and Altika altimetry data. Moreover, they agree well with the in situ data. According to the above results, the DJKR impoundment signal in 2014 was captured by GRACE, altimetry missions, and in situ data, indicating that GRACE and altimetry can be effectively used to monitor the surface water change in a small reservoir. Figure 6a shows the change in precipitation in the Hanjiang River basin (HRB) from the TRMM, the Xiantao hydrological control station, and the Danjiangkou city weather station. The Xiantao hydrological control station exhibits hydrological information for the entire HRB. The results show that the HRB precipitation data from the TRMM and the Xiantao hydrological control station are basically the same, and the correlation coefficients between them are over 0.8. They agree well with the precipitation data from the Danjiangkou city weather station, except for amplitude differences in the summer; thus, the data from the TRMM are sufficient to quantify the HRB precipitation. From the strong relationship between the monthly precipitations data from the TRMM and the Danjiangkou city weather station (e.g., high correlation coefficients), the rainfall situation for the entire HRB and DJKR are basically the same. Figure 6b shows the evapotranspiration in the Hanjiang River basin and the DJKR from the four GLDAS models, MODIS, and the Danjiangkou city weather station, and it indicates that evapotranspiration has a large uncertainty. As shown in Figure  6, the water in the HRB was sufficient. Precipitation and evapotranspiration were steady, without heavy rainfall and/or low evaporation in 2014.  The TWSA from GRACE and CLM4.0 data, water-level change from in situ data, and precipitation from the Danjiangkou city weather station are all compared in Figure 7. A time lag between precipitation and GRACE-TWSA is revealed in this figure; in other words, when precipitation increases or decreases, GRACE-TWSA responds by increasing or decreasing after several months. The maximum correlation of 0.62 between precipitation and GRACE-TWSA is at a two-month lag ( Figure S1). The results also show that the HRB precipitation was steady without heavy rainfall in 2014, but the human-driven surface water storage change shows a steep increase. Therefore, the signal for the steep water level increase is due to impoundment by humans, which can be verified by GRACE, altimetry, and in situ data.
the DJKR from the four GLDAS models, Moderate Resolution Imaging Spectroradiometer (MODIS), and weather station.
The TWSA from GRACE and CLM4.0 data, water-level change from in situ data, and precipitation from the Danjiangkou city weather station are all compared in Figure 7. A time lag between precipitation and GRACE-TWSA is revealed in this figure; in other words, when precipitation increases or decreases, GRACE-TWSA responds by increasing or decreasing after several months. The maximum correlation of 0.62 between precipitation and GRACE-TWSA is at a two-month lag ( Figure S1). The results also show that the HRB precipitation was steady without heavy rainfall in 2014, but the human-driven surface water storage change shows a steep increase. Therefore, the signal for the steep water level increase is due to impoundment by humans, which can be verified by GRACE, altimetry, and in situ data. Figure 7. Comparison of the changes in water storage, precipitation, and water level.

Fingerprints of the MRSNWDP in the DJKR
Impoundment would be expected to affect the water storage in the HRB and DJKR. In order to investigate the fingerprints in the DJKR, an analysis was carried out on the spatiotemporal changes in precipitation from the TRMM, TSWSA from GRACE and GLDAS data, as well as HSWSA from GRACE and CLM4.0 data.
As shown in Figure 8, the water storage increased in the upper HRB, but decreased in the lower HRB. This agrees with the goal of the SNWDP, which is the diversion of water to the north from the south in China to address the increasing reduction of water in northern China. The TSWSA was fairly steady, without a significant surplus or shortage, but the HSWSA decreased in the upper HRB and increased in the lower HRB.
To analyze the water diversion fingerprints due to the DJKR impoundment, each month's precipitation, TSWSA, and HSWSA between October 2013 and December 2015 were obtained ( Figure S2). The precipitation change and TSWSA were basically steady and water sufficient; this observation indicates that the TSWSA were mainly influenced by the rainfall in this period ( Figure  S2). However, regardless of the precipitation change, the HSWSA obviously decreased in the upper HRB and increased in the DJKR and east of it, especially between October 2014 and January 2015 ( Figure 9). The water diversion fingerprints ( Figure S2) show the progress of water resources in the upper HRB as water is impounded to the DJKR. The impoundment caused water in the upper region to decrease and water from the DJKR and the east of it to increase. The above results reveal that the water diversion fingerprints due to the DJKR impoundment can be captured by the GRACE mission.

Fingerprints of the MRSNWDP in the DJKR
Impoundment would be expected to affect the water storage in the HRB and DJKR. In order to investigate the fingerprints in the DJKR, an analysis was carried out on the spatiotemporal changes in precipitation from the TRMM, TSWSA from GRACE and GLDAS data, as well as HSWSA from GRACE and CLM4.0 data.
As shown in Figure 8, the water storage increased in the upper HRB, but decreased in the lower HRB. This agrees with the goal of the SNWDP, which is the diversion of water to the north from the south in China to address the increasing reduction of water in northern China. The TSWSA was fairly steady, without a significant surplus or shortage, but the HSWSA decreased in the upper HRB and increased in the lower HRB. To analyze the water diversion fingerprints due to the DJKR impoundment, each month's precipitation, TSWSA, and HSWSA between October 2013 and December 2015 were obtained ( Figure S2). The precipitation change and TSWSA were basically steady and water sufficient; this observation indicates that the TSWSA were mainly influenced by the rainfall in this period ( Figure S2). However, regardless of the precipitation change, the HSWSA obviously decreased in the upper HRB and increased in the DJKR and east of it, especially between October 2014 and January 2015 (Figure 9). The water diversion fingerprints ( Figure S2) show the progress of water resources in the upper HRB as water is impounded to the DJKR. The impoundment caused water in the upper region to decrease and water from the DJKR and the east of it to increase. The above results reveal that the water diversion fingerprints due to the DJKR impoundment can be captured by the GRACE mission.

Verification TWSA from ILMM
In this study, the results show that ILMM can be used to effectively infer the TWSA in small basins (such as DJKR); thus, with the ILMM, there is a high potential for the full exploitation of GRACE data for hydrology research. GRACE-TWSA of the DJKR should be verified by water balance equations [20]. Unfortunately, we do not have in situ runoff data for the DJKR. Moreover,

Verification TWSA from ILMM
In this study, the results show that ILMM can be used to effectively infer the TWSA in small basins (such as DJKR); thus, with the ILMM, there is a high potential for the full exploitation of GRACE data for hydrology research. GRACE-TWSA of the DJKR should be verified by water balance equations [20]. Unfortunately, we do not have in situ runoff data for the DJKR. Moreover, evapotranspiration is difficult to estimate and has a large uncertainty (Figure 6b) [69]. Therefore, the TWSA obtained from ILMM are compared with different methods, and the result shows that the TWSA from the different methods agree with each other, but the periods and amplitudes from the ILMM are the most pronounced. The impoundment signal in 2014 was only detected by using the ILMM (Figure 4). Additionally, five methods are used to extract the TWSA from GRACE data, which can also be used to verify the accuracy of different hydrological models in different regions.
Additionally, there are no effective methods that can fully verify the results of GRACE, and this is one of the key limitations in the application of GRACE. Presently, the most effective method for verifying the results from GRACE is comparing GRACE data with in situ hydrological station data. Although GRACE's main signal can be considered hydrological data after deducing other signals from the background model, the TWSA from GRACE represent the total surface mass change, which is different from hydrological data. In addition, there are errors in the background model and hydrological data, which can also affect the results. Therefore, it is still necessary to further study an effective method in order to verify the accuracy of the TWSA detected from GRACE data.

Combination of In Situ, GRACE, and Altimetry Data to Manage Reservoir Water Resources
Changes in the water level and storage in artificial reservoirs are important components of the water cycle and water resources. With the increasing use of reservoirs, the annual and seasonal flow will change, and water cycles will be affected. There is an urgent need to study reservoir water resource management. Compared with satellite remote sensing technology, in situ data can provide observational data with higher precision, and are more appropriate for monitoring an individual water body of importance. Although the satellite is important for monitoring water storage changes, its advantage is its consecutive data and ability to cover large areas. Therefore, the potential capacity of satellites (such as GRACE and altimetry missions) to monitor water resources in small-scale basins should be established. In situ and remote sensing data can be combined to manage global and regional water resources.
Many studies have indicated that GRACE and altimetry can be used to detect and manage water resources, especially in remote regions and in transboundary countries [7]. Here, we also confirm that the combination of GRACE remote sensing, altimetry, and model data can provide another effective method for understanding the hydrological processes in small reservoirs.
The combination of satellite gravimetry and multi-source altimetry provide a unique tool for global water resource management. With the development of space technology, regional and global hydrological models that are developed by means of space-based observations will aid in the management of water resources [20] and disrupt the deadlock caused by unavailable data and management opacity. Combining hydrological remote sensing observational data, hydrological modeling, and existing high-precision hydrological station observational data can provide a wide range of high-precision freshwater maps for a particular region or the entire world. This kind of scientific research is essential for more efficient, sustainable, and cross-boundary cooperative water resource management.

Conclusions
The surface water in rivers, lakes, and reservoirs is an important part of the environment and a key strategic resource for human development. Artificial reservoirs play an important role in water supply and flood control, and may become more critical with the increasing frequency and intensity of extreme weather events. However, few studies have focused on small artificial reservoirs. Here, we propose the ILMM to improve the potential of the GRACE mission for detecting the TWSA in the DJKR. The water diversion fingerprints due to DJKR impoundment are quantified and analyzed by using GRACE, multi-source altimetry, Landsat, land surface model, precipitation, and in situ data. The main findings are as follows: (1) Changes in the water level and area in the DJKR show a good linear correlation. The relationship between the changes in area and changes in the water level is determined, and it can be used to fill in and predict missing data, verify the TWSA from GRACE data, and benefit water resource management. (2) The ILMM can improve the spatial resolution and enable the use of GRACE data to detect the TWSA in small-scale basins (such as DJKR).
(3) GRACE and altimetry missions can be effectively used to monitor human-induced surface water changes, such as the impoundment of small artificial reservoirs. (4) The precipitation change and TSWSA in the HRB are basically steady and water sufficient.
According to GRACE and CLM4.0 data, regardless of the precipitation changes, HSWSA obviously decreased in the upper HRB and increased in the DJKR and to the east of it, which indicates that these are the human-induced TWSA. The GRACE mission can capture the phenomenon, i.e., the water diversion fingerprints due to the DJKR impoundment.
The GRACE mission has made significant and unique contributions to Earth science. This mission ended in October 2017. The GRACE Follow-On (GFO) mission was launched in May 2018, and the SWARM constellation [70] has the potential to determine the time-variable gravity field and the TWSA [71]. Therefore, we can combine GRACE, GFO, and SWARM data to detect Earth's mass transport. Additionally, the altimetry missions CryoSat-2, ICESat-2, and future missions (Water Cycle Observation Mission (WCOM) and Water and Ocean Topography (SWOT)) [72] can be used to monitor water resources, and all of these missions have the potential to complement or replace in situ datasets for global or regional water resource management.
Author Contributions: N.C. initiated this study. Z.L., X.S. and Z.W. carried out the data analyses and wrote the first draft; N.C. and G.C. finalized the data interpretations and the draft; F.L. helped with some of the plots and data analyses.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Lagrange Multiplier Method
Wahr et al. [40] proposed a simple Gaussian smoothing method to determine the terrestrial water storage anomalies (TWSA); however, the method was unable to isolate a specific basin. Additional technology is required to estimate basin TWSA by using GRACE data.
The vertical integration of water storage in any average region can be written as follows [66]: where λ is the geocenter longitude, and θ is the geocenter latitude. The accurate average kernel ϑ(θ, λ) is the function of the shape of the described regional (such as a river basin): where dΩ = sin λdλdθ is the solid angle element, and the integration in a given spherical domain Ω region (spherical area) is ϑ(θ, λ). From Equations (1) and (A1), we can obtain where a is the semi-major axis of the reference ellipsoid; ρ ave is the average earth density (5517 kg/m 3 ); k l is the Love number [59]; l, m are the degree and order; ∆C lm , ∆S lm are the residual spherical harmonic coefficient after subtracted the long-term mean of the Earth gravity field; ϑ c lm and ϑ s lm are the spherical harmonic coefficients of ϑ(θ, λ): ϑ(θ, λ) = 1 4π Using W(θ, λ) to replace the accurate average kernel ϑ(θ, λ) in Equation (A1), where TWSA region is the approximate regional average, and W can be expanded to obtain the following: ], (∆δ c lm , ∆δ s lm ) is the nominal error of GRACE data [34]. The Lagrange multiplier method fixes one type of error (satellite or leakage errors) to a specific value and minimizes the other type of error (leakage or satellite errors). A type of signal leakage is defined as the ratio of the variance of accurate and approximate average kernel difference to the variance of the accurate average kernel: When λ is obtained from Equation (A14), we can plug it into Equation (A12) to obtain W c lm and W s lm . Finally, we put them into the Equation (A8), and TWSA can be inferred. This is the Lagrange multiplier method, which is the fixed satellite measurement error to minimize the signal leakage error.