3D Active Source Seismic Imaging of the Alpine Fault Zone and the Whataroa Glacial Valley in New Zealand

The Alpine Fault zone in New Zealand marks a major transpressional plate boundary that is late in its typical earthquake cycle. Understanding the subsurface structures is crucial to understand the tectonic processes taking place. A unique seismic survey including 2D lines, a 3D array, and borehole recordings, has been performed in the Whataroa Valley and provides new insights into the Alpine Fault zone down to ∼2 km depth at the location of the Deep Fault Drilling Project (DFDP)‐2 drill site. Seismic images are obtained by focusing prestack depth migration approaches. Despite the challenging conditions for seismic imaging within a sediment filled glacial valley and steeply dipping valley flanks, several structures related to the valley itself as well as the tectonic fault system are imaged. A set of several reflectors dipping 40°–56° to the southeast are identified in a ∼600 m wide zone that is interpreted to be the minimum extent of the damage zone. Different approaches image one distinct reflector dipping at ∼40°, which is interpreted to be the main Alpine Fault reflector located only ∼100 m beneath the maximum drilled depth of the DFDP‐2B borehole. At shallower depths (z < 0.5 km), additional reflectors are identified as fault segments with generally steeper dips up to 56°. Additionally, a glacially over‐deepened trough with nearly horizontally layered sediments and a major fault (z < 0.5 km) are identified 0.5–1 km south of the DFDP‐2B borehole. Thus, a complex structural environment is seismically imaged and shows the complexity of the Alpine Fault at Whataroa.

The Alpine Fault Deep Fault Drilling Project (DFDP) consisted of two main drilling phases (Figure 1b, central part of the Alpine Fault) aiming to intersect with complementary depth sections of the Alpine Fault. The second phase (DFDP-2) took place in the Whataroa Valley and drilled the nearly 900 m long deviated DFDP-2B borehole to be used for multidisciplinary research (see overview by Sutherland et al. [2017]; Townend et al. [2017]; Toy et al. [2017] and technical report by Sutherland et al. [2015]). During the drilling, a fiber-optic cable was installed outside the borehole casing to be used for temperature measurements (Janku-Capova et al., 2018) and for seismic imaging using the distributed acoustic sensing (DAS) technique .
Around the drill site in the Whataroa Valley, several seismic experiments have been performed over the last two decades. At shallow depths (≤5 km), the reflectors are heterogeneous and difficult to identify but have similar dips ranging from 40° up to nearly 80° directly below the surface (King et al., 2020;Lay et al., 2016;Lepine, 2016;Lukács, 2017). At a greater depth of 15-30 km, reflectors dipping at 40°-60° to the southeast have been identified in active-source studies (Davey et al., 1998;Stern et al., 2007;Van Avendonk et al., 2004).
Velocity models have been derived on a small scale mainly by analyzing refracted waves from active-source data sets (Davey, 2010;Garrick & Hatherton, 1974;Lay et al., 2016Lay et al., , 2020Lukács, 2017) and show a 100-460 m thick sedimentary layers with slower P-wave velocities are identified above the higher-velocity basement. At a larger scale, earthquake records have been analyzed (e.g.., Boese et al., 2012;Chamberlain et al., 2017); Feenstra et al., 2016;Guo et al., 2017;Michailos et al., 201)) and generally show the Alpine Fault zone to have lower velocities in respect to surrounding rock (Bourguignon et al., 2015;Stern et al., 2007).  , deep fault drilling project drill sites, and active seismic data sets.

Seismic Data Set
In 2016, an extensive 3D combined surface-borehole seismic survey was conducted in the Whataora Valley at the DFDP-2 drill site (Lay et al., 2020;Townend et al., 2016) making use of the recently drilled DFDP-2B borehole. The main goal of the field experiment is to image the Alpine Fault zone in detail.
Throughout the survey, an EnviroVibe seismic source with a hold-down mass of 6,800 kg was used. A combination of receivers in the borehole (DAS-type heterodyne distributed vibration sensing (hDVS) and three-component conventional borehole tool by Sercel) and an array of seismic sensors distributed at the surface recorded the seismic wavefields. Figure 2a shows the complete survey layout. The DFDP-2 drill site is located in the central part of the survey and the Alpine Fault surface trace is shown as the frontal and dextral trace after Toy et al. (2017) and Langridge et al. (2018). In a first step, the seismic data were used to derive a 3D P-wave velocity model described in detail in Lay et al. (2020). In a second step, for imaging purposes, two subsets of the surface recordings named as VSP1 and VSP2 were used for the seismic imaging presented here. Table 1 summarizes the acquisition parameters important for the presented imaging. Full details are described in the technical report by Townend et al. (2016). In general, VSP1 consists of multi-azimuthal source lines combined with many fixed receiver locations whereas VSP2 aimed at a dense coverage with three-component receivers through a rolling mode of nodal receivers combined with fewer source locations explained in detail below.
For data subsets VSP1 and VSP2, a 20 s record length was used for a 16 s long source sweep followed by 4 s of listening time. The cosine-tapered up-sweep had a frequency range from 10-150 Hz. To enhance the signal-tonoise ratio, four sweeps were stacked at each source location.

of 21
Subset VSP1 consists of source lines with a variety of azimuths recorded in the borehole. The source points were located along three major profiles (dark green stars in Figure 2a). Additionally, geophone surface receivers recorded using an Aries system were located at the source lines parallel to the Whataroa River. The response of the single-component vertical 10 Hz geophones deployed along two 2D lines was recorded in a triggered mode with a sampling interval of 1 ms (Table 1). The total number of receivers used along both lines was 412, with 312 deployed along the western line and 100 along line 4000 further east. Both profiles cross at least the frontal trace of the inferred Alpine Fault surface traces.
Subset VSP2 consists of a dense 3D grid coverage with continuously recording three-component 4.5 Hz geophones connected to DSS Cube data loggers (Figure 2b). Given the vertical-force nature of the seismic source, only the vertical component has been used in the analysis. In total, 12 patches of 160 single nodal recorders (cubes) were deployed subsequently in a rolling mode with each patch recording all 71 source locations. The first array part is visible as the patch north of the Whataroa river in Figure 2. The following 11 patches lie south of the river. In total, 1916 receiver locations were used, each of them recorded all 71 source points. The whole 3D array is oriented so that the x-direction is approximately parallel to the valley (NNW-SSE direction) with a geophone spacing dx of 10 m while the y-direction is approximately parallel to the main strike of the Alpine Fault with a geophone spacing dy of 20 m (WSW-ENE). Local coordinates are defined accordingly (see maps in Figure 2a). The local x-axis is approximately perpendicular to the strike of the Alpine Fault. Local depth is expressed as depth below sea level, that is, with the positive z-axis pointing downwards. The average topography is about 100 m above sea level, that is, at z = −0.1 km. The main emphasis of the following analysis is on the Aries 2D line 2000 with receiver and source locations along the same profile and also on the cube 3D array.

Seismic Processing
A general goal of seismic reflection processing is to enhance reflected arrivals while minimizing noise (surface waves, refracted arrivals, etc.) through a variety of applied filters and methods, and then to construct subsurface images from these reflections. At first, conventional common midpoint (CMP) approaches were applied to the 2D Aries line 2000 by Hall et al. (2016). However, their results are challenging for interpretation but will be used to complement our seismic images. In contrast to the CMP approach, this paper focuses on the use of focusing prestack migration (PSDM) approaches for both VSP1 and VSP2.
In this work, the processing is divided into two sequential parts. First, the preprocessing part aims to improve the seismic records to reveal potential reflections and enhance the signal-to-noise ratio. Second, the imaging part is performed by using PSDM approaches incorporating a detailed 3D velocity model. Here, the term reflection is used when describing reflected arrivals in shot gathers in the time domain but the PSDM will also include scattering and diffractions. In contrast, the term reflector is used for coherent events in seismic depth sections obtained by migration, which relocates the recorded seismic reflections to their true reflector origin in the subsurface.

Seismic Preprocessing
Within the single-source gathers of the 3D cube array, the 3D data are sorted into xz-profiles (parallel to the river and thus perpendicular to the average strike of the Alpine Fault) and yz-profiles (crossing the river/valley, parallel to the average strike of the Alpine Fault). Tests of the output of different filters are analyzed on both yz-profiles (see Figures 3b and 3d) and xz-profiles (see Figures 3a and 3c).
Several filter types have been tested and applied to the complete 3D cube array data set. The final parameters are summarized in Table 2. The crucial filters improving the signal-to-noise ratio are spectral whitening, deconvolution, air blast attenuation and band-pass filtering. Additional tests with frequency wavenumber (fk) filters did not improve the 3D array data without introducing too many artifacts and were consequently not included in the final processing scheme. Finally, gain corrections were applied before the hand-picked top mute that removes the refracted first-arrivals. A similar preprocessing scheme was applied to the Aries 2D line 2000 by Backes (2020) who obtained better performance of the fk-filter due to the more regular spacial sampling. After the preprocessing sequence, the data were input into the prestack depth migration routine.  Table 2, this surface wave energy is clearly attenuated, as illustrated in Figures 3b and 3d. The air wave (gray arrow, AW) is well attenuated, too. Applying the preprocessing sequence including top mute enhances the signal quality and the signal strength overall so that the reflections are enhanced and clearly visible.

Optional Static Corrections
Static corrections were applied during the preprocessing as optional parameters ( Table 2). The tomostatics approach implemented as part of the software GeoTomo (Zhang & Toksöz, 1998) was used to determine the static corrections. The near-surface velocity model was calculated via turning-ray tomography using the complete 3D-VSP data set used for the tomographic inversion described in detail by Lay et al. (2020). This detailed near-surface velocity model was then used to compute the associated static corrections for source and receiver locations. This approach has advantages in areas with highly heterogeneous geology and in which the near-surface structure departs appreciably from a simple horizontally layered model and was thus chosen for the seismic data from the Whataroa Valley. The computed source and receiver statics (based on Lay et al. (2020) using GeoTomo by (Zhang & Toksöz, 1998)) were applied in the preprocessing sequence and, in general, make reflections appear more coherent. However, static corrections are particularly helpful for long, continuous reflections, which is not the case for the seismic data from a highly heterogeneous subsurface, such as a glacial valley. In this case, static corrections do not significantly improve the coherency of reflections. In general, as many filters as necessary but as few as possible should be applied. Thus, static corrections were tested but were not ultimately included in the final preprocessing sequence. Having only small elevation changes along the valley might also help not to need static corrections. Overall, the applied preprocessing summarized in Table 2 was able to generally improve the signal quality and to reveal several clear reflections in the source gathers.

Prestack Depth Migration
A crucial goal of seismic imaging, is to convert the recorded seismic wavefield into structural images of the subsurface. A widely used approach is Kirchhoff prestack depth migration (PSDM) based on ray tracing in the subsurface. Here, we apply the focusing Fresnel volume migration (FVM) (Buske et al., 2009)   quality for sparse data sets in particular with low coverage, limited aperture, and in the case of steeply-dipping reflectors.
Using FVM, we estimate the emergent angle for each time-sample at a given seismic trace using a slowness calculation after Hloušek and Buske (2016), in which neighboring traces are defined as a radial offset. The slantstack method applied to neighboring traces using the local velocity at the receiver yields inclination and azimuth. In this way, the total emergent angle is determined and used within FVM for back-propagation of the wavefield into the subsurface.

P-Wave Velocity Model
In prestack depth migration, the use of a sufficiently accurate velocity model is often crucial since the seismograms are transferred from the time domain to the depth domain. In the Whataroa Valley case, the migration is performed from the topography. Thus, the migration algorithm incorporates the near-surface velocity variations as well as the topographic differences in source and receiver locations so that there is no need for additional static corrections. Lay et al. (2020) derived a detailed 3D P-wave velocity model from the complete extended 3D-VSP data set. For imaging purposes, this velocity model is slightly simplified in order to minimize the potential of introducing artifacts through the usage of a full, more complex, 3D velocity model in areas that were not well constrained and were not covered by a sufficient number of rays during construction of the velocity model. This simplification includes the usage of a representative 2D slice along the x-direction that is used in the same way for all xz-profile across the entire area to create a 3D-volume that still includes the most important velocity variations in the near-surface. The resulting velocity model is a so called 2.5D model, whereby the chosen representative 2D slice (see Lay, 2021) is taken from the central part of the valley with the thickest sedimentary cover of more than 400 m. This approach ensures that relatively low velocities are preserved in the model. Imaged reflectors obtained by prestack depth migration are hence surely data-driven by seismic reflections and not introduced by nonphysical reflectors caused by focusing effects during the ray tracing.

Image Processing
The major goal of this study is to image the Alpine Fault zone. However, reflections from the fault zone are not clearly visible at first sight in the single-source gathers and further efforts are necessary to make sure that weak reflections are not missed during the preprocessing or imaging steps. For that purpose, we applied seismic image processing in the form of Gaussian filtering (see overview by Hall [2007]).
Before image processing, a 2D slice is created from the full 3D FVM seismic cube indicated in Figure 2a by summing all values in the y-direction through the seismic cube. The image processing is then applied to this slice. This means that structures present at all y-locations will stack constructively, but conversely, the 3D character in the y-direction is lost. The main goal of this procedure is to enhance potential fault-related structures that are nearly perpendicular to the local x-axis, which is the case for the Alpine Fault. As the reflectors obtained from this result are more smeared than in the pure migration image, both results are jointly interpreted. The Gaussian filter is implemented as a frequency filter in the frequency domain to emphasize structures depending on their wavelengths in the seismic image. In the Whataroa Valley case, the filter length of the Gaussian curve is chosen with a half-width of 100 m and longer wavelengths are suppressed.
Due to the applied procedure, smearing occurs and decreases the resolution and the ability to distinguish single reflectors. Therefore, this method is applied only at locations well into the basement where only the Alpine Faultand no 3D valley structures-are expected. Absolute amplitudes are plotted (black) instead of phases.

Parameter Tests
To further enhance weak Alpine Fault reflections, several parameter tests were undertaken during the migration process, as described and discussed in detail by Lay (2021). These tests include calculations with different velocity models and comparison with static corrections applied before prestack depth migration, which did not improve the image quality but confirmed the reliability of the imaged structures. Various approaches and parameters within the migration step were also tested, including the effects of offset and angle restrictions allowing us to emphasize reflectors from a certain direction. Different stacking routines to combine information from single-source gathers were tested as well. All these tests enhanced the images of different parts in the subsurface and were in that sense successful. However, a potential fault reflection itself could not be improved significantly by these different approaches so they were not applied for the following results. Figure 4 summarizes the FVM results of the cube 3D array (VSP2). On the left (Figures 4a-4d), results from the single-source gather at the zero-offset location close to the borehole are shown, with a prominent reflector marked by black arrows. Additional shallower reflectors at z = 0.3 km in Figure 4a are interpreted as top of the basement and are found to correspond to reflection R5 in the seismic section in Figure 3b.

FVM Images
A stack of all 71 migrated single-source gathers is shown on the right of Figures 4e-4h). Absolute amplitudes are stacked which proved to be a robust stacking method emphasizing particularly strong reflectors. Within the stack, more data are included so that strong features such as the top of the basement that are imaged in all single-source gathers are clearly visible and easier to interpret. The different color scale further accentuate these reflectors.
Strong signals from the top of the basement are visible at the western valley flank (gray arrows). During a detailed analysis (Lay, 2021), the reflector in the east (light gray arrow) could be assigned to reflection R7 in Figure 3d. In other words, the migration clearly reveals that the reflections observed to have an opposing slope in the seismic sections originate from the eastern side of the valley.
Although dominant structures causing strong reflections are well imaged in the stack, some details are lost due to the stacking as even neighboring source location images vary in details due to the true 3D nature of the valley and structure studied. In particular, the potential fault-related reflector from Figures 4a-4d is not visible in the stack (Figures 4e-4g). Depending on the image details of interest, both types of data representation will hence be used for the following interpretation.

Resolution
In general, seismic resolution decreases with target depth. This effect is increased for longer wavelengths due to the generally higher velocities at greater depths as well as the lower frequency content of the data at greater depths due to stronger attenuation of high frequencies. This effect is in the depth migrated image is shown in Figure 4a, in which the reflectors at z = 0.5 km are sharper than the strong reflector at z = 1 km.
The migration grid has a node spacing of 5 m in the x-,y-, and z-directions to allow for a detailed seismic image. The seismic data exhibit dominant frequencies of 10-40 Hz with significant energy from the higher frequencies of up to 60 Hz found at near offsets. With P-wave velocities in the range of 1,000-5,000 m/s, the resulting wavelengths λ are between approximately 25 and 200 m. Thus, the resulting spatial resolution is at best 8-50 m using the quarter wavelength criterion. The sweep source type results in a zero-phase Klauder source wavelet. Distinguishing side lobes from the next clear reflection is challenging, especially when the geology is complex. Thus, the spatial resolution is expected to be more in the range of 15-100 m with a better resolution at shallower depths (i.e., containing higher frequencies and having lower P-wave velocities).  Figure 5i. The Gaussian-filtered stack of all FVM cube 3D array data summed along the y-axis (Figure 5a), is considered here to be the most reliable image since it is obtained through the largest volume of input data, that is, the FVM stack of 71 source gathers each with ≥1900 receivers. Also, the summation along the y-axis enhances structures perpendicular to the x-axis such as the potential Alpine Fault related reflector. Although the summation causes smearing and thus lowers the resolution of the identified reflectors, it is used to further enhance the relatively weak single reflectors specifically for the Alpine Fault as discussed in Section 3.5. Figure 5b shows the single FVM-migrated source gather 10148 from the cube 3D array that is consistent with the zero-offset source location shown in Figure 4a. Figure 5c shows the result of Kirchhoff PSDM of the zero-offset data recorded by the hDVS fiber-optic cable. Note that the resulting reflectors are symmetrical around the borehole. A detailed analysis of the KPSDM data set has been described by Kleine (2020). Figure 5d shows an additional FVM image from source gather 2108, recorded by the 2D Aries line 2000 that has been preprocessed by Backes (2020). The additional reflectors (blue lines) form a planar zone parallel to the inferred fault that is at least 550 m in width (perpendicular to the main dip direction) with generally higher reflectivity. This zone could be interpreted as the minimum extent of the outer damage zone around the Alpine Fault at Whataroa further discussed below.  and is observed to be both inherited from previous structures and fault-related Simpson et al., 2020;Williams et al., 2018). Laboratory analyses of seismic P-wave velocities Jeppson & Tobin, 2020a, 2020bLi et al., 2020;Simpson et al., 2020) are not directly transferable to seismic scales due to their small-scale sensitivity (cm-range resolution for laboratory samples instead of tens-of-meter-range for seismic) but imply damage zones as thick as 1 km.

Interpretation of the Alpine Fault Damage Zone
Within the results of seismic imaging presented in Figure 5, several subparallel reflectors with dips of 40°-55° to the south-southeast are observed. Distinct reflectors are imaged from different source locations, that is, illuminations and only one of them (red) is imaged coherently in all images. Single reflections might be caused by faults or fractures illuminated by a certain source-receiver combination. A more detailed understanding of the nature of the reflection and its origin would be possible by analyzing the reflection coefficient as for example, Adam et al. (2020); Stern et al. (2007) have done for the Alpine Fault. However, this is challenging here as reflections observed in the presented data set are weak and discontinuous.
Nevertheless, the location orientation and content of the seismic reflectors are consistent with a zone of extensive fracturing also identified by Massiot et al. (2018), that is, as at least 600 m wide (outer) damage zone. Resolving single structures within the inner damage zone (<160 m after Williams et al. (2018)) interpreted from the observations in DFDP-2B and by fault-zone guided-waves (Eccles et al., 2015) is not possible given the resolution of surface seismic data. However, a zone with denser gouge filled fractures will likely show significant variations in seismic impedance and thus cause reflections. Thus, the total structures represented by the red reflector in Figure 5i are plausibly interpreted to correspond to an inner damage zone.

Comparison of Fault-Related Structures From the 3D Array, 2D Line and the P-Wave Velocity Model
Indicators of the main Alpine Fault from the long Aries 2D-line and the Gaussian-filtered stack of all cube 3D array FVM images (Figure 5e) are used as the basis for the interpretation illustrated in Figure 6. The reflector at z = 0.2 km imaged by the Aries data dips 56° to the SE (black arrow in Figure 6a). The strongest reflector from the Gaussian filtered stack analysis discussed in Figure 5 (red) at a depth of z = 0.7 km has a smaller dip of about 40°. The red dashed lines interpolate between both reflectors and could be interpreted as the principal slip zone of the Alpine Fault. The interpolated fault trace is consistent with both the frontal surface fault trace (Langridge et al., 2018) and the changes in velocity (Lay et al., 2020) included in Figure 6b, both for the top of the basement and the near-surface unconsolidated sediments. Interpolating between additional reflectors, lateral velocity changes and the dextral trace of the fault at the surface might define the extent of the damage zone. If so, the interpreted damage zone would then be nearly 600 m wide and exhibit reflectors with dip angles of 40°-56°. The dip angles of the imaged reflectors are well within the range of observed values from other analyses showing dips of 45°-65° both in the DFDP-2 borehole Townend et al., 2017;Toy et al., 2017) and extrapolated shallow trenches (Langridge et al., 2018).
Additional structures are visible in the summary of Figure 6. Shallow reflectors at z = 0.1-0.3 km between x = 1.8-2.9 km correlate well with the top of the basement extracted from the P-wave velocity model (Lay et al., 2020). The structures seem to form a 100-150 m-deep buried sedimentary basin at x = 2.5 km, which is also seen in gravity results (Jenkins et al., 2020). Whether the central strong reflector at x = 2.4-2.7 km gently dipping to the north-northwest originates solely from lithology or could also be fault-related or a combination of both cannot be determined without further investigations beyond the scope of this analysis (e.g., of the reflection coefficient). This reflector likely represents the top of the basement, which is slightly depressed at this location.  Hall et al. (2016) are projected to the local coordinate system used here (details in Lay (2021)) and plotted as a line drawing. The CMP results are based on a 2D stack and hence less sensitive to 3D reflector origins than the FVM results. The interpreted footwall damage zone might extend to even shallower depths as indicated by the orange dotted box. Other identified potential faults (brown dotted lines) add complexity to the shallow subsurface.

Combined Interpretation and Comparison With Geological and Geophysical Data
A comparison with results from gravity modeling (Jenkins et al., 2020) and geological interpretation  is shown in Figure 7e. The results from gravity modeling by Jenkins et al. (2020) are projected by ∼1 km onto the seismic line from further north. The profile from Toy et al. (2017) is located slightly further south. Nevertheless, both results show similarities to the seismic results. In particular, the depression of the top of the basement at x = 2.4 km, that is, a thicker sedimentary layer is also seen in the gravity measurements and models of Jenkins et al. (2020).
The fault segmentation of the Alpine Fault described by Norris and Cooper (1995) is also thought to cause en-echelon faults in the strike-slip segment. According to Toy et al. (2017), such serial partitioning might explain the second of their two interpreted fault traces (blue dashed line in Figure 7e). Thus, the red area confined from seismic images might alternatively be interpreted as a fault zone exhibiting serial partitioning rather than the damage zone. However, further evidence would be necessary to prove this hypothesis.

Shape of the Whataroa Glacial Valley
The shape of the buried glacial valley derived from the P-wave velocity model of Lay et al. (2020) correlates well with the seismic images, as shown in Figure 8. The migration velocity model used (see Section 3.4) does not include y-direction structures but only x-direction variations. Therefore, the good fit of the reflector to the western valley flank in Figures 8c-8e is not an artifact introduced by the velocity model. In other words, the correlation is a confirmation of results from the velocity model because both methods image the same structure but with somewhat independent means.
Within the valley itself, the isosurface of the P-wave velocity model and the representative seismic images (yz-slices, as in Figure 4g) show the morphology of the buried glacial valley. A typical U-shaped glacial morphology is readily discernible from the seismic images, and exhibits a steep western valley flank. The inferred 3D structure of a valley with relatively steep valley flanks can explain the seismic data very well.
The results are in good agreement with the findings of Jenkins et al. (2020); Thomas (2018) that the Whataroa Valley is an over-deepened glacial trough, as has also been identified seismically in other regions (e.g., Burschil et al., 2019). After the Last Glacial Maximum (27-20 ka), the glacier carving the Whataroa Valley retreated (Barrell, 2011;Suggate, 1990). During the ice retreat, a proglacial lake formed (Thomas, 2018) with rapid sedimentation rates and high lake levels preserving the steep flanks of the bedrock valley. The resulting sediments are interpreted to have been deposited in evolving environments: glacial, lacustrine or marine, and fluvial (Thomas, 2018). The steep valley flanks, particularly in the west, are also seen in gravity models by Jenkins et al. (2020) but might alternatively hint at a subvertical or reverse fault at the western flank. The marked reflector in Figure 8d at the eastern valley side may also be associated with distinct topographical features. Remnants of alluvial terraces were mentioned by Langridge et al. (2018) and are in good agreement with geological maps (Cox & Barrell, 2007). However, we think it is more likely that the higher reflectivity here is caused by heterogeneities or fractures reflecting within the basement.

Indicators of Additional Faults
Within this study, there are weak indicators for reflectors above the top of the basement, that is, indicators of sedimentary structures within the buried glacial valley. However, single-source seismic sections show only weak short reflections above the basement reflection. Thus, it is likely that these incoherent reflections are not constructively interfering during the migration and hence are not clearly visible in the seismic images. In contrast, applying the alternative CMP-processing scheme to the Aries 2D line (Hall et al., 2016) yields more detailed insights into the shallow sedimentary structures. A final Kirchhoff time migration and subsequent depth conversion were included in their processing.
A detailed comparison between CMP post-stack depth migration and FVM pre-stack depth migration results is undertaken by Lay (2021). Although some differences exist, the main features are clearly imaged in both resulting seismic images. The major difference is that the CMP stack shows more details in the shallow (sedimentary) part, which facilitates interpretation at these depths. However, the structural information of the 3D FVM images enables the geometry of structures can be derived, which is a great advantage in a structurally complex environment.
Combining both approaches takes advantage of the strengths of each method.  . This feature is weak in this case and will not be further interpreted. Second, the fault can cause reflections as seen for the Alpine Fault reflector itself (arrow 2, as in Figure 6a) that has been discussed in detail above. Third, a fault is also identifiable by offsets of reflectors with a clear structural discordance in between (arrow 3). We clearly analyze the 3D geometry of this southern structure below.

Internal Sedimentary Structures
Internal sedimentary structures imaged by seismic reflection analysis would give additional information on both, the glacial history of the Whataroa Valley and the shallow fault structures. However, Quaternary sediments form a challenging imaging environment, and previous studies aiming to understand the Alpine Fault (Davey, 2010;Lay et al., 2016;Stern et al., 2007) were unable to robustly resolve sedimentary structures within the Whataroa Valley. Several shallow seismic reflection profiles were analyzed with common-midpoint processing (CMP) approach and were able to identify a few sedimentary structures (Hall et al., 2016;King et al., 2020;Lepine, 2016;Lukács, 2017). Thus, the combination of both imaging approaches (FVM and CMP) provides further information on internal sedimentary structures in Figure 9.
The zoom in Figure 9b shows details of the sedimentary trough and fault 3. The interpreted fault is clearly visible as an offset between reflectors in both the FVM and CMP stack. Within the trough, a set of linear, approximately horizontal reflectors are visible. Several internal sedimentary layers are distinguishable with at least one distinct layer boundary and potentially up to four boundaries (i.e., five seismic sequences). These might be associated with glacio-fluvial or lacustrine sedimentation in the Whataroa Valley (King et al., 2020;Lepine, 2016;Lukács, 2017;Thomas, 2018). The deepest sedimentary package would then be less than 100 m thick. Possibly, these sediments are consistent with the deepest sediments described by Thomas (2018) that were interpreted to result from a depositional environment of ice contact to an ice distal lacustrine setting. Also, the bedrock of the DFDP-2B borehole  noted a higher amount of gravel close to the bedrock contact.

3D Analysis of the Southern Fault Three
On this xz-profile, the seismic unconformity forming the structural discordance is subvertical. A 3D plot of the phase consistent FVM stack for the cube 3D array is shown in Figure 10 with a view from the north into the valley clarifying the fault geometry. A strong reflector is visible on the presented depth slice (Figure 10b) with a change in seismic characteristics just north of where the arrow points to the potential fault trace. Identifying fault three in the yz-profile is more difficult due to the subvertical dips, but the most likely associated feature is marked in Figure 10d. Additional analysis of fault three in-depth slices is presented in Lay (2021) and illustrates the east-west strike.
The isosurface at v p = 2,600 m/s from the P-wave velocity model (Lay et al., 2020) correlates well with the structures identified in the seismic FVM image (Figures 10c and 10e). A central trough of thicker sediment layers is visible between x = 3-4 km. That trough is bordered by a steep western valley flank (right in Figure 10). Again, the good correlation between reflectors in the seismic image and the results of the P-wave velocity model is obvious. Moreover, this trough correlates with a low gravity anomaly (4 mGal) identified by Jenkins et al. (2020) within the Whataora Valley south of the DFDP-2 drill site. They interpreted this structure as a thick sedimentary layer infilling an over-deepened glacial valley.
Fault 3 might be associated with massive fracturing  in the hydrogeological active hanging wall (Janku-Capova et al., 2018;Menzies et al., 2016;Townend et al., 2017) but might also be associated with local tectonic stress. Upton et al. (2018) show that topography controls the local stress field and also affects fault  structures. Hence, alterations in the local stress field might cause deeper faults potentially loosely linked to the Alpine Fault.

Anisotropy
Seismic anisotropy is expected from the analysis of elastic properties of rock samples collected near the Alpine Fault Christensen & Okaya, 2007;Godfrey et al., 2000;Graham, 2020;Jeppson & Tobin, 2020b;Li et al., 2020;Simpson et al., 2020). This anisotropy has a great influence on seismic imaging as proven, for example, by modeling studies at the Alpine Fault conducted by Adam et al. (2020). Simon et al. (2019) successfully applied anisotropic Kirchhoff PSDM to data collected at a borehole in Sweden. However, no appropriate anisotropic migration velocity model could be derived for the Alpine Fault given the generally low signal-to-noise ratio. Transferring laboratory measurements to field results is challenging as shown by Jeppson and Tobin (2020b). Even determining anisotropy in the shallow region of the Alpine Fault is very challenging according to Godfrey et al. (2000). Zhang and Schmitt (2020) have demonstrated the anisotropic nature of the heavily fractured rock mass along the DFDP-2B borehole, but a complete understanding of this anisotropy is not yet available. This known anisotropy is not included in the imaging described here and as such will affect the locations and inferred dip of the observed reflectors (e.g., Isaac & Lawton, 1999;Simon et al., 2019). Anisotropy might thus explain some mismatch as seen between different parts of the data set when comparing single-source gathers and borehole recordings. In the future, an anisotropic velocity model might be able to further improve the match of various reflectors. The source of the reflectivity, too, remains enigmatic and may await future drilling campaigns for a fuller appreciation.

Conclusions
An extended seismic survey including a dense 3D array and VSP records at the DFDP-2B drill site in the Whataroa Valley has been analyzed in detail. Two main goals were pursued to characterize the drill site: (a) identifying glacial valley structures and (b) illuminating Alpine Fault-related structures at depth. The applied FVM technique is able to image and distinguish both the glacial valley and fault-related structures despite the challenging inhomogeneous near-surface environment.
An interpreted Alpine Fault reflector is imaged at a depth of ∼0.8 km dipping 40° to the southeast and about 100 m below the bottom of the DFDP-2B borehole. This reflector lies in the center of a zone containing several sub-parallel reflectors dipping ∼40°-60° to the southwest at depths between 0.2 and 1.2 km. This reflecting zone extends ∼550 m perpendicular to the reflectors and might be interpreted as the minimum extent of the outer damage zone. Overall, the potential hanging wall damage zone is also seen at shallower depths of as much z = 0.2 km, that is, about 0.3 km below the surface. A steepening appears to occur at these shallower depths to at least ∼56° at a depth between 0.2 and 0.5 km.
Additionally, a major east-west-trending near-vertical fault is identified about 1 km south of the drill site that appears to coincide with topographical structures. This fault might be explained by shear-related processes that are locally influenced by the valley structures, but needs to be analyzed in more detail.
In addition to the fault-related structures, the glacial Whataroa Valley itself has been revealed by the seismic data to have steep valley flanks, particularly in the west. This confirms the results of previous geological and geophysical studies and provides additional information as the analyzed data set is a rare example of a densely spaced 3D seismic array in this area. Earlier assumptions of side reflections from the valley flanks are now confirmed by seismic data.
About 0.5 km SSE of the drill site, a 350 m-thick trough is identified that is likely to be filled by sediments, which has also been detected by other geophysical methods such as gravity. Some indicators hint at distinguishable internal sedimentary layers with at least one distinct layer boundary and potentially up to four boundaries (i.e., five seismic sequences). Thus, the seismically imaged trough might reveal details of sedimentation related to the glaciation and postglaciation history of the Whataroa Valley.
FVM images have been required to directly image steeply dipping faults. However, the complementary use of CMP-based processing provides additional information about the shallow subsurface that we interpret in terms of glacial sedimentation sequences in the over-deepened glacial deep point south of the borehole.
Overall, the 3D seismic images provide new insight into tectonic and glacial features at the DFDP-2B drill site in the Whataroa Valley. Thus, comprehensive 3D seismic data are clearly needed to understand complex geological environments.
However, the newly imaged structures also show that the fault zone of the Alpine Fault is complex and has still not revealed all its details. Further analysis of seismic anisotropy and borehole data might be able to illuminate further aspects of the fault zone. lington, University of Otago, University of Auckland, University of Alberta, University of Calgary, Schlumberger and TU Bergakademie Freiberg). DRS was supported by NSERC discovery and Canada Research Chairs Program, the Sercel System was provided through grants from the Canadian Foundation for Innovation and the Alberta Government. Supplementary field equipment was provided by the Geophysical Instrument Pool Potsdam. The processing of the seismic data presented was performed with Seismic Unix (Cohen & Stockwell, 2002) and Promax. TU Bergakademie Freiberg acknowledges a software grant from Haliburton Landmark for this. Perceptually uniform color maps from Kovesi (2015) are used to avoid interpretation artifacts introduced by the color scale. Additionally, we used GMT (Wessel et al., 2013) for data handling and plotting. A full description of the field procedures to obtain the presented data is provided in the report by Townend et al. (2016). They used seismic data are made available in an archive by . We thank Associate Editor Gail Christeson and two anonymous reviewers who helped to improve the manuscript.