Diving-wave time-lapse delay for CO₂ thin layer detection

SUMMARY

We have derived an analytical approximate expression to estimate the delay in diving seismic waves due to thin layers of CO₂. The expression is valid for high frequencies and can be used to estimate the delay in diving waves at seismic frequencies for large separations between the source and receiver (offset). The approximation may be used to assess CO₂ detection limits using diving waves and to support survey planning for CO₂ monitoring and full-waveform inversion (FWI) cycle skipping analysis. In this study, we analyse the diving-wave response to a thin layer of CO₂ for band-limited data using acoustic finite-difference modelling, and compare the results against the analytical calculations. We find that the responses are offset-dependent and related to double- and single-leg interactions between the diving waves and the CO₂. To test the methods, we created a synthetic representation of the 2010 subsurface conditions for the top CO₂ layer at the Sleipner storage complex in the North Sea, by combining base and monitor post-stack seismic data with field velocity trends. Using the acoustic finite-difference method, we model pre-stack data that captures the complexity of field data and demonstrate the use of the diving-wave delay for CO₂ migration monitoring and CO₂ thin layer detection.

Numerical modelling, Time-series analysis, Body waves

INTRODUCTION

Detection of thin layers of CO₂ in saline aquifer formations is a key objective for carbon sequestration monitoring projects. After injection, CO₂ migrates upwards and laterally, forming an inverted cone whose shape is controlled by the interplay between viscous, capillary and gravity forces (Ringrose et al. 2021). Due to gravity and rock heterogeneity, the cones tend to have a low aspect ratio, developing into thin and elongated layers of CO₂. This behaviour can be observed at the Sleipner storage site, where most of the CO₂ layers are thinner than one-quarter of the propagated seismic wavelength (Cowton et al. 2016).

Time-lapse seismic reflection analysis has been the primary tool for monitoring the migration of CO₂ in field storage projects. The time-lapse seismic difference reveals both amplitude and traveltime changes, which are used to map the boundaries of the CO₂ accumulation. Bandwidth limitations and the strong tuning interference (Widess 1973) between the reflections from the top and the base of a thin layer of CO₂ limit the kinematic resolution of these layers using the reflection method. Combining several techniques, such as direct mapping of the CO₂ layer in seismic sections, amplitude/push-down time-shift analysis and spectral decomposition techniques, have narrowed the detection levels for shallow settings, such as the top CO₂ layer at Sleipner, down to ∼1 m (White et al. 2018).

Time-lapse analyses using transmitted energy offer a complementary tool to reflection seismic analyses. The refraction method has been applied in various studies to monitor changes in subsurface conditions over time, such as ground ice changes (Hilbich 2010) and shallow gas migration (Landrø et al. 2019). Zadeh & Landrø (2011) derived approximations for the refraction time-shift due to a 2-D anomaly with finite lateral extension in the lower half-space and studied the shift on field data to monitor an underground blowout. Haavik & Landrø (2014) used refraction time-shift analysis to monitor gas flow in and out of a tunnel valley. Landrø et al. (2021) provided an analytical approximation for the diving-wave time delay due to the presence of a thin gas layer in a constant velocity gradient medium with a water column. However, their scheme ignores the ray deflection in the low-velocity CO₂ layer, resulting in a large overestimation of the time delay. Time-lapse analyses using diving waves are not affected by the tuning effect in the thin CO₂ layers and allow the measurement of kinematic delays that can be used for the monitoring of CO₂ injection. Moreover, the diving waves can provide information about thin layers of CO₂ at low frequencies. They are also sensitive to the low velocities in the CO₂, which makes it possible to link their traveltime delay to the properties of the CO₂ layer. Diving waves also sample the medium- and long-wavelength components of the velocity structure of the subsurface (Shipp & Singh 2002), enabling successful application of nonlinear, ill-posed inverse problems such as full-waveform inversion (FWI, Shipp & Singh 2002; Virieux & Operto 2009; Kazei et al. 2013; Raknes et al. 2015). Despite the potential of the long-offset time-lapse method, its application to the Sleipner CO₂ storage complex is not possible due to the lack of long-offset recordings in the baseline seismic survey, and evaluating its use requires numerical modelling.

In this study, we present an analytical approximation to the diving-wave time-lapse delay due to a thin layer of CO₂ injected into a constant gradient aquifer. We perform numerical modelling of the delay for band-limited sources using the acoustic finite-difference method and assess the interactions between the diving waves and a thin layer of CO₂ using simple subsurface models. We compare our analytical delay against the results from the numerical modelling. Finally, we combine post-stack reflectivity data and field velocity trends to create a high-resolution representation of a hypothetical single-layer Sleipner scenario, which captures the complexity of field data. The models are used to generate realistic pre-stack synthetic data, which are then utilized to demonstrate the use of the diving wave delay for monitoring thin layers of CO₂.

DIVING-WAVE TIME-LAPSE DELAY DUE TO A THIN CO₂ LAYER

Consider an isotropic baseline background medium where the velocity v increases linearly with depth z:

$$\begin{eqnarray} v\left( z \right) = \ {v}_\mathrm{ o}\ + \ gz, \end{eqnarray}$$

(1)

where |${v}_\mathrm{ o}$| is the velocity at the surface and g the velocity gradient.

The traveltime in the background medium from the source to the receiver can be computed using the diving-wave moveout equation for constant gradient medium (Stovas & Alkhalifah 2014):

$$\begin{eqnarray} {T}_{{\mathrm{ base}}}\left( x \right) = \frac{2}{g}\log \left( {\frac{{gx}}{{2{v}_\mathrm{ o}}} + \sqrt {1 + \frac{{\ {g}^2{x}^2}}{{4v_\mathrm{ o}^2}}} \ } \right), \end{eqnarray}$$

(2)

where x is the lateral distance between the source and the receiver (offset).

We next assume CO₂ is injected into a subsurface aquifer, where both the aquifer and the seal follow the background velocity trend. The CO₂ forms an infinitely long, thin low-velocity layer at depth |${z}_1\ $| (see Fig. 1). The layer has half-thickness |$\mathrm{ d}z$| and velocity |${v}_g,$| with top and base located at |${z}_1-\mathrm{ d}z$| and |${z}_1\ + \ \mathrm{ d}z$|⁠, respectively. The CO₂ layer will cause a deflection of the diving rays from source to receiver in Fig. 1, resulting in two additional travel paths that delay the signal with respect to the baseline medium: a straight path within the CO₂ layer and a diving wave path in the underburden. For a common-offset scenario, the overburden travel path also changes.

Figure 1.

Common-offset travel paths for the base (B) and monitor (M) media. After injection, the CO₂ forms a thin low-velocity layer that deflects the rays, yielding a traveltime delay of the monitor with respect to the baseline medium.

The one-way offset and traveltime for this diving ray segment travelling from the top of the medium to the top of the CO₂ layer are computed using Stovas & Alkhalifah (2014) as:

$$\begin{eqnarray} {x}_{\mathrm{ ob}}\left( p \right) = \frac{1}{{pg}}\left( {\sqrt {1 - {p}^2v_\mathrm{ o}^2} - \sqrt {1 - {p}^2{{\left( {{v}_\mathrm{ o} + g\left( {{z}_1 - \mathrm{ d}z} \right)} \right)}}^2} } \right), \end{eqnarray}$$

(3)

$$\begin{eqnarray} {T}_{\mathrm{ ob}}\left( p \right) = \frac{1}{g}\log \left( {\frac{{\left( {{v}_\mathrm{ o} + g\left( {{z}_1 - \mathrm{ d}z} \right)} \right)\left( {1 + \sqrt {1 - {p}^2v_\mathrm{ o}^2} } \right)}}{{{v}_\mathrm{ o}\left( {1 + \sqrt {1 - {p}^2{{\left( {{v}_\mathrm{ o} + g\left( {{z}_1 - \mathrm{ d}z} \right)} \right)}}^2} } \right)}}} \right), \end{eqnarray}$$

(4)

where, for departing angle |${\theta }_i$|⁠, |$p\ = \ \mathrm{ sin}( {{\theta }_i} )/{v}_o$| is the horizontal slowness in the monitor medium.

The one-way offset and traveltime within the constant-velocity CO₂ layer are computed using the offset-traveltime parametric equations:

$$\begin{eqnarray} {x}_{\mathrm{ tl}}\left( p \right) = \frac{{2\mathrm{ d}zp{v}_g}}{{\sqrt {1 - {p}^2v_g^2} }}, \end{eqnarray}$$

(5)

$$\begin{eqnarray} {T}_{\mathrm{ tl}}\left( p \right) = \frac{{2\mathrm{ d}z}}{{{v}_g\sqrt {1 - {p}^2v_g^2} }}. \end{eqnarray}$$

(6)

The one-way offset and traveltime in the underburden, from |$z\ = \ {z}_1\ + \ \mathrm{ d}z$| down to the post-injection maximum penetration depth, are calculated using Stovas & Alkhalifah (2014):

$$\begin{eqnarray} {x}_{\mathrm{ ub}}\left( p \right) = \frac{{\sqrt {1 - {p}^2{{\left( {{v}_\mathrm{ o} + g\left( {{z}_1 + \mathrm{ d}z} \right)} \right)}}^2\ } }}{{pg}}, \end{eqnarray}$$

(7)

$$\begin{eqnarray} {T}_{\mathrm{ ub}}\left( p \right) = \frac{1}{g}\log \left( {\frac{{1 + \sqrt {1 - {p}^2{{\left( {{v}_\mathrm{ o} + g\left( {{z}_1 + \mathrm{ d}z} \right)} \right)}}^2\ } }}{{p\left( {{v}_\mathrm{ o} + g\left( {{z}_1 + \mathrm{ d}z} \right)} \right)}}} \right). \end{eqnarray}$$

(8)

The total two-way offset and traveltime after CO₂ injection are computed by adding the contributions from the three path segments:

$$\begin{eqnarray} {x}_{{\mathrm{ monitor}}}\left( p \right) = 2\left( {{x}_{\mathrm{ ob}} + {x}_{\mathrm{ tl}} + {x}_{\mathrm{ ub}}} \right), \end{eqnarray}$$

(9)

$$\begin{eqnarray} {T}_{{\mathrm{ monitor}}}\left( p \right) = 2\left( {{T}_{\mathrm{ ob}} + {T}_{\mathrm{ tl}} + {T}_{\mathrm{ ub}}} \right). \end{eqnarray}$$

(10)

Finally, the diving-wave time-lapse delay in the infinite frequency limit due to a layer of CO₂ can be computed analytically as a function of offset as:

$$\begin{eqnarray} \mathrm{ d}t\ = {T}_{{\mathrm{ monitor}}}\ -\ {T}_{{\mathrm{ base}}}, \end{eqnarray}$$

(11)

where |$\ {T}_{{\mathrm{ monitor}}}\ = \ {T}_{{\mathrm{ monitor}}}( p )$| and |${T}_{{\mathrm{ base}}}\ = \ {T}_{{\mathrm{ base}}}( {{x}_{{\mathrm{ monitor}}}( p )} )$|⁠. The delay in eq. (11) is non-zero for |$x\ > \ {x}_{{\mathrm{ monitor}}}( {p = 1/v( {{z}_1 + \mathrm{ d}z} )} ).$|

ANALYTICAL EXPRESSION AS A FUNCTION OF OFFSET

By expanding |$\mathrm{ d}t$| in eqs (4), (6) and (8) in series as a function of |$\mathrm{ d}z$|⁠, and neglecting second order and higher terms, we obtain

$$\begin{eqnarray} \mathrm{ d}t \approx \frac{{4\mathrm{ d}z\sqrt {1 - {p}^2v_g^2} }}{{{v}_g}} - \frac{{4\mathrm{ d}z\sqrt {1 - {p}^2\left( {{v}_\mathrm{ o} + g{z}_1} \right)_a^2} }}{{{v}_\mathrm{ o} + g{z}_1}}, \end{eqnarray}$$

(12)

where p > 0.

We approximate p to that in the background medium as:

$$\begin{eqnarray} p \approx \frac{2}{{\sqrt {4v_\mathrm{ o}^2 + {g}^2{x}^2} }}. \end{eqnarray}$$

(13)

Substituting eq. (13) in eq. (12) yields an approximation to the double-leg time-lapse diving-wave delay as a function of offset:

$$\begin{eqnarray} \mathrm{ d}t\left( x \right) \approx \frac{{4\mathrm{ d}z\sqrt {1 - \frac{{4v_g^2}}{{4v_\mathrm{ o}^2 + {g}^2{x}^2}}} }}{{{v}_g}} - \frac{{4\mathrm{ d}z\sqrt {\frac{{g\left( { - 8{v}_\mathrm{ o}{z}_1 + g\left( {{x}^2 - 4z_1^2} \right)} \right)}}{{4v_\mathrm{ o}^2 + {g}^2{x}^2}}} }}{{{v}_\mathrm{ o} + g{z}_1}}. \end{eqnarray}$$

(14)

Fig. 2 shows a comparison between the exact and approximate traveltime delays calculated using eqs (11) and (14), respectively. Model parameters are: |$\mathrm{ d}z$| = 5 m, |${z}_1$| = 718 m, |${v}_g$| = 1500 m s⁻¹, |${v}_o$| = 1657 m s⁻¹ and g = 0.73 s⁻¹. The approximation is valid for offsets greater than the offset at which the ray becomes horizontal in the background medium at |${z}_1$|⁠. The error is maximum at this offset (∼0.86 ms, ∼10.6 per cent for our model) and tends to zero as offset increases.

Figure 2.

Exact versus approximate diving-wave double-leg delay as a function of offset due to a thin and infinitely long layer of CO₂ within a constant gradient medium.

GENERALIZED DIVING-WAVE DELAY FOR 2-D CO₂ MODELS

The expression in (14) assumes an infinitely long CO₂ layer, and up- and downgoing interactions (ray segments) from the source and receiver sides between the diving wave and the CO₂ layer. Under the assumption of the same p between baseline and monitor media, it can be shown that separate contributions from the source and receiver sides, and contributions from multiple layers of CO₂, can be linearly added together as:

$$\begin{eqnarray} \mathrm{ d}t\left( x \right) \approx \mathop \sum \limits_{k = 1\ }^{{N}_l} \mathop \sum \limits_{i = 1}^2 \frac{2\mathrm{ d}z_{i}\sqrt{1-\frac{4v^2_{gi}}{4v^2_{\mathrm{ o}}+g^{2}x^{2}} }}{v_{gi}}- \frac{2\mathrm{ d}z_{i} \sqrt\frac{g\left(-8v_{\mathrm{ o}}z_{1i}+ \left(x^{2}-4z_{1i}^2 \right) \right)}{4_{\mathrm{ o}}+g^{2}x^{2}}}{v_\mathrm{ o}+ gz_{1i}}, \\ \end{eqnarray}$$

(15)

where |${N}_l\ $| is the number of CO₂ layers. The subscript i denotes the source and receiver side properties, and the sum over k gathers the contributions of individual CO₂ layers at different depths. This expression enables the handling of 2-D CO₂ models (⁠|${v}_g,\ {z}_1,\ \mathrm{ d}z$|⁠), single-leg interactions between the diving waves and the CO₂, and multiple layers of CO₂ at different depths. Examples of 2-D CO₂ models are a CO₂ layer of finite extent (dz = 0 where there is no CO₂), or a lateral variation in CO₂ saturation described by a variable |${v}_g$|_. For 2-D CO₂ geometries, a pre-calculation step to compute source and receiver side models for the CO₂ properties (⁠|${v}_{i},{z}_{1i},\mathrm{ d}{z}_{i} $|⁠) is required.

The depth of the ray in the overburden can be computed as (Margrave & Lamoureux 2003):

$$\begin{eqnarray} {z}_{\mathrm{ ob}}\left( {p,d} \right) = \frac{1}{{pg}}\left( {\sqrt {1 - {{\left( {pgd - \sqrt {1 - {p}^2v_\mathrm{ o}^2} } \right)}}^2} - p{v}_\mathrm{ o}} \right), \end{eqnarray}$$

(16)

where d is the lateral distance from the source (or receiver).

We next find d for which the diving ray hits the top of the CO₂ layer by:

$$\begin{eqnarray} \mathop {\min }\limits_d \ \ {z}_{\mathrm{ ob}}\left( d \right)\ - \ {z}_1\left( {{x}_{\mathrm{ src}}\ + \ d} \right)\ = \ 0\ , \end{eqnarray}$$

(17)

$$\begin{eqnarray} \mathop {\min }\limits_d \ \ {z}_{\mathrm{ ob}}\left( d \right)\ - \ {z}_1\left( {{x}_{\mathrm{ src}} + x - d} \right)\ = \ 0\ , \end{eqnarray}$$

(18)

where x_src is the lateral coordinate of the source.

Finally, the source models |${v}_{g1} = {v}_{g}( {{x}_{\mathrm{ src}} + d} ), {z}_{1\,1} = {z}_{1}( {{x}_{\mathrm{ src}} + d} ), \mathrm{ d}{z}_{1} = \mathrm{ d}z( {{x}_{\mathrm{ src}} + d} ), $| and receiver models |${v}_{g2} = {v}_{g}( {{x}_{\mathrm{ src}} + x - d} ), {z}_{1\,2} = {z}_{1}( {{x}_{\mathrm{ src}} + x - d} ), \mathrm{ d}{z}_{2} = \mathrm{ d}z( {{x}_{\mathrm{ src}} + x - d}) $|⁠, are inserted in eq. (15) for delay estimation.

RESULTS

Numerical modelling on a simple 2-D model

The analytical expression estimates the delay due to the diving waves interacting with the CO₂. The main contributions to this delay are the slower traveltimes through the CO₂ layer compared to the baseline, and the longer path in the underburden for the monitor, as shown in Fig. 1. In the field, we record not only the diving waves but also a series of events which may interfere with the diving waves, as shown in Fig. 3.

Figure 3.

Schematic representation of events recorded in the field for long offsets: diving waves (black), guided waves (white), CO₂ layer post-critical reflection (red) (only base reflection shown here), reverberations within the CO₂ layer (dashed blue), diving-wave internal multiples (orange) and diving-wave multiples underneath CO₂ layer (dark blue). Other events recorded but not shown here are: source ghosts, receiver ghosts, receiver side of source ghosts, CO₂ layer free-surface multiples, etc.

To understand the interactions between the diving waves and a thin layer of CO₂, and the interplay with other coherent events, we performed 2-D acoustic finite-difference modelling using simple subsurface models.

First, we analyse the diving wave responses using velocity models that exclude the water column and include an absorbing free surface to avoid free-surface reflections and multiples. This enables assessment of the diving wave responses to the CO₂ in the same conditions as for the analytical approximation. Our baseline velocity model consists of a constant velocity gradient medium with parameters |${v}_\mathrm{ o}\ $| = 1657 m s⁻¹ and g = 0.73 s⁻¹, as shown in Fig. 4(a). These trend parameters are derived from least-square fitting on smoothed RMO (residual moveout) Dix-converted velocity cubes from the overburden at the Sleipner CO₂ storage complex in the North Sea. The CO₂ layer parameters for the monitor scenario are: dz = 5 m, |${z}_1$| = 718 m and |${v}_g\ $|= 1500 m s⁻¹, which are the same as used for the analytical calculation. The CO₂ layer is 3 km long and the density of the CO₂ layer is 1000 kg m⁻³. The source and receiver depth is 10 m. The sources are fired at an interval of 25 m, yielding 50 m offset classes and 25 m spacing between common-mid-point (CMP) traces. The maximum separation between the source and receiver is 6 km.

Figure 4.

Simple 2-D velocity models: (a) without water column and (b) with water column. The baseline model is a constant gradient velocity trend with parameters |${v}_\mathrm{ o}$| = 1657 m s⁻¹ and g = 0.73 s⁻¹. After injection, a 10 m thin, 3 km long CO₂ layer is formed within a reservoir that follows the baseline trend. Source locations are shown by the red lines at the top of the models. The boat sails left to right.

Fig. 5 shows a 250 Hz peak frequency, far-offset section, showing the delay in the diving waves due to their interaction with the CO₂. We see five delay regions depending on the position of the source with respect to the CO₂ layer, and the type of interaction that occurs. When the diving wave does not interact with the CO₂, the arrival time is constant and equal to that in the background medium (regions 1 and 5). As the position of the source moves (left to right), we observe delays related to single-leg interactions from the source side, double-leg interactions from the source and receiver sides, and receiver-side interactions (regions 2, 3 and 4). Note also the complex diffraction patterns from the top and base of the CO₂ layer for the source and receiver sides.

Figure 5.

High-frequency common-offset seismic section showing the delay in the diving waves due to a thin layer of CO₂ within a constant gradient medium. The delay is caused by single-leg interactions from the source and receiver sides, and double-leg interactions between the diving waves and the CO₂. No interaction takes place on regions 1 and 5. The section is extracted from an offset near the onset of interaction with the CO_2.

Fig. 6 shows the diving-wave delay due to the CO₂ computed from trace-by-trace cross-correlation between the base and monitor surveys for offsets beyond the onset of interaction. The time-shift maps show a decrease in delay magnitude as a function of offset. Near the onset of interaction, we observe single-leg interactions from the source and receiver (regions 2 and 4, respectively) and double-leg interactions between the diving waves and the CO₂. As offset increases, the single-leg regions separate spatially and the double-leg region ceases to exist owing to the finite lateral extent of the CO2 layer. The asymmetry of the blue region in the plot is due to the one-sided recording spread.

Figure 6.

Cross-correlation time-shift map for seismic data modelled using a 250 Hz peak frequency source. The dashed line represents the location of the seismic section shown in Fig. 5.

Fig. 7 shows cross-correlation time-shift maps computed on data modelled using sources of varying peak frequency. The high-frequency delay (250 Hz) shows a strong similarity with the delay computed using infinite-frequency approximation (eq. 15, shown at the bottom of Fig. 7), both in terms of offset patterns and time-shift magnitude. For seismic frequencies (∼20 Hz) the offset patterns exhibit a smoother appearance. For these frequencies, the offsets near the onset of interaction show larger time-shifts than for the analytical approximation. This is explained by interference between the diving wave and the CO₂ layer reflection. For offsets near the onset of interaction, these events cannot be separated kinematically, causing a tuning response that increases the magnitude of the delay computed by cross-correlation. This response is further exacerbated at the low frequencies (5 Hz) where the interference is strongest. At far offsets, the diving waves separate kinematically from the CO₂ layer reflection, reducing the interference and yielding a delay closer to the analytical calculation.

Figure 7.

Cross-correlation time-shift maps for seismic data modelled using sources with varying peak frequency: (a) 5 Hz, (b) 20 Hz and (c) 250 Hz. (d) Time-shift computed using eq. (15). Dashed lines highlight the offset for the sections shown in Fig. 8.

Fig. 8 shows common-offset sections for the monitor data at the far-offset location highlighted in Fig. 7. At 250 Hz, the primary signal is well separated from internal CO₂ layer reverberations, while these interfere at lower frequencies. Diffractions are also visible. At low frequencies, the diffractions interfere with the diving-wave arrival, and the cross-correlation captures the time difference with respect to the diffraction, yielding an erroneous time-shift. Nevertheless, for far offsets, the delays at high and seismic frequencies, shown at the top of the plots, are similar.

Figure 8.

Far-offset sections for monitor data modelled using sources of varying peak frequency. The curves at the top represent the cross-correlation time-shift with respect to the baseline. For far offsets, the seismic delay is similar to the high-frequency delay.

In the field, the interface between water and air acts as a reflecting surface that increases the complexity of the acquired wavefield by adding ghosts and multiples. Therefore, we analyse next the diving wave responses to the CO₂ in the presence of a water column and free surface multiples. The linear trend parameters for the baseline are the same as in Fig. 4(a), but shifted to incorporate the water column. The baseline trend and the monitor velocity models are shown in Fig. 4(b). The water depth is set at 82 m, which is the average water depth in the Sleipner CO₂ storage complex. Fig. 9 shows cross-correlation delay maps for various frequencies. We see only minor differences compared to the time-shift maps in Fig. 7, where no water column and multiples are present. Fig. 10 shows the monitor data in the presence of the free surface. At high frequencies, we observe the primary signal, the source and receiver side ghosts interfering constructively, and the receiver side of the source ghost. We also notice the internal CO₂ layer reverberations, and the first-order water column multiple. Despite the interference between these events at the low frequencies, and the presence of the water column, the diving wave delay does not change significantly from before, and the analytical equation (even without a water layer) can explain the high-frequency and the far-offset seismic delays.

Figure 9.

Cross-correlation time-shift maps assuming a reflecting free surface and a water column, for data modelled using sources with varying peak frequency: (a) 5 Hz, (b) 20 Hz and (c) 250 Hz. (d) Time-shift computed using eq. (15). Dashed lines highlight the offset for the sections shown in Fig. 10.

Figure 10.

Far-offset sections for monitor data modelled using sources of varying peak frequency assuming a reflecting free surface and a water column. The curves at the top represent the cross-correlation time-shift with respect to the baseline. For far offsets, the seismic delay is similar to the high-frequency delay.

Diving-wave delay on realistic subsurface models: 2010 single-layer Sleipner

Velocity Models

To assess the delay in the diving waves due to a thin layer of CO₂ under realistic subsurface conditions, we built a realistic velocity model inspired by the Sleipner CO₂ storage complex in the North Sea. The model represents a single CO₂ layer scenario for the 2010 top-layer conditions in this field. The model is created by combining high-frequency velocity changes derived from post-stack seismic data from the baseline 1994 and the monitor 2010 seismic surveys, with smooth velocity trends from the field. The velocity trend is the RMO Dix-converted velocity volume from 2010. This field contains long- and mid-wavelength velocities associated to the slowdown in the CO₂ plume. The velocity cube has therefore been edited to suppress this slowdown, smoothed and then corrected for the velocity push-down, yielding a smooth velocity trend representative of the pre-injection conditions.

To estimate the high-frequency velocity changes, we assume that the seismic post-stack section is a bandpassed representation of the Earth's reflectivity, which is the derivative of the acoustic impedance (product of velocity and density). Thus, by performing trace-by-trace mathematical integration of the post-stack cubes we can estimate the velocity changes. We then integrate the 1994 post-stack data and estimate the velocity changes for the baseline. These are further edited to attenuate non-primary signal such as remnant first-order water column multiple in the near-surface. Finally, the edited velocity changes are combined with the smooth velocity trend to yield a baseline model, as shown in Fig. 11(a).

Figure 11.

Single CO₂ layer model inspired by the top CO₂ layer at Sleipner: (a) baseline, (b) monitor and (c) profile view for baseline (black) and monitor (yellow) with location highlighted by the yellow triangles.

The mathematical integration is also applied to the 2010 post-stack data to derive a velocity change for the top CO₂ layer. The band-limited reflectivity of the CO₂ layer varies laterally as a function of the CO₂ layer thickness due to tuning interference between the reflections from the top and base of the CO₂ layer (Widess 1973). Therefore, we incorporate only the structural information from the integrated seismic, and enforce a velocity change of 550 m s⁻¹. Such slowdown is in agreement with the rock physics analysis of Dupuy et al. (2017) for CO₂ replacing brine within the Utsira sand reservoir. Due to bandwidth limitations, the CO₂ layer inverted via mathematical integration is expected to be thicker than the true CO₂ thickness. Therefore, we performed additional thinning of the CO₂ layer via iterative smoothing and velocity clipping. Finally, we combined the velocity change in the CO₂ layer with the baseline model to create the monitor scenario, as shown in Fig. 11(b). The CO₂ layer is located at approximately 800 m in depth, with thicknesses around 10 m. The minimum thickness defined in the model is 4 m. The CO₂ layer velocity |${v}_g$| varies around 1500 m s⁻¹ with a minimum of 1480 m s⁻¹ (Fig. 12).

Figure 12.

CO₂ layer properties: (a) depth (⁠|${z}_1$|⁠), (b) thickness (2∗dz) and (c) velocity (⁠|${v}_g$|⁠). CO₂ thicknesses down to 4 m are defined in the model.

Numerical modelling

We have modelled a 3-D seismic survey using an acoustic VTI (vertical transverse isotropy) finite-difference method. The acquisition comprises a single gun, single cable geometry, with shots spaced 50 m in x and y, yielding offset classes of 100 m, and a trace spacing between CMP positions of 50 m. The extent of the shot carpet is 12 000 m by 2400 m, as shown in Fig. 13. The shot carpet extends further north than the shot carpet of the true field data. The idea behind this extension is to ensure far-offset coverage over and around the entire CO₂ layer. The dense shot carpet aims for fine spatial sampling in common-offset domain, to provide a detailed assessment of the diving-wave delay. We use a zero phase, 16 Hz peak frequency source, with a temporal length of about 400 ms, as shown in Fig. 13. At these frequencies, the diving waves have enough signal-to-noise ratio to enable reliable delay measurements. Modelling is done with a reflecting free surface to incorporate ghosts and multiples as in the field.

Figure 13.

Geometry and source wavelet for 3-D acoustic finite-difference modelling on synthetic Sleipner scenario: (a) shot carpet over CO₂ layer. (b) Source wavelet. (c) Source amplitude spectrum. The peak frequency is 16 Hz. Shot spacing is 50 m in x and y.

Fig. 14 shows a raw shot gather from the synthetic baseline survey with a comparison against the monitor for the far-offsets. A large discrepancy between baseline and monitor can be observed in the post-critical region between 4500 and 5000 m. This difference is interpreted as the onset of interaction between the diving waves and the CO₂. Here, the delay is likely dominated by double-leg interactions and the strong interference with the CO₂ post-critical reflection. As offset increases the discrepancy between base and monitor reduces. The time-shift reduction becomes evident for offsets beyond 5000 m. This is interpreted as the onset of single-leg interactions. Additional delay reduction can be attributed to a reduced interference between the CO₂ post-critical reflection, and the natural behaviour of delay reduction with offset.

Figure 14.

Diving-wave delay on shot gather: (a) baseline shot gather with zoom-in area highlighted by the rectangle. (b) Far-offset zoom-in for baseline monitor shot gather overlay. Note the difference between the base and monitor beyond 4500 m offset. Upper and lower double-headed arrows depict the diving and guided waves, respectively.

Fig. 15 shows 6000 m common-offset sections along and perpendicular to the main axis of the CO₂ layer. The upper panel shows the time-shift calculated from the cross-correlation between the base and the monitor data. The panel shows no shift far north and south, and far east and west of the CO₂ layer. This suggests no interaction between the diving waves and the CO₂. A delay of ∼3.5 ms is measured at CMP's over the CO₂ layer. The delay increases up to ∼10 ms further north and up to ∼5 ms to the south. Negative shifts are a consequence of seismic interference and complications in the cross-correlation in the presence of diffractions. A maximum delay of about 11 ms is measured over the northern part of the ridge. The large magnitudes to the north and south are attributed to single-leg interactions, which are expected to take place at such offsets. Reliable delay measurements for offsets near the onset of interaction are complicated due to interference.

Figure 15.

Diving-wave delay on common-offset (6000 m) sections for base and monitor. The upper panel shows the cross-correlation time-shift in the section. Note the largest delay magnitudes are measured at two different locations along the main axis of the CO₂ layer.

Analytical comparison

Fig. 16(a) shows a cross-correlation time-shift map for the 6000 m common-offset cube. The single-leg interactions between the diving wave and the CO₂ become evident. Even with a 400 ms long signal, and the complex interaction with other coherent events, such as post-critical reflections, ghosts, and multiples, we can still pick the delay for the thin layer of CO₂. It is worth remembering that the minimum thickness CO₂ of 4 m is captured by the cross-correlation. Other strong events such as the guided waves, arrive at later times enabling the picking of the delay in the diving waves. These events separate further as offset increases due to differences in velocity moveout.

Figure 16.

Diving-wave time-shift maps for common-offset 6000 m: (a) delay measured by base monitor cross-correlation. (b) Delay calculated using eq. (15). The outline of the CO₂ layer is shown in black. Largest delays measured north and south of the CO₂ layer suggest single-leg interactions between the diving waves and the CO₂. Red arrows refer to profile locations in Fig. 15.

We have also calculated the delay analytically for our irregular CO₂ layer geometry. The background trend parameters are inverted from linear regression over a single trace representing the average of the baseline model. We estimate |${v}_\mathrm{ o}$| = 1767 m s⁻¹ and g = 0.39 s⁻¹ from the seabed down to 1500 m, which is the approximate penetration (turning) depth at 6 km offset. Fig. 16(b) shows the result of using eq. (15) with the inverted background parameters and the CO₂ layer models in Fig. 12. At this offset, the calculation predicts single-leg interactions between the diving waves and the CO₂. The single-leg regions mimic the shape of the CO₂ layer. The analytical time-shift and the delay at the seismic frequencies show similar spatial patterns. Similar to the modelling results using the simple 2-D models, the delay maps appear smoother for lower frequencies. Also, the magnitudes for the analytical delay are smaller. This is primarily explained by interference with the CO₂ layer post-critical reflections in the low end of the spectrum. However, other factors such as 3-D wave propagation, and different horizontal slowness between the base and monitor contribute to the discrepancy.

DISCUSSION

Diving-wave delay estimation

We have derived an analytical expression that can estimate the time delay in the diving waves due to CO₂ layers in a constant velocity gradient, isotropic background medium, at high frequencies, considering a single-transmission and the absence of a water column. We discuss ahead some, but not all, of the technical challenges related to analytically estimating the diving-wave response to CO₂ using our expression.

Eq. (15) models a single transmission, that is, only the diving-wave response to the CO₂, while our numerical examples, gather the response of the entire wavefield, including reflections, multiscattered events, etc. For seismic frequencies, and for offsets near the onset of interaction, the interference between the diving waves and the CO₂ post-critical reflection is strong, as shown in Fig. 17. Diving waves have a velocity moveout that increases with offset, which enables separation between the diving waves and the CO₂ reflection as offset increases. This onset of separation depends on the frequency, hence making the measured delay frequency dependent. Beyond the onset of separation, the delay magnitudes approach the analytical time-shifts. Interference with other coherent events such as reverberations inside the CO₂ layer, ghosts and other forms of multiples play a minor role, as seen in the simple modelling examples in Figs 7 and 9.

Figure 17.

Band-limited diving-wave delay for 1-D CO₂ layer, as shown in Fig. 1. Left: 20 Hz shot gather highlighting strong interference between the diving waves and the CO₂ reflection indicated by the double-headed arrow. Right: seismic data and cross-correlation time-shift at 5 Hz (solid red) and 20 Hz (solid blue). The black arrow highlights the CO₂ post-critical reflection. Dashed red and dashed blue curves represent the anisotropic delay at 5 and 20 Hz, respectively. The various delays are similar to each other beyond the onset of separation.

The analytical expression has been derived assuming an isotropic background medium, while some degree of anisotropy is expected in real subsurface conditions. Fig. 17 shows the effect of adding a constant anisotropy (2 per cent delta and 4 per cent epsilon) to the simple 1-D model in Fig. 1. The anisotropy changes the traveltimes and the offset for the onset of interactions. Close to this offset, the differences between the isotropic and anisotropic delays are large. At 6 km offset, the discrepancy reduces to ∼3.6 and ∼2 per cent, for 5 and 20 Hz, respectively. Overall, the analytical approximation is more accurate for large offsets where frequency and anisotropy effects on the delay are smaller.

The presence of a water column is expected to be only a minor source of error in the delay estimation for shallow water environments such as the North Sea. In the examples assessed, with an 82 m water column, such as that in the Sleipner field, the water layer travel path differences between base and monitor are negligible, and its contribution to the delay is a second-order effect. However, for deeper water environments this may become a key source of discrepancy.

Important sources of variability between the numerical modelling and the analytical calculation in Fig. 16 are: changes in the velocity gradient with depth, 3-D wave propagation and CO₂ layer geometry, and time-shift picking challenges. In general, lithological variations, depth-varying compaction regimes (Walderhaug 1996; Lander & Walderhaug 1999), and overpressure (Marcussen et al. 2009), lead to a large velocity heterogeneity that can be hard to simplify with a constant gradient trend. The linear trend parameters |${v}_\mathrm{ o}$| and g, control the penetration depth of the diving waves, and the offsets at which the various interaction patterns occur. Explaining the delay for all offsets using a single gradient poses thus a limitation. For a given offset, the delay calculation may be improved by iterative calculation of |${v}_\mathrm{ o}$| and g via linear regression on the velocities down to the expected penetration depth. Our approximation assumes a 2-D CO₂ model, whereas waves propagate in the 3-D space and interact with 3-D CO₂ geometries. Finally, picking the delay in the presence of diffractions via base/monitor cross-correlation is also challenging and an important source of discrepancy for the comparison in Fig. 16. Manual picking may be an alternative, but it may overestimate the delay where strong interference takes place.

The analytical expression may thus be used for assessing CO₂ thin layer detection limits using diving waves. It may support survey planning for CO₂ migration monitoring by enabling quick assessment of the extent of the shot carpet and offsets required to capture the various interaction patterns between the diving waves and the CO₂. For multilayer configurations, where large delays are expected due to contributions from the various CO₂ layers, the expression may be used as a predictive cycle skipping tool, providing a lower bound for the time-shift magnitude for a given CO₂ configuration.

Diving-wave time-lapse analysis

Transmission time-lapse analyses using data modelled using realistic representations of the subsurface at Sleipner, showed that it is possible to capture the delay in the diving waves due to thin layers of CO₂ despite the complexity of the subsurface. Such long-offset analysis is not possible to perform using current field data due to differences in streamer length between the base and monitor. We next discuss our diving-wave analysis in relation to the field data challenges.

Our acoustic VTI finite-difference modelling using the single-layer 2010 Sleipner model captures the complexity of field data, particularly the strong interference with other coherent events such as ghosts, multiples, guided waves, post-critical reflections and diffractions. We assess the delay at far offsets where the diving waves become kinematically separated from high-amplitude events such as guided waves. Therefore, inaccuracies in the amplitude modelling of the latter using the acoustic formulation are not expected to influence the results. Viscous effects, such as anelastic attenuation and dispersive shift due to the waves travelling through the CO₂ layer are expected to occur in the field (Furre et al. 2017; Papageorgiou & Chapman 2021). However, this generally affects the high frequencies. In our example, and in general, diving waves are analysed at low frequencies < 20 Hz, which are only mildly affected by attenuation. Elastic effects, such as P- to S-wave conversions (Agudo et al. 2018), are expected to result in an additional amplitude loss in the transmitted wavefield for field data. However, the latter does not influence the kinematics of the P-diving wave used in the time-lapse analysis.

Injecting CO₂ into an aquifer may also introduce time-lapse effects in the overburden and below the reservoir. Aquifer dilatation and overburden compaction have been documented in the In Salah CO₂ storage complex, in Algeria, where reservoir expansion led to surface uplift (Ringrose et al. 2013; White et al. 2014). These effects are not modelled in our time-lapse analysis and their impact on detectability have not been studied. At Sleipner, the pressure rise is minimal due to high injectivity and large storage capacity (Eiken et al. 2011), and injection-induced time-lapse effects are negligible. However, these geomechanical effects may influence CO₂ detectability in low injectivity settings, where pressure rise and overburden compaction may counteract the delay in the diving waves.

Repeatability is also an important aspect of detectability and is considered to play an important role in transmission time-lapse analysis (Landrø 2015). In our study, baseline source and receiver positions are replicated for the monitor survey, while some discrepancy is expected for field data. A 2 m thin CO₂ layer yields a diving-wave delay of 2–5 ms, as shown in Fig. 18, with additional variability expected depending upon aquifer stiffness. Picking the time-shift on upsampled low-frequency data makes it possible to detect anomalies of a fraction of a millisecond, and thus sense CO₂ layers of thicknesses < 1 m. For field data, the detection limits are higher and are bounded by the levels of non-repeatable noise.

Figure 18.

5 Hz diving-wave delay as a function of offset for a 1-D CO₂ layer of varying thickness: 2 m (red), 4 m (blue), 10 m (green) and 20 m (yellow). Base model is the simple trend in Fig. 4(b). A 2 m CO₂ layer yields a low-frequency delay of 2–5 ms.

Acquiring seismic records that capture diving waves interacting with the CO₂ plumes may require increasing the length of the streamer for traditional marine tow-streamer acquisitions. This may lead to longer turning times in the acquisition racetrack patterns and to some increase in survey costs. Nevertheless, the use of sparse node systems and distributed acoustic sensing fibre optic cables is on the rise as effective low-cost solutions for CO₂ monitoring. These technologies enable the acquisition of long-offset seismic records and favour the application of diving-wave time-lapse analyses.

The method presented here offers an alternative to reflection time-lapse seismic, particularly for thin layers of CO₂, where tuning interference limits the resolving power of reflections. The diving waves were able to sense down to 4 m thin CO₂ layers, which is the minimum CO₂ thickness in our model, despite the limited source bandwidth and the long time span of the source. Direct picking of the signal in the common-offset volumes offers limited knowledge about the depth of the anomaly. Also, for far offsets, where single-leg interactions dominate, the anomalies may not correlate spatially with the location of the CO₂ layer in the seismic section. Advanced waveform imaging techniques such as FWI offer a scheme to exploit the diving-wave delay over the full bandwidth, and to locate the time-lapse anomalies at their true subsurface position (Landrø 2015). Its application combining reflections and transmissions for CO₂ thin layer detection is a subject for future research.

CONCLUSIONS

We have derived an expression that can approximate the delay in diving waves due to thin layers of CO₂, as a function of offset, for various CO₂ layer conditions. The expression is valid for high frequencies and can handle some degree of lateral variations in CO₂ properties, double- and single-leg interactions between diving waves and the CO₂, and multilayer CO₂ geometries. For seismic frequencies, the expression is only valid for far offsets, where the diving waves become kinematically separated from the CO₂ post-critical reflection due to moveout differences. This expression can be used to assess CO₂ detection limits using diving waves, provide insights into survey planning for CO₂ migration monitoring, and to conduct FWI cycle skipping analysis. Diving-wave time-lapse analyses offer a complementary method to reflection time-lapse analyses for CO₂ migration monitoring and CO₂ thin layer detection. Numerical modelling using our realistic representation of the single-layer 2010 Sleipner scenario showed that the diving waves can capture the delay due to a thin and complex CO₂ layer, sensing CO₂ thicknesses down to 4 m, despite limited bandwidth. At far offsets, the diving waves separate from other coherent events such as the guided waves, making it possible to measure their delay due to the interaction with the CO₂. The main limitations of mapping the delay directly from the seismic is that, at far offsets, where single-leg interactions dominate, the anomalies do not conform spatially with the CO₂ layer, and in general, it is challenging to precisely locate the depth of the anomaly. Advanced imaging techniques such as FWI, can use the diving-wave delay for all offsets, and over the full bandwidth, to locate the thin CO₂ layer response at its true subsurface position and its application is subject to future research.

Acknowledgement

This work was conducted as part of the Carbon Capture and Storage management program of the Centre for Geophysical Forecasting (CGF, https://www.ntnu.edu/cgf/). We would like to thank CGG for waveform modelling and computing resources, AkerBP and NTNU for seismic waveform modelling, Equinor and CO₂ Datashare Group for access to post-stack time-lapse surveys and velocity trends provided by the Sleipner license partners Equinor Energy AS, Vår Energy ASA, PGNiG Upstream Norway AS, KUFPEC Norway AS. Peng Zhao from CGG is thanked for valuable discussions.

DATA AVAILABILITY

The seismic post-stack data were derived from sources in the public domain: [https://co2datashare.org/dataset]. The velocity model is confidential and sharing is only possible upon approval from Equinor and the Sleipner License Partners. All remaining synthetic data produced in this study can be shared upon reasonable request to the corresponding author.

CONFLICT OF INTEREST

The authors acknowledge that there are no conflicts of interest recorded.

AUTHOR CONTRIBUTIONS

Ricardo Martinez: Formal analysis, investigation, methodology, project administration, writing – original draft, writing - review and editing. Vetle Vinje: Methodology, project administration, resources, supervision. Alexey Stovas: Conceptualization, formal analysis, supervision. Joachim Mispel: Supervision, methodology. Philip Ringrose: Conceptualization, supervision, writing - review and editing. Kenneth Duffaut: Supervision, methodology. Martin Landrø: Conceptualization, funding acquisition.

FUNDING

The work was funded by the Centre for Geophysical Forecasting (grant no. 309960).