Improved Imaging of the Large-Scale Structure of a Groundwater System With Airborne Electromagnetic Data

The but Abstract Working with airborne electromagnetic (AEM) data acquired in the Kaweah Subbasin in the Central Valley of California, USA, we developed a new approach for imaging the bedrock surface and the confining Corcoran Clay layer. Our approach, which incorporated the prior knowledge of the targets for the improved imaging, included multiple L 2 -norm and L p -norm inversions as well as an interpolation process. The major improvement in imaging the targets was made in the L p -norm inversion step by incorporating prior knowledge. For the Corcoran Clay, pairs of resistivity and driller's logs at two wells guided the selection of the primary resistivity model and were used to increase the accuracy of the estimated Corcoran Clay thickness. The bedrock surface was poorly constrained by well data in the existing groundwater model, appearing as a flat surface. Acquired AEM data covered most of the area, so we had higher confidence in the obtained bedrock surface at depths ranging from 160 There was relatively good agreement between the location of Corcoran Clay AEM

2 of 25 can be highly variable and challenging to quantify. In addition, the spatial density of the driller's logs might be too low to adequately capture the continuity of the large-scale subsurface features. Uncertainty in representing the large-scale structure can have a significant impact on the predictions obtained from a groundwater model (Nilsson et al., 2007;Refsgaard et al., 2012).
The airborne electromagnetic (AEM) method is potentially a highly effective means of imaging large-scale structure for the development of a groundwater model (Barfod et al., 2018;Foged et al., 2013;He et al., 2014;Kang et al., 2021;Knight et al., 2018;Sattel & Kgotlhang, 2004;Wynn, 2002). Inversion of the acquired AEM data recovers a resistivity model of the subsurface. When well data are available proximate to the AEM data, this resistivity model can be transformed to a model displaying the variation in lithology or sediment type. This observed variation can then be interpreted to identify the large-scale features of the groundwater system, for example, mapping out the lithologic units or sediment packages that define the major aquifers and aquitards. One challenge when using the AEM method for this application, however, is the limited spatial resolution of the method resulting in uncertainty in the recovered resistivity model and interpreted locations of the large-scale features.
For the improved imaging of the large-scale features, in this study, we developed a new approach to the AEM inversion workflow for the improved imaging of the large-scale structure. This approach was designed for the scenario where AEM data are being acquired in an area with no existing groundwater model and limited high-quality well data. Our specific objective was to achieve improved imaging of the subsurface so as to capture, as accurately as possible, the large-scale features required for developing a groundwater model primarily using the AEM data; but high-quality well data from a few locations were needed as additional data to guide the inversion process. We accomplished this through a "targeted inversion approach" where we defined as targets specific features which we wanted to accurately resolve in our study area-the bedrock surface (at the top of the bedrock) underlying the sediments of the aquifer system and the confining clay layer dividing the system into an upper and lower aquifer. While these two features were defined for our study area, the Kaweah Subbasin located in California's Central Valley, USA, they define part of the large-scale structure in many mountain-valley aquifer systems throughout the Central Valley and elsewhere in the world. Using the commonly adopted L 2 -norm approach to the inversion of AEM data Siemon et al., 2009;Viezzoli et al., 2008), a bedrock surface will appear as a smooth transition rather than a distinct boundary in the recovered resistivity model, introducing uncertainty in its location. Similarly, the lack of a distinct boundary at the upper and lower boundaries of a clay layer will introduce uncertainty in its location and can cause the thickness of the layer to be significantly overestimated. Errors in the depth to bedrock and the thickness of the confining layer, due to the uncertainty in the AEM inversion process, could significantly impact the accuracy of the results obtained from a groundwater model with boundaries derived from the AEM data. Further, these errors are accelerated by the limited resolution of the AEM data. Alternate approaches to the inversion of AEM data, such as the sharp inversion, which constrains the inversion to image sharp boundaries , could improve the accuracy in locating the bedrock surface. However, this approach would still fail to locate the depth and thickness of the confining layer given that the AEM data cannot independently resolve the thickness and resistivity of the layer (Geowissenschaften et al., 2007). Therefore, a new approach is needed to incorporate additional information about either the resistivity or thickness of the layer for accurate imaging of the confining layer.
Our targeted inversion approach utilized a multi-step inversion process which included both the standard L 2 -norm inversion and the L p -norm inversion (Fournier & Oldenburg, 2019), the latter allowing us to incorporate prior knowledge about the targets. To form 2D maps displaying the locations of the targets-the large-scale features, we integrated high-quality well data and the resistivity models recovered through L 2 -norm and L p -norm inversions. To demonstrate the efficacy of our approach, the locations of the targets that we obtained were compared with the locations within the existing groundwater model which had been constructed using available well data. Once the large-scale structure had been identified, this was used in a final inversion to obtain improved imaging of the smaller-scale structure. We use the term "smaller-scale" as defined by Nilsson et al. (2007), to refer to the scale at which there is spatial variation in the hydrogeology within a large-scale unit. Through the developed multi-step approach, we formalized a process, readily transferable to other locations, designed to obtain improved imaging of groundwater systems from AEM data. KANG ET AL.

Hydrogeology and Existing Groundwater Model in the Study Area
The Kaweah Subbasin is located in the southern part of the Central Valley, referred to as the San Joaquin Valley, with the eastern edge defined as the interface between the valley floor and the foothills of the Sierra Nevada Mountains. As part of their efforts to better understand their groundwater resources and meet the requirements of California's Sustainable Groundwater Management Act (SGMA, 2014), the local water agencies developed a groundwater model of the subbasin using 600 driller's logs and 50 oil and gas well logs, the latter including both lithology logs and resistivity logs (Fugro West, 2016). The model covers the subbasin, extending ∼2 km outside the subbasin boundary except along the eastern edge, ranges from a depth of ∼120 m in the east to ∼500 m in the west and contains three layers, of variable thickness, with a cell size of 150 m × 150 m. The groundwater model boundary is shown in Figure 1a.
A description of the hydrogeology of the subbasin is given in the report by Fugro West (2016). The key features defining the layers in the groundwater model are in a geologic cross-section (modified after Fugro West, 2016), the location of which is shown in Figure 1a with the cross-section in Figure 1b. The boundaries of the layers were fixed in advance and not adjusted during the groundwater modeling process. The sediments within the subbasin are divided into an upper and lower aquifer in the western half by a confining clay layer, the Corcoran Clay. The aquifers are described as Quaternary alluvium with numerous interbedded sands and clays. Generally, more clays are present in the lower aquifer than in the upper aquifer. The Corcoran Clay is described as Quaternary lacustrine and marsh deposits. The aquifers are underlain by bedrock in the east which abuts low permeability sediments to the west along a boundary interpreted to be the Rocky Hill Fault.
The groundwater model contains three layers: The upper aquifer which ranges from 100 to 150 m in thickness; the Corcoran Clay which ranges from 10 to 20 m in thickness, and the lower aquifer which ranges from 200 to 350 m in thickness. Note that the three layers are defined throughout the model. In the current groundwater model the extent of the Corcoran Clay was modified from the work of Page (1986) using available well data (Fugro West, 2016) and corresponds to the location where the Corcoran Clay is no longer continuous. Where the Corcoran Clay is absent, the three layers effectively become a single aquifer. As shown in Figure 1b, the base of the model in the east is defined by the top of the bedrock, then going west is interpreted to lie at the top of the package of low permeability sediments.
Both the top of the bedrock and the location Corcoran Clay are important large-scale features of the groundwater model, the former defining the position of a no-flux boundary and the latter controlling the hydraulic connectivity between the upper and lower aquifers; thus, the accuracy of their locations and the thickness of the Corcoran Clay directly impacts the accuracy of model predictions. The spatial coverage of the well logs used to constrain the depth to bedrock was poor (GKGSA, 2020), due to the fact that there were relatively few wells reaching the bedrock surface. While the spatial coverage of the well data used to map the location of the Corcoran Clay was reasonably good, there are concerns about the quality of the well data (e.g., a large locational error [∼2 km] for the driller's logs) and interpretations of continuity between wells, particularly at the edges of the Corcoran Clay where it becomes discontinuous.

The AEM Method
The AEM method uses electromagnetic induction to obtain information about the resistivity of the subsurface. Time-varying electric currents are injected through a transmitter loop suspended below an aircraft to generate induced currents in the subsurface. These induced currents will depend upon the resistivity of the subsurface and generate an induced voltage that can be measured at a receiver loop also carried by the aircraft. The aircraft moves continuously while the receiver loop is recording, with the raw voltages stacked at predetermined intervals to provide measurements referred to as AEM soundings or the observed AEM data.
Inversion of the observed AEM data is necessary to obtain a 3D resistivity model of the subsurface resistivity. While there are many approaches to modeling AEM data (Brodie & Sambridge, 2009;Farquharson et al., 2003;Yang & Oldenburg, 2012), in hydrogeologic applications when it can be assumed that the structure is layered, a spatially constrained inversion approach (Viezzoli et al., 2008) is most commonly used to recover, at each AEM sounding location, a vertical 1-D profile of resistivity fitting the observed AEM data. The spatial constraints, KANG ET AL.
10.1029/2021WR031439 4 of 25 favoring a smooth transition of resistivity values between adjacent sounding locations, are implemented in the inversion through a regularization function. The regularization function also includes a reference model through a smallness constraint (Oldenburg & Li, 2005). Conventionally, an L 2 -norm is used for the constraints in the regularization function (Oldenburg & Li, 2005;Tarantola & Valette, 1982;Tikhonov & Arsenin, 1977). In this study, in addition to an L 2 -norm we utilized an L p -norm which makes it possible to incorporate additional prior knowledge (Fournier & Oldenburg, 2019;Vignoli et al., 2015).
The AEM inversion is non-unique, that is, there are many resistivity models that can fit the data. A typical form of non-uniqueness is the presence of a depth, commonly referred to as the depth-of-investigation (DOI), below which resistivity changes in the model do not make a noticeable difference in the AEM response. We used Oldenburg and Li's (1999) approach performing two inversions with different homogeneous reference models to define the DOI. Another well-known form of non-uniqueness, of specific relevance to our study, is due to the fact that the AEM response is sensitive to the conductance of a layer, which is the thickness divided by the resistivity (Geowissenschaften et al., 2007). While in rare circumstances both the resistivity and thickness of a layer can be independently resolved, generally this is not possible. Given this limitation, other reliable information needs to be incorporated into the inversion.

Available Data
In 2018, ∼800 km of AEM data were collected over the Kaweah Subbasin using the SkyTEM 312 system (Sorensen & Auken, 2004). Figure 1a shows the AEM soundings along the flight lines of the survey, which were located to maximize coverage in the region of the groundwater model while avoiding urban areas (Kang et al., 2020). The acquisition of the AEM data was managed by, and the data processed by, Aqua Geo Frameworks (Asch et al., 2019).
Our objective was to develop an approach that would make it possible to accurately image the bedrock surface and confining clay layer using primarily the AEM data. In addition to the AEM data, we utilized data from six wells that could be considered, with confidence, to be of high quality. This included three driller's logs in the northeastern corner of the subbasin (locations shown in Figure 1a) that had been manually inspected (Steklova, personal communication, 2020). Criteria used to identify the high-quality well include the knowing the well location to within 100 m, the depth of the well where deeper is better, the presence of geophysical logging data, and detailed descriptions of the sediments. These provided point information about the depth to bedrock. We also had access to an additional three pairs of logs, resistivity (16-inch normal) and driller's logs from the same monitoring wells, referred to as Wells A, B, and C (locations shown in Figure 1a), considered by the local water agencies to be of very high quality. These pairs of logs, in the western part of the subbasin, provided information about the location of the Corcoran Clay. We took each pair of logs and, at the well location, defined the depth and thickness of the Corcoran Clay.

Methodology
The workflow for our targeted inversion is illustrated in the top section of Figure 2. The targeted inversion approach includes three steps; (a) L 2 -norm inversions to obtain a starting resistivity model, a domain of interest, a target-containing layer from AEM data, and a DOI; (b) L p -norm inversions to obtain the locations of the targets along with error, at each AEM sounding location, and (c) the interpolation process integrating the locations of the targets from the AEM data and the high-quality well data to obtain 2D maps displaying the locations of the two targets-the top of the bedrock and the Corcoran Clay layer.
In addition, we developed a subsequent approach that can extract additional information about the smaller-scale features of the groundwater system from the AEM data. The output of the targeted inversion approach-the location of large-scale structure-was used to constrain the inversion, and therefore it was referred to as structurally constrained inversion; a workflow of this approach is shown in a bottom panel of Figure 2. Applications to date have been in the inversion of ground-based or marine-based EM data and have used the large-scale structure interpreted from seismic sections available in the same area as the EM data (Brown et al., 2012;Key, 2009).
For all inversions of the AEM data, we used the regularization function, ( ), which can be written as 6 of 25 where the first term is the smallness constraint, and the two other terms are the spatial constraints, one acting in the radial direction, denoted by , and the other in the vertical direction, denoted by ; and ref indicate an unknown, to be recovered, resistivity model and a reference model, respectively (Cockett et al., 2015a;Oldenburg & Li, 2005). The unknown and reference resistivity models were defined at each AEM sounding location with resistivity values assigned in 39 vertical cells; the thickness of a cell at the ground surface was 3 m and increased at a constant rate resulting in a total depth of 550 m. Alpha values, , , , determine the relative importance of each term. Given that layered sediments are the dominant materials in the study area, was set to five times smaller than to indicate our preference for laterally continuous rather than vertically continuous structure; was fixed to 1. The value of was set to either include ( = 1) or not include ( = 0) the smallness constraint in the regularization function. Further details about these parameters can be found in Oldenburg and Li (2005).
The exponent acts on all the constraints. In terms of the smallness constraint, p influences the form of the resulting distribution of all the resistivity values in the 3D resistivity model. The conventional L 2 -norm inversion, which corresponds to setting = 2, favors a smooth Gaussian distribution, where the center of the distribution is the reference model. In the L p -norm the value of can vary from 0 to 2 (Fournier & Oldenburg, 2019), which has the impact of favoring distributions that vary from sparse (isolated peaks) at low values of , to Gaussian at p = 2. The extreme value of = 0 denotes a distribution that has the fewest possible number of peaks deviating from the reference model.
The value of also acts on the spatial constraints. Setting = 2 denotes a smooth transition in resistivity values in space (laterally and vertically). The extreme value of = 0 indicates a model with an abrupt transition − a large resistivity contrast, at the fewest number of interfaces between cells possible while still fitting the data. Assuming that the greatest resistivity contrast will be found at the interfaces (lateral and vertical) between hydrogeologic units, not within the units, this allows an abrupt transition only at an interface between subsurface units.
While the model space for the solution with the L 2 -norm is convex, making it relatively straightforward to solve the optimization problem, when ≤ 1 the model space for the L p -norm is non-convex making it challenging to find a solution. The L p -norm has been implemented for the spatial constraints in the inversion of AEM data (often referred to as a sharp inversion; Vignoli et al., 2015). In this study, we implemented the L p -norm for both the spatial and smallness constraints for the inversion of AEM data. To the best of our knowledge, this had not been previously done, and introduced an even greater challenge to finding a solution given the increased level of non-convexity, but provided the enhanced flexibility we desired for incorporating prior knowledge. A similar challenge was identified while implementing the L p -norm for inverting other types of geophysical data Portniaguine & Zhdanov, 1999;Vignoli et al., 2012). We used the strategy developed by Fournier and Oldenburg (2019) for the inversion of potential field data, which first finds an L 2 -norm solution then activates the L p -norm inversion.
The cell-based and face-based (where "face" refers to interface between the cells) weightings ( cell and face , respectively) capture the level of confidence in the smallness constraint and spatial constraints at a given cell or face; the greater the weighting, the higher the level of confidence. For each face of the resistivity model, a different value of the face-based weighting can be set; similarly, for each cell of the resistivity model a different value of the cell-based weighting (e.g., cell ) can be set. These weightings were used for the structurally constrained inversion; for other AEM inversions all these weightings were set to 1. The faced-based weightings were used to allow for a sharp resistivity contrast at the target. Assigning a smaller weighting to selected vertical or lateral faces allows the inversion to put a sharp resistivity contrast at the face. The cell-based weightings for the spatial constraints: cell , cell were used to minimize the resistivity variations in the vertical and lateral directions for cells corresponding to either the bedrock or the Corcoran Clay layer compared to the other cells. Assigning larger weightings to selected cells than to other cells allows smaller resistivity variations. The cell-based weighting for the smallness constraint was used to provide a higher level of confidence to selected cells of a reference model by assigning a larger weighting than that assigned to other cells.
To conduct our inversions we used SimPEG, a Python-based open-source geophysics software package (Cockett et al., 2015b;Heagy et al., 2017). Further details regarding the solution of the inverse problem can be found in Appendix A.

Step 1: L 2 -Norm Inversions
The parameters used in all inversions in our multi-step approach are given in Table 1. In Step 1, we carried out three L 2 -norm inversions to obtain a starting resistivity model, a domain of interest, a target-containing layer, and a DOI ( Figure 2). A domain of interest was obtained for each target; one for the top of the bedrock, bedrock , and the other for the Corcoran Clay layer, clay . These domains defined the general regions within which we 8 of 25 identified the presence of each target. Given that our purpose was to identify zones in which targets are present, an L 2 -norm inversion favoring a smooth transition of resistivity values in the lateral directions was most appropriate for this purpose, as it would limit the number of smaller-scale features included in the recovered resistivity model, making it easier to identify the large-scale features. We only used spatial constraints, setting = 0 in the first inversion (Inversion one in Table 1) so did not need to define a reference model. This setup (or setting a very small value of ) is considered as standard practice for the inversion of AEM data; a recovered resistivity model from this inversion was selected as a starting resistivity model. bedrock and clay were selected from the L 2 -norm inversion by visually identifying the boundaries of the targets (i.e., where there appeared to be the largest change in resistivity in the vicinity where the target was expected). We then added a buffer of a few kilometers to these boundaries to ensure that we captured the full extents of the targets within the defined domains. This manual approach was selected as there was no obvious numerical way (e.g., Caterina et al., 2013) to identify these boundaries, given the many other changes in resistivity within the models.
While the first L 2 -norm inversion was not sufficiently accurate for locating the Corcoran Clay, we used it to identify a depth interval that would include the Corcoran Clay that was used in developing reference models for the L p -norm inversions in the next step; this depth interval was referred to as a target-containing interval in the workflow shown in Figure 2. A target-containing interval was not needed for the bedrock surface. Identifying the top and base of the Corcoran Clay was done subjectively, by selecting the top and base of a continuous low-resistivity zone in the 3D resistivity model. Note that this "Corcoran-Clay-containing" depth interval was on the order of 80-160 m in thickness, so much larger than was expected for the thickness of the Corcoran Clay alone, based on driller's logs.
Two other L 2 -norm inversions (Inversions 2 and three in Table 1) were used, with the smallness constraint added ( = 1), to calculate the DOI (Oldenburg & Li, 1999). Adding the smallness constraint with a relatively high value regularizes the inversion to be close to the reference model. This impact of the smallness constraint increases with depth due to decreasing sensitivity of the AEM data. Oldenburg and Li used this concept to estimate a depth at which AEM data is insensitive to subsurface resistivity variations. Two separate inversions with two homogeneous reference models were required for this estimation; the two homogenous resistivity values used in this study were 17 Ωm and 25 Ωm .

Step 2: L p -Norm Inversions
The L 2 -norm inversions allowed us to identify the regions where our targets were present, but there was still considerable uncertainty in the exact locations of these targets. In this next step, in order to improve our ability to locate these targets, we incorporated prior knowledge specific to each target and iteratively performed L p -norm inversions within the domain of interest for each target until two resistivity models-a primary model (the better) and a secondary model-showing reasonable match with the high-quality well data were obtained. The locations Step 1: L 2 -norm inversion 2 (2, 2, 2) 17 Step 2: L p -norm inversion 9 of 25 of each target were obtained from the primary model; both primary and secondary models were used to provide a measure of error for the purpose of weighting in the next interpolation step.
Our prior knowledge was based on a general understanding of the hydrogeology of the study area (Faunt, 2010;GKGSA, 2020) and familiarly with the resistivity of geologic materials. In terms of locating the top of the bedrock, we used the fact that there is a large contrast between the resistivity of saturated sands and clays which make up the aquifer and the resistivity of granitic bedrock. In locating the Corcoran Clay, we used the fact that there is a large contrast between the resistivity of the coarser-grained sediments and clay. Well data also indicated that the thickness of the Corcoran Clay layer was generally much thinner than the over/underlying aquifers. Below we describe the methodology used to resolve each of the targets.

Top of the Bedrock
To locate the top of the bedrock, we performed two L p -norm inversions, setting = 0 (only the spatial constraints were used). Selecting the -values of the spatial constraints ( , ) was an iterative process, where a number of inversions were run with variable -values ranging from 0 to 2. We found all the resulting resistivity models showed good agreement with the driller's logs. From these models, we selected two end members; a) a primary model, where the top of the bedrock was easily identified due to the sharp resistivity contrast at the bedrock surface, and b) a secondary model, where the top of the bedrock was difficult to identify due to the smoothly varying resistivity values across this surface. The primary resistivity model used = 0 and = 0 and the secondary resistivity model used = 2 and = 2 (i.e., L 2 -norm; Inversions 4 and five in Table 1, respectively). We note that setting values of and to zero for the L p -norm inversion tailored the inversion to assign a large resistivity contrast only when required to fit the data.
In the primary resistivity model, we determined, at each sounding location, whether the bedrock surface was present or absent and, if present, estimated the depth to bedrock. At sounding locations where the bedrock surface was present, we also estimated the depth to bedrock from the secondary resistivity model. We assigned the top of the bedrock to the depth that divided the full depth range into two depth intervals-one above and one below, with minimal variance in the resistivity values of each depth interval. We explored alternatively using the minimal variance in the inverse of resistivity, or the log of resistivity, but found these approaches resulted in the same or very similar (± one vertical cell) bedrock surfaces.

Corcoran Clay Layer
Two L p -norm inversions were conducted to locate the Corcoran Clay so as to have a primary resistivity model and a secondary model. We defined the primary resistivity model as one that showed good agreement with the resistivity and driller's logs from Wells A and B; and the secondary resistivity model as one that showed good agreement with the resistivity and driller's logs from Well B. Both L p -norm inversions used the same spatial constraints but a different smallness constraint by changing the reference model.
In setting the spatial constraints, we used the prior knowledge that there should be a large resistivity contrast at the boundaries of the Corcoran Clay so, as was done in locating the top of the bedrock, -values for the spatial constraints were set to zero. In addition, we incorporated the fact that the Corcoran Clay is much thinner than either the overlying upper aquifer or underlying low aquifer. In a recovered resistivity model therefore, most of resistivity values should correspond to the aquifer materials and relatively few to the Corcoran Clay. To obtain such a sparse distribution of resistivity values, centered on the aquifer resistivity values, we added the smallness constraint to the regularization with s = 0. We note that these aquifer resistivity values are unknown, and must be estimated.
To define the two reference models used to obtain the primary and secondary models, we identified two regions: outside and inside the Corcoran-Clay-containing depth interval identified using the starting resistivity model in Step 1. In both reference models we set the resistivity values for the cells outside the Corcoran-Clay-containing layer equal to those in the starting resistivity model (from Step 1), but varied the resistivity value of the cells inside the Corcoran-Clay-containing layer. We allowed the homogeneous resistivity of the Corcoran-Clay-containing layer, clay ref , to range from 10 Ωm to 30 Ωm , iteratively solving for two resistivity values which, when used in the reference models, resulted in the recovery of the primary and secondary resistivity models. Given that the vertical extent of the Corcoran-Clay-containing layer likely included borders of the surrounding aquifers, the range of resistivity was determined by choosing 1st and 99th percentile of the resistivity values within the 10 of 25 Corcoran-Clay-containing layer. Resulting values of clay ref for the primary and secondary models were 30 Ωm and 20 Ωm , respectively (Inversions 6 and 7 in Table 1, respectively). In the two recovered resistivity models (primary and secondary), we determined at each sounding location whether the spatially continuous Corcoran Clay was present or absent, interpreting a large resistivity contrast at two interfaces as the top and base of the Corcoran Clay. The depth to the top interface was defined as the depth to the Corcoran Clay. For the Corcoran Clay thickness, we could have used the distance between the top and base, but wanted improved accuracy given the inability of the AEM method to independently resolve resistivity and thickness. Representing the Corcoran Clay, with total thickness cc and resistivity cc , as a series of layers, each of thickness with resistivity , what is captured in the AEM data is an integrated measurement that can be represented as applying an electric field oriented parallel to these layers resulting in the following relationship: At an AEM sounding, the Corcoran Clay, with total layer thickness AEM and total resistivity AEM , was recovered in the resistivity model as a layer of two or three resistivity cells each with thickness c,i and resistivity c,i where We wanted to accurately recover the true cc as opposed to AEM . We accomplished this through a calibration process using well data, referring to the determined thickness of the Corcoran Clay as calibrated and the determined resistivity of the Corcoran Clay as calibration .
where superscript A indicates the closest sounding location to Well A. Using the two soundings closest to Well A and the two closest to Well B yields: Rearranging this, we obtain The calibration process, resulting in a calibrated Corcoran Clay thickness, calibrated can be written as This calibration was applied to all values of Corcoran Clay thickness ( AEM ) obtained from the primary resistivity model to obtain our best estimates of Corcoran Clay thickness. The calibration resistivity of the primary model was 5 Ωm . The same procedure was applied for the secondary model, but the Corcoran Clay was only found to be present at the sounding location closest to Well B, not at the sounding location closest to Well A. Therefore, only Well B was used in the calibration, so the calibration resistivity became: calibration = wellB cc ∑ c,i c,i . The resulting calibration resistivity for the secondary model was 7 Ωm .

Step 3: Interpolation
The locations of the targets were obtained from Step 2 at each sounding location. In this step, the locations of the targets from the AEM data were integrated with the high-quality data from six wells to generate 2D maps of the targets on a uniform grid with a cell size of 150 m, which is the same as the groundwater model ( Figure 2). The accuracy of the depth to the bedrock surface and to the Corcoran Clay layer that we were able to obtain from AEM data was limited by the thickness of the resistivity cells. Therefore, we used the thickness of the resistivity cell at the target depth to provide a measure of error in the 2D maps. Due to the calibration process, the estimated Corcoran Clay thickness was not limited by the thickness of the resistivity cell. The fact that the calibration resistivity values were 5 Ωm and 7 Ωm indicated that there could be roughly 30% relative error in the estimate of the Corcoran Clay thickness.

Top of the Bedrock
The top of the bedrock was displayed in a 2D map as the depth to bedrock. For developing this map, the input data were: (a) depth to bedrock from the three driller's logs, (b) depth to bedrock estimates at all sounding locations in the primary and secondary resistivity models within bedrock where the bedrock surface was determined to be present in Step 2.
The depth to bedrock values from the primary model, aem target , were interpolated giving higher weighting to the values with lower errors; with error defined as the absolute difference between the estimates from the primary and secondary models. The depths extracted from the three driller's logs, well , were used as hard constraints. Given the degradation in the lateral resolution of AEM data with depth, an estimate of the depth to bedrock obtained at a sounding location was assigned to an area, the radius of which depended upon the depth, target . The definition of this radius is based upon the propagation of induced currents in a homogeneous subsurface (Nabighian, 1979), and is given by the following equation: where tx loop is the effective radius of the transmitter loop; tx loop was 11 m. Assigning depth values to the 2D map was done through an averaging function, avg [⋅] , which averages grid values of the 2D map, 2 , within the distance aem from each sounding location: The interpolation is an inverse operation of Equation 8, which can be written as where aem target indicates the measure of error. The solution of this inverse problem adopted the framework used for the inversion of the AEM data (Appendix A). Given that this inverse problem was non-linear due to the imposed hard constraints, providing a starting 2D map of the depth to bedrock was required. For this, we used an inverse distance weighting of aem target to populate all cells in the 2D grid. To incorporate our preference that the depth to the bedrock varies smoothly in the lateral direction, we only used the spatial constraints for a regularization function, which can be written as Both alpha values for the spatial constraints were set to 1, providing an equal weighting in the -and -directions.

Corcoran Clay Layer
The location of the Corcoran Clay was displayed in two 2D maps as the depth to the top of the Corcoran Clay and the Corcoran Clay thickness. As was done for mapping the depth to bedrock, we interpolated the weighted-by-error locations of the Corcoran Clay layer from the primary resistivity model, using the secondary model to estimate error. We included values for the depth to the Corcoran Clay only at sounding locations where the 12 of 25 Corcoran Clay was identified in the primary resistivity model from Step 2. We included values for the Corcoran Clay thickness at all sounding locations within clay with the thickness given as 0 at any sounding location where the Corcoran Clay was not identified. The three pairs of well data were used as hard constraints. The starting 2D map was generated by using the inverse distance weighting of the locations of the depth to the Corcoran Clay and Corcoran Clay thickness from the primary resistivity model.

Structurally-Constrained Inversion
After the locations of the targets were obtained from the targeted inversion approach, we applied the structurally constrained inversion (Figure 2) to recover the final resistivity model accurately imaging both the large-scale and smaller-scale structures. For this inversion (Inversion 8 in Table 1), we used the faced-based weightings, the cell-based weightings, and the reference model in the regularization function of the AEM inversion (Equation 1).
To allow for sharp resistivity contrast at the targets, we first found the closest vertical or lateral faces for each target and then assigned zero for the corresponding face-based weightings while the face-based weightings for other faces were set to 1. To minimize resistivity variations for cells corresponding to either the bedrock or the Corcoran Clay compared to the other cells, we assigned a large value (10) for the cell-based weightings corresponding to the spatial constraints (i.e., cell and cell ) for cells below the top of the bedrock and within the Corcoran Clay layer, while for other cells the cell-based weightings were set to 1.
Due to the large sensitivity of AEM data to conductive materials, it was easier for the inversion to update resistivity values around the Corcoran Clay layer, which altered the boundaries of the Corcoran Clay layer that we incorporated through the face-based weightings. Thus, it was necessary to provide additional constraint related to the resistivity information about the Corcoran Clay to preserve the boundaries from our targeted inversion approach. For this, we utilized the smallness constraint in the regularization function (Equation 1); this includes a reference model and a cell-based weighting. We generated a reference model composed of two regions: the Corcoran Clay layer and others. Due to the fixed vertical thickness of resistivity cells, it was not possible to directly use the depth to the Corcoran Clay and Corcoran Clay thickness obtained from our approach. The top and base of the Corcoran Clay were set to the closest faces resulting in small errors (∼2-5 m). Due to these errors, it was necessary to apply the same calibration idea (used in Step 2) when assigning resistivity values to cells corresponding to the Corcoran Clay layer. From the conductance of the Corcoran Clay layer obtained from the primary resistivity model ( Step 2), we calculated an equivalent resistivity value providing the same thickness of the Corcoran Clay at each sounding location. For other cells, we assigned a mean resistivity value, 16 Ωm , of the starting resistivity model (Inversion 1 in Table 1). The cell-based weighting for the smallness constraint was set to 1 for cells corresponding to the Corcoran Clay layer while a much smaller value (10 −3 ) was assigned for the other cells to minimize the impact of the reference model at those cells.

Results
The DOI, obtained in Step 1 using an L 2 -norm inversion, ranged from 240 to 380 m which covers the depth range within which both of the targets are present in the groundwater model. The domains of interest for the targets, obtained in Step 1 using an L 2 -norm inversion, are shown in Figure 3. As expected, bedrock was located along the eastern edge of the subbasin and clay in the western half. In Figure 3 we also show the locations within the domain of interest for each target where the target was determined, from the L p -norm inversion, to be absent; in the figure we show each tenth sounding (∼200-300 m spacing). As seen in the figure, there was a large AEM data gap at the center of the survey area due to the presence of urban areas. This resulted in a high level of uncertainty in using AEM data to map the extent of the targets in these areas. Shown in Figure 3 are the western boundary of the bedrock, and the northern and eastern boundaries of the Corcoran Clay from the groundwater model. The comparison with the results obtained from the AEM data will be discussed in a later section.
In Figures 4 and 5 we show, along two selected transects B-B' and C-C' (locations in Figure 3), vertical sections displaying the resistivity models recovered from the L 2 -norm and L p -norm inversions, respectively. All visualizations of recovered resistivity models are shown between the ground surface and the average DOI of 300 m. The section in Figure 4 focuses on the top of the bedrock, and Figure 4 on the Corcoran Clay with a larger resistivity range (5-300 Ωm ) shown in Figure 4 than in Figure 5 (8-100 Ωm ).

of 25
The sections in Figure 4 extend from B, at the western limit of bedrock , to B'. Included on the sections are three driller's logs which describe bedrock and the overlying sediments. In Figure 4a (L 2 -norm) we interpreted the most resistive unit to be the bedrock, but found a smooth transition from the bedrock to the overlying lower resistivity materials. The top of the bedrock, obtained using the minimal variance approach, falls within the depth range corresponding to this transition. Using the L p -norm inversion, with the recovered model shown in Figure 4b, we found an abrupt change in resistivity between the bedrock and overlying lower resistivity materials, the depth of which is in good agreement with what is seen in the driller's logs. Moving from east to west (B' to B), we see a sharp close-to-vertical interface separating resistive bedrock (to the east) from lower resistivity sediments (to the west). These lower resistivity sediments are identified, at shallower depths in the driller's log just east of this interface, as predominantly clay. This lateral location along B-B', which is referred to as B1, is marked as a vertical dashed line. All bedrock in regions westward of this location will be beneath the DOI so not seen in the AEM data. This location is interpreted as the western limit of the bedrock in this transect. This same vertical interface between bedrock and sediments was found for other flight lines within bedrock and defined as the western limit of the bedrock. Figure 5a is a vertical section extending from C to C' (on Figure 3) which displays the starting resistivity model recovered from the first L 2 -norm inversion. The interpreted location of the layer that contains the Corcoran Clay corresponds to the zone of lower resistivity varying in thickness from 80 to 160 m. In Figures 5b and 5c, which extend from C to the northern limit of clay , we show the primary and secondary resistivity models for the Corcoran Clay recovered from the L p -norm inversions. For both of these models the interpreted location of the Corcoran Clay lies within the depth range of the Corcoran Clay-containing layer located using the L 2 -norm inversion, but the thickness is significantly reduced. In the primary model close to the northern boundary of clay , we see a location where the continuous portion of the Corcoran Clay ends, which is referred to as C1; this is defined in this section as the northern boundary of the Corcoran Clay.
In Figure 6, we compare 1D profiles from the three resistivity models: The starting resistivity model, and the primary and secondary resistivity models for the Corcoran Clay with the resistivity logs and sediment type information at Wells A and B within the depth range of 0-250 m from the surface. As expected, for the recovered resistivity model from the L 2 -norm inversion we see a gradual transition in resistivity values around the top and 14 of 25 base of the Corcoran Clay, leading us to (subjectively) identify a Corcoran-Clay-containing layer. The depth to the top of this layer is about 60-80 m and its thickness ranges from 80 to 160 m. We note, that if only an L 2 -norm inversion were available, this Corcoran-Clay-containing layer would likely be interpreted as the location of the Corcoran Clay. In contrast to what is seen in the L 2 -norm resistivity model, the primary resistivity model (from the L p -norm inversion) shows a sharp transition in resistivity values across the top and base of the Corcoran Clay. The interpreted location of the Corcoran Clay agrees well with the driller's logs at Wells A and B, and corresponds to a depth interval where there are significant changes in the resistivity logs. The resistivity values in the model are lower than those in the logs by a factor of 1.5-2, likely to due to a bias toward sampling conductive materials in the AEM method. The secondary resistivity model shows good agreement with the logs from Well B, but fails to match the resistivity values in Well A. This was interpreted to be due either to the underestimation of the reference resistivity value, clay ref or to the fact that the resistivity logging measurement is more sensitive to resistive materials (Geowissenschaften et al., 2007).
We now present our results as 2D maps in Figures 7 and 8. While we have results, from the interpolation process in Step 3, throughout the entire domains of interest for each target, the uncertainty will increase as the distance from an AEM sounding increases. We thus elected to display only those results that were within 3 km of an AEM sounding. We compared the 2D maps to the existing groundwater model by calculating the difference in areas where both our 2D maps and the groundwater model contain the targets.
The 2D map displaying the depth to bedrock is shown in Figure 7a. The first feature to note is the solid red line which represents the western boundary of observed bedrock in the AEM data. This boundary was located by connecting the locations of the most-western soundings where bedrock was present. In some places, the location of this boundary is well-defined, constrained between two soundings that are close together; for example, one  Figure 3). The white gaps indicate locations where the AEM data were not collected due to our design of the AEM survey to maximize coverage while avoiding urban areas, and locations where the data were omitted due to a significant level of noise.

of 25
to the east where bedrock was present and to the west where bedrock absent. In a number of locations, however, the uncertainty in the location of this western boundary is relatively high, increasing as the distance to an AEM sounding increases. There are differences between our mapped western boundary and that defined in the groundwater model that cannot be explained by this uncertainty alone and are not due to a limitation in the DOI of the AEM data. The depth to bedrock in the groundwater model ranges from 100 to 170 m below the surface; this depth range is located above the DOI and is a depth range within which the sensitivity of the AEM data is relatively high. The depth to bedrock determined from the AEM data generally increases toward the west resulting in a thickening of the package of sediments overlying the bedrock. The depth to bedrock, on average 60 m, ranges from 15 m ± 4-160 m ± 14 m, where the error is due to the thickness of resistivity cells at these depths. In Figure 7b, we display the difference obtained when subtracting the depth to bedrock in the groundwater model from the depth that we determined. The depth that we determined was always less, with the difference averaging −80 m and ranging from −150 m to −5 m. Ωm (secondary resistivity model). Resistivity models are shown along C-C' in Figure 3. Horizontal black dashed lines indicate the top and base of the Corcoran-Clay-containing layer interpreted from the L 2 -norm inversion; a vertical dashed line indicates the northern limit of the Corcoran Clay. Vertical dotted line indicates clay and dashed line indicates a lateral location C1 (marked in Figure 3). The white gaps indicate locations where the AEM data were not collected due to our design of the AEM survey to maximize coverage while avoiding urban areas, and locations where the data were omitted due to a significant level of noise.

of 25
In Figures 8a and 8b, we present 2D maps displaying the thickness of and the depth to the Corcoran Clay, respectively. The blue solid line contours where the thickness (determined through interpolation in Step 3) is zero, which we define as the northern and eastern boundary of the continuous portion of the Corcoran Clay. As was the case in mapping the extent of the bedrock, the uncertainty in the location of this boundary increases as the distance to an AEM sounding increases. We again observe differences in the lateral extent of the Corcoran Clay interpreted from the AEM data and that shown in the groundwater model. The Corcoran Clay thickness, as interpreted from the AEM data, increases toward the southwest, with localized thinning or thickening. The thickness averages 17 m, ranging from 3 ± 1 m to 25 ± 7.5 m; with the relative error estimated from the calibration with the well data (∼30% error). The depth to the Corcoran Clay, averaging 100 m, also generally increases toward the southwest, ranging from 50 ± 6 m to 130 ± 12 m; with the error equal to the thickness of the resistivity cells.
In Figures 8c and 8d, we show the differences in thickness and depth-subtracting the values from the groundwater model from the values obtained with our approach. For the Corcoran Clay thickness (Figure 8c) we found relatively large differences of −15 m to −10 m in the region between the boundary of the Corcoran Clay in the groundwater model and the boundary that we determined. Elsewhere the difference ranged from −6 to 5 m with the absolute difference averaging 3 m; the absolute difference was used to take into account the presence of negative and positive differences. The differences in the depth to the Corcoran Clay (Figure 8d) averaged −15 m, ranging from −30 m to −10 m; we always found the Corcoran Clay at a shallower depth. The difference is greater along the eastern and western edges of the region than it is in the center.
In Figure 9a, we show a three-dimensional view of the final resistivity model recovered from the structurally constrained inversion; this covers the entire survey area. Also, we show two vertical sections of the final resistivity model at D-D and C-C', respectively, in Figures 9b and 9c. In the west, the Corcoran Clay divides the upper and lower aquifers, thinning out toward the east, where the upper and lower aquifers merge. In the eastern part of the basin, the resistive bedrock underlies low-resistivity sediments. As expected, we see sharp resistivity contrasts delineating the targets.

Comparison of the L 2 -Norm Inversion and the Targeted Inversion Approach
Let us first consider how the targeted inversion approach, when compared to using the L 2 -norm inversion alone, improved our ability to image the two targets of interest-the bedrock surface and the confining clay layer. As shown in Figure 4a, the top of the bedrock appeared in the L 2 -norm resistivity model as a smooth transition between the resistive bedrock and overlying low-resistivity sediments. As a result, delineating the top of the bedrock from this model was not straightforward, requiring the use of a variance-minimization algorithm. In contrast, the L p -norm resistivity model shown in Figure 4b revealed a sharp interface at the top of the bedrock, that we were able to recover by incorporating the prior knowledge of a large resistivity contrast between the 18 of 25 bedrock and the overlying sediments into the spatial constraints of the L p -norm. The good agreement with the three driller's logs gave us confidence in our approach.
In addition to using the spatial constraints in the L p -norm inversion to locate the Corcoran Clay, the smallness constraint was added so as to also incorporate the prior knowledge that the Corcoran Clay is much thinner than the over/underlying aquifers. As shown in Figure 5, what we interpreted as the Corcoran-Clay-containing layer in the recovered resistivity model from the L 2 -norm inversion was shallower and thicker than the interpreted Corcoran Clay in the driller's and resistivity logs from Wells A and B. In particular, the thickness of the Corcoran-Clay-containing layer was 4-5 times greater than the thickness of the Corcoran Clay layer seen in the well data. This is to be expected. The smoothing of the resistivity model in an L 2 -norm inversion will always result in an overestimation of the thickness of a clay layer in a resistivity background. The Corcoran Clay layer imaged in the primary model from the L p -norm inversion is much closer in thickness to the logs from Wells A and B ( Figure 6) demonstrating the value of the L p -norm inversion for imaging the Corcoran Clay layer. An important aspect of the L p -norm inversion was utilizing high-quality well data, without which we would have not able to select the primary resistivity model. We note that around 80 m depth at Well A, above the top of the Corcoran  When compared to the conventional approach of using an L 2 -norm inversion, we found that the targeted inversion approach, involving a combination of L 2 -norm and L p -norm inversions, yielded significant improvements in our ability to image the depth to bedrock and the confining clay layer in the study area. While requiring more time and input, in the form of prior knowledge, we were able to obtain more accurate information about the locations of the large-scale features.

Comparison of the Results Obtained From the AEM Data to the Existing Groundwater Model
The existing groundwater model from our study area was constructed using relatively dense well data. The fact that this model was not utilized in our targeted inversion approach provided us an opportunity to assess the efficacy of our approach by comparing the locations of the targets we obtained to those in the existing groundwater model.
We start by discussing the location of the bedrock surface, where we found a relatively large difference between our approach and the groundwater model. The location we determined for the boundary marking the western extent of the bedrock (within the depth range of the groundwater model and the AEM data) differed by ± 5 km from the boundary location in the groundwater model. In addition, the average difference between the depth estimates was about −80 m corresponding to a relative difference of roughly 130%. In contrast to the flat top of the bedrock in the groundwater model (as illustrated in Figure 1b), our mapping showed considerable topography (i.e., variation in depth), as would be expected given the erosional history of an exposed bedrock surface. Given the lack of well data providing information about the location of the bedrock surface in the groundwater model, it is very likely that neither the extent of the bedrock surface nor the depth to the bedrock are accurately captured in the groundwater model. With the AEM data, we were able to obtain accurate imaging of the bedrock surface-in those areas where AEM data are available. The limitation with the AEM data is the lack of coverage in some areas due to the spacing of the flight lines, and the presence of urban areas (where AEM data acquisition is not allowed); but the coverage is far superior to that of the well data. In comparing the two approaches to mapping the bedrock surface, we have much higher confidence in the location obtained through targeted inversion than that shown in the groundwater model.
In building the groundwater model, a greater number of wells were available in the region of the Corcoran Clay than in the region where the bedrock surface is present. It is, however, very difficult to accurately map out the Corcoran Clay using the well data given the abundance of other clay interbeds in the upper and lower aquifers. The Corcoran Clay is blue to gray or bluish green in color so is commonly referred to in driller's logs as "blue clay". Unfortunately, the color cannot conclusively identify the Corcoran Clay as many clays in the area that lack iron oxide are blue or gray. Further, even if the description in the driller's log does not say blue, the clay might be part of the Corcoran clay, as the driller may have decided that all clays in the area were blue so it was not worth repeating. Given all of this, mapping the Corcoran Clay with AEM data can provide a valuable complement to the interpretation based on the well data.
Let us first discuss what is shown as the northern and eastern boundary of the Corcoran Clay. As seen in Figure 8a, the boundary of the Corcoran Clay that we located is about 3-7 km inside (to the west of) the boundary in the groundwater model, except for a small area where our boundary extends to the east of the model boundary. In Figure 10, we compare the mapping of the Corcoran Clay-as either present or absent, in the AEM data and in the well data. In the well data we used the description "blue clay" within the depth range of 50-150 m in the driller's log as an indication that Corcoran Clay was present at that location. In the region where the boundary of the Corcoran Clay in the model lies to the east of what we determined, there are only three wells with "blue clay" and 18 without. Along the southern stretch of this region, our mapped boundary seems to show better agreement with the well data than the groundwater model boundary does. As evidence of the challenge of using well data to locate the Corcoran Clay in the groundwater model, it is interesting to note the 13 reports of "blue clay" in wells to the east of the boundary in the model. The thickness of the "blue clay" in these wells was on average 10 m, thick enough to be detected by the AEM data if it was laterally continuous.
One likely explanation for the difference in the boundary locations, where the AEM boundary lies inside that in the groundwater model, is the challenge of using well data to differentiate the continuous Corcoran Clay from discontinuous Corcoran Clay, or from the occurrence of other clays at a similar depth. In reviewing the resistivity models recovered from the AEM data, we see clay-rich zones in the region between the two boundaries, but these 21 of 25 zones are not continuous with the Corcoran Clay. An example of this is shown in Figure 5, where clay is present to the south of C1 but it is not continuous with the Corcoran Clay. It is possible that the Corcoran Clay becomes discontinuous in this region, but was mapped as part of the continuous Corcoran Clay in the groundwater model. It is also possible that these clay-rich zones are not part of the Corcoran Clay at all, but just identified as such in the groundwater model.
In the area where our determined boundary crosses over and goes to the east of the model boundary, our boundary is not well-constrained by the appearance/absence of "blue clay" in the driller's logs. While there are limited AEM data due to the presence of an urban area, the density of soundings and sensitivity of the AEM measurement to the presence of the Corcoran Clay, gives us a high level of confidence in using the interpolated AEM data to map the extent of the Corcoran Clay in this area.
The advantage we have in the AEM data, is the superior ability to see the parts of the continuous Corcoran Clay layer proximate to the AEM sounding locations compared to the sparse well data. As such, we tend to put greater confidence in the location of the boundary determined from the AEM data. The key limitation is the high level of uncertainty in areas with poor AEM data coverage. If the boundary location in these areas has a significant impact on the accuracy of the groundwater model, additional data such as well data and ground-based electromagnetic data should be acquired. In regions where the Corcoran Clay was identified in the AEM data and the groundwater model, we identified small to moderate differences in Corcoran Clay thickness. In the groundwater model the Corcoran Clay thickness ranged from 10 to 20 m; in our 2D map it ranged from 3 ± 1 m to 25 ± 7.5 m. The relative difference was ∼15%. The depth to the Corcoran Clay in the groundwater model ranged from 50 ± 90 m to 160 m ± 12 m; in our 2D map it ranged from 50 ± 6 m to 130 m ± 12 m, again resulting in a relative difference of ∼15%.
What defines the targeted inversion approach is the focus on imaging specific targets. We have done this by including the L 2 -norm inversion, the L p -norm inversion, and the interpolation process in our approach. This developed approach is transferrable to any other hydrogeologic setting where specific targets of interest can be defined. For example, in a coastal area where the change in salinity dominates the resistivity response, the inversion would be set up to recover the location of the freshwater-saltwater interface by setting the -values for the spatial constraints to be small (e.g., 0 and 1) to reflect the large resistivity contrast at the interface.

Further Use of the Results From the Targeted Inversion
While the targeted inversion approach provided accurate locations of the large-scale features of interest in this study area, there is additional information at a smaller scale that can be obtained about the groundwater system from the AEM data. This requires a high-quality resistivity model over the entire study area. To obtain such a model, we implemented a structurally-constrained inversion where information about the location of large-scale features is used to constrain the inversion. Obtaining such a model was the goal of the structurally constrained inversion where information about the location of large-scale features was used to constrain the inversion.
In Figure 9a, we show a three-dimensional view of the final resistivity model covering the entire subbasin. As expected, we see sharp resistivity contrasts delineating the targets. In addition to delineating the targets-as had been done with the targeted inversion approach, the final resistivity model recovered from the structurally constrained inversion also identified in the presence of a package of lower permeability sediments below the lower aquifer and information about smaller-scale structure within the upper and lower aquifers; this information can also be used to improve the existing groundwater model. Figure 1b, the existing groundwater model is defined with the base at the top of a package of "low-permeability sediments" throughout much of the Kaweah subbasin and at the bedrock surface along the eastern edge. When we review a vertical section through the final resistivity model in Figure 9b (along transect D-D' shown in Figure 9a) we see low-resistivity materials adjacent to the bedrock in the east (near D') that then continue to the west for ∼30 km extending to the DOI. We interpret these to be clay-rich sediments equivalent to the package of "low-permeability sediments" described in the groundwater model. The top of these sediments is above the base of the groundwater model along most of the transect. This suggests that either the base of the model needs moving to shallower depths or a new layer needs to be added to the lower aquifer. As can be seen in the 3D view of the resistivity model in Figure 9a, the package of low permeability sediments is present throughout the eastern half of the subbasin. This has important implications for groundwater flow within the subbasin and for recharge coming into the valley from the mountain block.

As shown in
The resistivity model recovered from the structurally constrained inversion provides more accurate information about the resistivity values in the upper and lower aquifers. In Figure 9c, we see a clearly imaged layer in the upper aquifer, just overlying the Corcoran Clay, with resistivity values higher than those found in the lower aquifer. This is to be expected given the reports of more clay in the lower aquifer than the upper aquifer (Fugro West, 2016). This resistivity feature was not clearly identified in the resistivity model recovered from either the L 2 -norm inversion (Figure 5a) or from the L p -norm inversion (Figure 5b), but was clearly identified in the model from the structurally constrained inversion (Figure 9c). This demonstrates the value of the structurally constrained inversion, where knowing the large-scale structure allows for improved imaging of the smaller-scale features.

Conclusions
Accurate groundwater models, to support groundwater science or management, require as input information about the large-scale structure of the groundwater system. We conclude that a multi-step targeted-inversion approach provides an effective way to extract the most accurate information from AEM data about the large-scale structure. Having defined the large-scale structure, this can then be used to improve the recovery of the smaller-scale resistivity structure.
By implementing an L p -norm inversion, we were able to incorporate prior knowledge. This made it possible to refine the recovered resistivity images so as to more accurately locate the targets. In implementing an L p -norm inversion in other areas, the needed prior knowledge should not be difficult to obtain if information is available about the types of geological material present and the expected changes in resistivity associated with the largescale features. The presence of high-quality well data would provide an additional source of information that could be incorporated into the inversion and would definitely result in improved imaging; but we were successful in this study with data from only six wells.
The targeted inversion approach includes multiple L 2 -norm and L p -norm inversions as well as an interpolation process, so it requires more computation and time than other approaches. We conclude, however, that this can easily be justified by the benefits obtained in terms of improved imaging. In addition, running multiple inversions 23 of 25 could readily be parallelizable, and many of the other processes in the approach (e.g., interpolation) could be automated.
In the adoption of AEM data for the development of groundwater models, the optimal approach is to first work with all existing well data and other sources of information to develop a groundwater model, and then determine where the AEM data could be most valuable in reducing uncertainty in the model. All available data, included the acquired AEM data, would then be integrated to develop a model that fits all sources of data. Such an undertaking would require a significant effort with hydrogeologists working closely with geophysicists.
In this study, we elected to use the existing groundwater model solely as a means of assessing our ability to extract information about the large-scale features from the AEM data. This allowed us to address the key question: what information can we obtain from AEM data to support the development of a groundwater model? This study has shown the information about large-scale structure that can be obtained from AEM data. There will inevitably be uncertainty in locating large-scale features, uncertainty that-as with well data, increases as the distance from the AEM data increases.
Given the growing use of the AEM method for groundwater science and management, we hope there will be continued development of our approach for the inversion of AEM data and integration with other forms of data, so that we can maximize the benefit of all available datasets. To accelerate this, we have publicly released the numerical codes used in this study through a Python-based open-source software, SimPEG (https://www.simpeg. xyz).
= relative error (%) × 0.01 × | | + floor (A4) The relative error takes into account potential random noise or motion noise that remained after the stacking process, the amplitude of which is proportional to the amplitude of the data. The floor is related to ambient noise or system noise, which may have a constant amplitude. The relative error and floor were set to 3% and 10 −15 V/A-m 4 for most of sounding locations except for the locations close to the eastern edge of the survey area; estimates for these values were obtained through an iterative inversion processes-running an inversion and checking how well the data could be fit without generating artifacts in the recovered resistivity model. We found that the signal-to-noise (S/N) ratio at the eastern edge was relatively low due to shallower bedrock surface and the impact of rapid topographic change on the motion noise. In addition, this rapid topographic change can introduce 1D modeling errors. Given this decreased level of S/N ratio and the potential 1D modeling error at the eastern, we assigned the greater level of data error: 10% and 10 −14 V/A-m 4 . The target misfit, * , was set to assuming the chi-square distribution of the data error (Oldenburg & Li, 2005).

Data Availability Statement
The AEM data used in this study are publicly available through the Stanford Digital Library: https://purl.stanford. edu/yc041mz9691. Funding for this research was provided to Rosemary Knight by the Gordon and Betty Moore Foundation (Grant No. GBMF6189) and the NASA Applied Sciences Water Resources Program (award no. 80NSSC19K1248). Additional funding for the acquisition of the AEM data was obtained from a Proposition 68 Sustainable Groundwater Management grant from the California Department of Water Resources to water agencies in the Kaweah subbasin. We would like to thank Aqua Geo Frameworks for providing the processed AEM data, Larry Dotson from Kaweah Delta Water Conservation District for providing the geophysical logs, and Klara Steklova for digitizing the driller's logs. We also wish to thank Aaron Fukuda, Michael Hagman, and Eric Osterling from the Kaweah Subbasin Groundwater Sustainability Agencies, as well as Derrik Williams (Montgomery and Associates) for providing, and assisting with our use of the groundwater model. We are grateful to the SimPEG community for developing and publicly releasing the SimPEG-EM1D module used in this study. Lastly, we would like to thank Andrew Binley and three anonymous reviewers for their insightful feedback on our manuscript. [Note to editors and reviewers-the following will be included once the paper has been accepted for publication].