A hydrogeologic conceptualization is critical to understand, manage, protect, and sustain groundwater resources, particularly in regions where data are sparse and accessibility is difficult. We used airborne electromagnetics (AEM), shallow seismic reflection and refraction, and downhole nuclear magnetic resonance (NMR) logs to improve our understanding of an arid groundwater system influenced by palaeovalleys. We found that there is a limited connection between the palaeovalley and fractured bedrock aquifers because they are separated by a spatially variable layer of saprolite, which is the layer of chemically altered rock on top of the fractured bedrock. The AEM data provided an estimate of the top of the saprolite but failed to effectively image the bottom. In contrast, the seismic data provided an estimate of the bottom of the saprolite but failed to image the top. This geophysical combination of electrical and seismic data allowed us to map saprolite thickness in detail along a 1.7 km long transect that runs perpendicular the main trunk of a well-defined palaeovalley. These data indicate that the palaeovalley is lined with a heterogeneous layer of saprolite (approximately 3–120 m thick) that is thickest near its edges. Despite the observed variability, only a small percentage of the fractured bedrock aquifer (8%–17%) appears to be in contact with the palaeovalley aquifer. Furthermore, the lack of an elastic boundary at the top of saprolite suggests that the porosity of the saprolite is similar to the palaeovalley sediments — an observation that is supported by the downhole NMR-derived water contents. The electrical change at the top of saprolite is caused by a combination of a decrease in total dissolved solids of the groundwater in the saprolite and a change in pore structure associated weathering in situ versus transported weathered materials. The presence of saprolite, which commonly behaves as an aquitard, may limit the groundwater exchange between the palaeovalley and bedrock aquifers, with implications for the regional groundwater resource potential.

Successful exploration and sustainable management of groundwater requires a comprehensive conceptualization of the local groundwater system: the set of hypotheses that describe the hydrostratigraphy, system boundaries, variations in hydraulic properties, and system fluxes in and out of the system. The uncertainty in conceptualization can be considerable, particularly in remote areas with limited access, primarily because these settings are also characterized by limited data. Geophysical methods can constrain aquifer geometry (Cardimona et al., 1998; Holzschuh, 2002; Chen et al., 2010; Enemark et al., 2020), identify contrasts in salinity (Franco et al., 2009; Viezzoli et al., 2010; Kirkegaard et al., 2011), and identify processes driving the exchange of fluxes (McClymont et al., 2011, 2012; Travelletti et al., 2012; Valois et al., 2016).

In ancient landscapes, such as those in arid and semiarid Australia, palaeovalleys often host productive aquifers (Magee, 2009). We define palaeovalleys as ancient valleys incised into underlying weathered bedrock, filled with alluvial or colluvial sediments (Scott and Pain, 2009; Graham et al., 2010; Clarke, 2014; Wilford et al., 2016). The sediments and the underlying weathered and fractured bedrock can form aquifers. A challenge in palaeovalley aquifer system conceptualization is defining the interaction between the palaeovalley sediments and the fractured bedrock aquifer. Understanding this connectivity is important for sustainable management, notably in understanding replenishment or recharge and for determining potential impact from pumping.

Palaeovalleys are well studied and exist worldwide (Blum et al., 2013). They have been mapped using digital elevation models and spaceborne data from the National Oceanic and Atmospheric Administration — Advanced Very High-Resolution Radiometer (Tapley, 1988), imaged by seismic refraction (Leslie et al., 2000; Deen and Gohl, 2002; Goodwin et al., 2017) and airborne electromagnetics (AEMs) (Wilford, 2009; Gilfedder and Munday, 2013; Ley-Cooper and Munday, 2013; Wilford and Thomas, 2013). Palaeovalleys have been drilled extensively, and their stratigraphy has been characterized (Hou et al., 2003a, 2003b, 2008; Hou, 2008; Alley et al., 2009). Palaeovalleys can host large sedimentary mineral deposits and act as critical groundwater resources in semiarid to arid climates (Anand and Paine, 2002; Hou et al., 2003b; Broekert and Sandiford, 2005).

At the base of a palaeovalley, an ancient weathering profile may separate a deeper fractured bedrock aquifer from the overlying palaeovalley sediment fill. A well-developed weathering profile has unaltered and unfractured bedrock at the bottom and chemically depleted materials on top (Taylor and Eggleton, 2001; Scott and Pain, 2009). A significant part of the weathering profile will be comprised of saprolite. We define saprolite as chemically weathered bedrock that retained some expression of the original fabric of the bedrock (Scott and Pain, 2009; Graham et al., 2010; Wilford et al., 2016). The thickness, spatial distribution, and hydrogeologic properties of this ancient weathering profile can influence the hydrogeologic interaction between groundwater in palaeovalley and fractured bedrock aquifer systems.

This study combines AEMs, seismic refraction and reflection, and borehole nuclear magnetic resonance (NMR) to quantify the thickness and spatial distribution of the weathering profile as well as its hydrogeologic properties relative to the overlying palaeovalley sediments. We posit that the combination of AEM, seismic, and NMR data can, in some geologic settings, robustly define hydrogeologic units in which boundaries are controlled either by a porosity change (e.g., due to lithology) or by changes in water quality (e.g., solute content). There is no guarantee that the elastic and electrical boundaries will occur at the same spatial location. The robustness of this approach stems from the fact that all of the three measurements are sensitive to porosity at varying magnitudes. The seismic and NMR responses are dominated by the effects of porosity, whereas an electromagnetic (EM) response will be mainly influenced by the electrical conductivity of the fluid contained in the material being sensed.

The multimethod geophysical approach used in this study revealed that the main electrical boundaries, identified by AEM, are different than the elastic boundaries identified through seismic refraction and reflection. It is inferred that the contrast in electrical conductivity identifies the top of saprolite because this is driven mainly by a change in the groundwater total dissolved solids (or solute content). In contrast, the seismic data identified the top of the fractured bedrock but failed to identify the top of the saprolite. We attribute this to the porosity in the saprolite being similar to that of the overlying palaeovalley sediments. This observation is supported by water contents in the borehole NMR. These observations suggest that, although porous, the saprolite layer dividing the palaeovalley and the fractured bedrock limits the connectivity between the two aquifers. This critical information addresses one of the key conceptualization challenges in such palaeovalley systems, namely, that the interaction between the palaeovalley sediments and the underlying fractured bedrock aquifers may be limited.

The site is located in the Musgrave Province in central Australia on the Anangu Pitjantjatjara Yankunytjatjara (APY) lands 5 km southeast of the settlement of Kaltjiti (Figure 1). The site is part of the Lindsey-East palaeovalley system (Hou et al., 2003b, 2007; Hou, 2008; Krapf et al., 2019). This region is remote, which make collecting hydrogeologic observations and drilling logistically challenging and the available data sparse. The existence of palaeovalleys in this region has been known previously (Hou et al., 2007; Magee, 2009), but a regional AEM survey (Ley-Cooper, 2017a, 2017b) mapped the geometry of the palaeovalley network in detail (Gilfedder and Munday, 2013; Ley-Cooper and Munday, 2013; Munday et al., 2013) (Figure 1). The palaeovalleys in this region are critical to understand regional groundwater resources (Magee, 2009; Munday et al., 2020).

The palaeovalleys are cut into a Proterozoic weathered crystalline bedrock comprised of a mix of the gneissic rocks of the Birksgate complex (1605 ± 6 Ma) and the granitic plutons of the Pitjantjatjara supersuite (1219 ± 12 Ma) intermixed with dikes of Alcurra Dolerite (1085 ± 2 Ma) (Geoscience Australia and Australian Stratigraphy Commission, 2017). The basement rocks underwent at least one phase of intensive deep weathering and erosion from the Late Cretaceous through the Cenozoic, producing 100 m plus deep weathering profiles (Conor, 1987; Krapf et al., 2019). In the Late Paleogene to Early Neogene, the network of palaeovalleys was carved into the weathered basement rocks (Krapf et al., 2019). In the Mid to Late Neogene, a shift to a warmer, wetter, subtropical climate caused the region to transition from an eroding landscape to an aggrading landscape, which filled the palaeovalleys with locally sourced alluvial and fluvial sediments from the Musgrave Ranges (Magee, 2009; Krapf et al., 2019; Munday et al., 2020). The contemporary landscape is characterized by a topographically flat surface contacting the base of the Musgrave Ranges that define the northernmost region of the APY lands (Pawley et al., 2014; Krapf et al., 2019) (Figure 1). Current climate for the study area is semiarid to arid with a hot, dry desert climate, short cool to cold winters and low unreliable rainfall (Varma, 2012; Munday et al., 2020).

A compilation of drilling, geochemistry, and geophysics was used to construct a hydrogeologic conceptual model of the APY region (Munday et al., 2020). The region is typically divided into two aquifer units: alluvial and colluvial sediments, and fractured bedrock. Recharge occurs in the Musgrave Ranges in the north and regional groundwater flows southward. There are no perennial streams, and runoff occurs rapidly in response to large and local sporadic rain events. Chloride mass balance measurements suggest episodic and localized groundwater recharge into the fractured bedrock aquifers that outcrop in the Musgrave Ranges, on average between 2 and 20 mm/year (Munday et al., 2020). Furthermore, episodic recharge is possible in localized ephemeral drainages traversing surficial alluvium and colluvium though it is smaller occurring at an estimated 0.5–10 mm/year (Munday et al., 2020). The low hydraulic gradient between the northernmost region and southernmost regions combined with the large distance between monitoring wells means that the local groundwater flow paths are difficult to observe.

In the current hydrogeologic conceptualization a key question remains: Are the fractured bedrock and palaeovalley aquifers hydrologically connected (Figure 2)? Regional magnetics, groundwater levels and salinity, and bore yields suggest that aquifers are compartmentalized. It is unlikely that groundwater can flow from the recharge region in the north and exit the region solely through the fractured bedrock aquifer (Munday et al., 2020); rather the water most likely discharges into the overlying sediments. The existing groundwater monitoring network has insufficient spatial coverage and resolution to characterize the level of hydraulic connection between both aquifers.

To understand the hydraulic connectivity between the palaeovalley and the fractured bedrock aquifers, we focus on mapping the spatial distribution and relative hydrogeologic properties of the weathering profile that divides them (Figure 2). At the top of a well-developed weathering profile is the mottled zone. Below the mottled zone is saprolite. Saprolite (Scott and Pain, 2009; Graham et al., 2010; Wilford et al., 2016) formed from crystalline rocks can be porous, upward of 30% (Emerson et al., 2000; Goodfellow et al., 2014; Holbrook et al., 2014; Flinchum et al., 2018a). Under the saprolite is saprock, which is a rock that is chemically altered but more physically related to the underlying fractured bedrock. The saprock transitions into fractured bedrock that could be locally weathered along fractures. The fracture density is assumed to be higher near the saprock boundary and gradates into unaltered and unfractured bedrock at depth.

For the purposes of this paper, we simplify the definition of the weathered/unweathered profile as comprising two layers based on their geophysical response (Figure 2). The first is the “saprolite layer” that includes the mottled zone and the saprock, in addition to the saprolite. It is overlain by transported cover (palaeovalley sediment fill). The saprolite is interpreted to be more resistive than the overlying palaeovalley sediments but has a similar seismic velocity. The second layer is the fractured bedrock layer. We do not define the top of unfractured bedrock because it is a gradational boundary. The top of fractured bedrock is defined by an increase in seismic velocity. The fractured bedrock has electrical conductivity similar to the saprolite. Although we focus on mapping the thickness of saprolite under a palaeovalley, the presence of saprolite is known to exist (from drilling) throughout the region and, where present, it divides the colluvial and alluvial sediment fill and fractured bedrock (Figure 2).


The South Australian Government funded two AEM surveys across the Musgrave Province to help manage South Australia’s groundwater resources (Ley-Cooper, 2017a, 2017b). We focused on the region where the AEM data were collected with a SkyTEM dual-moment 275/25 Hz time-domain EM system. The AEM data were inverted using a 1D layered earth inversion using AarhusInv inversion codes (Auken et al., 2014). Specifically, we used a smooth model laterally constrained inversion (LCI) (Auken and Christiansen, 2004; Auken et al., 2005). Each sounding was defined by a smoothly varying 30-layer 1D model with fixed vertical discretization. The topmost layer thickness was set to 2 m, and the depth to the bottom layer was 600 m. The thickness of each layer increased logarithmically with depth.

Approximately 17,000 line km of AEM data were collected in the regional survey with line spacing at approximately 2 km (Munday et al., 2013; Davis et al., 2020). Within this survey domain, a few select regions were flown with higher resolution line spacing of 250 m. The higher resolution data sets were acquired in areas of particular focus for groundwater resource investigations, and the seismic data and drillholes discussed here were acquired in one of these, specifically to the south and east of the town of Kaltjiti (Figure 1). Data for this subregion were reinverted with different regularization parameters, including layered or blocky model (four layers), a smooth 30-layer spatially constrained inversion, and a sharp inversion (Vignoli et al., 2012), which uses a focused regularization technique allowing for an accurate reconstruction of resistivity distributions while maintaining the capability to reproduce horizontal boundaries (details in Appendix  A). The use of different inversion approaches was undertaken to assess whether a similar conductivity structure was defined for the area of interest — the Lindsay East Palaeovalley located in the east of Kaltjiti township (Davis et al., 2020). These reinversions show a similar conductivity structure for a transect across the Lindsay East Palaeovalley (Davis et al., 2020). A separate and independent ground-based investigation validated the observed conductivity structure defined in the inverted AEM data (Davis et al., 2020).


A series of boreholes were drilled in August 2018 and are located approximately 20 m to the southwest of the main road through the study site (Figure 1). The first hole, DH-1a, is one of the deepest holes drilled for groundwater in the Musgrave Province and extends to a depth of 115 m, with core collected to a depth of 90 m (Keppel et al., 2018; Costar et al., 2019). This hole and a series of others (screened at different depths) targeted the Lindsay East palaeovalley sediment-fill sequence. This configuration allowed the groundwater to be sampled at different levels within the palaeovalley sediment package (Costar et al., 2019). A second hole, DH-1f, was drilled to a depth of 78 m. DH-1a and DH-1f were expected to encounter bedrock; however, competent bedrock was not encountered in either.

Seismic data collection

The seismic data were collected with a 48-channel Geometrics system with 10 Hz vertical geophones spaced at 6 m. The source was a truck-mounted 40 kg accelerated weight drop that struck a 20 × 20 × 2 cm steel plate. At each shot location, we stacked the records 10 times. The shot spacing was 6 m, or every geophone. In general, we attempted to include 24 offline shots before and after the last geophone (extending 138 m off each end), although time constraints made these offline shots variable. The 1.7-line km of data were collected in seven different rolls. A gap exists between roll 7 and roll 1 because of a narrowing of the road raised safety concerns (Figure 1). An additional 300 m of profile was collected on the easternmost portion of the palaeovalley. All of the distances on this additional profile are negative with respect to the original start of the profile. The seismic data were processed twice: once for seismic refraction tomography and then to generate a reflection image.

Seismic refraction

Seismic refraction tomography uses the first-arrival times of the seismic energy to estimate seismic velocity of the subsurface. First-arrival times were manually picked. Noisy traces were not picked. Because the acquisition geometry (i.e., the shot spacing and geophone spacing) was optimized for reflection, we had 7700 reciprocal measurements. The data quality was high, as indicated by 87% of the reciprocal traveltimes being less than 2 ms. All of the arrivals were assigned a data weight prior to inversion based on the distance from the shotpoint using their distance from the source (details in Appendix  B). The inversion was run in Python Geophysical Inversion Modeling Library (PyGIMLi) (Rücker et al., 2017). Given the large volume of traveltime picks (22,057), we ran a bootstrapping algorithm to quantify uncertainties in the velocity (details in Appendix  B). The final velocity model is the average of 20 inverted models after randomly removing 40% of the traveltime picks. The model fits are determined by the χ2 value — which is a weighted-type average that incorporates our uncertainty of each data point or by an rms error between the observed and modeled traveltimes. Additional details about model fits and specific inversion parameters can be found in Appendix  B.

Seismic reflection

We applied a traditional reflection processing flow to build the reflection image. All of the data were processed using Seismic Unix (Stockwell, 1999). Even with efforts to increase the signal-to-noise ratio, we elected to stack the data with gains applied prior to normal moveout (NMO) to bring out the structure, working under the assumption that we could not stack what we could not see. We elected to stack the data using a trace-normalized gain, but a comparison between different gains prior to stacking can be found in Appendix  C.

To process the reflection data, we followed a seven-step workflow. First, the geometry was verified and noisy traces were identified and removed; second, all traces were passed through a band-pass filter between 50 and 200 Hz; third, the shot gathers were sorted into super common-midpoint (CMP) gathers comprised of three CMP gathers (Figure 3); fourth, surface waves were removed using a surgical mute as a function of offset; fifth, an estimate of a single NMO velocity for each supergather in the data set was determined (Figure 3); and sixth, the NMO velocity function was smoothed out (Figure 3) using a low-pass filter and then we flattened and stacked the data using the smoothed NMO profile (Figure 3a). Finally, we applied a time-to-depth conversion at each stacked trace using the smoothed NMO velocity (Figure 3a). Although we believe the reflection image could be cleaned up further, for the application of defining the bottom of a palaeovalley, we believe that this basic processing flow is sufficient. Details pertaining to each processing step can be found in Appendix  C.

Nuclear magnetic resonance

NMR in a noninvasive geophysical technique that provides a way to quantify aquifer properties. The NMR measurement transmits an EM pulse that produces an oscillating magnetic field. The oscillating field perturbs an existing magnetic moment that is caused by the alignment of nuclear magnetic spin moments in hydrogen protons when they are immersed in a magnetic field. After the oscillating field is turned off, the magnetic moments produced by the hydrogen protons precess around the background magnetic field axis because they relax back to equilibrium (Bloch, 1946; Bloch et al., 1946; Torrey, 1956; Kleinberg, 1999).

During relaxation, the bulk magnetic moment generated by the protons produces a quantifiable time-varying and exponential decaying magnetic field. The decay rate of the signal is represented by a sum of exponential decays each quantified by a relaxation time (T2 in seconds or milliseconds). This sum of exponential decays is often expressed as a distribution in which the values in this distribution represent water relaxing from different pore sizes in a given geologic material (Brownstein and Tarr, 1979; Keating and Falzone, 2013; Falzone and Keating, 2016a). The T2 relaxation time has been linked to pore size (Kleinberg and Horsfield, 1990; Mohnke and Yaramanci, 2008; Grunewald and Knight, 2011; Falzone and Keating, 2016b) and hydraulic conductivity (Straley et al., 1987; Kenyon et al., 1988; Howard et al., 1993; Walsh, 2008; Weller et al., 2010; Dlubac et al., 2013; Ren et al., 2019). Theoretically, the maximum amplitude of the NMR signal is proportional to the number of protons excited, which, in hydrogeologic applications, is from groundwater. In reality, instrument dead times prevent us from measuring the maximum amplitude, and we rely on the multiexponential fit to estimate the amplitude of the signal. As a result, measurement noise can influence the water content estimates. Under fully saturated conditions, NMR can provide an estimate of porosity.

Boreholes DH-1a and DH-1f (Figure 1) were logged with the Vista-Clara Inc. Javelin JP350 borehole NMR tool (Walsh et al., 2010). Diameters defining the shells of investigation for the JD350 are 27, 30, 34, and 38 cm. The Javelin probe uses two recovery times: The first is optimized to characterize shorter relaxation times (Tr1 = 80 ms), and the second is optimized to characterize longer relaxation times (Tr2 = 200 ms). These values represent the time between each stack. For the shorter recovery time (Tr1), the echo time was 1.6 μs, and for the longer recover time (Tr2), the echo time was 4 μs. We stacked the shorter recovery time (Tr1) 24 times and the longer recover time (Tr2) 4 times. These parameters were selected to optimize the collection speed so that the full extent of all existing wells could be completed in the allocated field time. A radio-frequency noise cancellation device was used to eliminate near-surface noise, especially in the readings near the surface. Data processing, stacking, and filtering were conducted using software developed by Vista-Clara Inc. We express water as either being bound by clay (T2 < 3 ms), bound (3 ms < mean T2 < 33 ms), or mobile (T2 > 33 ms) (Timur, 1969; Walsh et al., 2013).


The conductivity structure exhibited by inverted AEM line profiles in the Musgrave Province tends to define or approximate a three-layer structure (Figure 4a). The first layer is generally resistive (conductivities <40 mS/m), and it is associated with aeolian sands and unsaturated sediments above the water table (Figure 4a). Below this layer is an electrically conductive layer (conductivities >40 mS/m), which varies in thickness. There is significant heterogeneity within this layer. In places, this layer exhibits zones with conductivities of >200 mS/m (Figure 4a). Commonly, this second layer is associated with the alluvial and colluvial sediment fill of the palaeovalleys and effectively defines the palaeovalley geometry (Figure 4a). Beneath the electrically conductive layer, there is a marked decrease in electrical conductivity — a resistive layer with conductivities of <40 mS/m. This boundary, dividing the second and third layers, occurs at between 40 and 50 mS/m.

To generate an approximation of the palaeovalley geometry, we extracted the depth to the 45 mS/m contour along all the flight lines in the study region (Figure 4a). The flight spacing in the area was 250 m. We interpolated the results to a 100 × 100 m grid using a continuous curvature with a spline-tension algorithm (Smith and Wessel, 1990). This interpolated surface was then deemed to represent the top of saprolite in the region mapped (Figure 4b). The main trunk of the Lindsay East palaeovalley runs north–south, and it is approximately 2 km wide. At its deepest point, the bottom of the palaeovalley is greater than 100 m (Figure 4c). The depth to the saprolite reveals two smaller tributary palaeovalleys that run approximately perpendicular to the main north–south-trending palaeovalley (Figure 5). Outside of the palaeovalleys, the ancient topography is characterized by low-relief undulating hills (Figure 4), consistent with a landscape of low relief.

We used the AEM surface to refine the location for the seismic profile and to select the locations to drill boreholes DH-1a and DH-1f (Figure 4c). Logistic and cultural constraints meant the seismic data had to be collected along main roads (Figure 1). We placed the 1.7-line km seismic profile along a road so that it would image the deepest part of the palaeovalley and the western edge of the valley. The borehole DH-1a was placed into the deepest part of the palaeovalley, and DH-1f was selected to intersect the Lindsay East palaeovalley at the middle of the hillslope that defines its western edge.

Drilling results

Drilling logs, based on cuttings, show distinct units (Figure 5). The key lithologic change is a 25 m thick clay sequence (from approximately 65 to 90 m). Palynological samples from this formation revealed that it was deposited in a subaqueous environment characterized by fluctuating freshwater to brackish conditions in a fluvial to lacustrine depositional environment (Krapf et al., 2019). This clay sequence is visible on the AEM profiles as a thick conductor that lies at the bottom of the palaeovalleys but truncates near the edges (Krapf et al., 2019) (Figure 4a). All of the other sediments are interpreted to be sourced from the Musgrave Ranges and have a similar mineralogic composition because the underlying bedrock and are believed to be mobilized weathered materials form the bedrock in the north.

In the palaeovalley fill (Figure 5), DH-1a produced 18 L/s (high for the region) and had low salinity (approximately 800 mg/L) (Costar et al., 2019). DH-1a sampled just above the saprolite boundary, and it had low yield (<2 L/s) with moderate quality (1200 mg/L). DH-1a encountered saprolite at approximately 97 m (Figure 5a), but it never penetrated into fractured bedrock. Groundwater was never sampled or pumped in the saprolite at DH-1a. DH-1f encountered mottled saprolite at a depth of 43 m. The drillers expected to find bedrock near 45 m based on the AEM data; they continued to 78 m with no indication of intersecting competent bedrock. Within the saprolite layer, DH-1f had low water yields (<1 L/s) and slightly lower salinity values (1000 mg/L) than water sampled in the palaeovalley sediments above the saprolite at DH-1a (Costar et al., 2019).

Seismic refraction

The final velocity model produces traveltimes that match observed traveltimes well with a χ2 value of 1.79 and an rms value of 2.79 ms (Figure 6). More details about the model fits can be found in Appendix  B. The refraction data are well constrained from 500 m onward. The transparent overlay represents areas where we have no ray coverage (Figure 5). Ray coverage through the model is good with the exception of the first 500 m of the line and the existence of the large data gap because of the curve in the road. Deeper parts of the model are not well constrained between 0 and 500 m because we were unable to pick far-offset data on roll one because of noise generated by drillers living on site.

There is a strong and well-constrained refraction indicated by an increase in velocity (Figure 6). We have highlighted the boundary using the 4500 ± 500 m/s (Figure 6) based on the velocity of outcropping granite (Holbrook et al., 2014; Flinchum et al., 2018b) and laboratory measurements on bedrock from a palaeovalley in the Lawler’s region in Western Australia (Emerson et al., 2000). Near this boundary, the vertical velocity gradient is the strongest, indicating a strong change in elastic properties, specifically seismic velocity. This boundary occurs 45 m below the final depth of DH-1f and 100 m below the depth where DH-1f encountered saprolite (Figure 6). The boundary becomes less constrained under DH-1a, but it is still deeper, or not present, at the location where DH-1a intersected saprolite.

Finally, the tomography model shows homogeneous velocities from the surface to the depth of the main refractor (approximately 4500 m/s). The velocities average just more than 2000 m/s without much lateral heterogeneity. The exceptions are small (20–40 m wide) high-velocity zones (>3000 m/s) just below the surface.

Seismic reflection

The final depth-converted seismic reflection section reveals a great amount of subsurface structure (Figure 7). A strong and continuous reflection captures the palaeovalley geometry — specifically the palaeovalley edges (Figure 7). Between −500 and −430 m, the strong reflector decreases in elevation by 60 m. On the other side of the palaeovalley, between 1250 and 1380 m, the elevation changes by 80 m. This strong continuous reflection also defines a lot of lateral heterogeneity throughout the bottom of the palaeovalley. Near DH-1a, the elevation of this reflector is 400 m (a depth of 120 m). Approximately 500 m along the profile, the elevation of this reflector is 435 m. At the location of DH-1f, the elevation of this reflector is 390 m, which is 10 m deeper than at DH-1a.

The reflector is continuous across the entire profile with the exception of the region between 0 and 250 m and a few breaks at approximately 520, 635, and 1350 m (Figure 7). The reflector is vertically homogeneous in the time domain, suggesting that our velocity depth conversion is influencing the structure. The reflector has its deepest point at an elevation of 360 m at 1250 m along the profile and comes to an elevation of 450 m on either end of the profile (Figure 7a).

There is a region between 0 and 250 m, where the strong reflection weakens (Figure 7). In this region, there are also many reflections below the major boundary. These reflections extend from approximately 400 m in elevation all the way down to 250 m in elevation. There are a few other regions between 1000 and 1250 m that show strong reflectors below the key boundary. The reflectors occur well below the final depths of DH-1a. These reflectors below the main reflector are strong and observable in the individual supergathers.

There are only a few locations where we observed reflections above the key boundary. The most prominent set of these reflectors starts at approximately 750 m and extends to just under 200–935 m. There is another set of strong reflectors between −300 and −250 m and again between 400 and 500 m. These are the only prominent reflections that we see above the bedrock boundary.

Nuclear magnetic resonance

Between DH-1a and DH-1f, there is more than 140 m of downhole NMR measurements (Figure 8). The standing water table in both boreholes is at 8 m. DH-1a and DH-1f have similar water contents and T2 distributions in the palaeovalley sediments (Figure 8). The average water content in the palaeovalley sediments at DH-1a is 0.18 ± 0.06 m3/m3, and the mean T2 value is 38.0 ± 41.7 ms. The large standard deviation on this measurement is caused by the multimodal nature of the T2 distributions in the palaeovalley sediments (Figure 8). The average water content in the palaeovalley sediments at DH-1f is 0.19 ± 0.05 m3/m3, and the mean T2 value is 22.8 ± 16.8 ms.

There are boundaries that can be identified in the borehole NMR logs. In DH-1a, at a depth 65 m, there is a clear change in water content and T2 distributions (Figure 8). The water content increases, the T2 distributions shorten, and the multimodality of the distribution transforms into a uniform distribution. This change corresponds to the silty clay sequence in the drilling logs (Keppel et al., 2018; Costar et al., 2019; Krapf et al., 2019) (Figure 5). In DH-1f, at a depth of 45 m, the water content increases slightly and the T2 values transition from multimodal to a more uniform distribution. This transition corresponds to the location where drilling logs mark the transition to mottled saprolite (Keppel et al., 2018; Costar et al., 2019; Munday et al., 2020) (Figure 5). Within the saprolite, there are two regions based on their water content and T2 values. The first is a region between 45 and 55 m represented by slightly higher water content (0.30 ± 0.05 m3/m3) and higher mean T2 values (20.0 ± 6.0 ms). The second is a region between 55 and 67 m that has slightly lower water content (0.23 ± 0.03 m3/m3) and slightly lower mean T2 values (17.4 ± 6.8 ms) (Figure 8).

In this section, we explore the data in three sections. First, we show that the elastic and electric boundaries occur at different depths and estimate uncertainty on the depth to the seismic reflection boundary. Second, we explore why the discrepancy between the electrical and elastic data exists. Finally, we use the discrepancy to estimate saprolite thickness across the palaeovalley and show that there is limited hydraulic connection between the fractured bedrock aquifer and the overlying palaeovalley sediments.

Time-to-depth uncertainty

With reflection imaging there is uncertainty associated with the time-to-depth conversion. An important question needs to be addressed prior to combining different geophysical and hydrologic observations: How reliable is our single-layer NMO velocity model when used to go from time to depth (Figure 3b)? The NMO velocities derived from the CMP supergathers created additional structure in the time-to-depth conversion (Figure 7). The two boreholes along the seismic profile never intersect the strong reflection boundary to validate the depth of the strong reflection. Instead, we must rely on the velocities from seismic refraction tomography (Figure 6).

The seismic refraction image was processed independently using only the arrival times of the data. During the refraction inversion, we gave the model adequate space to ensure that the rays turned on their own. We also selected inversion parameters that gave the inversion the ability fit sharp boundaries. These two criteria, combined with the simple gradational starting model, suggest that the location and geometry of the main refractor is data driven (Figure 6). We can compare the seismic refraction results with the reflection results by overlying the 4500 ± 500 m contours on the seismic reflection image. The match between the refraction boundary and the strong reflection is best in locations where we have the best ray coverage. These locations also have a high number of far-offset traveltime picks. A good example of this is between 1000 and 1500 m.

The general structural trends between the refraction and reflection data sets suggest that our NMO velocities are appropriate for converting from time to depth. The 4500 m/s velocity contour and the strong reflection intersect at an elevation of 380 m at 1200 m along the profile and again at an elevation of 435 m at 1400 m along Figure 7 suggesting that the rise in bedrock topography at the palaeovalley edge is consistent between the two methods. Because refraction tomography produces smooth velocity models (Zelt et al., 2013; Rawlinson and Spakman, 2016), the reflection image shows significantly more structure.

In addition to the similar structural trends in the inverted velocity model and the reflection image, the NMO velocities are consistent with the inverted refraction velocities. The NMO velocity represents the average velocity of all layers above a picked reflection. We expect the NMO velocities to be an average of velocities in the traveltime tomography model above bedrock (i.e., above 4500 m/s). This is visually consistent with the traveltime tomography velocity model (Figure 6). It is confirmed when the histograms of the NMO velocity function and all pixels in the velocity model are plotted on top of each other. The average NMO velocity from our picks is 2130 m/s with a standard deviation of 230 m/s (Figure 3b). A change in velocity of 230 m/s would put depth uncertainty on a reflector at ±11.5 m. Thus, a conservative estimate of error on the depth surface of the reflector is approximately 15 m. The most accurate depth estimate will be at the key reflector. Changes in velocity above and below this boundary are not accounted for in our simple velocity depth conversion.

The geophysical discrepancy

Our geophysical data set provides us with two geophysical parameters: (1) electrical conductivity from the AEM surveys and (2) seismic impedance, which is the ratio between the material density and the seismic velocity above and below an elastic boundary. We observe a discrepancy between these two geophysical properties across the palaeovalley system (Figure 9b). We are confident in the location of the strong reflector because the independently inverted traveltime tomography models show a strong change in velocity in same location. We are also confident in the location of the decrease in electrical conductivity in the AEM data because they has been inverted multiple ways, with consistent results being generated, and the derived conductivity structure is directly comparable with that independently derived using borehole inductive conductivity data and a ground-based TEM method that yields similar structure (Davis et al., 2020). Thus, it is unlikely that the locations of the boundaries are caused by errors or smoothing — they are data driven. The discrepancy is real; therefore, it contains information about the hydrogeologic structure.

When the interpreted saprolite surface, calculated by extracting the 45 mS/m contour in individual AEM transects (surface from Figure 4b), is overlaid on the seismic reflection profile, and it is always shallower (Figure 9b). The saprolite surface, interpreted by AEM, is well constrained across the seismic profile containing the intersection of nine north–south flight lines. The interpreted saprolite surface likely contains additional heterogeneity that was smoothed out by the LCI. At the palaeovalley center (at DH-1a), there is a conductive layer that is approximately 20 m thick. This conductor is not present near the valley edges, and it is associated with a clay-rich estuarine sedimentary sequence linked to the valley development and fill (Keppel et al., 2018; Krapf et al., 2019) (Figure 5).

The magnitude of the difference between the AEM boundary and the elastic boundary can be significant (Figure 9). It varies across the palaeovalley, with the largest differences occurring on the palaeovalley edges. Near DH-1f, the difference between these two boundaries is greater than 80 m. On the easternmost side of the palaeovalley (−500 m), the difference between the two boundaries is also significant, differing by 60 m. Furthermore, there is a lack of reflections between these two boundaries (Figure 9). These differences cannot be resolved by regularization or changing the inversion parameters in the AEM or by adjusting the NMO velocity model in the reflection data.

In general, the saprolite surface, as interpreted by AEM (Figure 4), appears to be associated with the upper saprolite when it is encountered in DH-1a and DH-1f (Figure 9). At DH-1a, the difference is less than 2 m, and, at DH-1f, the difference is 15 m. These differences may, in part, be explained by the footprint of the AEM system and the lateral averaging that occurs in the stitching between soundings. Based on the drilling results, we suggest that the strong electrical boundary defined in the AEM data and represented by a conductor above a resistor effectively determines the top of the saprolite. Although there is no evidence of a transition to fractured bedrock in the AEM data, this may occur in places. Typically, saprolite is often thought to be more electrically conductive because of the presence of clay, and/or the soluble salts it may contain within its pore structure. The lack of fractured bedrock at either DH-1a or DH-1f indicates that the saprolite is more resistive than the palaeovalley sediments above it. Thus, the material above the saprolite interface (the blue line in Figure 9) is more electrically conductive than the material below.

We do not expect a large mineralogical difference between the alluvial and colluvial sediments and the underlying saprolite because former were largely derived from the regolith itself. The biggest difference is that the palaeovalley sediments have been transported whereas the saprolite is weathered in situ (Hou et al., 2003b, 2008). Drill cuttings confirm that the alluvial/colluvial fill at the base of DH-1 is of similar character to the underlying saprolite, namely, they comprise kaolinitic clays and quartz sands with slightly rounded grains.

Without a strong mineralogical change, the elastic properties, specifically the seismic velocity, will be dominated by porosity. An increase in porosity will decrease the seismic velocity (Hashin and Shtrikman, 1963; Budiansky and Oconnell, 1976; Nur et al., 1998; Bachrach et al., 2000; Berryman et al., 2002). However, the rate of decrease depends on the rock type and the pore structure, as well as the rock-physics model (Mavko et al., 2009). Seismic velocity is also affected by saturation (Nur and Simmons, 1969; Gregory, 1976; Bachrach and Nur, 1998; Bachrach et al., 2000; Berryman et al., 2002), but it will not be significantly affected by the water quality, or more specifically the solute content of the water within the pore space. The water table at the study site was at 8 m, well above the transition from palaeovalley sediments to saprolite. Thus, the material above the strong reflector (the black line in Figure 9) and refractor (the gray line in Figure 8) is probably softer and more porous than the material below.

The lack of an elastic boundary at the top of the saprolite suggests that the porosity of the saprolite is close to, if not equal to, the porosity of the alluvial and colluvial sediments in the bottom of the palaeovalley. This is supported by the downhole NMR data at DH-1f (Figure 8). Here, we assume that the measured water content is a good proxy for porosity by assuming that all of the pores are filled with water and that they are within the detectable range of the NMR instrument. The NMR log at DH-1f shows that water content increased across the saprolite boundary (Figure 8c). The transition from palaeovalley sediments into saprolite was associated with shift from a bimodal distribution to a shorter and more multimodal distribution (Figure 8). Furthermore, the water content measured in the saprolite at DH-1f is consistent with measured porosity of saprolite developed on crystalline rocks between 0.25 and 0.50 m3/m3 (Johnston, 1987a, 1987b; Emerson et al., 2000; Lane and Pracilio, 2000; Munday et al., 2001; Rutherford et al., 2001; Holbrook et al., 2014, 2019; González-Álvarez et al., 2016; Flinchum et al., 2018a, 2019; Hayes et al., 2019).

The presence of a change in electrical properties at the saprolite interface is likely caused by two variables because the elastic data rule out changes in porosity and the drilling logs limit changes attributable to lithology. The first variable responsible for the electrical contrast is the groundwater electrical conductivity. If groundwater conductivity in the saprolite is less than the overlying palaeovalley sediments, we expect the overall conductivity in the former to drop. This change would not affect the NMR data or seismic data. The second variable is more complicated but is related to a change in pore structure, specifically an increase in tortuosity. An increase in tortuosity might produce a drop of electrical conductivity and a drop in hydraulic conductivity. With electrical conductivity measurements alone, even if the porosity, lithology, and saturation are known, it is difficult to determine which variable is causing the change. It is likely a combination of both variables.

Water sampled from the saprolite at DH-1f had low yields (<1 L/s) and salinities of approximately 1000 mg/L (Costar et al., 2019). Water sampled in the palaeovalley sediments above saprolite at DH-1a had low, but slightly higher yields (<2 L/s) and a slightly higher salinity of approximately 1200 mg/L (Costar et al., 2019). From this observation at least a small part of the decrease in conductivity in the saprolite is likely due to a decrease in fluid conductivity or contained soluble salt content of the water. Based on the total dissolved solids (TDS) and electrical conductivities reported by Keppel et al. (2018), a 200 mg/L difference in TDS would produce a decrease in electrical conductivity of approximately 50 mS/m. Unfortunately, without electrical conductivity measurements of water and rock samples and more thorough testing, it is difficult to determine how much of the observed drop is due to a decrease in total dissolved solids within the saprolite.

The slightly fresher water in the saprolite at DH-1f was an unexpected occurrence. In saprolite, groundwater and its constituent solutes infiltrate and move through the saprolite by processes of diffusion (or matrix flow) and dispersion where movement is influenced by macropores and preferred pathways (Johnston, 1987a, 1987b). Slower rates of groundwater movement are associated with diffusion, which encourages higher groundwater residence times and the storage of salt. High electrical conductivities have been reported from AEM data sets for saprolite and have been shown to be associated with high water solute content (Lane and Pracilio, 2000; Munday et al., 2001; Rutherford et al., 2001; González-Álvarez et al., 2016).

The electrical conductivity is sensitive to the pore structure, specifically to tortuosity, which is often represented as an exponent in Archie’s law (Archie, 1942; Robinson et al., 2012). An increase in tortuosity would theoretically produce a decrease in hydraulic conductivity and electrical conductivity. In NMR studies, the T2 decay time is often associated with a drop in hydraulic conductivity and usually represents materials with smaller and tighter pores (Keating and Falzone, 2013; Falzone and Keating, 2016b, 2016a; Osterman et al., 2016). Thus, at DH-1f, the shift in mean T2 values from a multimodal distribution with a mean T2 average of 38.0 ± 41.7 ms in the palaeovalley sediments to a more Gaussian distribution with a mean T2 of 20.0 ± 6.0 ms suggests a change in the pore structure, especially that the saprolite is missing the larger pores (Figure 8). This slight decrease in mean T2 decay times suggests a decrease in hydraulic conductivity.

Estimating saprolite thickness

We calculate the saprolite thickness by subtracting the elevation of the AEM depth to the top of saprolite and the seismic reflection boundary depth to bedrock (Figure 9). Our data set suggest that the saprolite thickness underneath the palaeovalley is tens to hundreds of meters thick (Figure 9a). This is greater but on the same order of magnitude of the 70 m of weathered profiles observed elsewhere in the region (Krapf et al., 2019). The AEM boundary shows a clear change in electrical properties at or near the known top of saprolite. There is no indication of an elastic boundary near the top of saprolite in either the refraction data or the reflection data (Figure 9).

Our data show that saprolite thickness is laterally heterogeneous and varies from a maximum value of 120 m to less than 2.5 m across the 1.7 km profile (Figure 9a). In selecting regions where the palaeovalley and bedrock aquifer are in contact, we look at locations where the saprolite thickness is less than 15 m because of the uncertainty associated with time-to-depth conversion in the reflection processing. There are only two places where the saprolite thickness is less than 15 m (Figure 9a). The first is at the beginning of the profile at 0 and 50 m. The second is a region between 475 and 580 m. Between these two locations, only 8% of the fractured bedrock aquifer might be in contact with the palaeovalley aquifer suggesting that the palaeovalley and bedrock aquifers are most likely physically separated by a layer of saprolite.

A third, more speculative, possibility for structural control on the interaction between the two aquifers might be the presence of a buried tor or large corestones (Figure 9). There are few places where we see reflections above the strong boundary. One of these places is between 750 and 925 m (Figures 6d and 8). Given the prevalence of tors and corestones in the Musgrave Ranges approximately 60 km north (Figure 8c and 8d), the lack of continuity but the presence of clear reflections, and the proximity just above the main reflector, we interpret this set of reflections as a buried tor or weathered corestones. We do not have enough detail to determine if these corestones or tors are connected to the underlying fractured bedrock. Given the high localized recharge in the Musgraves around outcrops these features could provide additional vertical hydraulic conductivity into the fractured bedrock aquifer. If we count these features as additional surface area where the two aquifers might be connection as a conservative estimate, it still means that only 17% of the entire seismic profile has locations where the two aquifers are in contact. Thus, it still likely that the two aquifers are physically separated.

This combination of electrical, elastic, and NMR data allowed us to show that the saprolite is highly variable and can be locally thick (>100 m). It is rare, or unlikely, to find a location where fractured bedrock is in direct contact with the overlying alluvial and colluvial sediments. We propose that in the best-case scenario only 17% of the entire 1.7 km of data may be in physical contact. The NMR data show a loss of longer T2 decay times suggesting a drop of hydraulic conductivity. Thus, the saprolite, despite being porous, will likely act as an aquitard separating the palaeovalley and bedrock aquifers. More work is required to determine exact values of hydraulic conductivity, but it seems likely that the saprolite has a lower hydraulic conductivity than the palaeovalley sediments.

The geophysical observations provided two pieces of evidence to suggest that the saprolite limits the hydrogeologic connection between the valley fill and the underlying fractured bedrock system (Figure 9). First, the thickness varied from more than 120 m to less than 3 m. Yet only 8%–17% of the region represents locations where bedrock could be in contact with the overlying sediments. Potentially, the saprolite lines the entire palaeovalley (Figure 10). Thus, the presence of thick saprolite will significantly reduce the chances of direct hydrogeologic connection. Combined with the loss of longer T2 decay times in the saprolite shown by the NMR data suggesting a drop in hydraulic conductivity, we show that the saprolite limits the connection between the fractured bedrock and palaeovalley aquifers. Second, the presence of two geophysical boundaries (electrical and elastic) is valuable information. Although we cannot quantify the porosity and hydraulic conductivity through these data directly, they suggest that the porosity in the saprolite is similar to the bottom of the palaeovalley fills but the pore structure and perhaps the permeability decreases. A pore structure change associated with an increase in tortuosity would reduce the hydraulic conductivity, which is consistent with the lower yields observed in DH-1f in saprolite (Costar et al., 2019) and an increase in bound water identified via NMR logs (Figure 7). From the combined geophysical observations, we believe that the fractured rock and alluvial systems can be treated as two separate aquifers in the hydrogeologic conceptualization (Figure 10).

When using geophysical data to inform a hydrogeologic conceptualization, it is critical to understand the relationships between the measurable geophysical properties and the hydrogeologic properties of interest (e.g., porosity and hydraulic conductivity). Here, we explored the relationship between electrical conductivity and elastic properties and properties that influence groundwater flow and storage, mainly the porosity. We observed a significant discrepancy, greater than 80 m in some places, between boundaries identified by electrical and elastic properties. Thus, on a field scale, our data analyses demonstrated that a change in one geophysical property does not guarantee a change in the other. Although the boundaries in each geophysical method may not match, having a combined geophysical data set provided another layer of observations that were valuable in the integration of geophysical data into a hydrogeologic conceptualization.

We exploited a discrepancy between key geophysical boundaries to extract qualitative hydrogeologic information about the saprolite. The lack of an elastic boundary suggested that the saprolite has similar porosity to the overlying palaeovalley sediments; this was supported by an increase in water content in the NMR data. The decrease in electrical conductivity at the saprolite boundary suggested a change in pore structure to smaller, tighter pores, probably resulting in lower hydraulic conductivity. This was supported by the loss of longer T2 decay times in the NMR data. The combination of the data allowed us to map the distribution of saprolite in detail and show that that only a small percentage of the palaeovalley is in direct contact with the bedrock aquifer. When combined, these observations suggested that the saprolite acts as an aquitard limiting the connection between the two aquifers. This addressed a key challenge in the existing hydrogeologic conceptualization. Thus, the addition of multiple geophysical data sets has led to a better understanding of the hydrogeologic system, which will result in better management and exploration of new groundwater resources.

The authors would like to thank the Goyder Institute for Water Research (a research alliance between the South Australian government through the Department for Environment and Water, CSIRO, Flinders University, the University of Adelaide, and the University of South Australia) for its support of the research described here, which was conducted under the auspices of its G-FLOWS 3 Project. CSIRO’s Deep Earth Imaging Future Science Platform and CSIRO Land and Water also acknowledged for critical support for the work undertaken. Andrea Viezzoli is acknowledged for helping with the AEM inversions.

Data associated with this research are available and can be accessed via the following URL: https://doi.org/10.25919/5dcc9dedd8a67.



All of the airborne data presented in this paper were collected with a time-domain EM system in which the profiles run north–south. The region that we focused on (Figure A-1) was covered by the heliborne SkyTEM dual-moment 275/25 Hz EM and magnetic system (Figure A-1). All of the data were inverted using the AarhusInv inversion program (Auken et al., 2015) and the Aarhus Workbench LCI algorithm (Auken et al., 2002, 2005). Each sounding was defined by a 30-layer model in which the topmost layer thickness was set to 2 m and the depth to the bottom layer was 600 m. In between the thickness of each layer increased logarithmically with depth. Each sounding was inverted with a multilayer, smooth-model model in which the thickness of the layers was held constant.

For the inset study region (Figure 1), the focus of this study, survey line spacing was 250 m intervals. Within this region, some profiles have been reinverted with different parameters such as four-layer models, sharp inversions, and 30-layer smooth models (Davis et al., 2020) (Figure A-2). These inversion results all show similar structures. A separate investigation was also undertaken as a campaign to validate the airborne EM profiles with NanoTEM (ground TDEM) data (Davis et al., 2020). The NanoTEM data also show similar structures to the airborne data. For more information on the processing and collection of the airborne data, the reader is referred to Ley-Cooper and Munday (2013), Ley-Cooper (2017a, 2017b), and Munday et al. (2018).



First-arrival times were manually picked in a Python-based program (Flinchum, 2021). Prior to picking, all data were trace normalized to enhance the first-arrival time. To pick far-offset picks, we would occasionally apply a band-pass filter between 10 and 150 Hz to reduce some high-frequency noise probably caused by wind. Any traces in which there was not a clear arrival time form the seismic source were not picked. Overall, the data set contains 22,057 first-arrivals picks (Figure B-1). The dense shot geometry that is required for seismic reflection (but not refraction) meant that we had reciprocal traveltimes for almost all picks (except offline shots) (Figure B-1). The reciprocal traveltime is defined as the difference between the time required to travel from the source to the receiver when the source location and receiver locations are swapped. Reciprocal traveltimes should be zero and are used to indicate the quality of the picks. The reciprocals show an extremely high-quality data set (Figure B-1).

We tried using the weights from reciprocal traveltimes, but there was no discernable trend. Furthermore, we do not have reciprocal values for shots not located at a geophone (e.g., offline shots). To drive the error weight matrix in the PyGIMLi inversion, a weight must be given to every pick (Rücker et al., 2017). To ensure an error is assigned to every pick, we applied a linear function as a function of offset (distance from the source). In this case, the error (in seconds) is defined by
where offset is the distance from the source to the receiver, and the spread length was 282 m. This means that weight errors start with a 1 ms error at zero offset and increasing in a way such that the error at the farthest offsets (47*6 m = 282 m) was equal to 2 ms (equation B-1). We felt justified with the small data weights at far offsets because of the high quality of the reciprocal traveltimes (Figure 3b). These data weights are used to provide a stopping criterion for their inversion. Each error is read into the inversion and used to weight the traveltimes in which the smallest errors are given the most weight. Furthermore, the weights are used to calculate the χ2 fit:
where tobs is the observed traveltime, tmodel (ms) is the modeled traveltime (ms), σ is the assigned picking error (ms), and N is the total number of picks.

The inversion of the data was done in PyGIMLi. PyGIMLi runs based on a shortest path algorithm and a deterministic Gauss-Newton inversion scheme (Gance et al., 2012). The mesh extended at a depth of 180 m. This is well below the perceived depth of channel. We selected this depth so that velocities would not be artificially increased to get the rays to turn on their own. Based on the presence of a strong refractor in rolls 4, 5, and 6, we minimized the vertical smoothing parameter. This should allow the model to fit a sharper bedrock boundary, consistent with the observed picks. The inversion was set to run no more than 15 iterations or until the χ2 value dropped below one. There is a trade-off between fitting the near-surface arrivals and the far-offset arrivals. If the data weights are too big, the model tends to fit the longer offset and the near offsets can be off by 5–7 ms. This is an interesting observation that might have a rippling effect in traveltime tomography models.

Given the high number of first-arrival picks (22,057), we ran a bootstrapping algorithm to quantify uncertainty in the velocities. The bootstrapping algorithm removed 40% (8822 picks), calculated the errors (equation S1), and then ran an inversion. These inversion results were saved, and then the original 22,057 picks were used again to randomly removed another 40% and reinverted. This process was done 20 times. We kept track of each resulting velocity model and χ2 value, and at the end of 20 runs we averaged all velocities and can calculate a standard deviation (Figure B-2). The velocities do not vary all that much, and the velocity model is robust (Figure B-2).

We calculated the forward model through the average velocity model from the bootstrapping results (Figure 5). The forward model produces the possible raypaths through the final velocity model (Figure B-3). Ray coverage was very good with the exception of the first roll, in which we had significant noise and lacked far-offset picks (Figure B-3a). In the final model, areas with no raypaths are slightly transparent (Figure 5). Investigation of the modeled versus observed traveltimes shows that the model explains most of the observed data (Figure B-3b and B-3c). In the residuals (Figure B-3b), it is clear that the model is having some difficulty fitting the nearest offset data and some of the furthest offset data (Figure B-3b). Changing the data weights can help the user improve the fits of the near or far offsets of the data. The vertical gradient of the velocity model also shows a strong elastic boundary co-located with the strong reflector (Figure B-4).



We applied a basic reflection processing flow to build the reflection image presented in this manuscript. All of the data were processed using Seismic Unix. We acknowledge that with more resources and more time the image could be cleaned up. As a result, we caution the reader about over interpreting detailed structures. Even with all efforts to increase the signal to noise, we elected to stack the data with gains applied prior to stacking to bring out the structure, working under the assumption that we cannot stack what we cannot see. Thus, when interpreting the data, it is important to remember that these are not true amplitudes.

We followed the following basic steps to produce the final time-domain stacked image (Figure C-1). First, noisy traces were identified and removed. Second, data were passed through a band-pass filter between 50 and 200 Hz. Third, data were sorted into super CMP gathers. A supergather is the combination of three CMP gathers: the CMP at the location and the two adjacent CMPs. This increased the maximum fold from 48 to 144. This greatly increases the signal to noise at the cost of lateral resolution. The supergathers also made the semblance plots a lot stronger to pick reflections (Figure C-2). The super gathers will reduce the horizontal resolution by smearing, but after experimentation, the reflections were best observed within the supergathers.

Surface waves were removed using a surgical mute as a function of offset. The mute had a slope of 500 m/s (Figure C-3). We used a surgical mute instead of an f-k filter because the 6 m spacing resulted in severe spatial aliasing. This removal of surface waves means that even though we increased the fold the signal to noise is not great in at short offsets. The most reliable data occur between 0.8 and 0.3 s where there are a lot of traces to flatten and stack.

We estimated a single NMO velocity for all CMP gathers in the data set. We choose a single NMO velocity to keep things simple, especially in the depth conversion. Initially, we were interested in defining the bottom of the palaeovalley, which was presumed to be a strong reflector. The NMO velocity is an average velocity of everything on top of the reflector; thus, this velocity will be most accurate for the chosen reflector. During the picking of the NMO velocity, we used semblance plots as an aid. The semblance plot was calculated using Seismic Unix on automated gain controlled data with a window of 0.1 s. During the picking of the NMO velocity, we tried to pick the strongest visible reflector. Given that we manually picked all NMO velocity, there is coherent velocity structures between each supergather (Figure C-3). The average NMO velocities peak with the average velocities from the seismic refraction (Figure C-4). The depth of the strong reflector also varied as a function of time, but for the most part it changed gradually and as something that was geologically realistic.

We smoothed out the NMO velocity (Figure 3) using a low-pass filter and then flattened and stacked the data using NMO. At this stage is where we experimented with gains. We applied a simple time-to-depth conversion at each stacked trace using the smoothed NMO velocity. We also experimented with different gains prior to stacking; however the strong reflection shows up no matter how the data were gained prior to stacking (Figure C-5).

Biographies and photographs of the authors are not available.

Freely available online through the SEG open-access option.