Cascading Dynamics of the Hydrologic Cycle in California Explored through Observations and Model Simulations

: As drought occurs in a region it can have cascading effects through the water cycle. In this study, we explore the temporal co ‐ evolution of various components of the hydrologic cycle in California from 2002 to 2018. We combine information from the Gravity Recovery and Climate Experiment (GRACE) satellites, the North American Land Data Assimilation System (NLDAS) suite of models, and the California Department of Water Resources (DWR) reservoir levels to analyze dynamics of Total Water Storage (TWS), soil moisture, snow pack, large reservoir storage, and ultimately, groundwater. For TWS, a trend of − 2 cm/yr is observed during the entire time period of our analysis; however, this rate increases to about − 5 cm/yr during drought periods (2006 − 2010 and 2012 − 2016). Results indicate that the majority of the loss in TWS is caused by groundwater depletion. Using proper error accounting, we are able to identify the start, the peak, and the ending of the drought periods for each individual water state variable in the study domain. We show that snow and soil moisture are impacted earlier and recover faster than surface water and groundwater. The annual and year ‐ to ‐ year dynamics shown in our results portray a clear cascading effect of the hydrologic cycle on the scale of 8 − 16 months


Introduction
In recent decades, the demands of fresh water for urban, agricultural, and environmental uses have exceeded the natural renewable supply in many regions of the world, especially for arid and semi-arid regions such as California [1,2]. To date, this gap between the available surface water supply and the growing water demand has caused significant stresses on fresh water systems globally. This has created a large reliance on groundwater, a slowly renewable resource conveniently stored at shallow subsurface aquifer layers. Because of continued extraction and insufficient recharge that is complicated by climate variability (i.e., droughts), depletion of subsurface aquifers has become a dire environmental consequence that is being observed globally in arid and semi-arid regions with highly variable precipitation [e.g., 3−12]. Generally, various aspects of the hydrologic cycle, such as snow pack and soil moisture, have been impacted by increasing temperatures and more variable precipitation. This overall shock to the water cycle can be seen in both observations and model simulations of total water storage for many regions of the world [7−12].
Despite the critical need for fresh water supplies, water resource availability is often poorly monitored. For example, groundwater monitoring networks do not exist at the same scope and scale for those that track surface water [13]. Conventionally, groundwater level change is monitored through in situ measurements at observational wells, which are sparse, very cost ineffective, and difficult to access. Furthermore, measurement records are often discontinuous. To overcome these challenges, various methods, such as geostatistical interpolation and hydrologic model simulation, have been proposed to estimate fresh water availability [14,15]. Even though hydrologic models may be capable of inferring water availability at high spatiotemporal resolution, the conventional land surface models, such as Noah-MP [16], omit consideration of human activities. This lack of representation could cause biases in simulations, for example by underestimating rates of groundwater loss [17]. In addition to groundwater, other sources of fresh water are hardly monitored or tracked, such as soil moisture or snow water. Although methods exist to gather data on these processes, such as airborne or satellite observations, in situ observations, or methods like data assimilation that fuse observations with model simulations [15,18,19], detailed knowledge about the state of these systems is still lacking. This causes a great deal of uncertainty in estimates of the different aspects of the hydrologic cycle.
In California, particularly in productive agricultural regions such as the Central Valley, has experienced periodic declines in total water storage in recent decades. A main reason for the loss of fresh water in this region is due to groundwater depletion. Annually, at least 40% or more of the Central Valley's demand for water is supplied from groundwater, with the majority is used in agriculture [11,20] Data from the United States Geological Survey (USGS) show that between 1962 and 2003, groundwater use in the Central Valley increased from 0.75 to 2.5 km 3 per year [21]. The groundwater extracted for irrigation often exceeds the natural recharge, leading to declines in the groundwater table [13,21−23]. Furthermore, the winter snow pack in the Sierra Nevada is usually a reliable source of fresh water for the drier spring and summer seasons. However, recent years have shown a large increase in temperatures and a slight decrease in precipitation, causing a depletion of the snowpack for many regions where snow is usually available [15]. Also, large surface water reservoirs have shown periods of depleted reserves during droughts, and many of the largest reservoirs were in critical condition during recent droughts [24]. The impacts from droughts are usually felt in one hydrologic process before it reaches another, for example, depletion of the snow pack occurs before groundwater reserves are relied on. In this study, we aim to capture this cascading effect of the hydrologic cycle for California, and to quantify its uncertainty by analyzing different sources of information that are available.
More recently, much effort has been focused on monitoring aspects of the hydrologic cycle using remote sensing techniques [e.g., 4,5,7,13,25−27]. For example, NASA's Gravity Recovery and Climate Experiment (GRACE) mission [28] has proven to be an extremely valuable tool for regional to global scale water cycle studies [7,28]. Satellite observations of Earth's gravity field from GRACE are processed routinely into estimates of surface mass change and can provide information about basinscale dynamics of groundwater. GRACE mass change estimates can be combined with other hydrologic information, such as model simulations or in situ observations, to infer groundwater changes [e.g., 7,[11][12][13]30]. For example, information on snow and soil moisture can be obtained from high-resolution land surface models, and reservoir storage information can be obtained by local water managers, and this information can be helpful to estimate groundwater dynamics derived from the GRACE satellites. In this way, various sources of information are fused to quantify the water cycle along with its uncertainty.
Total Water Storage (TWS) estimates from GRACE include all of the snow, ice, surface water, soil water, canopy water, and groundwater in a region, and when combined with auxiliary hydrologic datasets, TWS can be utilized to infer process information on groundwater, or other water sources. To investigate the various water processes in our study, we utilize data from the North American Land Data Assimilation System (NLDAS-2) for outputs of soil moisture and snow water equivalent [19,31]. For surface water, we use data from the California Department of Water Resources (DWR) on the 50 largest reservoirs in the state. We ignore smaller terms that may be negligible, such as ice or canopy water. This leaves us with the groundwater component. By combining all of these sources of information, we can investigate the temporal evolution of each water process and answer questions such as: 1) Is there a cascading effect in the hydrologic cycle that can be clearly extracted from the data? 2) Are there obvious differences in the timing of drought (beginning, peak and end of drought) for each process, and if so, by how many months (or years)? By utilizing the various sources of information available to us via satellite observations, model simulations, and ground observations, we are able to infer the cascading effect of the water cycle and can make quantitative assessments of the evolution of each hydrologic process.

GRACE Satellites
NASA's Gravity Recovery and Climate Experiment (GRACE) satellite mission provided a new means to observe fresh water variability from space [e.g., 5,7,13,29]. Since 2003, satellite observations of time variable gravity from the GRACE mission have provided monthly estimates of total mass change globally. These estimates have been widely used for hydrology and water resource monitoring.
Various recent studies have demonstrated that GRACE derived estimates of variations of total water storage (TWS) can provide fresh water availability estimates with sufficient accuracy [11,18,26]. These GRACE-based methods have been applied to regions such as Northern India [5,32], the Middle East [33,34], Northern China [35,36], and California [7,11,13,30], among others. In this study, the JPL RL05M.1 Version 2 GRACE mascon product, which can be downloaded at https://grace.jpl.nasa.gov/data/get-data/ [37,38], is used to examine water storage changes in California ( Figure 1) and its underlying aquifer system during the GRACE satellite period, from April 2003 through June 2017. Figure Figure 1, especially during droughts. Furthermore, the GRACE product also provides uncertainty estimates, which will be used in our overall assessment of groundwater uncertainty described in Section 2.5. In other words, the standard deviation of the TWS estimates at each time step, or TWS (t), is provided with the downloaded GRACE product.

NLDAS-Mosaic, Noah, and VIC
The North American Land Data Assimilation System (NLDAS-2) constructs quality-controlled, and spatially and temporally consistent, land-surface model (LSM) datasets from the best available observations and model outputs [19,31] Data sets are available on hourly, daily, monthly and climatological time scales. NLDAS products are used in modeling, research, and applications, such as drought and flood monitoring, watershed and water quality management, and case studies of extreme events. NLDAS integrates a large quantity of observation-based and model reanalysis data to drive offline (not coupled to the atmosphere) land surface models (LSMs), and is executed at 1/8thdegree grid spacing over North America.  Figure 2 shows temporal trends in SM and SWE for 2003-2018 simulated by the suite of NLDAS-2 models. As shown in Figure 2, there are large discrepancies between the various model simulations, which create uncertainties for our estimates of SM and SWE. We define the uncertainty for each process as the standard deviation (std) between the ensemble of three models at each time step ( Figure 3B,C):

Reservoir Data from the California DWR
The California Department of Water Resources (DWR) is a government agency that is responsible for managing and protecting California's water. The DWR provides the California Data Exchange Center (CDEC), which installs, maintains, and operates an extensive hydrologic data collection network including data on the largest reservoirs in California, which we use in this study. The data can be downloaded at http://cdec.water.ca.gov/reservoir.html. We extract the information for the 50 largest reservoirs in California for the period 2003-2018, and use this as our Surface Water (SW) portion of the hydrologic cycle for this study. As done in previous studies [e.g., 7], we apply 15% noise to the reservoir level data to obtain uncertainty estimates for surface water ( Figure 3D), i.e., SW (t) = 0.15 × ( SW (t) )

Data Processing to Calculate Groundwater
Each data product was re-gridded to match the GRACE product resolution (GRACE resolution = 0.50° × 0.50°). The original NLDAS simulations had a resolution = 0.125°, but these were re-gridded using bi-linear interpolation to match the gridding of the GRACE product. Once all products were on the same grid, we performed an area-weighted average that accounted for grid areas as a function of latitude. Then each of the data pieces was accumulated (summed) to obtain a temporal curve with uncertainties for each process (i.e., TWS(t), SM(t), and SWE(t)). The DWR reservoir data were in the form of point information, but we added up all these values to get one temporal data vector for California (i.e., SW(t)). Finally, once all the data was in the form of a temporal vector (with uncertainties), we were able to estimate the groundwater and its uncertainty.

Groundwater Estimates and Associated Uncertainties
The GRACE products estimate total mass change; therefore, the contributions of mass changes from other processes must be estimated before the groundwater mass change can be inferred [7]. Typically, this is done with a numerical land surface model. However, one standing issue with respect to GRACE is how much of the derived groundwater mass change is due to anthropogenic causes versus model error. Addressing this requires a rigorous and robust treatment of uncertainties that come with individual model corrections. Below, we outline the essential steps for deriving groundwater storage change and the approach we adopt to quantify the associated uncertainties.
The GRACE total water storage (TWS(t)) data were processed in a way that was able to quantify the uncertainty of the temporal dynamics of groundwater storage change aggregated for California ( Figure 3A). We used information for soil moisture (SM(t)), snow (SWE(t)) and surface water (SW(t)) aggregated for the GRACE analysis region as described in Sections 2.2 and 2.3, and then subtracted their effects from the TWS variations (TWS(t)) obtained from GRACE, where t indicates time. Our analysis was performed at monthly resolution. Therefore, the time-specific changes in groundwater (GW(t)) for a region are given by the following: We applied the Monte Carlo method and draw n = 1000 ensembles by introducing the uncertainty from each hydrologic process (i.e., TWS, SM, SWE, and SW). The calculation uses a normal (Gaussian) error with standard deviation , or N( 0 ,  ): The uncertainty of TWS, or TWS in Equation (4), is directly the error associated with the GRACE measurements, which is offered as part of the data product. The uncertainty information for SM and SWE, or SM and SWE in Equation (4), was obtained from the NLDAS-2 models. We characterized the uncertainty of SM and SWE as the standard deviation (std) between the ensemble of three models at each time step from the Mosaic, VIC, and Noah land surface models. The data for the reservoirs that represent the surface storage were obtained from the California Department of Water Resources (DWR). The uncertainty of the surface reservoir storage, or SW in Equation (4), was assigned as 15% of the magnitude of the storage, a conservative approach that has been employed in previous studies [e.g., 7].
After propagating all the sources of uncertainty, we were able to estimate the temporal groundwater dynamics for California with uncertainties (see Figure 3E). Note that there should theoretically be a correlation between the error components, since these processes are directly linked with one another. One way to account for this is to include these correlation terms (if they were known) in the non-diagonals of the error matrices. Although accounting for these correlations is beyond the scope of the current paper, we believe it could provide further insight on how these processes co-evolve together.

Wavelet Analysis
Many scientists have made use of the wavelet method in analyzing time series in geophysical or climate studies [39,40]. At present, there are many easy-to-use wavelet software packages for analyzing two time series together. Here, we used the cross wavelet transform and wavelet coherence (i.e., https://www.mathworks.com/matlabcentral/fileexchange/47985-cross-wavelet-and-waveletcoherence) for examining relationships in the time frequency domain between two time series [40]. Grinstead et al. demonstrated how phase angle statistics can be used to gain confidence in causal relationships and test mechanistic models of physical relationships between time series. As hypothesized, the annual and year-to-year dynamics of the hydrologic cycle should portray a clear cascading effect. To investigate this further, we applied a wavelet analysis, which shows the frequency (in months) with which there is a clear signature between the temporal correlation of two processes. We applied this test to show the wavelet analysis results for each pair of processes estimated in our study.  Figure 2, trends in SM and SWE simulated by the NLDAS-2 suite of models are shown. For SM (Figure 2A−C), there are significant differences in the simulated trends between the various models, but losses of water from SM are estimated to have been roughly -1 (+/-0.25) cm/yr, and as high as −3 cm/yr in some regions. For SWE ( Figure 2D−F), most models agree that there has been roughly -1 (+/− 0.10) cm/yr loss of fresh water from the snow pack.

Trends in the Water Cycle
After applying the method described in Section 2.5, we are able to estimate the groundwater anomalies for the period of the study. Figure 3 shows the timeline of each hydrologic process, along with their respective uncertainty bounds. Table 1 reports the signal of water loss from each process, for different time periods. These numbers are also displayed in Figure 4, where the processes that contribute to the trends seen in TWS are disentangled and displayed separately. In summary, these results indicate that the majority of the loss in TWS is caused by groundwater depletion. This is especially apparent when we look at the volume of water lost in the right column of Table 1, where groundwater makes up the majority of this signal. Also, in Figure 4E−G, we show that groundwater makes up 60%−90% of the signal seen in TWS loss (red histograms).

Cascading Effect of Hydrologic Cycle
In this section, we aim to capture the cascading nature of the hydrologic cycle. Figure 5A shows the mean annual cycle of the monthly averages for each process. We see that SWE and SM have peaks and also minimums that occur earlier than surface water and groundwater. For SWE and SM, the peaks occur around winter time, and the minimums around summer. For surface water and groundwater, the peaks occur around the spring, and the minimums occur in the fall. The cascading effect of the water cycle on an annual scale is therefore depicted here in these figures. Figure 5B-F illustrates the de-seasonalized data, which we use to estimate the start, the peak, and the ending of the drought periods for each hydrologic process. This processing step is done to change the raw simulations into ones that show deviation from the annual cycle. In other words, if a certain process has a low or high value relative to the same month in the historic record, this processing step captures that deviation. This also allows us to investigate all of the processes as anomalies, which matches the forms of TWS and GW. These deviations, shown with black dots in Figure 5B-F, are computed by subtracting the mean annual cycle shown in Figure 5A from the mean simulations shown in Figure 3. A 13-month running mean (blue dots) show the state of each water process relative to its climatology for that month while smoothing out high-frequency fluctuations. This allows us to understand whether a certain process has low or high levels of water storage relative to its past states, and thus to quantify the timing of the beginning, peak, and ending of the drought for each process. , which shows each water process normalized by its interannual variability. A 13-month running mean (blue dots) show the state of each water process relative to its climatology for that month. This allows us to investigate the beginning, the peak, and the ending of each drought period.
Our estimated start, peak, and ending of each drought period for each hydrologic process are displayed in Table 2. Also included are the peak magnitudes of volumetric water loss (in km 3 ) for each component. These results indicate that snow and soil moisture are impacted earlier than surface water and groundwater. By initial inspection using the values in Table 2, we can estimate that there is a lag of roughly 8-16 months between the impact on SM and SWE compared to that of SW and GW. The next section highlights another way to quantify this temporal lag correlation between the various processes through wavelet analysis.

Results from the Wavelet Analysis
A B C D E F As mentioned in the previous section, the annual and year-to-year dynamics shown in our results portray a clear cascading effect of the hydrologic cycle on the scale of 8−16 months. To investigate this further, we apply a wavelet analysis which shows the frequency at which there is a clear signature between the temporal correlations of two processes. Figure 6 shows the wavelet analysis results for each pair of processes that we estimate in our study. The time series used for each process in the wavelet analysis shown in Figure 6 are the mean state estimates from Figure 3 (blue lines). We also applied this analysis on the 13-month running means from Figure 5B−F (blue dotted lines), but do not show the wavelet analysis for these time series since there were no coherent relationships between them. This is possibly due to the loss of information when applying a running average.
The panel on the top right of Figure 6 shows that, for TWS and groundwater, there is a high power signal (yellow colors), suggesting there is high correlation at around 8−16 months between these two processes. Outside of this time window, e.g., at 4 or 32 months, there is a low power signal (blue colors), which indicates there is little correlation between the two processes. the frequency of correlation at around 8−16 months is observed between all pairs of hydrologic processes. This indicates that the evolution of the hydrologic cycle occurs through an 8-16 month window. This window presumably begins with the rise and fall of soil moisture and snow pack in the winter and spring, followed by the decline of surface water and groundwater. These sources of water will continue to be relied upon and depleted in the summer time until the next season with sufficient rain begins, which could be the next or even the following Fall. This is what describes the pulse of roughly 8−16 months for the hydrologic cycle in California.
The correlation seen in the 8-16 month time range is a combination of both natural as well as human factors. First, the period of the annual rainfall accumulation affects other processes, like the accumulation of snow, soil moisture, and reservoir storage. Once the rainfall is finished for the year and the natural processes take their course, such as snow melting and soil moisture evaporation, the next step of this cascading evolution is the human use of water. By releasing reservoir storage and pumping groundwater to meet the demands of the human population, the final aspect of this cascading hydrologic cycle is complete as water reserves are depleted. The cycle starts again with the new rainfall of the next year. This roughly explains the 8-16 month window that is captured with the wavelet analysis, which is a result of both natural and human factors.

Summary and Discussion
In this study, we explored the temporal co-evolution of various aspects of the hydrologic cycle in California from 2002-2018. We combined information from the GRACE satellites, the NLDAS suite of models, and the California DWR to analyze dynamics of Total Water Storage, soil moisture, snow pack, large reservoirs, and finally, groundwater. We calculated groundwater after propagating uncertainties from other hydrologic processes and factored them out of the TWS estimates from GRACE.
There were clear trends of losses in water volume for all hydrologic processes. For TWS, a trend of -2 cm/yr was observed during the entire time period (Figure 1), and −5 cm/yr during drought periods (2006−2010 and 2012−2016). We disentangled the various processes from this trend, and found that the majority of the loss in TWS is caused by groundwater depletion (Figures 2 and 4; Table  1). We quantified the start, the peak, and the ending of the drought periods for each hydrologic process, and showed that snow and soil moisture are impacted earlier than surface water and groundwater ( Figure 5; Table 2). The annual and year-to-year dynamics shown in our results portray a clear cascading effect of the hydrologic cycle on the scale of 8−16 months (Figures 5 and 6). Overall, this study displayed the cascading evolution of the hydrologic cycle in California. This timeline begins with the accumulation of soil moisture and snow pack in the winter and its release in the spring, followed by the decline of surface water and groundwater during the drier summer period. Studies like this can help scientists and water managers better define and characterize droughts and their impacts. Generally, drought is defined by the lack of precipitation or the depletion of visible reserves such as reservoir storage or low snowpack. This study showcased the lagged effect of a drought, which also includes the depletion of groundwater reserves. It can be a mistake to assume the end of a drought period is reached as soon as rainfall occurs, since other aspects of the hydrologic cycle may not be replenished for months or years after. In this study, we hoped to portray that all aspects of the hydrologic cycle should be considered when characterizing a drought, and the cascading effect of this cycle can be confidently quantified, as we did here.

Acknowledgments:
The research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). Copyright 2020. All data used in this study are publicly available. The GRACE data set is available at https://grace.jpl.nasa.gov/data/get-data/, the NLDAS model output data are available at https://disc.gsfc.nasa.gov/datasets?keywords=NLDAS, and the DWR dats can be found at http://cdec.water.ca.gov/reservoir.html. The software package for wavelet analysis can be found at https://www.mathworks.com/matlabcentral/fileexchange/47985-cross-wavelet-and-wavelet-coherence.

Conflicts of Interest:
The authors declare no conflict of interest