Seismic Velocity Structure Beneath the Tofino Forearc Basin Using Full Waveform Inversion

Abstract

Given the effects of steep dips and large lateral variations in seismic velocity beneath the Vancouver Island continental shelf, seismic processing and travel time inversion are inadequate to obtain a detailed velocity model of the subsurface. Therefore, seismic full waveform inversion is applied to multichannel seismic reflection data to obtain a high-resolution velocity model beneath the Tofino fore-arc basin under the continental shelf off Vancouver Island margin. Seismic velocities obtained in this study help in understanding the shallow shelf sediment structures, as well as the deeper structures associated with accreted terranes, such as Pacific Rim and Crescent terranes. Shallow high velocities, as large as ∼5 km/s, were modeled in the mid-shelf region at ∼1.5–2.0 km depth. These coincide with an anticlinal structure in the seismic data, and possibly indicate the shallowest occurrence of the volcanic Crescent terrane. In general, seismic velocities increase landward, indicating sediment over-consolidation related to the compressional regime associated with the ongoing subduction of the Juan de Fuca plate and the emplacement of Pacific Rim and Crescent accreted terranes. Seismic velocities show a sharp increase about 10 km west of Vancouver Island, possibly indicating an underlying transition to the Pacific Rim terrane. A prominent low velocity zone extending over 10 km is observed in the velocity model at 800–900 m below the seafloor. This possibly indicates the presence of a high porosity layer associated with lithology changes. Alternatively, this may indicate fluid over-pressure or over-pressured gas in this potential hydrocarbon environment.

1. Introduction

The Juan de Fuca plate is subducting beneath the North American plate at 46 mm/year in the northern Cascadia margin. As a result of subduction, sediments have been scraped off from the downgoing oceanic plate to form the accretionary wedge [1]. On southernmost Vancouver Island, three major northeast dipping faults separate the crust in to two narrow zones: the ocean basaltic Crecent terrane and the marine sedimentary Pacific Rim terrane (Figure 1). The origin of these terranes is still debated, but it is believed that these terranes were accreted to the margin in late Eocene time, between 55 and 42 Ma ago [1,2,3]. The Pacific Rim terrane is separated from the Wrangellia terrane by the landward dipping West coast fault. The Crescent terrane is separated from the Pacific Rim terrane by the landward dipping Tofino fault. The landward dipping Crescent thrust bounds the Crescent terrane to the west. The Tofino forearc basin lies above the accretionary wedge, the Crescent terrane and the Pacific Rim terrane (Figure 1). It has a width of ∼60 km and a length of ∼200 km, and it contains up to 4 km of sedimentary rocks [1]. Shell drilled six exploratory wells in the Tofino basin in the 1960s, which show up to 3 km of Miocene to Recent mudstones and minor sandstones [4]. Some of these wells also indicate Eocene sediments and volcanics [4,5]. The volcanics are the source of the Prometheus magnetic anomaly [6], with a magnetic anomaly high over the Crescent terrane and a low over the Pacific Rim terrane [1]. Beneath the eastern portion of the shelf, seismic reflection data indicate a localized thickening of the Tofino basin sediment section to more than 3 km, and this was interpreted as a fossil trench [1]. As a consequence of the great sediment thickness, a gravity anomaly low is observed over the fossil trench region [7].

Subduction of the Juan de Fuca plate beneath the North American plate has been the main tectonic process controlling the growth and deformation of basin structures. The Tertiary strata beneath the continental shelf are relatively undeformed and form a wedge which thickens seaward from the west coast of Vancouver Island to the outer shelf. The observed faults, folds and unconformities in the basin indicate that the deformation increases southward from Brooks peninsula at the north end of Vancouver Island [8]. Anticlines and thrust faults which occur beneath the shelf edge are the result of intracrustal compression associated with the subduction process along the margin [9]. Hyndman et al. [1] argue that these structures result from high pore pressures arising from horizontal tectonic shortening and fluid expulsion.

A series of multi-channel seismic (MCS) lines were collected across Tofino basin in 1985 and 1989 by the Geological Survey of Canada (Figure 1). Interpretation of seismic reflection profiles, traveltime inversion velocity models and biostratigraphic well data suggest that character of the basement has largely controlled deformation of the overlying Tofino basin sediments [5,11,12,13,14,15]. On migrated seismic sections, some thrust faults appear to penetrate close to the top of the subducting oceanic crust [6,16]. The amount of displacement along the faults and style of deformation varies along the margin, which indicates localized variations of pore fluid pressure. Two-dimensional traveltime tomography of MCS first-arrival times suggest that sediment deposition increased more rapidly in the later half of the Tofino basin history [10,12]. According to this study, the development of the basin is related to the reactivation of the terrane-bounding Tofino fault. P-wave velocities derived from tomography show that anticlinal folds in the accretionary wedge sediments indicate low velocities at the apex of the fold. In contrast, folded sediments over the Crescent Terrane exhibit a velocity pattern that more closely mimics the shape of the folds [12]. The lower velocities at the apex of the fold are interpreted as the result of fracturing of older (Late Pliocene), more lithified sediments that contain fluids derived from deeper sediments and the accretionary wedge as opposed to overlying younger (Pleistocene to Holocene) sediments where the low-velocity is due to their lower burial depth.

Recent studies showed that seismic full waveform inversion (FWI) provides more detailed information of the subsurface structure [15,17,18,19,20]. Kamei et al. [17] applied the frequency-domain acoustic FWI to long offset marine ocean bottom seismic (OBS) data to image the megasplay fault system of the seismogenic Nankai subduction zone off Japan. The velocity model from this study clearly shows a low-velocity zone extending laterally over 40 km associated with fluids driven from the downgoing subducting plate, and clearly delineates the accretionary prism from the downgoing Phillipine Sea Plate. In a different study, Smithyman et al. [19] applied 2.5-D frequency-domain visco-acoustic FWI to land vibroseis data collected along crooked roads to obtain velocity and attenuation models of the Nechako-Chilcotin plateau of British Columbia, Canada. Models from this study provide important information of the Cenozoic volcanic structures down to 3 km depth.

Here, the 2D frequency domain acoustic full waveform inversion of Pratt [21] is applied to one of the marine MCS reflection lines on the Vancouver Island continental shelf (89-06, Figure 1). A time migrated section of this line is shown in Figure 2. The processing details of this line are presented in Yelisetti [15]. Yuan et al. [22] applied 1D full waveform inversion to a small number of individual common-reflection-point gathers from several of the 1989 Geological Survey of Canada MCS lines located in the mid-slope region of the margin. Our study is the first 2D waveform tomography application to seismic data on this margin. We show that this technique can resolve significant structural details in MCS data that are not achievable with previously applied methods including 2D traveltime tomography of first arrivals. We also present the imaging challenges encountered while inverting structurally complex, noisy, short offset data devoid of low frequencies. However, the primary objective of this study is to quantitatively image the structure of the upper part of the Tofino fore-arc basin under the continental shelf to further our understanding of the margin. In particular, we aim (1) to obtain the detailed P-wave velocity structure of the Tofino basin sediments, (2) to understand the nature of deformation within Tofino basin sediments, and (3) to understand the relationship between Tofino basin sediments and the underlying Crescent terrane and Pacific Rim terrane.

2. Theory

The first step in the application of full waveform inversion is the determination of a very good starting velocity model. We obtained this model by traveltime tomography of first arrivals on individual shot gathers, similar to the method applied by Hayward and Calvert [12]. However, we utilized the tomographic techniques of Zelt and Smith [23] and Zelt and Barton [24], rather than that of Aldridge and Oldenburg [25] applied by Hayward and Calvert [12]. Thus, we first present an overview of the traveltime inversion methodology, and subsequently discuss some of the key theoretical considerations in the full waveform inversion method of Pratt [21].

2.1. Traveltime Tomography

In traveltime tomography, forward modeling and inversion are used to iteratively update the velocity model. Inversion tries to minimize the misfit between the measured times from forward modeling and observed times. However, the solution of an inverse problem is highly non-linear and is non-unique as the data contain errors and the problem is under-determined. This is controlled by regularizing the solution, which helps to stabilize the ill-conditioned and singular valued problem by using prior information. Prior information can reflect the types of solution desired, or can consist of estimates of model parameters or their derivatives. The final model is obtained by minimizing the objective misfit function ($\mathsf{\Phi }$) with respect to the model (m) as follows [24]:

$$\mathsf{\Phi }(m)=\delta t^T{C}_d^{-1}\delta t+\lambda [m^T{C}_h^{-1}m+S_z{m}^T{C}_v^{-1}m],$$

(1)

where $\delta t$ is data residual vector, $C_v$, $C_h$, and $C_d$ are vertical and horizontal roughening and data covariance matrices, respectively; $S_z$ controls the horizontal versus vertical smoothness. The trade-off parameter ($\lambda$) provides the model with the least structure, and also determines the relative importance of data and prior information. $\lambda$ is selected based on expected misfit for Gaussian errors to give misfit to data of ${\chi }^2=N$. This applies the prior information subject to fitting data to a statistically appropriate level.

2.2. Full Waveform Inversion

Waveform tomography represents a more advanced approach, in which the whole waveform is used in the reconstructions. In general, the approach of waveform tomography represents an attempt to fully account for the complex interaction of the propagating wave and the target.

In the present study, traveltime tomography uses only first arrivals, whereas waveform tomography uses reflections, refractions, and diffraction events. Therefore, traveltime tomography is fundamentally limited in resolution and waveform tomography restores the missing resolution of traveltime tomography [26,27]. The resolution is of the order of first Fresnel zone width in traveltime inversion, whereas it is of the order of wave length in waveform inversion. For a velocity of 2000 m/s and a frequency of 10 Hz, the wavelength is 200 m, and the first Fresnel zone width is about 866 m considering a maximum propagating distance of 3758 m in this study. The high resolution comes from the fact that it uses not only the first arrivals but also secondary arrivals.

This study uses the 2-D frequency domain acoustic waveform inversion approach of Pratt [21]. In the forward modeling, synthetic shot gathers are estimated using the following constant density acoustic wave equation in the frequency domain.

$${\nabla }^{2}u(x,\omega )+\frac{{\omega }^2}{c^2(x)}u(x,\omega )=-f(x,\omega ),$$

(2)

where u is the wavefield vector (pressure or displacement), and f is the source vector. $\omega$ is the angular frequency, and c is the wave velocity, which can be complex valued if the medium is attenuating.

The above acoustic wave equation is solved using finite difference method in frequency domain. The equation can be written in the matrix form as shown below.

$$Su=f,$$

(3)

or

$$u=S^{-1}f,$$

(4)

where S is an impedance matrix or differencing matrix.

The non-linear Born scattering equation which describes the interaction of the source wavefield with the scattering medium while traveling to the receiver is given by Wu and Toksöz [28]:

$$U(r_s,r_g)=-{\omega }^2\underset{V}{\int }\delta s(r)G(r,r_s)G(r,r_g)dr,$$

(5)

where $U(r_s,r_g)$ is scattered wavefield, $\omega$ is angular frequency, $\delta s$ is scattering potential, and G is free space Green’s function.

The residual wavefield ($u_{res}$) is estimated by substracting the observed wavefield (d) from the forward modeled wavefield (u).

$$u_{res}=u-d.$$

(6)

The data misfit (E) is defined as the $L_2$ norm of the data residuals.

$$E=\frac{1}{2}{u}_{res}^T{u}_{res},$$

(7)

where $u_{res}^T$ indicates the conjugate transpose of $u_{res}$. The inversion modeling minimizes the misfit between observed waveforms and forward modeled waveforms by iteratively updating the model (m) in the frequency domain. The model perturbation ($\delta m$) is related to the data misfit by the following equation (Newton method):

$$\delta m=-H^{-1}\nabla E,$$

(8)

where H is Hessian matrix which contains the second derivatives of the misfit function. The misfit function is given by the sum of squared data residual errors. The computation of the gradient of the misfit function does not involve any partial derivate calculations. It is computed at each iteration by multiplying the forward propagated wavefield with the backpropagated data residuals [29]. However, the Hessian matrix is very huge and often singular. Hence, it is very difficult to invert. The conjugate gradient method is one of the most effective methods for the iterative inversion. It ignores the Hessian and uses the gradient of the misfit function ($\nabla E$) directly to update the model as shown below.

$$m^{(k)}=m^{(k-1)}-{\alpha }^{(k-1)}\nabla E^{(k-1)},$$

(9)

where k is the number of iterations, and $\alpha$ is the step length which is given by

$$\alpha =\frac{|\nabla E_m{|}^2}{|J{\nabla }_m{E|}^2},$$

(10)

where $||$ represents the Euclidean length of the vectors, and the Fréchet derivative matrix or the sensitivity matrix J is given by

$$J=\frac{\partial u}{\partial m}=-S^{-1}\frac{\partial S}{\partial m}u.$$

(11)

3. Traveltime Tomography Method

3.1. Starting and Final Velocity Models

In this study, we used multichannel seismic data acquired in 1989 by Geological Survey of Canada using a 144-channel streamer with maximum offset of 3783 m [30]. The shot spacing was 50 m, and the airgun array (128.15 L) was towed at a depth of 10 m. The data were recorded for a time of 14 s at a sample interval of 4 ms. The data has excellent quality with high S/N ratio. No major problems were encountered during the data acquisition.

First arrival traveltimes are identified out to maximum offset. The apparent velocities of these arrivals are about $1.9$ km/s but increases up to $4.3$ km/s in some areas. An example of the data is shown in Figure 3. A starting model was carefully constructed using travel times from every 5th shot (every 250 m) along the line by ray tracing first in a block model [23] and then in a grid model [24]. The high-resolution grid model, which uses only first arrivals, is well-suited to accommodate a large number of shot gathers, since it is a much more automated process than the block model, where the structure must be specified layer-by-layer and reflector-by-reflector. The problem is that the automated approach in the high-resolution grid model is not very suitable when there are rapid lateral variations in velocity (Figure 3). In our initial application of the grid model, the inversion process would find the minimum traveltime residuals, but the residuals were still very large in the regions of large velocity change. By using the block model in a combined forward and inverse modeling approach, the flexibility in the model parametrization allowed us to obtain a better fit in these regions of high residuals (Figure 4a). The output velocities from the block model were used as input for the grid model with all the traveltime observations. The grid model uses node-based forward modeling and cell-based regularized inversion. Additional details about the inversion process and model parameterization are presented in Yelisetti [15].

A total of 139 shots with a spacing of 250 m were used for travel time inversion. The model length is $38.5$ km. Using the combination of forward modeling and damped-least squares inversion in a block model, we achieved a ${\chi }^2$ of $3.76$ and rms travel time residual of 12 ms, respectively (Figure 5a). Tomographic inversion with a uniform grid spacing of 25 m was used to obtain the final grid model (Figure 5b), which produced model traveltimes that agreed with the traveltime observations within the picking uncertainty of 6 ms. Velocity anomalies as large as $\pm 240$ m/s were observed in the velocity anomaly model that was obtained by subtracting the final grid model from the block model (Figure 5c).

Figure 4 shows the travel time residual plots using both the block model and the fine grid model. The large travel time residuals observed at 15–20 km model distance and at 26–36 km model distance (Figure 4c) from the block model were reduced using the fine grid model (Figure 4d). Observed and calculated travel times using these two models were plotted over model distance 24–38 km (Figure 4a,b). Arrows in Figure 4a,b indicate significant improvement in travel time residuals for the fine grid model relative to the block model.

3.2. Corrugation Tests

Corrugation tests were done to further assess the lateral resolution of the tomographic fine grid velocity model. We used a constant vertical velocity perturbation of $5\%$ to the final velocity model and constant corrugation widths of $0.5$ km, 1 km, and 2 km. Then, we traced the rays through the perturbed model and calculated the travel times using forward modeling. With the final velocity model as a starting model, these calculated times were used for inversion to obtain the corrugated velocity model. Perturbations were then calculated by subtracting the final velocity model from the corrugated velocity model. These tests showed that the model was able to resolve structures with a lateral resolution of ∼1 km to depths of 0.5–1 km (Figure 6). Anomalies of $0.5$ km width are resolved only in the near surface (∼250 m).

4. Waveform Inversion Strategy

Waveform inversion strategy used in this study closely follows that of Takougang and Calvert [31]. A typical shot gather showing the direct arrival, reflections and refractions is shown in Figure 7a. The corresponding frequency spectrum of the near-offset trace is shown in Figure 7b, which indicates the minimum available frequency for inversion of 6 Hz.

The starting velocity model derived from traveltime inversion was used in forward modeling to calculate modeled travel times. These are within $\frac{1}{2}$ cycle of the observed travel times (Figure 8), indicating that the initial velocity model is of sufficient quality to carry out inversion. The source and receiver ghosts were mitigated during the modeling by using a free surface boundary condition on the top of the model. Absorbing boundary conditions were used on the remaining three sides of the model. The waveform inversion process updates the starting velocity model at each iteration step by minimizing the misfit between the observed and synthetic waveforms. The forward propagated wavefield is multiplied with the back-propagated data residual to obtain the gradient of the misfit function [10].

Data preconditioning steps, such as low pass filtering, amplitude scaling, and time windowing, are required prior to inversion. The field data were low-pass filtered to 15 Hz to remove the high frequency noise, which is difficult to model during the inversion process. Amplitude scaling was performed to match synthetic data with raw data. The goal of preconditioning is to organize the field data in a form that is appropriate for waveform inversion. Therefore, any part of the data that is not predicted in the forward modeling using the 2-D acoustic wave equation should be corrected or removed. These features include shot-to-shot energy variations, shear waves, bad traces, coherent noise, and amplitude discrepancy. Furthermore, amplitude variations in real data may arise from acoustic attenuation effects, poor instrument coupling, noise contamination, or elastic effects. Thus, the real data must be pre-processed so that its amplitude versus offset behavior matches that of the modeled first arrivals. In order to match the amplitudes of both the data sets, we followed a method proposed by Brenders and Pratt [32], in which the observed data are multiplied by a constant computed from a linear regression analysis of amplitude variation with offset. After scaling, the observed data are transformed into frequency domain for input to inversion.

4.1. Time Windowing

A time window is commonly applied to the input seismic data before extracting the frequency domain data. The main goal of time windowing is to stabilize the inversion process by allowing only first arrivals, direct, and refracted energy from the seismic data in the initial stages of the inversion process [31]. Nonlinearity can be mitigated by excluding late arrivals, seafloor multiples and scattered energy that originate from outside the imaging plane. In addition, early arrivals helps in recovering the long wavelength or low frequency features of the model, which helps in stabilizing the inversion process. The time window size can be increased in the later stages of the inversion to incorporate late arrivals to image deeper structure (e.g., Reference [32]). The choice of time window will influence the inversion frequency interval. We used a window size of 2 s after the first arrival with maximum modeled time set to $2.5$ s at a sample interval of 4 ms.

4.2. Time Domain Source Signature

The source wavelet can be estimated using the procedures outlined in Pratt [21] if it is not known prior to the inversion. The traveltime inversion velocity model was used along with a spike wavelet to obtain the initial source wavelet. The source signature is usually updated while inverting for a new velocity model. The uniformity in estimated source signature (Figure 9a) from shot-to-shot indicates the quality of the starting velocity model. For our final inversion, we used the averaged source signature over all shots (Figure 9b).

4.3. Frequency Selection

There is no unique way to choose frequencies in waveform inversion. Frequencies are generally selected based on the Nyquist sampling theorem (e.g., Reference [33]). A frequency interval $\Delta f=\frac{1}{T_w}$ is needed to describe the data within the time window $T_w$. The spectrum of all frequencies within this interval for which the amplitudes are non-zero needs to be sampled. Frequencies within the spectrum are inverted sequentially rather than simultaneously, starting from low to high [14,15,17,31,34].

Over the frequency range from 6 to 14 Hz, our inversion used 21 frequencies at an interval of 0.4 Hz; for each frequency, there were 5 iterations [10]. These frequencies were chosen based on the Nyquist sampling criterion to avoid aliasing and time wrap-around effects. Attenuation is set to zero (Q = infinity) for velocity inversion. The grid spacing was $12.5$ m in both horizontal and vertical directions, which results in 3081 and 161 grid points in the horizontal and vertical directions, respectively. This fine grid interval was chosen based on the dominant frequency in the data, which is ∼30 Hz. Therefore, the smallest wavelength that can be modeled is 50 m, assuming a minimum model velocity of 1500 m/s. These parameters satisfy the criterion of four grid points per wavelength that is required to obtain $95\%$ accuracy in estimated numerical velocities.

The amplitude and phase of the source signature were iteratively updated during the inversion process. A 2D-bandpass filter was applied in the wavenumber domain to filter the gradient in both horizontal and vertical directions. We used values of 0–3 for $K_x$ and 1.5–14 for $K_z$. These numbers are calculated using the minimum model velocity and the dominant frequency present in the data.

5. Properties of Velocity Model Derived from Waveform Inversion

Our final velocity model in Figure 10b clearly resolves more detailed structures than the starting velocity model shown in Figure 10a. The velocity perturbations relative to the starting velocity model are clearly shown in the difference model obtained by subtracting the initial model from the final model (Figure 10c). Velocity anomalies as large as $\pm 150$ m/s are observed after the final iteration step.

Figure 11a indicates the variation in objective function at frequencies 6, 8, 10, 12, and 14 Hz. The relative reduction in objective function with respect to its previous update is shown in Figure 11b. Velocity perturbations as large as $\pm 150$ m/s are observed in the initial stages of the inversion process, which are reduced to $\pm 10$ m/s in the final stages of inversion (Figure 11c).

The final velocity model from full waveform inversion indicates a prominent low-velocity zone between 2–14 km model distance (Figure 10b), which is not observed in the traveltime inversion model (Figure 10a). To further assess this low-velocity zone, we computed synthetic shot gathers using the final velocity model. A sample of these shot gathers that pass through this low-velocity zone at model distance 6 km, 8 km, and 10 km are shown in Figure 12. These are compared to observed shot gathers and to synthetic shot gathers from the traveltime inversion model. Clearly, the raw shot gathers are very similar to those computed from the final velocity model; specifically, as highlighted in Figure 12, prominent reflections in the observed shot gathers are present in the final modeled shot gathers, but they are absent from gathers using the starting or traveltime inversion model. We confirmed that much of this reflectivity originates from the low-velocity zone by producing a model in which this low-velocity zone is removed from the final model, by smoothly interpolating structure from above and below the zone. Figure 13 shows the synthetic shot gathers with and without the low-velocity zone. The difference between these two sets of synthetic shot gathers are arrivals due to the low-velocity zone, and many of these correspond to reflections in the observed shot gathers. However, there is still a lot of reflectivity in the plot without the low-velocity zone in the near-offset depth range from 1.2–1.5 s (Figure 13b). Presumably, this comes from structure deeper than the low-velocity zone—as opposed to the reverberations from the low-velocity zone shown in the difference plot (Figure 13c).

To further assess our final velocity model, synthetic shot gathers are computed at model distance 14 km, 16 km, and 18 km, where a prominent high velocity zone is observed at depths below ∼1.0 km. These shot gathers are compared to observed shot gathers and synthetic shot gathers from traveltime inversion model (Figure 14). Again, there is a high degree of similarity between the observed shot gathers and those computed from the final velocity model. The improvement, relative to the initial velocity model from traveltime inversion, is mainly in the near-offset half of the gathers, where the starting model gathers typically have too few reflections, except for an overly prominent strong reflection indicated by yellow arrow in Figure 14.

The waveform inversion velocity model was used to obtain common offset gathers (Figure 15) to assess the quality of modeled data with respect to the observed data. The modeled common offset gather (Figure 15a) at 3308 m compares reasonably well with the corresponding observed common offset gather (Figure 15b). Particularly, an anticlinal structure is modeled between traces 40–65, which is also present in the observed data. Similarly, a sub-basin modeled between traces 65–120 matches reasonably well with the observed data. The quality of modeled data below $1.5$ s deteriorates, which could be associated with elastic effects that are not included in the acoustic inversion.

To see the overall quality of our model, we plotted the observed and modeled data in the frequency domain (Figure 16) corresponding to a frequency of 6 Hz. We can clearly follow the phase information in these two data sets, which is not readily observed in the difference plot (Figure 16c). This indicates that the final velocity model incorporates most of the phase information that is contained in the data. However, strong amplitude anomalies are observed at the far-offset (receiver numbers 0–50, Figure 16c), which could be associated with elastic effects, such as mode conversion, that are not included in the acoustic inversion. These amplitude anomalies could be improved using visco-elastic inversion.

6. Checker Board Tests

Checker board tests are used to assess the final velocity model. Five percent perturbation is applied to the final velocity model using a check size of 250 m and 500 m. For a velocity of 2500 m/s and a frequency of 10 Hz, the wavelength of 250 m is consistent with the minimum check size used in this study. Forward modeling and inversion are carried out using the perturbed velocity model. The inverted models are subtracted from the final velocity model. The difference between the final and inverted models are shown in Figure 17. The velocity perturbation results show that the final model is resolved well in the upper 1 km for both 250 m and 500 m check size. Moreover, the model is well resolved down to 2 km depth between the model distance of 12–21 km, where the high velocities are observed. This is a big improvement relative to the resolution from the corrugation tests for the traveltime model (Figure 6). The imaging depth also increased compared to the traveltime models as a result of the large time window used in waveform inversion.

7. Comparision with Drill Hole Data

Waveform inversion velocities are compared to sonic logs from Zeus I-65, Zeus D-14, and Pluto I-87 that were drilled by Shell in the 1960s (Figure 18). These wells are not exactly on the line, and were projected on to the seismic line by 18 km, 10 km, and 11 km, respectively. Zeus I-65 and Pluto I-87 were drilled on top of anticlinal structures. Zeus I-65 penetrated through the shallow part of the Crescent terrane, while Pluto I-87 penetrated through a thin gas-bearing sand at depth. Zeus D-14 penetrated through volcanic rocks on the flank of an anticlinal structure [4]. Since the original logs were not available, the digitized logs of Zeus I-65, Zeus D-14, and Pluto I-87 were assigned uncertainty values of $\pm 200$ m/s, $\pm 100$ m/s, and $\pm 150$ m/s, respectively, based on the wiggles in the sonic logs. This large error of $\pm 200$ m/s for a $0.5$ km thick layer with an average velocity of 2 km/s would only result in a two-way travel time error of approximately $\pm 0.05$ s. Even though the wells are projected on to the line, the inversion velocities match reasonably well with the sonic logs (Figure 18), down to the maximum modeled depth of ∼2 km, which indicates the reliability of our inversion results. Seismic velocities at Zeus-D14 well in the depth interval 1.5–2.0 km (Figure 18) appear to be higher at the projected location of the well than at the well itself, which is 10 km from the seismic line. This discrepancy could be related to sub-surface structural variation along the margin.

Figure 18. One-dimensional-velocity profiles through travel time inversion model (red) and full waveform inversion model (blue) along with sonic log velocties (green). Arrows in Figure 2 show locations of model velocity profiles that correspond to sonic logs. Solid black arrows indicate locations where Eocene volcanics were found in the Shell Canada exploratory wells from 1967, which have been re-evaluated by several studies, including Narayan et al. [11], Hayward and Calvert [12], and Johns et al. [5]. Lithology variations with depth from Zeus-I65, Zeus-D14 and Pluto-I87 are shown in (b), (c), and (d), respectively.

Figure 18. One-dimensional-velocity profiles through travel time inversion model (red) and full waveform inversion model (blue) along with sonic log velocties (green). Arrows in Figure 2 show locations of model velocity profiles that correspond to sonic logs. Solid black arrows indicate locations where Eocene volcanics were found in the Shell Canada exploratory wells from 1967, which have been re-evaluated by several studies, including Narayan et al. [11], Hayward and Calvert [12], and Johns et al. [5]. Lithology variations with depth from Zeus-I65, Zeus-D14 and Pluto-I87 are shown in (b), (c), and (d), respectively.

8. Results and Discussion

The final velocity model obtained from acoustic inversion is converted to two-way travel time, assuming vertical incidence rays, for direct comparison with the migrated seismic stack section along line 89-06 (Figure 19b). The magnetic anomaly along this line is also plotted (Figure 19a).

8.1. Velocity Structure

8.1.1. Low-Velocity Zone (LVZ)

Within the accreted sediments beneath the shelf, a prominent low-velocity zone is observed in the southwestern portion of the line, where the velocities are about 200 m/s lower than the surrounding velocities (Figure 10b and Figure 19b). The LVZ is located at depths of ∼600–1000 m below seafloor (mbsf) and extends laterally for ∼10 km. Superposition of the velocity model on the migrated seismic section indicates that this LVZ is associated with a strong seismic reflection (Figure 19b). A LVZ has not been observed in the previous studies [12] due to resolution limits imposed by traveltime inversion and conventional seismic reflection imaging. The observed low-velocity zone is likely a result of high pore pressure associated with fluid expulsion in the accreted wedge sediments under the continental shelf [10,35]. Hayward et al. [36] also observed a LVZ in association with high fluid pressures at the Barbados accretionary prism based on the vertical seismic profile observations. Shell exploratory wells in the Tofino basin also show high fluid pressure (lithostatic gradients) [1,4]. The low-velocity zone can also arise from over-pressured gas [4] or from a high porosity layer associated with lithology changes [10].

The LVZ identified here is associated with a possible gas migration pathway (indicated by a black arrow in Figure 19b) that reaches up to the seafloor at ∼14–15 km model distance. However, a coincident fault is not observed on the migrated seismic section. Another localized LVZ is identified on the northeast portion of the section at 750 m below the seafloor near 35–36 km model distance (Figure 10b). It is observed below a broad anticlinal feature and is associated with strong reflections as seen from the seismic section (Figure 19b). This feature may be associated with gas migration from below.

8.1.2. Crescent Terrane

A prominent high velocity feature is observed between 12–18 km model distance (Figure 10b). The velocities in this region increase from 2000 m/s to 2500 m/s at a shallow depth of 750 m. Below this depth, the velocities uniformly increase to about 5000 m/s at 2000 m depth. These high velocities likely correspond to the shallowest occurrence of the volcanic Crescent terrane [10]. A magnetic high is also observed over this region, which provides further evidence for the presence of Crescent volcanics (Figure 19a and Spence et al. [6]). However, the magnetic high occurs ∼5 km landward of the shallowest occurrence of the volcanic Crescent terrane. This indicates that the Crescent terrane is highly fractured where it is shallow or otherwise demagnetized. A diapiric feature or an anticlinal structure is observed between 17–22 km model distance (Figure 19b). This feature is located ∼2–3 km landward of the shallowest occurence of the volcanic Crescent terrane. This is probably associated with the seaward verging Tofino fault, which was interpreted as the top boundary of the Crescent terrane [6,10].

8.1.3. Pacific Rim Terrane

As seen in the final velocity model, in general, the shallow sediment velocities in the Tofino basin increase landward (Figure 19b). This is likely associated with the compressional environment associated with the ongoing subduction process [14]. Lateral variations in velocity model are better seen by subtracting the average 1D-velocity profile from the final FWI velocity model (Figure 19c). As a specific example, velocities in the northeastern portion of the model (between 32 and $38.5$ km model distance) at a depth of ∼0.2–0.5 km are approximately $0.6$ km/s higher than average velocities west of ∼30 km model distance. The shallow high velocities suggest a shallower depth for the transition to the Pacific Rim terrane.

9. Interpretation

In general, the modeled velocities can be summarized in to three layers. (i) Average velocities of up to 2 km/s in the shallow sediment Section (500–1200 m), represented by mostly red-yellow in the velocity color scale (Figure 19b). These may correspond to interbedded sediments with muddy siltstone/silty mudstone. Similar lithologies were found at equivalent depths in the wells. (ii) Below this, the velocities of up to $3.2$ km/s are observed in the depth interval that is represented by mostly greens in the velocity color scale, which corresponds to mudstone. (iii) The blues in the velocity color scale indicate >4 km/s which corresponds to Eocene volcanics.

Seismic velocities are >4 km/s below ∼2.0 km depth near Zeus-I65 well, which indicate the presence of Eocene volcanics at this depth (Figure 18b). However, at Zeus-D14, our seismic velocities from waveform inversion do not extend sufficiently deep to identify the top of massive volcanics, which may be observed at 2.2–2.5 km sub-surface depth based on the velocities observed at Zeus-I65.

Within the sediment section, seismic velocities at ∼1.0–2.0 km sub-surface depth near Pluto-I87 well indicate velocity undulations as large as $\pm 200$ m/s (Figure 18d). These could be related to alternate layers of fluid or gas saturation associated with bending-related faulting [37].

The thickness of the sediment fill increases towards fossil trench region and then decreases over the Pacific Rim terrane (Figure 19b). This fossil trench region [1] is bounded by the basement highs of the Crescent terrane and the Pacific Rim terrane. The sediment structures show steep dips at the western edge of the Pacific Rim terrane as opposed to gentle dips observed over the Crescent basement high. Comparision of this line with other seismic lines along the margin indicates that the fossil trench region has its greatest thickness along this line and that sediment deformation varies along the margin [6,12].

Sediment velocities closely follow the shape of the diapiric feature observed at 18–20 km model distance. This is in contrast to the south where the sediment velocities decrease in the core of the anticline on multichannel lines 85-01 and 89-02 [12], which indicates the structural variation along the margin. Furthermore, the sediment velocities are smoother and less deformed on the southwestern portion of the line where sediments overlie the accreted wedge, compared to the northeast where the sediments are more deformed (greater folding and faulting), indicating that the deformation increases landward.

10. Conclusions

In this study, frequency domain acoustic full waveform inversion was used to derive subsurface velocity model down to 2.0 km depth using controlled source seismic data from Vancouver Island continental shelf. This study shows that the waveform inversion is quiet effective in imaging geological features using short streamer data when proper preconditioning and inversion strategy are used. Specifically, when applied to marine multichannel seismic reflection data across the Tofino fore-arc basin beneath the Vancouver Island shelf, the inversion enables the identification of a prominent low-velocity zone at a depth of 800–900 m within the Tofino basin sediment section, extending over a 10 km lateral distance. The low velocity zone could be associated with (1) over-pressured gas, (2) a high porosity layer associated with lithology changes, and (3) fluid over-pressure. High seismic velocities (∼4–5 km/s) at depths of 1.5–2.0 km enable the identification of the top boundary of the Eocene volcanic Crescent terrane. In general, sediment seismic velocities increase landward, likely indicating the over-consolidation of shelf sediments. Furthermore, the sharp landward increase in sediment velocity at ∼10 km west of Vancouver Island indicates the underlying transition to the Pacific Rim terrane. The waveform inversion velocity modeling results from this study could be useful for future hydrocarbon drilling in this area.