Identification of fresh submarine groundwater off the coast of San Diego, USA, using electromagnetic methods

Climate change has a pronounced effect on water resources in many semiarid climates, causing populated areas such as San Diego County (USA), to become more vulnerable to water shortages in the coming decades. To prepare for decreased water supply, San Diego County is adopting policies to decrease water use and to develop additional local sources of water. One new local source of freshwater is produced by a desalination facility that purifies brackish groundwater from the coastal San Diego Formation. This formation has been studied extensively onshore, but little is known about the geology or groundwater quality offshore in the adjacent continental shelf. Because most groundwater systems are interconnected and complex, further analysis is needed to identify offshore geology, possible sequestration of freshwater in the shelf, and potential pathways for saltwater intrusion. This comprehensive understanding is important because seawater intrusion may limit use of the San Diego Formation and longevity of desalination facilities. Controlled-source electromagnetic methods are uniquely suited to detecting offshore groundwater as they are sensitive to changes in pore fluids such as the transition from fresh to brackish groundwater. This paper describes results from surface-towed electromagnetic surveys that mapped the pore-fluid salinity and possible fluid pathways in the continental shelf off the coast of San Diego. The results indicate a considerable volume of fresh-to-brackish groundwater sequestered in the shelf, both in continuous lenses and isolated pockets, that appear influenced by fault systems and shallow stratigraphy.


Introduction
Climate change significantly impacts the water resources available in semiarid climates and populations living within these regions are likely to become more vulnerable to water shortages. The southwestern United States is a semiarid region hosting several large population centers, including coastal southern California with the three highest population density counties in the state (Los Angeles, Orange, and San Diego). This region is dependent upon water from snowpack in north-ern California which is expected to experience a 48-65% loss by the end of this century and the Colorado River which is expected to see an additional 20% drop in flow by 2050 (Udall and Overpeck 2017). San Diego County is especially vulnerable to these predicted water disruptions as the county currently purchases 85-90% of its water from northern California and the Colorado River.
So that water supply agencies can adapt to the changing climate in semiarid regions, groundwater resources could be further developed to account for an expected decrease in overall water supplies. In coastal semiarid areas such as San Diego County, development of groundwater resources is complicated by saltwater intrusion and the potential presence of fresh groundwater reservoirs offshore. Groundwater in the San Diego region is further complicated by its location along the North American-Pacific Plate boundary where several strikeslip fault systems influence coastal geomorphology and subsurface geology.
The San Diego Formation (SDF) is a heavily faulted coastal aquifer that is currently supplying about 5% of the total water used in the San Diego area. The salinity of water extracted from the SDF generally ranges between slightly-saline to brackish (1,500-2,500 ppm), with some regions producing freshwater (defined as less than 500 ppm). The combination of increasing water demand and diminishing sources of imported water suggests that groundwater extraction from the SDF may increase, in particular to supply municipal water to communities in southwestern San Diego area (San Diego County Water Authority 2021). Groundwater extracted from the SDF could originate from naturally recharged precipitation, or alternatively, from aquifer storage and recovery, a method that involves artificial recharge of the SDF with excess water during periods of low water use followed by extraction of this water during periods of high-water use (Keller and Ward 2001).
To provide water agencies with the scientific understanding to evaluate these water-management scenarios, the United States Geological Survey (USGS)-in cooperation with the Sweetwater Authority, the City of San Diego, and the Otay Water District-began an extensive hydrogeology project in 2001 to map the geology and the water resources of the San Diego area (US Geological Survey 2020). A major part of this project focused on mapping the San Diego coastal aquifer, determining groundwater recharge and discharge rates, identifying groundwater flowpaths, and assessing groundwater quality. Installation of 16 USGS multiple-depth, monitoring-well sites enabled collection of three-dimensional (3D) geologic data, groundwater levels, and groundwater quality. The project study area, the generalized geology of the San Diego area, and the locations of USGS monitoring-well sites are shown in Fig. 1. Hydrogeologic understanding gained from these sites and from development and testing of geologic and hydrologic models helped to define the onshore coastal geology, groundwater recharge and discharge rates, and groundwater flowpaths.
Data collected from these onshore USGS monitoring sites from 2001 to 2020 indicate that fresh and brackish groundwater flows onshore from the continental shelf under San Diego Bay, driven by onshore groundwater extraction from the SDF (Anders et al. 2013). Additionally, noble gas isotope analyses of sampled well water indicate that groundwater age generally increases with proximity to the coastline (Seltzer et al. 2019). These observations, combined with the geologic model of the San Diego coastal aquifer, suggest that the SDF extends about 5 km west of the present shoreline and indicate that freshwater may have been sequestered in the SDF when the continental shelf was subaerially exposed approximately 20,000 years ago (Danskin 2012). Because most hydrogeologic information has been collected onshore, the location and extent of offshore groundwater are not well understood. Historically, electromagnetic (EM) surveying has been used on land to map groundwater resources with great success (Knight et al. 2018;McNeill 1988). Also, EM methods have been proposed by many authors to identify salinity of submarine groundwater (Cohen et al. 2010;Haroon et al. 2021;Micallef et al. 2020Micallef et al. , 2021Post et al. 2013) and several controlled-source electromagnetic (CSEM) surveys have been successful at identifying groundwater offshore, most notably offshore Martha's Vineyard (Gustafson et al. 2019) and Hawaii (Attias et al. 2021). The geologic setting of San Diego varies from these previous studies in that the presence of offshore groundwater has not yet been confirmed, and the SDF is heavily faulted, possibly isolating pockets of freshwater offshore. This current study aims to further the understanding of the onshore groundwater system by characterizing and mapping pore fluids within the offshore part of the SDF using CSEM methods.

The CSEM method
The marine CSEM method uses a human-made source of EM energy that passes through seawater and propagates into the seafloor and to receivers, deployed either on the seafloor or towed through the water, which measure the resulting electric fields. The recorded fields are proportional to electrical resistivity, making the method suitable for detecting fresh groundwater, which is significantly more resistive in comparison to seawater or the surrounding seafloor sediment. Additionally, as the porosity and cementation of the SDF are well studied onshore, the salinity of the submarine aquifer can be predicted with electrical conductivity data by applying Archie's Law. Using observations and data collected from USGS monitoring-well sites and hydrogeologic models, the porosity of the SDF is constrained to values between 30-45% with a mean value of 38% and the cementation exponent used in Archie's Law was determined to be 1.8. Of these two variables, porosity has the greatest effect on the predicted salinity of the pore fluids. The relationship between salinity and bulk resistivity given these values is shown in Fig. 2.
Because the water depths on the continental shelf are shallow and the SDF is not expected to extend more than 500 m below sea level, the small and low-power surfacetowed CSEM system of Sherman et al. (2017) can be used effectively to map the SDF. This system has the added benefit that it is inexpensive to use, can be hand-deployed, and requires a nonspecialized vessel to efficiently survey the shelf. The system, illustrated in Fig. 3, consists of an EM transmitter that generates an electric dipole, four electric field receivers, and a dorsal device that collects waterconductivity and water-depth data. The electric dipole, four receivers (referred to as porpoises), and dorsal are towed on the surface of the water behind a vessel traveling between 2 and 4 knots (approximately 1-2 m/s) on floating high-molecular-weight polyethylene rope. The frames of the receivers and dorsal are made from rigid plastic that is designed to slough off kelp or other ocean debris, hold a rigid dipole 0.67 m beneath the sea surface, and provide mounting for a vertical global positioning system (GPS) mast that doubles as a flashing beacon for better visibility to other vessels. To maximize sensitivity to the predicted geometry of the SDF, the receivers were separated by 100 m, creating a maximum source-receiver spacing of 400 m and a total towed array length of 430 m.
The EM transmitter operates on 110-240 VAC power and can output a GPS stabilized binary or ternary waveform of as much as 50 amps on an antenna typically 50 m long. For this survey, the transmitter outputs a 30-amp current-controlled 2.5-Hz fundamental waveform-D, described by Meyer et al. (2011), on a 10-m antenna. This configuration resulted in a 300-amp-m horizontal electric dipole moment, which is a signal strength well suited to the source-receiver spacing chosen for this survey. Waveform-D was chosen as it generates higher amplitude responses in the 1st, 3rd, 7th, and 13th harmonics than other binary and ternary waveforms, allowing the collection of data across a broad frequency range (Meyer et al. 2011).

Survey data quality
The data presented here were collected from September 2019 to September 2020 on 5 separate days of surveying aboard the Scripps Institution of Oceanography research vessel Bob and Betty Beyster. The cruises were designed to initially survey the shelf with a broad exploration method followed by increasingly focused surveys targeting resistive anomalies identified within the SDF from the initial results. The location of the survey lines presented here are shown in Fig. 4.
A total of 141 km of high-quality inline electric field response data were collected. For each day of survey, the data were inspected for quality before further processing. GPS data from the towed receivers were lost on the first day of data collection because of user error; however, redundant onboard GPS data were found to be sufficient and thus used. For the Fig. 2 Graph illustrating the predicted changes in bulk resistivity with changing pore fluid salinity using characteristic SDF attributes such as porosity values between 30 and 45%. Pore fluid salinity can be estimated from the resistivity models generated with controlled-source electromagnetic (CSEM) data using Archie's Law which describes the relationship between cementation factor, porosity, bulk resistivity, and pore fluid conductivity of a material Fig. 3 Schematic of the CSEM array used for the survey offshore San Diego County, California. The dipole is a 10-m horizontal electric dipole source. For the surveys discussed here, the dipole center is located 28.2 m from the GPS mast-mounted on the vessel following three survey days, all instruments functioned as designed.

Results
Amplitude and phases of the CSEM response functions were extracted from the collected raw time-series data using a method detailed by Meyer et al. (2011). To increase the signal-to-noise ratio, the resulting transfer function estimates were stacked using an arithmetic mean to obtain the transfer function estimates for every 30 s along with an error estimate. This method yielded high quality amplitude and phase response data for all four receivers as a function of position and frequency. Amplitude data for all four receivers at 2.5, 7.5, 13.5, and 32.5 Hz and phase data for 7.5 and 13.5 Hz were included in the inversion code. The amplitude data for these frequencies were well above the noise floor and sensitive to the known depths of the SDF. Phase was neglected at 2.5 Hz as there was little variability in the signal and addition of these data in the models resulted in an unbalanced misfit between frequencies due to the model fitting noise in the 2.5 Hz phase. Phase was also neglected at 32.5 Hz as there was significant impact from clock drift at this frequency.
The modeling software used in this study is the publicly available, goal-oriented, adaptive, finite-element two-dimensional (2D) MARE2DEM inversion and modeling code of Key (2016). This code uses Occam's Inversion, a method that regularizes the inversion to the smoothest resistivity model that fits the data to a specified misfit (Constable et al. 1987). CSEM data were scrutinized manually for obvious outliers and subjected to a 2% error floor before being included in the inversion as finite-length dipoles.
The starting models included the seawater as a fixed parameter, using conductivity data collected by the dorsal and available bathymetric data. Therefore, the free inversion regions were reduced to the area below the seafloor and set to a uniform starting resistivity of 1 Ωm. Inversion parameter grids were constructed using 30-m-wide quadrilateral cells that increased in height with depth to mimic the loss of resolution of the EM method. Due to the adaptive nature of the MARE2DEM code, the computation mesh was allowed to fine where necessary to produce accurate responses. The resolution depth was investigated by including a highly conductive (0.1 Ωm) or highly resistive (1,000 Ωm) layer as a base to the model at varying depths. From these tests, the limits of sensitivity are estimated to be at a depth between 400 to 500 m below sea level; thus, the final model interpretation and presentation is limited to these depths. Isotropic inversions fit the amplitude data, but these models produced biased phase residuals and a normalized root-mean-square of 1 was not achieved, suggesting the need to include anisotropy in the models (Sherman and Constable 2018). Thus, differing resistivity values in the vertical and horizontal directions were included in the models, and allowed for simultaneous fit of both the amplitude and phase data. These inversions resulted in anisotropic ratio (vertical resistivity/ horizontal resistivity) of up to 5 in the SDF, consistent with the geology of the SDF which has been observed to have a series of fining-upward stratigraphic sequences. The anisotropic resistivity inversions fit the data to a normalized rootmean-square of 1 with a two-percent error floor. See Fig. 5a for an example resistivity model and Fig. 5b for the data and data fits from this model, which are representative of the data quality and model fit of the surveys and resulting models discussed in this paper. The final resistivity models that focus on the main resistive anomalies offshore Imperial Beach, located in the southern coastal San Diego area, are shown in Fig. 6.

Discussion
The structural San Diego basin is a pull-apart sedimentary basin which underlies much of the San Diego area and was formed concurrently with deposition of the SDF during the middle-to-late Pliocene to Pleistocene (Kennedy and Peterson 1973. The San Diego basin is observed to increase in depth from north to south, with the SDF also thickening until both the basin and the SDF gradually shallow and flatten near the USA-Mexico border. The SDF also shallows and thins both to the east and west and is interpreted to be thickest beneath the city of Imperial Beach (Huntley et al. 1996). The western edge of the SDF is not well defined, but current hydrogeologic models which incorporate gravity, seismic, and borehole data indicate that the western boundary occurs approximately 5 km west of the present shoreline.
The 2D models of vertical resistivity contain numerous highly resistive anomalies offshore Imperial Beach. The depth to these resistive anomalies is approximately 50-80 m below sea level, which is consistent with the depth at which the SDF is encountered in nearby USGS monitoring-well sites. Additionally, Darigo (1984), using acoustic reflection methods, interpreted the locations of these resistive features to be the SDF (see Fig. 5a). As such, the locations of these resistive anomalies are interpreted to be within the SDF. Archie's law, the resistivity models, and known geologic parameters of the SDF were used to predict the salinity of the pore fluids within these resistive anomalies. The resistive anomalies (>10 Ωm) offshore Imperial Beach correspond to areas where pore fluids are fresh-to-slightly saline (<3,000 ppm). Assuming a porosity of 30% and a thickness of 100 m, consistent with model geometries, the volume of the mapped fresh-to-slightly-saline groundwater is approximately 390 million m 3 (103 billion gallons). The onshore volume of the SDF is estimated to be between 341 and 532 million m 3 (City of San Diego 2016). The estimated offshore volume of groundwater nearly doubles the total volume of the SDF which is consistent with predictions by Keller and Ward (2001).
The resistive anomalies appear to be north-south trending and are separated by down-dropped blocks controlled by segments of the Silver Strand fault zone. In previous seismic surveys, the down-dropped blocks have been observed to be collocated with Pleistocene and Quaternary paleochanlnels that were incised into the continental shelf during the last glacial maximum, a time when the continental shelf was subaerially exposed (Darigo 1984;Graves et al. 2021). Similar west (W) and east (E) paleochannels (PC) were found in this CSEM survey and are indicated by conductive U-shaped features labeled respectively as WPC and EPC on Figs. 5a and 6.
The paleochannels marked EPC are collocated with areas in which the SDF appears to be conductive, signifying saline pore fluids. This observation indicates, in these areas, that the paleochannels may have been incised into the lowerpermeability sediment that functionally trapped the fresh SDF pore fluids, allowing for fluid migration into and out of the SDF. As sea level rose, drowning the shelf, these incisions would allow for diffusion of saltwater into the SDF below the paleochannels. Saline pore fluids appear to have infilled the SDF below the paleochannels marked EPC. Conversely, the paleochannels marked WPC, which run parallel and west of Silver Stand section 1, also appear to have been incised into the SDF, but the resistive anomalies remain intact below these paleochannels. This preservation of fresh-to-slightly-saline pore fluids within the SDF, despite the presence of overlying Fig. 5 a Plot of 2D vertical resistivity model from one line of data, with data and data fits from the third receiver in the array (300-m transmitter to receiver offset). 2D vertical resistivity model with overlain geologic interpretations (marked in white) by Darigo (1984) are shown (a). The location of the vertical resistivity model is shown in Figs. 4 and 6 and the location of the geologic interpretations by Darigo (1984) are shown on Fig. 4. The resistivity model is run with an anisotropic penalty weight of 0.05 and is fit to a normalized root mean square of 1. Warm colors (red, yellow) indicate areas where pore fluids are fresh-to-slightly saline (<3,000 ppm total dissolved solids), assuming a porosity value between 30 and 45%. Cool colors (blue) indicate moderately-saline-to-highly-saline pore fluids. The white labels are from Darigo (1984): H Holocene sediment; P Pleistocene channel fill; Tmp undifferentiated Miocene and Pliocene (SDF). Conductive U-shaped features interpreted as paleochannels are labeled EPC (eastern paleochannel) and WPC (western paleochannel). b Plots the data and error (error bars) and modeled response (black lines) for all data used in the inversion from the third receiver in the array. The data quality and data fit presented here are representative of both the data quality of the other receivers in the array and the data fit achieved by the other inversion results discussed in this paper paleochannels, may result from locational differences in the SDF stratigraphy.
Depositionally, the SDF is a mixture of shallow marine and nonmarine sediment (Keller and Ward 2001). As a result, the SDF contains multiple fining-upward sequences, sandy marls, conglomerates, and clayey layers (Ellis 1919) which is evident in the anisotropic ratios generated from the final resistivity models. The fining-upward sequences and clayey layers may form a series of less permeable sediment, which could aid in preserving freshwater in the more permeable layers below. This vertical restriction to groundwater flow is illustrated in the confined response observed during regional-scale, multiple-well aquifer tests and single-well pumping tests (Brandt et al. 2020). The depositional character of the SDF may be one reason that the fresh-to-slightly-saline pore fluids are preserved below the paleochannels marked WPC, west of Silver Strand section 1. In these areas, the fining-upward sequences may form a series of caps to pore fluids within the SDF. The WPC paleochannels may have been incised into some, but not all, of these caps, which could leave parts of the SDF undisturbed.
An additional reason for irregularities in the highly resistive locations could be the discontinuous nature of the more permeable parts of the SDF, which are separated by clayey sediment (Huntley et al. 1996). This depositional pattern may explain the discontinuous nature of the resistors offshore, but it does not explain the boundaries being collocated with fault traces. Therefore, geometries and locations of the highly resistive features mapped via CSEM may result from a combination of stratigraphy, regional faulting, and lithification.
Finally, the Silver Strand fault zone is made up of a series of strike-slip, right-lateral faults within a generally transtensional regime (Maloney et al. 2016). These faults transect the offshore SDF, creating migration pathways for freshwater flushing. This phenomenon, in the form of freshwater springs along the fault zone, has been noted in historical records of Coronado Island (Keller and Ward 2001). The 2D resistivity models suggest upward fluid migration from the lower fresh-to-slightly-saline anomalies along the Silver Stand fault sections 1 and 3. These anomalies extend in some cases to the seafloor. The vertical orientation of, depth to the tops of, and collocation with known faults of these resistive features suggest the presence of freshwater flushing along Silver Stand fault sections 1 and 3. colors (red, yellow) indicate areas where pore fluids are fresh-to-slightly saline (<3,000 ppm total dissolved solids), assuming a porosity value between 30 and 45%. Cool colors (blue) indicate moderately-saline-tohighly-saline pore fluids. Conductive U-shaped features interpreted as paleochannels are labeled EPC (eastern paleochannel) and WPC (western paleochannel)

Conclusion
The use of CSEM offshore San Diego, California, USA, has resulted in a resistivity model identifying sections of the SDF with substantial volumes of fresh-to-slightlysaline pore fluids, sequestered offshore. These results also indicate that CSEM could be a useful and efficient tool for mapping the extent of other coastal aquifers across diverse geologic settings. The results, however, also illustrate the potential complications related to withdrawing water from the naturally charged SDF or using the SDF for aquifer storage and recovery. Possible freshwater flushing from these sections has been interpreted along the transtensional faults related to the Silver Strand fault zone. In this case, offshore faulting may act as conduits for saltwater intrusion into the onshore SDF. Conversely, offshore faulting may act as partial barriers to groundwater flow and inhibit movement of fresh or saline groundwater into the onshore SDF. Further research and data collection are needed to understand how well this offshore, submarine groundwater is connected hydraulically to the onshore SDF.