The Role of Atmospheric Rivers on Groundwater: Lessons Learned From an Extreme Wet Year

In the coastal regions of the western United States, atmospheric rivers (ARs) are associated with the largest precipitation generating storms and contribute up to half of annual precipitation, but the impact of ARs on the integrated hydrologic cycle, specifically on groundwater storage and hydrodynamics, is largely unknown. To better explore the hydrologic behavior of AR versus non‐AR event precipitation, we present a novel combination of two water tracking methods (one in the atmosphere and one in the subsurface) to explicitly track the full lifecycle of water parcels generated by ARs. Simulations of northern California's Cosumnes River watershed during the record wet 2017 water year are performed via the coupling of a high‐resolution regional climate model and a land surface‐groundwater model accounting for lateral groundwater flow. Despite ARs contributing more precipitation than non‐AR storms, we find less AR water is preferentially stored in aquifers by year end. Fractionally, ARs result in 300% less snow derived groundwater‐recharged compared to non‐AR precipitation. Rain‐on‐snow (RoS) plays an important role in AR‐driven discharge, where over 50% of total discharge from ARs snow is from RoS events. Finally, despite record‐breaking annual precipitation, simulated groundwater depletion occurs by year end due to estimates of groundwater pumping activities. The results from these simulations serve as a partial analogue of future hydrologic conditions where ARs are expected to intensify and provide a greater fraction of annual precipitation due to climate change.

(i.e., a higher fraction of precipitation which results in streamflow). ARs are also known to bring about warmer atmospheric conditions, resulting in higher snow levels and more rain than snow compared to non-AR precipitation events . This change in precipitation partitioning may directly impact groundwater recharge, because snowmelt more effectively infiltrates into the subsurface when compared to rainfall (Earman et al., 2006). Changes in precipitation phase partitioning may also lead to more rain-on-snow (RoS) events, also associated with AR conditions (Chen et al., 2019;Guan et al., 2016;Michaelis et al., 2022), potentially exacerbating the effects of less snowfall on groundwater recharge. Seasonally, increased occurrence of ARs may indicate a sharpening of the winter precipitation season, so that peak hydrologic water availability would occur earlier in the water year when potential evapotranspiration is lower, potentially resulting in less overall water flux back to the atmosphere.
Very few studies have investigated the role in which ARs play on hydrologic processes, and those which do solely focus on above-ground hydrology. Konrad and Dettinger (2017) found that high integrated vapor transport (IVT) conditions associated with ARs coincided with the highest streamflow events between 1945 and 2015. Chen et al. (2019) similarly found that ARs produced higher streamflow than non-AR events, where runoffs ratios were double that of non-AR runoff ratios if snowpack existed before the onset of the AR. They also found that RoS events during ARs are mostly driven by changes in air temperature (and to a lesser degree radiation changes) so that ARs suppress evapotranspiration.
Studies which assess changes to subsurface hydrology are even less prevalent in the literature, in part due to limitations in tracking methods, both in observational networks and numerically in model frameworks. This is despite observed relationships between antecedent soil moisture and AR-related discharge (Haleakala et al., 2022;Ralph et al., 2013). One study developed a water tracer model within the Noah LSM and found that a large fraction of precipitation from a simulated AR in the Pacific Northwest was retained in the soil after 6 months (Hu et al., 2018). However, the influence of ARs on deeper groundwater recharge rates and water storage remains a void in the hydro-meteorological and water resource communities. To the best of our knowledge, this study is the first to examine the role in which ARs contribute to groundwater dynamics.
Here, we combine two novel water tracking algorithms-above and below ground-to explore how ARs influence the integrated hydrologic cycle, with specific foci on groundwater recharge and subsurface flow. AR and non-AR precipitation is tracked across the entire atmosphere through bedrock continuum during an extreme wet year characterized by multiple ARs. We employ cutting-edge atmospheric and hydrologic simulations in high resolutions over an unimpaired watershed in northern California characterized by natural river flows. The novel coupling of above-and below-ground AR water tracking enables an examination of the full lifecycle of water. Information on AR contributions to groundwater (including water residence times, depth, and recharge timing), water source (rain vs. snow), and eventual exit pathway (evapotranspiration or discharge) are therefore all possible in this approach. This study builds on a wealth of developments in AR-tracking capabilities (Inda-Díaz et al., 2021;Rutz et al., 2019;Shields et al., 2018;Zhou et al., 2021). Similarly, subsurface tracking capabilities also have a rich history (Salamon et al., 2006) which were only recently extended to include surface processes (Bearup et al., 2016;de Rooij et al., 2013;Maxwell et al., 2019;Yang et al., 2021) with relatively few applications (Bearup et al., 2016;Maxwell et al., 2019;Wilusz et al., 2020;Yang et al., 2021), but have yet to be extended as done in this study to track individual storms or classes of storms. The cross-disciplinary approach employed here enables a physics-based understanding of the role in which atmospheric processes such as ARs have in shaping subsurface flow processes at local-and watershed-scales, relevant for decision makers and stakeholders who are adapting to rapidly changing environmental conditions . This study also serves as a first demonstration of individual storm tracking from atmosphere through bedrock to understand water partitioning behavior which may be important to understand in a no-analogue future climate.
In addition to demonstrating these modeling tools in a new, storm-based atmosphere to groundwater framework, we seek to study the role of ARs on hydrology and to gain key lessons learned by studying a recent, extreme winter in the context of the following science questions: • Does water derived from AR and non-AR precipitation partition differently? • Specifically, does the character of AR precipitation events impact groundwater storage differently than non-AR precipitation? Considerations include more prolonged and higher precipitation totals, impacts from rain-snow events, and different spatial distributions of precipitation. • And finally, what are the implications on groundwater storage dynamics and water resources?

Site Description and Selected Water Year
The Cosumnes watershed ( Figure 1a) located in northern California, spans the western slope of the Sierra Nevada and terminates in the Sacramento Central Valley, similar to many watersheds in California that supply water down gradient for domestic, agricultural, or industrial water users across the state. The Cosumnes River is unique in that it is one of the last major rivers in the western United States without a major dam, enabling the rare study of uninhibited natural river flows. The Cosumnes supports some of the last remaining wetland and riparian habitat in the Sacramento Central Valley as well as a burgeoning viticultural economy and other agricultural production. An integrated hydrologic model of the Cosumnes watershed has been developed with the ParFlow-CLM model, and used to explore a myriad of above-and below-ground watershed dynamics including drought whiplash (Maina, Siirila-Woodburn, Newcomer, et al., 2020), wildfires (Maina & Siirila-Woodburn, 2019) modeling impacts due meteorological forcing resolution (Maina, Siirila-Woodburn, & Vahmani, 2020), and projections of end-of-century hydrologic states (Maina et al., 2021). A series of model validation exercises have been conducted to test the model ability to represent key hydrologic processes. Model performance has been compared to point measurement of observed streamflow and groundwater levels (see Maina, Siirila-Woodburn, Newcomer, et al., 2020; Figure 3) over recent wet and dry periods in the region. We have also conducted spatially distributed comparisons of key hydrologic variables (evapotranspiration, soil moisture, and snow water equivalent) using a variety of remote sensing products including METRIC, SMAP, SNODAS, and a MODIS-derived snow product, ParBal (see Figure 4 of Maina, Siirila-Woodburn, Newcomer, et al., 2020). Further description of the ParFlow-CLM model is described in Section 3.2.
Water year 2017 (hereafter WY2017, which spans 1 October 2016 to 30 September 2017) was selected for this study as it was the most active AR season since, at least, the mid-twentieth century (Gershunov et al., 2017). In the Cosumnes, WY2017 resulted in the highest cumulative streamflow in a >100 years record (Figure 1c), providing an opportunity for "lessons learned" from an end-member precipitation year. Although a larger sample of ARs is most ideal to assess the statistical representation of individual storms, WY2017 provides an example of different storm types typical of AR conditions in the region, and the subsequent water partitioning associated with those storms.

Methods
A combination of high-fidelity, high-resolution models and tracking techniques are used here to study the movement of water across the bedrock-through-atmosphere interface. This approach (a conceptual model is shown in Figure 2) includes feature tracking in the atmosphere to identify ARs (described in Section 3.1) and discrete particle tracking of parcels of water in the subsurface (described in Section 3.2).

Atmospheric Simulation and AR Tracking
The Weather Research and Forecasting model (WRF;  is a state-of-the-art, fully compressible, non-hydrostatic, mesoscale numerical weather prediction model. A nested domain configuration of WRF (version 3.6.1) is used to downscale NLDAS-2-based forcings to 500 m resolution. Four domains with horizontal resolutions of 13.5, 4.5, 1.5, and 0.5 km are used with 30 vertical atmospheric levels. WRF is initialized using post-spin-up soil moisture from ParFlow-CLM. Other WRF initial conditions are based on NLDAS-2 as described in Maina, Siirila-Woodburn, and Vahmani (2020). The physics parametrizations used in WRF include the Dudhia scheme (Dudhia, 1989) for shortwave radiation, the Rapid Radiative Transfer Model (Mlawer et al., 1997) for longwave radiation, University of Washington Boundary Layer Scheme (Bretherton & Park, 2009) for the planetary boundary layer, the Morrison double-moment scheme (Morrison et al., 2009) for microphysics, and the Eta Similarity scheme (Monin & Obukhov, 1954) for the model surface layer. The Grell-Freitas cumulus parameterization (Grell & Freitas, 2014) is used only for domains d01 and d02. The adopted WRF configuration has been extensively validated in the literature (Vahmani & Jones, 2017;Vahmani et al., 2019). These studies show a reasonable performance for WRF over California, reproducing daily air temperature and evapotranspiration with RMSDs of 1.1°C and 0.74 mm/day, respectively. One-way coupling of the highest resolution (0.5 km) WRF output to ParFlow-CLM was recently done as part of an analysis to quantify forcing resolution biases on integrated hydrologic processes (Maina, Siirila-Woodburn, & Vahmani, 2020) and is extended for use here to track ARs. A comparison of spatially distributed, high-frequency precipitation and temperature in the Cosumnes was performed for the WY2017 forcing used in the study, and is detailed in Appendix 4 of Maina, Siirila-Woodburn, and Vahmani (2020), showing good agreement with observations. TempestExtremes is used to track landfalling ARs during WY17 in the outermost (d01) domain of the WRF simulations (Ullrich & Zarzycki, 2017;Ullrich et al., 2021). To track the ARs we calculate IVT at hourly resolution using meridional and zonal wind speeds and specific humidity across all 30 vertical levels. The AR tracking parameter settings, specifically SplineARs in TempestExtremes, are customized to the spatiotemporal resolution of the WRF simulations (informed by prior research in Rhoades, Jones, Srivastava, et al., 2020;Rhoades et al., 2021). Parameter settings for TempestExtremes are described in Table 1.

Integrated Hydrologic Flow and Particle Tracking
The integrated hydrologic model ParFlow coupled to the land surface model the Common Land Model (CLM) is used here to simulate the continuum of hydrologic conditions of a representative hillslope spanning from bedrock Figure 2. Conceptual image of the coupled atmosphere feature tracking scheme used to identify atmospheric rivers (ARs) and connection to subsurface particle tracking, including main processes considered.

Parameter
Variable to the lower atmosphere. ParFlow simulates variably saturated and fully saturated flow in three dimensions via the Richards' equation (Richards, 1931) and overland flow in two dimensions via the kinematic wave equation (Ashby & Falgout, 1996;Jones & Woodward, 2001;Kollet & Maxwell, 2006Maxwell & Miller, 2005). CLM simulates a coupled water energy balance at the surface layer of the domain by incorporating spatially distributed vegetative processes across specified land use types (Dai et al., 2003). CLM also includes a sophisticated, multi-layer snow model (which includes thermal processes, canopy interception/throughfall, sublimation, snowpack metamorphism and compaction, albedo decay, melt, and snow age) and is further described in detail by Ryken et al. (2020). Land cover is based on the 2011 NLCD Land Use and Land Cover database and is parameterized by the International Geosphere-Biosphere Program (IGBP) database (Homer et al., 2015;International Geosphere,Biosphere Programme, n.d.). Figure 1a shows the distribution of land use and land cover in the Cosumnes watershed; Figure 1b shows the elevation distribution and topography, spanning from the headwaters in the Sierra Nevada to sea level in the Sacramento Central Valley.
The coupled system is forced by eight meteorological variables at an hourly temporal resolution: downwelling long-and short-wave radiation, precipitation, 2-m surface air temperature, 2-m east-to-west (zonal) and north-to-south (meridional) surface winds, atmospheric pressure, and specific humidity, which are supplied via a one-way coupling with the WRF output (Section 3.1). A complete description of the Cosumnes ParFlow-CLM configuration, including input datasets and validation with observations and other remotely sensed datasets can be found in Maina, Siirila-Woodburn, Newcomer, et al. (2020). The spin-up condition to obtain the initial condition, which consists of a one-time spatial kriging to obtain a map of initial pressure field then a simulation period of precipitation minus evapotranspiration of average climate forcing, then followed be a final step of recursively running a nearly average transient water year (WY2019) recursively until convergence is reached within 1%, similar to Maina et al. (2021) and Maina, Siirila-Woodburn, Newcomer, et al. (2020).
The Lagrangian particle tracking model EcoSLIM is applied to track water particles starting from infiltration at the land surface. For this application, we tag particles as AR or non-AR precipitation based on the feature tracking done by TempestExtremes (see Section 3.1). EcoSLIM uses spatially distributed output from hydrologic models and forcing fields to insert, remove, and move particles (representative of parcels of water) through the flow domain (Maxwell et al., 2019). Evapotranspiration fluxes, subsurface pressure, subsurface saturation, and the three-dimensional subsurface flow velocity fields are provided by ParFlow-CLM. However, the approach used here is fully compatible with other integrated hydrologic models. Particles are introduced into the system via precipitation fluxes (either snow or rain) or initialized as groundwater. Particle mass is determined based on whether the particle was inserted following a precipitation event or initialized as groundwater; all particles are mass weighted. A full description of the code can be found in Maxwell et al. (2019). The model tracks the water source (i.e., snow, rain, or groundwater) of particles moving within the domain and how they exit the domain through either evapotranspiration or outflow. Corresponding residence times of all particles are also quantified, relying on an optimal local time-step of particle movement based on a prior version of the SLIM and SLIM-FAST codes (Maxwell, 2010)  The spin-up procedure for the particle tracking model follows after that of the ParFlow-CLM spinup, and consists of initializing one particle per cell starting with the pressure head distribution based on the final condition of ParFlow-CLM spin-up. The EcoSLIM spin-up begins with the final step of the ParFlow-CLM spin-up by recursively running the nearly average transient WY2019 until convergence has been reached, in this case around 35 years of simulation time. This ensures that both the velocity field (from ParFlow-CLM) and residence time distribution of the particles (from EcoSLIM) are converged. Following the spin-up, we then compare the evolution of particle mass that exit or remain in groundwater storage during WY2017, partitioned by water source (i.e., rain, snow, or groundwater). Storage of each cell either saturated (Storage Sat ) or unsaturated (Storage Unsat ), is calculated as follows: where is the cell porosity (−), dx, dy, dz are respectively the x-, y-, and z-cell resolutions (−), is the cell pressure head (m), and SS is the specific storage (m −1 ).
RoS events are identified when precipitation falls in the form of rain over cells that have snow water equivalent (SWE) greater than zero; any infiltrating particles over these regions are tracked in EcoSLIM as snow and further classified/tracked as RoS particles. As with all the dynamically added subsurface particles within EcoSLIM, the amount of RoS subsurface particles is determined by the net infiltration at the land surface, so that RoS particles are the result of the water-energy balance within the snowpack resulting in snowmelt, and not the precipitation rate alone.

AR Identification and Rain-Snow Partitioning
The cumulative precipitation total of WY2017, as well as the water partitioning between snow and rain, and AR events as identified via Tempest Extremes, and non-AR events, which is any other precipitation outside the identified AR windows, are shown in Figure 3. The majority of large precipitation events corresponding to sharp steps in the cumulative precipitation time series are identified as ARs (brown vertical boxes), although some exceptions exist. For example, note a precipitation event around WY day 180. WY2017 precipitation totals indicate that 52% of precipitation is provided by ARs. The majority of both AR and non-AR precipitation is in the form of rain and ARs contribute less snowfall than non-ARs: 19% of non-AR precipitation falls as snow, whereas only 13% of AR precipitation is snowfall. Across the 10 storms, all produce a mixture of snow and rain. We reserve any inferences about temporal trends with magnitude, rain-snow-partitioning, etc., given the storm sample size; however, an examination of storm-by-storm behavior is helpful to understand system dynamics.

Individual AR Characteristics and Impacts on Snowpack
Simulated maps of the AR-specific total rainfall and snowfall show an orographic dependence with greater precipitation occurring at higher elevations, in the headwaters of the Sierra Nevada (Figure 4). In the highest elevations of the watershed (1,500-2,500 m), all of the storms produce snowfall, and most (7 of 10) ARs produce a mixture of rainfall and snowfall over the duration of the storm. All ARs result in an increase in SWE in the headwaters, in some regions in excess of 300 mm during a single storm. Most of the ARs also result in RoS at intermediate elevations (between approximately 500-1,500 m), resulting in a clear elevation band of decreases in SWE. Of the 10 ARs, the only storms that do not result in a decrease in SWE are those that have negligible or nonexistent snowpacks prior to the onset of the storm (AR 1 and AR 10), and one smaller storm, which was also colder (AR 6).

Hydrologic Water Partitioning and Contributions to Groundwater
A basin-wide water budget partition of WY2017 precipitation is shown in Figure 5. The year-end fate of WY2017 precipitation is partitioned as either exited particles (evapotranspiration, ET, or discharge, Q), or into groundwater storage (both for individual storms and for water year totals; Figure 5). A time-series of basin-wide water partitioning is included in Figure S1 in Supporting Informatuon S1. Despite ARs contributing 52% of the annual precipitation, AR precipitation results in 15% less water going to annual groundwater storage than non-ARs (Figure 3b). In part, this difference can be attributed to the characteristically high rainfall rates of ARs and therefore greater potential for less infiltration and more runoff (either via infiltration excess overland flow or, if the storm duration is sufficient enough and/or the antecedent soil moisture conditions were wet enough, via saturation excess). Because ParFlow-CLM allows for both saturation excess and infiltration excess to be simulated, both runoff generation types are considered here. Interestingly, 93% of the AR contribution to groundwater is derived from rain, whereas groundwater contributions from non-ARs are 77% from rain and 23% from snow. Recall that annual snowfall amounts from ARs and non-AR storms are nearly equivalent (see Figure 3), so the differences in groundwater composition (rain vs. snow derived) cannot be attributed to snowfall-derived precipitation amounts. A possible explanation for the difference in groundwater contribution source is the high prevalence of RoS conditions during ARs (see Section 4.2). Figure 5b also shows the fraction of each water compartment that was derived from RoS events (stippled boxes). RoS largely governs discharge resulting from AR precipitation, where over half of discharge stems from AR-driven RoS events. In contrast, a negligible amount of non-AR discharge is driven by RoS events.
The fractional behavior of the 10 ARs ( Figure 5c) and annual totals (Figure 5d) help to differentiate and quantify contributions more explicitly when the volumes are normalized by precipitation amount. Fractionally, the behavior of AR and non-AR precipitation are very similar (Figure 5d), except for the contributions to groundwater from rain versus snow. Importantly, AR snow contributes significantly less to groundwater compared to non-AR snow: approximately 300% more non-AR snow resides in groundwater storage by year end compared to AR snow. These results suggest that the most significant difference in AR-sourced precipitation is on the groundwater recharge-snow relationship.

Dynamic Groundwater Storage
To better understand groundwater dynamics through WY2017 and between the vadose and saturated zone, Figure 6 shows the hourly relationship between watershed-averaged saturated and unsaturated storage, color-coded by month in (a) and by individual storm in (b). The cross-hairs of each diagram indicate the initial state of the storage at 1 October 2016, and the four quadrants indicate relative increases/decreases relative to that point in time.
By plotting the dynamic interaction of saturated and unsaturated storage, we are able to better understand overall groundwater storage in the presence (and absence) of ARs, as well as the system response through time. AR event-based ( Figure S2 in Supporting Informatuon S1) and monthly ( Figure S3 in Supporting Information S1) groundwater pressure head change maps are included in Supporting Information S1.
October groundwater storage (Figure 6a) shows the continued decline of the overall storage during the tail-end of summer baseflow conditions, but is quickly followed by sharp increases in winter and spring storage where the majority of storage variation occurs. During these wet months, relatively high groundwater storage increases from the initial state are evident. The timing of the 10 discrete AR events (Figure 6b; note separate color-legend) shows sharp and rapid changes in groundwater dynamics during two periods. These include periods when both the saturated and unsaturated storage are larger than the initial states (quadrant I), but also during long periods of higher saturated storage and lower unsaturated storage relative to the initial state (quadrant II). The latter occurs as near-surface variably saturated regions become completely saturated, effectively acting as a "trade-off" of storage as that residual storage becomes saturated. Lastly, as the final spring precipitation events occur and snow melts in April, the dynamic storage trend down toward quadrant III, where both saturated and unsaturated groundwater storage show decline relative to initial groundwater states. These simulations suggest that total groundwater storage will decrease by year-end relative to the beginning of the WY, despite WY2017 being the wettest on record for the Cosumnes watershed. This is consistent with Maina, Siirila-Woodburn, Newcomer, et al. (2020), which showed that groundwater pumping rates in the Sacramento Central Valley result in a net negative impact on watershed-average groundwater storage even during WY2017.

Conclusions
ARs play a central role in western United States precipitation, but to date there exists a lack of knowledge on how these storms impact the integrated hydrologic cycle, especially groundwater, which is crucial to long-term water resources and especially during periods of drought. With the use of coupled, atmosphere and integrated hydrologic models and above-, and below-ground water tracking techniques, we simulate a record wet year of an unimpaired northern California watershed. We used this model to draw lessons learned from how AR and non-AR water partitioning drives groundwater recharge. To the best of our knowledge, this study is the first to examine the role in which ARs contribute to groundwater dynamics. Our primary results include: 1. Despite ARs contributing over half of annual total precipitation, the contribution to end-of-year groundwater storage is less than the contribution made by non-AR precipitation. 2. Occurrences of RoS plays an important role in AR-driven discharge, where over half of AR discharge from snow comes from RoS events. 3. ARs result in 300% less snow-derived groundwater-recharge compared to non-ARs. 4. Despite WY2017 being the wettest year on record, groundwater pumping results in a net loss of groundwater storage by year-end.
Computational limitations associated with the full atmosphere-through-bedrock modeling and tracking algorithms, both which required supercomputing resources, limited the period of examination (a single water year; albeit an extreme year) that we could conduct in this study. Thus, our study represents a "proof of concept" for future work when computing resources enable us to conduct a longer set of simulations that better estimate the range of climate variability in California. WY2017 provides, at least, a partial analogue into the future given the higher occurrence and intensity of ARs that occurred. Antecedent soil moisture and groundwater conditions prior to each storm may impact the hydrologic response of each storm, especially the first storm of the season following a dry summer; here we focus on the hydrologic response of an observed extremely wet winter with a large sample size of ARs to understand their individually sequenced, and aggregate response on the hydrologic cycle to develop some key "lessons learned." Future studies could expand this study design to consider more scenarios of AR impacts on hydrology such as drought years with occasional extreme ARs (e.g., WY2015 or WY2018) or more near-average years (e.g., WY2009 or WY2016). Future studies could also consider the process-based hydrologic changes of ARs we have examined here but in conjunction with their societal and environmental impacts such as flooding (e.g., Konrad & Dettinger, 2017), levee breaks (e.g., Florsheim & Dettinger, 2015), and postfire debris flows (e.g., Oakley et al., 2017).
The results presented here suggest serious implications for groundwater recharge and water management given a future where annual precipitation is derived from larger, more intense, and warmer storms. Future studies may also examine AR conditions such as those from WY2017 but with synthetic modifications such as warmer air temperatures, different water demands and/or antecedent conditions, changes in precipitation phase (i.e., a low-to-no snow future; Siirila-Woodburn et al., 2021), and with changes in vegetation and/or land use (as projected by Earth system models in the coming decades) or different spring season characteristics (e.g., McEvoy & Hatchett, 2023). Isolating these factors of future climate conditions in the presence of ARs may help to inform watershed dynamics in the future as well as prioritizing investments intended to minimize flood risk and maximize water storage, including flood protection measures for unimpaired basins, water capture and recharge procedures such as managed aquifer recharge, and in basins with reservoirs, optimization of surface water reservoir and conveyance operations such as forecast informed reservoir operations.