A time lag exists between precipitation P falling and being converted into terrestrial water. The responses of terrestrial water storage (TWS) and its individual components to P over the global scale, which are vital for understanding the interactions and mechanisms between climatic variables and hydrological components, are not well constrained. In this study, relying on land surface models, we isolate five component storage anomalies from TWS anomalies (TWSA) derived from the Gravity Recovery and Climate Experiment mission (GRACE): canopy water storage anomalies (CWSA), surface water storage anomalies (SWSA), snow water equivalent anomalies (SWEA), soil moisture storage anomalies (SMSA), and groundwater storage anomalies (GWSA). The responses of TWSA and of the individual components of TWSA to P are then evaluated over 168 global basins. The lag between TWSA and P is quantified by calculating the correlation coefficients between GRACE-based TWSA and P for different time lags, then identifying the lag (measured in months) corresponding to the maximum correlation coefficient. A multivariate regression model is used to explore the relationship between climatic and basin characteristics and the lag between TWSA and P. Results show that the spatial distribution of TWSA trend presents a similar global pattern to that of P for the period January 2004–December 2013. TWSA is positively related to P over basins but with lags of variable duration. The lags are shorter in the low- and midlatitude basins (1–2 months) than those in the high-latitude basins (6–9 months). The spatial patterns of the maximum correlations and the corresponding lags between individual components of the TWSA and P are consistent with those of the GRACE-based analysis, except for SWEA (3–8 months) and CWSA (0 months). The lags between GWSA, SMSA, and SWSA to P can be arranged as GWSA > SMSA ≥ SWSA. Regression analysis results show that the lags between TWSA and P are related to the mean temperature, mean precipitation, mean latitude, mean longitude, mean elevation, and mean slope.
Terrestrial water storage (TWS) is a vital component of the global hydrological cycle and can be divided into surface water storage (SWS), soil moisture storage (SMS), canopy water storage (CWS), snow water equivalent (SWE), and groundwater storage (GWS) (Pan et al. 2017; Strassberg et al. 2009; Tregoning et al. 2012; see appendix for a list of acronyms used in this paper). Precipitation P is a major input to terrestrial water flux (Gao et al. 2014). In the water balance equation [d(TWSA)/dt = P − ET − R, where ET is evapotranspiration and R is runoff], it is the derivative of TWS anomalies (TWSA) that relates simultaneously to P, and therefore a time lag is expected (Eagleson 1978). Xu et al. (2018) also showed that when P is converted to TWS through water distribution, there is a theoretically delayed response of TWS to P. Determining the specific lag between TWS and P and the spatial variation of this lag are crucial to understanding interactions in the climate factors and the hydrological cycle.
Recently, the relationship between P and TWS has been investigated over a number of basins. Syed et al. (2008) demonstrated that terrestrial water storage change (TWSC) based on land surface model simulations in the Global Land Data Assimilation System (GLDAS) is strongly correlated with P in low-latitude areas. Based on TWSA from the Gravity Recovery and Climate Experiment (GRACE), Mo et al. (2016) found a stronger correlation between P and TWSA in southern than in northern China. Xu et al. (2018) identified a 2-month lag in a study of GRACE-based TWSA response to P in the Three Rivers source region of the Tibetan Plateau, similar to the 1–2-month lag between GRACE-based TWSA and P found by Soni and Syed (2015) for Indian river basins. Ndehedehe et al. (2016) determined a 1–2-month lag for TWSA responses to P for West Africa basins.
Previous studies have also found that the different components of TWS may have different responses to P. The maximum lag between TWSA and P in the Three Rivers source region of Tibet is 2 months, whereas soil moisture there can respond to P without delay (Xu et al. (2018). In southern Manitoba, Canada, Chen et al. (2002) found that the responses of groundwater level to P vary among different wells, with time lags ranging from 3 months to 3–4 years. Lorenzo-Lacruz et al. (2017) analyzed the relationship between a standardized groundwater index (SGI) and a 48-month standardized precipitation index (SPI) to investigate the relationship between groundwater storage and P variability in the Mediterranean region. They identified the responses of groundwater storage anomalies (GWSA) to P for different aquifers, including lags of 6 months for the Can Bajoca aquifer, 9 months for the Massanella red aquifer, and 46 months for the Estremera aquifer. Liesch and Ohmer (2016) established that the lags for both GWSA-P and SMSA-P were 1 month in Jordan.
In addition, several studies tried to reveal the reasons for the delayed response of TWSA to P. A study in Ethiopia (Awange et al. 2014) showed that aquifer properties have a strong influence on the lag between TWSA and P, such as 0-month time lag in karst-dominated regions and a lag of up to 6 months in unconsolidated sediment regions. Chen et al. (2002) established that the different recharge characteristics and permeability of sediments influence the response of groundwater level to P. Heterogeneous lithology and the percentage of clay in an aquifer were found by Lorenzo-Lacruz et al. (2017) to influence the response of GWSA and TWSA to P. For example, a longer lag between GWSA and P occurred mostly in areas where high-permeability rocks predominated, with lower correlation coefficients between P and aquifer storage being found in areas of higher clay content. Soni and Syed (2015) showed that the relationship between P and TWSA and the associated lag were influenced by such factors as climatic conditions, vegetation density, runoff pathlength, basin size, and the abundance of surface water bodies. However, despite these various studies, the spatial pattern of the responses of TWSA and its individual components to P, and the influence of climatic factors, vegetation behavior, and basin characteristics on the associated lags, remain unclear over global scale.
Traditionally, TWS is obtained from model simulations and in situ observations. Global Land Data Assimilation System (GLDAS) land surface models can be used to successfully simulate some components of TWS, such as SMS and SWE, but perform relatively poorly in capturing other components, such as GWS (Syed et al. 2008; Xu et al. 2018). To now, in situ measurements of TWS have not been obtained at a global scale. The GRACE satellite, launched in 2002 to monitor Earth’s time-variable gravity field, which reflects Earth’s surface mass redistribution and transport, provides unique data on TWS over large regions (Tapley et al. 2004; Wahr et al. 2004). These data can be used to accurately measure the vertically integrated water resources stored above and beneath Earth’s surface (Jin et al. 2010; Long et al. 2016; Rodell et al. 2004a; Syed et al. 2008; Wahr et al. 1998; Yang et al. 2013). Therefore, using TWS data from GRACE, the goals of this study are to 1) examine the responses of TWSA and its individual components to P for 168 global basins using both basin-average and grid-based approaches, 2) distinguish the differences in time lags between individual components of TWSA and P, and 3) explore the factors that may influence the lag between TWSA and P.
2. Materials and methods
a. Study regions
In this study, we analyze the TWSA–P relationship for 168 global basins. The criteria for basin selection followed those of Scanlon et al. (2016), which are 1) basins with areas > 40 000 km2 and 2) basins with radius > 200 km (Fig. 1). Applying these criteria, data for 168 basins covering 66 × 106 km2 are extracted from the Global Runoff Data Centre (GRDC; http://grdc.bafg.de) (Fig. 1).
1) GRACE data
Global monthly TWSA data retrieved from GRACE are obtained from the latest-release (RL05) mass concentration (mascon) solutions processed by Jet Propulsion Laboratory (JPL) processing centers at NASA and by the Center for Space Research (CSR) at the University of Texas at Austin. Based on the equal-area 3 × 3 spherical cap mascon function, the JPL mascon data uses a priori constraints, including altimetry data and forward models (GLDAS land surface models; Rodell et al. 2004a) for land, and Estimating the Circulation and Climate of the Ocean 2 for oceans (Menemenlis et al. 2008), to estimate the global gravity fields. These procedures can minimize the effect of measurement errors. The coarse 3° × 3° JPL mascon data (0.5° × 0.5° latitude–longitude grid) are then multiplied by the Community Land Model 4.0 (downsampled from 1° × 1° to 0.5° × 0.5°) provided by the Tellus website to reduce leakage errors introduced by the mascon basis function (Watkins et al. 2015). We obtain the JPL mascon data from https://grace.jpl.nasa.gov/data/get-data/jpl_global_mascons/. A detailed description of these JPL mascon data can be found in Wiese et al. (2016). The CSR mascons, calculated based on equal-area geodesic grids of approximately 120 km (1° at the equator), are processed by constraining the original GRACE level-1 data through the Tikhonov regularization method, which effectively suppressed the north–south stripe errors in the GRACE measurements. We obtain the CSR mascon data from http://www.csr.utexas.edu/grace/RL05_mascons.htm. A detailed description of the CSR mascon data is given in Save et al. (2016). Monthly anomalies of TWS are obtained relative to the baseline average for the period January 2004–December 2009. Here, the average of the adjusted months is applied to complement the missing months (Andrew et al. 2017; Long et al. 2015; Mo et al. 2016), and the mean TWSA from the GRACE JPL and GRACE CSR mascons is then calculated. Details of the TWSA dataset used in this study are reported in Table 1.
The global 0.5° monthly P datasets used in the present study are obtained from the Global Precipitation Climatology Centre (GPCC) (Schneider et al. 2014) and the Climatic Research Unit (CRU) Time Series Version 4.00 dataset (Harris et al. 2014), designated as PGPCC and PCRU, respectively. Details of these two datasets are given in Table 1.
3) Individual components of TWS
In this study, the soil moisture (SM), SWE, surface water (SW), and canopy water (CW), are obtained from NASA GLDAS land surface models (LSMs) (Feng et al. 2013; Beaudoing and 2015; Rodell et al. 2004b; https://disc.sci.gsfc.nasa.gov/). Owing to the different versions, structures, and parameters of these models, the LSMs have been previously divided into the Community Land Model (CLM) (Dai et al. 2001), Community Land Model 4.0 (CLM4.0) (Oleson et al. 2010), variable infiltration capacity (VIC) (Liang et al. 1994), Noah-1 and Noah-2.1 (Chen et al. 1996), and Mosaic (Koster and Suarez 1996). We obtain the individual components of TWS, at a spatial resolution of 1°, by averaging them from the six LSMs and then calculating the all individual component anomalies of TWSA by subtracting the average individual components of TWS from 2003 to 2014 (Chen et al. 2018). Details of the components included in the GLDAS LSMs are given in Table 2. Reservoir storage, considered as an important component of water storage (Castle et al. 2014; Famiglietti et al. 2011; Scanlon et al. 2012; Shamsudduha et al. 2017; Soni and Syed 2015; Strassberg et al. 2009), is disregarded in this study because of the lack of data on integrated reservoir storage for global basins.
4) Climatic and basin characteristics data
Data on basin properties, vegetation, and climate characteristics data are collected to investigate how catchments regulate the time lag between TWSA and P. Basin elevation (DEM) is derived from Global Multiresolution Terrain Elevation Data 2010 (GMTED2010) provided by the U.S. Geological Survey and the National Geospatial-Intelligence Agency (https://topotools.cr.usgs.gov/GMTED_viewer/viewer.htm) (Danielson and Gesch 2011). Slope data are obtained from the Food and Agriculture Organization of the United Nations (http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/en/). The P datasets are collected from the CRU and GPCC. Shortwave radiation data R and temperature data T at 0.5° spatial resolution for the period January 2003–December 2014 are obtained from the CRU and National Centers for Environmental Prediction (ftp://nacp.ornl.gov/synthesis/2009/frescati/temp/land_use_change/original/readme.htm). basin vegetation coverage (described by a monthly normalized difference vegetation index (NDVI) dataset with a spatial resolution of 0.0083°) is obtained from the Global Inventory Modeling and Mapping Studies group (https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/; Tucker et al. 2005). Details of all the above data are given in Table 1. These monthly datasets are aggregated to multiyear averages over basins scale.
c. Trend analysis
d. Time lag analysis
Pearson correlation analysis is used to explore the relationships between P and TWSA. For each basin, we obtain the correlation coefficients between P and TWSA for different time lags (0–12 months) by moving the TWSA time series forward 1 month at a time (Frappart et al. 2013; Ndehedehe et al. 2016; Soni and Syed 2015; Yang et al. 2014). We are then able to identify the maximum correlation coefficient and the lag time corresponding to that maximum correlation coefficient. Only positive and statistically significant correlations (p < 0.05) are included in the analysis.
e. Individual TWSA component anomalies
As mentioned in section 2b(3), the values of CWS, SWS, SMS, and SWE are averaged from six LSMs. We calculated the anomalies as
where SMSuυ is the SMS in month u of year υ and is the average SMS for the period January 2003 to December 2014 (Chen et al. 2017); the terms for the other components are defined similarly.
f. Quantifying groundwater storage anomalies
Theoretically, TWS is the sum of its individual components (Scanlon et al. 2016):
In addition to GWS, GLDAS-based TWS incorporates SWS, SMS, CWS, and SWE. Therefore, GWS can be obtained by differencing GLDAS-based TWS from GRACE TWS (Jin and Feng 2013; Rodell et al. 2009). GWSA can thus be calculated as
Jin and Feng (2013) found that GWS calculated from GRACE–GLDAS is robust for estimating global GWS. In addition, a good agreement has been found between GWS calculated from GRACE–GLDAS and in situ observations (J. Chen et al. 2014; T. Chen et al. 2014; X. Chen et al. 2014; Jin and Feng 2013; Long et al. 2016; Xiao et al. 2015).
g. TWSA component contribution ratios
We also determine the contributions of individual components (GWSA, SMSA, SWSA, and SWEA) to TWSA, which is crucial for understanding the effect of these components on land water changes. We apply the component contribution ratio (CCR) proposed by Kim et al. (2009) to calculate the contribution of individual components to TWSA for the studied global basins as follows:
where MAD is the mean absolute deviation of the individual component [; S refers to the individual storage components, for example, GWS, SMS, and SWS; St is the value of individual components in month t; and N is the number of months] and TV is the sum of the MAD of all components ().
h. Linking lags to climatic and basin characteristics
Linear regression analysis is used to evaluate the relationship between the lags of the response of TWSA to P and the climatic and basin characteristics. The climatic characteristics include: (i) basin mean temperature, (ii) basin mean precipitation, (iii) basin mean short radiation, and (iv) basin mean NDVI. The related basin characteristics include (i) mean latitude, (ii) mean longitude, (iii) mean elevation, (iv) mean slope, and (v) mean area. Table S1 in the online supplemental material lists the statistical measures of the climatic and basin characteristics for the 168 studied basins.
Multiple regression analysis is used to construct a statistical model of the relationship between climatic, vegetation, and basin characteristics and the lag between GRACE-based TWSA and P. Finally, for evaluating the applicability of the regression models, the normality, homoscedasticity, and independence of the residuals are checked, as detailed below.
i. Residuals analysis
1) Homoscedasticity test
When the variance of errors is the same for different values of the independent variables, homoscedasticity is indicated. A plot of the standardized residuals (the error) versus the regression standardized predicted value is used to check for homoscedasticity. The presence of homoscedasticity in the plot is indicated by the residuals being randomly scattered around 0 and the data scatter showing no clear pattern (Osbourne and Waters 2002).
2) Independence test
The test for independence of residuals involves checking that the residuals are uncorrelated (Osbourne and Waters 2002). The Durbin–Watson statistic is used to test for the presence of autocorrelation, with the value of the statistic ranging from 0 to 4. Values close to 2 indicate that the residuals are independent. Values near 0 indicate the presence of significant positive autocorrelation and those near 4 indicate the presence of significant negative autocorrelation (Altman and Bland 1995; Chan 2004; Durbin and Watson 1971).
3) Normality test
The regression model used is also based on the assumption that residuals follow a normal distribution (Altman and Bland 1995; Driscoll et al. 2000; Field 2013; Pallant 2001). Visual inspection of the distributions, namely, the frequency distribution (histogram) and the probability–probability (p–p) plot (Field 2013), are used for assessing normality. The histogram (of standardized residuals) provides a visual judgment of whether the distribution of the residuals is bell shaped (Thode 2002). The p–p diagram shows the cumulative probability of the expected values plotted against the cumulative probability of the observed values. Data plotting on or near a straight diagonal line indicate normally distributed residuals (Field 2013).
a. Spatial and temporal trends of precipitation and TWSA over global basins
Trend analysis is a widely used tool among hydrologists (Scanlon et al. 2018) and is used here to examine the covariation of TWSA and P. Figure 2 (top: PCRU) and Fig. S1 (PGPCC) show the spatial patterns of the trend of P for 168 global basins for the period January 2004 to December 2013. Results refer to that the trends of P based on the CRU and the GPCC are consistent with each other. However, there is no significant variation trend existing in P over global 168 basins. Among which, a nonsignificant decrease can be observed in South America (e.g., the Amazon and Orinoco basins), in southern North America (e.g., the Mississippi basin), in the Indus basin, and in southern Africa (e.g., the Orange and Zambezi basins). However, the increased trend for P is mainly in eastern Asia (e.g., the Yangtze and Amur basins), in northern Asia (e.g., the Lena, Yenisei, and Ob River basins), in eastern Australia (e.g., the Murray basin), and in northern Africa (e.g., the Nile and Congo basins).
Figure 2 (bottom) shows the spatial pattern of GRACE-based TWSA over for global basins. Most low-latitude basins show an increase (significant) in TWSA (26%), including the Nile, Congo, Murray, and Zambezi basins. Increased TWSA is observed in several basins in the midlatitudes of the Northern Hemisphere, such as the Mississippi and Yangtze basins. Decreases in TWSA are observed for 30.4% of global basins, which are seen mostly in basins at high and midlatitudes of the Northern Hemisphere, such as the Lena, Yenisei, and Ob basins, and those in northern India. These identified patterns of TWSA are consistent with those noted by some previous global or regional-scale investigations, such as those of Ahmed et al. (2014) and Scanlon et al. (2016).
Coherence between P and TWSA is found for most of the studied basins, but TWSA shows more instances of significant increase and decrease trends, implying potential influences of other factors on TWSA besides P. Interestingly, for some basins, such as the Ob, Lena, and Yenisei basins, increasing P corresponds to decreasing TWSA. However, in the Amazon basin, for example, P shows a decreasing trend and TWSA an increasing trend.
b. Response of TWSA to precipitation
Figure 3 shows the variation in correlation coefficient between TWSA and P with different lags (0–12 month) for 31 large basins (shown in Fig. 1). For some basins, the correlation coefficient transforms from positive to negative with the increasing of the lag months, particularly in low- and midlatitude basins, such as the Amazon basin (Fig. 3, basin label 1). Some basins change from negative to positive, which mainly occur in high latitudes, such as in the Mississippi basin (label 3). Similar patterns are shown for TWSA and PGPCC (Fig. S2).
The maximum correlation coefficients between GRACE-based TWSA and P, and the lags (measured in months) corresponding to those maximum correlation coefficients, are determined using both a grid-based approach (Figs. 4a,b, Figs. S3a,b) and a basin-average approach (Figs. 5a,b, Figs. S4a,b). The spatial patterns of the maximum correlation coefficient and the corresponding lag month are similar for both approaches. The correlation coefficients are significantly positive over several very large land areas (Fig. 4a and Fig. S3a), including northern South America, Africa, southern India, southern China, and most of Russia. At a basin scale, the larger correlation coefficients are observed for the Amazon basin (0.92 for PCRU and GRACE-based TWSA; 0.94 for PGPCC and GRACE-based TWSA), Nile basin (0.86 and 0.84), Niger basin (0.89 and 0.90), Yangtze basin (both 0.85), and Mekong basin (0.90 and 0.93) (Fig. 5a and Fig. S4a).
The lags between GRACE-based TWSA and P corresponding to the maximum correlation coefficients are shorter for basins in low latitudes than for those in mid and high latitudes. Generally, the GRACE-based TWSA at low latitudes responds to P with a 2–3-month lag but in the mid and high latitudes with a 7–9-month lag. For most basins in South America, Africa, Australia, China, India, and Australia, the lag is 1–2 months (Fig. 4b and Fig. S3b). At the basin scale, the lower lags are found in the Amazon (2 months), Congo (2 months), Nile (2 months), Yangtze (1 month), and Ganges (2 months) basins (Fig. 5b and Fig. S4b). The larger lags are observed mostly in the mid and high-latitude parts of the Northern Hemisphere and in southern South America, southern Africa, and southeastern Australia (Fig. 4b and Fig. S3b). At a basin scale, these include the Mississippi (9 months), Ob (8 months), Yenisei (8 months), Lena (8 months), and Murray (7 months) basins (Fig. 5b and Fig. S4b).
In addition, the maximum correlation and the corresponding lag of GLDAS-based TWSA and P are evaluated (Figs. 4c,d and 5c,d, Figs. S3c,d and S4c,d). We find that the maximum correlation coefficients are positive (Figs. 4c and 5c, Figs. S3c and S4c) and that the spatial pattern of the lags between GLDAS-based TWSA and P is similar to that of the lags between GRACE-based TWSA and P (Figs. 4d and 5d and Figs.S3d and S4d). However, the lag between GLDAS-based TWSA and P (about 7–8 months at high latitudes) is shorter than that between GRACE-based TWSA and P (about 9–10 months at high latitudes).
Furthermore, we calculate the difference of the lag months between GRACE-based TWSA to P and GLDAS-based TWSA to P by the grid-based approach (Fig. 6a and Fig. S5a) and basin-based approach (Fig. 6b and Fig. S5b), respectively. We find that the spatial patterns for both methods are similar. The differences of the lag months are about 1 month for most of the global areas, which are mainly in the Africa (Congo basin, Nile basin, Niger basin, etc.), in the northern South America (Amazon basin and Parana basin, etc.), in India (Ganges basin, etc.), in the south of the China (Mekong basin, etc.), and in the most of the north of Asia (Ob basin, Yenisei basin, and Lena basin, etc.). In addition, we find that the differences are larger in most of southern North America (it is mainly shown for the grid-based approach), in Australia (Murray basin), and in the east of northern Asia (it is also mainly shown for the grid-based approach). These possible reasons for the differences of the lag month between GRACE-based TWSA and GLDAS-based TWSA may be caused by 1) the lack of GWSA in GLDAS-based TWSA, 2) uncertainties of the GLDAS-based TWSA (Nie et al. 2016; Scanlon et al. 2018), and 3) errors of the GRACE TWSA (Scanlon et al. 2016, 2018).
c. Responses of the individual components of TWSA to precipitation
1) Separating the components of TWSA
Figure 7 shows the CCRs of GWSA, SMSA, SWEA, and SWSA for 31 large global basins. The results are similar to those of Felfelani et al. (2017), who determine the CCRs for subsurface water storage, river and reservoir water, and snow water storage with respect to TWS based on the HiGW-MAT model. In the present study, SWEA accounts for the highest proportion of TWSA in the northern high-latitude basins, including the Kolyma (48.6%), Mackenzie (50%), Yenisei (48.9%), Lena (43%), and Ob (42.2%) basins. In the midlatitude and subtropical basins, the CCRs of SWSA, GWSA, and SMSA are high, and the CCR of SWEA is low. In tropical basins, such as the Congo and Ganges basins, GWSA and SMSA dominate TWSA. It should be noted that because of the lack of global reservoir water, we consider only runoff and snowmelt as SWS, which may result in an underestimation of the CCR of SWSA, especially for subtropical basins and those basins with managed reservoirs, such as the Yangtze, Brahmaputra, and Ganges basins.
2) Responses of the individual components of TWSA to precipitation
Responses of the individual components of GRACE-based TWSA to P are determined by correlation analysis using both a grid-based approach (Fig. 8 and Fig. S6) and a basin-average approach (Fig. 9 and Fig. S7), respectively. Results show that the maximum correlation coefficients of GWSA (Figs. 8a and 9a, and Figs. S6a and S7a), SWSA (Figs. 8c and 9c, and Figs. S6c and S7c), SMSA (Figs. 8e and 9e, and Figs. S6e and S7e), and CWSA (Figs. 8g and 9g, and Figs. S6g and S7g) with P are significantly positive for most regions and basins. Moreover, higher correlation coefficients are found in the low and midlatitudes than in the high latitudes regardless of which method is used. For SWEA, the larger correlation coefficients are observed mainly in the mid and high latitudes (Figs. 8i and 9i, and Figs. S6i and S7i).
Generally, the lag months of SWSA and SMSA with P are shorter in the low and midlatitudes, usually ranging from 0 to 3 months, compared with the lags of 8–12 months found in high-latitude basins (Figs. 8d,f and 9d,f, and Figs. S6d,f and S7d,f). For GWSA, the longer lags are observed mainly in western and northern North America, most of central and northern Eurasia, southern South America, southern Africa, and southeastern Australia (Figs. 8b and 9b, and Figs. S6b and S7b). The lags between SWEA and P are longer in mid- and high-latitude basins than in low-latitude basins (Figs. 8j and 9j, and Figs. S6j and S7j). In particular, for the Middle East and Tibetan Plateau, the lags are about 10 months. As expected, the lags between CWSA and P are short or nonexistent for all basins, revealing that CWSA can respond immediately to changes in P (Figs. 8h and 9h, and Figs. S6h and S7h).
Furthermore, the lags of GWSA, SMSA, and CWSA with P are compared (Figs. S8 and S9) within global areas and basins. We found that the lag can be arranged as GWSA > SMSA ≥ SWSA in most basins and regions (Note that the yellow areas and basins with blue boundaries indicate that the lags of GWSA, SMSA, and SWSA with P are ordered as GWSA ≥ SMSA ≥ SWSA).
d. Correlation of climatic and basin characteristics with the TWSA–P lag
The data described in section 2b(4) are analyzed to produce measures of climatic and basin characteristics. The climatic characteristics are mean T, mean P, and mean R. The vegetation data is mean NDVI. The basin characteristics are mean latitude, mean DEM, and mean slope and basin area. The relationships between these characteristics and the lag months are evaluated by linear regression analysis for the 168 global studied basins (Fig. 10 and Fig. S9). According to the linear fit line and the significance level shown in the Fig. 10 and Fig. S9, we found that the basins mean T (p < 0.05) (Fig. 10a and Fig. S9a), mean R (p < 0.05) (Fig. 10b and Fig. S9b), mean P (p < 0.05) (Fig. 10c and Fig. S9c), mean DEM (p < 0.05) (Fig. 10d and Fig. S9d) show significantly negative correlations with the lag months of GRACE-based TWSA and P. The significant positive correlation is shown for the relationship of mean latitude (p < 0.05 in the north latitude) (Fig. 10e and Fig. S9e) and mean slope (p < 0.05) (Fig. 10f and Fig. S9f) with the lag months. However, mean NDVI (Fig. 10g and Fig. S8g) and mean basin area (Fig. 10h and Fig. S9h) are found to play negligible roles in determining the lag months between TWSA and P.
e. Regression model of the TWSA–P lag
Given the results of section 3d, the data for mean T, mean R, mean P, mean DEM, mean latitude, and mean slope are used to construct a statistical model to explain the lag months of the 168 studied basins. The following multiple linear regression models are established for the lag months of GRACE-based TWSA to PCRU (Lag_CRU) and the lag months of GRACE-based TWSA to PGPCC (Lag_GPCC), respectively:
where X1, X2, X3, X4, X5, and X6 are mean T, mean R, mean P, mean DEM, mean latitude, and mean slope, respectively. Equations (5) and (6) explain 66.0% and 63.9% of the observed variance in Lag_CRU and Lag_GPCC, respectively.
We also use tests for independence, normality, and homoscedasticity of the residuals as described in section 2i to check the applicability of the regression models. The Durbin–Watson test for independence yields a value of 1.7 for PCRU and 1.9 for PGPCC, which indicates that the residuals are independent. The histograms and p–p plots presented in Fig. S10 clearly show that the residuals satisfy the normality assumption. The scatterplots in Fig. S11 show no clear pattern, and the line of best fit is nearly horizontal (slope = 6.74 × 10−16 for PCRU and slope = 4.88 × 10−16 for PGPCC), which together indicate that the residuals are homoscedastic.
4. Discussion and conclusions
This study examines the response of GRACE-based TWSA and the individual components of GRACE-based TWSA to P using a basin-average approach for global 168 basins for the period January 2004–December 2013. We quantify the maximum correlation and lag between TWSA and P, and that between the individual components of TWSA and P, for the studied basins, with great implications for understanding the role of precipitation in the global hydrological cycle.
We investigate the spatial pattern of the trend of global P based on GPCC and CRU datasets and GRACE-based TWSA. Consistent with Scanlon et al. (2016), Ramillien et al. (2014), and Humphrey et al. (2016), increases in GRACE-based TWSA are identified for various basins including the Nile, Niger, and Yangtze basins. Humphrey et al. (2016) and Ahmed et al. (2014) explained the increase in GRACE-based TWSA in the Niger basin is related to the rising precipitation. However, in some basins, such as the Ob, Lena, and Yenisei basins, a nonsignificant increase in P is associated with a decrease in GRACE-based TWSA. Herein, according to the analysis of the period of August 2002–March 2008 by Muskett and Romanovsky (2009), it may attribute that the GRACE-based TWSA trends are related to groundwater, runoff, and snowmelt. For the Amazon basin as a whole, a nonsignificant decrease in P and an increase in GRACE-based TWSA are observed, consistent with the findings of Humphrey et al. (2016). The 2009 Amazon flood event may explain this phenomenon (Chen et al. 2010). The groundwater depletion in the northwest India likely causes the observed decrease in GRACE-based TWSA for that region (J. Chen et al. 2014; Humphrey et al. 2016; Rodell et al. 2009).
The maximum correlation coefficients between GRACE-based TWSA and P, as well as the corresponding lag months to this maximum value, are calculated using both a basin-average approach and a grid-based approach. Results show that the spatial patterns produced by these two approaches are similar, including the positive relationship between TWSA and P in most basins and regions and the common lag between TWSA and P. The lags are shorter in the low- and midlatitude basins, confirming the findings of previous studies in those regions, including the 1–2-month lags in India (Soni and Syed 2015), the ~2-month lag in the Amazon basin (Famiglietti et al. 2011), and the 2-month lags in Africa (Ndehedehe et al. 2016). Humphrey et al. (2016) argued that the particular water storage processes could cause the short lags in the tropical and subtropical region. Such as the transport processes of water within the extensive floodplains caused the lag of P to TWSA in the Amazon basin (Frappart et al. 2013). However, longer lags are found in the high latitudes, which are likely due to the effect of seasonal snowmelt (Humphrey et al. 2016; Trautmann et al. 2018).
We also analyzed the relationships between individual components of GRACE-based TWSA and P as well as the corresponding lag months. The relationship between SWSA and P is significantly positive. Moreover, the lags between SWSA and P are shorter in the low- and midlatitude basins (0–3 months), suggesting a more direct control of P on SWSA (Huang et al. 2019), whereas the lag months of SWSA to P are longer in the high latitude basins (8–12 months). Previous studies have argued that lags in high-latitude basins may not alone be explained by precipitation. Liesch and Ohmer (2016) referred to the idea that P and evapotranspiration (ET) were the main controls on SWS change and that ET relied mainly on T. Therefore, we analyze the relationships of T and P with SWSA (Figs. S12–S14). Results show that in low-latitude and coastal subarctic areas, P controls SWSA, with a shorter lag (0–1 month). In mid- and high latitudes, T dominates SWSA, and the lags are 3–6 months (Fig. S11). Moreover, we found that the response of SWSA to T mostly occurred in the spring and autumn seasons with different lag months. The reason may be that the lower SWSA occurs in the autumn, lagging the peak T by 4–6 months (Fig. S15) (Humphrey et al. 2016), and during spring, SWSA responds rapidly to T (yellow area in Fig. S15). In coastal subarctic areas (such as the Columbia and Colorado basins), owing to the temperature difference between land and ocean, the maximum P occurs during the autumn, which results in the response of SWSA to P (Humphrey et al. 2016).
Similar to the relationship between SWSA and P, the relationship between SMSA and P is significantly positive, and the lags are shorter (0–3 months) in the low- and midlatitude basins. However, the longer lag months are observed in high latitudes. Somorowska (2017) explained variation in SMS variability in Poland not only in terms of precipitation but also in ET, which depending on meteorological conditions, vegetation, and soil hydraulic properties, can also directly influence the surface soil layers. In the Songhua basin, northeastern China, the lag of SMSA mainly influence by precipitation, snow melting, and seasonal permafrost thawing processes (Chen et al. (2019).
For GWSA, the longer lags are found mainly in western and northern North America, most of central and northern Eurasia, southern South America, southern Africa, and southeastern Australia. The lags between GWSA and P are controlled by several factors. In the Songhua basin, Chen et al. (2019) found increasing in GWSA in 10–12 months, due to the mainly exploitation of wildly distribution of snow, ice, and seasonal permafrost. Under snowpack accumulation and melting processes, groundwater was recharged in 1–6 months. It can explained as the end of snowpack melting, the increasing of temperature, evapotranspiration and irrigation activity. From 6 to 10 months, the GWSA was mainly influenced by precipitation. However, the recharged of P to GWSA in these seasons is far less than the consuming of groundwater by human activities such as irrigation, groundwater exploitation and the evapotranspiration caused by climate factors, such as temperature and speed, etc. In addition, the soil characteristics can also affect the lag between GWSA and P. For example, Huang et al. (2019) found that the lag between GWSA and P was shorter in the high karstic region than in the low karstic region owing to the soil characters and Topographic condition.
SWEA is influenced by T (Butt and Bilal 2011; Chen et al. 2019). Lags between SWEA and P are longer (~10 months) and are mainly in the mid- and high latitudes particularly in the Middle East and Tibetan Plateau. As expected, the lags between CWSA and P are short or nonexistent for all studied basins, revealing that CWSA can respond immediately to P.
Furthermore, the lags of GWSA, SMSA, and SWSA with P can be arranged as GWSA > SMSA ≥ SWSA for most of the large basins. This pattern may be related to the distribution of P in the atmosphere and on the land surface. Oki and Kanae (2006) and Chen et al. (2018) stated that water evaporates from the ocean, converting from liquid to gas, and then changes to a liquid (or solid) and falls as precipitation to the ground. Some of this water flows as surface runoff, some soaks into the soil, and some becomes groundwater. These broad components of the hydrological cycle provide a first-order explanation of the differences in the lags of GWSA and SMSA with P.
Using multiple linear regression models, we also examined the relationships between climatic factors and basin characteristics and the lag between GRACE-based TWSA and P. Result shows that the mean T, mean P, mean R, mean latitude, mean DEM, and mean slope are related to the lag months. Temperature mainly influences the change of SWSA and SWSA (Liesch and Ohmer 2016). The influence of mean latitude and mean slope and DEM on the TWSA–P lag shows the control of basin attributes, but these attributes also in part reflect climatic characteristics, such as mean temperature and radiation (Parry et al. 2016). The result is consistent with Yang et al. (2017), who indicated that the basin’s mean DEM and latitude significantly control the lag of hydrological drought recovery to that of meteorological drought in Australia. In the present study, the basin mean T, mean P, mean latitude, mean R, mean DEM, and mean slope together explain 66.0% and 63.9% of the observed variance in Lag_PCRU and Lag_PGPCC, respectively.
The remaining unexplained degree of explanation may be related to variation in soil characteristics such as soil porosity and permeability (Huang et al. 2019) and irrigation activity (Chen et al. 2019), etc.
TWS plays an important role in climatological and hydrological processes. This study conducted the first subtle evaluation of the response of TWSA and the individual components of TWSA to P over global basins. Beyond the above findings, this study should help hydrological researchers to understand the interactions and mechanisms operating between climatic conditions and variation in TWS, as well as the nature and speed of the processes by which precipitation becomes surface water, soil water, and groundwater. Because of the complex processes involved in water circulation within soil, especially in high latitudes (given the frozen surfaces), soil water circulation is the most difficult process to be described by LSMs (Chen et al. 2019). Bao et al. (2017) found that the correlation between GLDAS-based soil moisture data and observed data was higher during the April–October nonfrozen period (compared with the November–March frozen period) over the Tibetan Plateau, which indicates that GLDAS soil moisture estimates have certain representativeness. However, uncertainties exist, including the change in soil moisture with temperature during freezing–thawing in the frozen period and the bias of values caused by GLDAS precipitation-forcing data in the nonfrozen period. In the Songhua basin, Chen et al. (2019) showed that SMSA estimates from different GLDAS LSMs have apparent discrepancies. Therefore, the relationship between SMSA and P can be influenced by uncertainties caused by the SMS dataset in different LSMs. In addition, uncertainties in SMS may then propagate to GWS. Consistent with our study, previous studies have often taken averages of SMS based on different GLDAS LSMs to reduce the uncertainties (Feng et al. 2013; Gong et al. 2018; Soni and Syed 2015).
This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant XDA20060402), the National Key Scientific Research and Development Program of China (Grant 2017YFA0603601), and the National Natural Science Foundation of China (Grant 41671083). The authors thank Min-Hui Lo for providing the GLDAS CLM4.0 dataset. This paper was supported by the GRACE RL05 Mascon Solutions data (https://grace.jpl.nasa.gov/data/get-data/) that was obtained from the JPL and the CSR. The global P dataset was produced by the GPCC (https://www.esrl.noaa.gov/psd/data/gridded/data.gpcc.html) and the Climatic Research Unit (CRU Time Series v4.00) datasets (http://www.cru.uea.ac.uk/data/). Six global land models were applied here, and the basins were obtained from the GRDC.
List of Acronyms
Component contribution ratio
Community Land Model
Climatic Research Unit
Center for Space Research
Canopy water storage
Canopy water storage anomalies
Climatic Research Unit, National Centers for Environment Prediction
Global Land Data Assimilation System
Global Precipitation Climatology Centre
Gravity Recovery and Climate Experiment
Groundwater storage anomalies
Jet Propulsion Laboratory
Normalized difference vegetation index
Standardized groundwater index
Soil moisture storage
Soil moisture storage anomalies
Standardized precipitation index
Snow water equivalent
Snow water equivalent anomalies
Surface water storage
Surface water storage anomalies
Total runoff integrating pathway
Terrestrial water storage
Terrestrial water storage anomalies
Terrestrial water storage change
Variable infiltration capacity
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JHM-D-18-0253.s1.