Supplement of Identifying recharge under subtle ephemeral features in a flat-lying semi-arid region using a combined geophysical approach

Abstract. Identifying and quantifying recharge processes linked to ephemeral surface water features is challenging due to their episodic nature. We use a combination of well-established near-surface geophysical methods to provide evidence of a surface and groundwater connection under a small ephemeral recharge feature in a flat, semi-arid region near Adelaide, Australia. We use a seismic survey to obtain P-wave velocity through travel-time tomography and S-wave velocity through the multichannel analysis of surface waves. The ratios between P-wave and S-wave velocities are used to calculate Poisson's ratio, which allow us to infer the position of the water table. Separate geophysical surveys were used to obtain electrical conductivity measurements from time-domain electromagnetics and water contents from downhole nuclear magnetic resonance. The geophysical observations provide evidence to support a groundwater mound underneath a subtle ephemeral surface water feature. Our results suggest that recharge is localized and that small-scale ephemeral features may play an important role in groundwater recharge. Furthermore, we show that a combined geophysical approach can provide a perspective that helps shape the hydrogeological conceptualization of a semi-arid region.



Travel-time Tomography Processing
To ensure that we collected the first arrival the data in the field we used a 10 ms buffer prior to the sledgehammer trigger (t0). The data were collected with 5 m geophone spacing and 5 m shot spacing. The high shot spacing density was selected because, although not shown here, the data were set up for a reflection survey (not discussed in this manuscript) which requires more shots to build common midpoint fold. The profile was running parallel to a well-traveled road connecting Adelaide to Port Wakefield. The highway is approximately 100 m south west of the line. This highway caused unpredictable low-frequency noise on various shots throughout the day ( Figure S1).
To reduce the effect of traffic noise we stacked each shot eight times. We also removed a noisy trace at 95 m. First arrivals were picked manually on the stacked shot files (Figures S1). Prior to picking each shot was trace normalized. We attempted to pick as many first-arrivals without the use of a band-pass filter to prevent any shifts caused by the filter at early times. After the arrivals were picked on the raw shot gathers we applied a band-pass filter between 10-120 Hz, which significantly reduced the low frequency traffic noise and allowed us to pick faroffset arrivals ( Figure S1). Despite the traffic noise, we were able to confidently pick 2014 travel-times (87% of the data) ( Figure S2a).
The travel-time picks were high-quality. The quality of the pick is based on the difference between reciprocal travel times. A reciprocal travel-time occurs when a shot at station x and a receiver at station y is picked, and then later in the data set the shot occurs at station y and the receiver at station x. Theoretically, the difference between these travel-times should be zero, or when plotted they should lie along a 1:1 line ( Figure S2c). In our case, because of the high density shot spacing, almost every pick has a reciprocal travel-time and all the reciprocal picks fell within 5 ms of each other ( Figure S2a).
We use the reciprocal values to drive the error weight matrix in the PyGIMLi inversion (Rücker et al., 2017). Thus, an error be assigned to every pick. In most cases we used the reciprocal error plus 1.2 ms (refwhy?). In locations where we did not have a reciprocal travel-time we applied a linear function as a function of offset (distance from source). In this case error (in seconds) was defined by Eq. S1 (ref). where offset is simply the distance from the source to the receiver, and the spread length was 235 m. We use Eq. S1 and the calculated reciprocal values to ensure that each pick has an assigned error ( Figure 2b). It also provides a way to determine if the data are over fit or under fit. These errors are read into the inversion and used to weight the travel-times during the inversion as well as calculate the χ2 fit: Where is the observed travel time, is the modeled travel time, is the picking error, and N is the total number of picks.
With the travel-time picks and associated errors collected we can invert the model. The mesh was selected so that there were at least three nodes between each geophone and the maximum triangle size was no more than 2 m. Based on the slopes of the first arrivals we determined the average velocity was around 2000 m/s. Thus, we used a gradient starting model with 400 m/s at the surface and 2000 m/s at 40 m depth. The depth of the model was selected to not limit the ray paths. This took some experimentation, but the final mesh extends a few meters below the lowest ray path. Smoothing parameters were manually selected through trial and error to achieve a feasible geologic model and still achieve a good model fit. Specific settings in the the PyGIMLi refraction inversion we set Lambda equal to 200 and the z-weight to 0.5 and the maximum iterations to 5.
Throughout the manuscript we interpret small changes in velocity. Although it is difficult to say exactly how accurate these velocities are, we did attempt to quantify velocity uncertainty. To quantify uncertainty, we took the 2014 travel time picks and randomly removed 30% (604 data points). We used this decimated data set to run an inversion. We would then remove another random 30% of the original data set and invert. This process was repeated 50 times. Each velocity model was saved. At the end all 50 models were averaged together. The average model was run through the forward model to get the travel-times and the error fits. This averaged model is the one shown in the manuscript (Figure 2a). This means that at each model parameter we have a standard deviation, indicating how much it changed ( Figure S3c). The standard deviations highlight areas of the highest uncertainty. The values never varied more than 200 m/s ( Figure s3c). We are thus confident in the P-wave velocity profile. The final model also outputs the ray coverage. This is the total number of rays that cross through a single model parameter ( Figure S3b). This is how we defined the masked region in the main velocity profile ( Figure 2). The final model had a χ2 of 0.82 and a RMS value of 1.73 ms. Both indicate that the modeled travel-times reasonably fit the observed travel times ( Figure s3a).

Surface wave Processing
To construct the 2D S-wave velocity profile shown in the manuscript ( Figure 3) we used an open source code written in Matlab called Surface Wave Inversion and Profiling (SWIP) (Pasquet and Bodet, 2017). In this section we provide our specific geometry and figures showing the high-quality nature of data to provide a sense of confidence in the final S-wave velocity profile ( Figure 3). SWIP uses spatial windowing and frequency stacking to extract local dispersion images representative of 1D slices of the investigated medium, thus detecting more lateral heterogeneity while increasing the signal-to-noise. In this case we extracted dispersion images across 8 traces ( Figure S4a). The assumption is that the velocity structure under these 8 traces is horizontally layered and then assigned to the midpoint of the traces. To avoid near-source effects we ensured that there was at least 25 m (5 traces) between the source and the closest geophone. We allowed the source to be up to 115 m away (24 traces) from the first geophone. The code sorts valid pairs of shots and traces through all of the seismic data, then computes dispersion images for each valid shot gathers, and finally stacks all dispersion images corresponding to a single midpoint. On average we stacked ~40 dispersion images per midpoint which generated high quality and easily identifiable dispersion curves ( Figure S4b and S4d). Even though we used 14 Hz geophones we had clear energy down to at least 10 Hz ( Figure S4b and S4d). The lower frequencies could have been generated from the traffic on the near-by highway.
Based on the shot and geophone geometry we extracted and picked dispersion curves for 41 midpoints every 5 m beginning at 17.5 and ending at 217.5. To show that the low velocity zone between 60-80 m was a clear signal we show two dispersion curves, one from 62.5 m and the other from 162.5 m ( Figure S4). During picking it was clear that there were two modes visible ( Figure S4). We picked both the fundamental and first higher order mode ( Figure S4). It has been shown that higher-order modes can stabilize the inversion and provide additional depth penetration (Xia et al., 2003).To visualize the dispersion picks, they are converted to wavelength using the phase velocity and frequency (λ=v/f). The picks can be plotted in space for each midpoint ( Figure S4e and S4f). The dispersion curve picks for both the fundamental and higher order first mode show a clear pattern and consistency from midpoint to midpoint ( Figure S4e and S4f). It should be noted here that the low velocity zone between 60 and 80 m is clear in both modes of the data, and thus should be reflected in the final S-wave velocity profile ( Figure S4e and S4f).
The inversion of the surface wave data consists in independent (not laterally constrained from midpoint to midpoint) 1D inversions. The dispersion curve picked at each midpoint are used to in a Monte Carlo inversion scheme. The inversion does three runs of 15000 models each. For each midpoint, the best fitting 1000 S-wave velocity models that replicate the dispersion curves are kept ( Figure S5). Thus, for each midpoint we have a suite of models that all converge around a narrow range of a 1D velocity profile ( Figure S5a and S5c). The final 1D velocity model is selected by averaging the velocity of top 1000 models. This 1D profile is then used to calculate the theoretical dispersion curves. The modeled dispersion curves fit within our error bars in most cases ( Figure S5b and S5d). The difference between our picks and modeled dispersion curve can then be used to calculate a misfit ( Figure  S5g). A misfit value below 0.5 is considered adequate. Furthermore, we can visually look at the difference between our picked dispersion curves and modeled dispersion curves ( Figure S5e and S5f) to see that the model is representative of the data. The model fits provide confidence in the final 2D velocity model in the manuscript (Figure 3).

Modeling Electrical Conductivity with Archie's Law
The electrical conductivity of the formation (the one that we measure) can be modelled using Archie's Law (Archie, 1941;Robinson et al., 2012).
In Eq. S3 φ is porosity, S is saturation, Rform is the formation resistivity, Rfluid is the fluid resistivity, and conductivity is the inverse of resistivity. The variable m is the cementation coefficient, which accounts for effects of permeability and tortuosity. This value varies between 1.2 and 4.4 (Lesmes and Friedman, 2005;Robinson et al., 2012). The variable n is an empirically defined coefficient but is commonly set at a value of 2 (Knight, 1991;Day-Lewis, 2005;Robinson et al., 2012). Using Archie's equation (Eq. 2) we can explore how the conductivity of the formation will respond as a function of fluid saturation and the conductivity of the fluid. It is often assumed that as the water saturation increases, the fluid conductivity will also increase. We used Archie's law to show that it is possible to observe a drop in electrical conductivity if the pores are being saturated with a more resistive fluid or if a less conductive fluid replaces the higher conductive fluid originally in place.
Using the fully saturated volumetric water content below the water table from the downhole NMR data provides an estimated porosity of 0.25 m3/m3 (Figure 5d). Assuming a soil density of 1.5 provides a porosity estimate of 0.23 m3/m3 from the soil samples collected, which is consistent with the NMR results (Figure 5a). Using equation 2, we calculated the electrical conductivity assuming a fixed porosity of 0.25 m3/m3, setting m to 1.3, setting n to 2, and varying the fluid conductivity from 500 to 4000 mS.m-1 and a saturation from 0 to 1 m3/m3 (Figure 7a).
There are three end-member cases that can be observed from the calculation of the formation conductivity ( Figure S7a). The first is filling the pores with water that has the same electrical conductivity as the groundwater ( Figure S7a). The result is that the formation conductivity rises exponentially ( Figure S7b) and this is consistent with the common interpretation that increasing the saturation also increases the electrical conductivity. The second case is filing the pores with saltier water; that is the fluid filling the pore-space has an increasingly higher electrical conductivity. In this case the formation conductivity also rises exponentially, but at an even faster rate than just filling the pores with water ( Figure S7b). The last end member case is filling the pores with fresher water. That is as the saturation increases the fluid conductivity decreases. In this case, it is possible to get a small drop in electrical conductivities ( Figure S7b). Furthermore, from this relationship it can be observed that at the same saturation level, the formation conductivity will reduce if the fluid in the pore space is replaced with a more resistive fluid ( Figure  S7a). This basic forward modelling      Figure S4e) and the modeled dispersion curve for the fundamental model. (f) The difference between the observed picks ( Figure S4e) and the molded dispersion curve for the first higher-order model. (g) The sum of the differences for both modes plotted spatially for each individual inversion.

Figure S6
Results from the downhole NMR sounding at the Geoprobe location (Figure 1). Panels a-c show the measured decays in black and the multi-exponential model fits in gray.  (a) The relationship between saturation, fluid conductivity, and formation conductivity calculated from Eq. 2, assuming a constant porosity of 0.25, m=1.3, and n =2. The solid black line represents a case where the pore space is being filled with water that gets less conductive (i.e. fresher). The dashed line is a case where water is being replaced with water that gets more conductive (i.e. more saline). The dotted/dashed line is the case of filling the pores with water of the same conductivity. (b) The formation conductivity as a function of saturation for the three end-member cases shown in panel a.