Generating High‐Fidelity Reflection Images Directly From Full‐Waveform Inversion: Hikurangi Subduction Zone Case Study

Full‐waveform inversion (FWI) can resolve subsurface physical properties to high resolutions, yet high‐performance computing resources have only recently made it practical to invert for high frequencies. A benefit of high‐frequency FWI is that recovered velocity models can be differentiated in space to produce high‐quality depth images (FWI images) of a comparable resolution to conventional reflection images.

• Full-waveform inversion (FWI) using high-frequency data produces high-resolution velocity models that can be directly utilized to produce high-fidelity depth images • FWI images are comparable to conventional reflection images and provide high-resolution physical property models to assist interpretation • As FWI is applied to raw data, there are few subjective processing steps, making images easily reproducible for time-lapse imaging 2 of 10 A common practice in the petroleum sector is to obtain FWI velocity models using relatively low-frequency inversions (<10 Hz) and using these to migrate processed seismic reflection data, which leads to an improvement in the coherence of the reflection image (e.g., Routh et al., 2016;Sirgue et al., 2010;Warner et al., 2013). In the past decade, high-performance computing resources have become more readily available in both industry and academia, making it feasible to increase the maximum inversion frequency, for example, 18 Hz in Arnulf et al. (2021), 31 Hz in Mispel et al. (2019), 40 Hz in Routh et al. (2017), and 100 Hz in Kalinicheva et al. (2020). Consequently, it is now possible to resolve the physical properties of the shallow crust on length scales of tens of meters. These greater resolutions mean that high-frequency-FWI can now be used to resolve the physical property contrasts responsible for generating reflections without performing conventional seismic processing, velocity-model building, and migration imaging. As a result, FWI can be used to generate accurate depth images directly (e.g., Kalinicheva et al., 2020;Routh et al., 2017;Sedova et al., 2019;Zhang et al., 2020) using a simple integrated workflow and minimal pre-processing of the input data.
Here, we present the methodology, and demonstrate the utility, of FWI for high-resolution subsurface imaging at the Hikurangi subduction margin, offshore New Zealand's North Island ( Figure 1). We applied high-frequency FWI to a legacy 2D data set acquired in 2005 and a near-coincident 2D profile, inline M1220, selected from the recently acquired 3D NZ3D data set. Both lines cross several International Ocean Discovery Program (IODP) drill holes. We use a multi-stage workflow that starts with building the macro velocity model and then moves onto resolving reflection horizons through appropriate selection of frequency and offset ranges. Finally, we differentiate our FWI velocity models vertically in space to obtain depth images of the subsurface. The resulting velocity models and depth images resolve ambiguities in previous seismic imaging, enabling more informed interpretations of the stratigraphy surrounding an IODP drill hole, located on the frontal accretionary wedge.

Case Study: Hikurangi Subduction Zone
The shallow subduction interface (<15 km depth) at the Hikurangi subduction margin, offshore Gisborne, New Zealand, demonstrates a wide range of seismic behaviors, from slow-slip events (SSE) through to tsunami earthquakes (Figures 1a and 1b) (Wallace & Beavan, 2010). The enigmatic seismic behaviors of this plate interface, combined with its shallow depth, has made it the focus of many geophysical studies. We have chosen to attempt high-frequency FWI along the well-studied 05CM-04 seismic profile (Figure 1b), as there are existing seismic pre-stack depth migrated (PSDM) images ( Figure 1c) (Barker et al., 2018) that we use for comparison, and a previous application of low-frequency FWI to these data was successful (Gray et al., 2019). Much of this profile has also been re-surveyed within the modern NZ3D seismic experiment (Figure 1b, Han et al., 2021), which allows us to test the performance of high-frequency FWI on both a legacy and modern data set. Additionally, the 05CM-04 profile passes through three IODP drill sites , allowing us to ground-truth our models and constrain stratigraphic interpretations.
Seismic imaging in this region is complicated by extreme topography and subsurface structure, anisotropy, and attenuation ( Figure 1c, Arai et al., 2020;Gray et al., 2019). There are regions where seismic horizons are difficult to track in the previous conventional seismic images, leading to ambiguity in their interpretation (Barker et al., 2018;Barnes et al., 2020). These challenging conditions, combined with the high level of geophysical interest, provide an excellent case study to demonstrate the FWI method and its imaging capabilities.
Here, we focus on imaging the shallow Pāpaku thrust fault (Figure 1c) within the accretionary prism, which was drilled during IODP expedition 372 and 375 . We compare the results of earlier reflection imaging from the 05CM seismic experiment with FWI models and imaging derived from the same data and FWI results obtained from the modern NZ3D seismic experiment.

05CM Seismic Experiment
05CM-04 is a legacy seismic reflection data set ( Figure 1b) acquired in 2005 (Barker et al., 2009). Seismic acquisition parameters included a 12 km streamer with 960 channels spaced 12.5 m apart, recording shots fired every 37.5 m from a 4,140 cubic inch airgun array and sampled at 2 ms (Barker et al., 2009).
Conventional seismic reflection processing of this data set involved band-pass filtering, gain recovery, swellnoise attenuation, seafloor reflection-based multiple elimination, signature deconvolution, radon multiple removal, and CDP binning at 12.5 m (Barker et al., 2018). Kirchhoff PSDM was performed iteratively using velocities picked from reflection semblance spectra, with updated P-wave velocity models used to improve subsequent images (Figure 1c).

NZ3D Seismic Experiment
The NZ3D experiment is a modern three-dimensional seismic reflection and refraction data set acquired in 2018 ( Figure 1b, Han et al., 2021). Acquisition parameters included four 6-km-long hydrophone streamers spaced 150 m apart, each with 468 channels spaced 12.5 m apart, recording shots from two 3,300-cubic-inch airgun arrays alternately fired every 25 m and sampled at 2 ms. Fifty-eight seismic lines were acquired, spaced 300 m apart on average. Ninety-seven ocean bottom seismometers (OBS), spaced 2 km apart, also recorded the seismic sources. Refracted arrivals recorded on the OBS array were used to perform 3D travel-time tomography, producing a 3D P-wave velocity model, which we have used for our starting model (Arai et al., 2020). Here, we invert inline M1220 of the NZ3D data set, which is coincident with the 05CM-04 profile.

IODP Expeditions 372 and 375
Site U1518 was drilled, logged, and cored during IODP Expedition 372 and 375 in 2018 ( Figure 1b, Saffer et al., 2019;Wallace et al., 2019). This site is located on the frontal accretionary wedge of the lower continental slope ( Figure 1c). The target of this drill site was the Pāpaku fault, a thrust fault that extends from the plate interface to the seafloor. It is believed that this fault system hosts significant inter-plate motion and could be one of the source locations for SSEs (Fagereng et al., 2019). Drill coring reached a depth of 494 m below seafloor (mbsf), successfully sampling the hanging wall, fault zone (323 mbsf), and footwall. Seismic data can be used to track lithologies and stratigraphic horizons away from the drill site and can aid the understanding of deformation and rupture processes. Interpretations are, however, reliant on the quality of the imaging.

Methodology
Our multi-stage workflow has three steps. In step 1, we perform low-frequency FWI and include long offset data to obtain a model of the macro velocity structure. We use this as the starting model for step 2 so that structures recovered will be positioned at the correct depth. In step 2, we perform high-frequency FWI on short-offset data to recover changes in P-wave velocity over shorter wavelengths. In step 3, we differentiate this model to obtain a reflectivity image. Further rationale on our approach can be found in the Supporting Information S1.
We tested various inversion strategies and applied identical schemes to the hydrophone streamer data from inline M1220 of the NZ3D data set and the 05CM-04 data set.
Before running FWI, we produced accurate source wavelets and starting models, the details of which can be found in the Supporting Information S1.

Step 1: Low-Frequency Inversion Scheme
We first checked that data predicted using our starting model ( Figure 2a) were not cycle skipped; they reproduce the observed wavefield to within half a seismic wavelength at low frequencies (<6.0 Hz) (Figures S5a and S6a). We downsampled the starting model to 12.5 m node spacings to decrease the compute time. Then, we performed 2D, time-domain, acoustic, isotropic FWI using the codes of Warner et al. (2013). Water density was assumed to be 1,000 kg m −3 , while subseafloor densities were determined using Gardner's relationship (Gardner et al., 1974) between P-wave velocity and density for velocities greater than 1,560 m s −1 . FWI was performed using the following inversion scheme: 1. Observed data were first limited to offset distances of 3,200 m, the offset to which the synthetic wavefield matched the observed wavefield on all shot gathers. 2. Data were time-windowed by 3 s after the first arrival. This windowing preserved refracted arrivals, which provide a firm control on velocity while removing later reflections. 3. FWI was then performed with 40 iterations using frequencies ≤3.0 Hz, followed by a further 40 iterations with frequencies ≤3.4 Hz. 4. (1-3) were repeated with increasing data offsets of 4,000, 5,000, and 6,000 m, with the output model of the previous inversion being used as the starting model for the new inversion. Refracted data at longer offsets have sampled greater depths in the subsurface, so this incremental layer-stripping scheme builds deeper velocity structure with each offset increment. 5. At 6,000 m of data offset (NZ3D hydrophone streamer length), we inverted with the upper frequency in each band increasing from 3.0 to 6.0 Hz (3.0, 3.4, 3.9, 4.5, 5.2, and 6.0 Hz) iterating 20 times per step.
Such an inversion scheme builds the velocity structure from the top down, which means that data misfits at the longer offsets (>3,200 m) were corrected as the shot-receiver offset was incrementally increased. The resulting low-frequency velocity model of the M1220 profile can be seen in Figure 2b (complete profile models: Figures S3b-S4b).
To check the validity of the models, we compared synthetic data with the observed data ( Figures S5 and S6) at multiple stages during the inversion. As the inversion progressed, the modeled and observed wavefield match improved at all offsets, indicating that the velocity model was moving in the right direction. We also checked how our velocity models compared with the existing PSDM images of Barker et al. (2018) and observed a strong correlation in the long-wavelength structures ( Figures S7 and S8). These quality assurances suggest that our final models accurately represent the subsurface velocity structure down to at least 500-600 m below the seabed, the maximum penetration of the refracted wavefield ( Figure S9).

Step 2: High-Frequency Inversion Scheme
The output of the low-frequency inversion scheme is a macro velocity model of the subsurface, which we use as the starting model for this high-frequency inversion scheme. The basic parameterization of the FWI remained consistent. The high-frequency inversion scheme was as follows: 1. Low-frequency models were resampled to 6.25 m node spacings, allowing FWI up to 64 Hz. 2. Observed data were limited to offsets <3,200 m to include only sub-vertical reflections and short offset refractions, as elastic and anisotropic effects in longer offset data caused smearing in the final reflection image. 3. Data were time windowed by 6 s after the first arrival to include deeper reflections. 4. FWI was performed, starting with frequencies up to 3.0 Hz and increasing in stages to a maximum frequency of 38 Hz. The maximum inversion frequency was increased by ∼15% for each new frequency band (e.g., 3.0, 3.4, 3.9 Hz, etc…). For each frequency band, the inversion was iterated 20 times.
The resulting velocity model of the M1220 profile, in the location of IODP site U1518, can be seen in Figure 2c (complete profile models: Figures S3c and S4c). These models present structures with much shorter wavelengths when compared to the low-frequency inversion result (Figure 2b). We repeated the same quality checks that we employed for the low-frequency models (Figures S10-S13) and generated common image gathers (CIGs) ( Figure S14).

Step 3: Generating Reflection Images
Reflections are generated by contrasts in acoustic impedance in the Earth's subsurface, dictated by changes in rock seismic velocity and density at depth. We performed a vertical differentiation along each trace of the FWI velocity models to highlight rapid variations in P-wave velocity that generate reflections. This differentiation of the velocity field represents an approximation of the subsurface impedance contrast (Huang et al., 2021). Figure 2D-2F shows the results of the vertical differentiation applied to the starting model, low-frequency, and high-frequency inversion models, respectively, of the NZ3D data set. In the smooth starting model (Figure 2d), the seafloor is the only interface with a velocity contrast capable of generating a reflection. Long-wavelength structures become observable within the differentiation of the low-frequency model (Figure 2e), including folded sedimentary layering, faulting, and a bottom-simulating reflector (BSR). Differentiation of the high-frequency FWI velocity model (Figure 2f) reveals an increase in the number of reflection horizons and significant sharpening of subsurface features, highlighting the importance of high frequencies in FWI imaging (FWII).

Results and Discussion
We now have three independent depth images of the subsurface in the region of IODP Site U1518 (Figures 3a-3c) (complete profiles: Figure S15). The first is the conventionally processed 05CM-04 PSDM image of Barker et al. (2018) (Figure 3a); the second is the FWI image derived directly from unprocessed raw field data using the legacy 05CM-04 data set ( Figure 3b); the third is the FWI image derived using inline M1220 of the modern raw 3D NZ3D seismic data set ( Figure 3c). The broad match between the structures resolved by the PSDM and FWI imaging method is remarkable. To first order, the Pāpaku fault and its primary splay faults and the principal sedimentary reflections are all imaged by both methods and datasets. Furthermore, the FWI images shown in Figures 3b and 3c demonstrate that high-frequency FWI can be used to produce high-quality images of reflectivity, from raw field data, with a comparable bandwidth to a PSDM (Figure 3a). Note that these FWI images are not made by explicitly migrating the reflection data using the FWI velocity model, they are extracted directly from the FWI velocity model itself.
The three images are similar, but the quality of the image and the ease in tracking particular reflection horizons is variable across the accretionary prism. For example, the Pāpuka fault is clear in all three images but is best imaged on the PSDM below ∼3.5 km depth and at shallower depths in the NZ3D FWII (Figures 3a-3f).
A suite of bright, subparallel, antiformal reflections is clearly imaged in the PSDM in the western half of Figure 3a but becomes challenging to track eastwards toward U1518 and the Pāpaku fault (Barker et al., 2018;Barnes et al., 2020). In our FWI velocity models and images, reflections in the hanging wall and footwall of the Pāpaku Fault are clearer (arrows i-v in Figures 3b and 3c). We see steeply dipping reflections that offset other reflections (e.g., iv in Figures 3b and 3c), leading us to believe they are likely to be splay thrust faults (Figure 4). Consequently, when we use the FWI images to track horizons away from the drill hole (Figures 3e and 3f), we interpret them differently from previous studies, such as the one shown in Figure 3d ( Barnes et al., 2020).
Our interpretation of folds and thrust faults close to the drill hole in the FWI images ( Figure 4b) results in the base of SU4 being significantly below the base of the hanging wall at U1518. This interpretation contrasts with that of Barnes et al. (2020), in which SU4 intercepts the fault at the location of U1518. It appears that our new interpretation is more consistent with biostratigraphic data obtained from a drill core at IODP Site U1518 (Figure 4f) (base SU4 has been dated at ca. 0.78 Ma at drill site U1520, and the sediments above the fault zone at U1518 are 0.6 Ma) (Crundwell & Woodhouse, 2021). In the biostratigraphy, there are two age inversions, one in the hanging wall of the Pāpaku fault (∼185 mbsf) and one across the Pāpaku fault (323 mbsf) (Figure 4f). Such age inversions are typically the result of reverse faulting juxtaposing older lithologies above younger ones. One of the reflections that we interpret as a splay thrust in the FWI image ( Figure 4b) crosses IODP Site U1518 at the same depth as the shallower age inversion.
Another difference between the imaging methods is the imaging of a bottom simulating reflector (BSR). The reflection is located ∼700 mbsf in the FWI images (i in Figures 3b and 3c) but is poorly imaged in the PSDM. The reflection amplitude of the BSR is relatively weak through the deformation front, suggesting it is probably associated with only minor changes in velocity and density. We assume that this event has been impaired by the PSDM processing, perhaps resulting from the poorer velocity model used for the migration. We note that this BSR is resolved in post-stack time migrated (PSTM) imaging of the NZ3D data set (Han et al., 2021).
An issue with comparing velocity models is that velocities obtained using isotropic FWI principally recover horizontal velocity, whereas near-normal incidence reflections recover moveout velocities, and sonic well logs measure vertical velocity (e.g., Christeson et al., 2018;Thomsen, 1986). In Figure 4c, we compare the starting (red) and FWI velocity (blue) models with velocities obtained from wireline logging (black) at Site U1518. FWI has significantly improved the fit with the down-hole velocity log compared to the starting model between 60 and 325 mbsf and shows the velocity decrease across the Pāpaku fault. When running the inversions, we found that recovered velocities close to the seabed changed depending on how we set the parameters that controlled the transition from water to rock density, and this may be responsible for the poor fit in the upper 60 mbsf. Below about 325 mbsf, the FWI velocity is between 50 and 100 m s −1 higher than the well log velocity, which can be explained by anisotropy since layered silts and clays usually are faster in the subhorizontal (bedding) direction.
A significant advantage of using FWI to generate reflectivity images is its simplicity. Many operations are performed separately in a conventional PSDM workflow, such as denoise, deghost, demultiple, zero-phase the data, building a velocity model, removing refractions, and migrating reflections to generate a final stacked image. In contrast, using FWI, these operations are all implicit within a single integrated algorithm. Another possible benefit of the FWI imaging method is its reproducibility, given that we recover similar results from two seismic datasets collected 15 years apart. These results support the use of FWI for time-lapse imaging and interpretation, for example, at carbon capture storage sites (Mispel et al., 2019). Variations in the velocity models over time could infer temporal changes in the presence of fluids or fluid pressure changes accompanying large earthquakes along the Hikurangi subducting margin.

Conclusions
Here, we have presented the application of FWI imaging as an effective alternative to established conventional reflection imaging methods. FWI imaging can improve the resolved structure of the shallow subsurface. At the same time, the production of high-resolution physical property models and depth images provides the user with multiple interpretative tools. By applying the high-frequency FWI method to two generations of seismic data, we hope to motivate researchers to reprocess existing seismic datasets; encourage future seismic experiments to be optimized for FWI; and demonstrate the reproducibility of FWI results for time-lapse imaging.