TimeLapse FullWaveform Inversion Using CrossCorrelationBased Dynamic Time Warping
https://doi.org/10.1007/s11804024004403

Abstract
Offshore carbon capture, utilization, and storage (OCCUS) is regarded as a crucial technology for mitigating greenhouse gas emissions. Quantitative monitoring maps of sealed carbon dioxide are necessary in a comprehensive OCCUS project. A potential highresolution method for the aforementioned purpose lies in the fullwaveform inversion (FWI) of timelapse seismic data. However, practical applications of FWI are severely restricted by the wellknown cycleskipping problem. A new timelapse FWI method using crosscorrelationbased dynamic time warping (CDTW) is proposed to detect changes in the subsurface property due to carbon dioxide (CO_{2}) injection and address the aforementioned issue. The proposed method, namely CDTW, which combines the advantages of crosscorrelation and dynamic time warping, is employed in the automatic estimation of the discrepancy between the seismic signals simulated using the baseline/initial model and those acquired. The proposed FWI method can then backproject the estimated discrepancy to the subsurface space domain, thereby facilitating retrieval of the induced subsurface property change by taking the difference between the inverted baseline and monitor models. Numerical results on pairs of signals prove that CDTW can obtain reliable shifts under amplitude modulation and noise contamination conditions. The performance of CDTW substantially outperforms that of the conventional dynamic time warping method. The proposed timelapse fullwaveform inversion (FWI) method is applied to the Frio2 CO_{2} storage model. The baseline and monitor models are inverted from the corresponding timelapse seismic data. The changes in velocity due to CO_{2} injection are reconstructed by the difference between the baseline and the monitor models.Article Highlights● Considering the cycleskipping problem is crucial to the effective application of fullwaveform inversion in an OCCUS project.● As a method for quantifying the difference in seismic waveforms induced by CO_{2} injection, crosscorrelationbased dynamic time warping is investigated.● The cycleskipping problem can be mitigated by the CDTWFWI method using the traveltime information.● The performance of CDTWFWI is investigated on the Frio2 CO_{2} storage model, and a comparison is made with the conventional FWI method. 
1 Introduction
Offshore carbon capture, utilization, and storage (OCCUS) is an important technology for greenhouse gas mitigation (Vandeweijer et al., 2011; Hansen et al., 2013; Furre et al., 2017; Ringrose and Meckel, 2019; Flohr et al., 2021; Roche et al., 2021). CO_{2} mitigation can be realized through OCCUS technology, which is a solution and technological system developed by coastal countries or regions through engineering methods. Compared to onshore CO_{2} storage technology, OCCUS technology offers several advantages, such as extensive storage sinks and enhanced security (Meckel et al., 2017; Eide et al., 2019; Suicmez, 2019; Li, 2022). CO_{2} injection changes the elastic parameters of the subsurface media, and therefore the seismic response. This makes seismic methods are crucial tools for monitoring CO_{2} injection. Seismic methods offer the potential to collect comprehensive data on elastic parameters and reservoir properties, including lithologies, porosity, and the state of pore fluids. These data can supply detailed information that is valuable for understanding the subsurface environment during CO_{2} injection (Bosch et al., 2010). Seismic methods have been extensively applied in various CO_{2} storage projects during the past two decades to monitor the migration of subsurface CO_{2} plume (Arts et al., 2004; Daley et al., 2007; Saito et al., 2008; Harris et al., 2016; Wang and Morozov, 2020; Li et al., 2023). The migration of CO_{2} plume is successfully identified through migrated images, velocity differences, or reservoir parameters using these applications.
With the advancement of OCCUS, a recognized need for the quantitative estimation of CO_{2} and detailed monitoring of its flow path emerges. Among various seismic methods available, fullwaveform inversion (FWI) is continuously gaining increasing attention due to its potential to realize the highest resolution (Tarantola, 1984; Virieux and Operto, 2009; Zhang et al., 2019). Recent studies (Queißer and Singh, 2013; Li et al., 2021; Hu et al., 2023) have shown that timelapse FWI of repeated seismic surveys can accurately estimate highresolution variations related to CO_{2} injection. A direct relationship is observed between these parameters and the reservoir properties of the CO_{2} storage site. Thus, timelapse FWI is considered a valuable tool for quantitative analysis and monitoring of CO_{2} flow in underground storage sites. Regardless of its success, FWI encounters challenges in practical OCCUS projects due to the cycleskipping problem associated with the conventional pointtopoint minustype objective function. The accuracy and efficiency of the inversion process can be restricted by this issue (Tarantola, 1984; Virieux and Operto, 2009). Furthermore, the difficulties encountered by FWI in practical OCCUS projects are further exacerbated by practical issues, such as limited aperture, absence of lowfrequency data, unknown seismic wavelets, and inaccurate initial models (Virieux and Operto, 2009).
Researchers have introduced different variants of FWI to address the challenges associated with the conventional approach (Bunks et al., 1995; Sirgue and Pratt, 2004; Djebbi and Alkhalifah, 2014; Métivier et al., 2018). One basic idea involves using portions or subsets of the full information contained in seismic data. For instance, the frequency continuation strategy can mitigate the cycleskipping problem by using the lowfrequency data firstly and feed in the highfrequency data successively (Bunks et al., 1995). Travel time is an important piece of information that supports seismic data analysis. This information represents the time taken for seismic waves to propagate from the sources to the receivers. The incorporation of traveltime information into the FWI frame facilitates the inversion of the elastic parameters, which are called skeletonized FWI or traveltimebased FWI in the literature (Luo and Schuster, 1991; Tromp et al., 2005; Ma and Hale, 2013; Lu et al., 2017). Traveltimebased FWI is extensively applied in geoscience applications due to the quasilinear behavior of its objective function and reduced sensitivity to noise (Dahlen et al., 2000; Montelli et al., 2004; Pyun et al., 2005).
An essential component of traveltimebased FWI is the determination of traveltime shifts between the simulated seismic signals and the observed seismic signals collected before/after the injection. The traveltime shifts expose the disparity between the compared signals, with larger shifts indicating greater differences between the current baseline/initial models and the postinjection models. Crosscorrelation and dynamic time warping (DTW) are two commonly used methods for traveltime shift estimation (Knapp and Carter, 1976; Hale, 2013). The small shift limitation of crosscorrelation and the instability and high sensitivity to noise of DTW seriously plague their applications in complex situations (Wang et al., 2023). Crosscorrelationbased dynamic time warping (CDTW) is a leading technique for estimating traveltime shifts. The CDTW method combines the benefits of the commonly used crosscorrelation and DTW methods (Wang et al., 2023), enabling automatic estimation of timevarying shifts with robustness.
A novel timelapse FWI method, in which CDTW is used to measure the discrepancy between simulated and observed seismic signals of baseline and monitor seismic data, is developed in the current study. The abbreviation timelapse CDTWFWI is adopted for the proposed method to enhance clarity. In the CDTWFWI method, the objective function is constructed using the L2 norm of the traveltime shifts estimated by CDTW. Limitedmemory Broyden–Fletcher– Goldfarb–Shanno (LBFGS) optimization technology is used to address the inverse problem (Nocedal and Wright, 2006). The traveltime shifts can be backprojected into the subsurface space domain through the proposed CDTWFWI method, and the images of the media can then be provided. For practical implementation, the adoption of the adjointstate method facilitates the efficient computation of the objective function gradient (Plessix, 2006). This study focuses on pressure seismic waves and aims to reconstruct the pressurewave velocity model.
The remainder of this paper is organized as follows. First, the theory of the CDTWFWI method is present in the methodology section. The CDTW and conventional DTW methods for estimating the traveltime shifts on pairs of signals are compared. The performance of the proposed CDTWFWI method for reconstructing velocity variations induced by CO_{2} injection using the Frio2 CO_{2} storage model is demonstrated in the "Numerical experiment on the Frio2 CO_{2} storage model" section. Finally, conclusions regarding the proposed method are drawn.
2 Methodology
The theory of the proposed CDTWFWI method, including the governing equation of seismic waves, the inverse problem of CDTWFWI, and the review of CDTW technology, is presented in this section.
2.1 Governing equation of seismic waves
Complex phenomena of seismic waves, such as reflection, scattering, wavemode conversion, and attenuation, can be attributed to the wave equation. This study focuses on the pressure waves and aims to reconstruct the pressurewave velocity. Thus, the constantdensity acoustic equation is assumed to govern the seismic waves:
$$ \frac{1}{v^2(\boldsymbol{x})} \frac{\partial^2 p(\boldsymbol{x}, t)}{\partial t^2}\nabla^2 p(\boldsymbol{x}, t)=f(t) \delta\left(\boldsymbol{x}_s\right) $$ (1) where p(x, t) is the pressure wavefield, which is a function of spatial coordinate x and the time coordinate t; v(x) is the pressurewave propagation velocity; ▽^{2} is the Laplacian operator; f(t) is the source term with x_{s} as the source excitation location; and δ is the Dirac function. The secondorder temporal and tenthorder spatial finitedifference scheme is employed to solve Equation (1).
2.2 Inverse problem of CDTWFWI
The objective function of conventional FWI is defined as L2 norm of pointtopoint difference between simulated and observed waveform:
$$ E=\frac{1}{2} \sum\limits_{s, r} \int_0^T\left(u_{\text {cal }}\left(\boldsymbol{x}_s, \boldsymbol{x}_r ; t\right)u_{\text {obs }}\left(\boldsymbol{x}_s, \boldsymbol{x}_r ; t\right)\right)^2 \mathrm{~d} t $$ (2) where s and r indicate the indexes of sources and receivers, respectively; T denotes the recording length of the acquisition system; u_{cal}(x_{s}, x_{r}; t) and u_{obs}(x_{s}, x_{r}; t) indicate sim‐ ulated waveforms using Equation (1) and observed waveform, respectively. The receiver is denoted by the variable x_{r}.
Localized optimization technologies can be used to address the inverse problem in Equation (2). The relationship between the pressure wavefield and the velocity is complex. Thus, the objective function in Equation (2) vibrates severely and has several local minima (Virieux and Operto, 2009). Therefore, the cycleskipping problem contributes to the trapping tendency of conventional FWI in local minima (Chi et al., 2015). The cycleskipping problem can be successfully mitigated by traveltime informationbased FWI. The objective function is formulated as:
$$ E=\frac{1}{2} \sum\limits_{s, r} \int_0^T \Delta t^2\left(\boldsymbol{x}_s, \boldsymbol{x}_r ; t\right) \mathrm{d} t $$ (3) where Δt (x_{s}, x_{r}; t) denotes the traveltime shift between the simulated and observed waveform.
Traveltimebased FWI mainly involves the accurate estimation of the time shifts. In this study, the time shifts estimated by CDTW (Wang et al., 2023) are incorporated into the traveltimebased FWI. The objective function of CDTWFWI is:
$$ E=\frac{1}{2} \sum\limits_{s, r} \int_0^T \Delta t_C^2\left(\boldsymbol{x}_s, \boldsymbol{x}_r ; t\right) \mathrm{d} t $$ (4) where Δt_{C}(x_{s}, x_{r}; t) indicates the time shift estimated by CDTW. The gradient of Equation (4) with respect to velocity is formulated as follows (Chi et al., 2015; Wang et al., 2023):
$$ g(\boldsymbol{x})=\frac{\partial E}{\partial v(\boldsymbol{x})}=\sum\limits_s \int_t \frac{2}{v^3(\boldsymbol{x})} \ddot{u}\left(\boldsymbol{x}, t ; \boldsymbol{x}_s\right) w\left(\boldsymbol{x}, t ; \boldsymbol{x}_s\right) \mathrm{d} t $$ (5) where the superscript two dots represent the secondorder temporal derivative; u (x, t; x_{s}) and w (x, t; x_{s}) indicate source and receiverside wavefields, respectively. The zerolag crosscorrelation between the source and receiverside wavefields is used to calculate the gradient g (x). The sourceside wavefield is calculated using Equation (1), and the receiverside wavefield is computed using the following equation:
$$ \frac{1}{v^2(\boldsymbol{x})} \frac{\partial^2 w\left(\boldsymbol{x}, t ; \boldsymbol{x}_s\right)}{\partial t^2}\nabla^2 w\left(\boldsymbol{x}, t ; \boldsymbol{x}_s\right)=f^*\left(\boldsymbol{x}_s, \boldsymbol{x}_r ; Tt\right) $$ (6) where f ^{∗}(x_{s}, x_{r}; T − t) denotes the adjointstate source (Plessix, 2006), which is formulated as follows (Chi et al., 2015; Wang et al., 2023):
$$ f\left(\boldsymbol{x}_s, \boldsymbol{x}_r, t\right)=\Delta t_C\left(\boldsymbol{x}_s, \boldsymbol{x}_r ; t\right) \frac{\dot{u}_{\text {cal }}\left(\boldsymbol{x}_s, \boldsymbol{x}_r ; t\right)}{\int \ddot{u}_{\mathrm{cal}}\left(\boldsymbol{x}_s, \boldsymbol{x}_r ; t\right) u_{\mathrm{cal}}\left(\boldsymbol{x}_s, \boldsymbol{x}_r ; t\right) \mathrm{d} t} $$ (7) The first or secondorder temporal derivatives are denoted by the superscript one or two dots, respectively. As shown in the negative t in Equation (6), the adjointstate source is backward propagated into the subsurface.
Considering the gradient, the velocity is updated by:
$$ v^{k+1}(\boldsymbol{x})=v^k(\boldsymbol{x})+\gamma^k o(g(\boldsymbol{x})) $$ (8) where k indicates the iteration number; v^{k}(x) and γ^{k} indicate the velocity and step length at the kth loop, respectively; ο is the preconditioning operator. The LBFGS technology is used in this study as the preconditioning operator. That is, LBFGS technology is employed to calculate the search direction ο (g (x)). The combination of Equations (3) – (8) facilitates the backprojection of the traveltime shifts in the data domain to the subsurface space domain using the CDTWFWI method. Therefore, the subsurface pressurewave velocity is obtained.
The estimation of traveltime shifts in the CDTWFWI method remains unclear. The estimation strategy of the shifts using the CDTW method will be introduced in the next section.
In timelapse CDTWFWI, the baseline and monitor models are sequentially inverted. First, the baseline velocity models are inverted using the seismic data obtained preinjection from an initial model. Second, starting from the inverted baseline model, the monitor velocity model is inverted using the seismic data obtained postinjection at a specific time. Finally, the velocity variations can be estimated by taking the difference between the inverted baseline and monitor models.
2.3 CDTW technology
CDTW is regarded as one of the stateoftheart technologies for traveltime shift estimation. The CDTW technology is briefly introduced and then compared with conventional DTW.
2.3.1 Brief review of DTW and CDTW methods
In the conventional DTW method, the shifts along the identified optimal path by the error matrix are used to facilitate the matching of the two input signals. The error matrix is defined as shown below (Hale, 2013):
$$ \boldsymbol{e}[i, l] \equiv(f[i]g[i+l])^2 $$ (9) where e indicates the error matrix; f and g denote the input signals. For instance, f and g indicate the simulated and observed waveforms, respectively, in the CDTWFWI method. Based on the error matrix, the time shifts are estimated using dynamic programming, which comprises two steps: accumulation and backtracking (Hale, 2013). The first step involves the computation of the distance matrix through accumulation:
$$ \boldsymbol{d}[0, l]=\boldsymbol{e}[0, l] $$ (10) $$ \begin{gathered} \boldsymbol{d}[i, l]=\boldsymbol{e}[i, l]+\min \left\{\begin{array}{l} \boldsymbol{d}[i1, l1] \\ \boldsymbol{d}[i1, l] \\ \boldsymbol{d}[i1, l+1] \end{array}\right. \\ i=1, 2, 3, \cdots, N1 \end{gathered} $$ (11) where d is the distance matrix, and N is the number of samples for the input signals. The second step involves the implementation of backtracking according to the following formulas:
$$ s[N1]=\underset{l}{\operatorname{argmin}} \boldsymbol{d}[N1, l] $$ (12) $$ \begin{aligned} s[i1] & =\underset{l \in\{s[i]1, s[i], s[i]+1\}}{\operatorname{argmin}} \boldsymbol{d}[i1, l] \\ i & =N1, N2, N3, \cdots, 1 \end{aligned} $$ (13) After the aforementioned two steps, the time shift s is automatically estimated. The pointtopoint error matrix contributes to instability and sensitivity to noise of the conventional DTW algorithm (Zhao and Itti, 2018; Wang et al., 2023), severely restricting their use in practical situations. Wang et al. (2023) introduced the CDTW method, in which the error matrix is defined by the local crosscorrelation, to address the aforementioned problem:
$$ \boldsymbol{e}[i, l] \equiv 1\frac{\int_{iK}^{i+K} f[j] \times g[j+l] \mathrm{d} j}{\sqrt{\int_{iK}^{i+K} f^2[j] \mathrm{d} j} \sqrt{\int_{iK}^{i+K} g^2[j+l] \mathrm{d} j}} $$ (14) where K denotes half the length of the local window. After computing the error matrix (Equation (16)), the timevarying shifts are estimated through the CDTW using accumulation and backtracking algorithms (Equations (10) – (15)). The computation of the error matrix in local windows via CDTW using crosscorrelation is shown in Equation (16). Hence, the combination of the error matrix constructed from crosscorrelation with the accumulation and backtracking programs from conventional DTW is realized using CDTW. An appropriate window length is necessary to capture the local structures of the signals. For additional details on CDTW technology, refer to Wang et al. (2023).
2.3.2 Comparison of conventional DTW and CDTW
The performance of DTW and CDTW is compared on pairs of signals with known time shifts. The use of detailed maps is necessary in OCCUS monitoring. The resolution of the inversion results can be enhanced by employing a highfrequency. Therefore, the accuracy of DTW and CDTW in measuring the discrepancy between the simulated signal and the observed signal obtained postinjection is analyzed using a 500 Hz Ricker wavelet.
The convolution model a (t) = W (t) ⊗ R (t) generates the original signal shown in Figure 1(a) (the red line). In this model, ⊗ is the convolution operator, and W (t) and R (t) are the wavelet and reflection coefficients, respectively. This example uses a 500 Hz Ricker function as the wavelet. Random function randn() generates the reflection coefficient in MATLAB. Known shifts (Figure 1(b)) are added to the original signal to obtain the perturbed signal (Figure 1(a), the blue line). The travel shifts induced by CO_{2} injection beneath the ocean bottom in OCCUS or other scenes can be simulated using this pair of signals. Distance matrices computed by DTW and CDTW in this noisefree situation exhibit clear optimal paths along which accurate time shifts are estimated (Figure 2).
DTW and CDTW provide accurate time shifts in the aforementioned example. However, the field signals display amplitude modulation and noise contamination, which will aggravate the challenges in shift estimation. The original signal (red line) and the amplitudemodulated version (blue line) are shown in Figure 3(a). Multiplying a timevarying factor (Figure 3(b)) to the perturbed signal (Figure 1(a), the blue line) facilitates the generation of the modulated signal. In this case, a clear optimal path is absent on the plot of the distance matrix by conventional DTW (Figure 4(a)). The estimated time shifts by conventional DTW degrade at times around 0.05 and 0.18 s. On the contrary, the amplitude modulation does not affect the CDTW method (Figure 4(b)). The computed distance matrix and estimated time shifts by CDTW are comparable to those in the noisefree case (Figure 2(b)).
The performance of conventional DTW and CDTW technologies in the noise contamination case is further compared. The randn() function generates random noise in MATLAB. The random noise is then added to the perturbed signal (Figure 1(a), the blue line). The signal/noise ratio of the noisecontaminated signal (Figure 5(a), the blue line) is 2 dB. The conventional DTW technology severely hinders the noise. No visible optimal path is observed on the distance matrix by conventional DTW. The estimated time shifts by DTW are far away from the actual one in the time range from 0 s to 0.13 s (Figure 5(b), the white circles). Notably, the time shifts reach the artificially set boundary (−20 ms), indicating that conventional DTW unsuccessfully performs in the noise contamination case. By contrast, the noise influence on the CDTW technology is minimal. The estimated time shifts by CDTW fit well with the true one (Figure 5(c), the white circles).
The performance of conventional DTW and CDTW technologies in the noisefree, amplitudemodulated, and noisecontaminated cases is compared in this section using highfrequency seismic signals. The presence of amplitude modulation and noise severely affects and contributes to the unstable conventional DTW, thus hindering its use in practical OCCUS applications. By contrast, the CDTW technology can resist amplitude modulation and noise presence. This technology demonstrates promising potential in practical OCCUS applications. The proposed CDTWFWI will be tested in the next section using numerical experiments.
3 Numerical experiment on the Frio2 CO_{2} storage model
The sources are typically placed near the ocean surface during seismic surveys for monitoring in OCCUS. Meanwhile the receivers are either placed near the ocean surface or down to the ocean bottom. In both cases, the seismic waves propagate through and are affected by the seawater volume. Data processing and quantitative monitoring become increasingly difficult due to this phenomenon. For instance, the images and inversion results are severely degraded by the freesurface multiples (Berkhout and Verschuur, 1997). An alternative to the aforementioned surface sources and receivers is downwell sources and receivers, which can mitigate the seawater influence and are closer to the monitoring target (Daley et al., 2011). Based on previous work (Hovorka et al., 2006; Daley et al., 2007), borehole seismic techniques can image the in situ change in seismic velocity due to the displacement of brine by supercritical CO_{2}. The effectiveness of the proposed CDTWFWI method is tested in this section using a crosswell acquisition configuration, in which sources are deployed in a well and receivers in another well far away at an appropriate distance, respectively.
The pressurewave velocity model is modified upon the Frio2 CO_{2} pilot experiment (Daley et al., 2007). In this smallscale injection experiment, approximately 300 t of supercritical CO_{2} are injected into the 17 mthick Blue sand reservoir at a depth of 1 650 m at the Frio site. The original horizontal and vertical grid intervals are 0.15 m. The original model is discretized on 467×434 grids. The original horizontal and vertical grid intervals are enlarged by an order of magnitude without varying the number of grids (i.e., the grid interval of the modified model is 1.5 m) (Figure (6)).
In the crosswell survey, the source and receiver wells are assumed to be located at x = 15 m and x = 685.5 m, respectively (Figure 6(a)). A total of 109 sources and 217 receivers are evenly deployed in the source and receiver wells, respectively, to record the seismic data. A Ricker wavelet, which has a dominant frequency of 250 Hz, is used as the source term to generate the seismic "observed" data using the pre and postinjection velocity model (Figure 6). The velocity inside the reservoir varies from 2 650 m/s to 2 765 m/s and is equal to 2 700 m/s outside the reservoir.
In the first stage, the seismic data acquired preinjection are inverted for the baseline model. A homogeneous velocity model of 2 700 m/s is used in this stage as the initial model. The velocity updates by CDTWFWI are shown in Figure 7(a). The updates fit reasonably well with the baseline model. Conventional FWI is conducted using the CDTWFWI model to obtain highresolution baseline velocity. The frequency band from 35 Hz to 500 Hz is used in FWI by a frequency continuation strategy (Bunks et al., 1995). Figure 7(b) shows that the inverted baseline velocity model has high quality.
Next, the monitor seismic data for the monitor velocity model are then inverted based on the inverted baseline model.
A crucial component for updating the velocity model is known as the kernel. Figure 8(a) shows the kernel of the first iteration, where the source and the receiver are located at z=199.5 m and z=600 m, respectively. The kernel is superimposed on Figure 6(a). The sensing area of the kernel is clearly depicted in Figure 8(a). This figure highlights the region illuminated by the kernel. Additionally, the travel time of the source–receiver pair will be affected by variations in velocity within the corresponding aquifer area. This finding indicates the feasibility of designing a costeffective acquisition system with a limited number of sources or receivers. Related work (Daley et al., 2007) has demonstrated that valuable insights into the migration of a CO_{2} plume can be obtained using a carefully arranged downhole source. Figure 8(b) displays the final gradient after summation over 109 sources in the initial iteration. This final gradient effectively emphasizes the lowvelocity area induced by CO_{2} injection.
The reconstructed pressurewave velocity variations by the CDTWFWI method are shown in Figure 9. The lowvelocity area matches well with the truth (Figure 6(b)). The reconstructed result (Figure 9) reveals the high CO_{2} saturation area around x=200 m and the narrow low CO_{2}saturation area. To some degree, the velocity is underestimated because the traveltime information is sensitive to low and intermediatewavenumber velocity variations and not sensitive to the highwavenumber velocity variations. The evolution of the normalized objective function versus iteration number (Figure 10) illustrates a considerable reduction in the seismic data discrepancy.
Conventional FWI was performed using different initial models to obtain a highresolution monitor model. A frequency continuation strategy (Bunks et al., 1995) is adopted from 35 Hz to 500 Hz. Based on the inverted baseline model (Figure 7(b)), the gradient of the first loop from conventional FWI shows an incorrect updating direction (Figure 8(c)). Thus, the inverted velocity model (Figure 11(b)) is far away from the true one (Figure 11(a)), indicating trapping of FWI starting from the inverted baseline model in a local minimum. By contrast, the gradient of the first loop from conventional FWI demonstrates satisfactory updating direction using the CDTWFWI model inverted from timelapse data (Figure 9). Starting from the velocity model reconstructed by CDTWFWI (Figure 9), a highresolution velocity model is provided by conventional FWI (Figure 11(c)). Figure 11(c) shows the inverted monitor model. If low frequencies (< 35 Hz) are available, then the low wavenumbers of the final result can be further improved. Quantitative estimation of reservoir parameters such as saturation can be established by combining rock physics with the inverted monitor model. The highresolution velocity variations induced by CO_{2} injection can be attributed to the difference between the inverted baseline and monitor models (Figure 11(d)).
Finally, Figure 12 reveals the reduction in data differences induced by CO_{2} injection using the proposed CDTWFWI method. A remarkable data difference is observed between the monitor shot gathers and the shot gathers simulated using the homogeneous model (Figure 12(a)). A substantial reduction in the discrepancy between the monitor data and the data simulated using the inverted baseline model is observed after the first stage inversion (Figure 12(b)). As shown in Figure 12(c), the computed waveform using the inverted monitor is close to the true monitor data, indicating the successful estimation of velocity changes.
4 Conclusions
The timelapse CDTWFWI method was applied to facilitate the inversion of seismic data and estimation of velocity models for CO_{2} sequestration monitoring. The cycleskipping problem was effectively addressed through this method, and the velocity variations due to CO_{2} injection were accurately reconstructed. Highresolution velocity models were obtained by combining CDTWFWI with conventional FWI methods. The obtained results on pairs of signals demonstrate that CDTW can accurately measure the data discrepancy due to CO_{2} injection, even in cases of amplitude modulation and noise contamination. Furthermore, the effectiveness of the CDTWFWI method was proven through numerical experiments on the Frio2 CO_{2} storage model. Thus, future research should focus on the application of the proposed method to field OCCUS seismic data.
Competing interest The authors have no competing interests to declare that are relevant to the content of this article. 
Arts R, Eiken O, Chadwick A, Zweigel P, Van der Meer L, Zinszner B (2004) Monitoring of CO_{2} injected at Sleipner using timelapse seismic data. Energy 29(910): 13831392. https://doi.org/10.1016/j.energy.2004.03.072 Berkhout AJ, Verschuur DJ (1997) Estimation of multiple scattering by iterative inversion, Part I: Theoretical considerations. Geophysics 62(5): 15861595. https://doi.org/10.1190/1.1444261 Bosch M, Mukerji T, Gonzalez EF (2010) Seismic inversion for reservoir properties combining statistical rock physics and geostatistics: A review. Geophysics 75(5): 75A16575A176. https://doi.org/10.1190/1.3478209 Bunks C, Saleck FM, Zaleski S, Chavent G (1995) Multiscale seismic waveform inversion. Geophysics 60(5): 14571473. https://doi.org/10.1190/1.1443880 Chi B, Dong L, Liu Y (2015) Correlationbased reflection fullwaveform inversion. Geophysics 80(4): R189R202. https://doi.org/10.1190/GEO20140345.1 Daley TM, AjoFranklin JB, Doughty C (2011) Constraining the reservoir model of an injected CO_{2} plume with crosswell CASSM at the FrioⅡ brine pilot. International Journal of Greenhouse Gas Control 5(4): 10221030. https://doi.org/10.1016/j.ijggc.2011.03.002 Daley TM, Solbau RD, AjoFranklin JB, Benson SM (2007) Continuous activesource seismic monitoring of CO_{2} injection in a brine aquifer. Geophysics 72(5): A57A61. https://doi.org/10.1190/1.2754716 Dahlen FA, Hung SH, Nolet G (2000) Fréchet kernels for finitefrequency traveltimes—I. Theory. Geophysical Journal International 141(1): 157174. https://doi.org/10.1046/j.1365246X.2000.00070.x Djebbi R, Alkhalifah T (2014) Traveltime sensitivity kernels for wave equation tomography using the unwrapped phase. Geophysical Journal International 197(2): 975986. https://doi.org/10.1093/gji/ggu025 Eide LI, Batum M, Dixon T, Elamin Z, Graue A, Hagen S, Hovorka S, Nazarian B, Nokleby PH, Olsen GI, Ringrose P, Vieira RAM (2019) Enabling largescale carbon capture, utilisation, and storage (CCUS) using offshore carbon dioxide (CO_{2}) infrastructure developments—A review. Energies 12(10): 1945. https://doi.org/10.3390/en12101945 Flohr A, Schaap A, Achterberg EP, Alendal G, Arundell M, Berndt C, Blackford J, Böttner C, Borisov SM, Brown R, Bull JM, Carter L, Chen B, Dale AW, Beer DD, Dean M, Deusner C, Dewar M, Durden JM, Elsen S, Fischer JP, Gana A, Gros J, Haeckel M, Hanz R, Holtappels M, Hosking B, Huvenne VAI, James RH, Koopmans D, Kossel E, Leighton TG, Li JH, Lichtschlag A, Linke P, Loucaides S, MartinezCabanas M, Matter JM, Mesher T, Monk S, Mowlem M, Oleynik A, Papadimitriou S, Paxton D, Pearce CR, Peel K, Roche K, Ruhl HA, Saleem U, Sands C, Saw K, Schmidt M, Sommer S, Strong SJ, Triest J, Ungerböck B, Walk J, White P, Widdicombe S, Wilson RE, Wright H, Wyatt J, Connelly D (2021) Towards improved monitoring of offshore carbon storage: a realworld field experiment detecting a controlled subseafloor CO_{2} release. International Journal of Greenhouse Gas Control 106: 103237. https://doi.org/10.1016/j.ijggc.2020.103237 Furre AK, Eiken O, Alnes H, Vevatne JN, Kiær AF (2017) 20 years of monitoring CO_{2}injection at Sleipner. Energy Procedia 114: 39163926. https://doi.org/10.1016/j.egypro.2017.03.1523 Hale D (2013) Dynamic warping of seismic images. Geophysics 78(2): S105S115. https://doi.org/10.1190/geo20120327.1 Hansen O, Gilding D, Nazarian B, Osdal B, Ringrose P, Kristoffersen JB, Eiken O, Hansen H (2013) Snohvit: The history of injecting and storing 1 Mt CO_{2} in the fluvial Tubåen Fm. Energy Procedia 37: 35653573. https://doi.org/10.1016/j.egypro.2013.06.249 Harris K, White D, Melanson D, Samson C, Daley TM (2016) Feasibility of timelapse VSP monitoring at the Aquistore CO_{2} storage site using a distributed acoustic sensing system. International Journal of Greenhouse Gas Control 50: 248260. https://doi.org/10.1016/j.ijggc.2016.04.016 Hovorka SD, Benson SM, Doughty C, Freifeld BM, Sakurai S, Daley TM, Kharaka YK, Holtz MH, Trautz RC, Nance HS, Myer LR, Knauss KG (2006) Measuring permanence of CO_{2} storage in saline formations: the Frio experiment. Environmental Geosciences 13(2): 105121. https://doi.org/10.1306/eg.11210505011 Hu Q, Grana D, Innanen KA (2023) Feasibility of seismic timelapse monitoring of CO_{2} with rock physics parametrized full waveform inversion. Geophysical Journal International 233(1): 402419. https://doi.org/10.1093/gji/ggac462 Knapp CH, Carter GC (1976) The generalized correlation method for estimation of time delay. IEEE Transactions on Acoustics, Speech and Signal Processing 24: 320327. https://doi.org/10.1109/TASSP.1976.1162830 Li D, Peng S, Guo Y, Lu Y, Cui X, Du W (2023) Reservoir multiparameter prediction method based on deep learning for CO_{2} geologic storage. Geophysics 88(1): M1M15. https://doi.org/10.1190/geo20210717.1 Li D, Peng S, Huang X, Guo Y, Lu Y, Cui X (2021) Timelapse full waveform inversion based on curvelet transform: Case study of CO_{2} storage monitoring. International Journal of Greenhouse Gas Control 110: 103417. https://doi.org/10.1016/j.ijggc.2021.103417 Li J (2022) Accelerate the offshore CCUS to carbonneutral China. Fundamental Research. https://doi.org/10.1016/j.fmre.2022.10.015 Lu K, Li J, Guo B, Fu L, Schuster G (2017) Tutorial for waveequation inversion of skeletonized data. Interpretation 5(3): SO1SO10. https://doi.org/10.1190/INT20160241.1 Luo Y, Schuster GT (1991) Waveequation traveltime inversion. Geophysics 56(5): 645653. https://doi.org/10.1190/1.1443081 Ma Y, Hale D (2013) Waveequation reflection traveltime inversion with dynamic warping and fullwaveform inversion. Geophysics 78(6): R223R233. https://doi.org/10.1190/geo20130004.1 Meckel TA, Trevino R, Hovorka SD (2017) Offshore CO_{2} storage resource assessment of the northern Gulf of Mexico. Energy Procedia 114: 47284734. https://doi.org/10.1016/j.egypro.2017.03.1609 Métivier L, Allain A, Brossier R, Mérigot Q, Oudet E, Virieux J (2018) Optimal transport for mitigating cycle skipping in fullwaveform inversion: A graphspace transform approach. Geophysics 83(5): R515R540. https://doi.org/10.1190/geo20170807.1 Montelli R, Nolet G, Dahlen FA, Masters G, Engdahl ER, Hung SH (2004) Finitefrequency tomography reveals a variety of plumes in the mantle. Science 303(5656): 338343. https://doi.org/10.1126/science.1092485 Nocedal J, Wright SJ (2006) Numerical optimization. Springer, New York, USA Plessix R (2006) A review of the adjointstate method for computing the gradient of a functional with geophysical applications. Geophysical Journal International 167(2): 495503. https://doi.org/10.1111/j.1365246X.2006.02978.x Pyun S, Shin C, Min DJ, Ha T (2005) Refraction traveltime tomography using damped monochromatic wavefield. Geophysics 70(2): U1U7. https://doi.org/10.1190/L1884829 Queißer M, Singh SC (2013) Localizing CO_{2} at sleipner—Seismic images versus Pwave velocities from waveform inversion. Geophysics 78(3): B131B146. https://doi.org/10.1190/geo20120216.1 Ringrose PS, Meckel TA (2019) Maturing global CO_{2} storage resources on offshore continental margins to achieve 2DS emissions reductions. Scientific Reports 9(1): 110. https://doi.org/10.1038/s4159801954363z Roche B, Bull JM, MarinMoreno H, Leighton TG, FalconSuarez IH, Tholen M, White PR, Provenzano G, Lichtschlag A, Li J, Faggetter M (2021) Timelapse imaging of CO_{2} migration within nearsurface sediments during a controlled subseabed release experiment. International Journal of Greenhouse Gas Control 109: 103363. https://doi.org/10.1016/j.ijggc.2021.103363 Saito H, Nobuoka D, Azuma H, Xue Z (2008) Time lapse cross well seismic tomography monitoring CO_{2} geological sequestration at Nagaoka pilot project site. Journal of MMIJ 124: 7886. https://doi.org/10.2473/journalofmmij.124.78 Sirgue L, Pratt RG (2004) Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies. Geophysics 69(1): 231248. https://doi.org/10.1190/1.1649391 Suicmez VS (2019) Feasibility study for carbon capture utilization and storage (CCUS) in the Danish North Sea. Journal of Natural Gas Science and Engineering 68: 102924. https://doi.org/10.1016/j.jngse.2019.102924 Tarantola A (1984) Inversion of seismic reflection data in the acoustic approximation. Geophysics 49(8): 12591266. https://doi.org/10.1190/1.1441754 Tromp J, Tape C, Liu Q (2005) Seismic tomography, adjoint methods, time reversal and bananadoughnut kernels. Geophysical Journal International 160(1): 195216. https://doi.org/10.1111/j.1365246X.2004.02453.x Vandeweijer V, Van Der Meer B, Hofstee C, Mulders F, D'Hoore D, Graven H (2011) Monitoring the CO_{2} injection site: K12B. Energy Procedia 4: 54715478. https://doi.org/10.1016/j.egypro.2011.02.532 Virieux J, Operto S (2009) An overview of fullwaveform inversion in exploration geophysics. Geophysics 74(6): WCC1WCC26. https://doi.org/10.1190/L3238367 Wang J, Dong L, Huang C, Wang Y (2023) Crosscorrelationbased dynamic time warping and its application in wave equation reflection traveltime inversion. Geophysics 88(6): R737R749. https://doi.org/10.1190/geo20230089.1 Wang Y, Morozov IB (2020) Timelapse acoustic impedance variations during CO_{2} injection in Weyburn oilfield, Canada. Geophysics 85(1): M1M13. https://doi.org/10.1190/geo20190221.1 Zhang Z, Alkhalifah T, Wu Z, Liu Y, He B, Oh J (2019) Normalized nonzerolag crosscorrelation elastic fullwaveform inversion. Geophysics 84(1): R1R10. https://doi.org/10.1190/geo20180082.1 Zhao J, Itti L (2018) ShapeDTW: Shape dynamic time warping. Pattern Recognition 74: 171184. https://doi.org/10.1016/j.patcog.2017.09.020