{'fraction_non_alphanumeric': 0.03229878553084055, 'fraction_numerical': 0.013842336656074523, 'mean_word_length': 5.033088235294118, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Recent global decline in endorheic basin water storages

Endorheic (hydrologically landlocked) basins spatially concur with arid/semi-arid climates. Given limited precipitation but high potential evaporation, their water storage is vulnerable to subtle flux perturbations, which are exacerbated by global warming and human activities. Increasing regional evidence suggests a probably recent net decline in endorheic water storage, but this remains unquantified at a global scale. By integrating satellite observations and hydrological modelling, we reveal that during 2002–2016 the global endorheic system experienced a widespread water loss of about 106.3 Gt yr−1, attributed to comparable losses in surface water, soil moisture and groundwater. This decadal decline, disparate from water storage fluctuations in exorheic basins, appears less sensitive to El Niño–Southern Oscillation-driven climate variability, which implies a possible response to longer-term climate conditions and human water management. In the mass-conserved hydrosphere, such an endorheic water loss not only exacerbates local water stress, but also imposes excess water on exorheic basins, leading to a potential sea level rise that matches the contribution of nearly half of the land glacier retreat (excluding Greenland and Antarctica). Given these dual ramifications, we suggest the necessity for long-term monitoring of water storage variation in the global endorheic system and the inclusion of its net contribution to future sea level budgeting. Hydrologically landlocked basins worldwide have experienced widespread decline in water storage over the past decade.


Supplementary Text
Supplementary Figures S1 to S21 Supplementary Tables S1 to S6

Validation of modeled snow water equivalent
Modeled monthly time series of snow water equivalent (SWE) is validated against the Snow Data Assimilation System (SNODAS) 1 data product (version 1, accessed from https://nsidc.org/data/g02158). SNODAS is a modeling and data assimilation system developed by the National Operational Hydrologic Remote Sensing Center (NOHRSC) of the NOAA National Weather Service. It provides some of the highest-quality estimates of snowpack properties, e.g., snow depth, SWE, temperature, and melt runoff, across the conterminous United States and southern Canada. Its spatial coverage includes ~60% of the total area of the endorheic North America (Fig. S19).
We processed the SNODAS SWE daily product (at 1 km spatial resolution) to monthly time series of average SWE (in mm) within the covered endorheic North America. The processed SNODAS SWE time series is compared with the region-average monthly SWE simulated by the seven hydrological models as described in Methods (also see Table S1).
Considering that the SNODAS product is not available until late 2003, our validation is performed on the complete annual cycles (12 years) from April, 2004 to the end of our study period (March, 2016). Consistent with our methods for TWS partitioning, the modeled SWE trend is fitted from the output ensemble of CLM, Mosaic, Noah, VIC, and WGHM models, and the 95% confidence interval (CI) is propagated from (i) SWE variations among these five models, (ii) the discrepancy between the ensembles of the five models and all seven hydrological models (including CLSM and PCR-GLOBWB) til the end of 2014, and (iii) the residuals of linear trend fitting. The CI for the SNODAS trend is inferred from the residuals of its trend fitting alone.
As shown in Fig. S21a, modeled full SWE time series generally follow that of the SNODAS product, although the models tend to underestimate the seasonal maxima by an average factor of ~2. The agreement between models and SNODAS appears to be overall improved after the removal of climatology (Fig. 21b)

Validation of modeled soil moisture
We validate modeled monthly soil moisture (SM) time series against in situ SM measurements through the Soil Climate Analysis Network (SCAN) (accessed at www.wcc.nrcs.usda.gov/scan) administrated by the National Water and Climate Center of the USDA Natural Resources Conservation Service. Typically, a SCAN station monitors volumetric moisture content at various soil depths, together with some other energy and meteorological parameters (e.g., air temperature, relative humidity, and liquid precipitation). The standard soil depths include 2, 4, 8, 20, and 40 inches or about 5, 10, 20, 50, and 100 cm (a total of 5 layers).
All of these stations happen to be located in the Great Basin. IDs of the 22 SCAN stations are labeled (in red) in Fig. S19, and their names and geographic coordinates are provided in Table   S6.
As described in Methods, our SM trends are calculated from the ensemble of the four GLDAS LSMs (CLM, Mosaic, Noah, and VIC). For each LSM, SM values are integrated from the simulated water contents throughout the entire soil compartments (typically covering surface, vadose, and other subsurface layers). The soil compartment properties, however, are model-specific (see Table S1). The maximum soil depth ranges from 1.9 m in VIC to 3.4 m in CLM, which all exceed the maximal depth measured by the SCAN stations (~1 m). The total soil layer number ranges from 3 in VIC and Mosaic to 10 in CLM, and the depths of individual layers are not always consistent with those of the SCAN measurements. These factors cause some challenges in the SM validation. As a solution, we compare the modeled and in situ volumetric SM contents (in %) instead of absolute water depths. For each station, we calculate (i) the correlation between deseasonalized in situ and model-ensemble monthly SM anomalies integrated through their total soil layers, and (ii) both the correlation and interannual trends in the SM anomalies for the soil layer of 0-10 cm. We focus on the top 10 cm soil because this layer is the most commonly available: it is shared by the SCAN stations and three of the four LSMs (CLM, Noah, and VIC). Table S6, out of the 19 SCAN stations where measurements are available for at least 5 consecutive years in all soil layers, 17 stations have in situ measurements that are significantly correlated (p < 0.05) with our model-ensemble anomalies. Out of the 22 stations where measurements of over 5 consecutive years are available for the top 10 cm soil, 21 stations have in situ measurements significantly correlated with the modeled anomalies, and the correlations remain significant for 20 of them when the significance level (α) is lowered to 0.01.

As shown in
Although the SCAN measurements provide an important in situ validation source, their measurements are made at each station site while our modeled SM contents represent the condition of the model pixel (0.25 degree or coarser) that includes the station site. Such spatial inconsistency causes inevitable discrepancies between modeled and in situ SM contents. Even though, the interannual trends in modeled and in situ time series for the top 10 cm soil are comparable for most of the 22 SCAN stations (Table S6 and  It is worth noting that GLDAS LSMs do not include human processes, so our calculated SM trends do not account for the impacts of groundwater-fed or surface-water-fed irrigations. Irrigated areas cover ~10% or much less of the arid/semiarid regions (according to the Global Map of Irrigation Areas (GMIA, www.fao.org/nr/water/aquastat/irrigationmap/index60.stm) 2,3 .
The purpose of irrigation is to supplement local SM and satisfy the crop's biophysical demands (e.g., for transpiration) during the growing season. Excess irrigation (water exceeding the crop's immediate consumption) can temporarily increase the local SM. If the soil is saturated, the excess water will be drained, evaporated, or percolate to the groundwater table. The remnant moisture supplement may eventually be consumed by the crop, rather than changing long-term SM. For these reasons, we assume that during the GRACE period, direct irrigation impacts on SM were regional and transient (e.g., limited to monthly to seasonal timescales), and did not substantially alter the inter-annual SM trends at the studied zonal to global scales. We here acknowledge that the accumulated impact of irrigation could be more substantial over the past century or several decades, due to rapid expansion of irrigated areas since ~1950s 4 . Existing studies also showed that altered water and energy budgets by irrigation could increase precipitation via surface-atmosphere coupling, which in turn affected SM [5][6][7] . However, such impacts on precipitation and energy budgets may have been integrated to the GLDAS climate forcing inputs.

Linkage of water storage changes to climate variability
Periodic variability of the climate system can strongly influence regional meteorological patterns and the associated TWS [8][9][10] . Recent studies showed that detrended interannual TWS variations during the GRACE era were correlated with the El Niño-Southern Oscillation (ENSO) over much of the globe, particularly tropical and subtropical regions (e.g., the Amazon and Southeast Asia). In these analyses, the calculated trend was usually removed from the deseasonalized TWS time series to isolate interannual signals, assuming that the trend was mainly induced by "secular" climate change and/or direct human impacts rather than interannual climate variability [11][12][13] . Although this is valid for many cases, cautions are needed when associating detrended TWS with climate variability during a relatively short period. For example, the first decade after 2000 (till 2012 or so) favored La Niña episodes more than El Niño, leading to excess precipitation on land and a hiatus in sea level rise 14 . This ENSO asymmetry, an effect of interannual climate variability, contributed to a net decadal increasing trend in the TWS (excluding ice sheets and glaciers) on the global continental surface 15,16 .
Given this caveat, we compare the intensity of climate variability with deseasonalized TWS anomalies (hereafter TWSA) without removal of the linear trend. As such, a strong correlation between TWSA and climate variability suggests the importance of climate variability on TWS variations despite the possible existence of secular factors (e.g., global warming and groundwater depletion). Similarly, a weak correlation, in combination with a strong TWS trend, implies a possible dominance of the effect of more secular factors over that of short-term climate variability. Importantly, we here do not aim at a quantitative separation of the effect of climate variability from that of secular factors, but to examine whether climate variability may considerably contaminate or determine our observed endorheic TWS decline.
We here examine five prominent modes of climate variability related to sea surface temperature or pressure anomalies in the Pacific, Atlantic, or Indian Oceans ( Fig. S12; Table   S4) Table S4).
For each climate index, we calculate its correlation with TWSA in both global exorheic/endorheic basins and the 10 regional endorheic zones. On a basin scale, the flux of TWS (i.e., ΔTWS between each time interval) reflects the residual between precipitation (the source), and evapotranspiration and runoff (the sinks). Numerous studies suggested that climate variability strongly influences precipitation anomalies in various regions 9,21-27 . As runoff and evapotranspiration depend on the availability of precipitation, their anomalies are also affected by climate variability 11 . However, the sink rates are not solely determined by contemporaneous weather conditions (e.g., precipitation inputs). They are often related to the status of historical land water storage (e.g., baseflow and lake/wetland water detention), which can result in evident time lags between climate index and TWS flux. When the climate index misrepresents the contribution of the sinks, we may observe the index to (lag-)correlate with TWSA more than ΔTWSA (changes between consecutive TWSA) 11,13 , given that TWS integrates ΔTWS and the previous TWS related to the sinks. For this reason, we explore the correlations between each climate index with both TWSA and ΔTWSA. To target interannual scales, time series of climate index and TWSA are smoothed by an average low-pass filter with a 13-month window. The Pearson correlation coefficients (r) are then computed between filtered climate index and TWSA/ΔTWSA time series for different time lags at a monthly step. We here constrain the lag within one year, following a range of ~3-12 months as applied/implied in some existing studies 11,13,21,28 . The pattern of TWSA with PDO is similar to that with ENSO ( Fig. S12c), except that the correlation between PDO and exorheic TWSA is more moderate (r = -0.54). This conforms to the fact that PDO is inherently related to ENSO at interannual timescales, but includes additional variability at decadal timescales 21  As described previously, phases in IOD and AMO have been mainly positive during our study period, although recent studies identified a strong negative IOD in late 2016 19 and a likely emerging cold phase in AMO 29 . Apparently, the limited GRACE era precludes a diagnostic comparison between their full-cycle variations and TWS observations. Based on our results In Table S4, we also report the maximum (lag-)correlations between each climate index and ΔTWSA in global exorheic/endorheic basins. The correlations are overall weaker than those with TWSA, but presents a consistent pattern that the global endorheic TWS is generally less correlated to interannual climate variability.
Although endorheic TWS is less sensitive to interannual variability on a global scale, it is important to note that stronger correlations exist in some regional endorheic zones (Table S4).
For example, in Dry Andes and Patagonia, TWSA is strongly correlated with ENSO and PDO (r ≥ 0.60) and moderately correlated with NAO (r = 0.55). TWSA in the Inner Tibet Plateau (ITP) is negatively correlated to PDO and ENSO, and to a lesser extent, positively correlated to IOD. In endorheic North America, TWSA shows moderate to strong correlations to all climate indices except IOD. Moderate linkages to multiple internanual indices (ENSO, NAO, and PDO) are also found in Australia and Great Rift Valley/Southern Africa (GRVSA). These results imply that climate variability, at least on interannual to subdecadal timescales, may exert important influences on the above-mentioned endorheic regions, which is generally consistent with existing literature 11,15,22,27,[30][31][32][33] . However, these regions are limited in area (accounting for 31% of the global endorheic area), and include some of the endorheic zones where TWS declines are less dramatic (e.g., Australia and North America), or TWS experienced increasing trends (e.g., ITP and GRVSA). As a result, they end up with a net but uncertain TWS increase of 5.5 (±6.1) Gt yr -1 (Table 1). In comparison, the remaining regions in Sahara/Arabia and Central Eurasia (except ITP), the dominant contributors to the global endorheic TWS decline, are overall less correlated to interannual climate indices, which likely suggests substantial impacts from more secular factors such as human water managements. Again, this is largely consistent with recent literature 15,33,34 , and also corroborated by our finding that 45% of the net TWS decline in Sahara/Arabia and Central Eurasia is attributed to groundwater depletion (Table 1).
In addition to correlations, we perform another experiment to test whether the times series of each climate index may produce a linear trend which is similar to that of endorheic TWSA. Given any endorheic zone, we simply scale the time series of each index to that of TWSA using a polynomial regression. We then calculate the linear trend and its 95% CI for each scaled index series, and compare them with the trend and CI for the GRACE-observed TWSA.
Except endorheic Australia, we found that no scaled climate index, for any endorheic zone, produces a trend that is statistically similar to the observation (that is, the difference between scaled and observed trends insignificant to the CIs).
As major findings of this paper, we show that the recent decadal decline in endorheic TWS is substantial, and it stems from comparable losses in all three primary water storage components (surface water, soil moisture, and groundwater). This partitioning, together with zonal comparisons with climate variables (e.g., precipitation and temperature in Fig. S13), implies that the endorheic TWS decline is a complicated result of climate variability and forced responses (e.g., to anthropogenic climate change and direct human activities). Except the partitioning among water storage components, a definite decoupling between the contributions of climate variability and forced responses is beyond the scope of this study. Therefore, our reported TWS trend may not imply a continuous secular signal that can be extrapolated into the past or future. We also acknowledge that the current GRACE era prohibits a secure linkage of available TWS observations to any lower-frequency climate variability or secular trend. Nevertheless, our exploration does indicate that compared to exorheic basins, TWS in the continental interiors (endorheic basins) is overall less sensitive to some prominent modes of climate variability at interannual timescales. This comparison highlights the important influences of human water management and longer-term climate conditions (e.g., (multi)decadal variability and anthropogenic warming), in addition to interannual variability of the natural climate system.
It is entirely possible that the recent decadal trend in global endorheic TWS may be altered by future climate variability. But the core message of our findings undeniably conveys that, even though the total water mass in the arid endorheic system may be limited, its change at least at the studied decadal timescale could be so substantial that it largely outpaced that in the entire exorheic landmass and signified a potential impact on the sea level budget. From this considerable decadal decline, we hereby suggest to the geoscience community the importance of a continued monitoring and understanding of TWS changes in the global endorheic system. A meaningful start is to look into how different water storage components contribute, as we demonstrate in this work.   Water levels are paired with mapped areas with a temporal tolerance of 8 days. Levels are averaged if more than one level measurements are available available within any 8-day window. Horizontal error bars represent the scales of level errors reported in the source altimetry data. The labeled hypsometric equation is fitted from the plotted level-area pairs, so it may not be suitable for extrapolating or representing hypsometry of the same lake beyond the plotted levelarea ranges. For the same reason, we do not account for mapping-induced volume uncertainties arising from below the minimum water level during the study period (see Methods section "Estimating lake storage changes"). (b) Monthly time series of estimated lake storage anomaly (gigatons of water) (i.e., departures from the period mean). When no altimetry level is available for an 8-day window, water areas mapped from high-quality MODIS images (if available) are used to inverse water levels using the calibrated hypsometric curve. All measured/inversed water levels are used to derive the 8-day volume anomaly time series, which is further averaged to generate the monthly anomalies. (c) Monthly time series of deseasonalized lake storage anomaly. Vertical error bars and shades represent 95% confidence intervals (two standard deviation uncertainties) for water storage anomalies and best-fit linear trends, respectively. See "Methods" for more detailed uncertainty propagation.          For each enumeration region, the left panel illustrates the scatter plot between the net change in equivalent water thickness (EWT) and total precipitation for each year (blue dots). The relationship is evaluated on an annual basis in order to remove the correlation due to seasonal impacts. Considering that the initial GRACE observation was available in April 2002, we define a full annual cycle as starting from April and ending in March of the next year. r shows the Pearson correlation coefficient between annual net EWT changes (integrated variations within each annual cycle) and total precipitation, and asterisks indicate significant correlations with α = 0.05. Positive correlations are observed in the global endorheic system and most continental zones except a weak negative correlation in Sahara and Arabia (g). To test the pattern consistency, we recalculate the relationship for each region by using net EWT changes and total precipitation in each full calendar year (January to December for 2003-2015; black squares). The general relationships remain similar, except that in Sahara and Arabia, the previous negative correlation becomes weakly positive, implying a likely minor impact of precipitation on its recent TWS changes (which is consistent with the dominant groundwater contribution seen in Fig. 4). For each region, the red horizontal line approximates the period-average precipitation level on which EWT change is 0. Theoretically, this precipitation level meets the demand for evapotranspiration, and maintains the overall water budget equilibrium in the endorheic region. Precipitation below or near the equilibrium line may lead to possible water budget deficits. Right panels illustrate monthly temperatures and 95% uncertainties from April 2002 to March 2016. Quantified temperature trends confirm that warming has prevailed in each of the endorheic zones, with an average rate of 0.06 K yr -1 exceeding the global average of approximately 0.02 K yr -1 39 .       Table S6.     Table S5. Surface water storage changes and contributions of different storage components in global and zonal endorheic basins. Uncertainties represent 95% CIs (two standard deviations). Table S6. Validation of monthly soil moisture time series from GLDAS model ensemble using SCAN in situ soil moisture measurements in endorheic North America. Validation includes Pearson correlation coefficient (r) and interannual trend (% yr -1 ) between deseasonlized monthly time series of modeled and in situ soil moisture contents (%), for a total of 22 SCAN stations (see Fig. S19 for locations) where measurement records exceed 5 consecutive years. Trend uncertainties represent 95% CIs (two standard deviations). The GLDAS model ensemble for total soil layers includes VIC, Noah, CLM and Mosaic, and the 0-10 cm ensemble only includes VIC, Noah, and CLM as the Mosaic model does not contain the 0-10 cm soil layer.
Note: SCAN in situ soil moisture data sets can be accessed at www.wcc.nrcs.usda.gov/scan. The symbols ** and * represent pvalues that are less than 0.01 and 0.05, respectively.