Modeling managed aquifer recharge processes in a highly heterogeneous, semi-confined aquifer system

Widespread groundwater overdraft in alluvial aquifer systems like the Central Valley (CV) in California, USA, has increased interest in managed aquifer recharge (MAR). Like most clastic sedimentary basins, recharge to the productive semi-confined CV aquifer system remains a challenge due to the presence of nearly ubiquitous, multiple confining units (silt and clay) that limit recharge pathways. Previous studies suggest the presence of interconnected networks of coarse-texture sand and gravel deposits that bypass regional confining units over a small fraction of the CV near the American and Cosumnes rivers. Here, variably saturated infiltration and recharge processes were simulated across a domain that includes high-resolution representation of the heterogeneous alluvial geologic architecture in this area. Results show that recharge potential is highly dependent on subsurface geologic architecture, with a nearly 2 order-of-magnitude range of recharge across the domain. Where interconnected coarse-texture recharge pathways occur, results show that these features can (1) accommodate rapid, high-volume MAR and (2) propagate widespread and rapid pressure responses over multi-kilometer distances in the semi-confined aquifer system. For all MAR simulations, results show that the majority of MAR is accommodated by filling unsaturated-zone (UZ) pore volume. Results also show that coarse-texture UZ facies (where present) accommodate the majority of MAR volume during early time, but fine-texture facies ultimately accommodate the majority of the total MAR volume, even for coarse-dominated sites. These findings highlight the large variability of MAR potential across the landscape and demonstrate the importance of fine-texture facies for accommodating MAR in alluvial aquifer systems.


Your article is published under the Creative
Commons Attribution license which allows users to read, copy, distribute and make derivative works, as long as the author of the original work is cited. You may selfarchive this article on your own website, an institutional repository or funder's repository and make it publicly available immediately.

Introduction
Civilizations have typically obtained water from natural and constructed surface-water resources throughout most of human history. Only during the last 50-70 years has a significant quantity of water for humans been obtained through pumping from wells (Alley et al. 2002). During this short time, alarming levels of groundwater depletion have been observed in many regions, especially in semi-arid and arid areas that rely heavily on groundwater pumping from clastic sedimentary basins (Konikow 2013;Taylor et al. 2013;Wada et al. 2010Wada et al. , 2012. Groundwater has commonly been a source of high-quality freshwater and an important safeguard against uncertain inter-annual and inter-decadal shortfalls in precipitation and surface-water supplies (Hanson et al. 2012). However, overdraft of this important resource has only accelerated during the twenty-first century (Wada et al. 2011) and is further threatened by future climate uncertainty (Milly et al. 2008;Mirchi et al. 2013). Despite continued unsustainable groundwater abstraction in many areas, water policy efforts continue to respond to near-term crises and fail to anticipate long-term future conditions (Karl et al. 2009).
Large inter-annual variability of precipitation and streamflow (Dettinger et al. 2011) and heavy reliance on groundwater pumping for agricultural irrigation (Scanlon et al. 2012) has created an especially precarious scenario for water management in California's Central Valley (CV) aquifer system, which has had overdraft conditions in its southern portion for decades (Brush et al. 2013;Faunt et al. 2009;Konikow 2013;Scanlon et al. 2012). Historically abundant snowmelt runoff from the Sierra Nevada provides an estimated 54% of water for CV crops, on average (Faunt et al. 2009), with the remainder provided by direct precipitation and groundwater. Average temperatures in California (CA) are expected to increase by 1.5-4.5°C by the end of the twentyfirst century (Cayan et al. 2008), which will decrease the proportion of precipitation as snow and initiate earlier spring snowmelt and runoff (Cayan et al. 2008(Cayan et al. , 2010Vicuña and Dracup 2007), increase evapotranspiration (ET) and decrease late-season baseflow (Hayhoe et al. 2007;Huntington and Niswonger 2012), and likely increase the likelihood of cooccurring flooding and water-shortages in the same water year (Knowles et al. 2006;Swain et al. 2018). Future population growth and land-use change in CA, USA, and will likely increase drought risk (Barnett et al. 2008;Cayan et al. 2010), and elevate competition for existing water resources (Gleick 2000). These stressors are not unique to CA, but are symptomatic of increasing vulnerability of water resources worldwide (Stewart et al. 2005;Vicuña et al. 2011), especially in snowmelt-fed semi-arid and arid regions (Konikow 2013).
To reverse the negative effects of overexploitation of groundwater resources, groundwater must transition from being treated mainly as an extractive resource to one in which recharge and subsurface storage are pursued more aggressively. This remains a challenge because unlike surface-water reservoirs that are typically replenished on annual timescales, the clastic sedimentary aquifer systems are replenished on much longer time scales (Taylor et al. 2013), especially if no particular effort is devoted to augmenting groundwater recharge. Managed aquifer recharge (MAR) has been used for decades to supplement natural recharge and to strategically store surface water in groundwater aquifer systems for future water supply (Bouwer 2002), often as part of a conjunctive use (Bredehoeft and Young 1983) and/or water market framework (Israel and Lund 1995), including in select groundwater basins in CA (e.g., Asano 2016; Kletzing 1987), elsewhere in the Southwestern US (e.g., Jacobs and Holway 2004), and globally (e.g., Dillon et al. 2019). There is interest in expanding use of MAR in CA, both to offset overdraft and to hedge against future decreases in snowpack-water storage and changes in the timing and volume of surface-water availability. Established MAR projects commonly use dedicated infiltration basins located over locally coarse-texture geological deposits to increase recharge, but increasingly, MAR on agricultural fields during nongrowing seasons (Ag-MAR) has been proposed as an alternative to infiltration basins (Dahlke et al. 2018;Harter and Dahlke 2014;Niswonger et al. 2017). Studies have noted that even during periods of water scarcity, wet-season high-magnitude streamflows (HMF) can often provide ample unmanaged surface water for MAR during nongrowing seasons in CA (Beganskas and Fisher 2017;Kocis and Dahlke 2017) and elsewhere (Chinnasamy et al. 2018). The largest quantity of winter HMF in CA are in the CV, where most of the statewide groundwater overdraft occurs. HMF in the CV typically occur during episodic 5-7 day windows (Kocis and Dahlke 2017), and can be quite large (3.2 km 3 annual average during years with HMF). However, augmenting groundwater recharge with ephemeral HMF remains a challenge in the CV because the sedimentary aquifer system is composed of mostly silt and clay sediments (Faunt et al. 2009) that form nearly ubiquitous, multiple confining layers that create semiconfined conditions and limit infiltration rates over most of the landscape.
Geologic heterogeneity is ubiquitous across scales and strongly affects movement of water and solutes through the subsurface; however, seldom are enough subsurface data available to represent heterogeneous features explicitly in models (De Marsily et al. 2005;Koltermann and Gorelick 1996). As a result, development of groundwater models often involves simplifying or up-scaling heterogeneity to enable adequate representation of regional groundwater flows for purposes of regional water resources management. For example, in models of sedimentary aquifer systems (e.g., Phillips and Belitz 1991;Fogg 1986) where the aquifer sediments amount to 20-50% of the aquifer systems (the remainder being aquitard sediments), effects of the ubiquitous aquitard beds are approximated by regionally reducing the vertical hydraulic conductivity (K v ) relative to the horizontal hydraulic conductivity (K h ) by several orders of magnitude in order to match both horizontal and vertical hydraulic gradients. This approach is capable of producing good approximations of regional flows and groundwater budget components, but tends to smooth any local variations in recharge and vertical flow in the sedimentary connected network (e.g., Fogg et al. 2000;Fleckenstein et al. 2006). Stochastic methods like transitionprobability-based-indicator geostatistics provide opportunity for realistically representing subsurface heterogeneity of the major aquifer and aquitard facies while honoring measured data (Carle and Fogg 1996;. Results from studies implementing these methods show strong influence of subsurface heterogeneity on groundwater/surface-water interactions (Engdahl et al. 2010;Fleckenstein et al. 2006;Liu 2014).
A number of studies have measured infiltration and recharge processes in porous media in field settings (e.g., Batlle-Aguilar and Cook 2012; Bresciani et al. 2018) and laboratory settings (e.g., Fichtner et al. 2019). Other work has demonstrated the importance of geologic heterogeneity on infiltration processes and stream-aquifer interactions, including the often-significant contribution of focused stream leakage in arid and semi-arid areas (e.g., de Vries and Simmers 2002;Bresciani et al. 2018;Irvine et al. 2012), including the CV (Fleckenstein et al. 2006). Geologic heterogeneity is also important for natural recharge processes and MAR in karst systems (Hartmann et al. 2017;Xanke et al. 2017). In heterogeneous clastic sedimentary systems, even a small fraction of permeable hydrofacies in correlated random media tend to be interconnected in three dimensions (Fogg et al. 2000;Harter 2005), especially in the absence of spatially persistent geologic unconformities, providing potential recharge pathways to semi-confined aquifer systems. Interconnected, highly permeable sand and gravel deposits have been shown to occur in select locations in the southern CV that are potentially conducive to considerably higher rates of regional recharge than would be possible over the rest of the landscape (Weissmann et al. 2004). Studies have suggested the presence of these features in the northern CV (Meirovitz 2010;Shlemon 1967), compelling further study of MAR dynamics in this system. Several studies have focused on synthesizing MAR suitability characteristics from a combination of spatial data, including remote-sensed imagery, geologic maps, and soil surveys to identify favorable surface site characteristics for MAR in CA (O'Geen et al. 2015), and elsewhere (Adham et al. 2010;Ghayoumian et al. 2007). These studies provide a valuable initial survey of site suitability, but do not account for deeper subsurface geologic heterogeneity that has been shown to be important for recharge (e.g., Weissmann et al. 2004).
This research aims to explicitly simulate variably saturated water-flow dynamics in a highly resolved representation of complex subsurface geologic heterogeneity of the CV that includes both interconnected, highly permeable sand and gravel deposits and more typical silt-and clay-dominated sediments. Additionally, the goal of this research is to (1) gain insight into infiltration and recharge phenomena that are challenging to observe and have not been included in regional-scale groundwater models, and (2) guide MAR strategies for regions reliant on diminishing snowpack water storage and in overdrafted groundwater basins.

Domain extent and local hydrostratigraphy
The model domain ( Fig. 1) covers about 1,640 km 2 of lowelevation alluvial fan topography typical of the CV, CA. The climate of the area is Mediterranean, with 75% of annual precipitation occurring between November and March. The domain comprises the lower portions of the American and Cosumnes River watersheds on the east side of the northern CV, CA (i.e., the Sacramento Valley), where previous studies have shown the presence of massive, interconnected highly permeable sand and gravel deposits embedded in a matrix of fine-grained sediments (Meirovitz 2010;Shlemon 1967).
Such localized stratigraphic features have been shown to be potentially conducive to considerably higher rates of regional recharge than would be possible over the rest of the landscape (Weissmann et al. 2004), and likely occur throughout the CV but are mostly still undiscovered (e.g., Weissmann et al. 2005).
The hydrostratigraphy of the domain area includes the overlapping American River and Cosumnes River alluvial fans, each of which have markedly different depositional characteristics. The northern portion of the domain is dominated by the American River fan, which drains a large catchment (>4,000 km 2 ) that extends to the Sierra Nevada crest->3,000 m above mean sea level (amsl). Because of its high catchment elevation, deposition of the American River fan was significantly influenced by cyclic Plio-Pleistocene climate variation and episodic Sierra Nevada glaciation, which deposited relatively coarse-grained sediments in paleochannels (Shlemon 1967;Ford and Chee 1974). Some deep coarse-texture paleochannels, called incised valley fill (IVF; Weissmann et al. 2004Weissmann et al. , 2005, have been identified in the domain area (Meirovitz 2010). Incised valley fill are the result of a combination of interglacial river incision and subsequent glacially derived, dominantly coarse-grained sediment deposition in an aquifer system that otherwise contains mostly horizontally stratified fine-grained deposits. Plio-Pleistocene Sierra Nevada glacial cyclicity resulted in multiple IVF depositional events, which, in places, overlap and interconnect within the fan sequence from land surface into the deeper aquifer system to provide a relatively high permeability conduit for groundwater recharge that can bypass local confining units to depths >30 m (Meirovitz 2010). Incised valley fill associated with recent glaciation are especially important for recharge because they extend from near-surface, mainly just below the surface soil profile, downward into the subsurface. The most recent major glacial advance affecting the American River watershed, referred to as the Recess Peak, occurred approximately 13,000 years ago (Moore and Moring 2013). Importantly, all the major rivers draining highelevation catchments on the west side of the Sierra Nevada likely have IVF deposits, but few have been documented: e.g., the American River (Meirovitz 2010), Tuolumne River (Weissmann et al. 2005), and Kings River Weissmann et al. 2004). The southern part of the domain is intersected by the Cosumnes River, which is the last major undammed river in CA, and drains a smaller (900 km 2 ) lower-elevation catchment, which was much less influenced by glaciation, and hence lacks a shallow IVF deposits. As a result, the Cosumnes River fan is typically much finer in texture. Recent research suggests that some American River paleochannels have been shown to flow across and underneath the current path of the Cosumnes River in the Holocene and Quaternary due to a more southwest course of the ancestral American River when sea level was lower, creating complex, overlapping stratigraphy and multi-scale geologic heterogeneity (Meirovitz 2010).
Previous studies have identified the main groundwaterbearing units within the American-Cosumnes fan as the Laguna-Riverbank complex, and Mehrten formation (Ford and Chee 1974;Fleckenstein et al. 2006). The lithology of both units is an assemblage of granitic sand, silt, and clay, along with gravel bodies that are each up to 100 m thick in the domain area. Whereas the Laguna-Riverbank complex is Pleistocene/Pliocene in age and mostly granitic in origin, the underlying Mehrten formation is Miocene in age and primarily volcanic in origin. Studies have identified the presence of low-permeability paleosols in portions of the domain area, which can be difficult to differentiate from fine-grained silts and clays and can act as laterally persistent confining units, except when intersected by IVF (Niswonger and Fogg 2008).
The complex overlapping fan network and cross-cutting IVF in the domain area results in an aquifer system that, at shallow depths, is unconfined (and sometimes perched) or semi-confined and increasingly confined with depth (Fleckenstein et al. 2006;Liu 2014;Niswonger and Fogg 2008). Groundwater pumping in the region typically occurs in the deeper (>30 m depth) semi-confined/confined portion of the aquifer system (Liu 2014).

Hydrofacies model development
A transition probability Markov-chain geostatistical (TPROGS) model (Carle andFogg 1996, 1997;Carle 1999) was used to simulate the distribution of subsurface hydrofacies in the domain area. Development of the facies model with TPROGS is described in detail by Meirovitz (2010). Roughly 1,200 well logs, soil surveys, geologic cross-sections, and mapped paleochannels were used as conditioning data for the model. Subsurface sediments were parsed into four textural categories: gravel, sand, muddy sand, and mud (undifferentiated silt and/or clay)- Table 1 (Fleckenstein et al. 2004;Meirovitz 2010). The key parameters for developing the model include the proportion of each hydrofacies, the mean lengths of each hydrofacies along the principal directions, and the embedded transition probabilities to represent cross-correlation between different facies. Differences in depositional environment for the American and Cosumnes River fans resulted in important differences in these parameters for each region, requiring development of separate facies models for the northern and southern portions of the domain (i.e., for American and Cosumnes River fans, respectively). These models were subsequently combined to represent the sequence stratigraphy of overlapping fans and resultant multi-scale heterogeneity (Fig. 2). Meirovitz (2010) provides a detailed procedure for developing and combining the models. This process is summarized in section 'Hydrofacies model development with TPROGS'.
The model uses an orthogonal grid with the x-and ydirections rotated 17.85°counterclockwise from the cardinal directions, and with the z-direction oriented vertically. A grid cell size of 200 × 200 × 1 m was chosen to reflect the hydrofacies mean lengths in the x-, y-, and z-directions, respectively. The total domain size is 36.2 km × 45.4 km × 265 m in the x-, y-, and z-directions, respectively (i.e., 181 × 227 × 265 cells). Cells located above the land surface elevation are designated as inactive in the model, resulting in about 7.3 million active cells in the domain.
In these equations, S s is specific storage

Boundary conditions
Domain boundaries were chosen with the interest of simplifying the assignment of boundary conditions for the flow model. Model boundary conditions are discussed in greater detail in Liu (2014). The eastern boundary roughly coincides with the Sierra Nevada foothills, and the northern, southern, and western boundaries roughly coincide with the American River, Dry Creek, and the Sacramento River, respectively ( Fig. 1). A specified head boundary condition was applied for the eastern boundary, where heads were estimated from local monitoring well data (Liu 2014). A general head boundary condition was applied to the western boundary, which roughly coincides with the Sacramento River and Sacramento-San Joaquin Delta along the northwestern, and southwestern portions of the western boundary, respectively. A general head value of 0 m amsl was set 1 km beyond the western boundary to approximate these features. Because the groundwater flow direction is generally from east to west, no-flow boundary conditions were applied along the northern, southern, and  bottom boundaries. Combinations of specified-flux and specified-head upper boundary conditions were assigned both for model spin-up and for the recharge simulations, and are described in greater detail in subsequent sections.

Hydrofacies hydraulic properties
Hydraulic properties for each facies category (Table 2) were calibrated manually to simulate observed well hydrographs in both the shallow and deeper aquifer system; the calibration process is described in detail in Liu (2014) and is summarized in section 'Hydrologic model calibration'. To simulate the deeper aquifer system where geologic conditioning data are limited, a separate, uniform facies designation ('deep aquifer'; Table 2) with upscaled hydraulic properties was designated for the lower portion of the domain (Liu 2014). The anisotropy ratio (R) of vertical K s to horizontal K s was considered equal to unity (i.e., no anisotropy) for all facies except for the deep aquifer facies because it was assumed that the highresolution juxtaposition tendencies of facies, which is preserved by the geostatistical method, represents the true anisotropy of the system . A value of R = 0.003 was specified for the deep aquifer facies to reflect upscaled alluvial aquifer properties calculated by model simulation experiments of Liu (2014). Calibrated hydraulic properties are consistent with the range of literature values for the CV, CA, and for similar alluvial systems (Anderson et al. 2015;Botros et al. 2009;Fleckenstein et al. 2004;Frei et al. 2009;Maserjian 1993;Niswonger and Fogg 2008;Sager 2012).

Model spin-up
A 16-year simulation period reflecting water years 1970-1985 was used to bring the simulated hydrology into dynamic equilibrium and generate a realistic water-table configuration and vertical groundwater gradients within the domain is described in detail in Liu (2014). An initial potentiometric surface was specified using interpolated groundwater level data. Monthly estimated urban and agricultural groundwater pumping rates were applied as specified fluxes representing wells screened in lower portions of the domain that coincide with typical screened intervals of municipal and agricultural pumping wells in the region. Dominant sources of recharge for the region include stream recharge from the American River, Cosumnes River, and Deer Creek, and from deep percolation of agricultural and urban return flows, while recharge from precipitation is minimal. Estimated monthly urban and agricultural recharge volumes were applied as specified-flux boundary condition across the top of the domain based on groundwater model development by the Sacramento County Groundwater Agency (RMC 2011). Weekly estimates of spatially distributed river stage for the American River, Cosumnes River, and Deer Creek were applied as specified heads along land surface cells coincident with each of these features. Each of these water budget components were adjusted manually along with hydraulic properties as part of the model calibration process and are described in detail by (Liu 2014) and summarized in section 'Hydrologic model calibration'. An additional 1-year spin-up period was simulated in which a uniform 1 mm day −1 specified-flux boundary condition was applied across the domain to equilibrate soil moisture conditions in the near-surface UZ cells. This additional period was necessary to facilitate model convergence during the initial time steps of subsequent recharge simulations. The spin-up period was shown to be sufficient to equilibrate soil-moisture storage according to qualitative metrics described by Ajami et al. (2014).

Recharge simulation experimental design
Site selection Five 5.76 km 2 recharge sites ( Fig. 2), each encompassing 144 upper-boundary cells, were chosen to approximate hypothetical MAR infiltration basins. The size of each site was chosen to reflect a regional-scale MAR site. In CA alone, MAR sites range from individual infiltration basins over several hectares  Beganskas and Fisher 2017) to large networks of basins >25 km 2 in size (e.g., Kern Water Bank Authority 2018). Site locations were chosen ad hoc to maximize a range of site characteristics favorable for MAR, including the initial depth-to-water, proportion of coarse-texture hydrofacies (i.e., sand and gravel) at the land surface and in the unsaturatedzone below each site, and effective vertical K s of UZ coarsetexture hydrofacies, as approximated by the geometric mean of K s for vertically coincident UZ cells (Fogg et al. 2000; Table 3; Fig. 3). All sites were chosen to minimize interaction with lateral boundary conditions. Sites 1-5 range from characteristically most favorable (site 1) to least favorable (site 5) for recharge. Vertical cutaways of hydrofacies distributions are shown in Fig. 4, which illustrate the relative proportions, and local degree of interconnection of coarse-texture facies beneath each recharge site.

Recharge scenarios
Two sets of numerical experiments (i.e., recharge scenarios) were simulated for each of the five sites: (1) a 180day recharge scenario was simulated for each site to evaluate the system response to recharge stress over a prolonged recharge period, and (2) a 3-year simulation was run for each site, wherein recharge was simulated for the initial 90 days of each year, followed by 275 days without recharge to evaluate the effects of successive, ephemeral recharge stresses on system response. The 90day recharge periods chosen for the 3-year simulations correspond with the January-March peak wet season for the Cosumnes River catchment (Booth et al. 2006). For all recharge simulations, surface ponding was approximated by a specified head boundary condition representing 10 cm of ponding at each recharge site, with no recharge specified for the remaining upper-boundary cells. Finally, two additional 180-day and 3-year simulations were run in which no recharge was specified for all upper-boundary cells, i.e., as no-recharge scenarios. Results were output at 2-day and 5-day intervals for the 180-day and 3-year simulations, respectively.

Model superposition and post-processing
The spatiotemporal distribution of ψ, and total subsurface water storage (TSS) from recharge stress was isolated from other model stimuli with a model superposition approach. At each time step, the spatial distribution of ψ, and TSS for the norecharge simulation was subtracted from the distribution of ψ and TSS for each recharge simulation. In this way, perturbations in ψ and TSS from the recharge stress are isolated from other model stimuli at each model cell, including transient model response to regional boundary condition effects. For each recharge simulation, the pressure response of the recharge pulse was quantified by determining the number of model cells with ψ perturbations above 1-cm, 10-cm, and 1m thresholds. To isolate system response in both saturatedand unsaturated-zones, ψ and TSS perturbation responses were further delineated according to whether the response in a given model cell occurred above or below the initial water table (i.e., S w < 1 or S w = 1 for each cell, respectively).

Domain-wide pressure and recharge response
Temporal variation of recharge rates and differences in mean recharge rates and volumes Results from the five 180-day recharge simulations show large temporal variation of recharge rates at each site, typically with a rapid decay from initially high rates (Fig. 5a). Sites 1 and 2 showed very high early-time recharge rates (130 and 51 cm day −1 , on average, during the first 10 days, respectively) relative to both late-time rates at these sites and to rates at other sites. This behavior reflects the high proportion of near-surface, coarse-texture facies in the unsaturated zone beneath those sites (Table 3), which can rapidly accommodate large recharge volumes. For all sites, the average recharge rate during the first 30 days was >4× greater than for last 30 days (Table 4). Recharge rate also varied widely between sites, with a 62× greater 6-month average recharge rate for site 1 (18.5 cm day −1 ) than for site 5 (0.3 cm day −1 ). Similarly, the cumulative change in subsurface storage (i.e., cumulative recharge volume) at the end of the 180-day simulation was 63× greater for site 1 than for site 5 (Fig. 5b). For all sites, >50% of the total recharge volume for the 180-day simulation is achieved during the first 50 days of simulation time. Pressure response at 1-cm, 10-cm, and 1-m thresholds The pressure response of the recharge pulse was quantified by determining the number of cells with ψ perturbations above 1cm, 10-cm, and 1-m thresholds, at each timestep for all five sites (Fig. 6). Results show that at the end of the 180-day simulation, site 1 has 2.4×, 5.1×, 9.5× more cells perturbed at 1-cm, 10-cm, and 1-m thresholds, respectively, than for sites 2-5 combined. Differences between site 1 and the remaining sites were less for lower (1-cm) thresholds and greater for higher (1-m) thresholds. For site 1, >907,000 cells were perturbed at a 1-cm threshold at the end of the 180-day simulation, which is~12% of the active domain volume. For site 5, only 496 cells were perturbed at a 1-cm threshold (i.e., >3 orders-of-magnitude fewer than for site 1) because propagation of the recharge pulse did not extend beyond the uppermost several meters of unsaturated facies during the 180-day simulation window. This response reflects the fact that the uppermost surface layer of site 5 is comprised entirely of mud facies, which impede recharge.
Pressure propagation distances Figure 7a shows the maximum distance of ψ perturbation at a 1-cm threshold from the center of each recharge site for each timestep. Results show the maximum distance of ψ perturbation propagation is >4 km at the end of the 180-day simulation for sites 1-4, and > 10 km for site 1. For all sites, half of the maximum distance at 1-cm threshold was achieved during the first 40 days. Figure 7b shows the spatial distribution of ψ perturbation propagation for each site at 1-cm, 10-cm, and 1m thresholds in plan view.

Influence of geologic heterogeneity on domain-wide response
Visualizations of the ψ perturbation response of each recharge pulse throughout the domain after 180 days are shown in Fig. 8. The relatively large area of influence of ψ perturbation response for site 1 (Fig. 8a) reflects the large proportion of UZ coarse texture facies (0.88) and the relatively greater average effective vertical K s for UZ facies (42.2 m day −1 ) compared to other sites (Table 3). The complex, non-uniform ψ perturbation response at each site reflects the influence of subsurface hydrofacies heterogeneity on recharge pulse propagation. For example, Fig. 8f highlights the influence of a 'choke point' at site 2, wherein the recharge pathway is reduced to a narrow network of interconnected high-K facies at~10-m depth, despite the upper-most surface layers being comprised entirely of gravel. Minimal propagation of ψ perturbation observed for site 5 (Fig. 8g) reflects the absence of any coarse-texture facies at land surface (Table 3), which impedes infiltration, and confines the response to the uppermost several meters of nearsurface facies.

Unsaturated-vs. saturated-zone behavior
Pressure response and change-in-storage above and below the initial water table Further analysis of the recharge response included parsing the ψ and TSS perturbations according to whether the response occurred in the unsaturated zone or in saturated zone. This designation was delineated by the initial water-table configuration at t = 0 of the simulation. For example, a response in any cells above the initial water table was designated as an unsaturated-zone response, even though cells directly below each recharge site likely become saturated as the wetting front advances from land surface to the water table. Results show that the majority of the total change in storage occurs by filling previously unsaturated cells, but the majority of the ψ perturbation response is propagated through the saturated system. Figure 9a shows the arrival of the wetting front at the water table as a rapid (but still small) increase in the proportion of the change-in-storage response below the initial water table. Wetting front arrival at the water table occurs at t = < 2, 10, 28, and 16 days for sites 1-4, respectively. For site 5, the wetting front did not reach the water table during the 180-day simulation, and was not included in this analysis. At the end of each simulation, the proportion of change-in-storage that occurred below the initial water table ranged from <0.01 to 0.06 for sites 2-4 and was 0.23 for site 1. For all sites, the majority of change-instorage volume accommodated by filling unsaturated pore space which is primarily controlled by ϕ and is substantially greater than S s (Table 2), which controls storage response in the saturated zone. Figure 10b also shows a rapid increase in Fig. 4 Vertical cutaways of the subsurface hydrofacies configuration below each recharge site, where a includes all hydrofacies, b includes coarse-texture hydrofacies (i.e., gravel and sand), and c includes finetexture hydrofacies (i.e., muddy sand and mud) and deep aquifer Legend the proportion ψ perturbation that occurred below the initial water table upon the arrival of the wetting front at the water table. Unlike change-in-storage, however, the proportion of the ψ perturbation response occurring below the initial water table was >0.90 for sites 1-4 at the end of each 180-day simulation. Results show that while the majority of recharge volume is accommodated by filling partially-saturated pore volume, the majority of the ψ perturbation response occurs by propagating through the semi-confined (and saturated) aquifer system.

By hydrofacies
Recharge response for each 180-day simulation was parsed according to whether change-in-storage was accommodated by coarse-texture (i.e., gravel, sand) or fine-texture (i.e., muddy sand, mud) facies. Figure 10 shows the proportion of the change-in-storage accommodated by fine-texture facies for each simulation, both cumulatively over the 180-day simulation ( Fig. 10a) and during each 2-day time step (Fig. 10b). Results show that the cumulative change in storage accommodated by fine-texture facies at the end of the 180-day simulation was >0.35 for all sites, and >0.50 for sites 3-5. Interestingly, the change-in-storage accommodated by finetexture facies at each time step for all sites was >0.50 after t = 90 days. Results show that fine-texture facies accommodate a substantial portion of the recharge volume, especially during late-time, even for coarse-dominated sites (i.e., sites 1-2). These findings demonstrate that while fine-texture facies may not facilitate widespread or rapid transmission of water volume or pressure response, these facies do in fact contribute substantially to the overall aquifer storage due to the relatively greater ϕ and S s for fine-texture facies than for coarse-texture facies ( Table 2). Results suggest that these results are minimally influenced by parameter choice (see section 'Parameter sensitivity on storage accommodated by fine-texture facies'). These findings are consistent with other research that has shown that simplified aquifer-system conceptualizations which treat fine-texture facies or confining units as noncontributing to overall aquifer storage are often inadequate (Konikow and Neuzil 2007).

Domain-wide recharge rates and volumes
Recharge rate and cumulative change-in-storage during each 90-day recharge period in the 3-year simulation were compared ( Fig. 11a,b, respectively). Results show nearly 2× greater average recharge rate during the 1st 90-day recharge period (year 1) as compared to the 3rd 90-day recharge period (year 3; Table 5). Results show the cumulative recharge   Legend volume for site 1 is 3.4× greater than for the next largest (site 2) and 74× larger than for the smallest (site 5). The 90-day cumulative recharge volume for each 90-day recharge period similarly showed an average 2× greater volume during year 1 than during year 3 (Fig. 11a). This response reflects both groundwater mounding (i.e., local increases in groundwater levels) and increases in S w in UZ cells below each recharge site after successive recharge periods, which reduce the unsaturated-zone pore volume that can be filled during subsequent recharge periods. Additionally, results were parsed according to whether the change-in-storage occurred in coarse-or fine-texture facies (Fig. 11c). For sites 1-2, results show most recharge volume is accommodated by coarse-texture facies during the initial 90-day recharge period, after which the volume decreases as infiltrated water is then slowly accommodated by fine-texture facies during the recovery period. This behavior is repeated for subsequent recharge periods, but the relative contribution of fine texture facies increases steadily with time. For all sites, the contribution of fine-texture facies to the cumulative recharge volume exceeds that for coarse-texture facies at the end of each 3-year simulation.

Proportion split between coarse and fine hydrofacies
The proportion of recharge volume accommodated by finetexture facies for each site during the 3-year simulation is shown in Fig. 12. Results for sites 1-3 show similar behavior, wherein the proportion of recharge volume in fine-texture facies generally increases with time, except during the second and third 90-day recharge periods (i.e., during years 2 and 3, respectively), during which time the proportion decreases due to rapid recharge into coarse-texture facies. In general, the proportion of recharge volume accommodated by fine-texture facies increases more rapidly during nonrecharge periods. For sites 1-3, each recharge inundation is initially accommodated by filling partially saturated, coarse-texture facies, which creates a large vertical pressure stress that is then transmitted rapidly through interconnected networks of coarse-texture facies into the saturated zone. Subsequently, the pressurization of the interconnected network of coarse-texture facies then 'bleeds into' adjacent fine-texture facies. For sites 4-5, which are dominated by locally high proportions of finetexture facies, the proportion of recharge volume accommodated by fine-texture facies generally decreases with time, reflecting spreading and interception of recharge by adjacent coarse facies. For all sites, the proportion accommodated by fines exceeds 0.50 between years 2 and 3. These results suggest that fine-texture facies are the largest subsurface reservoir for this system, even for recharge areas with coarse-dominated interconnected alluvial aquifers. However, much more time is required to add or remove water from this reservoir than for coarse-texture aquifer systems, which can be readily pumped and/or recharged over much shorter timescales.
Pressure propagation distances Figure 13a shows that the maximum distance of ψ perturbation at a 1-cm threshold from the center of each recharge site for each time step. Results show that the maximum distance of ψ perturbation propagation is >9 km at the end of the 3-year simulation for sites 1-4, and >16 km for site 1. For all sites, half of the maximum distance was achieved during the first year. For site 1, this distance was achieved within the first 90 days. Figure 13b shows the  Legend maximum lateral extent of ψ perturbation propagation at a 1-cm threshold at the end of years 1, 2, and 3.

Discussion
Increasing water scarcity has accelerated overdraft of groundwater resources in CA (Scanlon et al. 2012) and globally (Wada et al. 2011) and has created newfound interest in replenishing overdrafted aquifer systems with MAR. Studies have shown that thoughtful implementation of MAR practices can increase sustainability of groundwater resources where each 90-day recharge period is outlined with a gray background Fig. 12 Cumulative proportion of change-in-storage that is accommodated by fine hydrofacies (i.e., muddy sand and mud) during the 3-year simulation period for each recharge simulation. Each 90-day recharge period is outlined with a gray background  year 1 year 1 year 2 year 3 (Niswonger et al. 2017). However, widespread implementation of MAR is often impeded by institutional barriers like transfer of water rights and water accounting uncertainty (Asano 2016), along with infrastructure limitations, including water conveyance and land acquisition costs (Gailey 2018), and water quality concerns (Hartog and Stuyfzand 2017). Past research has shown the importance of heterogeneity on recharge processes in clastic sedimentary aquifer systems (e.g., Fleckenstein et al. 2006;Irvine et al. 2012). Results presented here further demonstrate that subsurface sedimentary geologic architecture is an important consideration for infiltration and recharge processes, and especially so when considering MAR effectiveness. Results show a highly transient recharge response that is consistent with field experiments (Batlle-Aguilar and Cook 2012) and a wide range of recharge rates that are possible across the landscape, including rapid, highvolume recharge in select areas where sand and gravel IVF outcrop at land surface. Siting MAR over IVF deposits has the potential for outsized recharge rates compared to the rest of the landscape. For example, these results suggest that the hypothetical cumulative recharge volume for site 1 during the 180-day scenario (~2.0 × 10 8 m 3 ) is >1/10th of the estimated annual groundwater overdraft for the CV, CA, during 2003-2010 (Famiglietti et al. 2011). Interannual variability of precipitation in CA is large, and is typically derived from just a few storm events (Dettinger et al. 2011), and excess surface water available for recharge in CA is flashy and typically falls within a short (<10 day) window (Kocis and Dahlke 2017), so identifying sites that can accommodate large volumes of recharge during a short period of time is paramount. Conversely, results show that laterally continuous fine-grained facies can impede MAR rates considerably, and may limit MAR feasibility over much of the landscape. Results show that vertical interconnection of coarse-texture facies is important for MAR, rather than just the upper-most facies designation. This suggests that GISderived surface metrics of recharge potential, while valuable, are not fully diagnostic of MAR potential. Rather, these data products should be considered as an important component of a thorough site evaluation that includes investigation of deeper subsurface geologic architecture.
Importantly, IVF deposits have been identified in several CA river fans, including the American, Tuolumne, and Kings rivers, and likely occur in other major river fans that drain high-elevation, glacially influenced catchments on the west side of the Sierra Nevada (Meirovitz 2010;Weissman et al. 2005;Weissmann et al. 2004). In addition, similar coarse-texture, glacially influenced IVF have been identified in other fans in the Western USA (e.g., Pierce and Scott 1983). The authors believe that results presented here can foster renewed effort on the part of the hydrologic sciences community to identify and catalog locations with IVF for MAR. These results demonstrate that the recharge potential for these features is sufficiently strong that they could be considered for special land use prioritization such as a recharge preserve.
One may conclude that the relationship between proportion of coarse-texture facies and recharge potential illustrated here can be deduced from first principles, and is thus trivial. Certainly, in a general sense, this finding is obvious; however, to the authors' knowledge, no study has used a 3D variably Fig. 13 a Maximum distance of pressure head (ψ) perturbation (1-cm threshold) from the center of recharge sites 1-4 for each 3-year simulation. Each 90-day recharge period is outlined with a gray background. Site 5 not shown due to minimal ψ propagation. b Pressure head (ψ) perturbation snapshots (plan view) for sites 1-5 for each 3-year simulation. Orange, red, and dark-blue outlines are the maximum lateral extent of ψ perturbation at 1-cm thresholds, at the end of years 1, 2, and 3, respectively. Square black boxes represent the lateral extent of the recharge site for each simulation Legend saturated water flow code to explicitly simulate MAR dynamics through a highly heterogeneous domain. This approach couples a detailed representation of subsurface geology with physically realistic water flow physics to elucidate important processes that can (1) help improve representation of recharge processes in coarse-resolution, management-focused groundwater models, (2) help prioritize site investigation and data collection for proposed MAR projects, and (3) inform management entities to the potential impacts of MAR.
Results illustrate an important dichotomy between change in storage and pressure response in the aquifer system. Indeed, results show that a pressure response can be registered in wells screened in the semi-confined aquifer system several kilometers from the origin of the recharge stress. Of course, the increase in pressure is not related to physical water from the recharge site entering that well, and these results illustrate this important concept, which can be a challenging to convey to laypersons, and has important implications for water rights and water management. For example, recently-passed groundwater management legislation in CA requires the creation of local groundwater sustainability agencies (GSAs) to limit both the "chronic lowering of groundwater levels" and "significant and unreasonable reductions in groundwater storage" (Kiparsky et al. 2016). Results presented here demonstrate that MAR can help mitigate both of these impacts. While the benefit of physical change in storage occurs locally, the increase in groundwater heads can be regionally beneficial, potentially benefitting adjacent jurisdictions outside of the immediate GSA.
These results also suggest that while networks of interconnected coarse-texture facies provide a conduit for rapid infiltration and widespread pressure response, the fine-texture facies accommodate a substantial fraction of the total recharge volume. This finding is consistent with other work showing the importance of fine-texture facies storage in aquifer systems (e.g., Konikow and Neuzil 2007), and challenges common aquifer-system conceptual frameworks, wherein finetexture facies are often considered a non or minimally contributing component of the aquifer system. These findings suggest that fine-texture facies may in fact be the largest reservoir in this alluvial aquifer system. Preliminary sensitivity analyses indicate this response is fairly robust to parameter uncertainty (see section 'Parameter sensitivity on storage accommodated by fine-texture facies'). Importantly, these findings support conceptual models of groundwater flow and storage in alluvial aquifer systems that include fine-texture facies as an important storage reservoir (e.g., Konikow and Neuzil 2007). In essence, the connected network of coarse-texture facies provide for relatively fast flow and recharge phenomena, while the finetexture facies end up accommodating most of the changes in storage but on longer time scales. The storage depletion and replenishment can be viewed as a two-stage process, in which rapid declines in storage occur in the coarse-texture aquifer network followed by slow drainage (leakage) from the fines. Conversely, during storage augmentation, the immediate increases occur in the aquifer network on time scales of days to months, followed by much slower but pervasive increases in storage in the fine-texture facies on time scales of months to years. From a whole-watershed perspective, one can deduce fine-texture facies to be the largest (but least accessible) reservoir within this system, followed by coarse-texture facies, and finally surface-water reservoirs, which are the most readily accessible and replenishable. In general, these results are somewhat reminiscent of the leaky aquifer analytical model development by Neuman and Witherspoon (1972), who pointed out an investigation bias toward the hydrology of aquifers and suggested that future work should focus on the aquiferaquitard complexes that compose aquifer systems.
The authors acknowledge some limitations to the approach-for example, the TPROGS technique for developing the geologic domain is informed by ample conditioning data; however, the approach is inherently stochastic, which can limit the robustness of facies prediction in areas of the domain with sparse conditioning data. In addition, only a single TPROGS realization was used for these simulations, and the authors acknowledge that a more rigorous ensemble approach could provide greater insight into potential facies distributions within the domain. In addition, the model spin-up included a domain-wide 1 mm day −1 recharge flux during the final year to facilitate model convergence for subsequent MAR simulations. This boundary condition may be unrealistic with respect to recharge rates reported for this semi-arid area (25-275 mm year −1 ; Fleckenstein et al. 2006) and may contribute to some overestimation of antecedent soil-moisture conditions in the uppermost model cells and influence recharge rates for subsequent MAR simulations. Importantly, this study should not be treated as a thorough site investigation for MAR in this region. Rather, the authors present these findings as a proof-of-concept to demonstrate the influence of geologic heterogeneity on MAR dynamics in a hypothetical but physically realistic domain. Moreover, the simulations presented here do not consider several surface processes that influence real-world MAR feasibility and dynamics, including topographic site limitations, evaporative losses, and clogging effects (Bouwer 2002). In addition, the authors acknowledge that a more detailed investigation of the role of geologic heterogeneity on K upscaling and a rigorous uncertainty or sensitivity analyses of hydraulic properties would permit broader interpretation of these findings, and thus warrants further study. Despite these limitations, these findings have implications for understanding MAR dynamics and for assessing MAR feasibility in clastic alluvial groundwater basins in CA and globally.

Conclusions
This research explores variably-saturated water flow dynamics in a highly heterogeneous geologic domain that reflects the complex alluvial geologic architecture on the east side of the Northern Central Valley, CA. The research objectives are to inform MAR implementation in CA and elsewhere by (1) highlighting the role of subsurface geology for recharge dynamics and (2) identifying important recharge phenomena that are not easily observed or simulated by typically coarse-resolution regional groundwater models. The approach uses the variably saturated water-flow code, ParFlow, to simulate recharge over a range of configurations of unconsolidated alluvial geology, including over sand and gravel IVF deposits that interconnect from land surface to the deeper semi-confined aquifer system. Two sets of recharge scenarios were simulated at five sites to evaluate system response to both prolonged, and shorter successive recharge stress.
Results show a large (nearly 2 order-of-magnitude) range of cumulative recharge volumes between sites that is dependent primarily on the configuration of subsurface geologic facies. Recharge rates were highly variable in time, with all sites showing relatively high initial rates (e.g., >100 and >50 cm day −1 during the first 10 days of simulation time for sites 1 and 2, respectively) followed by rapid decay to a quasi-constant recharge rate. Results demonstrate that the overall subsurface geologic architecture, rather than just the upper-most soils or facies designation, is important for recharge. All sites showed >4× reduction in average recharge rates for the last 30 days compared to the first 30 days of the 180-day simulation, and a~2× reduction in cumulative recharge for year 3 compared to year 1 for the multi-year simulations. This behavior reflects the effect of groundwater mounding that limits rapid filling of coarse-texture UZ facies.
Results suggest that the majority of recharge volume is accommodated by filling unsaturated-zone facies, but where there is sufficient hydraulic communication between land surface and the deeper aquifer system, the majority of the pressure response is propagated through the saturated aquifer system once the recharge wetting front intersects the water table. Results show that if there is sufficient hydraulic connection between the recharge site and the semi-confined aquifer system, the recharge pressure response can be widespread and rapid, propagating over several kilometers over a period of days or weeks. These results provide a valuable illustration of two physically distinct benefits of recharge: (1) local increases in groundwater storage and (2) the possibility of a more widespread re-pressurization effect in the regional aquifer system. The distinction between these responses has important implications for water rights, groundwater management regulations, and other water-policy issues.
Results also suggest that while the majority of water volume and pressure response is transmitted through coarsetexture facies, the majority of the recharge volume is eventually stored in fine-texture facies, even for sites that have disproportionally large fractions of coarse-texture facies. This result suggests that fine-texture facies are the largest, albeit least accessible, reservoir for recharge in this system. This finding has important implications for aquifer conceptualization, because fine-texture facies are often considered as aquitards (or aquicludes) that do not appreciably participate as part of the overall aquifer system. cross sections) to calculate lateral mean lengths and transition probabilities. Ten statistically similar domain realizations were generated by Meirovitz (2010). Only small variations in volumetric proportions and regional connectivity of facies were noted between realizations. For these reasons, and because of the significant time and computational resources required for model spin-up, as described in section 'Model spin-up', a single realization was chosen for the recharge simulations described herein (Fig. 2).

Hydrologic model calibration
The calibration effort by Liu (2014) was aimed at constraining model uncertainty and avoiding equifinality of parameter sets. In summary, estimates of water budget terms and aquifer properties were first gathered from previous studies (e.g., RMC 2011; Sager 2012) and used as preliminary model parameter values. Select parameters were adjusted manually within expected ranges to match transient groundwater-level observations within the domain area. Because coarse-texture aquifer properties were considered relatively well constrained by aquifer test data (RMC 2011), particular effort was made to estimate fine-texture facies K s by adjusting these values to match transient groundwater-level observations in the semiconfined and confined portions of the aquifer system where large, seasonal fluctuations in heads are very sensitive to K s of the fine-texture aquitard facies. Aquifer system configuration was considered adequately constrained by the ample well log data used to develop the geologic model (>1,000 well logs; Meirovitz 2010), wherein strong lateral connectivity of high-K facies was observed virtually throughout the domain area. A number of previous studies conducted within the area and in similar alluvial aquifer systems helped in the estimation of regional water budget terms as well as the expected ranges of aquifer properties for each facies type (Anderson et al. 2015;Botros et al. 2009;Fleckenstein et al. 2004;Frei et al. 2009;Maserjian 1993;Niswonger and Fogg 2008;Sager 2012). Specifically, regional agricultural pumping and recharge volumes were adjusted moderately by Liu (2014) to approximate average groundwater levels and changes in storage during a 16-year calibration period. Then, K s values of the finetexture facies (Table 2) were adjusted to approximate the appropriate seasonal fluctuations in heads which ranged systematically from small fluctuations at the water table to much larger fluctuations in deeper intervals, owing to greater confinement. Because the geologic model is a single stochastic realization, some local-scale deviation between simulated and observed heads was expected. Nevertheless, the strong lateral connectivity of the high-K facies in this model would produce very similar system behavior among different realizations (e.g., LaBolle and Fogg 2001;Fogg et al. 2000), obviating the need for multiple realizations in a study of this type.
Thirty-nine monitoring well locations distributed roughly evenly throughout the domain, including both 'shallow' (<20 m) semi-confined and 'deep' (>30 m) confined aquifer levels, were chosen to assess correspondence with measured water levels. Goodness of fit was evaluated both with root mean square error (RMSE) of water level residuals and with Pearson's r of simulated and observed water levels. Average RMSE and Pearson's r were 1.6 m and 0.97, respectively, for the calibrated model. During the calibration process, monthly regional agricultural pumping rates were adjusted within ±15% of initial literature values. K s values were adjusted within ±65%, on average, from initial values. Adjustments of K s were small relative the five order-of-magnitude range of K s for facies designations within the domain.
Parameter sensitivity on storage accommodated by fine-texture facies Seven additional numerical experiments were run to evaluate the influence of parameter choice on the proportion of storage accommodated by fine-texture facies for site 3. Six runs in which parameters describing hydraulic properties of fine-texture muddy sand and mud facies were perturbed and compared with a baseline run with the original model parameters (see Table 6). Each numerical experiment included a 90-day simulation with an initial 30day recharge period in which 10 cm of ponded water was simulated over the site, followed by a 60-day recovery period in which no recharge was specified over the entire domain. An additional 90-day no-recharge simulation was performed with identical parameters, and the outputs from the recharge and no-recharge runs were differenced following procedures described in section 'Model superposition and post-processing'. The proportion of storage accommodated by fine-texture facies was compared against the baseline run at the end of each 90-day simulation. Shorter-duration 90-day simulations were chosen to represent the recharge and fines-storage phenomena while minimizing computational expense. Similarly, parameter sensitivities were evaluated at a single site; site 3 was chosen because its recharge potential was balanced between the other high-and low-recharge sites (sites 1-2 and sites 4-5, respectively). Parameter perturbations were 10% of a parameter range proposed by Sager (2012). Results for the baseline run show that 60% of the recharge volume was accommodated by fine-textured facies at the end of the 90-day simulation, and ranged between 59 and 64% for each of the six runs in which parameters were perturbed (Table 6). Results indicate that the high proportion of storage accommodated by fine-textured facies is minimally influenced by parameter choice and is fairly robust to parameter uncertainty.
Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.