Managed aquifer recharge site assessment with electromagnetic imaging: Identification of recharge flow paths

Surface spreading recharge, the intentional flooding of the ground surface to replenish a groundwater system, is one approach used to mitigate groundwater overdraft in California's Central Valley (CV). Choosing appropriate sites for surface spreading recharge, in regard to the sites’ ability to convey water from the ground surface to the desired recharge depth, can be challenging because of limited knowledge of the properties of the subsurface. In this study, we present an approach for using a towed time‐domain electromagnetic (tTEM) imaging method to develop three‐dimensional (3D) models of sediment type, map potential flow paths through the subsurface, and evaluate sites for surface spreading. We began with tTEM data from seven sites in the CV along with an existing resistivity‐to‐sediment type transform. We leveraged geostatistical methods to generate multiple 3D models of binary (flow and no‐flow) sediment type from the tTEM data. We then developed two metrics to assess the quality of sites for recharge: (a) the depth to the shallowest no‐flow unit beneath each point at the surface and (b) preferential flow paths lengths measured by finding the shortest distance through connected flow units between surface points and the desired depth. We explored how these metrics can be used to identify optimal areas within a site, then developed a way to compare and assess the relative suitability of each site using the decay in the number of vertical flow paths as a function of depth. Our methods can be used to rapidly identify potential sites for surface spreading recharge.

passage of the Sustainable Groundwater Management Act (2014a, 2014b, 2014c), there is now considerable interest in finding ways to supplement natural recharge processes through managed aquifer recharge.
A detailed overview of various recharge methods can be found in Bouwer (2002). One of these methods, which we refer to as surface spreading recharge, involves the use of water to flood an available area on the ground surface, such as an agricultural field during nongrowing seasons, open land that has been converted to a recharge basin, or a natural or restored floodplain. The objective is to have the applied water percolate through the subsurface to the water table, which marks the top of aquifer and intended storage location, thus replenishing the groundwater system and storing groundwater for future use in dry years. However, this method involves certain risks that are related to the conditions and hydrogeologic properties of the selected site: 1. Slow infiltration of the water could lead to ponding at the surface or lingering of water in the root zone, both of which may weaken the roots and bases of trees, leaving them susceptible to damage in strong winds (Dahlke et al., 2018;O'Geen et al., 2015;Ganot & Dahlke, 2021). 2. Slow infiltration and subsequent ponding can lead to increased evaporative loss. 3. If the spatial distribution of the sediment types favors lateral over vertical movement, the water may be unable to reach the intended storage location. If the goal is to store groundwater for future use or reduce subsidence, it is important that the water reaches the intended location. Furthermore, barriers inhibiting downward flow can result in perched water tables or water mounds that will slow infiltration rates (Bouwer, 2002). 4. When contaminants are a concern, either as result of agricultural practices or the chemistry of the recharge water, recharge can have a negative effect on water quality (Bouwer, 2002).
The first risk applies to agricultural fields, whereas the latter three apply to any site under consideration for surface spreading recharge. Given the above, one of the key questions that needs to be addressed before starting the design and operation of a surface spreading project is where to recharge to reduce these risks. In this study, our focus was on the identification and assessment of recharge flow paths, with the objective being to maximize the efficiency with which water can move from the surface to the water table. Related to this, Bouwer (2002) described two major criteria for the hydrogeologic properties of the vadose zone that should be satisfied: (a) fast infiltration in the shallow subsurface and (b) downward, uninhibited flow connections to the intended storage location. In this study, we developed a methodology to aid in the selection of a recharge area-an entire site, or part

Core Ideas
• A geophysical method, tTEM, can provide critical information about potential recharge sites. • The tTEM resistivity can be transformed to sediment type using cone penetrometer testing. • An algorithm designed to find connected coarsegrained material can locate recharge pathways. • A scheme for selection of recharge sites uses the length of recharge pathways as a key metric. • Multiple models of sediment type reveal uncertainty in site assessment for recharge.
of a site-that satisfies these criteria. The effect of surface spreading on water quality will also influence the selection of a recharge area. While our methodology can provide information relevant for assessing potential changes in water quality, far more information about subsurface conditions and properties would be required. Adequately addressing water quality issues is thus beyond the scope of the developed methodology. O'Geen et al. (2015) developed a map of the recharge suitability of agricultural regions of the CV using the Soil Agricultural Groundwater Banking Index (SAGBI), which is a scoring system used to rank regional areas for surface spreading recharge by assessing the hydrogeologic properties of the soils along with other surface features, such as topography and the chemicals used in farming. Two of the five factors used in the scoring, root-zone residence time and deep percolation, were designed to predict where water would quickly infiltrate the surface soils and move past the root zone to deeper sediments, respectively. These two factors are determined using estimates of the hydraulic conductivity (K) of the soil horizons specific to the region, where high K is better for recharge. These soil horizons and assigned K estimates are obtained from the USDA-NRCS digital soil survey (USDA-NRCS, 2014) and, in general, are only valid estimates up to the first 2 m below the ground surface. While the SAGBI has great utility for assessing the first 1 or 2 m, root zones may extend deeper (up to 4 m) and water tables in the CV are often tens of meters below the surface, so knowledge of the sediment types at greater depths within the vadose zone is important for predicting where there are uninhibited flow connections in the shallow subsurface and to the aquifer. An additional limitation in the use of the SAGBI to assess specific sites is the fact that it was developed from soil surveys consisting of map units ranging in size from 0.02 to 2 km 2 (5-500 acres), and while many units contain multiple soil profiles, the SAGBI used only the most abundant type in each map unit. Therefore, the K estimates may not best represent the local conditions within any individual site, which are often on the order of <0.5 km 2 (tens of acres). In fact, the USDA-NRCS explicitly recommends onsite investigations when planning the use of a small, local area (USDA-NRCS, n.d.).
One approach to assessing the suitability of a site for surface spreading recharge would be to conduct the modeling required to simulate flow in the vadose zone. Such modeling can be complex and time-consuming and requires the accurate estimate of parameters-infiltration rates and hydraulic conductivity-that are challenging to obtain, especially in the absence of nearby wells. However, after the initial saturation, flow patterns become similar to those observed in the saturated zone, where water travels in preferential flow paths that are dominated by networks of sediments with high K, while occasionally "jumping" through smaller regions of low K (e.g., Bianchi et al., 2011;Rizzo & de Barros, 2017). These findings suggest that, upon initial wetting of the vadose zone, regions dominated by coarse-grained sediments (sand and gravel) that have high K are the conduits for flow, whereas areas dominated by fine-grained sediments (silt and clay) that have low K inhibit flow. If it is possible to resolve the spatial distribution of sediment type in terms of the amount of coarsegrained material, then it may be feasible to predict where there will be predominantly vertical preferential flow paths both in the first few meters for fast infiltration and connecting the ground surface to the intended storage location, using sediment type as a proxy for K. What this requires is a model of sediment type with sufficiently high spatial coverage and resolution to evaluate the suitability of sites for surface spreading recharge.
The CV hydrologic model (CVHM) (Faunt, 2009) provides information about the unconsolidated sediment type in the CV. The CVHM includes an estimated fraction of coarse-grained material, which we refer to as coarse fraction, throughout the CV, extending to a maximum depth of ∼800 m with cells of lateral dimension 1,600 by 1,600 m and vertical dimension of 15 m. Because of its widespread spatial coverage, Jankowski et al. (2018) used the CVHM to augment SAGBI with a factor for scoring deep sediment type, where deep refers to the zone between the soil (described with SAGBI) and the water table. At each CVHM cell, they calculated the harmonic mean of coarse fraction in this zone, with higher fractions receiving better scores in the assessment. With a recharge site typically smaller than a single cell in the CVHM, all that would be obtained from CVHM is an average value of coarse fraction in 15-m-depth intervals. In order to identify recharge pathways from the surface to the water table, we need to resolve a finer-scale spatial distribution of coarse-grained material, that is, on the scale of a few meters, rather than kilometers, per cell.
Over the past 10 yr, there has been growing interest in the use of geophysical methods as a cost-efficient way of characterizing sediment type to support managed aquifer recharge operations, including surface spreading recharge. Studies have used a variety of geophysical methods including surface and borehole electrical resistivity tomography, time domain electromagnetics, and nuclear magnetic resonance (Cook et al., 1992;Gottschalk et al., 2017;Maliva et al., 2015;Mawer et al., 2013;Parsekian et al., 2014;Sendrós et al., 2020;Lawrie et al., 2013). Two of the most significant impediments with the use of many of these methods are the time it takes to acquire data and the challenge of efficiently and accurately deriving the information needed about the subsurface to assess the suitability of a site for recharge.
In our research, we are using a new geophysical system, towed time-domain electromagnetics (tTEM), which is a time-domain electromagnetic system  towed behind an all-terrain vehicle that has been found to be a highly effective approach for subsurface characterization at the scale of an agricultural field. The tTEM system contains a pair of rectangular transmitter and receiver loops of size 2 by 4 m and 1 by 1 m, respectively, spaced with a 9-m offset. Time-varying currents are injected into the transmitter loop; this induces eddy currents in the subsurface with the resulting signals measured at the receiver loop. Through inversion, a model of the electrical resistivity of the subsurface is recovered from the measured signals.
The tTEM system was used at seven sites (four recharge basins, two nut tree orchards, and one open field) in the Tulare Irrigation District (TID) in the San Joaquin Valley of California (Behroozmand et al., 2019). The relatively small size of the tTEM system meant that data were easily acquired by driving the all-terrain vehicle between rows of trees in almond [Prunus dulcis (Mill.) D. A. Webb] and pistachio (Pistacia vera L.) orchards. Recovered through inversion of the acquired data were one-dimensional (1D) models of the electrical resistivity of the subsurface. Horizontal spacing of the 1D resistivity models was ∼10 m along each acquisition line, with lines separated by 10-20 m depending on the site. Although the 1D vertical resistivity models are continuous in depth, the noncontinuous horizontal spacing between tTEM measurements results in sparsely filled, three-dimensional (3D) resistivity models when all are combined to represent the full site. Vertical discretization of the resistivity models begins at ∼1 m at the surface and gradually increases to 6-12 m by the depth of investigation, which was between 50 and 70 m; the depth of investigation is considered to be the depth to which the estimates of electrical resistivity are reliable and depends on the subsurface electrical resistivity. Using measurements of sediment behavior, acquired using cone-penetrometer testing, the resistivity data in the vadose zone from one of the sites were transformed to the fraction of coarse-grain-dominated material, referred to as coarse-dominated material (i.e., material consisting mainly of sands and gravels), by empirically determining a resistivity-to-sediment-type transform (Goebel & Knight, 2021).

F I G U R E 1
Model of the fraction of coarse-dominated material at a site in the Tulare Irrigation District. The preferential flow path depicted by the dashed white line was manually selected using the spatial distribution of units with a high fraction of coarse-dominated material. Figure modified from Goebel and Knight (2021) In the study by Goebel and Knight (2021), the sparse 3D model displaying the spatial variation in the fraction of coarse-dominated material was qualitatively interpreted to determine where there appeared to be preferential flow paths from the surface to the water table at a site. An example is shown in Figure 1 (modified from Goebel and Knight [2021]), where the authors identified a potential connected pathway of coarse-dominated material through the subsurface. This qualitative interpretation required manually scrolling through the data to identify these pathways and relied on visual interpolation of sediment type in vertical sections lacking tTEM data. While valuable for illustrative purposes, this approach did not provide a way to easily locate or quantify the length of these pathways nor did it consider the uncertainty associated with identifying continuous pathways in a sparsely filled model.
In this study, we advanced the use of tTEM data for assessing potential recharge sites by developing and demonstrating a more quantitative methodology, most of which can be automated and efficiently used, for identifying preferential flow paths that connect the surface and the water table and quantifying their length. We assessed the quality of the flow paths by considering first, their ability to rapidly transmit water beyond the root zone and, second, their distance to reach the water table. In this context, shorter preferential flow paths with no shallow low-K barriers were deemed to be of higher quality, as these will transmit water most efficiently to the intended storage location. This method allowed us to identify areas within a site on which recharge efforts should be concentrated and evaluate the relative suitability of potential sites for surface spreading recharge. While a binary system (high-K and low-K) is a simplified representation of the subsurface, it is a useful model for rapidly identifying preferred sites based on our metrics. These sites can then be selected for further tests (e.g., infiltration rates) rather than performing these time-consuming and often expensive assessments at all available sites.
Our methodology of finding high-quality preferential flow paths consisted of two parts. We began with resistivity models derived from the inversion of tTEM data and used a geostatistical simulation to transform them into many 3D models of sediment type, classified as either high-K or low-K, which could represent the subsurface (described in Methods I section). By developing numerous models with varied input parameters for each simulation run we were able to capture the uncertainty resulting from the assumptions we needed to make in their development.
We then assessed the quality of preferential flow paths in the produced models by locating those that had no shallow barriers and provided the shortest connections of high-K units between the surface and the water table (described in Methods II section). This allowed us to evaluate the preferential flow paths at each site and compare the different sites to each other.

DATA AND STUDY AREA
The tTEM data used in this study were acquired in October 2017 by Behroozmand et al. (2019) at seven sites in the TID (Figure 2). Four of the sites-Swall-North, Swall-South, Martin-West, and Martin-East-have been used as active recharge basins. Another cluster of two sites consists of an almond grove, Almond-West, and a smaller field across the street, Almond-East, which is under consideration for surface spreading recharge. The final site is a large pistachio orchard given the designation of Nichols. The land surfaces at all of the sites are close to level in terms of elevation.  (Faunt, 2009). The sites are consistently smaller than the CVHM cell size We defined the depth to the water table at each site using the depth-to-groundwater map created by the TID (2017). At sites Almond-West and Almond-East, the depth obtained from this map agreed well with a detailed analysis conducted by Goebel and Knight (2021) to identify the water table depth using the tTEM data.
A regional assessment of our study area indicates that the top 30 m of the subsurface consists predominantly of young, continental alluvium distributed in interbedded deposits of clay, silt, silty sand, and gravelly sand (Mid-Kaweah Groundwater Sustainability Agency, 2019). Below these deposits is a transition zone to less-permeable, clay-rich sediments that would greatly inhibit the flow of recharge water. Thus, there is significant spatial heterogeneity in sediments that needs to be assessed in order to characterize the recharge potential of sites in this area.

MATERIALS AND METHODS
The methods we developed can be split into two parts: (a) Methods I, the development of multiple models of subsurface sediment type and (b) Methods II, the assessment of subsurface models for surface spreading. The full Methods I workflow is visually depicted in Supplemental Figure S1.

3.1.1
Recovering uniformly gridded 3D resistivity models from tTEM data We began this study with the resistivity models from the study by Behroozmand et al. (2019). The resistivity model at each site, recovered through inversion of the tTEM data, is nonunique, meaning that many resistivity models could fit the data equally well. We accounted for the uncertainty in the inversion by generating 100 resistivity models at each site using a posterior sampling (PS) technique (Kang et al., 2021;Rue, 2001), which used the recovered resistivity model as an input. This technique uses a linear approximation and assumes a Gaussian posterior distribution. We refer to these 100 models as the PS resistivity models.
Next, we resampled these PS resistivity models onto uniform 5-× 5-× 1-m 3D grids, leaving a 25-m buffer zone extending beyond the bounds of the resistivity models, which were contained within the edges of the site. This resulted in 100 sparse uniform grids of resistivity-one for each of the 100 PS resistivity models. We assigned resistivity values only to those voxels for which there was a resistivity value within 5 m during the resampling, leaving all other voxels, including the buffer zone, empty.

F I G U R E 3
Transform from resistivity to categorical probabilities. We started with (a) the transform from resistivity to fraction coarse-grain-dominated material (FRAC CD ) as obtained in Goebel and Knight (2021). From these distributions, we derived (b) the transform to categorical probabilities by comparing the number of counts at or above the high-K threshold [white dashed line in (a)] with the total counts at each resistivity bin. Resistivity values within the window of uncertainty in (b) may represent either a high-K or a low-K unit. We applied the transform in (b) to each resistivity value 3.1.2 Transforming resistivity to probabilities of sediment type The next step for each site was to transform the uniformly gridded resistivity models to subsurface models in which each voxel containing a resistivity value was assigned a probability of being high-K material (P high-K ). We began with the transform developed by Goebel and Knight (2021), which defines distributions of the fraction of coarse-dominated material (FRAC CD ) that correspond to integer resistivity values as shown in Figure 3a. While this relationship was developed using sediment samples at the Almond-West site, we assumed that it was applicable to the other six sites where the sediments are similar in terms of composition and depositional environment. If a highly accurate subsurface model had been required, the procedures followed by Goebel and Knight (2021), using cone-penetrometer testing, could have been repeated at a site so as to develop a site-specific transform. In order to define P high-K , we needed to select a cutoff between values of FRAC CD that are low-K vs. high-K. Setting an appropriate cutoff should take advantage of any information about the properties of the site and will be dependent on the objective of the investigation and the risk being managed. The number of voxels identified as high-K in the workflow would increase as the FRAC CD value selected for the cutoff is decreased. If it is important not to overestimate the quality of preferential flow paths, so as to reduce the risk of surface ponding, then higher values of FRAC CD should be used. Alternatively, to reduce the risk of overlooking a suit-able site by underestimating the quality of preferential flow paths, lower values of FRAC CD should be used.
Lacking information about the way in which K varied with FRAC CD at any of the sites, we referred to the work of Fogg et al. (2000), which showed that a volume of sediments in the CV requires a minimum of 18% coarse-grained material to have a high-K pathway for flow. We thus set FRAC CD = 0.18 as the cutoff in our study area noting that this value may be readily adjusted to fit the needs of the study.
After selecting the cutoff between high-K and low-K, we transformed all of the resistivity values to P high-K using the probability vs. resistivity relationships shown in Figure 3b; these were extracted from the FRAC CD distributions in Figure 3a by relating the number of counts above the FRAC CD cutoff to the total number of counts in each resistivity bin. This step resulted in 100 sparse P high-K models, which we reduced to a single mean P high-K model by calculating the arithmetic mean value at each populated cell.

3.1.3
Generating fully populated models of binary sediment type The final stage in our workflow was to generate fully populated models of binary-sediment-type at each site, which could be used to find the preferential flow paths. We sought to generate many binary-sediment-type models to capture the uncertainty resulting from varying the assumptions about spatial correlation of sediment type. To accomplish this, we used a geostatistical module known as the co-sequential indicator simulation (co-sisim) in Stanford's Geostatistical Modeling Software (SGeMS) (Remy et al., 2009), which outputs fully populated categorical models (e.g., of binary-sediment-type) using sparse input data and various spatial parameters that describe the expected spatial correlation and proportion of the output categories. In our workflow, we used the sparse mean P high-K model as the sparse input data; for the spatial parameters, we defined (a) the expected total (global) proportion of high-K units in the fully populated model (G high-K ) and (b) a variogram model, which designates the spatial correlation (range) of sediment type. We characterized the variogram model by its mutually orthogonal maximum horizontal, minimum horizontal, and vertical ranges, and the angle between the axis of the maximum horizontal range and north. The parameters we used to define the model grid for each site (G high-K , the experimental variogram model and depth to water table) are listed in Supplemental Table S1.
To generate multiple models, we ran co-sisim 200 separate times. During each run, we used the same sparse input data but varied the values of the spatial parameters, as we were able to estimate them only roughly from the available data. After first estimating the five spatial parameters, we developed uniform distributions for each. These distributions allowed us to randomly select values for these parameters during each simulation run, resulting in 200 different models of binarysediment-type. These steps, founded in Monte Carlo principles, follow the same basic structure as other geostatistical modeling studies that attempt to account for the uncertainty in simulation parameters (e.g., Soutter & Musy, 1998;Christensen, 2004;Caers & Hoffman, 2006;Wang et al., 2020). Supplemental Text S1 describes in detail the steps we used to obtain the initial estimates of G high-K and the variogram model, develop the distributions for each spatial parameter, and randomly generate a parameter set during each run of co-sisim.

Methods II: Assessment of subsurface models for surface spreading recharge
The binary-sediment-type models at each of our seven sites provided the data needed to evaluate preferential flow paths and identify optimal areas for surface spreading recharge within each site. At each site, we evaluated each of the 200 models using two metrics, described below, and selected the optimal rectangular recharge areas based on the arithmetic mean of the outputs. We chose rectangular areas for this study for practical reasons; it is easier to select and analyze rectangular areas on gridded data, and if only a portion of a site were selected for recharge, barriers could be positioned to follow the row crops, which are gridded as well.

3.2.1
Defining metrics and parameters for selection of recharge areas In order to avoid ponding and root damage, water must infiltrate past the root zone as quickly as possible. While we cannot accurately quantify infiltration rates with our models, we can identify regions of the site that have high-K units beyond the root zone, thus serving as indicative of fast infiltration rates because the presence of low-K units slows infiltration. Consequently, we developed Metric 1 to quantify the depth to the shallowest low-K unit beneath each surface pixel, here referred to as DEPTH low-K and demonstrated in Figure 4a. Values for DEPTH low-K will range from 0 m to the depth of the water table, where a larger value is better for recharge.
In addition to Metric 1, we needed a metric that identified flow paths with short connections to the water table. We developed Metric 2 to quantify the lengths of preferential flow paths, which we defined as the routes through high-K units that minimized the distance between points at the surface and the water table. A two-dimensional schematic showing examples of preferential flow paths (which will extend into 3D) is given in Figure 4b. An ideal recharge area will have short preferential flow paths where the minimum length would be equal to the depth to the water table. It is important to note that these paths are not necessarily the routes through which water would flow, rather, we considered the path length to be a proxy for the time to reach the water table (i.e., shorter is faster).
Quantifying both metrics allowed us to identify optimal areas for recharge within an individual site. We first used Metric 1 to locate areas with deep DEPTH low-K , then used Metric 2 to additionally prioritize regions with short preferential flow paths. We chose thresholds that defined suitably "deep" and "short," and then selected a rectangular optimal area that satisfied both thresholds. The DEPTH low-K threshold should be chosen by considering the depth of the root zone for the crops at the sites and should err on the side of being greater than deemed necessary if there is a high level of concern about ponding and root damage. We elected to use DEPTH low-K threshold as 8 m, which is twice the maximum depth of the root zone for the almond and pistachio trees in our study area. In practice, the choice of the DEPTH low-K threshold will also depend on the volume of recharge water, the surface area over which the water is applied, and the infiltration rates.
To select a threshold for Metric 2, we first defined a normalized path length (NPL) as a preferential flow path length divided by the local depth to water table. The NPL is thus a measure of the tortuosity of the path where a value of 1 is the minimum and deemed "best." In this study, we selected a NPL threshold of 2 as the cutoff for short paths. A higher NPL corresponds to predominantly lateral flow, thus increased residence time, which may result in perched water tables that F I G U R E 4 Conceptual illustration of metrics used to assess the models of fully populated binary-sediment-type. (a) Metric 1, the depth to the shallowest low-K unit beneath each pixel (DEPTH low-K ). We consider only pixels with DEPTH low-K exceeding the selected threshold (gray horizontal line) as suitable for recharge. (b) The normalized path lengths (NPL) from each pixel to the water table. For both metrics, we only analyzed surface pixels within the bounds of the site noting that preferential paths may extend beyond these bounds once they infiltrate past the surface. Though pixels eight and nine have the deepest DEPTH low-K , there are no connections to the water table, thus they have the longest NPL. In our model, pixels are 5 by 5 by 1 m (illustration is not to scale) ultimately slow infiltration rates (Bouwer, 2002). We note that the tolerance for lateral flow will depend on each individual case.

Calculation of metric maps
For each model of binary-sediment-type at a site, we calculated both metrics at every 5-× 5-m surface pixel within the lateral bound of the site (25 m inward from the edges of the 3D grid), displaying the result in plan-view maps. This resulted in 200 maps per metric. We obtained a single representative map for each metric by calculating the arithmetic mean value at each pixel across the 200 maps. Metric 1 involved a simple calculation of the depth to the shallowest low-K unit in each vertical column. For Metric 2, we calculated the length of the preferential path lengths using the Python package scikit-fmm (Furtney, 2015), referred to here as the "eikonal solver," which uses the fast-marching method (Sethian, 1996) to solve the eikonal equation, an approximation of the wave equation. This equation is a nonlinear partial differential that can be used to find the shortest distance (or travel time) between points and surfaces. When calculating distance, as done in our workflow, the simplifying assumption made was that there were only flow-and no-flow units. In our models of binary-sediment-type, high-K units were defined as flow units and low-K units were defined as no-flow units.
For each binary-sediment-type model, we configured the eikonal solver to find the shortest distance between each pixel at the surface to the water table. Once a path went below the surface, it could extend beyond the bounds of the site and into the 25-m buffer zone. Additional constraints included no-flow boundary conditions on the vertical sides of the 3D grid and an added layer of flow units at the surface to mimic the ability of water to move laterally across the surface of the site. We normalized the resulting maps of preferential flow path lengths to display the NPL at each pixel.
We used the representative metric maps to identify pixels with mean DEPTH low-K ≥ 8 m and mean NPL ≤ 2. We refer to the contours that contain these pixels as the 8-m-DEPTH low-K contour and 2-NPL contour, respectively, and overlaid both on the site's representative maps. At each site, we used these maps to manually select an optimal rectangular area for recharge. These optimal areas could range in size from a small rectangular area to the entire site. We sought to select optimal areas within which all metric values would have mean DEPTH low-K ≥ 8 m and mean NPL ≤ 2. In practice, some optimal areas included values outside these thresholds because of the challenge of selecting a rectangular area with irregularly distributed metrics.

Comparison of areas within a site
Identifying an optimal area for recharge can involve assessment of a single site to select an area within the site. This might be done, for example, if there was low tolerance for the risk of ponding (e.g., an agricultural field). A selected area could be separated with flood barricades from the rest of the site.
We demonstrate the use of the two metrics for selection of the optimal area for surface spreading within the Almond-West site. Figure 5 shows the (a) mean DEPTH low-K and (b) mean NPL maps calculated from the 200 model results and our manually selected optimal recharge areas. It is clear that Almond-West displays high variability in DEPTH low-K and NPL, suggesting significant spatial heterogeneity and underscoring the importance of resolving sediment type at high resolution. The CVHM texture model, for example, can only show that the entire almond grove has a harmonic mean value of percentage coarse-grained material of 58% between the Vadose Zone Journal F I G U R E 5 Metric maps for Almond-West site. (a) Mean depth to the first low-K unit (DEPTH low-K ) and (b) mean normalized path lengths (NPL) are the average metric results calculated from the 200 generated models surface to a depth of 45 m. The representation we get from tTEM is a much more informative view of sediment variation and how that may affect recharge efforts.
In selecting the optimal area for surface spreading at Almond-West, we first consider Metric 1. The mean DEPTH low-K map indicates that only ∼50% of the site may have sufficiently deep low-K units for fast infiltration during recharge, with areas outside of the 8-m-DEPTH low-K contour considered to be at risk of slow infiltration and subsequent ponding. When we then consider Metric 2, the southern half of the field has no NPL ≤ 2, despite having sufficiently deep DEPTH low-K . Only the northern half of the site has any areas with NPL ≤ 2. An additional consideration is the small size of the two eastern-most areas. This suggests that a high number of the preferential flow paths merge along these routes, which might affect infiltration rates. We used the overlap between the two metrics to delineate the optimal area shown in Figure 5, which largely satisfies both recharge criteria.
Similar figures for each of the other sites can be found in Supplemental Figures S2-S7. When there were no, or extremely small, overlapping contours, as is the case for Martin-West and Martin East (Supplemental Figures S6 and  S7), we prioritized DEPTH low-K while attempting to include the local minimum of NPL, noting that these sites would not be ideal for recharge under our criteria.

Comparison of multiple sites
We now consider the case where a number of sites are under consideration for recharge. Our approach can also be used to obtain a relative ranking of the sites, but it is difficult to draw any conclusions based on the metric maps alone because of variations in the size of sites and water table depths, except for scenarios in which we can eliminate sites from consideration because of little-to-no overlap in the metric contours (Martin-West and Martin-East; Supplemental Figures S6 and S7, respectively). We developed a way of quantitatively comparing the sites by calculating the number and percentage of vertical high-K paths as a function of depth from the ground surface to the water table. The total number of possible vertical paths is the number of pixels in the DEPTH low-K map, where the value at each pixel denotes the expected depth at which the first low-K unit will be encountered, indicating that the path vertically connecting the surface to the pixel's DEPTH low-K consists only of high-K units. Moving from the surface to the local water table at 1-m depth slices, we defined the number of vertical paths at each depth as the number of pixels that had not yet exceeded their DEPTH low-K value. The same exercise can be done for any area within a site, where the total number of possible vertical paths is the number of pixels within the selected area. A decay in the number of vertical paths as a function of depth indicated that vertical flow paths initiating from multiple points at the surface were merging. Thus, a low rate of decay indicated a better site. Using both the number and percentage of vertical paths, where the percentage is the number of vertical paths at some depth normalized by the number at the surface, allowed us to identify, respectively, sites with (a) the largest areas well-suited for recharge and (b) the areas that retain the highest proportion of vertical paths with depth, independent of the site size. Figure 6 presents this quantitative comparison of the seven sites. For each site, we assessed both the whole site (Figure 6a,b) and the optimal area within Vadose Zone Journal F I G U R E 6 Quantitative comparison of vertical paths vs. depth between all seven sites. We compare the number and percentage of vertical paths in two areas: Row 1, the whole site; Row 2, the optimal area; where percentage is the number of vertical paths normalized by the number at zero depth. The selection of an optimal area typically reduces the decay rate of vertical paths especially in the shallow (<8 m) portion of the subsurface the site (Figure 6c,d), which was chosen using the approach described in the preceding section. Here, to demonstrate the utility of this approach to site comparison, we will present a few important observations that can be made.
As seen in Figure 6a and c, some fields have many more vertical paths because of their larger starting surface areas. In each plot, the number of vertical paths at 0-m depth indicates the ground surface area of the site or selected optimal areas. The Nichols site has the largest surface area (1.5-24× the size of the other sites) and consistently has more vertical paths with depth than all other sites despite having significant decay beginning at ∼10 m. The same is true when comparing the selected optimal areas for the sites. Nichols has a large subsurface volume through which water can quickly move past the root zone, and it has a larger optimal area than any other site. Were there a desire to select a single site based on the criteria we are using and under the assumption that the recharge operations were not limited by water availability, Nichols would be ranked as the best of the sites.
Another key observation, seen in Figure 6b and d, is that some sites retain a much higher percentage of vertical paths as depth increases regardless of the site size. A scenario in which this might be important would be one where the objective was to spread water over an entire site, or its optimal recharge area, and to maximize the value of the site (or the optimal area) per unit area covered with water, where value is defined as the retention of vertical paths. If we considered only the number of vertical paths, as noted previously, Nichols would seem to be the best choice. Looking at the percentage vertical paths, however, Swall-North would emerge as comparable; for example, despite its significantly smaller size, the whole site retains ∼5% vertical paths to the depth of its water table at 24 m and ∼10% within its optimal recharge area. Nichols, on the other hand, loses 100% of its vertical paths by the depth of its water table at 41 m, even in the optimal area.
Assuming that all of these sites are farmland where ponding is of high concern, the number and percentage of vertical paths available from the surface to the water table in both the entire site and the optimal area would lead us to rank Nichols, Swall-North, and Almond-West as the best sites for recharge. Although Martin-West was comparable in terms of the results shown in Figure 6d, this site showed little overlap between contours in its metric maps (Supplemental Figure S7) and thus should not be selected for recharge using our criteria. We note that the top three sites are also the largest. When searching for quality recharge areas, investigating larger fields increases the likelihood of finding an optimal area of sufficient size for recharge.

The value of generating multiple sediment-type models
The results from the previous section were derived from averaging across 200 generated sediment-type models for each site; therefore, the metric maps ( Figure 5) and vertical-path decay curves ( Figure 6) represent expected values. However, we could also use the metric maps and vertical-path decay curves from each individual model to capture the uncertainty in these curves as a result of differences in the models. In Figure 7 we show (as gray lines) the decay in the percentage of vertical paths with depth for each model within each site's optimal area. The variation in the vertical-path decay curves displays the differences that could be present at a site and can F I G U R E 7 Vertical path decay curves in the optimal areas. The individual model results highlight the variability around the various average representations of the sites be used to tailor the assessment of the site so as to avoid specific undesired outcomes. For example, if avoiding shallow ponding was a high priority, a choice could be made to reject sites where any model shows a significant decay in the number of vertical flow paths near the surface. Alternatively, if shallow ponding were of little concern, such as for reclamation of a floodplain, sites with many models that retain a high percentage of vertical paths to the water table (e.g., Swall-North and Swall-South) could be selected for further analysis regardless of the level of uncertainty in percentage of vertical paths in the shallow subsurface.
In Figure 7 we show for each optimal area the mean of the individual vertical path decay curves (blue line) and the result from Figure 6d (red line). If the 200 binary-sedimenttype models for the optimal area at a site are similar, the red and blue curves will be similar. Alternatively, the red and blue curves will be significantly different when there is high variability in the individual models. For example, the individual vertical path decay curves from Almond-West are tightly clustered, and its blue and red curves show good agreement, particularly beyond a depth of 20 m. Swall-South, on the other hand, has a large spread in its individual curves. Consequently, its blue and red curves differ in magnitude and shape through-out the subsurface, suggesting we have little certainty about the subsurface and its effect on recharge at this site. The difference between the red and blue curves represents the extent to which the spatial structure varies from the average. In contributing models of the subsurface to be used by others in planning and implementing recharge operations, it is important to highlight the level of uncertainty in terms of the spatial distribution of sediments in the subsurface or how water will flow during recharge. Plots such as those in Figure 7 can be used to assess and demonstrate this uncertainty and guide informed decisions for locating recharge sites.

CONCLUSIONS
In this paper, we have demonstrated the use of tTEM data to select the optimal area within a site for surface spreading recharge based on the vertical distance between the surface and the shallowest low-K unit underlying each point at the surface (Metric 1) and the length of preferential flow paths between each point at the surface and the water table (Metric 2). Even with our simplified representation of the subsurface using a binary system, our results show high spatial variability in the distribution of binary-sediment-type at each site and, consequently, high variability in our metric assessments, underscoring the need for high-resolution geophysical imaging. The comparison of sites using these metrics becomes challenging because of differences in the size of the sites and the depth to the water table. We therefore developed an additional way to assess the relative suitability of each site using the decay in the number of vertical flow paths as a function of depth. A limitation of our method is that we reduced the sediment types in the subsurface to a binary system (high-K and low-K). Although this is a simplified representation of the subsurface, it allowed us to demonstrate the use of the geostatistical methods and the eikonal solver. Our methods could be extended to create fully populated models of the fraction of coarse-grain-dominated material. With available field data that makes it possible to determine a relationship between sediment type and hydraulic properties (e.g., infiltration rate and K), we could then reconfigure the eikonal solver to determine actual travel times within the vadose zone.
We have developed methods that can be used with electromagnetic imaging data to rapidly identify promising sites for surface spreading recharge with respect to their ability to efficiently transmit water from the surface to the water table. Further studies could then be conducted at these sites including, for example, infiltration tests, modeling and sampling to better quantify flow rates and flow paths, and an assessment of potential effects on water quality. We wish to emphasize that the factors deemed most important in the selection of a site or an area within a site will vary depending on the goals and risk tolerance of the recharge operation. With the concerns of groundwater depletion in the CV, new approaches that can reduce the time and expense involved are needed to support the growing commitment to surface spreading for recharge. Electromagnetic imaging techniques with the complementary methods for data analysis have an important role to play.

S U P P L E M E N T A L M A T E R I A L
The supplemental material includes further details regarding the geostatistical methods used in our workflow. In particular, we describe the development of the empirical variogram model, the distributions for each random input variable, and how we sampled the parameter distributions during each geostatistical run. The supplemental table provides the variogram model and grid parameters for each site. We also include figures displaying the workflow for obtaining the binary-sediment-type models and the metric map results for the remaining six sites (Almond-East, Nichols, Swall-North, Swall-South, Martin-West, and Martin-East).

A C K N O W L E D G M E N T S
This research was supported by funding to R. Knight from the Stanford Woods Institute for the Environment (REIP 2019 grant award) with additional support from the Gordon and Betty Moore Foundation (GBMF6189). We wish to thank the landowners who provided access to the sites and Aaron Fukuda, General Manager of Tulare Irrigation District, who played a critical role in facilitating the acquisition of the field data.

C O N F L I C T O F I N T E R E S T
The authors declare no conflicts of interest.