Mapping Aquifer Systems with Airborne Electromagnetics in the Central Valley of California

The passage of the Sustainable Groundwater Management Act in California has highlighted a need for cost‐effective ways to acquire the data used in building conceptual models of the aquifer systems in the Central Valley of California. One approach would be the regional implementation of the airborne electromagnetic (AEM) method. We acquired 104 line‐kilometers of data in the Tulare Irrigation District, in the Central Valley, to determine the depth of investigation (DOI) of the AEM method, given the abundance of electrically conductive clays, and to assess the usefulness of the method for mapping the hydrostratigraphy. The data were high quality providing, through inversion of the data, models displaying the variation in electrical resistivity to a depth of approximately 500 m. In order to transform the resistivity models to interpreted sections displaying lithology, we established the relationship between resistivity and lithology using collocated lithology logs (from drillers' logs) and AEM data. We modeled the AEM response and employed a bootstrapping approach to solve for the range of values in the resistivity model corresponding to sand and gravel, mixed coarse and fine, and clay in the unsaturated and saturated regions. The comparison between the resulting interpretation and an existing cross section demonstrates that AEM can be an effective method for mapping the large‐scale hydrostratigraphy of aquifer systems in the Central Valley. The methods employed and developed in this study have widespread application in the use of the AEM method for groundwater management in similar geologic settings.


Introduction
To ensure a secure food supply, many agricultural regions of the world require reliable sources of surface water and/or groundwater to meet irrigation needs. The Central Valley of California is one such region. Covering 50,000 km 2 , bounded by the Sierra Nevada to the east and the Coast Ranges to the west, the valley yields a third of the produce grown in the United States valued at $17 billion dollars per year (U.S. Geological Survey 2016). Much of the irrigation water in the valley has historically been taken from surface water-the rivers, lakes, and reservoirs replenished by winter storms and the spring melting of the snow pack in the Sierra Nevada. This water is provided through a series of federal, state and private irrigation canals, dams and lakes. In times of drought, most recently in the periods 2007 to 2009 and 2012 to 2016, when surface water deliveries were substantially reduced, the only way to meet irrigation needs has been through extensive pumping of groundwater. This has exacerbated an already serious problem in the Central Valley, where some areas have experienced declining water levels for several decades. Drought conditions have led to further lowering of groundwater levels, due to increased pumping and decreased recharge, causing wells to go dry and, in some areas, significant subsidence of the ground surface. The overdraft has been so significant, that there are now approximately140 million acre-feet (MAF) of unused groundwater storage space in the Central Valley, a value calculated based on the estimated difference between predevelopment and current conditions (The Nature Conservancy 2016). In contrast, the total surface water storage capacity in California is 42 MAF.
The alluvial sedimentary geology of the Central Valley is typically composed of more than 50% to NGWA.org Vol. 56, No. 6-Groundwater-November-December 2018 (pages 893-908) 70% fine-grained deposits dominated by silt and clay beds. These silts and clays make the groundwater basin confined to semi-confined, significantly impeding vertical groundwater flow. Within this geologic system are networks of sand and gravel that both constitute the aquifer system and provide pathways for recharge. The sustainable management of the groundwater resources of the Central Valley, now required with the passage of the Sustainable Groundwater Management Act (SGMA) by the California Legislature in 2014, highlights the need to better understand the hydrostratigraphy of the subsurface to improve the conceptual models of the aquifer systems and provide the information needed to inform management decisions. Mapping the locations of the sand and gravel networks would make it possible to select the best locations for surface spreading techniques so that recharge could be dramatically increased, and repressurization of the confined aquifer networks could be accomplished.
Mapping of the aquifer systems could also help assess the vulnerability of an area to subsidence associated with groundwater withdrawal, as subsidence in the Central Valley occurs primarily in areas where there are numerous clays layers interlayered with the sands and gravels (Faunt and Sneed 2015). In addition, accurate conceptual models could guide the siting of expensive monitoring wells.
The key question is then: How do we map the aquifer systems in the Central Valley in a cost-effective way? We propose regional implementation of the airborne electromagnetic (AEM) method to obtain information needed about the hydrostratigraphy of the aquifer systems. There are various AEM methods available, all of which use an airborne platform to move geophysical sensors over an area to map out subsurface variations in electrical resistivity to depths of approximately 500 m (depending on the electrical resistivity). Vertical resolution is on the order of meters to tens of meters. Lateral resolution along the flight lines is typically 30 m; between the flight lines it is determined by the line spacing, which is set according to the objectives of the survey and can range from 250 m to 30 km. Given the relationship between electrical resistivity and the properties of the subsurface geological materials, AEM data can be used, along with existing lithologic information available from drillers' logs, to develop a model of the large-scale (tens to hundreds of meters) hydrostratigraphic packages that define the architecture of the aquifer systems. The AEM method has been used throughout the world for groundwater exploration and aquifer mapping (e.g., Sattel and Kgotlhang 2004;Podgorski et al. 2013). The various methods that have been used to develop the hydrostratigraphic model are described in a recent paper by Christensen et al. (2017) and include both knowledge-driven cognitive approaches (e.g., Jørgensen et al. 2013) and geostatistical approaches (e.g., Marker et al. 2015).
We conducted a pilot study to assess the use of AEM to map the aquifer systems in the Central Valley of California. There had been one private survey completed, so this represented an opportunity to acquire the first publicly available dataset in the valley. We elected to use the SkyTEM system. An important factor in this selection was the ability of the system to provide the resolution required within the study area; studies have shown that the SkyTEM system can provide good resolution in similar environments (Bedrosian et al. 2015). The SkyTEM system, originally developed by Aarhus University (Sørensen and Auken 2004), was specifically designed with a dual moment transmitter that allows for near-surface as well as deep penetration. SkyTEM has been used over the past 15 years in Denmark to acquire 60,000 line kilometers of data covering 15,000 km 2 . Together with the airborne sensor, processing software was developed and a national data repository established that enables access to processed data for water authorities, water suppliers, municipalities, and others (Thomsen et al. 2004;Auken et al. 2009;Møller et al. 2009).
The demand for acquiring these AEM data in Denmark has been driven by groundwater legislation that requires that all municipalities characterize and manage the groundwater systems. AEM has been found to be the most cost-effective way to get the subsurface coverage needed for compliance with this regulatory requirement; a significant benefit being that there is no need for land access. In California, SGMA requires that local authorities assess the state of their groundwater basins and develop plans for the sustainable management of groundwater. The acquisition of AEM data throughout California could play a central role in providing critical information that can inform the development of groundwater models, guide recharge efforts, assess geologic controls on observed subsidence, and aid in siting of monitoring wells.
The location selected for our study, shown in Figure 1, is in the Tulare Irrigation District in the San Joaquin Valley that, with the Sacramento Valley to the north, makes up the Central Valley. There have been extensive, chronic groundwater level declines and associated problems such as subsidence in this area between 2007 and 2009 (Farr and Liu 2014;Smith et al. 2017). The geology in this area is typical of that found throughout the Central Valley, so AEM performance at this location should be representative of how it would perform elsewhere in the valley. The first question we posed was related to the quality of the acquired data: how effective would the AEM method be in imaging beneath the shallow, electrically conductive clays commonly found in the Central Valley? The second question we posed was related to the interpretation of the resistivity models obtained through inversion of the acquired AEM data. In Denmark, the 15 years of working with AEM data have involved extensive analysis of resistivity data from geophysical measurements (AEM, ground-based and borehole) and lithologic information from wells to build a national atlas that can be used derive lithology information from resistivity measurements (Møller et al. 2009;Christiansen et al. 2014;Barfod et al. 2016). Nothing of this sort currently exists in California. The question we posed: Could we work with the resistivity models and the limited well data to obtain information about the architecture of the aquifer systems, mapping out key hydrostratigraphic packages? Answering these two questions, related to the ability to acquire and then interpret high quality data, was an essential first step in determining whether the AEM method could be reliably used for mapping the aquifer systems in the Central Valley. If so, this would be a new way of obtaining critical information needed for the implementation of SGMA, and could significantly transform the approach to groundwater management in California. While the focus of this study was the Central Valley of California, the findings and the new methodologies developed will have much broader impacts, advancing the use of AEM for imaging alluvial aquifer systems throughout the world. District; we derived this image from Interferometric Synthetic Aperture Radar (InSAR) data. InSAR is a method that can provide the change in ground elevation over time with accuracy on the order of millimeters to centimeters (Madsen and Zebker 1998;Rosen et al. 2000). The change in phase measured between two satellite passes is used to form an interferogram from which the change in elevation is calculated. To obtain the subsidence map in Figure 2, we used data from the Advanced Land Observing Satellite (ALOS), which has an L-band radar system with a wavelength of approximately 24 cm. This relatively long wavelength is more effective in agricultural areas, because it is not as easily decorrelated by vegetation. We processed 82 interferograms over the time frame 2007 to 2010 and used the small baseline subset method (Berardino et al. 2002) to calculate the relative motion over time. We then used this time series to solve for the mean velocity during the study period. As can be seen in Figure 2, subsidence reached a maximum value of 26 cm/year between 2007 and 2010.

Description of Study Area
Our flight lines, shown in Figure 2, covered 104 line-kilometers and were selected to go from the center of the subsidence bowl to an area with little, to no, subsidence. Survey flight planning included adherence to U.S. Federal Aviation Administration (U.S. FAA) regulations regarding flying and towing cargo over infrastructure, including not flying over buildings and large highways and flying well above power lines. The lines were spaced approximately 1 km apart in order to conduct a small-scale reconnaissance of the area of interest.
There are numerous wells in the study area, many with drillers' logs that describe lithology; these were obtained from the California Department of Water Resources. The majority of the wells (∼80%) are less than 100 m deep, but we were able to use lithology (drillers') logs from 12 wells (shown in Figure 2) within approximately 500 m of the flight lines that reach at least 150 m. The lithology logs describe a shallow sand and gravel aquifer, overlying a clay layer, 0 to 20 m thick, referred to as the Corcoran Clay. Beneath the Corcoran Clay, the limited well data suggest another sand and gravel aquifer unit with interlayered clays. There is one resistivity log within 1 km of a flight line, from well 20S23E14, which shows the location of the Corcoran Clay, as well as the alternating sand to clay nature of the upper and lower aquifer.
In addition to the lithologic logs, we had a report with cross sections, provided by the Kaweah Delta Water Conservation District (Fugro West 2007), which described upper and lower aquifer units that are separated by the Corcoran Clay. The upper aquifer is described as Quaternary older alluvium (oxidized) and the lower aquifer is described as Quaternary older alluvium (reduced). Both aquifers are interpreted to have numerous interbedded silts and clays, with the lower aquifer interpreted to have more fine-grained material than the upper aquifer (Fugro West 2007). The Corcoran Clay, described in the report as Quaternary lacustrine and marsh deposits, pinches out in the eastern part of the section. At the base of the lower aquifer is an impermeable unit described as Pliocene and Pleistocene (questionable) continental deposits.

Acquiring the AEM Resistivity Model: Methodology and Results
The AEM method has been used for many years to map geology (Palacky 1981) and, in the last 10 or so years, has been widely used to map groundwater prospects. The theory behind the method is described in Ward and Hohmann (1988). In the SkyTEM system used in this study, all of the hardware required for data acquisition is suspended beneath a helicopter, which moves over the land surface at approximately 80 to 100 km/h with the frame hanging approximately 30 m above the ground surface. Current flowing in a transmitter loop generates a primary magnetic field. The termination of current causes a time-varying decay in the produced magnetic field which causes eddy currents to flow at various depths beneath the land surface. The less electrically resistive the region, the stronger the current and the more slowly the current decays. The eddy currents generate their own secondary magnetic fields which are measured at the receiver mounted on the transmitter loop. The strength of the field and the time dependence contains information about how current flows through the ground. The acquired data are inverted to obtain a model of the spatial variation in electrical resistivity of the subsurface material.
Airborne data were acquired using the SkyTEM 508 system on October 27, 2015. Although data were acquired approximately every 2.5 m along each of the flight lines, the dual moment mode, in which the data were acquired and averaged, resulted in an effective EM sounding spacing of 30 m along each line.
The acquired data were inverted using the Aarhus Workbench (HydroGeophysics Group 2011). The procedure applied to the data processing, including noise assessment, follows that detailed in Auken et al. (2009). We then performed laterally constrained inversions ) where correlation along the flight line is imposed on the inversion. This was followed by spatially constrained inversions, defined in Viezzoli et al. (2008) as a methodology to impose correlation across survey lines. The values of the parameters used to enforce spatial coherency, from equation 12 in Viezzoli et al., were A = 1.5, B = 30 m, b = 0.75. These were determined based on both prior geological knowledge and empirically, that is, testing a few variations while analyzing data misfit. For the starting model we used a resistivity everywhere of 30 m. While no lithology logs or electrical logs were used as constraints in the inversion (the electrical log was too far, >150 m, from the flight lines), they were used postinversion for qualitative analysis and comparison with the AEM inversion results. Figure 3 presents the resulting resistivity model for Line 3 (the position of which is shown in Figure 2). The Y position in this figure is the distance along the survey line. Gaps in the image are due to electromagnetic coupling of the AEM acquisition system with ground interference such as power lines and cathodically protected pipelines. The affected soundings were filtered and then manually edited out of the data set prior to inversion. Each pixel in the resulting resistivity model has an assigned resistivity value. The size of a pixel corresponds to the spatial resolution, both horizontal and vertical, at that location. Horizontal resolution along the flight line is based on the AEM sounding separation (about 3 m) and the number of soundings averaged during processing (10 soundings on average). In this study the processed sounding separation was 30 m. Large gaps due to decoupling of electromagnetic noise will locally negatively impact the horizontal resolution. The ability to resolve features decreases with depth per the fundamental physics of the technique. As depth increases, typically each pixel in the resistivity model, for this acquisition system, increases in thickness by approximately 1.1 times the thickness of the previous layer. In the resistivity model shown in Figure 3, the pixels at the top of the section have a vertical dimension of 3 m, increasing to approximately 14 m at 100 m depth, and approximately 48 m at 400 m depth. This approach, of inverting for electrical resistivity in a multilayer model with layers of fixed thicknesses, is commonly referred to as smooth model. We elected to use this type of model because the many layers are needed to accommodate the complexity of the geologic setting (e.g., the discontinuous clay layers). In addition, this approach allows us to explore, relatively easily, possible subtle variations in the electrical resistivity structure. The depth of investigation (DOI), as applied here, plays a critical role in data interpretation. The DOI is a numerical estimation of the depth below which the resolution of a model, obtained from the numerical inversion, diminishes. This loss of model resolution is primarily due to the limitations of the method including: the equipment and methodology used for data acquisition, the noise in the data, and the resistivity structure (the geology) being investigated (Christiansen and Auken 2012). The Aarhus Workbench uses inverted model sensitivities (components of the inversion Jacobian matrix) determined during the last inversion step, along with an estimate of data uncertainty, to calculate the DOI's (Christiansen and Auken 2012). The Jacobian matrix is the sensitivity of the data to perturbations at the current location in model-space. As such, the Jacobian is also a representation of the sensitivity of changes in the model to the observed data. The last step in this process is calculation of the cumulative sensitivities, produced by summing up the individual thickness-normalized sensitivities, starting with the bottom layer. Two cumulative-sensitivity thresholds are selected to represent what are known as the "upper" and "lower" depths of investigation. We selected cumulativesensitivity thresholds of 0.6 and 0.1 for calculation of the upper and lower DOI's. In Figure 3 the upper estimate of the DOI is represented as a dotted line and the lower estimate of the DOI coincides on this particular line with the base of the displayed profile section.
In Figure 3 notice that the upper DOI ranges in elevation from a maximum of about −300 m down to an elevation of about −375 m. This is not a large range at these elevations. It is a bit shallower over the resistive zone on the southern end (the left side) of Line 3 and a little deeper on the northern end (right side) of Line 3. Typically, the DOI is deeper in resistive material and shallower in conductive materials. The key to interpreting the upper DOI in Figure 3 is to examine the resistivity values in the first 200 m. On the south end of the line there is more conductive material between elevations 0 to −100 m. On the northern end of the line, there is more resistive material between elevations of +75 m and −50 to −100 m. These differences in the first 200 m of the subsurface are influencing the character of the upper DOI. The lower DOI, at the bottom of the section along this line, displays an irregular character along the northern end of Line 3, indicating that the model's sensitivity to the resistivity values at depth in those locations is quite low.
For the depth of the water table (shown in Figure 3) we used the interpolation of water level data, measured between September 1, 2015 and November 25, 2015, provided by the California Department of Water Resources, 2016, Groundwater Information Center Interactive Map Application: https://gis.water.ca.gov/app/gicima/. The depth of the water table was estimated to vary in the study area between about 47 and 60 m. Given that the distance between wells used for the interpolation was on the order of kilometers, variations in the water table at the subkilometer scale will not be captured in these data. Most of the highest resistivity values in the profile section presented in Figure 3 are found in the region above the estimated water table. This can be seen in a comparison of the histogram of all the AEM inversion resistivity values from the region above the water table (Figure 4a) with the histogram of values in the region below the water table (Figure 4b).
Below the water table, resistivity values range from 6 to 43 m. Above the water table, resistivity values range from 8 to 150 m; we display in the figure the counts out to 60 m as only 2% of the data had resistivity values greater than 50 m. The difference in resistivity values above and below the water table becomes important in the interpretation of the resistivity model.

The Resistivity-Lithology Relationship: Methodology and Results
The relationship between resistivity and lithology is at the core of the use of any geophysical electromagnetic or resistivity method to map out the variation in subsurface lithology. The extensive literature studying the resistivitylithology relationship is reviewed by Knight and Endres (2005) and briefly summarized here, highlighting what we expect to be the dominant factors linking resistivity and lithology in our study area.
In sediments and sedimentary rocks with water in the pore space, the primary mechanism for electrical conduction is typically ionic conduction through the pore water. As a result, the electrical resistivity tends to decrease as the volume of water-filled porosity (equivalent to volumetric water content) increases and will also decrease as the salinity of the water increases. In the main aquifers of the Tulare Irrigation District, there are no reports of significant variation in salinity of the pore water so, in this study, we assumed that pore water chemistry does not affect electrical resistivity. Further work is needed to incorporate water quality into the interpretation. study area, if we consider only electrical conduction due to ionic conduction through the pore water, we would expect to see, as observed, higher resistivity values above the water table where the sediments are unsaturated. We would also expect to see resistivity decrease as the porosity of the material increases.
A second mechanism that contributes to electrical conduction in sediments is surface conduction, which occurs due to the presence of a high concentration of ions associated with the electrical double layer at the solid/water interface. While surface conduction can occur in any material, it is most significant when clays are present due to their high surface area. The presence of surface conduction causes electrical resistivity to decrease as the surface-area-to-volume ratio increases; that is, resistivity will decrease as the grain size decreases and as clay content increases. For the purposes of our study, if surface conduction contributes to the measured electrical resistivity, we would expect to observe the highest resistivity values in the gravels and sands, with resistivity decreasing in silts, and further decreasing in the clays.
Given the absence of any database comparing resistivity and lithology data in California, we began our analysis of the resistivity-lithology relationship with a review of all available geological data and then a comparison with the resistivity values in the AEM resistivity model. This established method, referred to as Method 1, involves the steps outlined in Jørgensen et al. (2003) and Høyer et al. (2015), and involves the development of a conceptual geological model followed by comparing lithological and borehole geophysical logs to the AEM resistivity model. A  West 2007). The lithology logs presented in that report, showing various units of clays, sands, and gravels, in combination with the 12 lithology logs we had available, provided key lithology information as well as location information of representative lithologies to which the resistivity models could be correlated. Based on these resistivity-lithology correlations, Table 1 was developed. Note, in defining the lithology categories, we used the terms seen in the Fugro West report and the drillers' logs. The ranges in the resistivity values presented in Table 1 encompass the information from the AEM data, published cross sections, and lithology logs but are limited by the relatively small number of lithology logs close to the survey flight lines, the shallow depth of some of the lithology logs, and by the difficulty in separating out the resistivities above the water table from below the water table. We therefore decided to explore a second approach.
Our second approach, referred to as Method 2, involved the development of a new methodology that used the 12 lithology logs and the AEM resistivity model to quantitatively solve for the ranges of resistivity values that correspond to defined lithologies in the study area. At the core of this method is a key point: What we wanted was the relationship between the AEM-measured resistivity and lithology, not between resistivity measured by some other method (e.g., resistivity logging) and lithology; so we used the AEM resistivity values. One of the limitations in our approach, common to any other approach that could be used in this study area, was the lack of information about lithology below approximately 150 m. We assumed that the established resistivity-lithology transform was valid at greater depths but this is an issue that requires further study.
We conducted a separate analysis for the regions above and below the water table; a change in water content will result in a large change in the resistivity of a material so it is important to develop the resistivitylithology relationship in a way that accounts for this. Based on a review of the lithologic logs, and given the sparsity of independent lithologic information, we defined three lithologies below the water table: sand and gravel, mixed fine and coarse, clay. There was limited mention of mixed fine and coarse materials in the shallower section, so we reduced the lithology categories in the region above the water table to consider only sand and gravel, and clay. We converted the descriptions in all lithology logs to these categories. The intervals that we classified as "mixed fine and coarse" had variable descriptions: silt, sandy clay, silty sand and fine sand. We then systematically worked through the 12 lithology logs, from ground level to the deepest layer described, and identified the pixels in the AEM resistivity model that were closest to the location of the layers described in the lithology logs. A schematic of the basic steps in our approach, Method 2, is shown in Figure 5. On the left we have the layers described in a lithology log; on the right we have the closest AEM pixel. What is shown here is not an actual log but given as an example to explain our methodology using the three lithologies defined below the water table. Starting at the top of the lithology log, each layer has an assigned lithology (sand and gravel, mixed fine and coarse, or clay) and a thickness t i where i corresponds to the number of the layer. In the logs, the thickness of each described lithologic layer typically ranges from 1 to 2 m. We assigned to each layer the resistivity corresponding to the lithology: ρ sg , ρ mixed , and ρ clay for the resistivity of sand and gravel, mixed fine and coarse, and clay, respectively. It is important to note that while we have used a single variable to represent the resistivity of each lithology, there will always be a range of resistivity values due to spatial variation in the water content, composition, and pore structure.
On the right in Figure 5 is a pixel from the AEM resistivity model that is closest to the well for which we have the lithology log. In this study, we worked with all of the depth intervals in the 12 lithology logs and found that the closest AEM pixels were located on average 250 m away from the wells, with the separation distances ranging from 50 m to 1 km. We made the reasonable assumption that the lithology sampled in the AEM measurement was that shown in the lithology log. As shown in Figure 5, the AEM measurement does not have the vertical resolution to resolve the resistivity structure at the scale of the individual layers. What is derived from the AEM measurement is a larger-scale resistivity value referred to as ρ AEM . The vertical dimension of the AEM pixel, denoted in the figure as t AEM , is the vertical resolution of the AEM measurement at that depth. This varied from approximately 3 m at the surface to approximately 10 m at the depth of the deepest layers used in this analysis. We set up a relationship between ρ AEM and the resistivity values in the corresponding layers. The physics of the AEM measurement results in a form of averaging of the resistivity values in the individual layers that can be described, to first order, by the following relationship between the layer resistivity values and ρ AEM : where for n layers, i , as defined above, refers to the layer and ρ i is the layer resistivity which will be ρ sg , ρ mixed , or ρ clay . This relationship can be derived by representing NGWA.org R. Knight et al. Groundwater 56, no. 6: 893-908 each layer by a resistor with resistance R i = (ρ i L)/(t i W ), where L is the length of each layer in the direction parallel to the orientation of the field lines and t i W is the crosssectional area of each layer perpendicular to the field lines. The orientation of the electric field lines during the AEM measurement is such that the total measured resistance, R AEM = (ρ AEM L)/(t AEM W ), can be estimated by adding the "layer resistors" in parallel. This assumes that the field lines are parallel to the layering; this is a reasonable approximation.
Assuming constant values of resistivity for ρ sg , ρ mixed , and ρ clay , the above equation can be re-written for every depth interval where we have layers described in a lithology log and a nearby AEM resistivity value, as follows: ( Working with all of the layers described in the 12 lithology logs, we had 75 such equations with the three unknowns, ρ sg , ρ mixed , and ρ clay , below the water table and 74 equations and two unknowns, ρ sg , and ρ clay , above. We randomly sampled (with replacement) the system of equations 1000 times, where the sample size at each of the 1000 iterations was equal to the parent sample size of 75. At each iteration, we used this random sample of equations to solve for ρ sg , ρ mixed , and ρ clay , producing a distribution of 1000 possible values of resistivity for each lithology. This method, known as bootstrapping (Efron and Tibshirani 1994), yielded a distribution of resistivity values for each lithology. The distribution obtained through the bootstrap analysis shows the uncertainty in determining the relationship between the resistivity values in the derived resistivity model, and lithology. We note that, in this case, the uncertainty includes data uncertainty and model (also referred to as epistemic) uncertainty. Our model, given in Equation 1, represents the resistivity of each lithologic unit with a single variable, yet we know that this variable has a range in values due to spatial variation in properties as noted above.
For the region below the water table, the distribution of resistivity values found for the three lithologies are shown in Figure 6:11 m < ρ clay < 18 m, 12 m < ρ mixed < 22 m, and 17 m < ρ sg < 34 m. Because we were limited in our analysis to working with 12 lithology logs, we did not sample the complete range of resistivity values for the three lithologic units. The histogram of resistivity values below the water table (Figure 4b) contains resistivity values that range from 6 to 43 m, yet the lowest clay resistivity value in Figure 6 is 11 m and the highest sand/gravel resistivity value is 34 m. Incorporating the resistivity values determined through inversion of the AEM data with those derived from the bootstrap method, we therefore set the lower bound on the resistivity of clay to be 6 m and the  Table 2: 6 m < ρ clay < 18 m, 12 m < ρ mixed < 22 m, and 17 m < ρ sg < 43 m.
The results of the bootstrap analysis for the region above the water table are shown in Figure 7: 22 m < ρ clay < 31 m and 25 m < ρ sg < 35 m. As was the case in our analysis of data from below the water table, the 12 lithology logs did not sample the full range in resistivity values in the region about the water table. In the AEM resistivity model, values were found that ranged from 8 to 150 m. We therefore again incorporated the AEM inversion results and so reduced the lower bound on the range of resistivity values for the clay and increased the upper bound on the range of resistivity values for the sand and gravel. Table 2 provides the final results: 8 m < ρ clay < 31 m and 25 m < ρ sg < 150 m. Reviewing all of the results compiled in Table 2, it is very clear the importance of determining two independent resistivity-lithology relationships, one valid for the region above the water table and one for the region below the water. Using the relationship between resistivity and lithology determined in the region below the water table would lead us to interpret large areas in the unsaturated zone as sand and gravel, when in fact they are likely to be unsaturated clays.

Transforming The Resistivity Model to Lithology: Methodology and Results
We used the results obtained from Method 2 to transform the AEM resistivity model to lithology. As can be seen in Figures 6 and 7   be interpreted, with a high degree of confidence, to be a specific lithology. For example, above the water table, ρ AEM > 31 m can be defined as an unsaturated sand and gravel, and ρ AEM < 25 m can be defined as clay, but for measured resistivity values in the range 25 m ≤ ρ AEM ≤ 31 m we cannot determine whether the lithology is clay or sand and gravel. Beneath the water table similar uncertainty exists. We can interpret as sand and gravel those areas where ρ AEM > 22 m, and interpret as clay those areas where ρ AEM < 12 m, but there are areas where uncertainty exists in terms of defining lithology. We elected to honor this uncertainty by displaying in Figure 8a and 8b the results in terms of the AEM resistivity values, showing on the color bar the correspondence between resistivity and lithology in the regions above and below the water table, using the two independent resistivity-lithology relationships that we have defined. We note that in displaying the data we have elected to bin all resistivity values above 50 m (which includes only 2% of the data) so as to expand the color bar. Figure 8a shows the results from Line 3 (location in Figure 2); Figure 8b is a fence diagram displaying all of the AEM survey results. We interpret, and label in Figure 8a, the relatively thick, continuous conductive layer at approximately 100 m depth on the left side of the section to be the Corcoran Clay; the conductive feature thins and terminates on the right side of the section. The Corcoran Clay is known, from drillers' logs, to be present at this depth in this area, becoming discontinuous in the northeast (the right side of Figure 8a). In addition to this presentation of the results we also transformed the resistivity model to lithology in two other ways that communicate information about the lithologic variation in the subsurface mapped with the AEM method. The three images in Figure 9 display the probability of sand and gravel (Figure 9a), mixed fine and coarse (Figure 9b), and clay (Figure 9c) at any given location. We created these images by taking the distributions shown in Figures 6 and 7 as indicative of the probabilities of finding a lithology at a location. The final way in which we display our results is given in Figure 10, which shows the most probable lithology at any location along Line 3 as a single profile plot, and Figure 11, which displays the full interpreted data set as a fence diagram. The conductive feature labeled in Figure 8a as Corcoran Clay, appears as the thick clay unit in Figures 10 and 11, at a depth of approximately 100 m, present in all sections in the southwestern part of the study area.

Discussion
The motivation for this study was to determine whether the AEM method could be used for mapping the aquifer systems in the Central Valley, recognizing that the findings here would have implications for using the AEM method in characterizing other alluvial aquifers. Let us consider the first question we addressed, related to the imaging capability of the method: How effective would the AEM method be in imaging beneath the shallow, electrically conductive clays commonly found in the Central Valley? Our concern had been that electrically conductive clays would limit the DOI. We obtained AEM data of excellent quality to a depth of approximately 400 m along all of the flight lines. Neither the thick Corcoran Clay, nor the numerous fine clay layers described in the lithology logs, negatively impacted the penetration depth of the measurement. The calculation of the DOIs determined an upper DOI between 300 and 400 m and the lower DOI at approximately 500 m. Somewhere between the upper and lower DOI's, and usually deeper than the lower DOI, accurate resolution of the true electrical resistivity of the subsurface region being sampled is reduced. At this depth the imaging, through numerical inversion, transitions from accurately resolving the values of electrical resistivity to detecting changes in the electrical resistivity. Note that while detection indicates that the true magnitude of the electrical resistivity may not be accurately resolved; there is still some sensitivity to changes. Somewhere deeper than the lower DOI is generally the line of demarcation between resolution and detection. It is different in every data set due to differences in data noise and the electrical resistivity structure. The lower DOI is usually taken as a depth to which there is very high confidence in the resolved resistivity values. The lower DOI usually represents a depth at which the transition occurs from resolution to detection. As an example of applying the concept of DOI in determining the level of confidence in a resistivity model, consider the high resistivity unit seen in Figure 8a on the southwestern end of Line 3 and in the fence diagram in Figure 8b. This unit falls between the upper DOI (the dotted line) and the lower DOI, approximately coinciding with the base of the displayed model. The fact that the high resistivity unit begins at approximately the depth of the upper DOI does not mean that it should be questioned as a potential geologic unit. Anything below the lower DOI is in the transition range between resolved and detected. So while the resistivity of the high resistivity unit might not be accurately determined (i.e., the resistivity might not be exactly 40 m) below the lower DOI, it will have a detected range of approximately 30 to 60 m; this provides a high level of confidence in its identification as a high resistivity unit. The finding in this study, that we can accurately resolve the resistivity model to a depth of approximately 500 m and detect changes in resistivity to greater depths, is an important result for evaluating and planning AEM surveys for subsurface mapping in other parts of the Central Valley.
In addition to determining to what depth we can resolve/detect changes in resistivity, it is also important to consider how well we can capture vertical changes in resistivity and the implications for estimating the thicknesses of various units. The Corcoran Clay is the main confining unit between the upper and lower portions of the aquifer system, and can be seen in the resistivity sections in Figure 8 and in the interpreted sections in Figures 10 and 11. In order to determine how accurately the AEM method can determine the thickness of the Corcoran Clay, we generated a simple model of the electrical resistivity variation across the clay, and then modeled the response of the SkyTEM system and inverted the data, to compare the inverted clay thickness to the "true" clay thickness.
In Figure 12 we show, as the red line, a simplified representation of the true electrical resistivity variation across the Corcoran Clay, using as "true" the resistivity values from an electrical log for well 20S23E14, which is located about 300 m from the flight lines. Above, below, and within the Corcoran Clay, we set the resistivity to constant values, equal to the average values seen in these zones in the electrical log. We generated synthetic AEM data by modeling the response of the SkyTEM 508 system to the variation in resistivity values. We then inverted the synthetic data, using the same inversion process as was applied to the full data set in this study, so that the layer thickness was set to 3 m at the surface and then increased by 10% with depth (so that each layer thickness was 1.1 times that of the previous layer). This yielded the resistivity variation that would be seen in the AEM resistivity model, shown as the blue line in Figure 12.
While the inverted resistivity values in Figure 12 do not perfectly match the true resistivity values, they do  capture the general trend, with the Corcoran Clay clearly imaged as being conductive. The basic guideline for electromagnetic modeling is that there will be sensitivity to a layer whose thickness is at least 10% of the depth to the layer. The resistivity values may not be highly accurate, but there will be evidence that a conductive zone is present.
It is important to note the smoothing seen in the inversion result that masks the abrupt changes in resistivity. There are constraints in the inversion routine, on the allowed change in resistivity between adjacent layers, so that the inverted AEM data show a more gradual change in resistivity than is present. As a result, the inverted Corcoran Clay thickness is greater than the true thickness. This is likely to be a general result, suggesting that in all of our sections displaying the interpreted variation in lithology, the thickness of the Corcoran Clay will be overestimated if it is thinner than 10% of the depth to the middle of the layer in which it is indicated. Sharper vertical boundaries in the resistivity models could be achieved with different inversion strategies such as the minimum gradient support or with a few-layered inversion that solves for the resistivity and thickness of each layer (Vignoli et al. 2017). Adding a-priori information to the inversion could also improve the accuracy of the spatial extent of the recovered conductive layer (e.g., Sapia et al. 2014).
Once a resistivity model was obtained, we addressed the second question: Could we use the resistivity models and the limited well data to obtain information about the architecture of the aquifer systems, mapping out largescale hydrostratigraphic packages? The first step involves transforming the resistivity values to lithology. We used two methods. Method 1 is an established approach that provided us with a range of resistivity values for each lithology. Method 2 was designed to specifically link the SkyTEM-measured resistivity to lithology. As part of this method, we further refined the range of values by separating the region above and below the water table and using a bootstrapping method to obtain the distribution of resistivity values for each lithology.
Comparing first the resistivity range determined for the sands and gravels, we find that both methods gave very similar results, with the total range from Method 1 being 19 to 150 m, and from Method 2 being 17 to 150 m. The resistivity range found for clay below the water table with Method 2 (6 to 18 m) includes close to the full range of values found for clay, silty clay, clayey silt, interbedded sand/clay/silt with Method 1, a very reasonable result, again indicating very good agreement between the two approaches. Higher resistivity values were found for clay above the water table using Method 2-resistivity values as high as 31 m. In Method 2 we defined a category "mixed fine and coarse" below the water table, finding resistivity values that overlapped with those for clay and sands/gravels. We conclude that Method 2 provided us with resistivity values that agree well with those determined using the established approach of Method 1; and provided further information about the impact of saturation state (above or below the water table) on transforming resistivity values to lithology.
An important feature in the use of Method 2 is the ability to capture a distribution of resistivity values for each lithology. This allowed us to display the probability of the occurrence of each lithology ( Figure 9) and our "best guess" (Figures 10 and 11) that transforms each resistivity value to the lithology most likely to occur. Working with distributions and providing a lithologic interpretation in terms of probabilities of occurrence is a first step towards capturing and communicating to endusers uncertainty in the interpretation of AEM data. A recent paper used the dataset presented here to develop new ways to use color wheels to display this uncertainty (Nordin et al. 2016).
Let us now consider whether the resulting variation in lithology, derived from the AEM resistivity models, allowed us to map out the large-scale hydrostratigraphy. We compare our results to a cross section in the area of the flight lines, provided by the Kaweah Delta Water Conservation District (Fugro West 2007). A simplified version of the portion of the cross-section closest to our flight lines (B-B in the referenced report), labeled according to the interpretation in the report, is provided in the lower half of Figure 13; the location of this portion of the cross section is shown in Figure 2, labeled A-A .
In the upper half of Figure 13, is the interpreted section of Line 3 that covers the same depth range as the cross section (Line 3 extended to a depth of ∼550 m below the ground surface, while the cross-section extends to a depth of ∼400 m). We have overlain on the Line 3 section the interfaces and units marked on the cross section at the depths given in the cross section. As can be seen in the location map in Figure 2, the cross section and Line 3 are not coincident. They cross in one location (the 15 km mark on Line 3) and then diverge so as to be approximately 10 km apart at the start of Line 3. We attempted to use information from the report to determine the strike and dip of the various units and obtain a projection of the cross section onto the same plane as Line 3, but the well control and quality of the data did not warrant this level of analysis. In many places in the report the mapped interfaces are labeled with questions marks.
There are four key units shown in the cross section that we wanted to compare to our interpretation of Line 3: the Corcoran Clay, the upper and lower aquifer, and the impermeable unit at the base of the section. We found good correspondence between the AEM results and the cross section in mapping the average depth of the Corcoran Clay. The thickness of the clay was determined, through our AEM measurements and inversion, to be approximately 50 m in the southwest, pinching out in the northeast. This is thicker than the average value of approximately 20 m obtained in this area from interpolation of electrical logs. We attribute this difference to the inability of the AEM data to quickly recover to the high resistivity values at the base of the clay, as discussed above. We therefore conclude that while we can map the presence and approximate depth of the Corcoran Clay with the AEM method, we will tend to overestimate its thickness when it is thinner than 10% of the depth to the middle of the model layer.
The cross section and Fugro West report identify an upper and lower aquifer, differentiated based on oxidation state. Reddish coloring in the cuttings is noted in the drillers' logs from the region referred to as upper aquifer; this is presumably due to the presence of iron oxide. The AEM method cannot differentiate between oxidized and reduced sediments unless there are resistivity contrasts, so the interface between these two units is not seen in the interpreted sections. What we do see in all the lines, as shown in Line 3 in Figure 13, are the interlayered coarse and fine sediments that are described as characteristic of both the upper and lower aquifer.
It is important to note that the AEM method used in this study is capable of resolving packages of interlayered materials, but cannot resolve individual thin layers. There will be an averaging of resistivity values. To demonstrate this, we show in Figure 14 the result of a modeling exercise. We constructed a one-dimensional model of sand and clay layers of varying thickness, using our determined resistivity values to assign a resistivity value of 25 m to the sand, and a value of 12 m to the clay. We modeled the response of the SkyTEM system to these layers, then inverted the synthetic data to obtain a resistivity model, and transformed that to lithology using our established relationship between resistivity and lithology. The result shows the complicated nature of the link between what is present and what is captured in the AEM sounding. Given that the vertical resolution degrades with depth, the ability to resolve specific layers depends on the thickness and the depth. Closest to the surface, where the vertical resolution is the best, we recover a 3 m thick clay layer. At greater depths, when the layer is sufficiently thick, the lithology is accurately identified. There are many cases, however, where alternating thin layers of sand and gravel, and clay are identified as mixed fine and coarse due to the averaging of the resistivity values. The sections that are predominantly sand and gravel, with thin layers of clay, appear as sand and gravel, and mixed fine and coarse. The sections that are predominantly clay, with thin layers of sand and gravel, appear as clay and mixed fine and coarse. We are thus able to differentiate sections that are predominantly coarse-grained from those that are predominantly fine-grained, but cannot map in detail the fine structure of lithology variation.
A key feature in the cross section in Figure 13 is the so-called impermeable zone at depth, seen in the lower right corner, which defines the base of the lower aquifer. In the right half of Line 3, our interpretation shows thick clay at depth, defining the base of the lower aquifer. While the transition to clay in the interpreted SkyTEM section, occurs below the top of the impermeable zone in the cross section, we note that the Fugro West report indicates that the top of this zone is not well constrained. There are only three electrical logs, and no drillers' logs reaching this depth. The electrical logs are of questionable quality and do not show a clearly interpretable transition from high resistivity to low resistivity. The top of the impermeable zone could easily be more than 50 m higher or lower.
We note that the impermeable zone in the cross section is shown as extending across the full length of the section, while we show no thick clay at depth in the left section of Line 3. The cross section in Figure 13 covers a more limited depth range than the AEM data, due to the absence of wells at greater depths. There are in fact no wells that support the presence of the impermeable layer at depth in the left half of the cross section; the extension of the impermeable layer is simply an interpretation. Seen in Figure 8a, in the lower left corner of the complete depth-section for Line 3, is a region of high electrical resistivity. Our interpretation of this high resistivity unit suggests a large area of sand and gravel, as shown in Figure 10.
The fence diagram in Figure 11, with the interpretation of all the lines of AEM data, provides a more complete view of the subsurface. Throughout the area we see the mix of lithologic units in the top 150 to 200 m. There are very distinct differences at greater depths, with the lower approximately 300 m in the northeastern half of the area imaged as a thick clay unit which is completely absent in the southwestern half. There we see a  nearly continuous layer of highly resistive materials interpreted as sand and gravel. While there are no lithology logs that reach this depth that could be used to support this interpretation, the recent drilling of deep wells in this area has reported finding sand and gravel (A. Fukuda, personal communication, 2016). The ability to map the hydrostratigraphy, below the current depth range of wells in an area, is one obvious benefit of utilizing the AEM method to acquire the data needed for improved conceptual models of aquifer systems. Ongoing research is now focused on integrating all well data with the AEM data to build a conceptual model for this area, and evaluate the extent to which incorporating the AEM data enhances the accuracy and usefulness of the groundwater models used to support local groundwater management decisions.
The motivation for conducting this study was to assess the value of the AEM method as a cost-effective way to map the aquifer systems of the Central Valley. In terms of mapping the aquifer systems, wells are an established method and can provide detailed information, at one location. While the vertical resolution of AEM data can never match that of a well, even abundant well data yield little information in the horizontal directions, where difficult to identify features such as windows through aquitards and coarse-grained channel deposits can radically change groundwater flow and transport velocities. There is considerable value in mapping out the large-scale hydrostratigraphy and the under-sampled, lateral variations in aquifer system properties with the AEM method. In addition, the many wells that have been drilled in the Central Valley tend to be shallow, so do not provide information about the deeper parts of the aquifer system. It is important to acknowledge, however, that the presence of wells is a necessary part of the analysis and interpretation of AEM data; but even with sparse well data, as was the case in this study, valuable information about the subsurface hydrostratigraphy can be obtained from the AEM data.
In terms of the costs of acquiring subsurface data, wells can be expensive to drill, especially deeper wells, with the cost of a monitoring well typically exceeding $100,000. A conservative estimate of the cost of an AEM survey, for data acquisition, analysis and interpretation is $450 per line km, with an additional approximately $20,000 for mobilization and demobilization of the system. A ground-based geophysical method could be employed at a few locations for less than the cost of mobilization of the AEM system. But the coverage that can be obtained with the AEM method could not be achieved with any other method for similar costs. Given the cost of drilling wells, especially monitoring wells, it could therefore be highly cost-effective to use AEM data to obtain a large-scale image of an area, and then use the results to guide the selection of drilling locations. The complementary use of well data and geophysical data, where each has its costs and benefits, is an effective approach to obtaining the information required to support the management of groundwater resources.

Conclusions
This study has allowed us to assess the viability of using the AEM method to map the aquifer systems of the Central Valley; and, more generally, to map similar aquifer systems in sedimentary basins containing significant fines. We conclude that the regional implementation of the AEM method could provide critical information about the hydrostratigraphy of the aquifer systems needed for groundwater management. What we were able to derive from the AEM data about the large-scale structure, far exceeds what is possible using traditional methods based on the drilling of wells.
We found it possible to image the electrical resistivity to a depth of approximately 500 m; given the geology of the Central Valley we would expect to find similar imaging depths at most other locations throughout the valley. This covers the relevant depth range for current groundwater management and provides information about the deeper regions of the aquifer system not currently sampled by wells. In addition, the lateral spatial resolution seen in the AEM data could never be obtained with well data, and is needed to reveal the large-scale heterogeneity that should be captured in groundwater models.
We developed a new methodology to transform the resistivity model to lithology which can be applied throughout the Central Valley, and could be widely adopted as a new approach to the interpretation of AEM data. The key limitation in this approach will always be the shallow depth of most of the lithology logs, resulting in a resistivity-lithology relationship that is established and thus valid at shallow depths, but only assumed to be valid at greater depths. We are currently exploring ways to correct for the effect of depth on this relationship. The approach yields a distribution of resistivity values, positioning us to be able to quantify and account for this source of uncertainty in using the AEM data to generate conceptual models. Our resulting interpretation of lithology was consistent with other lithologic information from the area, while noting that the vertical resolution of the AEM method resulted in an overestimation of the thickness of the Corcoran Clay unit and an inability to detect relatively thin layers at depth.
AEM imaging to depths of approximately 500 m can provide critical information about the distribution and connectivity of hydrostratigraphic packages that would support the development of conceptual models and would reveal permeable pathways that could be used to recharge groundwater at shallow and deeper levels. The acquired information would directly support specific management actions such as the selection of sites for surface spreading recharge, and the siting of monitoring wells. In this study we used the SkyTEM 508 system as we were interested in maximizing the depth of imaging while maintaining reasonably high resolution in the top approximately 100 m. If the goal in a project were to more accurately resolve the top 50 to 100 m to assess recharge potential, other AEM systems could be used.
This study was motivated by our interest in exploring the use of the AEM method to address the critical need for subsurface data, in order to implement new groundwater legislation in California. We conclude that the acquisition of AEM data throughout the Central Valley would provide the hydrogeological framework needed to support the establishment of sustainable groundwater management in California.