Airborne geophysical method images fast paths for managed recharge of California’s groundwater

Given the substantial groundwater level declines in the Central Valley of California, there is an urgent need to supplement the recharge of the groundwater systems by implementing managed aquifer recharge. With approximately 170 km3 (140 million acre-feet) of available groundwater storage space, water deemed to be excess during wet years could be spread on the ground surface at selected locations allowing it to move downward to recharge the underlying aquifer system. Along the eastern edge of the Central Valley there are large paleovalleys that can act as fast paths expediting the downward movement of water. These paleovalleys, incised and then filled with coarse-grained materials—sand, gravel, cobbles—at the end of the last glacial period, are referred to as incised valley fill (IVF) deposits. An IVF deposit has been mapped at one location in the Kings River alluvial fan, with others proposed to exist in the fans of major rivers. If located, these deposits would be optimal sites for managed recharge. In this study, we assessed the use of a helicopter-deployed geophysical method to efficiently locate IVF deposits throughout the Central Valley. We acquired 542 line-kilometers of airborne electromagnetic (AEM) data in the Kings River alluvial fan, with dense line-spacing over the Kings River IVF deposit which had been mapped as ∼2 km wide, extending over 20 km into the Central Valley, from the ground surface to a depth of 30 m. The IVF deposit was unambiguously imaged in the AEM data as an extensive linear feature that was more electrically resistive than the surrounding materials due to the high percentage of coarse-grained sediments. This study provides the evidence to support the rapid adoption of the AEM method to locate IVF deposits along the eastern edge of the Central Valley. These deposits provide valuable natural infrastructure for recharging California’s groundwater.


Introduction
California's Central Valley is one of the most productive agricultural regions in the world, providing 8% of the agriculture output of the United States (USGS 2022). The historical lack of groundwater management together with extended periods of drought and unrelenting demand for irrigation water have caused considerable overdraft in large portions of the valley. The result has been a significant decline in groundwater levels, now more than 150 m below the ground surface in parts of the valley (California Department of Water Resources (CDWR) 2022b). This decline is raising serious concerns about the associated impacts-non-sustainable depletion and eventual exhaustion of the resource, drying up of the relatively shallow wells that serve as the domestic drinking water supply for vast numbers of homeowners and disadvantaged communities (Pauloo et al 2020), subsidence of the ground surface (Zektser et al 2004, Faunt et al 2016, Jeanne et al 2019, Lees et al 2022, and degradation of water quality , Levy et al 2021, Pauloo et al 2021. Climate change predictions for California include an increased intensity of drought during dry years but also include an increased intensity of storms during wet years. In order to implement sustainable groundwater management, which is now a legislative requirement in the state (CDWR 2014), many water agencies are turning to the concept of managed aquifer recharge where excess surface water is captured during the wet years and used to recharge the groundwater systems in an effort to recover from, and prepare for, the droughts. Owing to the Central Valley groundwater level declines, there is approximately 170 km 3 (140 million acre-feet) of subsurface storage space available for storing additional recharged groundwater (The Nature Conservancy 2016). For comparison, the total surface water storage capacity in California is approximately 52 km 3 (42 million acre-feet).
One approach to managed recharge is to spread the water on the ground surface, for example in designated recharge basins, so as to recharge the underlying groundwater. But the heterogeneous groundwater systems of the Central Valley consist predominantly of relatively low-permeability silt and clay beds that limit the amount of water than can move from the surface to deeper parts of the systems (Weissmann and Fogg 1999, Fogg et al 2000, Maples et al 2019. Along the eastern edge of the Central Valley, however, there are geologic features that offer prime connections to move water from the ground surface-large paleochannels, referred to as incised valley fill (IVF) deposits (Weissmann and Fogg 1999). These deposits were originally formed when valleys were incised at the end of the last glacial period. The sediment supply produced by the glaciation combined with the persistent, high velocity flows during the glacial melt resulted in the filling of the incised valleys with very coarse-grained material including cobbles and gravel. The IVF deposits can act as fast paths, creating both lateral and vertical connections into the deeper, multi-channel aquifer systems. This will rapidly move water, over a time scale of weeks to years, from the surface to groundwater storage. These deposits, once located, would provide some of the natural infrastructure needed to develop a state-wide program to recharge California's groundwater. Weissmann et al (2005) have suggested that IVF deposits are present in the alluvial fans of a number of the major rivers coming into the Central Valley from the adjacent Sierra Nevada. In some areas, the IVF deposits are expected to be well preserved due to continued tectonic subsidence and compaction-driven past subsidence. The high priority in terms of recharge are shallow IVF deposits that are well-connected to the ground surface where water for recharge could be applied. The question we addressed: How can we efficiently find these deposits?
The extensive sand, gravel, and cobble units presumed to be present in IVF deposits are captured in driller's logs, but it would be challenging and timeconsuming to accurately delineate the extent of IVF deposits throughout the Central Valley using well data alone. In this study we explored the use of a helicopter-deployed airborne electromagnetic (AEM) method (Legault 2015) which can acquire approximately 200-400 line kilometers of data in a day, imaging the electrical resistivity to depths of hundreds of meters (depending on the electrical resistivity of the subsurface). As reviewed by Siemon et al (2009) the AEM method has been extensively used over the past 40 years for groundwater exploration. There has been a rapid increase in use over the last two decades due to the development of new AEM systems, with AEM data now acquired not just for exploration, but for detailed characterization of groundwater systems (e.g. Abraham et al 2012, Foged et al 2014, Christensen et al 2017, Minsley et al 2021. Of specific relevance to our study are examples of locating features similar in composition and geometry to IVF deposits: paleochannels (Chandra et al 2020(Chandra et al , 2021 and paleovalleys, also referred to as buried valleys (Gabriel et al 2003, Jørgensen et al 2003, Eberle and Siemon 2006, Best et al 2019. Also of relevance, given our study area in the Central Valley, are publications describing the successful adoption of the AEM method to map the large-scale structure of groundwater systems in the valley , Kang et al 2021, 2022 by transforming the measured geophysical property, electrical resistivity, to obtain images displaying the spatial variation in sediment type. As part of a larger study of the transition zone between the valley and the foothills of the Sierra Nevada in the Kaweah and Kings groundwater basins in the Central Valley, we acquired 542 line-kilometers of AEM data in the Kings River alluvial fan, focused on an area identified by Weissmann et al (2002) as containing an IVF deposit. If we were successful in this area in imaging the IVF deposit with the AEM method, the same could be done along the entire eastern side of the Central Valley, providing an efficient way to identify fast paths for recharging the groundwater systems of the valley.

Description of study area
The study area is a portion of the Kings River alluvial fan (as defined by Faunt (2009)) in the Kings subbasin in the southern part of the Central Valley. The area is shown in figure 1 along with the locations of the AEM flight lines. The geology of the subbasin has been described in publications by Page and LeBlanc (1969) and Weissmann et al (2002), and in the Central Kings Groundwater Sustainability Plan (Central Kings GSA 2019). The composition of the subbasin consists of older alluvium, mostly interbedded layers of silty or sandy clays, clay lenses, clayey and silty sands, gravels, cobbles, and boulders; younger alluvium with a similar composition; and flood-basin deposits. More than half of the groundwater system is composed of clay-and/or silt-dominated facies. The Kings River had a major impact on the geomorphologic landscape of the Kings subbasin forming one of the largest alluvial fans in the Central Valley (Burow et al 1997). The adjacent foothills are part of the basement complex of the Sierra Nevada, which is composed primarily of granite, with schist or serpentinite present north of the Kings River.
The genesis and presence of the IVF deposits was investigated in detail by Weissmann and Fogg (1999) based on lateral paleosol discontinuities observed in driller's logs. Evidence of one IVF deposit in the Kings River alluvial fan was then found in soil maps (taken from Huntington (1971)) and in cores from four wells, two of which were drilled directly into the IVF deposit. This IVF deposit, extending from the ground surface to 30 m below the surface, contains a thick, coarse sand and cobble unit between the depths of 22 m and 30 m, is at least 20 km long, and has an average width of 2 km (Weissmann et al 2004). The location of this IVF deposit was the target for our acquisition of AEM data and is shown in figure 1. The soil texture mapping used in the study of Weissmann et al (2002) also suggested other IVF deposits radiating from an intersection point located northeast of Sanger. As detailed in Weissmann's work on sequence stratigraphy of the fan, each of the multiple sequences that comprise the alluvial fan sediment column is expected to have an IVF deposit, so one would expect these IVF deposits to be stacked at the location where the Kings River emanates from the incised uplands.
The mapping of sediment texture in this area, as the percentage of coarse-grained sediments in 50-foot depth intervals, was conducted as part of developing the Central Valley Hydrologic Model (CVHM) (Faunt 2009). The CVHM texture data, shown over the depth interval of 0-30.5 m in figure 1, provide information on the distribution of sediment texture in the Kings River alluvial fan, but generally cannot resolve smaller-scale features such as the IVF deposits of interest in our work.

Compilation of well data and flight line planning
We had descriptions of the subsurface sediments at 270 locations from the work of Weissmann et al (2004) who had compiled existing driller's logs. In addition, we selected 550 driller's logs from the California Department of Water Resources Map Application (CDWR 2022a) that met the following criteria: had the global positioning system information required to determine the well location; had depths greater than 20 m, with preference given to deeper logs; and had one geologic description including a color and at least one texture description. (The last criterion was defined by Laudon and Belitz (1991) for assessing the quality of driller's logs). These logs were processed by transferring the descriptions into a digital, machine-readable format and the descriptions then automatically classified using keywords. For this study we worked with three major categories-sediments, weathered bedrock, and unweathered bedrock, with the sediments further divided into the two categories sand/gravel/ cobbles and clay. The category 'clay' was assigned to any description that included mention of clay. The locations of the wells corresponding to the driller's logs used in this study are shown in figure 1.
As can be seen in figure 1, we positioned closelyspaced flight lines (∼400 m spacing) over the identified IVF deposit. We placed additional flight lines in the area where the presence of other IVF deposits had been suggested. AEM data cannot be acquired over populated areas, so this impacted the location of the lines in the study area. To finalize the flight plan, flight lines were adjusted whenever feasible to be within 100 m of the digitized driller's logs. This was to allow use of the logs to support the interpretation of the resistivity model recovered from the AEM data.

AEM data acquisition, processing, and inversion
The AEM data were acquired in early December 2020 using the SkyTEM 304 system (Sørensen and Auken 2004), shown in figure 2 flying by the Pine Flat Dam on the Kings River. This is a time-domain system with a dual-moment transmitter which enables recording from ∼10 microseconds to ∼10 milliseconds thus providing good data quality from both shallow and deep regions of the subsurface. The raw AEM data were processed using Aarhus Workbench (Auken et al 2015) and stacked to sounding locations ∼35 m apart along each flight line. Processed data  and estimated data errors were exported and used as inputs for the subsequent inversion process. We used a spatially-constrained inversion (Viezzoli et al 2009)  to recover a 1D vertical resistivity profile at every AEM sounding location. In setting up each resistivity profile, a relatively fine discretization of 30 layers was chosen with the first layer being 3 m thick, and the last one being 28 m thick, summing to 280 m in total along the depth axis. A uniform resistivity model was adopted for both the reference and initial model; no additional a-prior information was incorporated. More information about the inversion method can be found in Oldenburg and Li (2005) and Cockett et al (2015). We determined the depth below which we had limited confidence in the recovered resistivity values, referred to as the depth of investigation (DOI), by conducting two separate inversions with two different reference models (100 and 200 Ω m) following the approach of Oldenburg and Li (1999). While there are other approaches to estimate DOI using a single model (e.g. Spies 1989, Christiansen andAuken 2012) we prefer to define the DOI by examining two models and identifying the depth above which there is consistency in the resistivity structure. All of the 1D vertical resistivity profiles were combined to produce, what we refer to as, a 3D resistivity model.

Obtaining the relationship between resistivity
and sediment type or lithology AEM data were acquired over regions of the subsurface composed of sediments and bedrock. Previous AEM studies in the Central Valley have shown that the location of the water table needs to be accounted for when interpreting the recovered resistivity model to determine the type of geological material that is present at a location (Knight et al 2018, Dewar andKnight 2020); this is because of the impact of saturation state on resistivity. Therefore, during the same month as the AEM survey, 73 groundwater monitoring wells in close proximity to the flight lines were resampled for water levels, in addition to the standard fall monitoring. This allowed us to define two zones in our resistivity model-above the water table and below the water table.
Exposed at the surface in the foothills is weathered and unweathered bedrock. Starting at locations with AEM data in the foothills, where the bedrock surface was at the top of the resistivity model, we followed the bedrock surface westward, out into the valley. The surface, which corresponded to the top of a very resistive zone, fell to depths below the DOI in the valley. The resistivity model was compared to the driller's logs in an attempt to determine the resistivity values that corresponded to weathered and unweathered bedrock.
In order to interpret the rest of the resistivity model in terms of sediment type, we used the method presented in Knight et al (2018). Application of this method involved building data pairs of the two data types: the resistivity values in the resistivity model and descriptions of sediment type, classified as either clay or sand/gravel/cobble and taken from the processed driller's logs. Each data pair consisted of a 'co-located' resistivity value and sediment type, where 'co-located' was defined as corresponding to the same depth interval with lateral separation distance between the two data types of less than 100 m. We had available 107 processed driller's logs less than 100 m from an AEM sounding which resulted in 396 data pairs above the water table and 373 below the water table. The method uses as input all resistivity-sedimenttype data pairs, assumes layering of the sediments, represents the AEM measurement as the application of an electric field parallel to sediment layers, and then solves using bootstrapping for the distribution of resistivity values corresponding to each sediment type. This approach assumes that there is not a variation in the total dissolved solids (TDS) in the water in the sediments that could also impact the resistivity. An evaluation of TDS measurements within the study area over the last 10 years suggests that this is a valid assumption (California State Water Resources Control Board 2022).

Interpretation of the resistivity model
The resistivity model, recovered through inversion of the AEM data acquired in the Kings River alluvial fan, is displayed in figure 3 as a 3D fence diagram. The DOI ranged from 100 m to 280 m, with an average depth of 130 m. Also shown is the resistivity model from the Kaweah subbasin, presented in a previous study , where the average DOI was 300 m. This deeper DOI is due to the fact that a different AEM system was used that has a larger moment. In order to facilitate comparison, we display both resistivity models to a depth of 130 m. The water table in the study area was found to vary from ∼7 m to ∼58 m below the ground surface, generally increasing in depth towards the west.
The bedrock surface was easily identified in the AEM data as the most resistive material, present at the ground surface in the foothills, then dipping over a short distance to depths below the DOI of the AEM data. We were unable to differentiate between weathered and unweathered bedrock due to the limited number of driller's logs and the observed overlap in the resistivity ranges of the two materials. Therefore, what we map as the bedrock surface is the top of undifferentiated weathered and/or unweathered bedrock.
Interpretation of the rest of the resistivity model was based on the relationships determined between resistivity and sediment type, shown in the zone above the water table in figure 4(a) and in the zone below the water table in figure 4(b). Also shown in these figures is the full distribution of observed resistivity values above the interpreted bedrock surface. The implemented method of Knight et al (2018) solves for the resistivity values of sediment layers that are typically much thinner than the vertical resolution of the AEM measurement; thus, it is not unexpected to see a large number of resistivity values falling in the range where the two resistivity distributions meet. Using the histograms in figures 4(a) and (b), we defined a cutoff between the two sediment types that assigns the more likely sediment type (higher value in the histogram) to a resistivity value. The final relationship that we determined between resistivity and sediment type is shown in figure 4(c). In each of the zones, above the water table and below the water table, all resistivity values from the cutoff to the maximum observed in the survey area were classified as sand/gravel/cobble, and all resistivity values from the cutoff to the minimum observed in the survey area were classified as clay/mixed. There can be uncertainty introduced into the classification due to uncertainty in the location of the water table. This was presumed to have minimal impact on our interpretation given the spatial density of the water table measurements. Not shown in figure 4(c) are the values corresponding to un/weathered bedrock. These ranged from 50-1200 Ωm, thus overlapping with the resistivity range corresponding to sand/gravel/cobble. As described above, mapping the bedrock surface utilized not just the presence of a sharp transition to high resistivity values, but also the spatial continuity of the surface as we traced it from exposure in the foothills to underlying the sediments in the valley.
We now provide a high-level interpretation of the 3D resistivity models displayed in figure 3. In the models from both this study and the previous study, we see in the foothills the pink color used to represent the highly resistive bedrock. As we move out into the valley, we see very low resistivity values, blue in color, suggesting clay-rich material and also see large areas with intermediate resistivity values, green in color, which most likely correspond to primarily fine-grained silt and clay with interlayered sand and gravel. What is most striking in figure 3 is the large package of more resistive material in the northern part of the Kings River alluvial fan, colors ranging from yellow through red, extending from the base of the Sierra foothills, out into the valley. While the resistivity model of the Kaweah subbasin shows the lower resistivity values associated with the finer-grained materials which make up much of the groundwater systems of the Central Valley, the resistivity model from the Kings River alluvial fan shows an area of more resistive, i.e. coarser-grained, materials. Within the alluvial fan, we see an extensive, highly resistive, linear feature in the area of the IVF deposit.

AEM imaging of the IVF deposit
In figure 5(a) we present a plan view of the integrated resistivity values in a shallow layer, corresponding to the depth interval that had been identified by Weissmann et al (2004) as the location of the IVF deposit. The shallow layer was defined as extending from the surface to a depth of 31 m or to the top of the bedrock, whichever was less. What is presented as integrated resistivity values was obtained by summing the products of resistivity and thickness over a depth interval and then normalizing by the total thickness of the depth interval. Areas of high integrated resistivity values indicate a high percentage of coarse-grained materials. A 20 km long, linear region of high resistivity values coincides with the mapped extent of the IVF deposit from Weissmann et al (2004) and continues towards the northeast. Our interpretation is that we are imaging the IVF deposit as an extensive, linear, highly resistive feature. Within this feature, most of the integrated resistivity values are over 100 Ωm, with a maximum resistivity value of ∼300 Ωm. We note that the added 100 Ωm contour in figure 5(a) agrees very well with the boundaries of the IVF deposit as mapped by Weissmann et al (2004). This level of contrast with the integrated resistivity values in the surrounding alluvial fan makes it possible to resolve the IVF deposit with the AEM data. The highly resistive feature, interpreted as the IVF deposit, dips gently towards the southwest and has a maximum width of about 3 km. As shown in the plan view in figure 5(b), where we display the integrated resistivity values from 41 m to 54 m in depth, the highly resistive zone is not visible below 41 m. We interpret the zone with lower resistivity values as underlying the base of the IVF deposit, in agreement with the observations of Weissman et al (2005).
Seen in figure 5(a) are other areas of high resistivity (>100 Ωm) to the southeast of the known IVF deposit. These could be the other IVF deposits in the area suggested by Weissmann et al (2002). There is also an area of high resistivity at the apex of the alluvial fan, in line with the trend of the main IVF deposit. This most likely corresponds to the stacking of IVF deposits where the paleovalleys from multiple glacial episodes converge.
We display a view of the resistivity model along the length of the imaged IVF deposit from west to east, A-A ′ , in figure 6(a) and two views across the deposit from south to north with B-B ′ in figure 6(b) and C-C ′ in figure 6(c); the locations of the sections are shown in figure 6(d). We also display the driller's logs that are within 100 m of the sections. The feature that we interpret as the IVF deposit coincides with high resistivity values in the location identified by Weissmann et al (2004) as an IVF deposit with descriptions of predominantly sand or gravel and a thick interval (>5 m) of cobbles. The potential value of the imaged IVF deposit for managed recharge is clearly evident in figure 6(a). The ∼30 m thick feature extends over 35 km from a surface exposure at the fan apex (at around 280 000 m in the easting direction) out into the valley, coming to the surface at other locations in the valley. Water applied to the ground surface at any of these surface locations would move rapidly into the subsurface, recharging the groundwater system. If there was concern about the presence of clay layers at the surface that could impede recharge but were too thin to be resolved in the AEM data, shallow boreholes could be drilled or cone penetrometer testing conducted to eliminate this possibility. The other highly resistive features seen close to the surface in the plan view in figure 5(a) are also seen in figures 6(b) and (c) over the same shallow depth interval as the mapped IVF deposit. These features could be IVF deposits but, regardless of origin, are large areas of high resistivity, indicative of coarse-grained materials very close the ground surface. Such features should also be prioritized as sites for managed recharge.
Beneath the surface layer containing the imaged IVF deposit and the other resistive features, we see ∼10-20 m of lower resistivity material, indicating the presence of fine-grained materials. Below this, there are areas of higher resistivity values, indicative of coarse-grained materials. It is possible that these are older IVF deposits consistent with the sequence stratigraphic characterizations of alluvial fans along the eastern side of the Central Valley depicted in Weissman et al (2005). Along the base of the sections, low-resistivity, clay-rich regions are seen. We do not see, in any of the sections, clear evidence of the water table. As was found in a study of AEM data from the Central Valley (Dewar and Knight 2020) the spatial heterogeneity in sediment type, and in the corresponding resistivity values, masks any change in resistivity due to the presence of the water table.  (c) is across the IVF deposit (C-C ′ ) with locations of sections shown in (d). The 'IVF boundary' corresponds to the extent of the IVF deposit given by Weissmann et al (2004). The black dashed line indicates the water table. Driller's logs within 100 m are shown. Intervals described as containing cobbles are denoted as solid white circles.

Conclusions
With the climate change predictions for California of increasing severity of floods and droughts, there is an urgent need to capture excess surface water during the flood years and use it to recharge the groundwater systems. A significant challenge is finding the optimal locations for recharge, where water spread on the ground surface can efficiently move downwards to recharge deeper sections of the systems. For surfacespreading recharge to be successful in the Central Valley, there is a need to locate sites where pathways of coarse-grained materials can move the water downwards, through the dominantly fine-grained materials that comprise much of the groundwater systems of the Central Valley.
IVF deposits, formed at the end of the last glacial period, offer the ideal natural infrastructure for implementing surface-spreading managed aquifer recharge along the eastern side of the Central Valley. The youngest of the IVF deposits will surface close to the valley edge, and then extend tens of kilometers out into the valley. These fast paths of coarse-grained materials are likely to be present in the alluvial fans of a number of rivers that enter the valley from the Sierras.
In this study, by acquiring AEM data with the dual-moment SkyTEM 304 system, we were able to resolve the distinct signature of the IVF deposit in the Kings River alluvial fan. The IVF deposit appeared in the AEM resistivity model as an extensive, linear feature of high resistivity values, allowing us to map its location within the alluvial fan. It is reasonable to presume that the other IVF deposits will have the same characteristics that were responsible for the distinct AEM signature of the Kings River IVF deposit: a linear feature extending tens of kilometers into the valley and a fill of sand, gravel, and cobbles that corresponds to high resistivity values. We therefore conclude that the AEM system employed in this study can be used to efficiently find IVF deposits along the entire eastern edge of the Central Valley due to their distinct appearance in the AEM resistivity model. By mapping the extent of the IVF deposits, the areas where they reach the ground surface can be identified as potential locations for surface-spreading managed recharge. These massive fast paths of coarse-grained materials are highly valuable natural infrastructure for recharging California's groundwater.
This study was conducted in California's Central Valley, but has significant implications for other mountain-bounded alluvial aquifer systems worldwide where past glaciation is likely to have formed similar IVF deposits. By mapping the location of IVF deposits, countries can put in place programs for recharge that integrate natural infrastructure into a traditional dependence on engineered infrastructure. Natural pathways of coarse-grained materials that can facilitate much-needed large-scale recharge efforts are present below the ground surface. The AEM method provides an efficient and cost-effective way to find them.