Recharge site assessment through the integration of surface geophysics and cone penetrometer testing

The ability to identify, at potential managed aquifer recharge sites, the presence of connected pathways of hydraulically conductive sediments from the ground surface to the water table could help minimize costs and risks associated with recharge operations. A spatially dense dataset had previously been acquired in an almond [Prunus dulcis (Mill.) D.A. Webb] grove in Tulare, CA, using tTEM, a towed transient electromagnetic (tTEM) geophysical method. In order to interpret reliable information about sediment type from the tTEM data, a transform from the tTEM‐derived property, electrical resistivity, to sediment type is required. The uncertainty associated with derived models of sediment type can be significantly reduced if a site‐ and dataset‐specific transform is used. Cone penetrometer testing (CPT) was conducted at five locations, strategically selected based on a review of the tTEM data. Co‐located measurements of sediment type, derived from the CPT, and electrical resistivity, derived from the tTEM data, were used to create a resistivity‐to‐sediment‐type transform, with sediment type classified as either coarse‐grain‐dominated (sand and gravel) or fine‐grain‐dominated (silt and clay) material. The transform captured the uncertainty associated with variable water salinity and content, the resolution of the tTEM data, and other components of the tTEM measurement workflow. Using this transform, models of sediment type were generated for the unsaturated zone at the site. Within these models are features, which we interpret as potential recharge pathways, corresponding to high fractions of coarse‐grain‐dominated material amongst regions of fine‐grain‐dominated material. The workflow developed at this site can provide a framework for using tTEM and CPT for recharge site assessment.

groundwater levels have been in decline for decades. In 2014, the California Legislature passed the Sustainable Groundwater Management Act (SGMA). The SGMA requires that groundwater be managed sustainably, defining sustainability as avoiding six key undesirable results, most of which come directly from the overdraft of groundwater resources. In order to meet the sustainability mandate imposed by the state, and to ensure groundwater supply in the face of future droughts, local water management agencies must find a way to balance the supply and demand of groundwater within their districts. One method of reaching this balance, being considered by many management agencies, is to increase groundwater recharge in times of deemed surface water excess, allowing the aquifer systems to serve as a storage location for water. This process is commonly referred to as managed aquifer recharge (MAR).
There are many ways to facilitate groundwater recharge, with varying degrees of intervention with the natural recharge process. These range from applying water to the ground surface and allowing it to infiltrate, to active injection of water into the subsurface below the water table (Asano, 1985). Here we focus on one recharge method, where water is applied at the ground surface and allowed to infiltrate, either within a dedicated recharge pond or on agricultural fields, the latter commonly referred to as Ag-MAR.
There are a number of factors that must be considered when assessing the effectiveness of a potential site for recharge. These include, but are not limited to, the properties of the surface and subsurface materials, topography, water conveyance infrastructure, and water availability. Although all of these factors are important, this study is focused on developing a methodology for identifying sites most suitable given their subsurface properties. In order to identify such recharge sites, the subsurface of each site must be characterized to determine the presence of interconnected pathways of sediments with high hydraulic conductivity, which can serve to efficiently move water, that has infiltrated through the surface layer, down to the water table. Generally, coarse-grained material, such as sands and gravels, tend to have high hydraulic conductivity, whereas fine-grained material, such as clay, tend to have much lower hydraulic conductivity. Often, soil maps are used to characterize the sediments at a potential recharge site (O'Geen et al., 2015). As these only contain information about the shallow soil layer, they can be used to assess initial infiltration but cannot provide any information about sediment type below this shallow. Where more information about the subsurface is desired, the next step is often to conduct drilling or coring. Although this approach can be used to map changes in sediment type with depth, it is costly, invasive, and only yields information about the subsurface at the specific points at which measurements are made. Depending on the level of heterogeneity in the subsurface, even with this type of data, there can still be significant uncertainty in understanding the potential for recharge at a site. What is needed are more cost-

Core Ideas
• tTEM data can be useful for sediment type mapping the subsurface at potential recharge sites. • The quality of a tTEM derived sediment type model depends on the resistivity transform used. • A site-and measurement-specific transform can be created using collocated CPT and tTEM data. • Using a high-quality transform, potential recharge pathways can be located in the subsurface.
effective methods for subsurface characterization at a spatial resolution that capture the locations of recharge pathways at the scale of a recharge site. A number of geophysical instruments have proven useful for recharge site characterization, as they can cost-effectively yield information about spatially varying subsurface conditions (Cook et al., 1992;Gottschalk et al., 2017;Lawrie et al., 2012;Maliva et al., 2015;Mawer et al., 2013Mawer et al., , 2016Parsekian et al., 2014;Sendrós et al., 2013;Walsh et al., 2014). One such instrument, recently developed in Denmark, is the tTEM system. The tTEM system is a transient electromagnetic (TEM) system, which is towed (the "t" in tTEM) behind an all-terrain vehicle (ATV), allowing for rapid data acquisition . Fortuitously, the system is of a size that fits between the rows of tree crops which are prevalent in the Central Valley of California, (USA), so was used in a recent study to assess potential recharge sites in the valley (Behroozmand et al., 2019). The data acquired with this system can be used to generate detailed models of the electrical resistivity of the subsurface down to 60+ m, spanning the depth intervals of greatest interest for recharge operations in much of the Central Valley. Maps of subsurface resistivity are useful for recharge site characterization and assessment because changes in resistivity generally correspond to changes in sediment type: coarse-grained material (such as sands and gravels) tend to have much higher resistivity than clay-rich material. The challenge with using measurements of resistivity to map sediment type, however, is that changes in resistivity can result from changes not just in sediment type, but also in water content and pore fluid salinity. Because of this, additional information is needed to establish a resistivity-to-sediment-type transform that can be used to interpret the resistivity models derived from the tTEM data.
There are many approaches to creating the resistivity-tosediment-type transform. As recently reviewed by Knight et al. (2021), this can range from applying empirical relationships derived in the laboratory at the core scale, to the development of more complicated, and likely more accurate, transforms at the field scale; the latter can incorporate spatial variability in both geophysical and material properties, uncertainty in the underlying data, and issues associated with the measurement workflow in acquiring, processing, and inversion of the geophysical data. In many cases, these field-scale transforms rely on independent measurements of sediment type described in drillers' logs from wells (Christiansen et al., 2014;Foged et al., 2014;Høyer et al., 2015 ). One challenge with using well data to build the transform at a potential recharge site is that there are typically few, if any, wells close enough to make a reliable transform.
Cone penetrometer testing (CPT) is one method of acquiring an independent measure of sediment type that can be used to assist with the interpretation of tTEM data. This method uses data collected from a cone-tipped rod pushed into the subsurface to generate soil behavior type (SBT) logs. This approach has proven successful for characterizing the shallow subsurface for a wide variety of geotechnical and groundwater applications (Coutinho & Mayne, 2013;Robertson, 2016;U.S. Department of Transportation, 2017;Vienken et al., 2012).
In this study, we developed a methodology for transforming resistivity models to sediment type in the vadose zone, constructing the field-scale transform using tTEM-derived resistivity and CPT as the measure of sediment type. Our methodology directly incorporates the tTEM-derived resistivity into the development of the transform. This provides a transform that takes into account the measurement workflow (data acquisition, processing, inversion) used to obtain the resistivity models. We tested the utility of our transform-development method for identifying potential recharge pathways, using data from a site in the Central Valley, where we had both tTEM-derived resistivity and CPT datasets. The method, however, is applicable to any site where there is a need to image the location of potential recharge pathways in unconsolidated material in the top ∼60 m beneath the ground surface. , and the deepest was 52 m below the surface (in 2014). In an effort to combat negative impacts associated with declining water levels, and to work towards becoming compliant with the new mandates of SGMA, the Tulare Irrigation District operates a number of MAR projects, with plans for expanding these operations in the future (Mid-Kaweah GSA, 2019). Given both the urgent need to increase recharge in the region, and the interest of the local water management in facilitating Ag-MAR, the study area is an ideal location to investigate the application of geophysical methods to recharge site assessment. Understanding the geologic setting of our study area is key to understanding the nature of potential recharge pathways and interpreting data collected at the site. The following summarizes material contained in the groundwater sustainability plan recently prepared by the Mid-Kaweah Groundwater Sustainability Agency. The ∼50,000 km 2 of the Central Valley are largely filled with sediments eroded from the Sierra Nevada Mountains to the east, and the Coast Ranges to the west. The valley is filled with continental, deltaic, and marine sediments, which can be up to 16 km thick. In the San Joaquin Valley, the uppermost sediments are of continental origin. Alluvial fans are the dominant type of deposits preserved along the valley margins, especially along the eastern edge, with river and flood basin deposits filling in between. In the vicinity of our study area, geologic studies have mapped younger alluvium, consisting of interbedded deposits of gravelly sand, silty sand, silt, and clay, in the top ∼30 m of the subsurface. Below this, there is a transition to older alluvium, which is composed of highly oxidized deposits that tend to be more clay-rich, and less permeable than the overlying sediments. This geologic understanding of the region suggests a high degree of spatial heterogeneity in the subsurface sediments, underscoring the importance of characterizing potential recharge sites in our study area.

Description of the study area
It is important to note that recharge through the vadose zone, if the sediments are very dry, involves first a phase of initial wetting. This is a complex process, in which the contributions of different sediment types to governing fluid flow may vary greatly over time. For the purposes of this study, and any discussion of recharge pathways, we assumed that the sediments had undergone this initial wetting. We presumed that coarse-grained material would enhance recharge, and fine-grained material would inhibit recharge in the vadose zone. This characterization is consistent with the work of others on recharge site assessment in the Central Valley (Maples et al., 2019).

2.2
Introduction to use of the tTEM system for mapping sediment type A detailed description of the tTEM system can be found in Auken et al. (2019). A tTEM measurement is made by passing a time varying current through a transmitter loop. This induces a primary magnetic field, which decays when the transmitter current is turned off. The decaying primary magnetic field induces eddy currents in the electrically conductive portions of the subsurface. These eddy currents generate a secondary magnetic field, the time decay of which is recorded by the receiver as dB/dt (magnetic field decay over time), with the earlier time measurements corresponding to signal coming from the shallower subsurface, and the later time measurements corresponding to signal coming from deeper in the subsurface. The system uses a high and low moment measurement, in order to maximize both resolution of the shallow subsurface and depth of imaging. Hundreds of measurements are made every second, with the system typically towed at speeds from 15 to 20 km h −1 . Processing includes removing noisy measurements and averaging and stacking measurements over a predetermined window along the acquisition line, resulting in soundings (stacked measurements of the magnetic field decay) at regular intervals (generally ∼10 m). Typically, over 1,000 low moment and high moment measurements are included in each sounding. The magnetic field decay in the soundings is presented as single values within a number of short time windows, referred to as time gates. Typically, there are 15 time gates in the low moment measurement, and 23 time gates in the high moment measurement.
Inversion of the acquired data involves solving for the resistivity model at each sounding location which produces simulated tTEM data that agree with the acquired tTEM data. There are a number of different approaches that can be taken to the inversion of these data, ranging from simple one-dimensional (1D) inversions, to full three-dimensional (3D) inversions. The 1D spatially constrained inversion (Viezzoli et al., 2008) was implemented in this study. This is a standard, computationally efficient approach for the inversion of tTEM datasets for hydrogeologic applications. In this approach, the inversion of each individual sounding uses a 1D forward model, which assumes changes in resistivity only occur in the vertical direction. This 1D inversion is, however, constrained laterally by the model parameters of surrounding soundings, thus producing laterally smoothed results between neighboring resistivity models. The resulting resistivity models are discretized so as to represent the vertical resolution and are typically truncated at the depth beyond which estimates of resistivity are deemed to be less reliable, referred to as the depth of investigation (DOI).
In hydrogeologic studies, often the last step in working with any geophysical data is to transform the models of geophysical properties, recovered through the inversion of the acquired data, in this case resistivity, to the material property of interest, in this case sediment type. As recommended by Knight et al. (2021), in this study, we used a transform specifically developed for field-scale application at the field site of interest. Of particular relevance to this study, Knight et al. (2018Knight et al. ( , 2021 and Dewar and Knight (2020) noted the critical importance, when working in California, of having a measure of the top of the saturated zone (TSZ) when transforming fieldscale resistivity data, because of the large thickness of the vadose zone in many places and the impact of water content on resistivity. To account for this we implemented the method developed by Dewer and Knight (2020) for locating the TSZ.
In this study, we developed and tested an approach to transforming tTEM data to models of sediment type. This approach relies on co-located measurements of resistivity and sediment type, accounts for the difference in scale in these measurements, accounts for the tTEM measurement workflow, and seeks to capture uncertainty in the transform. We focus specifically on capturing uncertainty associated with variable lithology, saturation, water salinity, and measurement scale. We note that there is also uncertainty associated with the inversion of the tTEM data, the inclusion of which is an important area of future research.

2.3
Introduction to use of the CPT method for determining sediment properties Cone penetrometer testing measurements are made by pushing an instrumented, cone-tipped rod into the subsurface. An introduction to CPT measurements can be found in the paper by Robertson and Cabal (2015). Two key measurements are made while the rod is pushed: the normal force against the cone tip and the shear force against the side of the rod above the tip. Normalizing these two force measurements by area yields the tip resistance (q c ) and sleeve resistance (f s ). An empirical relationship exists between these field measurements and SBT, which classifies sediments in terms of their physical and mechanical properties. This relationship can be further refined by accounting for the variable overburden pressure, resulting in a classification called the normalized SBT (SBTn). The empirical relationship used in this study follows the updated classification from Robertson (2010) and is shown in Figure 1. In making this updated relationship, Robertson reviewed empirical relationships from over 33 studies, with CPT sites around the world. Of key importance in our study is the difference seen between claydominated sediments, which tend to have a low tip resistance coupled with high sleeve resistance, and sand-or graveldominated sediments, which tend to have a high tip resistance paired with a low sleeve resistance.
The cones used to make these measurements can range in size but are commonly 35-42 mm in diameter. Measurements are made continuously while the rod is pushed, with data reported in sample intervals generally ranging from 2 to 10 cm. The depth to which these cones can be pushed depends on the subsurface material properties and cone size but, with the right equipment and soil conditions, can exceed 100 m. Given the size of the cone, and the methods by which measurements are made, the produced logs only represent the sediments in a small volume around the cone.

Vadose Zone Journal
F I G U R E 1 Normalized soil behavior type (SBTn) classification, showing the relationship between normalized cone resistance and friction ratio. Modified after Robertson (2010) Other co-located measurements can be made by including additional instruments on the rod. An example is the measurement of electrical resistivity made by injecting current between two electrodes and measuring the change in voltage between another two electrodes. The electrodes are generally spaced 10-150 mm apart, resulting in a measurement that is sensitive to material a small distance from the rod. The exact distance depends on the electrical resistivity of the surrounding material and the geometry of the electrodes but is typically assumed to be, at most, 30% of the length of the electrode array. We refer to the produced logs of electrical resistivity as cone penetrometer testing electrical resistivity (eCPT) logs.
It is important to note that although the CPT method often results in a description of sediments similar to what would be obtained through direct sampling, it does not involve any direct sampling of sediments. For this reason, depending on the application, it may be important to validate the classification from this method with other datasets. In this study, we justified the use of this method for sediment typing using both the eCPT logs, and experience with this method in other parts of the Central Valley, where co-located CPT and Geoprobe cores show very similar classifications of sediments (G. Osterman, personal communication, 2020). We do note, however, that this is a potential source of uncertainty in the final resistivityto-sediment type transform.

DATASETS AND METHODS
This study used two datasets: tTEM data acquired in October 2017, and CPT data acquired in October 2019. Figure 2 shows the locations of these two datasets within the almond grove. In this figure, the location of the 1D resistivity models resulting from inversion of the tTEM data are shown in blue, and the locations of the CPT datasets are shown in yellow.

tTEM resistivity models
The resistivity models were obtained in the study conducted by Behroozmand et al. (2019). All details of data acquisition, processing, and inversion are given in that study and thus are only briefly summarized here. In October 2017, 43 km (commonly referred to as line-km) of tTEM data were acquired in the almond grove, over ∼6 h. These data were collected along rows spaced ∼6.8 m apart. These data were stacked during processing, resulting in a sounding every ∼10 m along each row, and then inverted in Aarhus Workbench (Aarhus Geosoftware) using a 1D spatially constrained (SCI) approach. The data were inverted using a 30-layer smooth model, where the layer boundaries were fixed. Layer thickness in the resulting resistivity model begins at 1 m at the ground surface and increases with depth to reach a thickness of 6.8 m at ∼60 m, which was determined to be the average DOI. The DOI for each model was estimated using the approach outlined in Christiansen and Auken (2012), where a global sensitivity threshold, applied to the recalculated sensitivity matrix, is used to determine where there is no longer sufficient sensitivity of the electromagnetic method to subsurface changes in resistivity. At this depth, it is assumed that reliable estimates of resistivity cannot be recovered through inversion of the data.

F I G U R E 2
Map showing the study area, located just east of Tulare, CA. The locations of the two datasets collected at the sites are shown, with the resistivity models resulting from inversion of the towed transient electromagnetic (tTEM) data in blue and cone penetrometer testing (CPT) datasets in yellow F I G U R E 3 Two views of the one-dimensional (1D) resistivity models resulting from inversion of the towed transient electromagnetic (tTEM) data, underlain by a satellite image of the study area. Warm colors (reds) correspond to high resistivity values; cool colors (blues) correspond to low resistivity values. Panel a shows the full extent of these models, starting at the ground surface and extending to the depth of investigation (DOI, between 40-85 m below ground surface [mbgs]). Panel b shows the models from 10 mbgs to the DOI Figure 3 shows two views of the 1D resistivity models resulting from inversion of the tTEM data, conducted by Behroozmand et al. (2019). In these images, each model has been clipped to the DOI at that location; this ranges from 40 to 85 m below the ground surface (mbgs) across the study area. Variation in the DOI across the site is a function of both variable subsurface resistivity, as well as noise levels across the field. At this site, the DOI became shallower as we approached the western edge of the field where there is a transmission line. Repeat surveys would be expected to have similar DOIs; however, some variation could occur if there were changes in noise sources or subsurface resistivity. The images are colored according to the resistivity values, with warm colors (reds) corresponding to high resistivity values, and cool colors (blues) corresponding to low resistivity values. The color scale is clipped at 50 Ω-m, to highlight the variation in resis-tivity within the range most pertinent to changes in sediment type, but resistivity values in these models reach as high as 130 Ω-m. However, 80% of the modeled resistivity values are contained within the range displayed in Figure 3. Figure 3a shows these models beginning at the ground surface, whereas Figure 3b shows these models beginning at 10 mbgs.
These two images show significant spatial heterogeneity in the subsurface resistivity as mapped by the tTEM, both vertically and laterally. In particular, at 10 mbgs, shown as the top of the image in Figure 2b, there are distinct regions that are dominated by high resistivity material, and regions dominated by low resistivity material. These resistivity models suggest that within this relatively small site, there may be complex patterns of sediment type. This complexity could play an important role in determining how water applied at the surface would infiltrate and percolate down to the water table.
F I G U R E 4 Electrical resistivity (eCPT), nearest towed transient electromagnetic (tTEM)-derived resistivity models (from both smooth and sharp inversions), and normalized soil behavior type logs (SBTn) for two cone penetrometer testing (CPT) locations (CPT1 and CPT4) within the study area In order to interpret these images in terms of sediment type, however, a resistivity-to-sediment-type transform is required.
Although this study implemented the smooth SCI inversion approach, it is important to note that, in any use of EM methods, the form of the resistivity models recovered from the data will be affected by the selected inversion approach. Each approach has its advantages and disadvantages, which we discuss by considering two scenarios encountered at the field site. Let us first consider the scenario where there is an interface at which there is a significant change in resistivity. In selecting a smooth inversion, what would likely appear in the recovered resistivity model is a gradual transition from the resistivity above the interface to that below, as the regularization enforces a smooth transition. In selecting a sharp inversion, the recovered resistivity model would likely display one or more abrupt changes in resistivity as a way to transition between the resistivity values above and below the interface. Given the resolution of the tTEM measurement, there would always be uncertainty in the depth locations at which these abrupt changes occurred (unless their locations were constrained).
An additional scenario to be considered is one where there is a thin clay layer, with a thickness below the resolution of the tTEM measurement. The use of a smooth inversion would likely result in a small change in resistivity at the depth corresponding to the clay layer. The use of a sharp inversion would likely result in little to no change in resistivity at the depth corresponding to the clay layer, as the regularization penalizes small changes in resistivity. In this study, we wanted to create a transform that could capture information about variation in sediment type below the resolution of the tTEM measurement, so we wanted to ensure that the resistivity value would be sensitive to such thin layers. As a result, we elected to use a smooth inversion approach. We note, however, that the methodology developed in this study could be applied to the resistivity models recovered through any inversion approach.
For illustrative purposes, we have included in Figure 4 below the results of both a smooth inversion and a sharp inversion. In this figure, we see that although the two approaches resulted in very similar resistivity models, there are subtle differences. In particular, we note that the gradual transition in resistivity values seen in the models recovered from a smooth inversion-between 2 and 30 m in CPT1 and between 5 and 10 m in CPT4-both appear as a single abrupt change in resistivity in the models recovered from a sharp inversion.

CPT
In October 2019, CPT testing was conducted at five locations, in order to acquire five SBTn and eCPT logs (both acquired at the same time and place with the one CPT rod). These logs were collected over 2 d in the almond grove using a 20-metric ton track-mounted rig. Both of these logs were continuous, with sampling intervals of 0.1 m for the SBTn logs and 0.05 m for the eCPT logs. The CPT testing locations were selected so as to sample a range of the resistivity values seen in the tTEM models, and to intersect key continuous features observed in the resistivity models. The CPT rod was pushed at each point until refusal: the point at which the force exerted by the pushing rig was insufficient to overcome the tip resistance and side friction against the rod. The shallowest refusal occurred at 19 mbgs, and the deepest occurred at 32 mbgs, with an average depth of 25.6 mbgs. Figure 4 shows examples of the SBTn and eCPT logs from two locations at the site, as well as two resistivity models from the tTEM measurement locations closest to each CPT point. One of the resistivity models was derived from a smooth inversion of the tTEM data, and one was from a sharp inversion. The SBTn logs show that the subsurface predominantly consists of SBTs 3-7, which include clays, silts, sands, and gravely sand. In comparing the SBTn and eCPT logs, we see a correspondence that matches our expectation as to how the electrical resistivity will change with sediment type: the sands and gravels exhibit the highest resistivity values, with the clays exhibiting much lower resistivity values. We note that there are significant differences between the eCPT logs and tTEM resistivity models, with the tTEM resistivity models often looking like a smoothed version of the eCPT log.
There are a number of reasons for the differences. The first is that the eCPT tool was optimized by the manufacturer (as described in the system specifications) for lower resistivity values (i.e., the tool is most sensitive at values ≥20 Ω-m). This is below the range of resistivity values observed by the tool throughout most of our field site. As a result, although we have confidence in the depth and relative changes in resistivity values observed in the eCPT logs, we do not expect the measured magnitude of the resistivity values to be accurate.
Another reason for the differences in the eCPT and tTEM resistivity logs is the differences in the physics of the two forms of measurement, one being direct current, and the other being electromagnetic. The tTEM measurement has worse resolution and a larger support volume, both of which will result in averaging, and thus smoothing, of variation in resistivity values. The other cause of the significantly lower resistivity values seen in the tTEM data in the resistive sections is due to an inability, common to all electromagnetic methods, to resolve differences in resistivity once the resistivity is above a certain threshold; this is commonly referred to as a saturation of the resistivity.
For these reasons given above, the eCPT data in this study were used in a qualitative way to assist with defining the classification system for the SBT data but were not directly used in the development of the resistivity-to-sediment-type transform. Given our need to transform tTEM resistivity values to sediment type, we used the tTEM resistivity values to develop the transform.

Estimation of depth to the TSZ
With the goal of this work being to locate connected recharge pathways from the surface to the water table, the first step in working with the tTEM resistivity models was to determine the location of the water table and clip the tTEM resistivity models at this depth. Water table measurements are taken in wells in the vicinity of the study area twice a year; however, the elevation of the water table can vary significantly on short timescales as a result of pumping. Because of this, we elected not to use the water table measurements from wells; instead, we used the tTEM data to estimate the depth to the TSZ using the method of Dewar and Knight (2020). The TSZ is located above the water table, as it includes the capillary fringe. We presumed that the TSZ would be close to the water table at the field site due to the presence of thin layers of coarse sediment, which are expected to limit the extent of the capillary fringe. In addition, given the purpose of the study, mapping permeable pathways to the TSZ would be as useful as mapping pathways to the water table.
The Dewar and Knight (2020) approach estimates the TSZ as the depth at which there is a significant change in statistical properties of the distribution of electrical resistivity values, within a prescribed depth search window. This search window is critical, as variation in lithology can have a large impact on resistivity distributions so the depth at which lithology changes could therefore be misidentified as the TSZ. Using the well-based water table measurements, we defined the search window to be between 40 and 60 mbgs; the minimum and maximum water table depths measured for the 5 yr prior to the tTEM acquisition.
To execute the Dewar and Knight method, all of the tTEM inverted resistivity models from the study area were subsampled to 1-m depth intervals and were then collectively used to create a single resistivity distribution for each depth interval. The difference between the 25th and 75th percentile of each of these distributions was calculated, and the TSZ was assigned to the depth at which the largest change in this statistical property occurred, within the depth search window. Dewar and Knight tested a variety of different statistical properties and found that this one was best able to identify the TSZ, when applied to time domain airborne electromagnetic data from a region that encompassed, and extended past, our study area, validating against spatially and temporally co-located water table elevation measurements. Using this approach, the TSZ was estimated to be at 45 mbgs. This estimation is in very good agreement with the water table depth of ∼44 mbgs, obtained from interpolation of the fall 2017 well-based water table measurements. Figure 5 shows the variation of the statistic used to select the TSZ in our study area. In the figure, we see how important it is to constrain this selection by setting a reasonable search window, as without it a depth closer to 30 m would have been selected. The large peak at 30 m in the variation in the selected statistic, suggests that there is a significant change in the subsurface resistivity at this depth, which is likely due to a change lithology rather than saturation state. A depth of 30 m was the depth of the TSZ used by Behroozmand et al. (2019) in the original interpretation of the tTEM data. In that study, a more qualitative approach was used to estimate the depth to the TSZ from observed changes in the resistivity distributions. Using

F I G U R E 5
Plot showing the change in the difference between the 25th and 75th percentiles of the distributions of the resistivity values at each 1-m depth interval in our study area, the search window used to look for the top of the saturated zone, and the depths selected as the top of the saturated zone made when constrained and not constrained within the search window the TSZ estimation of 45 mbgs, we moved forward with analysis of the data above the TSZ.

Development of a resistivity-to-sediment-type transform
As reviewed by Knight and Endres (2005), changes in the resistivity of the subsurface can be the result of changes in water content, pore fluid salinity, and/or sediment type. Without detailed information on the first two factors, it is impossible to determine a single transform from resistivity to sediment type. For this reason, we adopted an approach that acknowledges this uncertainty. Our workflow yielded a distribution of the volume fraction of coarse-dominated sediments for each tTEM-scale resistivity value; the form of the distribution capturing the uncertainty in our transform. This allowed us to produce a 1D model of the most probable fraction of coarse-dominated sediments at every location, called the maximum-likelihood model, and also allowed us to generate many equally viable models displaying the fraction of coarse-dominated sediments.
Using each of the five SBTn logs and their nearest tTEM 1D resistivity model, we developed a transform from resistivity to sediment type applicable to the vadose zone at the site. In creating this transform, we made two key assumptions. First, we assumed that the locations where we acquired SBTn logs (and therefore the locations from which we were using data to create our transform) were representative of the potential variation in lithology, water content, and pore fluid salinity at the study area. Given that the CPT sites were specifically selected to sample the different resistivity zones seen in the tTEM data, and that the study area is small, we deemed this to be a reasonable assumption. The second assumption was that the subsurface structure observed in the SBTn logs was a reasonable representation of the subsurface structure within the volume sampled by the tTEM. This assumption represents one of the most significant challenges in creating the fieldscale resistivity transforms and will be discussed further in the Section 4. These two assumptions, in addition to the inherent uncertainty in the quality of field data, represent uncertainty that is not directly accounted for in the transform. Figure 6 shows the workflow for creating our resistivity transform where, in the two columns, data or derived products are shown as right-justified and green, and steps in the workflow are shown as left-justified and blue. Each component of the workflow is discussed in detail below.
The first step in creating a transform using the CPT data was to convert the SBTs, recorded in the SBTn logs, into descriptions of sediment type using a classification scheme that best fit the goals of this study. We desired a classification scheme that differentiated between sediments in terms of (a) their ability to provide a recharge pathway and (b) their resistivity values (i.e., we wanted to use resistivity to identify zones that were either conducive to, or an impediment to, recharge). The recorded SBTs, for each 0.1-m depth interval in all five logs, were grouped using k-medians clustering, where the data resided in a 10-dimensional space, with nine dimensions being a binary representation of each CPT SBT, and the tenth being the corresponding resistivity values measured during the CPT logging. This analysis suggested two clusters: (a) SBTn Classifications 6 (clean sand to silty sand) and 7 (gravely sand to sand), and (b) all other SBTn classifications. We therefore reclassified the SBTn data along the divide dictated by the clustering analysis by defining a sediment type classification scheme with two categories of sediment type: dominated by coarse-grained sediments, referred to as coarsedominated, for SBTn Classifications 6 and 7; and dominated by fine-grained sediments, referred to as fine-dominated for all other classifications. This data-driven approach resulted in a classification scheme that was conservative in what was designated as coarse-grained material, with only two of the seven SBTn categories classified as coarse dominated. This classification scheme is well suited to the problem of recharge site assessment, where it is preferable to be cautious in designating a region as coarse dominated. Overestimation of the fraction F I G U R E 6 Flowchart illustrating the workflow used to create a resistivity transform. Boxes right-justified and green highlight data or derived products, and boxes left-justified and blue highlight steps in the workflow. Frac CD refers to the volume fraction of coarse-dominated material. tTEM, towed transient electromagnetic; CPT, cone penetrometer testing; SBTn, normalized soil behavior type of coarse-dominated material could lead to slower infiltration than expected, prolonged ponding of water and, in the case of Ag-MAR, damage to crops. The result from this step was five sediment type logs with vertical resolution of 0.1 m.
To determine the mapping between the tTEM resistivity values and sediment type, we first used the bootstrapping approach introduced by Knight et al. (2018) to solve for the range of possible resistivity values for regions of coarse-or fine-dominated sediments in the vadose zone. This approach was selected as it allowed us to capture the range in the resistivity values for each of these sediment types, which in turn allowed us to estimate the uncertainty in our mapping of sediment type.
We first identified the tTEM-derived resistivity model closest to each CPT point. For each depth interval within these resistivity models for which for which there were SBTn sediment type data, we paired the sediment type data and the resistivity value in that interval. In total, this resulted in 51 depth intervals where we had both resistivity and sediment type information. In all cases, the nearest tTEM resistivity model was no more than 5 m away from the corresponding sediment type log. After the derivation from Knight et al. (2018), we assumed that the physics of the tTEM method could be approximated as applying an electric field parallel to the sediment layers, allowing us to describe the bulk resistivity of a package of sediments containing layers of variable resistivity as shown in Equation 1: (1) where ρ tTEM is the resistivity of a package of sediments at the tTEM scale, Frac CD is the volume fraction of coarsedominated sediments within the package of sediments, and ρ CD and ρ FD. are the resistivity of coarse-and fine-dominated material, respectively. Using this approach, we were left with 51 equations with two unknowns: the resistivity of coarse-dominated sediments (ρ CD ), and the resistivity of finedominated sediments (ρ FD ). The Frac CD in these equations spanned the full spectrum from 0 to 100% coarse-dominated sediments. We randomly sampled this system of equations, with replacement, 1,000 times, where the size of each sample set was equal to the parent sample size (51). In each iteration, we used the system of equations to solve for ρ CD and ρ FD . By doing so, we were left with a distribution of 1,000 possible resistivity values for each of the sediment type categories. This range represents uncertainty both in the data itself and in the true variation in resistivity that will exist at a site for any given sediment type category. In the range of resistivity values, we are seeing the impact of the variable sediment type, water content, and pore fluid salinity that is likely occurring within each of our two sediment-type categories. Figure 7 shows the calculated distributions of resistivity values for the two sediment type categories. In this figure, we see higher resistivity values for the coarse-dominated material than for the fine-dominated material, consistent with what we would expect for these materials.
The distributions shown in Figure 7 would allow us to transform each cell in the tTEM 1D resistivity models to one sediment-type category, coarse dominated or fine dominated, if we expected sediment type to be homogeneous at the scale of the model cells. We know, however, that there is most likely to be variation in sediment type below the resolution of the tTEM model cells. We would ideally like to estimate Frac CD from the tTEM resistivity measurement, using the distributions of values for ρ CD and ρ FD we solved for with the bootstrapping method. As a way to relate the tTEM resistivity values to Frac CD , we calculated the tTEM resistivity ρ tTEM of a layered system, using Equation 1, for 0 < Frac CD < 1, increasing Frac CD in increments of 0.01. For each value of Frac CD , we calculated ρ tTEM 10,000 times, randomly sampling the distributions shown in Figure 6 to assign values to ρ CD and ρ FD . This resulted in 1,000,000 pairs of Frac CD and ρ tTEM. Using these, we then created a distribution of all possible Frac CD values for each integer resistivity value, by aggregating all of the calculated pairs within ±0.5 Ω-m of each integer resistivity value.
The results of this can be visualized in two different ways in Figure 8. In Figure 8a, the distribution of Frac CD is shown for every resistivity value, where color corresponds to the count of pairs in each pixel. In Figure 8b, these distributions are plotted for three different resistivity values: 35, 45, and 55 Ω-m. In both of these representations, we see that an increase in resistivity corresponds to an increase in the fraction of coarsedominated material. We also see that for any given resistivity value there can be a broad range of Frac CD values, but there tends to be a peak in these distributions.
There are a number of approaches we could have taken to transforming the tTEM resistivity models to models parameterized in terms of Frac CD , using the derived distributions of Frac CD as a function of resistivity. We could have chosen to make many different realizations of Frac CD for each resistivity model, by repeatedly sampling from the distributions that correspond to the resistivity values in the tTEM model cells; this would have displayed the uncertainty in Frac CD . In this approach, for example, a model cell with a resistivity of 55 Ω-m would have been assigned a Frac CD ranging from 55 to 100%. Alternativity, we could have made a single, maximum-likelihood model, by assigning to each model cell the most probable Frac CD that corresponded to its resistivity value. In this approach, a model cell with a resistivity of 55 Ω-m would have been assigned a Frac CD of 82%. For the purpose of this study, we moved forward with the latter approach, knowing that, if needed, we could have returned to the transform and added more details to our assessment of the variation in sediment type in the study area. Figure 9 shows a visual representation of the most probable Frac CD for each tTEM-scale resistivity value. This was used to make the maximum-likelihood models. This figure shows a gradual transition from ∼100% fine-dominated material near 30 Ω-m to ∼100% coarse-dominated material near 65 Ω-m. In calculating the resistivity values corresponding to FracCD, there were very few values below 20 Ω-m and above 85 Ω-m, giving us little confidence in the Frac CD associated with resistivity values beyond these limits. We therefore interpreted all resistivity values below 20 Ωm and above 85 Ω-m to be 0 and 100% coarse dominated, respectively.

Transform of tTEM data to sediment type
Let us first compare each of the original SBTn logs with the closest maximum-likelihood Frac CD 1D model, obtained through transforming the tTEM resistivity model derived at each tTEM measurement location. This allowed us to evaluate the performance of our transform. Going into this comparison, our expectation was that there would be some differences between the SBTn and the maximum-likelihood Frac CD models, largely due to three key factors. First, small-scale layering captured in the SBTn logs would not be explicitly resolved in the maximum-likelihood Frac CD models, due to the difference in vertical resolution between the CPT and tTEM methods. Second, in locations where there was horizontal variation in the subsurface, the two methods may be sampling different material due to the differences in the horizontal sampling volumes, causing a lack of agreement in the derived sediment types. Last, sharp transitions captured in the SBTn logs would not be resolved in the maximum-likelihood Frac CD models, due to the smoothing inherent in the tTEM inversion approach used in this study. Figure 10 compares each of the SBTn logs to the nearest maximum-likelihood Frac CD model. The separation distance between the two was, in all cases, <5 m. Recall that in setting up our sediment-type classification scheme, for use in developing the transform, the SBTn classes of gravelly sand to sand and clean sand to silty sand were grouped in the coarsedominated category, with all other SBTn classes grouped in the fine-dominated category. The reclassified sediment-type logs were then upscaled to the tTEM vertical resolution prior to creating the resistivity transform. Thus, in reviewing the results displayed in Figure 10, perfect agreement between an SBTn and maximum-likelihood Frac CD model would show the colors orange (clean sand to silty sand) and yellow (gravelly sand to sand) in the SBTn log matching the color yellow (coarse-dominated) in the maximum-likelihood Frac CD model, and all other colors in the SBTn log matching the color blue (fine-dominated) in the maximum-likelihood Frac CD model. In Figure 10, in addition to displaying the maximumlikelihood Frac CD , there is also a white bar showing the interquartile range for the distribution of Frac CD from the transform for each tTEM interval. The maximum-likelihood Frac CD generally falls within this interquartile range; however, for distributions for which the maximum-likelihood is far removed from the median, this is not the case; we see this occurring at low resistivity values. Additionally, an interquartile range was not defined at locations where the resistivity value was either below 20 Ω-m or above 85 Ω-m, the bounds beyond which we assigned coarse-dominated fractions of 0 and 1, respectively.
The success of the transform at capturing large-scale structure can be seen at CPT 1. At a depth of 23 m where the SBTn log shows a transition to a thick package of material that correspond to the fine-dominated category, the Frac CD model also shows 100% fine-dominated. In the depth interval above 23 m, shown in the SBTn logs as composed primarily of material that correspond to the coarse-dominated category, the Frac CD model shows, in most places, high values of Frac CD . There are thin layers of material corresponding to the fine-dominated category, seen in the SBTn log at 10, 15, and 20 m, that have affected the tTEM measurement (seen in the reduction of Frac CD ) but are below the vertical resolution of the measurement. Although these thin layers, which likely would affect the migration of water, are not explicitly resolved, their presence can still be inferred by the reduction of Frac CD at these depths. Another place where differences are seen is the 5-m zone of coarse-dominated material in the SBTn log between 15 and 20 m that appears in the Frac CD model as containing 20-40% fine-dominated material. This is likely due to the influence of the clay layers above and below this unit, as well as regions of fine-dominated material affecting the tTEM measurement but outside the horizontal sampling volume of the SBTn log.
We see similar agreement between the SBTn logs and the Frac CD models at CPT 3, 4, and 5. Keeping in mind the dif-ferences in vertical resolution and horizontal sampling volume, we concluded that the tTEM data were providing useful information about the spatial variation in sediment type. Many of the transitions between fine-dominated and coarsedominated were seen to be gradual (the "staircase structures" in the Frac CD models) due to the smoothing of the tTEM inversion.
There was one location, CPT 2, where the match between the SBTn and Frac CD model was significantly poorer than at the other locations. We interpreted this to be an example of mismatch as a result of the difference in horizontal sampling between the two measurement methods. We further investigated this by considering the resistivity data surrounding the CPT point. Figure 11 shows a roughly north-south running profile connecting CPT 2 and CPT 3, with the nearest tTEM resistivity models projected onto the profile. The models are colored according to their resistivity values, and the locations and extents of the CPT logs are marked in black. In this figure, we see that there is limited horizontal variation in resistivity within the subsurface surrounding CPT 3. In contrast, we see that CPT 2 is located near the edge of a prominent layer at 10-m depth, which becomes thicker and more conductive to the south. Although this feature appears to be continuous and unbroken in the resistivity section, the SBTn log at CPT 2 shows that the upper 20 m consists of many thin interbedded layers of different material, which would have been averaged and smoothed in the tTEM measurement workflow. Combining these two observations, we interpreted the region surrounding CPT 2, around 10-m depth in particular, to be one of potential significant horizontal variation, where a F I G U R E 1 1 Roughly north-south profile along the line connecting cone penetrometer testing (CPT) 2 and CPT 3, with all t towed transient electromagnetic (tTEM) models within 10 m of the profile projected onto the profile. The location and extent of the two CPT datasets are marked in black, and the tTEM resistivity models are colored according to their resistivity F I G U R E 1 2 Two views of the sediment type models resulting from transformation of the inverted towed transient electromagnetic (tTEM) resistivity models, underlain by a satellite image of the study area. The color scale spans from 0% coarse-dominated material (blue) to 100% coarse-dominated material (red). Panel a shows the full extent of these transformed models, between the ground surface and the top of the saturated zone. Panel b shows the transformed models from 10 m below ground surface to the top of the saturated zone transition in dominant sediment type is occurring. We suspect that this is the reason the SBTn log, which reflects the subsurface just at that point, is so different from the Frac CD model, which is reflective of a larger area surrounding the CPT point.
Our comparison of the SBT n logs with the nearest maximum-likelihood Frac CD models showed that, within the expected limitations, we could use our transform to create models of sediment type from the tTEM data, which captured the large-scale structure in the subsurface. Although this gave us confidence to move forward with applying our transform to the full dataset, there is also an opportunity for future research to explore the impact of using resistivity models derived through the use of different inversion approaches. The final step in this study was to explore ways to use the maximum-likelihood Frac CD models to map out pathways with a high fraction of coarse-dominated material, from the surface to the TSZ, for the purpose of recharge site characterization and assessment.

4.2
Evaluation of recharge pathways in sediment type models Figure 12 shows two views of the derived sediment type models resulting from transformation of the tTEM resistivity models to maximum-likelihood Frac CD models. Figure 12a shows the models starting from the ground surface and extending down to the TSZ, whereas Figure 12b shows the models from 10 mbgs to the TSZ. The color scale spans from 0% coarse-dominated material (blue) to 100% coarse-dominated material (red). These show the models from the same vantage as the resistivity models shown in Figure 2. In this figure, we see that the sediment type mapped with this approach is highly variable, both laterally and with depth. Although much of the surface is covered in coarse-dominated material (Figure 12a), just 10 m below the surface (Figure 12b), there appear to be channels of coarse-dominated material surrounded by laterally extensive regions of fine-dominated material. This figure demonstrations that there can be F I G U R E 1 3 A view of the sediment type models resulting from transformation of the inverted towed transient electromagnetic (tTEM) resistivity models sliced halfway up the almond grove, underlain by a satellite image of the study area. The color scale spans from 0% coarse-dominated material (blue) to 100% coarse-dominated material (red). Markers 1-3 point to regions of interest discussed in the text. Dashed white line shows permeable pathway from the surface to the top of the saturated zone significant spatial heterogeneity in sediment type, the property controlling the movement of water in the subsurface, and that the approach outlined in this study can be used to reveal this heterogeneity.
With this 3D dataset, it is possible to scroll through the sediment type models to determine whether specific regions within the field have connected pathways that would allow recharge water to reach the TSZ. Figure 13 shows a portion of the models, with those from the southern half of the field removed. Markers 1-3 on the figure note three locations underlain with very different sediment type models. Let us first consider what is present immediately below the three markers. At Marker 1, fine-dominated material extends from the surface all the way to the TSZ. At Marker 2, there is a thick package of coarse-dominated material starting at the surface, and then a sharp transition to fine-dominated material at 30 mbgs. At Marker 3, there is a layer of fine-dominated material starting at the surface, followed by a layer of coarsedominated material and mixed material extending down to the TSZ. Were 1D logs with information about sediment type to be collected, either with a drilling or CPT system, at the locations of these markers, there would be no evidence of a direct connected pathway of coarse-dominated material from the surface to the TSZ, and it would be impossible to determine if there was an indirect pathway. Because of the spatial density of the tTEM measurements, however, we can clearly see a large, connected pathway (noted with the white dashed line) of coarse-dominated material originating below Marker 2, and connecting with the TSZ below Marker 3.
As demonstrated through Figures 10-13, our transform methodology, when applied to tTEM data, was able to produce models of subsurface sediment type that honored the structure seen in the CPT logs and was useful for mapping out pathways of coarse-dominated material from the ground surface to the TSZ.

Assessment of the impact of the number of control points on the development of the resistivity transform
Developed and demonstrated in this study is a methodology for the construction of a site-and measurement-specific resistivity-to-sediment-type transform that can be used to generate sediment type models from tTEM data. The construction of this transform required an independent measure of sediment type co-located with tTEM data. In this study, we were fortunate to have five SBTn logs within our small site. Additionally, we had the tTEM resistivity models to use when siting the CPT points, allowing us to guarantee that we sampled from locations representing a range of different subsurface conditions. Given that datasets of this density and quality may not be commonly available at other sites, we investigated the potential impacts of not having as much data available when constructing a resistivity-to-sediment-type transform. Specifically, we explored the potential impact of having fewer SBTn logs with which to create a transform. We repeated the process described in this study for constructing the transform using only two, three, or four SBTn logs. In each case, we tested every possible combination of SBTn logs, sampling from the set of five logs from our site, to see how much the transform might vary depending on which specific logs were used. This resulted in 10 transforms created using only two logs, 10 transforms created using only three logs, and five transforms created using only four logs. Figure 14 compares these different transforms, displayed as the most probable Frac CD for each resistivity value, to the transform obtained in this study using all five SBTn logs.
In Figure 14, we can see that reducing the number of SBTn logs used to create the transform can result in very different transforms, depending on which combination of SBTn logs were used. Although the use of four logs resulted in Using Four CPT pts F I G U R E 1 4 Comparison of the most probable volume fraction of coarse-dominated material (Frac CD ) for each towed transient electromagnetic (tTEM) resistivity found using only two (blue), three (green), or four (red) normalized soil behavior type (SBTn) logs, and the relationship established using five logs (black). CPT, cone penetrometer testing reasonably close clustering of the estimates from the constructed transforms, the reduction to three and then two logs resulted in large differences. These differences would then be propagated into the derived maximum-likelihood Frac CD models, potentially resulting in very different sediment type models, depending on the tTEM resistivity values being transformed. To assess the impact of this, we took the two end members from each set of transforms, defined as the ones that would, on average, predict the highest and lowest values of Frac CD , and used them to transform all of the tTEM resistivity models above the TSZ to maximum-likelihood Frac CD models.
In the maximum-likelihood Frac CD models created using the five-log transform, presented in this study, an estimated 34% of the vadose zone is coarse dominated. The Frac CD models, created using the end members of the two-log transform, contain estimates of 21 and 50% coarse-dominated sediments in the vadose zone. This large range shows that, given the complexity of the spatial variation in sediment type in our study area, it may be risky to use a two-log transform; there would be very different predictions of Frac CD in the subsurface depending on which two SBTn logs were used. Whether the risk is acceptable is tied to the specific questions being answered, or decisions being made.
When we performed the same exercise using the end members of the three-log transform, the bounding estimates of coarse-dominated sediments are revised to 22 and 45%, respectively. With the four-log transform, these are narrowed to 32 and 45%, respectively; predictions of Frac CD are much closer to the five-log transform than obtained with the use of two SBTn logs. Although it would be difficult to know, a priori, the number of logs needed to construct a reliable transform, this analysis serves to demonstrate how important it is to have the data used in developing the transform be a sufficiently representative sample of the subsurface.

CONCLUSIONS
Within California's Central Valley, there is increasing interest in, and a need for, MAR projects, as water managers seek to balance groundwater supply and demand. In many cases, the ability to characterize a site in terms of the presence, or absence, of good pathways for recharge from the ground surface to the water table would allow managers to minimize costs and risks associated with MAR projects. The tTEM method can be effective for assessing potential recharge sites, but if a high level of confidence is required in the derived maps of sediment type, the uncertainty associated with tTEMderived models of sediment type needs to be reduced. We have developed a methodology, integrating the acquisition of tTEM and CPT data, which can be used to construct a site-and measurement-specific transform between resistivity and sediment type. We found that eCPT, though useful in interpreting SBTn logs, was not appropriate for direct use in developing a transform for tTEM data given the differences in the measurement workflow between these different datasets. It was the SBTn logs, acquired using CPT, that provided an ideal data set with which to develop the transform required to map tTEM resistivity models into models of sediment type. We suggest that CPT be adopted, where viable, as the optimal method to obtain point-location measurements of sediment type. The SBTn logs can be categorized in a way that defines two sediment types, coarse-grain-dominated (sands and gravels) and fine-grain-dominated (silt and clay), that (a) differ in terms of hydraulic properties and (b) can be differentiated based on their electrical resistivity values. We found it was important to have acquired the CPT data after the acquisition of the tTEM data, so that we were able to site CPT measurements that strategically sampled locations representing the range of resistivity values and spatial heterogeneity at the site. This was critical, as insufficient sampling could lead