Seismic P Wave Velocity Model From 3‐D Surface and Borehole Seismic Data at the Alpine Fault DFDP‐2 Drill Site (Whataroa, New Zealand)

The New Zealand Alpine Fault is a major plate boundary that is expected to be close to rupture, allowing a unique study of fault properties prior to a future earthquake. Here we present 3‐D seismic data from the DFDP‐2 drill site in Whataroa to constrain valley structures that were obscured in previous 2‐D seismic data. The new data consist of a 3‐D extended vertical seismic profiling (VSP) survey using three‐component and fiber optic receivers in the DFDP‐2B borehole and a variety of receivers deployed at the surface. The data set enables us to derive a detailed 3‐D P wave velocity model by first‐arrival traveltime tomography. We identify a 100–460 m thick sediment layer (mean velocity 2,200 ± 400 m/s) above the basement (mean velocity 4,200 ± 500 m/s). Particularly on the western valley side, a region of high velocities rises steeply to the surface and mimics the topography. We interpret this to be the infilled flank of the glacial valley that has been eroded into the basement. In general, the 3‐D structures revealed by the velocity model on the hanging wall of the Alpine Fault correlate well with the surface topography and borehole findings. As a reliable velocity model is not only valuable in itself but also crucial for static corrections and migration algorithms, the Whataroa Valley P wave velocity model we have derived will be of great importance for ongoing seismic imaging. Our results highlight the importance of 3‐D seismic data for investigating glacial valley structures in general and the Alpine Fault and adjacent structures in particular.


Geological Setting
The Alpine Fault is a major plate-bounding fault that accommodates ∼80% of the relative motion between the Australian and Pacific plates in the western South Island of New Zealand (Little et al., 2002;Norris & Cooper, 2001;Sutherland et al., 2007). Recent paleoseismological research reveals that the central and southern portions of the Alpine Fault in proximity to the Whataroa Valley rupture in M W 7-8 earthquakes remarkably consistently with average recurrence intervals of 263 ± 68 and 291 ± 23 years, respectively (Cochran et al., 2017;Howarth et al., 2018). The last major earthquake on the Alpine Fault occurred in 1717 A.D., and the conditional probability of a M W 8 earthquake in the next 50 years exceeds 27-29% (Biasi et al., 2015;Cochran et al., 2017;Howarth et al., 2018). A corollary is that the Alpine Fault is now very late in its typical interseismic phase.

The DFDP
In the last decade, the Whataroa Valley and surrounding areas have also been the focus of extensive study under the auspices of the Deep Fault Drilling Project (DFDP), a multidisciplinary investigation of the structure, present-day hydrogeological and mechanical state, and earthquake-generating characteristics of the fault Townend et al., 2009;Figure 1c). DFDP has to date involved two phases of scientific drilling targeting the Alpine Fault at depths of 100-150 m (DFDP-1, in 2011) and 1,000-1,500 m (DFDP-2, in 2014). During DFDP-1 drilling at Gaunt Creek (Sutherland et al., 2012), continuous lithologic sections extending from the hanging wall across the principal slip zone and into the footwall were sampled in core and wireline logging data sets Toy et al., 2015). Fault zone-guided waves also indicate a low-velocity zone extending to depths of ∼8 km (Eccles et al., 2015).
The second phase of drilling, DFDP-2, was intended to intersect the Alpine Fault at a depth of ∼1 km to determine temperature, stress, and fluid pressure regimes Townend et al., 2017;Toy et al., 2017). The borehole is nearly vertical from the surface to about 400 m depth but is deviated by as much as 44 • to the northwest at greater depths Sutherland et al., 2015). Due to this deviation, the true vertical depth of the 893 m long borehole is ∼820 m and the bottom of the borehole lies 720 m below sea level. All depth and height values in the following are referring to depth below sea level (bsl), that is, with the positive z axis pointing downward.
Within the Whataroa Valley, an extremely high geothermal gradient was found indicating active hydrothermal flow Janku-Capova et al., 2018;Sutherland et al., 2017;Townend et al., 2017). A detailed analysis of the rocks along the borehole based on drill cuttings was undertaken by Toy et al. (2017). A major boundary between sediments and basement occurs at a depth of 145 m. Within the basement, the grade of mylonitization increases with depth and Toy et al. (2017) identified layers of schist, protomylonites, and mylonites.
Technical problems encountered during the casing of the DFDP-2B borehole brought the drilling program to a premature halt, before the principal slip zone of the Alpine Fault had been reached . On the basis of the borehole stratigraphy, the borehole is thought to have terminated 200−400 m above the principal slip zone. As a consequence of the casing damage, only the top 400 m of the borehole remains accessible. However, a fiber optic cable was successfully installed along the length of the borehole. The fiber has since been used to conduct repeated temperature surveys  and seismic imaging using a version of distributed acoustic sensing known as distributed heterodyne vibration sensing .
Here we use the fiber optic cable as a continuous seismic sensor in conjunction with a conventional borehole geophone array and a dense grid of surface sensors (10−20 m spacing) to image the detailed velocity structure in the vicinity of the Alpine Fault. The resulting 3-D survey geometry enables us to overcome some of the effects of complex subsurface structure, including out-of-plane effects arising from the steep valley walls and topography, and to thus better image the valley-filling sedimentary sequence and features associated with the Alpine Fault itself.  Langridge et al., 2018;Norris et al., 2013) and the DFDP drill sites (yellow circles). The Whataroa98 (blue), WhataDUSIE (gray), and newly acquired 3-D VSP (red) seismic data sets are represented. Detailed survey area shown in Figure 2 is marked.

Previous Geophysical Studies of the Whataroa Valley and Surroundings
The Whataroa Valley has been the site of detailed geophysical study for several decades. Active-source (e.g., SIGHT: Davey et al., 1998;Stern et al., 2007;Van Avendonk et al., 2004) and passive-source seismic studies (e.g., Bourguignon et al., 2015;Guo et al., 2017;McNab, 2017;Michailos et al., 2019, and references therein) have targeted the Alpine Fault at various depths within the study area, and generally provided information on P wave structure only. At depths greater than 5 km, the Alpine Fault zone is generally identified by lower velocities in comparison to the surrounding schist (e.g., Bourguignon et al., 2015;Davey et al., 1998;Scherwath et al., 2003;Smith et al., 1995;Stern et al., 2001Stern et al., , 2007Van Avendonk et al., 2004). At shallower depths within the Whataroa Valley itself, Garrick and Hatherton (1974) first derived seismic P wave velocities from a ∼4 km long 2-D refraction profile. Based on simple layered refraction analysis, they identified an approximately 100−250 m thick sedimentary layer (v p = 2,100-2,200 m/s) above the basement (v p = 4,800-5,200 m/s). A larger-scale P wave velocity model for the uppermost 3 km was derived by Davey (2010) from the Whataroa98 2-D reflection seismic data set. Davey (2010) also identified a thin sedimentary layer of only a few 100 m thickness above the basement within the valley. Similar results were obtained by Lepine (2016) who analyzed several 2-D seismic reflection lines and derived a sedimentary layer of 100-500 m thickness.
The WhataDUSIE 2-D profile was acquired in 2011 (Lay et al., 2016;Lukács, 2017). The P wave velocity model derived from first-arrival traveltimes of the Whataroa98 and WhataDUSIE 2-D data (Lay et al., 2016) suggested a deep 400-600 m thick sedimentary layer (v p = 2,300 m/s) on top of the basement (v p = 3,500-5,500 m/s). In order to satisfy all of the observed traveltimes, particularly early first arrivals and fading-out of refractions from the Whataroa98 refraction data, the inversion requires a low-velocity zone (with a downward decrease in velocity of ∼1,000 m/s) at a depth of z = 800-2,000 m. It is difficult to find a realistic geological explanation for this anomaly. It probably results from the limitations of the 2-D survey design within the more complicated 3-D geometrical structure of the valley. The anisotropy of the underlying schists and mylonites (Christensen & Okaya, 2007;Godfrey et al., 2000) may also have some effect.
In addition to the seismic surveys mentioned above, the Whataroa Valley has been studied using gravimetric methods focused on understanding the subsurface geometry of glacially eroded channels. Following preliminary work by Brikke (2007), Davy et al. (2013) (2017) using ambient noise analysis and shear wave velocity modeling near the drill site. . Survey layout of the 3-D VSP experiment located around the DFDP-2 drill site. The local coordinate system used for the analysis is illustrated. The locations of receivers and sources are marked for the experiments VSP1 and VSP2. The Alpine Fault trace is formed by two distinct lines Toy et al., 2017) where the northern line marks the frontal trace and the southern line marks the dextral trace.

Seismic Data Acquisition
Generally, a 3-D seismic survey is much more expensive and complex to do than a 2-D seismic survey. Thus, many surveys have only the limited information available from a more or less straight 2-D-line. However, for complex geological circumstances and rugged mountainous topography, 3-D surveys can provide better and also more reliable information. A previous 2-D seismic data analysis within the glacially eroded Whataroa Valley (Lay et al., 2016) left questions unanswered concerning the underlying geological structure; particularly concerning the expected effects of the 3-D nature of the Whataroa glacial valley on both the velocity model and seismic images.
Here, we use the 3-D-seismic data collected with surface sensors in conjunction with Vertical Seismic Profiling (VSP) data to push forward the geophysical and geological understanding of the Whataroa Valley. For this, we derive a detailed 3-D P wave velocity model from first-arrival traveltime tomography.

Survey Design
The extended 3-D VSP data set was acquired in January and Feburary 2016 (see technical report by Townend et al., 2016). Different source and receiver combinations were used to address a range of complementary technical objectives ( Figure 2). In particular, the seismic acquisition was designed to yield a robust 3-D tomographic P wave velocity model and a detailed seismic image, given the constraints imposed by topography and the presence of a large river.
Throughout the survey, a vibration source (IVI EnviroVibe) with 6,800 kg peak force was used with frequencies sweeping from 10-200 Hz (Table 1). An accelerated weight drop with a mass of 250 kg positioned near the top of the borehole was also used to provide impulsive shots with which to determine the orientations of the borehole geophones.
A variety of borehole and surface receivers were used to simultaneously record the seismic signals (Table 1). All source and receiver locations were accurately determined with differential GNSS.
Within the DFDP-2B borehole, the fiber optic cable recorded along the complete length of the deviated DFDP-2B borehole. The heterodyne Distributed Vibration Sensing (hDVS) system developed by Schlumberger (Hartog et al., 2014) was used to record longitudinal strains along the fiber optic cable with an averaging length (i.e., gauge length) of 20 m. Data were extracted every 2 m along the borehole at a 1 ms sampling period. Depth calibration of the measurement positions was performed using a freezing test on the fiber optic cable.
In addition to the fiber optic system, a three-component wall-locking sonde system (Sercel SlimWave) with four levels of 15 Hz geophones at 15 m spacing sampled at 0.5 ms was used in the accessible 400 m of the borehole. This sonde was moved successively to achieve a 1 m spacing for a zero-offset VSP profile and a coarser spacing later in the survey.
Several different receiver systems were used at the surface (Table 1). A total of 412 one-component 10 Hz SM24 I/O geophones were deployed with an Aries reflection seismic system along two fixed lines subparallel to the Whataroa River and sampled at 1 ms. In addition, 160 autonomous three-component seismic recorders (cubes) acquired data continuously at a sampling period of 2.5 ms (Omnirec DATA-CUBE with 4.5 Hz Oyo-Geospace geophones), using two different acquisition geometries. The cubes were first distributed sparsely across the whole valley area with an average spacing of 35 m (red crosses in Figure 2; see section 2.2 VSP1). Next, they were successively repositioned to densely cover the central part of the valley in a grid with an average spacing of 10 to 20 min a rolling mode (yellow circles in Figure 2, see section 2.2 VSP2).
Finally, five three-component Reftek stations (Reftek 130 recorders with Geospace HS1-3d 2 Hz seismometers, 2 ms sampling interval) and 36 one-component cube receivers (DATA-CUBE, 4.5 Hz Oyo-Geospace geophone, 1.25 ms sampling interval) recorded continuously throughout the ∼10 day duration of the experiment. Data from these sensors were intended to be used primarily to study earthquakes, but have also played a minor role in the analysis presented here.

Principal Data Sets
Different combinations of shots and receivers were used to acquire data sets for specific purposes, which we introduce here. All presented data sets were used in this study to derive a detailed 3-D P wave velocity model.

VSP0: Zero-Offset VSP
A zero-offset VSP survey was undertaken to obtain a detailed image of the subsurface surrounding the borehole and to be able to directly correlate with borehole logging results (e.g., the sonic log results of Jeppson, 2017;Townend et al., 2017) and reflection seismic images. Additionally, we are able to compare both methods (DVS and conventional three-component recordings), extending the preliminary analysis conducted by Constantinou et al. (2016). The source point was positioned 23 m southeast of the top of the borehole. For the zero-offset VSP experiment, a 16 s sweep from 10-200 Hz(see Table 1) was used. The three-component Sercel tool was shifted to allow the wavefield to be obtained at 311 points from z = 1to 299 m depth (

Journal of Geophysical Research: Solid Earth
10.1029/2019JB018519

VSP1: Multiazimuth Source Lines and Extended VSP Experiment
Following the zero-offset VSP experiment, multiazimuthal and multioffset source lines were used to obtain better 3-D sampling of the subsurface structure. A stack of four sweeps, each 16 s long and spanning frequencies of 10-150 Hz (Table 1), was used. The respective source lines are marked by green stars in Figure 2. From the central borehole location, they extend as far as possible within the limited valley area along three main source lines, two parallel to and one across the valley.
The vertical-component Aries receivers were deployed on the two along-valley source lines with 10 m spacing (purple triangles in Figure 2). Additional three-component receivers (cubes) covered the accessible valley area during VSP1 (red crosses in Figure 2).
Making use of the multiazimuthal source lines, the three-component Sercel tool recorded a walkaway VSP. The tool occupied different depths (see the supporting information and details in Townend et al., 2016). For each offset source location the number of recorded channels was between 4 and 39, recording at depths between z = 70 and 295 m.
The hDVS system also recorded each shot but it was found that at least 10 sweeps needed to be stacked to obtain acceptable signal quality due to the lower sensitivity of the hDVS/fiber system in comparison to geophones. This could not be accommodated consistently within the limited time frame available for the field work. A compromise using stacks of 20-50 sweeps at lower frequencies of 10-60 Hz from a subset of locations was trialed. Ultimately, 10 shots from source locations in addition to the zero-offset source location were recorded by the hDVS and are used in this study for the velocity analysis.

VSP2: Dense Coverage of Three-Component Surface Geophones
For the second phase of the experiment (VSP2), the three-component receivers (cubes) were moved successively along the valley. The 160 cubes were deployed in a nearly rectangular array with a long axis subparallel to the strike of the Alpine Fault ( Figure 2) and 10 m (NNW-SSE) and 20 m (WSE-ENE) grid spacing. The whole array was moved 12 times to record reflections from the Alpine Fault zone over a broad depth range. Altogether, data were recorded at 1,916 receiver locations for each of 71 source locations (Table 1). Thus, the detailed 3C array densely covered an area within the Whataroa Valley of approximately 1,800 m inline along the river (i.e., perpendicular to the fault strike) and 600 m crossline perpendicular to the river (i.e., parallel to the fault strike).
The VSP2 source used the same 16 s sweep as in the VSP1 experiment with frequencies of 10-150 Hz ( Table 1). The 71 source locations formed a loop (cyan stars in Figure 2) and covered the central part of the valley.

Data Preprocessing
Some of the sensors (i.e., three-component cubes) recorded in a continuous mode. The GNSS-shot-time was later used to extract data in an appropriately short time window. The corresponding shot time was recorded using a one-component DSS cube (1.25 ms sampling interval). The recorded raw vibroseis data were subsequently correlated with the source signal of the sweep in order to extract the Earth's response function. For this procedure, the pilot sweep was used. After correlation, the data were stacked in order to increase the signal-to-noise ratio. For the three-component geophones (Sercel slimwave tool and Omnirec three-component cubes), vertical component data were used.
All acquired data satisfy the sampling theorem so that even for the largest sampling interval of 2.5 ms (cubes), the highest excited frequencies (200 Hz for the zero offset) do not exceed the Nyquist frequency of 200 Hz. The data were homogenized for picking, that is, resampled to the lowest sampling interval of 0.5 ms.
Examples of the recorded seismic data are shown in Figure 3.

Modeling Strategy and Resolution
3-D first-arrival traveltime tomography was applied to data sets VSP0, VSP1, and VSP2 with receivers at the surface and in the borehole. The inversion algorithm is based on the work of Zhang and Toksöz (1998). Ray tracing through a regular gridded model was conducted to calculate traveltimes. The minimized objective function includes the misfit between the picked and calculated first-arrival traveltimes as integrated slowness, the misfit of the gradients of traveltimes as apparent slowness and a regularization constraining the model roughness using a smoothing factor tau (see equation (1) in Zhang & Toksöz, 1998). Tests with different grid sizes, starting velocity models and smoothing parameters were undertaken (Bodenburg, 2017; see also the supporting information illustrating the results of various tests). Finally, a constant smoothing parameter of = 0.5 and 40 iterations were used for all inversions. The final results obtained by Lay et al. (2016) were expanded to 3-D with constant velocities along the y axis to span the width of the valley to produce a starting model. In addition to the parameter tests and analysis of the ray coverage, the inversion was first run separately for different subsets of the data as a variation of the jackknife strategy (Rawlinson et al., 2014). All obtained results (see examples in the supporting information) confirm the robustness of the main features.
Spatial resolution is strongly dependent on seismic wavelength. Most energy occurs in the range of 10-40 Hz with a dominant frequency of 25 Hz. With velocities in the range of 1,000-5,000 m/s, respective wavelengths are between 40 and 200 m. Hence, both vertical and horizontal spatial resolution cannot be expected to be better than 20 m. Consequently, a grid with cell sizes of 20 min the x, y, and z directions was adopted. This grid also derived the best results showing convergence during inversion and yielding the smallest RMS-misfits (Bodenburg, 2017). Moreover, this grid size is reasonable for a tomographic first-arrival traveltime inversion given a geophone spacing along lines of >10 m but much larger across lines.
Although the survey layout and resolution capacity is different for varying receivers (summary in Table 1), the overall goal of our analysis is to derive a common velocity model for all available data. Thus, the model uses the coarsest model resolution of 20 m at best from cubes. The vertical resolution along the borehole can be expected to be a lot better with a geophone spacing as low as 1 m for the Sercel tool. This is taken into account by separate analysis, for example, by zero-offset velocity profiles from Constantinou et al. (2016).
In the obtained blended velocity model, small-scale heterogeneities will not be resolved but will contribute to a mismatch between picked and modeled traveltimes. However, the strength of the blended model is to explain the overall 3-D velocity structures in the Whataroa Valley.
During the inversion, no a priori information except the starting model was used. The use of a layer-based approach was tested but not applied for the final models as the boundaries could not be reliably determined. Thus, the tomographic inversion tends to smooth velocity contrasts in order to minimize the RMS-errors because no boundaries are introduced as a priori information. Due to this smearing effect, high velocity contrasts, such as these expected at major geological boundaries like the sediment-basement-contact, will occur as (steep) velocity gradients in the final P wave velocity model. In the following, both borehole and surface data sets are analyzed separately but the complete data set was used for the final interpretation.

First-Arrival Traveltime Picking
First-arrival traveltimes (first maximum) were picked for all shots and receivers. As we have a sweep source signal, the first-arrival signals resemble the Klauder wavelet. Thus, the maximum of the wavelet was picked (see insets in Figures 3e and 3f). The wavelet was not always perfect so we picked as consistently as possible.
Only reliable picks that showed consistency with neighboring shots were used. Although there is a dependency of the first-arrival traveltimes on attenuation and dispersion, Molyneux and Schmitt (2000) show a discrepancy of only 3% for the group velocity determined at the predominant frequency from laboratory analysis. After quality analysis, 225,487 out of 364,625 picks were used for the tomographic inversion. Only 10,436 of these, that is, less than 5% , were recorded within the borehole. The horizontal source-receiver offsets range between 0 and 3,740 m.
First, we carry out an inversion using the first-arrival picks of only the borehole receivers. We then use the resulting model to calculate forward-modeled picks. In this step we include all picks (i.e., those from surface receivers as well) to compare with later inversion results. Second, we do an inversion using only the first-arrival picks of the surface receivers and again calculate the complete set of forward modeled picks. Finally, we do the inversion with all available picks and no weighting between borehole and receiver picks.  Borehole data from source location 2 are plotted against the absolute depth (below sea level) with Figure 3e showing the hDVS data and Figure 3f showing the vertical component of the Sercel data. Insets of Figures 3e  and 3f show where first-arrivals were picked. Note that the signal-to-noise-ratio in Figure 3f is better than that in Figure 3e, presumably due to the smaller offset. Already in this comparison, the good fit between the two types of recordings is visible as discussed by Constantinou et al. (2016) in more detail. Due to loose casing, the shallow part of the borehole shows strong ringing up to a depth of about 70 m and was not used for velocity analysis. For deeper parts, first-arrivals can be picked for both data types. Polarity of the hDVS data was reversed to make it consistent with the Sercel data. Interestingly, there is a clear boundary occurring at a depth of z = 140-170 m indicated by both changes in the slope of first-arrival traveltimes and reflections. At these depths, the drilling ) observed a major change from sedimentary to basement rocks and there is also a major change in casing diameters .

First-Arrival Traveltime Residuals
Traveltimes overlain on real seismic data in Figure 3 are analyzed in more detail with only traveltimes shown in Figures 4a-4d and their respective residual times (picked minus modeled time) in Figures 4e-4h. Except for outliers, the traveltime residuals of surface receivers are close to zero (< ±5 ms) for most offsets (cubes in Figure 4e and Aries in Figure 4f). The traveltimes/residuals for all data (blue) and data only recorded at the surface (turquoise) are nearly identical because of the overwhelming proportion of picks from surface receivers relative to the borehole receivers. Larger residuals (>40 ms) occur at larger (>800 m) and smaller (<100 m) offsets in Figure 4e. For the latter, they might be explained by a heterogeneous near surface, that cannot be sufficiently matched by the cell size of 20 m. In particular in Figures 4b and 4f, a directional dependency is evident so that two branches exist for offsets greater than 400 m.
Looking at the traveltimes from borehole receivers (Figures 4g and 4h), the general trend between real and forward-modeled picks is good for deeper areas (<5 ms for z > 170 m). At shallower depths (<170 m), the larger misfit (>10-20 ms) might again be explained by local geological heterogeneities around the borehole not resolvable by the 20 mgrid. As expected, using only surface receivers (turquoise curve) to calculate the velocity model yields the weakest fit, especially in the deeper part of the borehole where the surface receivers that determined the velocity model were far above the points used to compare to the model. In contrast, using borehole receivers only (olive curve) produces the best fit but also the complete data set (blue curve) shows a good fit that is only ∼5 ms different from the best fit.
A histogram counting the number of picks within a certain residual time window (Figure 4i) illustrates the residuals for the complete data set (i.e., ∼225,500 values). The different spreading of borehole residuals (olive) is obvious in contrast to the residuals from the almost identical surface recorders (turquoise) and complete data set (blue). For the model using all data (blue), the relative amount of data is quantified as well. Thirty-nine percent (i.e., >85,000 picks) of the residuals have an absolute residual as small as 5 ms.
Of the residuals 67% are within ±10 ms. These values are good in comparison to the pick accuracy lying between 0.5 ms for the best conditions and ∼10 ms for unfavorable conditions (no clear correlated data, weak signal-to-noise ratio, difficult coherency to neighboring traces especially for cubes in VSP1).
Total misfits (i.e., root-mean-square, RMS, values) serve as a quality measure and can be calculated between all available picked traveltimes and their respective forward modeled traveltimes for the portions of the data set. Using only the borehole receivers, we obtain an overall misfit of 29.9 ms with a standard deviation of ±36.7 ms. This relatively high misfit arises because the model is based on only 5% of the picks (i.e., the picks from borehole receivers) and thus cannot sufficiently fit all picks resulting in a nonnormal statistical distribution (Figure 4i). Using only the surface receivers, we obtain a misfit of 12.6 ± 26.3 ms. This is considerably better because surface-recorded data makes up more than 95% of the picks. However, the combined inversion yields an even better fit of 12.0 ± 26.0 ms. So we choose the resulting model as our preferred output that fits both the surface and borehole data.

Composite P Wave Velocity Model
The P wave velocity model with the respective ray coverage obtained using all available first-arrival traveltimes is shown in Figure 5. Areas without ray coverage are blanked out. The resulting velocity model shows ray coverage to a depth of nearly 800 m below sea level (i.e., about 900 m below the surface).
The depth xy -slice (Figure 5a) is located 50 m below the topography (mean central valley surface level at z = −100 m bsl). It shows a relatively homogeneous velocity of v p = 1,500-2,500 m/s for the major part of the covered area. Most interestingly, there are higher local velocities in the west at x = 2,900 m and y = 1, 100 m with velocities of more than 4,000 m/s coinciding with a good ray coverage in Figure 5c. Additionally, there is an anomaly with lower local P wave velocities of about 1,200 m/s at x = 2,700 and y = 900 m. This structure is less reliable due to a low ray coverage.
In the vertical xz slice (Figures 5b and 5d) the differences in ray coverage become obvious. The upper sediments (<100 m) are well covered in the far northwestern part by the densely spaced Aries receiver lines. In the most southeastern part, only a few cube receivers provide control. Closer to the borehole, the ray coverage is more complete on the southeastern side. This is a result of the survey geometry as there are more source locations in the southeast than in the northwest. Lower ray coverage around the borehole at z = 0-200 m is partly explained by 3-D effects (see supporting information) but also by the fact, that about 65% of the borehole first arrivals are from z > 200 m. Additionally, the Sercel borehole receivers only record to z = 300 m and only the hDVS sensor can record deeper. Unfortunately, as already noted due to the lower sensitivity of the hDVS, only 10 offset shots were recorded with sufficient quality of recognizable first arrivals for velocity model building. Thus, the ray coverage is lower for z > 300 m and will result in less constrained results.
Nevertheless, already from this image we can see that there are lower P wave velocities of 2,200-2,600 m/s in the upper few hundred meters. Underneath this sedimentary layer, higher basement velocities of 3,500-5,000 m/s dominate. Local anomalies occur within the sediment layer with velocities as low as 1,200 m/s at x = 1,700 and 2,700 m.
In general, the data confirm the expected sedimentary layer on top of the basement. The respective velocities with means of 2200±400 m/s and 4200±500 m/s lie within a suitable range for the expected rock types (Quarternary sediments and schists/mylonites, e.g., Carpenter et al., 2014;Christensen & Okaya, 2007;Schön, 2015) as well as previous studies (Davey, 2010;Garrick & Hatherton, 1974;Lay et al., 2016;Lepine, 2016). Furthermore, the data show that the velocity model within the Whataroa Valley is variable in 3-D.

P Wave Velocity Structure
In order to further investigate the variation in P wave velocities we have divided the 3-D model ( Figure 5) into a series of subvolumes to represent reliable velocity variations ( Figure 6, corresponding ray coverage in the supporting information Figure S10). Areas without ray coverage are not included. The mean P wave velocities are calculated by horizontally averaging the velocity in each subvolume onto 2-D panels that extend down the valley (Figures 6a1-6a4 on the left from north to south) and across the valley (Figures 6b1-6b4 on the right from east to west). A small sketch in each part of the figure shows the locations of the panels.
Two isovelocity surfaces are marked as dashed lines on each panel: v p = 1,500 m/s (gray) and v p = 2,600 m/s (black). These surfaces are interpreted on the basis of the well tie (see section 4.3) to mark a change in sedimentary layers and the top of the basement, respectively. In all slices of the P wave velocity model, the general structure shows a lower-velocity sedimentary layer on top of the higher-velocity basement. In general, the eastern part of the model is relatively homogeneous and thicker subvolumes are formed (e.g., Figures 6a1 and 6b2). The largest variations in velocity appear as we approach the western edge of the valley (Figures 6a3 and 6b2). In particular, the top of the basement (black dashed line) is rougher and the velocity at depth is higher.

Sedimentary Features
A layer of low P wave velocities (mean of 2,200 ± 400 m/s for the complete area) is associated with Quaternary sediments. The lower boundary is defined by the isovelocity surface of v p = 2,600 m/s (black dashed line in Figure 6). Within the sediments, shallow near-surface basins occur with velocities as low as 1,050-1,500 m/s in the upper 50 m (gray dashed lines in Figure 6). The deepest low-velocity basin occurs at x = 2,700 m ( Figure 6a3) and has a thickness of 100 m. These basins are located close to the current location of the Whataroa River (see Figure 2). Thus, these basins may have arisen from recent fluvial sedimentation bringing in coarser, less consolidated material such as gravel. A top layer consisting of gravel and sand instead of silt is also supported by the local geology derived from the borehole  and excavations of trenches designed for paleoseismology studies .
In particular, there might also be a correlation with Alpine Fault-related structures. Figure 6a3 shows the presence of two basins with lower P wave velocities (between x = 1,600-2,000 and x = 2,400-2,900 m). The southern boundary of the northern sedimentary basin (x = 2,000 m) coincides with the location of the surface (frontal) trace of the Alpine Fault . The area that separates the two basins has higher velocities over a width of 400 m wide and is coincident with the gap between the mapped frontal and projected dextral fault (see Figure 2Langridge et al., 2018;Toy et al., 2017).

Basement Features
High P wave velocities (between 2,600 and 5,000 m/s; mean of 4,200 ± 500 m/s for the complete area) are associated with the basement. The estimated depth to the top of the basement (black dashed lines in Figure 6) shows significant 3-D variation along and across the valley. Sediment thickness varies between 100 and 460 m (at z = −90 in Figure 6b2 to z = 360 m in Figure 6a1 with topography at z = −190 to z = −100 m, respectively). We assume that the basement steepens on the western side of the valley, especially adjacent to the surface location of the borehole. The high velocities present at shallow depths indicate that we are imaging the western flank of the valley. On the eastern side, we cannot identify a similar pattern as the receiver and source arrays do not extend sufficiently to cover the limit of the valley. The floor of the glacial valley, as represented by the 2,600 m/s isovelocity surface, undulates between 60 and 350 m bsl along the current Whataroa Valley (Figure 6a2). Figure 7 shows a comparison of the inverted P wave velocity model with the results from the borehole Jeppson, 2017;Townend et al., 2017;Toy et al., 2017). From the 3-D P wave velocity model (Figures 5 and 6) the velocity values are extracted along the borehole trajectory.

P Wave Velocity Profiles From Tomography
In Figure 7a, the P wave velocities extracted from the tomographic model are shown. As above in section 3, three different scenarios are plotted using (1) all available receivers for the inversion (red curve), (2) only borehole receivers (olive curve), and (3) only surface receivers (turquoise curve). In all models a sharp velocity contrast, such as the silt to schist contact (Figure 7b), will manifest itself as a change in gradient of the velocity. Smearing occurs in the P wave velocity tomography because of ray tracing and grid-based inversion algorithms using a 20 m grid spacing.
Looking at the model with only surface receivers (turquoise curve), the top section (z < −35 m) shows a relatively steep gradient with velocities of 1,200-2,200 m/s. At greater depth, the gradient decreases and the velocity rises to 5,000 m/s at the base of the model with changes in the gradient at depths of around 300 and 650 m.
In contrast, the model with only borehole receivers (olive curve) shows a constant but steep gradient in the shallow part (<100 m) but shows more variations at greater depths. Between 220 and 350 m depths the velocity is relatively constant at 3,800 m/s. As depth increases, the velocity increases to 5,000 m/s at 500 m and stays constant until 650 m depth, where there is a slight decrease.
Considering all receivers, the red curve shows similar character to the turquoise surface receiver model in the top 100 m (Figure 7a). From 50 to 200 m the combined model is dominated by the borehole receivers. Below a depth of 400 m, the velocity is effectively constant to the bottom of the model at v p = 5,300 m/s.

10.1029/2019JB018519
Generally, below 240 m greater variations within the combined model (red) and between different models occur and are caused by low ray coverage at these depths (see Figure 5d). However, in the red curve we see a sharper definition of the boundaries that occur as changes in the gradient.
Additionally, anisotropy can be expected for P wave velocities as foliated schist is highly anisotropic (e.g., Allen et al., 2017;Christensen & Okaya, 2007;Godfrey et al., 2000;Li et al., 2018;Okaya et al., 1995) and could explain traveltime differences between differing main ray travel paths for borehole and surface receivers. In a similar experiment in Sweden, Simon et al. (2017) were able to determine Thomsen parameters and thus an anisotropic velocity model (vertical transverse isotropy, VTI) making use of the different dominant ray paths in surface and borehole analysis. However, the Alpine Fault in the Whataroa Valley dips at ∼50 • (Lay et al., 2016;Toy et al., 2017) and the fabric dip is between 45 • and 55 • Townend et al., 2017) so that the simplest anisotropy approach would be a tilted transverse isotropy case. This tilt makes the derivation of anisotropy parameters more challenging and it is not possible to only use the described differences in traveltimes observed by surface and borehole receivers. Although challenging, three-component and borehole data will be investigated further to gain more information and possibly derive anisotropy parameters. Simultaneously, rock samples are analyzed in detail to quantify the geometry of the anisotropy (Jeppson, 2017;Li et al., 2018). Figure 7b shows a simplified geological column along the borehole (from Toy et al., 2017 and. No core was collected from the borehole, but cuttings were analyzed throughout the drilling . The surficial sediments can be divided into three units, gravels from surface (−95 m) to a depth of −55 m (yellow), sands to a depth of −35 m(light orange), and silt to a depth of 145 m (dark orange).

Geological Interpretation
Close to the top of the metamorphic basement a second thin gravel layer was encountered. The contact at 145 m depth is a major change in all petrophysical properties as the upper part of the basement is a 40 m thick unit of schist. The basement consists mainly of protomylonites derived from schist and the grade of mylonitization increases with depth reaching true mylonites at ∼700 m .
Comparing this geological description to the tomographic P wave velocity model in Figure 7a (red curve) shows that major boundaries at −55 and 145 m correlate well with steep gradients in v p . First, within the well-resolved shallow part of the P wave velocity model, the change between sand and silt (at −35 m) can clearly be identified as a step with a steep velocity gradient underlain by roughly constant velocities of 2,000-2,200 m/s. The shallower change between gravel and sand at a depth of −55 m is not so pronounced, likely due to the sand layer being the same thickness as the vertical cell size in the tomographic model.
Second, directly at the silt/schist boundary a steep velocity gradient starts in the P wave velocity model (Figure 7a, red curve). The velocity at this inflection point is approximately 2,600 m/s. We have used the 2,600 m/s boundary to indicate the top of basement ( Figure 6).
There might also be a correlation at depths >700 m, interpreted to be close to the Alpine Fault principal slip zone, where a slight decrease in velocities is predicted by the tomographic model (Figure 7a). At the Alpine Fault, the mylonites in the hanging wall are more fractured so that the P wave velocity should decrease. However, these values are derived from only a few ray paths and are close to the maximum penetration depth of the inversion model.

P Wave Velocity Profiles From Tomography, Zero-Offset VSP, and Sonic Log
The final combined models from this study (red curve) and from the previous 2-D study (orange curve) are presented in Figure 7c compared with other estimates of the P wave velocity in the well. In the upper 100 m, both models coincide well as this area is dominated by surface receivers and is well resolved by a high density of ray paths. For deeper areas (>300 m), the two models differ but are within the same range of 4,000-5,500 m/s.
The results from the zero-offset velocity analysis (adapted from Constantinou et al., 2016) are also shown in Figure 7c. These models are derived from the same first breaks used in the 3-D tomographic inversion but are computed assuming a vertically incident downgoing wave from the surface source location at the well head. The plot shows models for the conventional Sercel borehole tool (green curve) and the hDVS fiber optic cable (blue curve) separately. In general, the correlation with the tomographic model (red curve) is very good. However, within the upper part of the basement, both zero-offset methods show significantly higher velocities than the 3-D tomography. These unrealistic velocities exceeding 6,000 m/s may be related to poorly resolved first arrivals caused by refraction effects, noise due to casing, and interference with upgoing waveforms . In the shallow, the zero-offset VSP data show a short section of realistic velocities of 1,500-2,500 m/s, which correlates well with the 3-D tomography model, and is consistent with silts.
Results obtained by Jeppson (2017) from the full-waveform sonic log acquired during the drilling are shown in the background of Figure 7c (gray curve). The full-waveform sonic data are complicated with weak first arrivals and low signal-to-noise ratios so the derived P wave velocity is scattered with no clear trends. The overall velocity ranges from 4,000 to 5,500 m/s which correlates with the velocities derived from the VSP and the 3-D tomographic models (Figure 7c). Despite their different characteristics, the overall P wave velocities along the DFDP-2B borehole show remarkable similarities. Figure

Comparison With the Previous Tomographic Velocity Model
To extend on section 4.3.3, both tomographic velocity models from the 2-D (Lay et al., 2016) and 3-D seismic data (this study) are compared in more detail in Figure 9. Within the presented 3-D study, more than 225,000 valid first-arrival traveltime picks from surface and borehole receivers for more than 400 shot locations were used for the 3-D tomographic inversion.
In contrast, for the combined 2-D WhataDUSIE and Whataroa98 tomography, significantly fewer first-arrival traveltime picks were used, that is, about 58,000 in total (Lay et al., 2016). Receivers and shots were located in the central part of the valley and were aligned along the 2-D profile. Although the tomography itself was done in 3-D, the 3-D effects of rays traveling through the faster valley flanks could not be resolved. As a result, a 2-D slice was extracted along the profile that was used as the final P wave velocity model. The upper 1 km of the P wave velocity model from Lay et al. (2016) is shown in Figure 9a. Most of the velocity information is derived from the long-offset Whataroa98 data (Davey, 2010). Within the area marked by the dashed line, additional information from the WhataDUSIE data set is incorporated. The location of the surface trace of the Alpine Fault and the borehole trajectory are included. The model shows a sedimentary layer (greenish colors) above the basement (orange to red). Interestingly, there is a layer of high P wave velocities at a depth of about 500 m. This layer was discussed in detail by Lay et al. (2016).
With the new more detailed 3-D velocity model, we can create a similar summary 2-D slice (see Figure 9b). As seen in Figure 6, with the help of the 3-D data, we can distinguish between contributions from the valley flanks and the central part of the valley. Hence, we use the eastern part of the model (y = 0-600 m) to create a representative model, constructed from the means of available velocities, for the central part of the valley with thick sediments. Most interestingly, no high-velocity layer within the basement is observed. Thus, we conclude that the previously measured high-velocity zone is caused by a 3-D effect. Rays traveling through the fast valley flank (especially in the west) and rays traveling through the slower central sediment layers could not be distinguished by the 2-D seismic data. Thus, both effects were mixed and caused the apparent high-velocity layer.
The difference between the new (Figure 9b) and the old (Figure 9a) velocity models is shown in Figure 9c. Within the shallow subsurface (upper 100-200 m), the differences are minor. Thus, we can conclude that both models create a reliable near-surface velocity model. However, we can assume that the new 3-D model is slightly better because there are more first-arrival traveltime picks available. Also, the new extracted profile (Figure 9b) is constructed from a more regular distribution of raypaths (averaged over a width of 600 m) rather than unbalanced raypaths from the crooked line geometry used previously (Figure 9a). Generally, the top of the basement (velocities exceeding 2,600 m/s) lies slightly deeper in the new model with the only exception being in the southeast at x = 4, 500 m, where there is a region with a shallower top of basement. At greater depths, major differences as high as 1,000 m/s occur. They are predominantly caused by the absence of the high-velocity layer in the new velocity model. In contrast, high-velocity valley flanks are identified in the 3-D model that could not be imaged in a central 2-D slice.
Hence, this study shows for the first time how 3-D structures within the Whataroa Valley have a great influence on the P wave velocity model. Furthermore, these structures can only be imaged successfully by 3-D seismic data.

Subglacial Valley Structures
Our observations of the 3-D structures in the P wave velocity model correspond to the origin of the Whataroa Valley as glacially formed and molded (Barrell, 2011;Lukács, 2017;Suggate, 1990). During the Last Glacial Maximum (27,000-20,000 years ago), the valley glacier reached its maximum extent, carving out the present-day valley and extending out on the coastal plain to the northwest (Barrell, 2011;Suggate, 1990). Afterward, the valley was successively filled by glacial, lacustrine or marine, and fluvial sediments (Thomas, 2018).
Across the valley, we can identify the typical U shape of a glacial valley with thicker sediment layers in the central part of the valley and thinner sediment layers at the valley flanks. Along the valley from north to south, additional variations in the depth of the basement undulate and become shallower to the south at x = 4, 500 m (Figures 6a1 and 6a2). These undulations are smooth but highly variable in 3-D and thus might be caused by former glacial erosion. Similar features are also observed by other seismic studies in glacial valleys (e.g., Ahmad et al., 2009;Büker et al., 2000Büker et al., , 2010Burschil et al., 2018Burschil et al., , 2019Lukács et al., 2018;Ogunsuyi & Schmitt, 2010).
The final 3-D P wave velocity model reliably fits the geological expectations of a glacially formed valley and gives us further insight into the 3-D structures of the Whataroa Valley. The sediment thickness of 100-460 m correlates not only with other seismic studies (Davey, 2010;Garrick & Hatherton, 1974;Lay et al., 2016;Lepine, 2016;Lukács, 2017;McNab, 2017) but also with gravity results (Jenkins et al., 2019).

Alpine Fault-Related Features
The Alpine Fault is expected to be a region of relatively low P wave velocity (Bourguignon et al., 2015;Carpenter et al., 2014;Davey et al., 1998;Stern et al., 2001Stern et al., , 2007Townend et al., 2013;Van Avendonk et al., 2004). However, there is no indication of a low-velocity anomaly corresponding to the Alpine Fault in our final P wave velocity model. The location of the surface trace of the Alpine Fault is only at the northwestern margin of the studied area ( Figure 6) but is expected to dip to the southeast. The maximum depth range of traveltime paths yielding velocity information is at z = 600 m. Thus, we assume that we either cannot capture the Alpine Fault as a low-velocity structure due to the limitations of the presented P wave velocity model or, more likely, the Alpine Fault is simply not developed as a strong low-velocity zone in the upper 700 m.
Shallow sediments may show some direct evidence of the fault. Two shallow basins with decreased velocities are imaged that are presumably filled by young, unconsolidated sediments. The border of the northern basin and the surface trace of the Alpine Fault coincide (Figure 6a3). Furthermore, Figure 6a3 shows an abrupt offset in the top of the basement at x = 2, 300 m. This feature might be related to fault structures, but it occurs in an area with insufficient ray coverage and would thus need further investigation. However, the identified velocity structures may suggest that the recent Alpine Fault location is also captured in seismically imageable sedimentary processes as seen in surface trenches by Langridge et al. (2018).

Conclusions
The newly acquired extended 3-D VSP data set includes a variety of surface and borehole (fiber optic cable and three-component geophones) receivers. From this 3-D VSP data we derived a reliable 3-D P wave velocity model at the DFDP-2 drill site within the Whataroa Valley. The applied first-arrival traveltime tomography made use of nearly 225,500 first-arrival traveltime picks.
Only by using a 3-D seismic data set have we been able to adequately constrain 3-D valley structures. Thus, we recommend the use of 3-D seismic data within complex geological environments such as glacial valleys that might otherwise contaminate 2-D data with out-of-plane 3-D signals. In the case of the Whataroa Valley, we successfully identified the 3-D structures.
Generally, two main facies are distinguished in the velocity model consisting of a sediment layer covering the basement. The sediment layer is 100-460 m thick and shows P wave velocities of 1,050-2,600 m/s with a mean of 2,200 ± 400 m/s. Lower velocities are found in the shallow subsurface and form up to 100 m thick basin-like structures of likely unconsolidated sediments. This fits well to the gravel, sand, and silt layers identified from samples and borehole findings, as well as other geophysical work.
The basement below consists of Alpine schist and (proto)mylonites, which is consistent with the P wave velocities obtained here of 4,000-5,000 m/s with a mean of 4,200 ± 500 m/s. Additionally, the top of the basement is defined as the beginning of a steep velocity gradient that starts at v p = 2, 600 m/s. The top of the basement depicts the preinfill 3-D structure of the Whataroa Valley. Along the valley, there are slight undulations interpreted to result from variable glacial carving. Across the valley, we can identify the typical shape of a glacially carved valley later filled with sediments. Within the central part of the valley, the sediments form a thick layer (up to 460 m thick) whereas at the flanks they become thinner (as thin as about 100 m).
High velocities are observed in shallow areas in the west that are interpreted to image the steep valley flanks. There is a correlation between the P wave velocities marking the basement and the topography with western steep mountain flanks above the surface. Also, the identified 3-D structures are consistent with the geological history of the Whataroa Valley as a glacially formed and molded valley.
In spite of the high resolution of the 3-D survey, no low-velocity region was clearly obtained in the vicinity of the Alpine Fault. However, there are weaker features such as velocity contrast in the sedimentary layer possibly indicating fault-related structures.
The derived P wave velocity model is valuable in its own but can also serve as a reliable basis for further (reflection) seismic data analysis. As the 3-D P wave velocity model enables us to distinguish seismic reflections/structures from glacial valley structures and reflections from the Alpine Fault, we expect a detailed seismic image of the Alpine Fault at the DFDP-2 drill site.