A Multi‐Hydrogeophysical Study of a Watershed at Kaiwi Coast (Oʻahu, Hawaiʻi), using Seismic Ambient Noise Surface Wave Tomography and Self‐Potential Data

We study a watershed containing an old stream valley at the Kaiwi Coast on Oʻahu, Hawaiʻi. The study site has undergone significant hydrogeological changes throughout its history, and at present displays complex interactions of a variety of geological formations. We highlight an innovative and efficient combination of two hydrogeophysical methodologies for studying groundwater systems: ambient noise surface wave tomography (ANSWT), and self‐potential (SP) measurements. We collected (1) 5 days of continuous seismic data, using a total of 54 autonomous seismic nodes, exploiting both ocean wave and traffic noise sources and (2) 386 SP measurements with an average 20 m spacing, in the same area. We showcase how these two data sets complement each other in a joint‐interpretation framework. The SP data display the distribution of groundwater in the region, as well as zones of higher infiltration and/or flow rates. Using ANSWT, we identify the structural contact between the older Koʻolau basalt formation and the younger post‐erosional Honolulu volcanics, as well as the distribution of sedimentary valley fill and coastal deposits. The joint interpretation of the SP and seismic data clearly illustrates that groundwater flow occurs in the identified paleo‐channels at the erosional surface of the basaltic bedrock. Furthermore, it highlights that the basaltic bedrock of the ridges likely forms an important groundwater flow unit. Some flow may occur in the shallower valley fill material. Complementing SP data with seismic data enables us to place groundwater flow inferences, as made from SP data, at depth.

Furthermore, at the coast, sedimentary deposits can form a low-permeability confining cap-rock (Figure 1), which overlies higher permeability volcanic formations and impedes the seaward discharge of freshwater (Oki, 2005). In addition to these deposits, some valleys were filled in by marine and terrestrial sediments during a period of higher sea level than the present-day scenario (Oki, 2005). Other geological features at a variety of spatial scales can have a major impact on the flow properties and distribution of groundwater as well, such as the presence of low-permeability dike intrusions and hydraulically conductive lava tubes.
These heterogeneities and strong gradients in hydraulic conductivity influence the amount and characteristics of groundwater flow in the bedrock. This results, for example, in a delicate trade-off between groundwater that flows predominantly through basaltic bedrock and the basaltic ridges of the ridge-valley complex, including its recharge into deeper aquifer systems, versus flow through shallow alluvium (Figure 1). For example, on the leeward side of Oʻahu, the older alluvium is of hydrological importance due to its aforementioned relatively low hydraulic conductivity. As a consequence, ridges composed of basaltic bedrock with much higher hydraulic conductivity, form the dominant flow units on this side of the island (Figure 1). However, that scenario differs on the windward side of Oʻahu, where Koʻolau rift zone dike swarm intrusions are abundant (Takasaki & Mink, 1985;Walker, 1986Walker, , 1987, making the basaltic bedrock less hydraulically conductive, allowing stream valleys to become important groundwater flow units (Mathioudakis, 2018).
Despite having access to these various methodologies, from a geophysical perspective, Hawaiian hydrogeological systems are extremely challenging to image and characterize. Different methodologies suffer from different practical challenges. For example, the quality of active seismic imaging in basaltic formations (providing primarily information on the geological structure and rock properties) suffers from a variety of challenges, including high seismic velocities, strong velocity contrasts with, for example, overlaying sediments, anelastic attenuation, and strong internal scattering due to the presence of fractures and faults (Behera & Sen, 2014;Planke et al., 1999). Electrical methods, such as electrical resistivity tomography (providing primarily information on the pore-fluid distribution and its properties), can occasionally suffer from both source (difficulties injecting the current) and receiver-side coupling issues due to high contact resistance, negatively impacting the measurements (Wilkinson et al., 2012). Magnetotelluric measurements encounter several practical challenges as well, such as the presence of cultural noise on densely inhabited islands, large conductor effects from the surrounding ocean, high contact resistance, steep topography, as well as difficulties burying the magnetometers and electrode dipoles in basaltic bedrock (Hoversten et al., 2003). As a consequence, understanding the complex subsurface hydrology of the islands and its role in watershed functioning remains a challenging task.
Given that each geophysical method has its own benefits and drawbacks, and informs on different aspects of the aquifer system, multi-geophysical (multi-physics) approaches that combine different geophysical methods that complement each other, have become increasingly popular in recent years (Colombo & De Stefano, 2007;Jardani et al., 2013;Robert et al., 2011). Multi-geophysical data can be integrated through, for example, a joint-inversion framework, and/or a joint-interpretation framework. The common objective is to try and overcome the limitations of single-method studies and to gain a more complete and detailed understanding of the subsurface through the use of complementary data sets.
In this study, we combine, for the first time to our knowledge, the SP and seismic ambient noise surface wave tomography (ANSWT) methodologies in an integrated multi-geophysical interpretation framework, aimed at understanding the present-day hydrology of the old Makapuʻu stream valley at Kaiwi Coast, Oʻahu, Hawaiʻi. The study site itself, a state park, is not a primary target for groundwater exploration. However, it is a non-urban, accessible area on the leeward side of the island of Oʻahu, that serves as an ideal analog study site for other valley-ridge complexes on leeward Oʻahu. As such, with this study, we aim to increase our understanding of the main groundwater flow units and hydrogeological characteristics of the valley-ridge complexes in watersheds on the leeward side of Oʻahu.
Both the SP and ANSWT methods (although innovative for water resource applications) have, individually, been proven to be robust, efficient, and successful methodologies for studying volcanic island environments (Brenguier et al., 2007;Finizola et al., 2004;Mordret et al., 2015;Revil et al., 2004;Soueid Ahmed et al., 2018;Zlotnicki & Nishida, 2003). The SP methodology provides primarily information on preferential groundwater flow (which may coincide with major structural features), whereas ANSWT offers structural information on the main geological units by providing 3D images of primarily seismic shear wave velocity. Using a joint-interpretation framework, we analyze the SP and ANSWT data sets that we have collected in the old Makapuʻu stream valley and coastal zone. This study site provides an accessible case study that can be used as an analog for other watersheds on Oʻahu and on other islands.

Study Site
Our study site is an old watershed located in the southeasternmost part of the Koʻolau volcano. The Koʻolau volcano is the youngest of the two volcanoes that comprise the island of Oʻahu: its subaerial shield flows are estimated to be ∼2.9 million-2.1 million years old. Subsequently, these shield flows were eroded for a considerable time, after which rejuvenated activity occurred. This rejuvenation volcanism has lasted until perhaps as recently as 30,000-50,000 years ago (Haskins & Garcia, 2004;Rowland & Garcia, 2004). Figure 2a shows a geological map of the study area. Figure 2b displays a satellite image with a map of the geophysical acquisition. Our study site involves both the old Makapuʻu stream valley and the coastal plain directly neighboring the Makapuʻu ridge (the easternmost ridge of the island of Oʻahu, indicated as Ridge 3 in Figure 2b), at the Kaiwi coast. It covers multiple valley-ridge complexes, indicated in Figure 2b. In the north, the stream valley is abruptly cut-off, showing steep cliffs. Moore (1964) provided evidence that much of the northeast flank of the Koʻolau volcano had collapsed as a debris avalanche onto the ocean floor. Further evidence of the occurrence of giant avalanches at Hawaiian volcanoes was provided by side-scan sonar data collected around the eight main Hawaiian islands (Moore et al., 1989). The debris avalanche associated with the Koʻolau volcano is considered one of the largest landslides on Earth, with blocks up to 40 km long and 2 km tall having slid as far as 100 km (Moore & Clague, 2002). Any avalanche scar is anticipated to be seaward of the most seaward remnant of the Koʻolau volcano, placing it at least GROBBE ET AL.  a kilometer oceanward of the current coastline (Rowland & Garcia, 2004). Immediately after the massive landslide, the cliff could have been around 2,000 m high. However, due to its likely high instability, smaller avalanches and rapid marine and stream erosion would have caused this cliff to be eroded back, landward (southwest). This process is referred to as "parallel retreat" and has been documented at a number of locations on Earth (Widdowson & Cos, 1996). Marine erosion has likely further contributed to the development of the cliff into its current location and shape. The cliff cut-off on the northern side of the stream valley has clearly impacted watershed functioning. The former recharge zones at higher elevation have been removed due to the debris avalanche and subsequent erosion, effectively removing the source region of groundwater flow. In Figure 2b, at the north of the Makapuʻu valley, we can notice the low elevation at the top of the study site. Additionally, the stream valley is located on the leeward side of the island, which has drier conditions than the windward side, and the southeastern tip of Oʻahu is on the dry side of the leeward precipitation range.
In terms of the geology of the study site (Figure 2a), most of the Koʻolau lava deposits are aphyric, but a few are very rich in olivines (≥20 volume %). The distributions of the olivines have been used to calculate lava viscosities, and lavas such as these at Makapuʻu yield values of 10 3 -10 4 Pa s. The ridges on either side of the valley are primarily composed of shield-building flows (Rowland & Garcia, 2004). The valley itself consists of sedimentary valley-fill deposits that are underlain by the same basaltic shield-building flows (Koʻolau basalt). This interface may have undergone mechanical weathering or chemical alteration into a saprolite formation.
In addition to these deposits, an alkali olivine basalt flow (Kalama flow) was erupted from a scoria cone in Kalama Valley westward of the study area, around 35,000-80,000 years ago (Gramlich et al., 1971). It flowed southward and subsequently spread out of the valley; northeastward along the coastline toward our study site, while its extension was stopped by Koko Crater to the west.
These flow deposits show characteristics of both Pāhoehoe (flows having smooth, hummocky, glassy surfaces surrounding highly vesicular interiors) and 'A 'ā (flows with a dense, massive, discontinuous central phase, the 'A 'ā core, bounded by spiny, fragmented lava breccia called clinker) (Lau & Mink, 2006). A single Pāhoehoe unit is highly porous, but has little inherent permeability since the vesicles may not be connected. However, as a whole it may be highly permeable because of secondary structural features such as unconformably matched flow units, lava tubes, and cooling joints (Lau & Mink, 2006). There tend to be relatively few vesicles in 'A 'ā units, and if present they are typically large and irregular in shape. Practically all permeability in 'A 'ā is the result of structural features, such as openings in clinker beds and vertical cooling joints (Lau & Mink, 2006). In terms of hydrological characteristics, Pāhoehoe and 'A 'ā cannot be treated separately; a single unit is of limited spatial extent, and Pāhoehoe and 'A 'ā units can alternate each other forming together a representative hydrogeological volume (Lau & Mink, 2006).
Around some of the tidepools near the shoreline, crude columnar cooling joints are evident. These Honolulu volcanic deposits reached the eastern limit of our study area, almost up to the Makapuʻu ridge (Figure 2a). In recent years, the flat coastal plain at the bottom and west of the stream valley was almost completely inundated by a tsunami in 1946. Wave runup heights were more than 6 m in this area (Macdonald et al., 1983), and the water reached inland essentially to where the highway is located today. This event has left additional sedimentary deposits on the plain.
From a hydrogeological standpoint, we can summarize our study area to first order as follows: an old stream valley containing sedimentary valley fill deposits combined with underlying, potentially weathered lava flows (saprolite). Both formations typically have a much lower hydraulic conductivity than unweathered basaltic bedrock and hence can form potential barriers to groundwater flow. Basaltic bedrock ridges are located on both sides of the stream valley, typically having higher hydraulic conductivity, with the west side of the stream valley consisting of another valley-ridge complex. Furthermore, at the mouth of the stream valley, and westward, a coastal plain is located with Koʻolau basalt overlain with a more recent lava flow (Honolulu volcanics), followed by tsunami and sedimentary deposits. In this study, we aim to identify any major groundwater flow paths in this old stream valley and the coastal plain area, using a joint acquisition and joint interpretation of SP and ambient seismic noise geophysical data.

Methods
We have acquired both ambient seismic noise data and SP data covering the whole study area (Figure 2b: red and blue dots for seismic, and yellow dots for SP). A total of 54 seismic nodes were deployed (29 1-component and 25 3-component), combined with a total of 386 SP measurements. These data were acquired, processed, and inverted as part of the first edition (summer 2019) of a newly established graduate-level summer fieldschool program at the University of Hawaiʻi at Mānoa called "Hydrogeophysics in Volcanic Environments."

Self-Potential
The SP method measures the naturally occurring difference in electrical potential between two electrodes coupled to the ground. Several transport processes can occur in porous media, including coupled processes. With coupled processes, we mean fluxes that are the result of an actuating gradient of a different type than the flux phenomenon itself. There exists an underlying symmetry for the phenomenological coefficients between actuating force and resulting flux, the Onsager reciprocal relations (Onsager, 1931a(Onsager, , 1931b. SP signals can be the result of a variety of these coupled flow phenomena . Often, the electrical double layer (EDL) that exists at the surface between the solid particles (e.g., grain surfaces in a rock formation) and the pore water  plays an important role in the generation of SP signals. As a consequence of the presence of an EDL, electrokinetic phenomena can occur, such as the streaming potential phenomenon, where a hydraulic gradient creates an electric field (Jouniaux et al., 2009;, or the reciprocal electroosmotic phenomenon, where an electric field drives fluid flow. Various electrochemical effects can induce SP signals as well, such as the coupling between chemical gradients and electric current density (diffusion of ions, Fick's law) (Jouniaux et al., 2009;Minsley, 2007), or signals associated with redox processes, for example, related to ore bodies (mineral exploration) and contaminant plumes (Jouniaux et al., 2009;Minsley, 2007). Thermoelectric effects (such as the Seebeck effect) can also occur, where direct conversion from temperature gradients to electric voltages takes place. The EDL has its own contribution to the total thermoelectric coupling (Leinov et al., 2010;Revil et al., 2016).
In addition to the various mechanisms that can generate SP signals, one also has to be aware of variations in electrical conductivity across the study site. Zones of higher electrical conductivity typically tend to attenuate electrical signals, and as such can influence the observed SP signals. However, there are also examples (e.g., at the summit of a hydrothermally active volcano) where a shallow conductor contacting a deep conductive layer provides a current path between regions with different SP values, resulting in an increase of the SP-value around the summit (Ishido, 2004). Ideally, one can complement SP measurements with electrical resistivity tomography surveys, or other geophysical methodologies that inform on the subsurface electrical conductivity structures. In this current work, we have not added such a survey, primarily due to the risk of induced fires in the state park that suffered from dry summer conditions. However, in this study, we do complement the SP signals with another type of passive geophysical measurement, ANSWT, that informs us on the structural geology of the subsurface and as such can help to identify significant heterogeneities in the rock formations.
The fact that naturally occurring subsurface fluid flow forms one of the major sources for measurable SP signals means that SP signals can provide information about subsurface fluid flow paths, thereby enabling us to identify preferential flow paths. In this manuscript, when using the term preferential flow, we do not mean localized discrete flow paths, but rather refer to regions of higher groundwater flow as compared to its surroundings. Many literature studies have demonstrated that SP signals, especially when used in combination with other sources of prior information and/or data types, can provide valuable insight into the local hydrogeology, preferential flow paths, and fluid flow dynamics (such as the Darcy velocity), as well as macroscale structural features of volcanic systems, regardless of their type, and on a variety of spatial and temporal scales (Fagerlund & Heinson, 2003;Finizola et al., 2010;Jardani et al., 2007).
The typical SP data acquisition uses a fixed electrode (the reference electrode) and a moving electrode, separated by a wire of a certain length. Using these two non-polarizable electrodes, combined with a voltmeter, the relative electrical potential difference between the two acquisition points (defined by the electrode positions) will be measured, with the reference electrode artificially chosen to be equal to 0 mV. After a measurement, the moving electrode moves to the next position and again measures the electrical potential difference between that point and the fixed reference electrode. Following this acquisition scheme, we have collected a total of 386 SP measurements at our study site (yellow dots in Figure 2b). The data have been collected at a spacing of 20 m, using 300 m long wires, resulting in individual section lengths of 300 m. Each time we arrived at the end of the section, the moving electrode at its last position became the new reference station for the next section, until the end of the wire was reached again. Again, an artificial value of 0 mV was assigned to this new reference electrode, relative to which all other measurements of that section were made. Two important observations have to be made from this procedure: (1) all SP data points are relative measurements, and we can only infer relative trends that indicate the presence of subsurface anomalies and (2) since we obtain data with several arbitrarily assigned 0 mV values for each new section, the different sections must be properly corrected and combined to obtain continuity of the data set along each profile (the so-called reference corrections). Subsequently, several SP profiles have to be properly combined and corrected for any drift that has occurred during the duration of acquisition (the so-called closure correction), to end up with a properly corrected, global data set for the whole study area. To arrive at this global data set, we have applied a step-by-step SP data processing scheme to the raw data that involves both reference and closure corrections (Barde-Cabusson, Finizola, & Grobbe, 2021;.

Ambient Noise Surface Wave Tomography
Seismic interferometry is a process that enables Green's function retrieval between pairs of recording stations from the correlation of the ambient seismic noise recorded at these stations (Wapenaar & Fokkema, 2006). The retrieved Green function describes the Earth's impulse response that would be observed at one of these recording stations as if there was an impulsive source at the other. Figure 3 schematically illustrates this Green function retrieval process. The correlation of a long time series of ambient seismic noise (black traces) recorded at two stations (seismic nodes) and originating from noise sources along the coast and in urban areas, leads to the retrieved seismic waveform (the Green function) shown in red. This Green function is dominated by surface waves, predominantly Rayleigh waves on the vertical components, traveling between these two seismic stations (let us call these stations A and B). The negative times correspond to the waves traveling from B to A, and the positive arrival times correspond to the waves traveling from A to B. These Rayleigh waves have a depth sensitivity that depends on their frequency content (Aki & Richards, 1980;Socco et al., 2010). The lower the frequency, the larger the wavelength and therefore, the deeper the sensitivity. Because seismic velocities vary with depth in general, surface waves at different frequencies will travel with a different velocity; this is referred to as dispersion. The effect of the frequency-dependent depth sensitivity of surface waves is used to perform seismic tomography of the subsurface, resulting in an image of the three-dimensional distribution of the seismic wave velocities. For Rayleigh waves, these retrieved velocities correspond predominantly to the shear wave velocity.
In this work, we carried out ANSWT making use of 54 MagseisFairfield seismic nodes, 29 1-component and 25 3-component stations, distributed over our study area (Figure 2b). The red dots indicate the deployment locations of the 1-component seismic nodes (i.e., measuring only the vertical component particle velocity), and the blue dots indicate the locations of the 3-component seismic nodes (i.e., measuring both the horizontal component as well as the vertical component of particle velocity). We have continuously recorded ambient seismic noise for five consecutive days, between June 07, 2019 and June 11, 2019. After converting the native data from the seismic nodal system into a SeisComP Data Structure (SDS) filled with daily MiniSEED files for each sensor and each component, we used the open-source software package MSNoise (Lecocq et al., 2014) to process the noise data and compute the cross-correlations between the various station pairs. SEED stands for "Standard for the Exchange of Earthquake Data," and is a data format intended primarily for the archival and exchange of seismological time series data and related metadata. MiniSEED is a stripped-down version of SEED and contains only waveform data; no station and channel metadata is included. The computed cross-correlations resulted in the retrieval of various Green functions between the different station pairs. For this work, we only used the vertical component measurements (i.e., measuring the vertical component particle velocity movement of the Earth due to seismic wave propagation). The spectrograms from two stations (KA.115, located far from the highway and KA.332, located close to the highway, Figure 2b) are shown in the Figures S1 and S2. This array provided a total of 1431 pairs of stations and, correspondingly, 1431 vertical-to-vertical (ZZ) correlations.
Theoretically, for optimal, symmetric, noise-based Green function retrieval, we require a uniform distribution of noise sources, 360° around the seismic stations (Wapenaar & Fokkema, 2006). Any deviation from this distribution will result in a suboptimal retrieval of the Green functions between the station pairs, forming a potential source of uncertainty in the final inferred seismic tomography images. To assess the distribution of the ambient noise sources around our field array, we performed a plane-wave beamforming analysis (Gerstoft & Tanimoto, 2007)   . Schematic cross correlation-based Green's function retrieval. The correlation of the seismic noise recorded at two stations A and B (black traces) results in a seismogram that represents the surface wave propagating between the two stations (red trace; A to B corresponds to the positive times, B to A corresponds to the negative times). This retrieved Green function, for this particular pair of stations, is laterally sensitive to the area indicated by the red ellipse. By cross-correlating different pairs of stations, different areas of the study site will be covered. The depth sensitivities of these surface waves are schematically displayed by the black curves in the cross section. The solid curve illustrates a typical short-period sensitivity; the dashed curve represents the sensitivity curve for a longer-period. The incoming seismic noise (black-wiggled arrows) is generated by the ocean and along the coast for the longer periods (above 1 s), and in urban areas, along roads (traffic noise) and trails for the shorter periods (below 1 s).
( Figure 4). To mitigate the effect of potentially high-amplitude signals, we normalized the filtered data for every 1 h window by its maximum amplitude. The 1 h long beamforming diagrams were then averaged for the whole recording period. The lower frequency bands correspond primarily to ocean noise, whereas noise from traffic and urban areas typically shows up in the higher frequency bands (∼5-10 Hz). However, in Figure 4, we see that for all frequency bands, the noise is mostly coming from the direction of the ocean and the coast lines (east and northeast), in agreement with the dominant swell in Oahu (S. Vitousek & Fletcher, 2008). At higher frequencies (>3.5 Hz), the energy of the noise seems to be correlated with the topography of the coastline, being stronger for sharp cliffs. At 0.25-0.5 Hz, due to the large wavelengths, our array does not provide the resolution anymore to determine from where, and with which velocity, the seismic waves are coming. For all bandwidths, there is little noise coming from the direction of the mountains and the center of the island (north-west), which would correspond more with the urban and traffic noise source directions. The apparent velocities (inverse of the slowness) are around 1 km/s and decrease with increasing frequency.
As expected from the noise distribution analysis, the correlations are fairly symmetrical for frequencies between ∼1 and 2.5 Hz ( Figure 5, left panel) and more asymmetric for higher frequencies >5 Hz ( Figure 5, right panel). The signal-to-noise ratio of the correlations is high for frequencies below 10 Hz. The aperture GROBBE ET AL.  . Plane wave beamforming analysis to assess the noise-source distribution in our study area. The warmer the colors (yellow-red), the higher the intensity of the noise from that direction. The polar plots show the north in the central-top, east central-right, south central-bottom, and west central-left of the diagram. Different panels display the beamforming analysis for different frequency bandwidths. The values on the axes correspond to the apparent slowness values (1/apparent velocity). We can observe that for all bandwidths, most of the noise comes from the ocean, east and northeast direction. of the array (total extent) is about 1 km. This aperture limits the lowest frequency (largest wavelength) that we can use to ∼1 Hz. Below 1 Hz, the wavelengths of the Rayleigh waves are too large for the inter-station spacing to enable accurate measurements of the dispersion curves.
We automatically picked the group velocity dispersion curves using a frequency-time analysis method (Ritzwoller & Levshin, 1998) with user supervision and manual rejection of the low-quality curves (Mordret et al., 2015). In total, we managed to pick 1055 dispersion curves between 1 and 10 Hz (Figure S3). The number of dispersion curves is not equal for each frequency and decreases with decreasing frequency, reducing the lateral resolution of the tomography at lower frequency. The group velocity dispersion curves were then regionalized into group-velocity maps at discrete frequencies using a straightray linearized inversion scheme (Mordret et al., 2013(Mordret et al., , 2015 on a Cartesian grid of 33 × 33 m grid cells. Each map was damped to the average velocity at the inverted frequency where the ray density per cell is lower than 10. A Gaussian smoothing was also used to regularize the inversion and avoid overfitting the data. The correlation length of the smoothing window was 45 m. The group velocity maps are shown in Figure S4. Each cell of the set of group velocity maps for every frequency describes a local group velocity dispersion curve. These dispersion curves can be inverted into 1D shear wave velocity versus depth profiles. For the depth inversion, we followed the approach of Mordret et al. (2015) by parameterizing the generic 1D shear wave velocity profile versus depth with a power-law backbone velocity profile overlain by four cubic splines to account for higher or lower velocity layers at depth. This parameterization ensures a smooth variation of the velocity with depth. The 1D velocity profiles were discretized with 30 layers above a half-space. We also inverted for the velocity of the half-space. In total, we inverted for seven parameters. The inversion was carried out using a Neighborhood Algorithm (Mordret et al., 2015) which helped to efficiently explore the parameter space by refining the random sampling around the most promising parameters, leading to the best fit to the data. Among the sampled models (8,000 for each geographical point of the grid), the best ones were chosen based on the misfit between the local dispersion curve extracted from the group velocity maps, and the synthetic dispersion curve forward modeled from the sampled 1D shear wave velocity model. Examples of the local inversion of group velocity dispersion curves into depth versus 1D shear wave velocity models are shown in Figures S6-S10. Their locations are shown in Figure S5  Hz, the right panel shows the bandwidth 5-10 Hz. Ideally (for a perfectly uniform noise distribution), the retrieved Green functions are symmetric around t = 0 s (horizontal axes). In conjunction with Figure 4, we can observe that for the bandwidth 1-2.5 Hz (estimated 1-5 Hz), the retrieved waveforms are reasonably symmetric, whereas we can observe an asymmetry for the bandwidth 5-10 Hz because the dominant high-frequency noise sources are mainly located East of the array. the average 1D models were then laterally interpolated to generate the final 3D shear wave velocity model, which is described in the next section. Figure 6 shows the final corrected global SP data set. It represents a map of SP-anomalies at our study site. As described in the Methods section, the SP data set has been corrected for reference and closure. The SP values are thus relative to a unique 0 mV reference taken at the sea, by convention (Barde-Cabusson, Finizola, & Grobbe, 2021). The blue colors represent relative negative SP values, the red colors represent relative positive SP values. The white dots indicate our SP measurement locations. Remember that we interpret our SP signals as being the result of the electrokinetic streaming potential, where a hydraulic gradient creates an electrical field. Negative SP anomalies are associated with flow away from the measurement electrodes, for example, fluid infiltration. A hydraulic gradient that causes groundwater flow always has a downward com-GROBBE ET AL.

Self-Potential
10.1029/2020WR029057 11 of 21 Figure 6. Self-Potential (SP) map of the study area, with hydrological interpretations, overlain on a 1 m Digital Elevation Model (DEM) from Hawaii Statewide GIS Program, https://geoportal.hawaii.gov/. The colorbar indicates the relative SP values. White dots represent each SP measurement, white arrows represent the inferred major groundwater flow paths, and SGD stands for submarine groundwater discharge. ponent (infiltration) and as such will show up as a negative anomaly. We can clearly observe a pronounced zone of negative SP values that we interpret as a zone of groundwater flow.
We have to bear in mind that all SP measurements are relative measurements, and that one should therefore focus on the variations in the observed SP signals as opposed to the absolute value of the SP anomalies. However, given that our data set has been referenced to the sea, which is considered a 0 mV equipotential surface, the area associated with large negative SP values as observed on our map can be considered as an area of major infiltration occurring along the path of underground water flowing from land to ocean. Thus, we associate the relative negative values (or lower SP values) with zones of preferential flow versus SP values higher than the 0 mV reference, representing zones of little to no flow. We can observe positive SP values with a maximum of 15 mV (mostly about 10 mV on the whole map), which rules out the presence of major uprising flows. Those values are likely associated with background noise of our study area and/or local shallow uprising flows created by small-scale obstacles in the ground. The resulting SP signal is not significant with respect to the observed negative values associated with infiltration processes (with SP values down to approximately −55 mV). Similar interpretations of the relative SP data in terms of fluid flow can be found in the literature, such as for hydrothermal systems where up-welling results in positive SP values, and infiltration and downward fluid motion in negative values (Revil et al., 2008), or in studies of groundwater flow (Revil et al., 2005).  inverted down to 1 km depth, the resulting models are only reasonably resolved (with uncertainties smaller than 150 m/s) down to about 400 m below the surface. The lateral resolution is about 150 m at 40 m depth and decreases to about 500 m at 300 m depth. It is mainly controlled by the density of rays available at each frequency.

Ambient Noise Surface Wave Tomography
On the map at 20 m depth, we clearly see a smearing of the anomalies in the south-west/north-east direction, because at high frequency, we were only able to pick data for pairs oriented in this direction, which is also the direction of the main source of high-frequency noise (Figure 4). Deeper, the observed structures are less affected by this problem because the azimuthal distribution of rays is more uniform at lower frequencies.
We observe four main areas of interest in the seismic model from north to south and from east to west. First, we observe a shallow low-velocity zone at the northern tip of the study area that attenuates below 160 m and becomes a high-velocity anomaly. Then, just south of it we see around 80-240 m depth an east-west trending high-velocity anomaly that reaches the southern tip of Ridge 2 to the west. South of this previous anomaly, is a large prominent low-velocity zone that reaches the coastline at shallow depths and penetrates toward the west at deeper depths. Finally, to the west of the studied area, we observe a broad high-velocity anomaly in the shallower part of the model that tends to remain visible only in the south-west corner of the studied area below 80 m depth and that vanishes at depths greater than 200 m.
To gain a better sense of the geometry of the structures at depth, vertical cross sections through the model are shown in Figure 8. They allow us to see the general increase of velocity with depth, as well as the geometry of a prominent interface around 100-200 m depth, represented by the isovelocity 1,400 m/s (Figure 8). The topography of this interface is shown in the lower half of Figure 9. The parameterization used for the depth inversion imposes smooth variations of the shear wave velocity with depth. Therefore, sharp interfaces are represented by fast but smooth increases of velocity rather than step-like discontinuities. This choice of parameterization also results in a vertical smearing of the structures in the model and an overestimation of the thickness of these structures.

Hydrogeological Joint-Interpretation
From the SP map, we can infer several preferential flow paths as indicated in Figures 6 and 9 (as discussed, we look at relative differences and can associate relatively negative values with zones of preferential flow). A major fresh groundwater flow path feeds from the Ridge 1-Ridge 2 area into a brackish water mangrove GROBBE ET AL.  bay at the coast. This is indicated as Submarine Groundwater Discharge (SGD) in Figure 6. Furthermore, we can observe some minor flow paths to the south-east near Ridge 3.
Let us now combine this knowledge with the seismic information. In general, seismic velocities are sensitive to different properties of the subsurface: they can represent variations in lithology, with certain rock types being intrinsically faster or slower than others. They can also represent differences in alteration of the rocks, water content, degree of fracturing, clay content, porosity, or pore aspect ratios. These factors should be kept in mind when interpreting seismic tomographic images, as one single seismic velocity value could correspond to different aspects of the nature of the subsurface.
The shear wave velocities larger than 1,300-1,400 m/s in the 3D model likely represent the Koʻolau basaltic bedrock of the old volcanic edifice, with the highest velocities corresponding to remnants (or subsurface extent) of ridges. The high-velocity anomaly linking Ridges 2 and 3 appears to be an actual extension of the old edifice bedrock below the surface and is a topographic high of the bedrock topography.
The shape of the bedrock topography at a depth around 100 m below the surface (isovelocity of 1,400 m/s) is probably an expression of the erosion of the ancient Koʻolau basaltic edifice during the Pleistocene and represents an unconformity (Figure 2). This erosional surface may consist of saprolite, which typically has a very low hydraulic conductivity. This also means that this surface can act as a seal, keeping groundwater above it separated from water below it.
The low-velocity anomaly filling this paleo-channel basin corresponds to Pliocene sediments and eruptive products from the Honolulu Volcanics group (Stearns & Vaksvik, 1935   fill can reflect their alteration due to preferential ground water flow paths or simply water-saturated formations (resulting in lower S-wave velocities), or a combination of the two. Especially given the fact that the study site is located close to the coast, we expect salt water-saturated formations to be present. The low-velocity anomaly at the northern tip of the study area is a small sedimentary basin (Figure 2a and section CC" in Figure 8), whose creation could have been facilitated by the transverse high-velocity structure to the south of it.
The data integration (seismic structural information with SP fluid flow information) highlights that the basaltic bedrock of the ridges forms the dominant groundwater flow unit, as opposed to the valley fill material, and that groundwater flow may occur at deeper levels spatially coinciding with the minor valley in the western ridge complex. Furthermore, the joint interpretation of the SP and seismic data strongly favors the main groundwater flow taking place in the deeper parts of the identified paleo-channels (Figure 9) of the erosional surface at an isovelocity of 1,400 m/s, at roughly 100-200 m depth. Complementing SP data with seismic data thus enables us to place the groundwater flow inferences, as made from the SP data, at depth, albeit with associated methodological uncertainties.

Discussion
It is important to realize that the SP method is a potential field method, meaning that the anomalies in electrical potential as measured at the surface are theoretically the result of an integration of the effects of subsurface fluid flow processes over a whole range of depths. This also means that, without any other (prior) information (geological, geophysical, borehole, or well data, etc.), it is difficult to reliably relate the observed anomalies to processes occurring at a certain depth in the subsurface. Matched bandpass filters (Phillips, 2001;Van Noten et al., 2015) are one way to estimate the depth of potential field sources directly form the data, but the method require a fairly regular spatial sampling of the field to remain unbiased. The field constraints at our study site prevented us to regularly cover the area with SP measurements therefore hampering an accurate estimation of the SP signal sources at depth. In this work, we have therefore integrated SP data with ANSWT data, to provide valuable additional information about the subsurface hydrogeological properties, enabling us to relate the observed SP anomalies to fluid flow processes occurring at depth. The results of this work, and the associated hydrogeological interpretation as presented in the previous section, have clearly demonstrated the benefits of using these two complementary data types. However, it is important to realize that even when using such a joint interpretation framework, there exists a certain degree of non-uniqueness and uncertainty, and alternative hydrogeological interpretations are possible.
The non-uniqueness and uncertainty are caused by a variety of aspects and their primary source is often methodological in nature. Besides possible measurement uncertainties and model uncertainties, various assumptions and choices are made in the processing of the acquired data. The SP data processing sequence (Barde-Cabusson, Finizola, & Grobbe, 2021) is relatively straightforward and applies several corrections to the data, such as reference and closure corrections, to mitigate the effects of changes in measurement conditions during acquisition, as well as other data acquisition inaccuracies, as much as possible. To go from point SP measurements along several profiles to a full SP map as presented in this work, an interpolation step is required, with intrinsic uncertainty. The ANSWT data processing is more complex. Potential sources of uncertainty are, for example, imperfect azimuthal noise distributions for theoretical Green function retrieval, various filtering and smoothing procedures and assumptions, user-supervised dispersion curve picking, and two inversion steps with inherent non-uniqueness: a straight-ray linearized inversion to create group-velocity maps from group velocity dispersion curves, and the inversion of these maps into 1D shear-wave velocity-depth profiles. These steps are followed by a lateral interpolation step to generate the final 3D shear wave velocity model, again with intrinsic uncertainty. One way to decrease both uncertainty and non-uniqueness is to incorporate several seismic attributes in the inversion. Here we only used the vertical component of the seismic data to obtain the highest lateral resolution possible, using every stations. But it is possible to use the horizontal components as well to measure, for example, Love-wave dispersion (Chmiel et al., 2019), Horizontal-to-vertical spectral ratios (Spica et al., 2018), or Rayleigh wave ellipticity (Berg et al., 2018). The Love wave information would help to constrain the velocity structure while Horizontal-to-vertical spectral ratio or the ellipticity information would better constrain the depth of the layer interfaces in the subsurface. We will devote future studies to assess the potential benefits of adding horizontal components information for hydrogeological applications.
A second source of uncertainty is an inherent consequence of the physical information content and sensitivity of the methodologies used. Earlier, we have already discussed the variety of physical mechanisms that can cause and influence SP signals, such as thermoelectric effects, electroosmosis, and electrochemical processes. This highlights that care has to be taken when interpreting SP signals. In this work, we interpret the SP signals to be predominantly the result of the electrokinetic streaming potential, where a hydraulic gradient creates an electric field, and as such we relate the observed signals to subsurface fluid flow phenomena. This assumption is very common when dealing with SP signals (Jouniaux & Ishido, 2012). We can make several observations that justify this assumption. First of all, no dominant thermoelectric effects are to be expected since we are working on a nonactive volcano (Koʻolau volcano), and recent geothermal assessments indicate that no significant heat flow is to be expected (Ito et al., 2017). Electroosmosis is not expected, since we do not inject any electrical current in the area. There might be some electrochemical processes that can influence the SP signals, especially in the coastal area where salt water underlies fresh water at relatively shallow depths, with a mixing zone of brackish water in between. This might produce SP effects that remain difficult to quantify. If we assume these effects to primarily cause random perturbations in this area, these effects should result in an elevation of the noise level of the SP signals, and might explain some of the shorter wavelength features that are observable in the SP map. Alternatively, they may also affect the total background SP signal over the whole study area in a more or less uniform fashion. Either way, it is highly unlikely that the observed major negative SP anomalies that we relate to subsurface flow can be explained by such electrochemical processes. No known ore bodies or contamination sites are present in the area, so redox processes are not expected to be a dominant process. As discussed, the SP signals may be influenced by subsurface variations in electrical conductivity related to, for example, signal attenuation. Nevertheless, several of our field observations, such as negative SP anomalies coinciding with visibly wetter areas, with the presence of surface water streams, or with vegetation patterns that are indicative of the presence of fresh water or SGD (like mangroves in the coastal zone), show that the observed SP signals correlate well with fluid flow. Moreover, by complementing the SP data with ANSWT results, we further enhance our confidence that the electrokinetic streaming potential is the dominant SP-generating mechanism for our study site.
Another source of uncertainty is related to the fact that SP data, as mentioned, provides information about electrical potential differences as measured at the surface, due to subsurface fluid flow effects integrated over a whole range of depths. It is therefore difficult to pinpoint with absolute certainty the depth of the dominant fluid flow processes causing the anomaly, and very likely, the anomalies are the result of several co-existing processes. Nevertheless, combined with other types of information, we can significantly narrow down the range of potential interpretations, as has been successfully demonstrated with other data types than ambient noise seismic in many previously published SP studies (see earlier references). The ANSWT data provides information on the seismic velocity of the medium, which is a function of the hydrogeological properties of the rock. Hence, we can infer information about the structural geometry and properties of the medium, and to some degree information about the degree of fluid saturation. However, different rock properties have different effects on the seismic velocity. Consequently, the relationships between rock properties (such as porosity, permeability, degree of fluid saturation, etc.) and seismic velocities are often complicated and nonlinear, impeding the direct inference of, for example, hydraulic conductivity. What we propose, and have successfully demonstrated in this manuscript, is to relate the inferred structural geological features with inferred zones of preferential fluid flow, to identify clear correlations between fluid flow and hydrogeological structures controlling these flow patterns.
With the above information, let us revisit the results of Figure 9 and discuss an example of potential non-uniqueness in our hydrogeological interpretation. Based on prior literature information (that dominant groundwater flow occurs in the basaltic bedrock and ridges on the leeward side of Oʻahu) combined with the integration of the SP and ANSWT results, we have identified that the basaltic bedrock of the ridges indeed likely forms an important groundwater flow unit in our study site. Furthermore, we have correlated the negative SP anomaly (indicating zones of higher flow/infiltration) that covers the zone between Ridges 1 and 2 (see Figure 2 for the ridge numbers) with the location and geometry of the identified paleochannel at the erosional surface of the basaltic bedrock, from which we infer the dominant groundwater flow patterns at depth. However, given the aforementioned uncertainties and non-uniqueness, particularly the fact that SP anomalies are the net effect of fluid flow properties over a range of depths, we cannot exclude any fluid flow in the shallower regions of the sedimentary infill between Ridges 1 and 2, as well as parts of the main valley fill between Ridges 2 and 3 and in the coastal zone. Also, if the identified erosional surface of the Koʻolau volcano indeed acts as a seal for fluid flow over the whole area, separating groundwater flow above this interface from groundwater flow below it, it may cause a separation of fresh and salt water, without the presence of a brackish water zone. Direct well measurements would be required to confidently make such an interpretation.
In short, the main cause of uncertainty is that both the SP and ANSWT data are indirect measurements of the subsurface with a variety of intrinsic methodological assumptions and limitations. Ideally, one would like to constrain the SP fluid flow inference and seismic inversion steps with direct measurements from, for example, borehole data or water well head-level information. Unfortunately, there are no available water wells in this study site to provide a ground-truth calibration to our model. Nevertheless, despite these uncertainties and sources of non-uniqueness, which are inherent to any (hydro)geophysical methodology, we can have confidence in our most likely hydrogeological interpretation. Both the SP and the ANSWT methodologies are well-known, widely accepted, and proven methodologies to infer fluid flow patterns and seismic velocity structures, respectively (see, e.g., earlier references). Furthermore, there is a clear correlation between the fluid flow patterns and the geometry and location of the paleo-channels in our study site. Additionally, we have field-based observational evidence for the rough location of discharge of fresh groundwater at the coast; it feeds into a brackish water coastal mangrove zone. Finally, it is well known from numerous literature studies (see introduction) that for watersheds on the leeward side of the island of Oʻahu, the dominant aquifer flow occurs in basaltic ridges and bedrock, that have relatively high hydraulic conductivity compared to the valley fill material that has lower hydraulic conductivities, a pattern that is consistent with the acquired SP and ANSWT data for our study site. As the next steps, similar joint SP-AN-SWT studies can be carried out in other valley-ridge complexes, both on the leeward and windward sides of Oʻahu, to investigate the similarities and differences between valley-ridge systems. Ideally, a study site with available ground-truth data, such as hydrological well data, combined with hydrological modeling results, should be prioritized, to enable validation and verification of the geophysical inferences.

Conclusions
Using a joint interpretation framework of ANSWT data and SP data, we have identified the regions of dominant groundwater flow in the old stream valley at the Kaiwi Coast, Oʻahu, Hawaiʻi. The SP data display the distribution of groundwater in the region and zones of higher infiltration and/or flow rates. Using ANSWT, we have identified the hydrogeological erosional contact between the older Koʻolau volcanic basaltic formation and the younger post-erosional Honolulu volcanics, as well as the distribution of valley fill alluvium, Honolulu volcanics, and sedimentary water-saturated deposits near the coast. The integration of SP and seismic data suggests that significant groundwater flow occurs in the deeper parts of the identified paleo-channels on the erosional surface between the Koʻolau basaltic bedrock and the overlying Honolulu volcanic and sedimentary formations, roughly at 100-200 m depth. Additionally, our results are consistent with the literature-based hypothesis that, for watersheds on the leeward side of the island of Oʻahu, the dominant aquifer flow occurs in basaltic ridges and bedrock with higher hydraulic conductivity than the valley fill material. However, due to the inherent limitations of the hydrogeophysical methods, and the lack of direct measurements such as well as data in the area, our study cannot exclude other regions of potential groundwater flow, for example, in the shallower regions of the valley fill material. Combining structural and saturation sensitive seismic ambient noise data with pore fluid (groundwater) flow-sensitive SP data opens up a way to place flow inferences, as obtained from SP data, at depth, thereby providing complementary insight into complex subsurface hydrology. It is a passive geophysical, integrated methodology with low environmental impact that will be valuable for studying any watershed system on Earth.