TY - JOUR
AU - Lumley,, D.
AB - Abstract Seismic monitoring provides valuable information regarding the time-varying changes in subsurface physical properties caused by natural or man-made processes. However, the resulting changes in the earth's subsurface properties are often small both in terms of magnitude and spatial extent, leading to minimal time-lapse differences in seismic amplitude or traveltime. In order to better extract information from the time-lapse data, we show that exploiting the full seismic waveform information can be critical. In this study, we develop and test methods of full waveform inversion that estimate an optimal subsurface model of time-varying elastic properties in order to fit the observed time-lapse seismic data with predicted waveforms based on numerical solutions of the wave equation. Time-lapse full waveform inversion is nonlinear and non-unique, and depends on the knowledge of the baseline velocity model before a change, and (non-)repeatability of earthquake source and sensor parameters, and of ambient and cultural noise. We propose to use repeating earthquake data sets acquired with permanent arrays of seismic sensors to enhance the repeatability of source and sensor parameters. We further develop and test time-lapse parallel, double-difference and bootstrapping inversion strategies to mitigate the dependence on the baseline velocity model. The parallel approach uses a time-invariant full waveform inversion method to estimate velocity models independently of the different source event times. The double-difference approach directly estimates velocity changes from time-lapse waveform differences, requiring excellent repeatability. The bootstrapping approach inverts for velocity models sequentially in time, implicitly constraining the time-lapse inversions, while relaxing an explicit requirement for high data repeatability. We assume that prior to the time-lapse inversion, we can estimate the true source locations and the origin time of the events, and also we can also obtain a reasonably accurate baseline velocity model. Extensive synthetic tests using a realistic velocity model developed from a real project area demonstrate the potential of full waveform inversion to estimate velocity changes from dense surface arrays of seismic stations recording a small number of repeating events. Analysis of sensitivity kernels suggests that positioning sensors at large distances allows for a stable recovery of velocity changes near the event locations by illuminating the inversion area with a wide aperture of angles. We show how full waveform inversion maps the errors in the baseline velocity model and the non-repeatable noise into the estimates of time-lapse velocity changes. Among the three time-lapse inversion methods, parallel inversion is most affected by non-repeatability factors, and is thus the least robust and most contaminated by artefacts. In contrast, the double-difference and bootstrapping methods result in more accurate time-lapse inversions. As the non-repeatability of both sources and noise increases, the bootstrapping method provides more robust and accurate results than the double-difference method. Inverse theory, Waveform inversion, Body waves, Seismic tomography, Wave propagation 1 INTRODUCTION Elastic properties of the earth change over time, due to natural processes such as tectonic stress loading and release (e.g. Schaff & Beroza 2004; Nakata & Snieder 2012) and volcanic activities (e.g. Brenguier et al. 2014), or due to anthropogenic activities such as fluid injection or extraction (e.g. Lumley 2001; Zoback 2010; Maxwell et al. 2010), geothermal energy (e.g. Gunasekera et al. 2003; Martínez-Garzón et al. 2013) and CO2 geosequestration (e.g. Arts et al. 2004; Daley et al. 2007). Associated changes in the stress/strain state often cause (microseismic) earthquakes. Thus, passive seismic monitoring has substantially contributed to understanding and evaluating time-lapse changes of subsurface properties (e.g. Schaff & Beroza 2004; Maxwell et al. 2010; Hotovec-Ellis et al. 2015; Ben-Zion et al. 2015). Conventionally, most subsurface changes have been qualitatively interpreted from the spatial distribution and source parameters of the seismic events (Maxwell et al. 2010; Baltay et al. 2011; Nishikawa & Ide 2015). However, because elastic properties are quantitatively related to material and fluid properties via rock physics and geomechanical relationships (Mavko et al. 1998; Saffer & Tobin 2011; Tsuji et al. 2014), accurate quantitative estimation of time-lapse changes in elastic properties supplements and significantly enhances understanding of the in situ changes (Gunasekera et al. 2003; Schaff & Beroza 2004; Peng & Ben-Zion 2005; Sawazaki et al. 2009; Brenguier et al. 2014; Hotovec-Ellis et al. 2015; Nakata et al. 2015). Table 1. List of mathematical symbols used in this paper. Superscripts (0) and (1) indicate that the variable is associated with recording time T(0) and T(1) respectively. Subscripts p, d, b indicate that the variable is associated with parallel, double-difference, and bootstrapping inversions. Subscripttrue indicate that the model is true. ‖ · ‖ indicates the L2 norm. Symbol . Description . T(0), T(1) Calendar recording time m Model m(0), m(1) Model at T(0) and T(1) mbg Time-invariant (background) component of model Δm(0), Δm(1) Time-varying component of models m(0) and m(1) Δm(1, 0) Time-lapse model change between T(0) and T(1) given by m(1) − m(0) mo Initial model δm Time invariant model error given by m − mtrue δΔm(1, 0) Time-lapse model change error given by |$\Delta \mathbf {m}^{(1,0)} - \Delta \mathbf {m}^{(1,0)}_{\rm true}$| δtl Component of δΔm(1, 0) associated with time-variant components of the models (Δm(0) and Δm(1)) δbg Component of δΔm(1, 0) associated with a time-invariant component of the models (mbg) d, u Observed and predicted data for time-invariant inversion d(0), d(1) Observed data at T(0) and T(1) u(0), u(1) Predicted data at T(0) and T(1) Δd(1, 0) Observed time-lapse data Δu(1, 0) Predicted time-lapse data δd Residual for time-invariant inversion given by u − d δd(0), δd(1) Residuals at T(0), T(1) given by δd(0) = δu(0) − δd(0), and δd(1) = δu(1) − δd(1), xs, xr Source and sensor parameters for time-invariant inversion and perfectly repeating case |$\mathbf {x}_s^{(0)},\mathbf {x}_s^{(1)}$|⁠, |$\mathbf {x}_r^{(0)}, \mathbf {x}_r^{(1)}$| Source and sensor parameters at T(0) and T(1) |$\hat{\mathbf {x}}_s$|⁠, |$\hat{\mathbf {x}}_r$| Assumed source and sensor parameters for double-difference inversion |$\mathcal {F}$| Forward modelling operator E Objective function K Sensitivity kernel G Green's function ω Real-valued angular frequency s Real-valued Laplace constant Ω, Ωl Complex-valued frequency Symbol . Description . T(0), T(1) Calendar recording time m Model m(0), m(1) Model at T(0) and T(1) mbg Time-invariant (background) component of model Δm(0), Δm(1) Time-varying component of models m(0) and m(1) Δm(1, 0) Time-lapse model change between T(0) and T(1) given by m(1) − m(0) mo Initial model δm Time invariant model error given by m − mtrue δΔm(1, 0) Time-lapse model change error given by |$\Delta \mathbf {m}^{(1,0)} - \Delta \mathbf {m}^{(1,0)}_{\rm true}$| δtl Component of δΔm(1, 0) associated with time-variant components of the models (Δm(0) and Δm(1)) δbg Component of δΔm(1, 0) associated with a time-invariant component of the models (mbg) d, u Observed and predicted data for time-invariant inversion d(0), d(1) Observed data at T(0) and T(1) u(0), u(1) Predicted data at T(0) and T(1) Δd(1, 0) Observed time-lapse data Δu(1, 0) Predicted time-lapse data δd Residual for time-invariant inversion given by u − d δd(0), δd(1) Residuals at T(0), T(1) given by δd(0) = δu(0) − δd(0), and δd(1) = δu(1) − δd(1), xs, xr Source and sensor parameters for time-invariant inversion and perfectly repeating case |$\mathbf {x}_s^{(0)},\mathbf {x}_s^{(1)}$|⁠, |$\mathbf {x}_r^{(0)}, \mathbf {x}_r^{(1)}$| Source and sensor parameters at T(0) and T(1) |$\hat{\mathbf {x}}_s$|⁠, |$\hat{\mathbf {x}}_r$| Assumed source and sensor parameters for double-difference inversion |$\mathcal {F}$| Forward modelling operator E Objective function K Sensitivity kernel G Green's function ω Real-valued angular frequency s Real-valued Laplace constant Ω, Ωl Complex-valued frequency Open in new tab Table 1. List of mathematical symbols used in this paper. Superscripts (0) and (1) indicate that the variable is associated with recording time T(0) and T(1) respectively. Subscripts p, d, b indicate that the variable is associated with parallel, double-difference, and bootstrapping inversions. Subscripttrue indicate that the model is true. ‖ · ‖ indicates the L2 norm. Symbol . Description . T(0), T(1) Calendar recording time m Model m(0), m(1) Model at T(0) and T(1) mbg Time-invariant (background) component of model Δm(0), Δm(1) Time-varying component of models m(0) and m(1) Δm(1, 0) Time-lapse model change between T(0) and T(1) given by m(1) − m(0) mo Initial model δm Time invariant model error given by m − mtrue δΔm(1, 0) Time-lapse model change error given by |$\Delta \mathbf {m}^{(1,0)} - \Delta \mathbf {m}^{(1,0)}_{\rm true}$| δtl Component of δΔm(1, 0) associated with time-variant components of the models (Δm(0) and Δm(1)) δbg Component of δΔm(1, 0) associated with a time-invariant component of the models (mbg) d, u Observed and predicted data for time-invariant inversion d(0), d(1) Observed data at T(0) and T(1) u(0), u(1) Predicted data at T(0) and T(1) Δd(1, 0) Observed time-lapse data Δu(1, 0) Predicted time-lapse data δd Residual for time-invariant inversion given by u − d δd(0), δd(1) Residuals at T(0), T(1) given by δd(0) = δu(0) − δd(0), and δd(1) = δu(1) − δd(1), xs, xr Source and sensor parameters for time-invariant inversion and perfectly repeating case |$\mathbf {x}_s^{(0)},\mathbf {x}_s^{(1)}$|⁠, |$\mathbf {x}_r^{(0)}, \mathbf {x}_r^{(1)}$| Source and sensor parameters at T(0) and T(1) |$\hat{\mathbf {x}}_s$|⁠, |$\hat{\mathbf {x}}_r$| Assumed source and sensor parameters for double-difference inversion |$\mathcal {F}$| Forward modelling operator E Objective function K Sensitivity kernel G Green's function ω Real-valued angular frequency s Real-valued Laplace constant Ω, Ωl Complex-valued frequency Symbol . Description . T(0), T(1) Calendar recording time m Model m(0), m(1) Model at T(0) and T(1) mbg Time-invariant (background) component of model Δm(0), Δm(1) Time-varying component of models m(0) and m(1) Δm(1, 0) Time-lapse model change between T(0) and T(1) given by m(1) − m(0) mo Initial model δm Time invariant model error given by m − mtrue δΔm(1, 0) Time-lapse model change error given by |$\Delta \mathbf {m}^{(1,0)} - \Delta \mathbf {m}^{(1,0)}_{\rm true}$| δtl Component of δΔm(1, 0) associated with time-variant components of the models (Δm(0) and Δm(1)) δbg Component of δΔm(1, 0) associated with a time-invariant component of the models (mbg) d, u Observed and predicted data for time-invariant inversion d(0), d(1) Observed data at T(0) and T(1) u(0), u(1) Predicted data at T(0) and T(1) Δd(1, 0) Observed time-lapse data Δu(1, 0) Predicted time-lapse data δd Residual for time-invariant inversion given by u − d δd(0), δd(1) Residuals at T(0), T(1) given by δd(0) = δu(0) − δd(0), and δd(1) = δu(1) − δd(1), xs, xr Source and sensor parameters for time-invariant inversion and perfectly repeating case |$\mathbf {x}_s^{(0)},\mathbf {x}_s^{(1)}$|⁠, |$\mathbf {x}_r^{(0)}, \mathbf {x}_r^{(1)}$| Source and sensor parameters at T(0) and T(1) |$\hat{\mathbf {x}}_s$|⁠, |$\hat{\mathbf {x}}_r$| Assumed source and sensor parameters for double-difference inversion |$\mathcal {F}$| Forward modelling operator E Objective function K Sensitivity kernel G Green's function ω Real-valued angular frequency s Real-valued Laplace constant Ω, Ωl Complex-valued frequency Open in new tab Estimating the changes in elastic properties is not a trivial exercise, and is sometimes more difficult than estimating a static velocity model, since the time-lapse changes tend to be small in terms of magnitudes and occur in a localized area especially when monitoring anthropogenic activities such as fluid injection or extraction. In addition, there are various other time-varying factors that affect seismic monitoring data including ambient and cultural noise, near-surface conditions, earthquake locations and source mechanisms, sensor locations, sensitivity and couplings. Such undesirable non-repeatability can be strong enough to obscure ‘true’ small changes in seismic data associated with the time-lapse signal within an area of interest (Lumley 2001). Therefore effects of non-repeatability on estimation of subsurface property changes need to be carefully reduced by the design of monitoring systems (e.g. using (semi-)permanent arrays as in Shulakova et al. 2014), data pre-processing (Rickett & Lumley 2001; Roach et al. 2015), and model estimation techniques. Repeating earthquakes (also called doublets or multiplets) occur at approximately the same locations with similar source mechanisms (Waldhauser & Ellsworth 2000). The high repeatability of the source parameters makes their seismic waveforms preferable for estimating changes in elastic properties (Schaff & Beroza 2004; Tkalčić et al. 2013; Hotovec-Ellis et al. 2015). Interferometric analysis of coda waves of repeating earthquakes has been extensively employed to extract small traveltime shifts in waveforms over time and then to estimate average changes in elastic parameters between the source focal point and sensor locations (Schaff & Beroza 2004; Peng & Ben-Zion 2005; Sawazaki et al. 2009; Hotovec-Ellis et al. 2015). Then an underlying tectonic or anthropogenic process is qualitatively inferred from a set of the average changes obtained for various earthquake-sensor pairs. In this study, we explore methods for obtaining a high-resolution distribution of subsurface changes directly from the main arrivals of repeating earthquake seismic data using an inversion-based wave-equation approach. Traveltime tomography has been employed to estimate changes in seismic velocity by performing a traveltime inversion procedure repeatedly over time (e.g. before and after changes in the subsurface) (Gunasekera et al. 2003), but has resulted in low resolution due to reliance on traveltime information only, and the underlying high frequency assumption of ray theory (Williamson & Worthington 1993). However, higher resolution information about subsurface changes is available by using the full frequency-dependent waveforms instead of only traveltimes, which can be achieved by using a wave-equation based approach. Full waveform inversion (Tarantola 1984) represents a set of methods to find a subsurface model that fits waveforms (not merely first arrival times) based on the numerical solution of a wave equation at finite frequencies, and is rapidly emerging as a new technique to estimate the subsurface model at a much higher spatial resolution than traveltime tomography (e.g. Zhu & Tromp 2013; Fichtner et al. 2013; Kamei et al. 2013). Full waveform inversion has been successfully applied to estimate time-invariant static elastic models at crustal and regional scales from earthquake data (Tape et al. 2009; Zhu & Tromp 2013; Fichtner et al. 2013), and in the upper crust from controlled-source data (Pratt 1999; Kamei et al. 2013; Operto et al. 2015). The applications of full waveform inversion to time-lapse changes in earth properties have been limited to controlled-source seismic monitoring data (e.g. Asnaashari et al. 2014; Raknes & Arntsen 2014), but only a handful of applications to field data sets have been reported (Queißer & Singh 2012; Zhang et al. 2013). Developing a robust time-lapse full waveform inversion method is challenging. First of all, full waveform inversion is a strongly nonlinear and ill-conditioned inverse problem, even when estimating a time-invariant subsurface model. Thus an accurate starting model is required for global convergence when using local optimization methods. To further reduce the non-uniqueness, most full waveform inversion applications adopt a multiscale approach (Bunks et al. 1995), and also focus on fitting main arrivals, thus excluding late-arriving scattered, reflected, diffracted and coda waves from the inversion. Time-lapse full waveform inversion faces additional difficulties in detecting and estimating small changes in the model, which can be easily obscured by the aforementioned data noise and other non-repeatability factors. Inversion strategies specifically designed for time-lapse full waveform inversion have been under development for controlled source time-lapse surveys (e.g. Asnaashari et al. 2014). Finally, the spatial distribution of repeating earthquakes tends to be sparse, increasing the ill-conditionedness of the inversion. In this paper, we investigate the feasibility of full waveform inversion methods to estimate time-lapse changes in velocities from repeating earthquake data recorded by dense surface arrays, which are now being increasingly deployed (e.g. Lin et al. 2013; Ben-Zion et al. 2015). We consider a case where fluid injection or withdrawal from a reservoir (e.g. of hydrocarbon, geothermal energy, groundwater, etc.) causes stress and fluid changes in the shallow part of the upper crust, and induces repeating microseismic events. We perform synthetic tests using a realistic complex velocity model based on a real reservoir and project area, and use a small number of repeating seismic events to model sparse distributions of repeating earthquakes. In this paper, we assume that the medium is 2-D acoustic and the source mechanism of the events is isotropic, and estimate the time-lapse changes in the P-wave velocities. This approach allows us to focus on evaluating the strengths and weaknesses of various time-lapse full waveform inversion strategies, and to avoid additional complications such as the non-uniqueness arising from multiparameter problems (Kamei & Pratt 2013) and complex source mechanisms. To improve the convergence of full waveform inversion, we assume a reasonably accurate initial baseline velocity model is available, and emphasize fitting main (early) P-wave arrivals. We develop and test time-lapse inversion strategies originally developed for controlled-source seismic data, and evaluate the sensitivity of full waveform inversion to the prior knowledge of the baseline (pre-change) velocity model, non-repeatable noise, and repeatability of earthquake source parameters. In the remainder of this paper, we first describe various methods of time-lapse full waveform inversion in Section 2. In Section 3, we introduce our realistic complex earth model and our inversion parameters, followed by an analysis of sensitivity kernels in Section 4, and the inversion tests using the true baseline velocity model in Section 5. In Section 6, we describe a series of inversion tests conducted using inaccurate initial velocity models. We discuss our experimental results in Section 7, and conclude in Section 8. 2 METHOD Full waveform inversion is a solution to a nonlinear inverse problem, and estimates a set of model parameters that represent the elastic properties of the earth, by fitting model-predicted waveforms to observed seismic data. The goal of the inverse problem is to minimize an objective function defined based on differences between observed and predicted waveforms. The predicted waveforms are given by the numerical solution to a wave equation, and the nonlinear objective function is optimized by local gradient methods. In this section, we first review the basic theory of full waveform inversion using local gradient methods, and then describe three inversion methods developed specifically for time-lapse applications. 2.1 Time-invariant full waveform inversion First suppose that the earth's subsurface elastic properties are invariant over time, and that seismic waveforms are excited by a set of earthquakes and recorded at an array of sensors. We would like to estimate an earth model that can generate seismic data which optimally fits the observed seismic waveforms. We discretize the elastic properties of the earth, and denote them as the model parameter vector m[x]. We also denote the observed waveforms as d[xr, xs], and the predicted waveforms as u[m; xr, xs]. We then denote a set of earthquake source parameters as |$\mathbf {x}_s=[\mathbf {x}_{s_1},\mathbf {x}_{s_2},\ldots ]$|⁠, and a set of receiver (sensor) parameters as |$\mathbf {x}_r=[\mathbf {x}_{r_1},\mathbf {x}_{r_2},\ldots ]$|⁠, where |$\mathbf {x}_{s_i}$| indicates the source parameters of the ith earthquake (e.g. location, focal mechanism etc.), and |$\mathbf {x}_{r_j}$| indicates the parameters of the jth sensor (e.g. location, coupling etc.). Note that the data can be defined in the time or frequency domains, and we drop this explicit dependence from the notation. The observed waveform data consist of signal, s, generated by a seismic event with source parameters xs, and cultural and ambient noise, n, at receivers represented by xr, \begin{equation} \mathbf {d}[\mathbf {x}_{r},\mathbf {x}_{s}] = \mathbf {s}[\mathbf {x}_r,\mathbf {x}_s] + \mathbf {n}[\mathbf {x}_r]. \end{equation} (1) Predicted (simulated, modelled) seismic data are associated with the earth model parameter m as, \begin{equation} \mathbf {u}[\mathbf {m}; \mathbf {x}_r,\mathbf {x}_s] = \mathcal {F}\lbrace \mathbf {m};\mathbf {x}_r,\mathbf {x}_{s}\rbrace , \end{equation} (2) where |$\mathcal {F}\lbrace \mathbf {m};\mathbf {x}_r,\mathbf {x}_s\rbrace$| is the forward modelling operator acting on m with source and receiver variables xr and xs and given by the numerical solution to a wave equation. The operator |$\mathcal {F}$| is nonlinear with respect to the elastic properties of the earth m. We define an objective function E from the L2 norm of the data residuals such that \begin{equation} E= \left\Vert \delta {\mathbf {d}}\right\Vert ^2, \end{equation} (3) where ‖ · ‖ is the L2 norm, and δd is the data residual, \begin{eqnarray} \delta \mathbf {d} = \mathbf {u}[\mathbf {m};\mathbf {x}_{r},\mathbf {x}_{s}] - \mathbf {d}[\mathbf {x}_{r},\mathbf {x}_{s}], \end{eqnarray} (4) \begin{eqnarray} \quad = \mathbf {u}[\mathbf {m};\mathbf {x}_{r},\mathbf {x}_{s}] - \mathbf {s}[\mathbf {x}_{r},\mathbf {x}_{s}] - \mathbf {n}[\mathbf {x}_r]. \end{eqnarray} (5) Given a starting model mo, we minimize the objective function in eq. (3) by a local gradient method: at each iteration, we update the model based on the gradient of the objective function derived by the adjoint method, and then re-compute the predicted data using the updated model to account for the nonlinearity. We then denote an optimal model as m[d; mo, xs, xr]. The initial model needs to be accurate enough to avoid cycle skipping (Sirgue & Pratt 2004), and the model estimate can also be easily contaminated by artefacts arising from factors affecting the ill-conditionedness of the inversion including data noise (n), and source-sensor geometries. We define the model error δm as \begin{equation} \delta \mathbf {m}= \mathbf {m}[\mathbf {d};\mathbf {m}_o,\mathbf {x}_s,\mathbf {x}_r] - \mathbf {m}_{\text{true}}, \end{equation} (6) such that \begin{equation} \mathbf {m}[\mathbf {d};\mathbf {m}_o,\mathbf {x}_s,\mathbf {x}_r] = \mathbf {m}_{\text{true}}+ \delta \mathbf {m}. \end{equation} (7) where δm denotes the difference between the true and estimated models (we refer to δm as the model error). Note that the accuracy of the forward modelling operator to approximate the real physics of wave propagation in the earth (modelling error) influences the inversion result (Tarantola & Valette 1982; Kamei et al. 2013), but in this paper we assume that the forward modelling operator is correct. We further assume that the medium is acoustic, and conduct the modelling and inversion for P-wave velocities only. However we have opted to keep the formulation more general, and describe our specific acoustic implementation in Appendix A. 2.2 Time-lapse full waveform inversion We now further assume that elastic properties of the earth change over time, and that we have acquired seismic records for repeating earthquakes occurring at multiple calendar times; T(k) (k = 0, 1, …). Sensor and earthquake source parameters may be different in general at each recorded time, and we denote these as |$\mathbf {x}_r^{(k)}$| and |$\mathbf {x}_s^{(k)}$| respectively. We then denote the model parameters at T(k) as m(k) that consists of the time-invariant (background) model component mbg and time-varying component Δm(k); \begin{equation} \mathbf {m}^{(k)} = \mathbf {m}^{bg}+ \Delta \mathbf {m}^{(k)}. \end{equation} (8) The time-varying Δm(k) represents both subsurface changes of our interests (e.g. velocity changes due to injection/extraction of fluid in a reservoir) and undesired changes (e.g. seasonal variations in a weathering layer, when we are interested in reservoir changes). Although the undesired model changes can be substantial (Lumley 2001), we assume that Δm(k) solely consists of the earth changes of our interests. The observed data are defined as |$\mathbf {d}^{(k)}[\mathbf {x}_r^{(k)},\mathbf {x}_s^{(k)}]$| that can be decomposed into data, s(k), and noise, n(k), as \begin{equation} \mathbf {d}^{(k)}\left[\mathbf {x}_r^{(k)},\mathbf {x}_s^{(k)}\right] = \mathbf {s}^{(k)}\left[\mathbf {x}_r^{(k)},\mathbf {x}_s^{(k)}\right] + \mathbf {n}^{(k)} \left[\mathbf {x}_r^{(k)}\right]. \end{equation} (9) Then the predicted waveforms u(k) are associated with the earth model at T(k) by \begin{equation} \mathbf {u}^{(k)} \left[\mathbf {m}^{(k)}; \mathbf {x}_r^{(k)},\mathbf {x}_s^{(k)}\right]= \mathcal {F}\lbrace \mathbf {m}^{(k)};\mathbf {x}_r^{(k)},\mathbf {x}_s^{(k)}\rbrace . \end{equation} (10) Now we focus on a set of two seismic records for earthquake sources at T(0) and T(1) to simplify discussions hereafter. However the following descriptions are readily expanded for the general case when repeating earthquakes occur more than twice. We define time-lapse changes in model parameters Δm(1, 0) between T(0) and T(1) as \begin{equation} \Delta \mathbf {m}^{(1,0)} = \mathbf {m}^{(1)} - \mathbf {m}^{(0)}. \end{equation} (11) Note that as the forward modelling operator is nonlinear, \begin{eqnarray} {\mathcal {F}\left\lbrace \mathbf {m}^{(1)}; \mathbf {x}_s^{(1)},\mathbf {x}_r^{(1)}\right\rbrace - \mathcal {F}\left\lbrace \mathbf {m}^{(0)};\mathbf {x}_s^{(0)},\mathbf {x}_r^{(0)}\right\rbrace}\nonumber\\ {\quad\ne \mathcal {F}\left\lbrace \Delta \mathbf {m}^{(1,0)};\mathbf {x}_s^{(0)},\mathbf {x}_r^{(0)}\right\rbrace ,} \end{eqnarray} (12) \begin{eqnarray} \quad\!\! \ne \mathcal {F}\left\lbrace \Delta \mathbf {m}^{(1,0)};\mathbf {x}_s^{(1)},\mathbf {x}_r^{(1)}\right\rbrace , \end{eqnarray} (13) even when |$\mathbf {x}_s^{(1)}=\mathbf {x}_s^{(0)}$| (i.e. perfectly repeating events), and |$\mathbf {x}_r^{(1)}=\mathbf {x}_r^{(0)}$| (i.e. perfect sensor repeatability). Suppose that we have some estimate of the subsurface model obtained before T(0) and denote it as mo. Then, our goal is to estimate Δm(1, 0) from the observed data, d(0) and d(1), and the initial model mo. This time-lapse inversion can be more challenging than the inversion for the time-invariant (static) model described in Section 2.1, since the true model change |$\Delta \mathbf {m}^{(1,0)}_{\rm true}$| tends to be a small percentage of the background velocity model, and is often distributed within a small spatial area (relative to a seismic wavelength). We define the error in the model change δΔm(1, 0) as \begin{equation} \delta \Delta \mathbf {m}^{(1,0)}= \Delta \mathbf {m}^{(1,0)} - \Delta \mathbf {m}^{(1,0)}_{\text{true}}. \end{equation} (14) which can be decomposed into two terms from eq. (8); \begin{equation} \delta \Delta \mathbf {m}^{(1,0)} \equiv \delta ^{tl} + \delta ^{bg}, \end{equation} (15) where δtl is the inversion error associated with errors in time-lapse components (⁠|$\delta \Delta \mathbf {m}^{(0)}=\Delta \mathbf {m}^{(0)}-\Delta \mathbf {m}^{(0)}_{\rm true}$| and |$\delta \Delta \mathbf {m}^{(1)}=\Delta \mathbf {m}^{(1)}-\Delta \mathbf {m}^{(1)}_{\rm true}$|⁠), and δbg is the inversion error associated with the error in a time-invariant (background) component (⁠|$\delta \mathbf {m}^{bg}=\mathbf {m}^{bg}-\mathbf {m}^{bg}_{\rm true}$|⁠). (Note we use ‘δ’ to indicate model errors and data residual, and ‘Δ’ to indicate time-lapse model and data.) As the given model mo is not precise in general, δbg is non-zero, and tends to span a larger area than the true time-lapse changes. Therefore the time-lapse inversion needs to reduce both δtl and δbg to very small values, in order to differentiate the desired time-lapse changes Δm(1, 0) from the errors (artefacts). Due to the non-uniqueness of the inversion, accomplishing this task is non-trivial. Time-lapse full waveform inversion tends to succeed in retrieving the time-lapse model changes δm when δtl (or Δm) is larger than δbg (δtl > δbg) before the inversion. In contrast when δtl < δbg, time-lapse full waveform inversion can become very challenging and may even fail completely. The time-lapse inversion is further complicated by non-repeatability of data especially in terms of data noise, source and receiver parameters. To illustrate this, we analyse the observed time-lapse waveform difference defined as \begin{eqnarray} \Delta \mathbf {d}^{(1,0)}\left[\mathbf {x}_r^{(1)}, \mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right] = \mathbf {d}^{(1)}\left[\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]- \mathbf {d}^{(0)} \left[\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!,\nonumber\\ \end{eqnarray}(16) and by using eq. (9) we obtain \begin{eqnarray} \Delta \mathbf {d}^{(1,0)}\left[\mathbf {x}_r^{(1)}, \mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]= \Delta \mathbf {s}^{(1,0)}\left[\mathbf {x}_r^{(1)}, \mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\nonumber\\ &&+\,\, \Delta \mathbf {n}^{(1,0)}\left[\mathbf {x}_r^{(1)}, \mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right], \end{eqnarray} (17) where Δs(1, 0) and Δn(1, 0) are the time-lapse change in signal and noise respectively given by \begin{eqnarray} \Delta \mathbf {s}^{(1,0)}\left[\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}, \mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right] & =& \mathbf {s}^{(1)}\left[\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right] - \mathbf {s}^{(0)}\left[\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right] \!, \end{eqnarray} (18) \begin{eqnarray} \Delta \mathbf {n}^{(1,0)}\left[\mathbf {x}_r^{(1)}, \mathbf {x}_r^{(0)}\right] &= & \mathbf {n}^{(1)}\left[\mathbf {x}_r^{(1)}\right] - \mathbf {n}^{(0)}\left[\mathbf {x}_r^{(0)}\right] . \end{eqnarray} (19) The time-lapse change of signal Δs(1, 0) depends solely on the time-lapse changes of the earth if the source and sensor parameters are perfectly repeatable. As Δm(1, 0) tends to be small, such ‘true’ time-lapse signal is commonly small when compared to that of the (absolute) signals, and Δs(1, 0) can be easily contaminated by the non-repeatability of source parameters (e.g. location, origin time, wavelet, focal mechanism) and sensor parameters (e.g. coupling, sensitivity, frequency response). The effects of data noise can also be significant. Since the cultural and ambient noise is not repeatable, the subtraction in eq. (19) does not cancel out non-repeatable noise. Thus the magnitude of time-lapse noise change is potentially similar to or larger than that of the small time-lapse signal change, suggesting a strong influence on the estimation of the time-lapse model change. We also note that we do not know exact earthquake source parameters, and errors in the estimated parameters can also result in errors in the estimated time-lapse model changes. Therefore in order to estimate small localized time-lapse model parameter change in a robust and stable manner, full waveform inversion methods need to accommodate issues arising from an inaccurate initial model, non-repeatable noise, non-repeatable source and receiver parameters, errors in pre-estimated source parameters, undesired subsurface changes, and forward modelling errors. In this study, we focus on the first three issues (an initial model, non-repeatable noise, and non-repeatable source and receiver parameters), and examine three methods of full waveform inversion proposed to optimally utilize time-lapse controlled-source seismic data sets (Lumley & Shragge 2013; Asnaashari et al. 2014; Raknes & Arntsen 2014). We use such methods referred to as the parallel, double-difference and bootstrapping methods, and describe and examine these methods in the rest of this section. 2.3 Parallel inversion The most straightforward way to estimate time-lapse changes in the model parameters is to run a set of two inversions in a parallel and independent manner. We invert for a model m(0) at time T(0) (or m(1) at T(1)), by starting from mo and fitting the observed data at d(0) (or d(1)). Using the L2 norm, the objective functions for parallel inversion at T(0) and T(1) are defined as \begin{eqnarray} E_p^{(0)}=\left\Vert \delta {\mathbf {d}^{(0)}_p}\right\Vert ^2, \end{eqnarray} (20) \begin{eqnarray} E_p^{(1)}=\left\Vert \delta {\mathbf {d}^{(1)}_p}\right\Vert ^2, \end{eqnarray} (21) where |$\delta \mathbf {d}^{(0)}_p$| and |$\delta \mathbf {d}_p^{(1)}$| are the data residual for parallel inversion defined as \begin{eqnarray} \delta \mathbf {d}^{(0)}_p = \mathbf {u}^{(0)}\left[\mathbf {m}^{(0)};\mathbf {x}_r^{(0)}, \mathbf {x}_s^{(0)}\right] - \mathbf {d}^{(0)}\left[\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right], \end{eqnarray} (22) \begin{eqnarray} \delta \mathbf {d}^{(1)}_p = \mathbf {u}^{(1)}\left[\mathbf {m}^{(1)};\mathbf {x}_r^{(1)}, \mathbf {x}_s^{(1)}\right] - \mathbf {d}^{(1)}\left[\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]. \end{eqnarray} (23) (Note the subscript ‘p’ refers to parallel inversion.) By minimizing the objective function |$E_p^{(0)}$|⁠, we obtain an estimated model at T(0), |$\mathbf {m}_{p}^{(0)}[\mathbf {d}^{(0)};\mathbf {m}_o,\mathbf {x}_r^{(0)}, \mathbf {x}_s^{(0)}]$|⁠. We also estimate a model at T(1), |$\mathbf {m}_{p}^{(1)}[\mathbf {d}^{(1)};\mathbf {m}_o,\mathbf {x}_r^{(1)}, \mathbf {x}_s^{(1)}]$| by optimizing |$E_p^{(1)}$|⁠. Both inversions are nonlinear inversion methods, and a local gradient method is adopted as in time-invariant full waveform inversion techniques. The velocity change Δmp is estimated by taking subtraction of the two inversion results; \begin{eqnarray} \Delta \mathbf {m}_{p}^{(1,0)}=\mathbf {m}_{p}^{(1)}\left[\mathbf {d}^{(1)};\mathbf {m}_o,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]-\mathbf {m}_{p}^{(0)}\left[\mathbf {d}^{(0)};\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right].\nonumber\\ \end{eqnarray} (24) Following eqs (14)–(15), we can rewrite eq. (24) as \begin{equation} \Delta \mathbf {m}_p^{(1,0)}=\Delta \mathbf {m}_{\text{true}}^{(1,0)}+\delta _p^{tl} + \delta _p^{bg}, \end{equation} (25) where |$\delta _p^{tl}$| is the difference in inversion errors of the time-lapse component of the models |$\mathbf {m}_p^{(0)}$| and |$\mathbf {m}_p^{(1)}$|⁠, and |$\delta _p^{bg}$| is the difference in inversion errors of their time-invariant component (see Appendix B for derivation). Since the above formulation does not require perfect repeatability, the parallel inversion approach appears advantageous in case of non-repeatable data, and one may consider that |$\delta _p^{tl}$| and |$\delta _p^{bg}$| be independent of repeatability. However, the parallel approach does not explicitly take the advantages of properties of repeating earthquakes being similar, nor takes into account the consistency of time-invariant subsurface structures mbg in m(0) and m(1). Therefore eq. (24) does not necessarily cancel out the errors in estimated time-invariant components of |$\mathbf {m}_p^{(0)}$| and |$\mathbf {m}_p^{(1)}$| that are induced e.g. by acquisition footprint, and noise. Thus |$\delta _p^{bg}\ne 0$| in general, and its magnitude increases as the non-repeatability increases and also as the accuracy of the initial model degrades. For a similar reason, the magnitude of |$\delta _p^{tl}$| also increases with the data non-repeatability. 2.4 Double-difference inversion A double-difference method was originally proposed to constrain the location of seismic events based on the traveltime differences between events (Waldhauser & Ellsworth 2000), and was then expanded to jointly estimate (absolute) time-invariant velocity models by traveltime tomography (Zhang & Thurber 2003, 2006). The approach assumes that the multiple events share the same propagation paths, and attributes any time shifts in waveforms to the combination of relative locations of the events and relative velocity anomalies near the events. Using double-difference time has led to improved event-location and velocity-anomaly estimates. The double-difference approach has been adopted also for full waveform inversion of time-lapse controlled-source data (e.g. Watanabe et al. 2004; Zhang & Huang 2013; Asnaashari et al. 2014). Double-difference full waveform inversion measures changes in waveforms of controlled sources located at the same locations repeatedly over time, and attributes waveform changes to time-lapse variations in elastic properties (mainly P-wave velocity). The approach again assumes the perfect repeatability in acquisition parameters, and the shared propagation paths outside the changed region. In this paper, we formulate the double-difference method for a general case where the source and sensor repeatability is not perfect. Then we need to assume (virtual) repeating source and sensor parameters |$\hat{\mathbf {x}}_s$| and |$\hat{\mathbf {x}}_r$| that can be related to actual source and receiver parameters via \begin{eqnarray} \hat{\mathbf {x}}_s = \mathbf {x}_s^{(0)}+\delta \mathbf {x}_s^{(0)} = \mathbf {x}_s^{(1)}+\delta \mathbf {x}_s^{(1)} , \end{eqnarray} (26) \begin{eqnarray} \hat{\mathbf {x}}_r = \mathbf {x}_r^{(0)}+\delta \mathbf {x}_r^{(0)} = \mathbf {x}_s^{(1)}+\delta \mathbf {x}_s^{(1)}, \end{eqnarray} (27) where |$\delta \mathbf {x}_s^{(0)}$| and |$\delta \mathbf {x}_s^{(1)}$| are differences between actual and assumed source parameters at T(0) and T(1), and |$\delta \mathbf {x}_r^{(0)}$| and |$\delta \mathbf {x}_r^{(1)}$| are differences between actual and assumed sensor parameters at T(0) and T(1). Then following the definition of the observed time-lapse waveforms in eq. (16), we define the time-lapse residuals as \begin{eqnarray} \delta \Delta \mathbf {d}^{(1,0)} \!=\! \Delta \mathbf {u}^{(1,0)}\left[\mathbf {m}^{(1)}, \mathbf {m}^{(0)}; \hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right]- \Delta \mathbf {d}^{(1,0)}\left[\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right],\!\!\!\!\!\!\nonumber\\ \end{eqnarray} (28) where Δu(1, 0) is the predicted time-lapse difference data, and \begin{eqnarray} \Delta \mathbf {u}^{(1,0)}\left[\mathbf {m}^{(1)}, \mathbf {m}^{(0)};\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right] = \mathbf {u}^{(1)}\left[\mathbf {m}^{(1)};\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right] - \mathbf {u}^{(0)}\left[\mathbf {m}^{(0)};\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right].\!\!\!\!\nonumber\\ \end{eqnarray} (29) We rewrite eq. (29) using Δm(1, 0) = m(1) − m(0) so that \begin{eqnarray} \Delta \mathbf {u}^{(1,0)}\left[\Delta \mathbf {m}^{(1,0)}, \mathbf {m}^{(0)};\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right] = \mathbf {u}^{(1)}\left[\mathbf {m}^{(0)}+\Delta \mathbf {m}^{(1,0)};\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right] \nonumber\\ -\, \mathbf {u}^{(0)}\left[\mathbf {m}^{(0)};\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right]. \end{eqnarray} (30) We then construct an objective function Ed as the L2 norm of the time-lapse residuals as \begin{equation} E_{d}= \left\Vert \delta {\Delta \mathbf {d}^{(1,0)}}\right\Vert ^2. \end{equation} (31) (Hereafter, the subscript ‘d’ refers to double-difference inversion.) Instead of minimizing Ed for both m(0) and m(1) or for both m(0) and Δm(1, 0), we minimize Ed solely for a time-lapse model change Δm(1, 0) by assuming that a given model mo is an accurate estimation of m(0). Replacing m(0) with mo in eq. (30) obtains \begin{eqnarray} {\Delta \mathbf {u}^{(1,0)}\left[\Delta \mathbf {m}^{(1,0)}, \mathbf {m}_o; \hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right]}\nonumber\\ {\quad = \mathbf {u}^{(1)}\left[\mathbf {m}_o+\Delta \mathbf {m}^{(1,0)}; \hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right] -\,\mathbf {u}^{(0)}\left[\mathbf {m}_o;\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right].} \end{eqnarray} (32) We implement the nonlinear double-difference method following Watanabe et al. (2004), Zhang & Huang (2013) and Asnaashari et al. (2014), and review the method in Appendix C. Note that our double-difference implementation remains fully nonlinear, since the assumption of mo being accurate does not lead to a linearization of the inverse problem, and also since we do not assume |$\Delta \mathbf {u}^{(1,0)} = \mathcal {F}\lbrace \Delta \mathbf {m}^{(1,0)}\rbrace$|⁠. The estimated model change |$\Delta \mathbf {m}_d^{(1,0)}[\Delta \mathbf {d}^{(1,0)}; \mathbf {m}_o,\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}]$| contains errors associated with time-lapse components |$\delta _d^{tl}$| and with a background model |$\delta _d^{bg}$| \begin{equation} \Delta \mathbf {m}^{(1,0)}_d= \Delta \mathbf {m}_{\text{true}}^{(1,0)}+\delta _d^{tl}+\delta _d^{bg}. \end{equation} (33) Eqs (31) and (32) show that the double-difference method is solving for Δm(1, 0) from Δd(1, 0), and hence the estimated model change |$\Delta \mathbf {m}_d^{(1,0)}$| can be potentially much more accurate than the value |$\Delta \mathbf {m}_p^{(1,0)}$| estimated by the parallel method. We obtain the best double-difference result, when we satisfy the two assumptions; the perfectly repeating source and receiver parameters, and the perfect knowledge of |$\mathbf {m}_{\text{true}}^{(0)}$| (i.e. |$\mathbf {m}_o=\mathbf {m}_{\text{true}}^{(0)}$|⁠). In such a case, the background error term |$\delta _d^{bg}$| is nil, and the time-lapse error term |$\delta _d^{tl}$| is minimal. Note that |$\delta _d^{tl}$| is still not zero, due to limitations in acquisition parameters, undesired time-lapse changes, and non-repeatable data noise (eq. 19). Degrading source and sensor repeatability contaminates the time-lapse observed waveforms Δd(1, 0) (eq. 18). As the predicted time-lapse waveforms are computed for incorrect source and sensor parameters |$\hat{\mathbf {x}}_s,\hat{\mathbf {x}}_r$| (eq. 29), the time-lapse residual δΔd can be significantly affected by the repeatability issues, resulting in increasing errors |$\delta _d^{tl}$| in the time-lapse model change estimates Δm(1, 0). Furthermore, we do not have perfect subsurface information. Thus the error in mo causes errors in propagation paths, and then can be mapped to |$\Delta \mathbf {m}^{(1,0)}_d$| as |$\delta _d^{bg}$|⁠. 2.5 Bootstrapping inversion Instead of directly minimizing the time-lapse data change Δd(1, 0) as in the double-difference method, we may sequentially invert for the two sets of the data in the order of recording time, and implicitly constrain time-lapse data changes without assuming neither perfect repeatability nor an accurate starting model. We refer to this as a bootstrapping approach, which has also been referred to as the sequential approach (Asnaashari et al. 2014). We start the inversion process by inverting the observed waveforms at T(0) to estimate the model m0 through minimizing the objective function \begin{equation} E^{(0)}_{b}=\left\Vert \delta \mathbf {d}{^{(0)}_b}\right\Vert ^2, \end{equation} (34) where δd(0) is the data residual at T(0) defined as \begin{equation} \delta \mathbf {d}^{(0)}_b = \mathbf {u}\left[\mathbf {m}^{(0)}; \mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right] - \mathbf {d}\left[\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]. \end{equation} (35) (Hereafter, the subscript ‘b’ refers to bootstrapping inversion.) The objective function and the data residual defined in eqs (34) and (35) are defined as in the same way as |$E^{(0)}_p$| and |$\delta \mathbf {d}^{(0)}_p$| in parallel inversion (eqs 20 and 22). By iteratively minimizing the objective function |$E^{(0)}_{b}$| with a local gradient method, we update the velocity model starting from the given initial baseline model mo, and obtain the new baseline model estimate at T(0) denoted by |$\mathbf {m}_b^{(0)}[\mathbf {d}^{(0)}; \mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}]$|⁠. Instead of starting the inversion from the initial model mo (as we have done for parallel inversion), we use the estimated model |$\mathbf {m}_b^{(0)}$| as a starting model to estimate m(1), and fit d(1) by minimizing the objective function \begin{equation} E^{(1)}_{b}=\left\Vert \delta \mathbf {d}{^{(1)}_b}\right\Vert ^2 , \end{equation} (36) where the residual |$\delta \mathbf {d}^{(1)}_b$| is defined as \begin{equation} \delta \mathbf {d}^{(1)}_b = \mathbf {u}\left[\mathbf {m}^{(1)}; \mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right] - \mathbf {d}\left[\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]. \end{equation} (37) Optimizing |$E^{(1)}_{b}$| with a local gradient method yields the estimation of the model at T(1) as |$\mathbf {m}_{b}^{(1)}[\mathbf {d}^{(1)}; \mathbf {m}^{(0)}_{b},\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}]$|⁠. The change in the model parameters is given by subtracting the estimated model at T(0) from the estimated model at T(1); \begin{eqnarray} \Delta \mathbf {m}_{b}^{(1,0)} = \mathbf {m}_{b}^{(1)}\left[\mathbf {d}^{(1)}; \mathbf {m}^{(0)}_{b},\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]-\mathbf {m}_{b}^{(0)}\left[\mathbf {d}^{(0)}; \mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right].\!\! \nonumber\\ \end{eqnarray} (38) Note that Δmb is equivalent to the total model update obtained in the second-stage inversion, as the second term |$\mathbf {m}^{(0)}_{b}$| on the right hand side of the equation is the starting model for the inversion. As in eq. (15), we can write \begin{equation} \Delta \mathbf {m}_{b}^{(1,0)} = \Delta \mathbf {m}_{\text{true}}^{(1,0)} + \delta _b^{tl} + \delta _b^{bg}, \end{equation} (39) where |$\delta _b^{tl}$| and |$\delta _b^{bg}$| are the differences of time-lapse and static background velocity errors between the two inversion results respectively. (See Appendix D for derivation of |$\delta _b^{tl}$| and |$\delta _b^{bg}$|⁠.) As in parallel inversion, the bootstrapping approach does not explicitly depend on the repeatability of the events or the recording parameters. However in contrast to the parallel approach, we start the inversion of the model at T(1) from the final inversion result at T(0), and thus the inversion can potentially concentrate largely on reducing the time-lapse data changes during the second stage of inversion. Therefore, we expect both of |$\delta _b^{tl}$| and |$\delta _b^{bg}$| to be smaller than those terms in parallel inversion (⁠|$\delta _p^{tl}$| and |$\delta _p^{bg}$|⁠). We can also expect a similar performance to the double-difference method, without enforcing the repeatability of the source and sensor parameters. Note that the bootstrapping inversion method is not still free from data repeatability issues, and is affected by the nonlinearity of both first- and second-stage inversions including the quality of the initial model. 3 SETTINGS OF INVERSION TESTS In order to examine the feasibility of time-lapse full waveform inversion and to compare different inversion approaches, we consider a hypothetical scenario of monitoring a producing reservoir using a dense array of sensors at the surface. We assume that production has caused changes in pore pressure, resulting in stress-induced velocity perturbations and repeating microseismic events occurring at times T(0) and T(1) near the reservoir during production. In order to reduce the complexity of the inverse problem, we assume that the medium is acoustic, and the source mechanism of all the events is isotropic. We focus on evaluating the potential for full waveform inversion to estimate time-lapse changes in P-wave velocity in the reservoir. We conduct acoustic full waveform inversion in the Laplace–Fourier domain following Kamei et al. (2013). Laplace–Fourier full waveform inversion (Shin & Cha 2009) can be cast as an extension of the frequency-domain full waveform inversion approach of Pratt (1999), and uses a Fourier transform of time-damped wavefields as data to explicitly control amount of later arrivals used for the inversion (Kamei et al. 2013). The time-damping is characterized by a decay function exp [−st], where s is the damping variable (also referred to as the Laplace variable). Adopting a large s leads to significant time-damping and an emphasis on inverting early arrival waveforms. See Appendix A for the brief review of the method. In this section, we first describe the P-wave velocity model, the generation of synthetic ‘observed’ data, and the basic parameters of inversion tests. 3.1 Model and data generation We use a realistic synthetic model based on an existing reservoir. The baseline P-wave velocity model at T(0) is displayed in Fig. 1(a). The model contains a wide range of complexity, including two large normal horst-graben faults and detailed small-scale heterogeneities. The reservoir is located at a depth of 4 km between bounding faults within thin layers that are approximately 50 m thickness. Sensors are densely located at the surface at 25 m intervals. Figure 1. Open in new tabDownload slide (a) True baseline velocity model (Model Gtrue) at T(0) and (b) true time-lapse velocity changes between T(0) and T(1). (c,d) Smoothed versions of the true velocity model at T(0) in (a) used as an initial estimate mo of a baseline velocity model in the inversion tests conducted in Section 6. Stars in (a) and (b) indicate the locations of repeating events SA, SB and SC. The Gaussian filter is applied to the velocity model with a standard deviation of (c) 100 m and (d) 300 m. We refer to the model in (c) as Model Gres1, and the model in (d) as Model Gres2. Figure 1. Open in new tabDownload slide (a) True baseline velocity model (Model Gtrue) at T(0) and (b) true time-lapse velocity changes between T(0) and T(1). (c,d) Smoothed versions of the true velocity model at T(0) in (a) used as an initial estimate mo of a baseline velocity model in the inversion tests conducted in Section 6. Stars in (a) and (b) indicate the locations of repeating events SA, SB and SC. The Gaussian filter is applied to the velocity model with a standard deviation of (c) 100 m and (d) 300 m. We refer to the model in (c) as Model Gres1, and the model in (d) as Model Gres2. We assume that the P-wave velocity has increased in the two thin layers (we refer to these layers as the ‘top’ and ‘bottom’ layers) due to the pore pressure changes in the reservoir during production. Between T(0) and T(1), values of velocity in the top layer increased by 15 per cent of the baseline velocity model, and velocities in the bottom layer increased by 7 per cent (see Fig. 1b). We then assume that three repeating isotropic-source seismic events occurred at a depth of 5 km; labelled SA along the fault at the distance of 5.5 km, SB at 7.75 km, and SC at 10 km (marked by the stars in Figs 1a–b). We also assume that the source parameters of the repeating events are similar and consistent over time. We denote the events at T(0) as SA(0), SB(0) and SC(0), and the events at T(1) as SA(1), SB(1), and SC(1). We compute synthetic data, d(0) and d(1), for the events using a Laplace–Fourier-domain finite-difference acoustic modelling code (Jo et al. 1996; Kamei et al. 2013) and applying absorbing boundary condition at all edges. We simulate the wavefields using a source wavelet with a dominant frequency at 5 Hz and frequency range between 0 and 15 Hz. The density model is generated from the P-wave velocity model by using Gardner's relationship (Gardner et al. 1974). We then convert the Laplace–Fourier-domain wavefields into the time domain by using the discrete Laplace–Fourier transform, and display the synthetic seismic records for events at SA in Figs 2(a) and (b). In Fig. 2(c), the time-lapse differences of the data Δd(1, 0) = d(1) − d(0) are shown. In Fig. 2(d), we display traces of the observed and time-lapse data shown in Figs 2(a)–(c) for sensors located at the horizontal distances of 4, 8 and 12 km to examine time and amplitude changes in detail. Due to the small time-lapse change in velocity, seismic data at T(0) and T(1) are very similar to each other, and the time-lapse data differences are small: the maximum amplitude of time-lapse waveforms (Δd(1, 0)) is less than 10 per cent of those of the observed records (d(0) and d(1)). We observe the clearest time-lapse changes in the first arrival waveforms recorded at sensors on the right-hand side of the model as small time shifts (the largest shift is 4 ms), and as a decay in amplitudes. Changes in later scattered arrivals (time-lapse coda) are also evident, though their magnitudes are smaller than those in the first arrivals. Figure 2. Open in new tabDownload slide Observed waveforms for (a) SA(0) (occurred at T(0)), (b) SA(1) (occurred at T(1)) and (c) the corresponding observed time-lapse waveforms recorded for the repeating events at SA, Δd = d(1) − d(0). The time-lapse waveforms are scaled by a factor of 10 for display purpose, and then (a–c) are displayed in the same colour scale. (d) Seismic traces recorded at sensors located at distances of 4, 8 and 12 km for events SA(0) and SA(1) in (a–c). The red line indicates the data observed at T(0) (same as (a)), the blue line indicates the data observed at T(1) (same as (b)), and the green line indicates time-lapse waveforms (same as (c), but without scaling). Figure 2. Open in new tabDownload slide Observed waveforms for (a) SA(0) (occurred at T(0)), (b) SA(1) (occurred at T(1)) and (c) the corresponding observed time-lapse waveforms recorded for the repeating events at SA, Δd = d(1) − d(0). The time-lapse waveforms are scaled by a factor of 10 for display purpose, and then (a–c) are displayed in the same colour scale. (d) Seismic traces recorded at sensors located at distances of 4, 8 and 12 km for events SA(0) and SA(1) in (a–c). The red line indicates the data observed at T(0) (same as (a)), the blue line indicates the data observed at T(1) (same as (b)), and the green line indicates time-lapse waveforms (same as (c), but without scaling). 3.2 Basic inversion parameters We perform a series of inversion tests in order to examine the capability of full waveform inversion to estimate the time-lapse velocity change from seismic waveforms excited by three repeating earthquakes (SA, SB, and SC). During the inversion, we use three input velocity models; the true (baseline) velocity model at T(0) (Fig. 1a), and two smoothed versions of the true baseline model (Figs 1c and d). Smoothing is performed in the slowness domain by applying a Gaussian filter of a standard deviation of 100 and 300 m. We refer to the true velocity model as Model Gtrue (Fig. 1a), the finer velocity model obtained by the 100-m Gaussian filter as Model Gres1 (Fig. 1c), and the smoother velocity model as Model Gres2 (Fig. 1d). The size of the narrower Gaussian filter corresponds to a half seismic wavelength at the reservoir depths at the dominant frequency of the source wavelet. The subsurface model can be obtained at such a spatial resolution by full waveform inversion of controlled-source seismic survey data (e.g. conducted before the production started), as the expected resolution of full waveform inversion is approximately a half wavelength (Wu & Toksoz 1987). In all inversion tests, we assume that we know the correct location and origin time of the events by performing a separate prior estimation problem. This is reasonable, as we can expect to accurately estimate these parameters from our relatively good initial velocity models. We minimize the L2 residuals of the Laplace–Fourier wavefields at frequencies between 2 and 15 Hz at intervals of 0.5 Hz and at a damping constant of 1 s−1. In order to compensate for excessive amplitude loss at far offsets by the damping function, we precondition the wavefields by applying an offset-weighting function exp [sD/c] where D is the offset, and c = 3000 m s−1 is the approximate average of the velocity models (Brenders & Pratt 2007). We display the equivalent preconditioned wavefields in the time-domain after the application of offset-weighting in Fig. 3 to examine the amount of information available for the inversion. Note that Fig. 3 is displayed in different scaling than Fig. 2. We confirm the suppression of arrivals later than 1 s from the first arrivals, and expect the inversion to fit the waveforms up to 1 s after the first arrivals where most of the significant time-lapse data changes are observed. During our noise-free inversion tests, instead of applying the Laplace–Fourier transform to the time-domain wavefields in Fig. 2, we invert for the Laplace–Fourier domain wavefields directly calculated in our Laplace–Fourier modelling program. This is to keep the data as noise-free as possible by avoiding conversion errors and instead focusing on comparisons of the three inversion methods. In contrast, for the non-repeatable data noise tests in Section 6.2, the noise-added data are generated using the time-domain data in Fig. 2, and transformed into the Laplace–Fourier domain to include the conversion errors. Figure 3. Open in new tabDownload slide Observed time-damped (a) waveforms for SA(0), and (b) time-lapse waveforms for SA. A time-damping function of exp [−(t − D/c)s] is applied to the raw waveforms in Fig. 2, where D is the offset, and c is a representative velocity value (3 km s−1 is used). The time-lapse waveforms are scaled by a factor of 10 for display purpose, and shown in the same colour scale as in (a). Note that the waveforms are shown in a different scale from that used in Fig. 2. Figure 3. Open in new tabDownload slide Observed time-damped (a) waveforms for SA(0), and (b) time-lapse waveforms for SA. A time-damping function of exp [−(t − D/c)s] is applied to the raw waveforms in Fig. 2, where D is the offset, and c is a representative velocity value (3 km s−1 is used). The time-lapse waveforms are scaled by a factor of 10 for display purpose, and shown in the same colour scale as in (a). Note that the waveforms are shown in a different scale from that used in Fig. 2. As the low-wavenumber spectral components of Models Gres1 and Gres2 are similar to those of the true model Gtrue, we do not observe severe cycle skipping in most of data at our frequency range (see discussions in each inversion tests). Thus we invert for all frequency components simultaneously, and do not employ a multiscale approach. The gradient of the objective function is filtered by an elliptic wavenumber filter to suppress undesired high-wavenumber oscillations (Sirgue & Pratt 2004). We select the elliptic wavenumber filter based on the theoretical resolution of full waveform inversion obtained by a Born approximation analysis |$k_o=\frac{2\omega }{c}$| (Wu & Toksoz 1987). We again set c = 3000 m s−1, and aspect ratio of the filter as |$\frac{1}{4}$| to enhance the horizontal continuity. The high-cut vertical wavenumber is |$\frac{3}{4} k_o$| tapered to zero at ko. We adopt a target-oriented approach, and mask the gradient to restrict the model update within a rectangular area of 12 km horizontally by 5 km vertically around the reservoir (Shragge & Lumley 2013). Combining the time-damping and target-oriented approach potentially accelerates the inversion convergence, by concentrating on both the model and data space most affected by the subsurface time-lapse velocity changes. After each iteration, we update the source wavelet using the linear inversion method described in Pratt (1999), and the density model using Gardner's relationship. 4 SENSITIVITY KERNEL Fréchet sensitivity kernels and wavepaths (Woodward 1992) indicate the illumination area of full waveform inversion, and provide insights into the inversion method's model resolution and behaviour. We can derive a monochromatic (source-free) sensitivity kernel with respect to P-wave slowness as \begin{equation} K(\mathbf {x}, \mathbf {x}_{r_j},\mathbf {x}_{s_i}) = \frac{(\omega +is)^2}{c(\mathbf {x})} G(\mathbf {x},\mathbf {x}_{s_i}) G(\mathbf {x}, \mathbf {x}_{r_j}), \end{equation} (40) where |$G(\mathbf {x},\mathbf {x}_{s_i})$| is Green's function at x excited by an earthquake source at |$\mathbf {x}_{s_i}$|⁠, c(x) is the velocity distribution, ω is the angular frequency, and s is the Laplace constant (see Appendix A for derivation). The gradient of the objective function with respect to P-wave slowness is then written using K as \begin{eqnarray} \nabla E = \Re \left\lbrace \sum _{i,j,l} K(\mathbf {x},\mathbf {x}_{r_j},\mathbf {x}_{s_i};\Omega _l) \delta d^*(\mathbf {x}_{r_j},\mathbf {x}_{s_i};\Omega _l) f(\mathbf {x}_{s_i};\Omega _l)\right\rbrace ,\nonumber\\ \end{eqnarray} (41) where ℜ is the real part of a complex function, * is the complex conjugate, and |$f(\mathbf {x}_{s_i})$| is the source function of the earthquake at |$\mathbf {x}_{s_i}$|⁠. In Fig. 4, we display monochromatic sensitivity kernels at 5 Hz for the sensors located at horizontal distances of 3, 7, 11 and 14 km and the event SA computed with the true velocity model at T(0). The central portion of the sensitivity kernels, depicted in dark grey colour, corresponds to the first Fresnel zone that is sensitive to first-arrival waveforms. The angle of the sensitivity kernels varies significantly over a range of sensors (or aperture from the event location). Sensitivity kernels are nearly vertical for sensors at a narrow aperture (Fig. 4b), and gradually tilt as the aperture between the event location and a sensor increases (Figs 4a and c). For sensors at the furthest offset, sensitivity kernels cover a large portion of the reservoir, as seismic waves travel in the horizontal direction within the high-velocity layers at the deep part of the model (Fig. 4d). These sensitivity kernels computed for different sensor locations overlap near the event SA in a variety of angles, which indicates that the reservoir can be reasonably well illuminated, and we expect that the velocity structures in the reservoir will be well constrained and imaged in high resolution. The sensitivity kernels become sparser, and the illumination becomes degraded outside the reservoir further from the event location. The observations suggest that having a few earthquakes may not be suitable to invert for the entire velocity model, but may be sufficient to recover local changes in the reservoir close to the event locations. Figure 4. Open in new tabDownload slide Real part of sensitivity kernels (wavepaths) at 5 Hz computed for the events at SA and a sensor at a distance of (a) 3, (b) 7, (c) 11 and (d) 14 km. Figure 4. Open in new tabDownload slide Real part of sensitivity kernels (wavepaths) at 5 Hz computed for the events at SA and a sensor at a distance of (a) 3, (b) 7, (c) 11 and (d) 14 km. 5 ESTIMATING TIME-LAPSE VELOCITY CHANGES BY USING THE TRUE BASELINE VELOCITY MODEL We first conduct full waveform inversion by assuming that we know the correct velocity model at T(0) to examine the maximum capability of full waveform inversion for estimating time-lapse velocity changes in the reservoir from a limited number of repeating earthquakes. By starting from the exact baseline velocity model (i.e. δbg = 0), we do not need to use the three time-lapse inversion methods described in Section 2, because we can directly invert the data set d(1) at T(1). We use the inversion parameters described in Section 3.2 to estimate |$\Delta \mathbf {m}=\mathbf {m}^{(1)} - \mathbf {m}_{\text{true}}^{(0)}$|⁠, the changes in the velocity model. The inversion converges well: after 200 iterations, the objective function is reduced to 0.07 per cent of that of the starting model. We display the velocity change model in Fig. 5(a). Fig. 5(b) shows the true velocity change smoothed by the Gaussian filter with a standard deviation of 50 m in order to compare the results at the level of resolution of full waveform inversion. We observe that full waveform inversion recovers the velocity change very well. The estimated model clearly delineates the top and bottom layers, though inversion artefacts exist and thus δΔm(1, 0) = δtl ≠ 0. In order to evaluate waveform fitting, we synthesize the seismic record u(1) at T(1) from the final velocity model, and display them for SA(1) in Fig. 6(a), and the corresponding residual waveforms in Fig. 6(b). Note that we display waveforms and residuals without damping to examine waveform fitting of the entire record including later arrivals. The inverted model reproduces the early part of the waveforms well, fitting both phase and amplitudes up to 1 s from the first arrivals, as expected from the selected time-damping value of s. The experiments suggest that full waveform inversion can accurately estimate detailed time-lapse velocity changes from the seismic data of a sparse set of repeating earthquakes. Figure 5. Open in new tabDownload slide (a) Inverted velocity changes when we assume that we know the correct velocity model at T(0). (b) Smoothed true velocity change model obtained by applying a Gaussian spatial filter to the original true velocity change model displayed in Fig. 1(b), so that its spatial resolution is similar to the inversion result in (a). Note that the colour scale is different from that used in Fig. 1(b). Figure 5. Open in new tabDownload slide (a) Inverted velocity changes when we assume that we know the correct velocity model at T(0). (b) Smoothed true velocity change model obtained by applying a Gaussian spatial filter to the original true velocity change model displayed in Fig. 1(b), so that its spatial resolution is similar to the inversion result in (a). Note that the colour scale is different from that used in Fig. 1(b). Figure 6. Open in new tabDownload slide (a) Predicted and (b) residual waveforms for the event SA(1) occurred at T(1) computed from the inverted model in Fig. 5(a). The inversion is performed by assuming that we know the correct velocity model at T(0). The corresponding observed waveforms are shown in Fig. 2(b). Note that no time-damping is applied to the waveforms. Waveforms are displayed with the same colour scale as in Fig. 2. The residual waveforms in (b) are scaled by a factor of 10 for display purposes, but too small to recognize. Figure 6. Open in new tabDownload slide (a) Predicted and (b) residual waveforms for the event SA(1) occurred at T(1) computed from the inverted model in Fig. 5(a). The inversion is performed by assuming that we know the correct velocity model at T(0). The corresponding observed waveforms are shown in Fig. 2(b). Note that no time-damping is applied to the waveforms. Waveforms are displayed with the same colour scale as in Fig. 2. The residual waveforms in (b) are scaled by a factor of 10 for display purposes, but too small to recognize. 6 INVERSION OF REPEATING EARTHQUAKE DATA USING A SMOOTH VELOCITY MODEL An inexact initial velocity model potentially affects the performance of time-lapse full waveform inversion due to δbg as discussed in Section 2. In this section, we evaluate the potential of parallel, double-difference, and bootstrapping inversion methods to recover the time-lapse velocity change, by starting inversions from the imperfect models Gres1 and Gres2. We start our investigation with noise-free data (Section 6.1), followed by testing non-repeatable noise effects (Section 6.2), and non-repeatable source parameters in Section 6.3. We finish our synthetic tests by comparing results from different initial models in Section 6.4. 6.1 Inversion of noise-free data First we estimate velocity changes from the noise-free data sets (i.e. d(0) = s(0), d(1) = s(1)) by using Model Gres1 in Fig. 1(c), as an initial model, and also by assuming perfectly repeating source and receiver parameters (i.e. |$\mathbf {x}_s^{(0)}=\mathbf {x}_s^{(1)}$| and |$\mathbf {x}_r^{(0)}=\mathbf {x}_r^{(1)}$|⁠). We analyse the effects of an imperfect initial model, and compare the three methods (parallel, double-difference and bootstrapping). 6.1.1 Effects of velocity error on waveforms To understand the effects of an imperfect velocity model, we first evaluate discrepancies in waveforms between the observed and synthetic records at T(0). We compute synthetic waveforms u(0) from Model Gres1, and display the synthetic record for the event SA(0) in Fig. 7(a). We display the differences between the observed and synthetic waveforms, (i.e. residual waveforms, δd(0)[mo] = u(0)[mo] − d(0)), in Fig. 7(b). To provide a detailed examination of waveforms near the first arrivals, Fig. 7(c) presents overlays of the observed, synthetic and residual waveforms recorded at the sensors at horizontal distances of 4, 8, and 12 km. As long-wavelength components of the smoothed baseline model largely agree with those of the true model at T(0), the phases of the synthetic waveforms coincide well with those of the observed waveforms near the first arrivals. However, the amplitudes show slight discrepancies. The synthetic waveforms do not contain later scattered arrivals due to the lack of high-wavenumber model components, but we can expect using a damping variable s reduces the effect of missing later arrivals on time-lapse full waveform inversion. When we compare the residual waveforms δd(0)[mo] in Fig. 7(b) with the observed time-lapse waveforms Δd(1, 0) in Fig. 2(c) (displayed in the same scale), the former are much larger in amplitudes than the latter, suggesting potential contamination of velocity change estimation by the background velocity errors (δbg ≠ 0). Figure 7. Open in new tabDownload slide (a) Predicted and (b) residual waveforms for the event SA(0) at T(0) computed from the smoothed velocity model Gres1. Waveforms are displayed with the same colour scale as in Fig. 2, while the residual waveforms in (b) are scaled by a factor of 10 for display purposes. (c) Seismic traces of (a and b) recorded at sensors located at distances of 4, 8 and 12 km. The red line indicates the observed data at T(0) (same as Fig. 2a), the blue line indicates the synthetic waveforms from Gres1 (same as (a)), and the green line indicates residual waveforms (same as (b) but without scaling). The magnitudes of the waveform residuals at T(0) are larger than the time-lapse waveforms in Fig. 2(c). No time-damping is applied to the waveforms. Figure 7. Open in new tabDownload slide (a) Predicted and (b) residual waveforms for the event SA(0) at T(0) computed from the smoothed velocity model Gres1. Waveforms are displayed with the same colour scale as in Fig. 2, while the residual waveforms in (b) are scaled by a factor of 10 for display purposes. (c) Seismic traces of (a and b) recorded at sensors located at distances of 4, 8 and 12 km. The red line indicates the observed data at T(0) (same as Fig. 2a), the blue line indicates the synthetic waveforms from Gres1 (same as (a)), and the green line indicates residual waveforms (same as (b) but without scaling). The magnitudes of the waveform residuals at T(0) are larger than the time-lapse waveforms in Fig. 2(c). No time-damping is applied to the waveforms. 6.1.2 Update of baseline velocity model Because the data residual δd[mo] arising from the discrepancies between Model Gres1 and the true velocity model at T(0) is greater than the time-lapse signal δs = Δd(1, 0), we first aim to reduce the data residuals by updating the model through the inversion of the waveforms at T(0). We use the inversion parameters described in Section 3.2, except that we do not mask the gradient, as model errors span the entire model. We iteratively update the velocity model for 200 times, and display the resulting velocity model m(0) in Fig. 8(a). The objective function is monotonically reduced to 2.5 per cent of the value computed from the initial model (Model Gres1). We then synthesize seismic waveforms u(0)[m(0)] from the inverted model, and display the predicted and residual waveforms in Figs 8(b)–(d). The predicted waveforms agree well with the observed waveforms, and the residual waveforms δd(0)[m(0)] are small when compared to those from the initial baseline model (Fig. 7b). However, the resulting velocity model is strongly contaminated by artefacts (δm(0)), which contrasts sharply with the previous stable inversion test in Section 5 that only needed to recover the local velocity changes near the events. Now that the background model errors span the entire model, the sparse source distribution does not provide sufficient illumination, and hinders the recovery of a clean velocity model as expected from the sensitivity analysis in Section 4. In the following section, we refer to the updated baseline velocity model as Model G|$^{\prime }_{\mathrm{res1}}$|⁠. Figure 8. Open in new tabDownload slide Baseline velocity model inversion: (a) Inverted baseline velocity model (Model G|$^{\prime }_\mathrm{res1}$|⁠) by using Model Gres1 as an initial model. (b) Predicted and (c) residual waveforms for SA(0) computed from the inverted velocity model in (a). Waveforms in (b and c) are displayed with the same colour scale as in Fig. 2, while the residual waveforms in (c) are scaled by a factor of 10 for display purposes. (d) Seismic traces of SA(0) recorded at sensors located at distances of 4, 8 and 12 km. The red lines indicate the observed data at T(0) (same as Fig. 2a), the blue lines indicate the predicted data (same as (a)), and the green lines indicate residuals (same as (b)). The blue lines are overlain by the red lines and thus not visible, since the velocity model reproduces the first 1 s very accurately as found in (b). No time-damping is applied to the waveforms. Figure 8. Open in new tabDownload slide Baseline velocity model inversion: (a) Inverted baseline velocity model (Model G|$^{\prime }_\mathrm{res1}$|⁠) by using Model Gres1 as an initial model. (b) Predicted and (c) residual waveforms for SA(0) computed from the inverted velocity model in (a). Waveforms in (b and c) are displayed with the same colour scale as in Fig. 2, while the residual waveforms in (c) are scaled by a factor of 10 for display purposes. (d) Seismic traces of SA(0) recorded at sensors located at distances of 4, 8 and 12 km. The red lines indicate the observed data at T(0) (same as Fig. 2a), the blue lines indicate the predicted data (same as (a)), and the green lines indicate residuals (same as (b)). The blue lines are overlain by the red lines and thus not visible, since the velocity model reproduces the first 1 s very accurately as found in (b). No time-damping is applied to the waveforms. 6.1.3 Parallel inversion As the first of the three time-lapse inversion methods proposed in Section 2, we evaluate the parallel inversion method described in Section 2.3. We follow the inversion procedure described in Section 3.2 during 200 iterations of individual full-waveform-inversion runs. The objective functions (⁠|$E_{p}^{(0)}$| and |$E_p^{(1)}$|⁠) monotonically decrease during both runs; we display the history of the objective functions in Fig. 9(a) where the solid and dashed lines indicate |$E_{p}^{(0)}$| and |$E_p^{(1)}$| respectively. The two curves follow a similar trend, and are reduced to approximately 11 per cent of the initial value after 200 iterations. The final values of the objective functions are larger than that in Section 6.1.2, because the previous inversion is conducted without applying gradient masking. In order to evaluate the convergence of the time-lapse component, we compute the residuals of the time-lapse waveform difference δΔd(1, 0) in eq. (28) after each iteration, and display the history of the L2 norm of the time-lapse residuals ‖δΔd(1, 0)‖2 in Fig. 9(b). Note that we apply the damping function exp (−st) with the same value of s as used in the inversion to validate the performance of the inversion. The time-domain residuals generally decrease, and ‖δΔd(1, 0)‖2 becomes less than 4 per cent of the pre-inversion value. However, the convergence is not monotonic, since the time-lapse residuals are not used to constrain the inversion. Figure 9. Open in new tabDownload slide Parallel inversion: History of (a) the L2 norm of waveform residuals (⁠|$E_p^{(0)}=\Vert \delta \mathbf {d}^{(0)}\Vert ^2$| and |$E_p^{(1)}=\Vert \delta \mathbf {d}^{(1)}\Vert ^2$|⁠) and (b) the L2 norm of time-lapse residuals (‖δΔd(1, 0)‖2) during parallel inversion using Model Gres1 as an initial model. The damping parameter s = 1 s−1 (the same value used for the inversions) is applied to compute residuals. |$E_p^{(0)}$| and |$E_p^{(1)}$| are used as the objective function of the inversions. The left axes indicate the absolute value, and the right axes indicate the value relative to that before the inversion. In (a), the solid line indicates |$E_p^{(0)}$|⁠, the objective function during the inversion of the seismic data d(0) at T(0), and the dashed line indicates |$E_p^{(1)}$|⁠, the objective function during the inversion of the data d(1) at T(1). Figure 9. Open in new tabDownload slide Parallel inversion: History of (a) the L2 norm of waveform residuals (⁠|$E_p^{(0)}=\Vert \delta \mathbf {d}^{(0)}\Vert ^2$| and |$E_p^{(1)}=\Vert \delta \mathbf {d}^{(1)}\Vert ^2$|⁠) and (b) the L2 norm of time-lapse residuals (‖δΔd(1, 0)‖2) during parallel inversion using Model Gres1 as an initial model. The damping parameter s = 1 s−1 (the same value used for the inversions) is applied to compute residuals. |$E_p^{(0)}$| and |$E_p^{(1)}$| are used as the objective function of the inversions. The left axes indicate the absolute value, and the right axes indicate the value relative to that before the inversion. In (a), the solid line indicates |$E_p^{(0)}$|⁠, the objective function during the inversion of the seismic data d(0) at T(0), and the dashed line indicates |$E_p^{(1)}$|⁠, the objective function during the inversion of the data d(1) at T(1). The estimated (absolute) velocity models |$\mathbf {m}_p^{(0)}$| and |$\mathbf {m}_p^{(1)}$| at T(0) and T(1) (not shown) resemble the estimated Model G|$^{\prime }_{\mathrm{res1}}$| displayed in Fig. 8(a), and do not appear geologically plausible. The artefacts are caused by errors in the time-invariant component of the initial model, and manifested by the limited and sparse event distribution. Following eq. (24), we estimate the time-lapse velocity change |$\Delta \mathbf {m}_p^{(1,0)}=\mathbf {m}_p^{(1)}-\mathbf {m}_p^{(0)}$|⁠, and display the model change after 50 and 200 iterations in Figs 10(a)–(b). In spite of the poor recovery of the absolute velocity models, the parallel inversion results indicate a velocity increase in the top layer, suggesting most of the artefacts in |$\mathbf {m}_p^{(0)}$| and |$\mathbf {m}_p^{(1)}$| associated with the background model errors δmbg are cancelled out. When compared to the tests with the true baseline velocity model in Section 5, the estimated velocity change contains more artefacts, and the inversion underestimates the velocity increase. With an increased number of iterations, the parallel inversion approach gradually delineates the top layer, but at the same time induces stronger artefacts especially at depths below 3 km (i.e. δmbg no longer cancels out and δbg becomes pronounced.). Differentiating between artefacts and the actual velocity features become difficult after 200 iterations (Fig. 10b), confirming that the two non-unique inversion runs do not necessarily converge to an optimal time-lapse velocity change estimation by performing more iterations. Figure 10. Open in new tabDownload slide Inverted velocity changes obtained by the (a,b) parallel, (c–f) double-difference and (g,h) bootstrapping inversion methods after (a,c,e,g) 50 and (b,d,f,h) 200 iterations. Model Gres1 is used as an initial baseline velocity model. Double-difference inversion is conducted with (c,d) the original Model Gres1 and (e,f) the updated baseline Model G|$^{\prime }_\mathrm{res1}$| in Fig. 8(a). Figure 10. Open in new tabDownload slide Inverted velocity changes obtained by the (a,b) parallel, (c–f) double-difference and (g,h) bootstrapping inversion methods after (a,c,e,g) 50 and (b,d,f,h) 200 iterations. Model Gres1 is used as an initial baseline velocity model. Double-difference inversion is conducted with (c,d) the original Model Gres1 and (e,f) the updated baseline Model G|$^{\prime }_\mathrm{res1}$| in Fig. 8(a). We perform the further validation of the parallel inversion method based on the velocity models after 50 iterations. First we synthesize the waveforms |$\mathbf {u}_p^{(0)}[\mathbf {m}_p^{(0)}]$| and |$\mathbf {u}_p^{(1)}[\mathbf {m}_p^{(1)}]$| from the estimated velocity models. Note that we do not apply the time damping function to evaluate the waveforms for the entire recording time. We display the predicted waveforms for SA(0) in Fig. 11(a), and the residual waveforms |$\delta \mathbf {d}_p^{(0)}[\mathbf {m}_p^{(0)}]=\mathbf {u}^{(0)}[\mathbf {m}_p^{(0)}] - \mathbf {d}^{(0)}$| in Fig. 11(b). We observe that the residuals are reduced by the inversion, but are still prominent even near the first arrivals. We consider that the application of the gradient mask results in insufficient reduction in the data residuals, since the inversion cannot reduce the background velocity errors outside the targeted area. Figure 11. Open in new tabDownload slide Parallel inversion: (a) Predicted and (b) residual waveforms for the event SA(0) computed from the velocity model obtained by parallel inversion after 50 iterations. Corresponding time-lapse (c) waveforms and (d) residuals. Parallel inversion is conducted by using Model Gres1 as an initial velocity model. The estimated velocity change is displayed in Fig. 10(d). Waveforms in (a–d) are displayed with the same colour scale as in Fig. 2, while the residual and time-lapse waveforms in (b,c) are scaled by a factor of 10, and the time-lapse residuals are by 20 for display purposes. Figure 11. Open in new tabDownload slide Parallel inversion: (a) Predicted and (b) residual waveforms for the event SA(0) computed from the velocity model obtained by parallel inversion after 50 iterations. Corresponding time-lapse (c) waveforms and (d) residuals. Parallel inversion is conducted by using Model Gres1 as an initial velocity model. The estimated velocity change is displayed in Fig. 10(d). Waveforms in (a–d) are displayed with the same colour scale as in Fig. 2, while the residual and time-lapse waveforms in (b,c) are scaled by a factor of 10, and the time-lapse residuals are by 20 for display purposes. We then compute the predicted time-lapse waveform change |$\Delta \mathbf {u}_p^{(1,0)}=\mathbf {u}_p^{(1)}-\mathbf {u}_p^{(0)}$|⁠. We show |$\Delta \mathbf {u}_p^{(1,0)}$| after 50 iterations for the repeating events at SA in Fig. 11(c), and the corresponding residual |$\delta \Delta \mathbf {d}_p^{(1,0)}=\Delta \mathbf {u}_p^{(1,0)} -\Delta \mathbf {d}^{(1,0)}$| in Fig. 11(d). The predicted time-lapse waveform change agrees well with the observed waveform change near the first arrivals, validating the estimated time-lapse velocity changes. Note the subsequent later scattered arrivals are not well reproduced due to the time-damping applied during the inversion. 6.1.4 Double-difference inversion Next we examine the double-difference method that focuses on minimizing the time-lapse residuals. As described in Section 2.4, using an appropriate baseline velocity model is a key to properly eliminate the propagation path effect outside the reservoir. We can use either the original velocity model Gres1 in Fig. 1, or the updated velocity model G|$^{\prime }_{\mathrm{res1}}$| in Fig. 8 that fits the record better but is not geologically plausible. We perform two sets of double-difference inversions using the models Gres1 and G|$^{\prime }_{\mathrm{res1}}$| as a baseline velocity model, respectively. We display the history of the objective functions Ed during 200 iterations in Fig. 12, and the estimated velocity changes Δmd after 50 and 200 iterations in Figs 10(c)–(f). Figure 12. Open in new tabDownload slide Double difference Inversion: History of the objective functions (Ed = ‖δΔd(1, 0)‖2) during the double-difference method. The damping parameter s = 1 s−1 (the same value used for the inversions) is applied. The left axes indicate the absolute value, and the right axes indicate the value relative to that before the inversion. The dashed line indicates the history during the inversion using the original initial velocity model (Model Gres1) and the solid line indicates the history during the inversion using the inverted velocity model (Model G|$^{\prime }_\mathrm{res1}$|⁠) as an input to double-difference inversion. Figure 12. Open in new tabDownload slide Double difference Inversion: History of the objective functions (Ed = ‖δΔd(1, 0)‖2) during the double-difference method. The damping parameter s = 1 s−1 (the same value used for the inversions) is applied. The left axes indicate the absolute value, and the right axes indicate the value relative to that before the inversion. The dashed line indicates the history during the inversion using the original initial velocity model (Model Gres1) and the solid line indicates the history during the inversion using the inverted velocity model (Model G|$^{\prime }_\mathrm{res1}$|⁠) as an input to double-difference inversion. The objective functions Ed in eq. (31) decrease monotonically, in contrast to the non-monotonic behaviour of the time-lapse misfit observed in the parallel inversion method (Fig. 9). The inversion using the estimated model G|$^{\prime }_{\mathrm{res1}}$| converges more rapidly than the inversion using the original velocity model Gres1. The objective functions after 200 iterations are 1.8 and 1.1 per cents of those computed before the inversion when using the models Gres1 and G|$^{\prime }_{\mathrm{res1}}$| respectively. Both double-difference inversions recover the top layer well. The amount of artefacts increases as the inversions progress, but are more suppressed when compared to the parallel inversion results, suggesting the effectiveness of explicitly constraining the time-lapse data changes. Double-difference inversion using Model Gres1 induces stronger and more coherent artefacts (δbg) than the inversion with Model G|$^{\prime }_{\mathrm{res1}}$|⁠, making identifying the geological layers with velocity change more difficult. Model G|$^{\prime }_{\mathrm{res1}}$| predicts wave propagation paths more accurately than Model Gres1, as indicated by the waveform fit. The observation supports the importance of the input model to the double-difference model as suggested in Section 2.4. We synthesize waveforms from the inversion results, and show the predicted and residual time-lapse waveforms (⁠|$\Delta \mathbf {d}_d^{(1,0)}$| and |$\delta \Delta \mathbf {d}_d^{(1,0)}$|⁠) after 50 iterations in Fig. 13. We observe clear improvements in the waveform fit, when compared to the parallel inversion result. The inverted time-lapse model changes predict the phases and amplitudes of the time-lapse waveforms very accurately near the first arrivals. The inversion using the estimated model G|$^{\prime }_{\mathrm{res1}}$| reproduces later arrivals better, especially after 0.3 s of the first arrivals. The observation further confirms the superior performance of the inversion using the updated baseline model G|$^{\prime }_{\mathrm{res1}}$| to that using the original model Gres1. Figure 13. Open in new tabDownload slide Double-difference inversion: (a,c) Predicted and (b,d) residual time-lapse waveforms computed from the inverted velocity changes in Figs 10(c) and (e) obtained by the double-difference method after 50 iterations. The velocity models are obtained by using (a,b) the original initial velocity model (Model Gres1) and (c,d) the inverted velocity model (Model G|$^{\prime }_\mathrm{res1}$|⁠) in Fig. 8 as the velocity model at T(0). No time-damping is applied to the waveforms. Waveforms are displayed with the same colour scale as in Fig. 2, while the residual and time-lapse waveforms in (b,c) are scaled by a factor of 10, and the time-lapse residuals are by 20 for display purposes. Figure 13. Open in new tabDownload slide Double-difference inversion: (a,c) Predicted and (b,d) residual time-lapse waveforms computed from the inverted velocity changes in Figs 10(c) and (e) obtained by the double-difference method after 50 iterations. The velocity models are obtained by using (a,b) the original initial velocity model (Model Gres1) and (c,d) the inverted velocity model (Model G|$^{\prime }_\mathrm{res1}$|⁠) in Fig. 8 as the velocity model at T(0). No time-damping is applied to the waveforms. Waveforms are displayed with the same colour scale as in Fig. 2, while the residual and time-lapse waveforms in (b,c) are scaled by a factor of 10, and the time-lapse residuals are by 20 for display purposes. 6.1.5 Bootstrapping inversion As the last of the three time-lapse inversion methods, we examine bootstrapping inversion that aims to implicitly minimize the residuals by sequentially inverting for the velocity models (see Section 2.5). We use the initial baseline velocity inversion in Section 6.1.2 at the first stage of the method, and treat the obtained Model G|$^{\prime }_{\mathrm{res1}}$| as |$\mathbf {m}_b^{(0)}$|⁠. As the data residual |$\delta \mathbf {d}_b[\mathbf {m}_b^{(0)}]$| shown in Fig. 8(c) is much smaller than the time-lapse data Δd(1, 0), we expect the second stage of the inversion to focus on estimating time-lapse component Δm(1, 0) rather than updating the remaining background model error. We estimate the velocity model at T(1) following the inversion parameters in Section 2, and Fig. 14(a) presents the history of |$E_b^{(1)}$| during 200 iterations of the second-stage inversion. As in the parallel inversion tests, we calculate the history of the L2 norm of the time-lapse residuals (⁠|$\Vert \delta \Delta \mathbf {d}_b^{(1,0)}\Vert$|⁠), and display it in Fig. 14(b). By subtracting Model G|$^{\prime }_{\mathrm{res1}}$| from the inverted velocity model |$\mathbf {m}_b^{(1)}$| at T(1), we calculate the time-lapse velocity change |$\Delta \mathbf {m}_b^{(1,0)}$|⁠, and display the estimated changes after 50 and 200 iterations in Figs 10(g)–(h). Figure 14. Open in new tabDownload slide Bootstrapping inversion: History of (a) the objective function (⁠|$E^{(1)}_b=\Vert \delta \mathbf {d}^{(1)}\Vert ^2$|⁠) and (b) L2 norm of time-lapse residuals (‖δΔd(1, 0)‖2) during the second stage of the bootstrapping method. Model Gres1 is used as the initial velocity model for the bootstrapping method, and then Model G|$^{\prime }_{\mathrm{res1}}$| is used as the input for the second-stage inversion. The damping parameter s = 1 s−1 (the same value used for the inversions) is used. The left axes indicate the absolute value, and the right axes indicate the value relative to that before the inversion. Relative L2 norms on the right axes appear larger than those in Fig. 9, as computed using the starting value of the second-stage inversion as a reference, while the left axis is set to be consistent with those in Fig. 9. Therefore the relative L2 norm in (a) range between 0 and 4. Figure 14. Open in new tabDownload slide Bootstrapping inversion: History of (a) the objective function (⁠|$E^{(1)}_b=\Vert \delta \mathbf {d}^{(1)}\Vert ^2$|⁠) and (b) L2 norm of time-lapse residuals (‖δΔd(1, 0)‖2) during the second stage of the bootstrapping method. Model Gres1 is used as the initial velocity model for the bootstrapping method, and then Model G|$^{\prime }_{\mathrm{res1}}$| is used as the input for the second-stage inversion. The damping parameter s = 1 s−1 (the same value used for the inversions) is used. The left axes indicate the absolute value, and the right axes indicate the value relative to that before the inversion. Relative L2 norms on the right axes appear larger than those in Fig. 9, as computed using the starting value of the second-stage inversion as a reference, while the left axis is set to be consistent with those in Fig. 9. Therefore the relative L2 norm in (a) range between 0 and 4. The objective function |$E^{(1)}_b$| monotonically decreases to approximately 10 per cent of the value computed at the start of the second-stage inversion, which is the value nearly equivalent to the parallel inversion result. However the history of the time-lapse data misfit is not oscillatory as observed in Fig. 9(b). The time-lapse misfit nearly monotonically decreases during the first 50 iterations, reaching a similar value to double-difference inversion, before starting to diverge. The estimated time-lapse velocity change well delineates the velocity increase in the top layer after 50 iterations (Fig. 10g). The estimated velocity change after 50 iterations is similar to the models estimated by the double-difference method in Figs 10(c) and (e), and hence much improved from the parallel inversion result shown in Fig. 10(a). The quality of the estimated model gradually degrades with the increasing number of iterations, as expected from the divergence of the time-lapse residuals. The bootstrapping model after 200 iterations (Fig. 10h) is more heavily contaminated by artefacts than the double-difference models (Figs 10d and f). We synthesize waveforms, and display the predicted time-lapse waveforms and residuals after 50 iterations in Fig. 15. We observe good fit in early arrival waveforms, validating our inversion results. The magnitude of the time-lapse data residuals |$\delta \Delta \mathbf {d}_b^{(1,0)}$| is similar to that of the double-difference results, but smaller than that of the parallel inversion result, as expected from the quality of the estimated model changes. We also note that after 50 iterations, the magnitude of the time-lapse data residual |$\delta \Delta \mathbf {d}_b^{(1,0)}$| is comparable to that of the data residual |$\delta \mathbf {d}[\mathbf {m}_b^{(0)}]$| associated with the remaining (background) model errors in Model G|$^{\prime }_{\mathrm{res1}}$| and shown in Fig. 8(c). Thus, the first 50 iterations are driven to implicitly minimize |$\delta \Delta \mathbf {d}_b^{(1,0)}$| that dominate the residual δd(1), the inversion is well behaved, and δbg < δtl. However, during the next 150 iterations the inversion tries to reduce both |$\delta \Delta \mathbf {d}_b^{(1,0)}$| and |$\delta \mathbf {d}[\mathbf {m}_b^{(0)}]$|⁠, and δbg increases. Then, the later stage of bootstrapping inversion behaves similar to parallel inversion, and increases the level of artefacts and the time-lapse misfit |$\Vert \delta \Delta \mathbf {d}_b^{(1,0)}\Vert ^2$| where δbg > δtl. Figure 15. Open in new tabDownload slide Bootstrapping inversion: Time-lapse (a) predicted and (b) residual waveforms computed from the inverted velocity changes after 50 iterations in Fig. 10(g) obtained by the bootstrapping method. Model Gres1 is used as the initial velocity model. Note that no time-damping is applied to the waveforms. Waveforms are displayed with the same colour scale as in Fig. 2, while the time-lapse waveforms in (a) are scaled by a factor of 10, and the time-lapse residuals in (b) are by 20 for display purposes. Figure 15. Open in new tabDownload slide Bootstrapping inversion: Time-lapse (a) predicted and (b) residual waveforms computed from the inverted velocity changes after 50 iterations in Fig. 10(g) obtained by the bootstrapping method. Model Gres1 is used as the initial velocity model. Note that no time-damping is applied to the waveforms. Waveforms are displayed with the same colour scale as in Fig. 2, while the time-lapse waveforms in (a) are scaled by a factor of 10, and the time-lapse residuals in (b) are by 20 for display purposes. 6.2 Effects of non-repeatable noise So far we have conducted the inversion tests using noise-free data. In order to evaluate the sensitivity of full waveform inversion to non-repeatable noise, we create Gaussian noise with zero mean and standard deviation corresponding to 5 and 10 per cent of the maximum amplitude of the noise-free records, and add a different realization of the noise (n(0) and n(1)) to the time-domain wavefields (s(0) and s(1)) at T(0) and T(1) displayed in Fig. 2 to create noise-added records (d(0) = s(0) + n(0) and d(1) = s(1) + n(1)). We display the 5-per cent noise-added records in Fig. 16(a), and the 10-per cent noise-record in Fig. 16(d). We also show time-lapse waveform change Δd(1, 0) of the both noise-added seismic records in Figs 16(b) and (e). In Figs 16(c) and (f), we display traces of the noise-added data , and their differences (i.e. the non-repeatable noise). Signal-to-noise ratios of the time-lapse wavefields are much lower than the recorded waveforms: The standard deviations of the time-lapse noise Δn(1, 0) of the 5- and 10-per cent noise-added seismic records are approximately 50 and 100 per cent of the maximum amplitude of the noise-free time-lapse waveform change Δs(1, 0). This suggests a strong impact of non-repeatable noise on time-lapse full waveform inversion, as indicated in Section 2. Figure 16. Open in new tabDownload slide Effects of non-repeatable noise: Noise-added (a,d) observed and (b,e) time-lapse waveforms for the event SA(0). (a,b) 5 and (c,d) 10 per cent uncorrelated Gaussian noise are added to noise-free seismic records, and noise at T(0) and T(1) are generated using different realizations. Waveforms in (a,b,d,e) are displayed with the same colour scale as in Fig. 2, while the time-lapse waveforms in (b,e) are scaled by a factor of 10 for display purposes. (c,f) Seismic traces of (c) 5 and (d) 10 per cent noise-added data recorded at sensors located at distances of 4, 8 and 12 km for the repeating events at SA. The red lines indicate the noise-added observed waveforms at T(0) (same as (c) in (a) or (f) in (d)), the blue lines indicate the noise-added observed waveforms at T(1), and the green lines indicate the noise-added observed time-lapse waveforms (same as (c) in (b) or (f) in (d), but without scaling). Note that signal-to-noise ratio of time-lapse waveforms in (b,e) are much worse than the recorded waveforms, and the time-lapse signals are difficult to distinguish. No time-damping is applied to the waveforms. Figure 16. Open in new tabDownload slide Effects of non-repeatable noise: Noise-added (a,d) observed and (b,e) time-lapse waveforms for the event SA(0). (a,b) 5 and (c,d) 10 per cent uncorrelated Gaussian noise are added to noise-free seismic records, and noise at T(0) and T(1) are generated using different realizations. Waveforms in (a,b,d,e) are displayed with the same colour scale as in Fig. 2, while the time-lapse waveforms in (b,e) are scaled by a factor of 10 for display purposes. (c,f) Seismic traces of (c) 5 and (d) 10 per cent noise-added data recorded at sensors located at distances of 4, 8 and 12 km for the repeating events at SA. The red lines indicate the noise-added observed waveforms at T(0) (same as (c) in (a) or (f) in (d)), the blue lines indicate the noise-added observed waveforms at T(1), and the green lines indicate the noise-added observed time-lapse waveforms (same as (c) in (b) or (f) in (d), but without scaling). Note that signal-to-noise ratio of time-lapse waveforms in (b,e) are much worse than the recorded waveforms, and the time-lapse signals are difficult to distinguish. No time-damping is applied to the waveforms. We conduct the parallel, double-difference and bootstrapping inversions by using the same strategies employed in Section 6.1. The noise-added time-domain data sets require pre-processing (the other data sets are created in the Laplace–Fourier domain by the forward modelling routine used for the inversion): we mute the noise before the first arrivals, in order to avoid excessive amplification of the noise before the first arrivals by the time-damping function exp [−st] during the inversion. Instead of manual picking, traveltimes of the first arrivals are estimated from the true model Gtrue by a fast marching eikonal solver available in the Madagascar package (Fomel et al. 2013). Then we transform the data into the Laplace–Fourier domain by the discrete Laplace–Fourier transform, and perform our inversion experiments. As conducted in Section 6.1.2, we update the initial baseline velocity model Gres1 to fit the noise-added data at T(0) with the non-target-oriented approach. After 200 model updates, the objective function is reduced to 16 per cent of the initial value. The double-difference inversion and the second step of the bootstrapping method are then conducted using this model as a starting model. We perform 50 iterations per full waveform inversion run. The objective functions steadily decrease, but the amount of reduction is much smaller than the noise-free inversions: With the 5-per-cent noise-added data, parallel inversion reduces the objective functions (⁠|$E_{p}^{(0)}$| and |$E_{p}^{(1)}$|⁠) by 26 percentiles, while double-difference and bootstrapping inversions decrease their respective objective functions (Ed and |$E^{(1)}_{b}$|⁠) by 29 and 42 percentiles. When conducting the inversions using the 10 per cent noise-added data, the relative reduction in the objective functions is approximately a half of the inversions of the 5-per-cent noise-added data. We display the final inverted velocity changes in Fig. 17. Adding Gaussian non-repeatable noise degrades the quality of the inverted models as expected from the signal-to-noise ratios of time-lapse waveforms. The parallel inversion method fails to recover a reasonable model of velocity changes regardless of the signal-to-noise ratios. The velocity increase in the top layer can be distinguished well by the double-difference and bootstrapping methods especially with the 5-per cent noise-added data. The bootstrapping method estimates the magnitude of the velocity increase slightly better than the double-difference method, but the discrepancies are marginal. The inversions of the 10-per cent noise data suffer from the increased level of artefacts that appear random, but the velocity increase in the top layer is delineated as coherent features. Figure 17. Open in new tabDownload slide Effects of non-repeatable noise: Velocity change estimated from seismic records with (a,c,e) 5 per cent and (b,d,f) 10 per cent uncorrelated non-repeatable noise using (a,b) parallel, (c,d) double-difference and (e,f) bootstrapping inversion methods. Model Gres1 is used as the initial baseline velocity model. Figure 17. Open in new tabDownload slide Effects of non-repeatable noise: Velocity change estimated from seismic records with (a,c,e) 5 per cent and (b,d,f) 10 per cent uncorrelated non-repeatable noise using (a,b) parallel, (c,d) double-difference and (e,f) bootstrapping inversion methods. Model Gres1 is used as the initial baseline velocity model. 6.3 Effects of changes in locations of repeating events We evaluate the feasibility of the parallel, double-difference and bootstrapping methods to handle non-repeatable source parameters. We examine effects of changes in depths of event locations, and further thorough investigations into other parameters including source wavelets and mechanisms, and also non-repeatable receiver parameters are the subject of ongoing research. In this paper, we assume that the source locations change over time, but we know the locations accurately. To conduct the evaluation of the effect of |$\mathbf {x}_s^{(0)}\ne \mathbf {x}_s^{(1)}$|⁠, we assume that the receiver parameters are repeatable (⁠|$\mathbf {x}_r=\mathbf {x}_r^{(0)}=\mathbf {x}_r^{(1)}$|⁠), and shift all three sources downwards by a constant value. The vertical shift is either 12.5 or 25.0 m. The 12.5-m shift corresponds to approximately 2 per cent of the wavelength at the dominant frequency, and causes a maximum of 2.5-ms delay in the first P-wave arrivals. The 25.0-m shift doubles the first arrival delay. The amount of the delays are of the same order as the actual time advances due to the velocity increase but in the opposite sign, that is, −5 ms. To see the effects further, we calculate the ‘observed’ time-lapse data change |$\Delta \mathbf {d}^{(1,0)} [\mathbf {x}_r,\mathbf {x}_s^{(0)},\mathbf {x}_s^{(1)}]$| by subtracting the observed record |$\mathbf {d}^{(0)}[\mathbf {x}_r,\mathbf {x}_s^{(0)}]$| at T(0) from the record |$\mathbf {d}^{(1)}[\mathbf {x}_r,\mathbf {x}_s^{(1)}]$| at T(1) as in eq. (16), and display the pseudo observed time-lapse data change for the repeating events at SA in Fig. 18. As described in Section 2, the pseudo time-lapse waveform changes are quite different from the true time-lapse waveforms depicted in Fig. 2 due to the non-repeatable source locations. Longer distances between the events and the sensors at T(1) delay the first arrivals of the record d(1) at all sensors when compared to the perfectly repeating case. As a result, the advance in first arrivals observed for the perfectly repeating data has diminished at the sensors at distances between 6 and 17 km, and arrival times are now delayed at sensors at distances between 0 and 6 km. Figure 18. Open in new tabDownload slide Effects of repeatability of source event locations: (a) ‘Time-lapse’ waveforms computed by subtracting the record of SA(0) from that of SA(1). The depths of the event SA(1) are (a,c,d) 5.0125 km (12.5 m deeper than the depth of SA(0)). Waveforms are scaled by a factor of 10, and displayed with the same colour scale as in Fig. 2. (b) Seismic traces recorded at sensors located at distances of 4, 8 and 12 km for the repeating events SA(0) and SA(1). The red line indicates the data observed at T(0), the blue line indicates the data at T(1), and the green line indicates ‘time-lapse’ waveforms (same as (a), but without scaling). The corresponding waveforms from the perfectly repeating events are displayed in Fig. 2. Note that no time-damping is applied to the waveforms. Figure 18. Open in new tabDownload slide Effects of repeatability of source event locations: (a) ‘Time-lapse’ waveforms computed by subtracting the record of SA(0) from that of SA(1). The depths of the event SA(1) are (a,c,d) 5.0125 km (12.5 m deeper than the depth of SA(0)). Waveforms are scaled by a factor of 10, and displayed with the same colour scale as in Fig. 2. (b) Seismic traces recorded at sensors located at distances of 4, 8 and 12 km for the repeating events SA(0) and SA(1). The red line indicates the data observed at T(0), the blue line indicates the data at T(1), and the green line indicates ‘time-lapse’ waveforms (same as (a), but without scaling). The corresponding waveforms from the perfectly repeating events are displayed in Fig. 2. Note that no time-damping is applied to the waveforms. Using Model Gres1 as an initial model, we conduct parallel, double-difference and bootstrapping inversions following the strategies used for the perfectly repeating events (Section 6.1). As the data at T(0) are the same as that of the perfectly repeating data tests in Section 6.1, we use the inverted velocity model G|$^{\prime }_{\mathrm{res1}}$| in Fig. 8 as an input for double-difference inversion, and the second stage of bootstrapping inversion. The double-difference method is performed by assuming that seismic events all occurred at SA(0), SB(0), and SC(0) (i.e. |$\hat{\mathbf {x}}_s=\mathbf {x}_s^{(0)}$| in eq. 26), and thus affected by location errors as seen in Fig. 18. The objective functions converge well, decreased to values nearly identical to those during the corresponding inversions of the seismic records of the perfectly repeating events in Section 6.1. We display the estimated velocity changes in Fig. 19. The imperfect repeatability of the event locations degrades the quality of the velocity change estimation of the all three inversion methods, and the effects become more significant as the vertical shift of the event locations increases. Parallel and double-difference inversions are severely affected by the discrepancies in the event location between T(0) and T(1) (Figs 19a–d). In the parallel inversion results shown in Figs 19(a) and (b), while we may be able to interpret the velocity increase in the top layer, the inverted velocity changes are more heavily contaminated by artefacts, when compared to the parallel inversion results obtained by using the perfectly repeating events in Fig. 10(a). The double-difference results (Figs 19c and d) delineate the top layer distinguishably, even though the assumption is not completely satisfied, and the ‘time-lapse’ waveform changes noticeably different from the true time-lapse wavefields. However. we observe velocity increase in the wide area of the model at depths greater than 3 km, which is introduced to compensate for the incorrect depth assumed for the event at T(1) (shallower than the actual depth). Additional artefacts are also present especially in the vicinity of the event locations. The bootstrapping method provides the best velocity change estimates (Figs 19e and f), least contaminated by the artefacts, as the method does not heavily rely on the repeatability of the event. However the magnitude of the inverted velocity changes becomes further underestimated when compared to the inversion tests using the perfectly repeating events (Fig. 10g). Figure 19. Open in new tabDownload slide Effects of repeatability of source event locations: the depths of the events at T(1) (SA(1), SB(1) and SC(1)) are (a,c,d) 5.0125 km (12.5 m deeper than the depth of SA(0), SB(0) and SC(0)) and (b,d,f) 5.025 km (25 m deeper than the depth of SA(0), SB(0) and SC(0)). (a,b) Parallel, (c,d) double-difference and (e,f) bootstrapping inversion results. Model Gres1 is used as the initial baseline velocity model. Note that double-difference inversion is conducted assuming that the events at both T(0) and T(1)are excited at the depth of 5 km. Figure 19. Open in new tabDownload slide Effects of repeatability of source event locations: the depths of the events at T(1) (SA(1), SB(1) and SC(1)) are (a,c,d) 5.0125 km (12.5 m deeper than the depth of SA(0), SB(0) and SC(0)) and (b,d,f) 5.025 km (25 m deeper than the depth of SA(0), SB(0) and SC(0)). (a,b) Parallel, (c,d) double-difference and (e,f) bootstrapping inversion results. Model Gres1 is used as the initial baseline velocity model. Note that double-difference inversion is conducted assuming that the events at both T(0) and T(1)are excited at the depth of 5 km. 6.4 Effects of quality of an initial velocity model Finally, we examine the effects of quality of an initial baseline velocity model by using Model Gres2 displayed in Fig. 1(d). As in Section 6.1.1, we first analyse the discrepancies between synthetic and predicted waveforms at T(0). Fig. 20 presents a synthetic seismic record u(0)[mo] at T(0) from Model Gres2, and its difference δd(0) from the observed data. When compared to those from Model Gres1 (Figs 2a and b), we observe larger amplitude discrepancies in the waveforms between observed and synthetic records. In addition, the first arrivals show clear phase shifts especially at far offsets. This indicates larger low-frequency residuals than Model Gres1, and cycle skipping at high frequencies. Also, the synthetic waveforms do not contain scattered arrivals observed in Figs 2(a) and (b), due to the greater smoothness of the model. Figure 20. Open in new tabDownload slide Effects of error in initial velocity models: (a) Predicted and (b) residual waveforms for SA(0) computed from Model Gres2 (used as a starting model for the inversion tests in Section 6.4). Waveforms in (a,b) are displayed with the same colour scale as in Fig. 2, while the residual and time-lapse waveforms in (b) are scaled by a factor of 10 for display purposes. (c) Seismic traces of (a,b) recorded at sensors located at distances of 4, 8 and 12 km. The red lines indicate the data observed at T(0) (same as Fig. 2a), the blue lines indicate the synthetic waveforms from Gres2 (same as (a)), and the green lines indicate residual waveforms (same as (b), but without scaling). No time-damping is applied to the waveforms. Figure 20. Open in new tabDownload slide Effects of error in initial velocity models: (a) Predicted and (b) residual waveforms for SA(0) computed from Model Gres2 (used as a starting model for the inversion tests in Section 6.4). Waveforms in (a,b) are displayed with the same colour scale as in Fig. 2, while the residual and time-lapse waveforms in (b) are scaled by a factor of 10 for display purposes. (c) Seismic traces of (a,b) recorded at sensors located at distances of 4, 8 and 12 km. The red lines indicate the data observed at T(0) (same as Fig. 2a), the blue lines indicate the synthetic waveforms from Gres2 (same as (a)), and the green lines indicate residual waveforms (same as (b), but without scaling). No time-damping is applied to the waveforms. Adopting the inversion parameters of the experiments for Model Gres1 in Section 6.1, we perform inversion tests to compare the parallel, double-difference, and bootstrapping methods. The inversion runs do not converge as well as in the Model-Gres1 inversion tests. The objective functions are reduced to 2–10 per cent of the value of the start of the inversion runs, both during the initial baseline velocity inversion, and during the parallel, double-difference and bootstrapping inversion runs. Effects of the poor convergence are apparent in the inverted time-lapse velocity changes displayed in Fig. 21: Parallel inversion does not estimate the time-lapse velocity change well. The top layer is imaged when using the double-difference and bootstrapping methods, but the magnitudes of the estimated velocity changes are much smaller than the true amount of velocity increase in the top layer and also than the inverted velocity changes from Model Gres1. The inversion results also contain substantial amount of both coherent and incoherent artefacts mapped due to the poor illumination to recover the background model. We observe an artificial coherent velocity increase close to the location of the bottom layer, and it is not trivial to distinguish the artefacts from the true velocity increase in the top layer. Figure 21. Open in new tabDownload slide Effects of error in initial velocity models: the inverted velocity changes when using Model Gres2 as an initial baseline velocity model. The velocity changes are estimated by (a) parallel, (b) double-difference, and (c) bootstrapping inversion methods. Figure 21. Open in new tabDownload slide Effects of error in initial velocity models: the inverted velocity changes when using Model Gres2 as an initial baseline velocity model. The velocity changes are estimated by (a) parallel, (b) double-difference, and (c) bootstrapping inversion methods. 7 DISCUSSION Our series of computational experiments demonstrates that full waveform inversion can accurately estimate small-scale velocity changes from dense surface sensor-array data sets using a limited number of sparse repeating seismic events. Our realistic test earth model is based on an actual reservoir, and contains a significant amount of heterogeneity, which generates realistic complex wavefields. Our time-lapse earth model features a the velocity increase over time of approximately 7–15 per cent of the background velocity model confined in two thin reservoir layers approximately 50 meters thick which corresponds to 0.83λ where λ is the wavelength at the dominant frequency. The small-scale time-lapse velocity changes cause a maximum time shift of only 5 ms in the first-arrival waveforms, approximately 0.025T where T is the dominant period. Such small data changes push the limits of full waveform inversion, especially with respect to its non-uniqueness, when confronted with complex wavefields, non-ideal source distributions, imperfect knowledge of the baseline velocity model, non-repeatable data noise, and non-repeatability of seismic events and sensor parameters (Section 2). We examine three time-lapse inversion strategies: the parallel, double-difference, and bootstrapping methods, in order to develop a robust inversion method capable of handling these challenging but realistic seismic monitoring conditions. We assume that reasonably accurate prior information is available (e.g. by separate data processing and analysis of relevant active-source seismic data sets) to constrain long- and mid-wavelength components of the time-lapse velocity models. This approach is plausible since an increasing number of passive monitoring projects are conducted along with controlled-source seismic surveys (e.g. Ben-Zion et al. 2015). Our computational experiments demonstrate that full waveform inversion can perform well even when the distribution of seismic events is sparse and not ideal. The vertical resolution of the image is approximately 100 m, much smaller than the size of the first Fresnel volume at the dominant frequency, |$\sqrt{\lambda L}\sim 2.4$| km where λ is the dominant wavelength, and L is the propagation distance. Our observation confirms the superior spatial resolution of full waveform inversion compared to that of traveltime tomography for the time-lapse case, as originally demonstrated for time-invariant cases (e.g. Kamei et al. 2013). This is primarily because full waveform inversion uses all waveform information including traveltimes, amplitudes, phase and waveform shapes rather than relying just on traveltimes, and properly handles seismic wave propagation at finite frequencies. The analysis of finite-frequency sensitivity kernels in Section 4 strongly suggests the importance of deploying dense sensors at wide aperture angles, and the advantage of having repeating earthquake locations occur near the localized velocity changes, which is often the case for induced microseismic events associated with manmade activities. Wide and dense sensor aperture provides excellent illumination at the vicinity of focal points (i.e. inside the reservoir), resulting in good recovery of the time-lapse velocity changes located just above the focal points, even with a small number of repeating earthquakes (Sections 5 and 6). However, the illumination is poor outside the reservoir (further from the seismic source), and we can no longer update the entire model in a geologically plausible manner (Section 6.1.2). For our repeating earthquake scenario, the effects of an inexact input baseline velocity model can be very substantial, since the background model error usually spans a large area. The associated artefacts (δbg) are then easily mapped as errors in the time-lapse velocity change estimation, and these velocity errors can be larger than the residual time-lapse velocity change error (δtl). A target-oriented approach, although accelerating the convergence of the inversion, potentially increases the artefacts by mapping velocity errors originating outside the target area into the target area. Using time damping can mitigate the artefacts by emphasizing first arrivals, and suppressing later scattered arrivals caused by small-scale heterogeneities outside the target area. Note that more careful time windowing may be required if a strong distinct arrival (e.g. reflection) is present but the corresponding geologic structure (e.g. a reflector outside of the target area) is not well represented in the initial model. Suppressing the inversion artefacts further is challenging, but all the three inversion methods show promising performances. Parallel inversion can reduce (or cancel out) the artefacts δbg by the subtraction of the two independently inverted models |$\mathbf {m}_p^{(0)}$| and |$\mathbf {m}_p^{(1)}$|⁠, implicitly relying on the assumption that the artefacts in |$\mathbf {m}_p^{(0)}$| and |$\mathbf {m}_p^{(1)}$| are similar or equal. Our computational experiments demonstrate that the parallel approach works well both when data shows very high signal-to-noise ratio, minimal non-repeatable noise, and highly repeatable source and sensor parameters , and when the initial velocity model is reasonably accurate as described in Section 2.3. Double-difference inversion mitigates the influence of the input baseline model error by directly using the time-lapse data difference Δd(1, 0), and focusing directly on updating velocity change estimates Δm(1, 0). While we confirm that the input baseline model does not have to be exact (relaxing the required condition), an additional baseline inversion is preferable prior to the application of double-difference inversion, so that the wave propagation paths are very accurate and δbg ∼ 0 (Section 6.1.4). As described in Section 2.4, data non-repeatability strongly affects the quality of the inverted time-lapse velocity change. In particular, double-difference inversion is vulnerable to source non-repeatability: small vertical shifts in the event locations are mapped erroneously as a coherent velocity increase. Bootstrapping (sequential) inversion accommodates both the input baseline model error and the non-repeatable source and sensor parameters by implicitly constraining the time-lapse data change through a two-stage inversion approach (Section 2.5). The first stage serves to reduce the model error δmbg, so that the residual for the second-stage inversion |$\delta \mathbf {d}^{(1)}_b$| is dominantly focused on estimating the time-lapse model change Δm(1, 0). The second stage of bootstrapping inversion acts similarly to double-difference inversion but without assuming data repeatability. Hence bootstrapping inversion can provide equivalent results to double-difference inversion when the data are perfectly repeatable, and more robust and accurate results when the repeatability becomes degraded. In the case that δbg > δtl, bootstrapping inversion can become unstable in a similar manner to parallel inversion, as exhibited in the late stage of the noise-free inversion test in Section 6.1.5). 8 CONCLUSIONS We demonstrate that full waveform inversion can accurately estimate a small time-lapse velocity change from surface sensor-array data recording a small set of repeated earthquakes. The spatial resolution of the velocity estimates is at sub-wavelength scale and is thus much higher than conventional ray-based time-lapse traveltime tomography. We identify important factors for successful estimation of time-lapse velocity changes using full waveform inversion as i) matching finite-frequency waveform data based on the wave equation, ii) having repeating events near velocity changes, iii) deploying seismic sensors over a wide aperture distance range, and iv) using an appropriate full waveform inversion strategy method that can accommodate initial baseline velocity model errors, non-repeatable data noise, and non-repeatable source and sensor parameters. Initial velocity errors and data/acquisition non-repeatability can create significant errors in the time-lapse inversion velocity change estimates. We develop, test and evaluate the parallel, double-difference and bootstrapping inversion methods for passive seismic recording scenarios involving repeating (microseismic) earthquake sources, which were initially developed for time-lapse controlled-source seismic monitoring. The parallel inversion method is the least stable, and most susceptible to the aforementioned issues. The double-difference and bootstrapping inversion methods generally provide much more robust and accurate inversion results than parallel inversion. The bootstrapping inversion method gives superior results compared to the double-difference method when seismic acquisition and noise conditions become increasingly non-repeatable. Acknowledgments We thank Don Sherlock and Joe Stefani of Chevron for their contributions in helping to develop the realistic velocity model used in this study. We thank the sponsors of the UWA:RM (University of Western Australia Reservoir Management) research consortium for providing partial financial support to conduct the research. We thank Jeffrey Shragge and Erdinc Saygin for their assistance in improving the manuscript before submission. We thank Jean Virieux and one anonymous reviewer for their review comments that have helped to improve our manuscript. REFERENCES Arts R. Eiken O. Chadwick A. Zweigel P. van der Meer L. Zinszner B., 2004 . Monitoring of CO2 injected at Sleipner using time-lapse seismic data , Energy 29 9–10 , 1383 – 1392 . Google Scholar Crossref Search ADS WorldCat Asnaashari A. Brossier R.J. Garambois S. Audebert F. Thore P. Virieux J., 2014 . Time-lapse seismic imaging using regularized full-waveform inversion with a prior model: which strategy? , Geophys. Prospect. 63 1 , 78 – 98 . Google Scholar Crossref Search ADS WorldCat Baltay A. Ide S. Prieto G.A. Beroza G.C., 2011 . Variability in earthquake stress drop and apparent stress , Geophys. Res. Lett. 38 6 , L06303, doi:10.1029/2011GL046698 . OpenURL Placeholder Text WorldCat Ben-Zion Y. et al. , 2015 . Basic data features and results from a spatially dense seismic array on the San Jacinto fault zone , Geophys. J. Int. 202 1 , 370 – 380 . Google Scholar Crossref Search ADS WorldCat Brenders A.J. Pratt R.G., 2007 . Full waveform tomography for lithospheric imaging: results from a blind test in a realistic crustal model , Geophys. J. Int. 168 1 , 133 – 151 . Google Scholar Crossref Search ADS WorldCat Brenguier F. Campillo M. Takeda T. Aoki Y. Shapiro N.M. Briand X. Emoto K. Miyake H., 2014 . Mapping pressurized volcanic fluids from induced crustal seismic velocity drops , Science 345 6192 , 80 – 82 . Google Scholar Crossref Search ADS PubMed WorldCat Bunks C. Saleck F.M. Zaleski S. Chavent G., 1995 . Multiscale seismic waveform inversion , Geophysics 60 5 , 1457 – 1473 . Google Scholar Crossref Search ADS WorldCat Daley T.M. Myer L.R. Peterson J.E. Majer E.L. Hoversten G.M., 2007 . Time-lapse crosswell seismic and VSP monitoring of injected CO2 in a brine aquifer , Environ. Geol. 54 8 , 1657 – 1665 . Google Scholar Crossref Search ADS WorldCat Fichtner A. Saygin E. Taymaz T. Cupillard P. Capdeville Y. Trampert J., 2013 . The deep structure of the North Anatolian Fault Zone , Earth planet. Sci. Lett. 373 C , 109 – 117 . Google Scholar Crossref Search ADS WorldCat Fomel S. Sava P. Vlad I. Liu Y. Bashkardin V., 2013 . Madagascar: open-source software project for multidimensional data analysis and reproducible computational experiments , J. Open Res. Softw. 1 1 , e8 , doi:org/10.5334/jors.ag . Google Scholar Crossref Search ADS WorldCat Gardner G.H.F. Gardner L.W. Gregory A.R., 1974 . Formation velocity and density—the diagnostic basics for stratigraphic traps , Geophysics 39 6 , 770 – 780 . Google Scholar Crossref Search ADS WorldCat Gunasekera R.C. Foulger G.R. Julian B.R., 2003 . Reservoir depletion at The Geysers geothermal area, California, shown by four-dimensional seismic tomography , J. geophys. Res. 108 B3 , 2134 , doi:10.1029/2001JB000638 . Google Scholar Crossref Search ADS WorldCat Hotovec-Ellis A.J. Vidale J.E. Gomberg J. Thelen W. Moran S.C., 2015 . Changes in seismic velocity during the first 14 months of the 2004–2008 eruption of Mount St. Helens, Washington , J. geophys. Res. 120 9 , 6226 – 6240 . Google Scholar Crossref Search ADS WorldCat Jo C.-H. Shin C. Suh J.H., 1996 . An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator , Geophysics 61 2 , 529 – 537 . Google Scholar Crossref Search ADS WorldCat Kamei R. Pratt R.G., 2013 . Inversion strategies for visco-acoustic waveform inversion , Geophys. J. Int. 194 2 , 859 – 884 . Google Scholar Crossref Search ADS WorldCat Kamei R. Pratt R.G. Tsuji T., 2013 . On acoustic waveform tomography of wide-angle OBS data–strategies for pre-conditioning and inversion , Geophys. J. Int. 194 2 , 1250 – 1280 . Google Scholar Crossref Search ADS WorldCat Kamei R. Pratt R.G. Tsuji T., 2014 . Misfit functionals in Laplace-Fourier domain waveform inversion, with application to wide-angle ocean bottom seismograph data , Geophys. Prospect. 62 5 , 1054 – 1074 . Google Scholar Crossref Search ADS WorldCat Lin F.-C. Li D. Clayton R.W. Hollis D., 2013 . High-resolution 3D shallow crustal structure in Long Beach, California: application of ambient noise tomography on a dense seismic array , Geophysics 78 4 , Q45 – Q56 . Google Scholar Crossref Search ADS WorldCat Lumley D.E. , 2001 . Time-lapse seismic reservoir monitoring , Geophysics 66 1 , 50 – 53 . Google Scholar Crossref Search ADS WorldCat Lumley D.E. Shragge J., 2013 . Advanced concepts in active and passive seismic monitoring using full wavefield techniques , ASEG Extended Abstracts 2013 1 , 1 – 4 . Google Scholar Crossref Search ADS WorldCat Mallick S. Frazer L.N., 1987 . Practical aspects of reflectivity modeling , Geophysics 52 10 , 1355 – 1364 . Google Scholar Crossref Search ADS WorldCat Martínez-Garzón P. Bohnhoff M. Kwiatek G. Dresen G., 2013 . Stress tensor changes related to fluid injection at The Geysers geothermal field, California , Geophys. Res. Lett. 40 11 , 2596 – 2601 . Google Scholar Crossref Search ADS WorldCat Mavko G. Mukerji T. Dvorkin J., 1998 . The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media Cambridge Univ. Press . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Maxwell S.C. Rutledge J. Jones R. Fehler M., 2010 . Petroleum reservoir characterization using downhole microseismic monitoring , Geophysics 75 5 , 75A129 – 75A137 . Google Scholar Crossref Search ADS WorldCat Nakata N. Snieder R., 2012 . Time-lapse change in anisotropy in Japan's near surface after the 2011 Tohoku-Oki earthquake , Geophys. Res. Lett. 39 11 , L11313 , doi:10.1029/2012GL051979 . Google Scholar Crossref Search ADS WorldCat Nakata N. Chang J.P. Lawrence J.F. Boue P., 2015 . Body wave extraction and tomography at Long Beach, California, with ambient-noise interferometry , J. geophys. Res. 120 1159 – 1173 . Google Scholar Crossref Search ADS WorldCat Nishikawa T. Ide S., 2015 . Background seismicity rate at subduction zones linked to slab-bending-related hydration , Geophys. Res. Lett. 42 17 , 7081 – 7089 . Google Scholar Crossref Search ADS WorldCat Operto S. Miniussi A. Brossier R.J. Combe L. Metivier L. Monteiller V. Ribodetti A. Virieux J., 2015 . Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation , Geophys. J. Int. 202 2 , 1362 – 1391 . Google Scholar Crossref Search ADS WorldCat Peng Z. Ben-Zion Y., 2005 . Spatiotemporal variations of crustal anisotropy from similar events in aftershocks of the 1999 M7.4 İzmit and M7.1 Düzce, Turkey, earthquake sequences , Geophys. J. Int. 160 3 , 1027 – 1043 . Google Scholar Crossref Search ADS WorldCat Polak E. Ribiére G., 1969 . Note sur la convergence de mèthodes de directions conjugueées , Revue Fr. Inf. Rech. Oper. 16-R1 35 – 43 . OpenURL Placeholder Text WorldCat Pratt R.G. , 1999 . Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale models , Geophysics 64 3 , 888 – 901 . Google Scholar Crossref Search ADS WorldCat Pratt R.G. Shin C. Hicks G., 1998 . Gauss-Newton and full Newton methods in frequency-space seismic waveform inversions , Geophys. J. Int. 133 2 , 341 – 362 . Google Scholar Crossref Search ADS WorldCat Queißer M. Singh S.C., 2012 . Full waveform inversion in the time lapse mode applied to CO2 storage at Sleipner , Geophys. Prospect. 61 3 , 537 – 555 . Google Scholar Crossref Search ADS WorldCat Raknes E.B. Arntsen B., 2014 . Time-lapse full-waveform inversion of limited-offset seismic data using a local migration regularization , Geophysics 79 3 , WA117 – WA128 . Google Scholar Crossref Search ADS WorldCat Rickett J.E. Lumley D.E., 2001 . Cross-equalization data processing for time-lapse seismic reservoir monitoring: a case study from the Gulf of Mexico , Geophysics 66 4 , 1015 – 1025 . Google Scholar Crossref Search ADS WorldCat Roach L.A.N. White D.J. Roberts B., 2015 . Assessment of 4D seismic repeatability and CO2 detection limits using a sparse permanent land array at the Aquistore CO2 storage site , Geophysics 80 2 , WA1 – WA13 . Google Scholar Crossref Search ADS WorldCat Roecker S. Baker B. McLaughlin J., 2010 . A finite-difference algorithm for full waveform teleseismic tomography , Geophys. J. Int. 181 2 , 1017 – 1040 . OpenURL Placeholder Text WorldCat Saffer D.M. Tobin H.J., 2011 . Hydrogeology and mechanics of subduction zone forearcs: fluid flow and pore pressure , Annu. Rev. Earth planet. Sci. 39 1 , 157 – 186 . Google Scholar Crossref Search ADS WorldCat Sawazaki K. Sato H. Nakahara H. Nishimura T., 2009 . Time-lapse changes of seismic velocity in the shallow ground caused by strong ground motion shock of the 2000 Western-Tottori Earthquake, Japan, as revealed from coda deconvolution analysis , Bull. seism. Soc. Am. 99 1 , 352 – 366 . Google Scholar Crossref Search ADS WorldCat Schaff D.P. Beroza G.C., 2004 . Coseismic and postseismic velocity changes measured by repeating earthquakes , J. geophys. Res. 109 B10 , B10302 , doi:10.1029/2004JB003011 . Google Scholar Crossref Search ADS WorldCat Shin C. Cha Y.H., 2009 . Waveform inversion in the Laplace-Fourier domains , Geophys. J. Int. 177 3 , 1067 – 1079 . Google Scholar Crossref Search ADS WorldCat Shragge J. Lumley D.E., 2013 . Time-lapse wave-equation migration velocity analysis , Geophysics 78 2 , S69 – S79 . Google Scholar Crossref Search ADS WorldCat Shulakova V. Pevzner R. Dupuis J.C. Urosevic M. Tertyshnikov K. Lumley D.E. Gurevich B., 2014 . Burying receivers for improved time-lapse seismic repeatability: CO2CRC Otway field experiment , Geophys. Prospect. 63 1 , 55 – 69 . Google Scholar Crossref Search ADS WorldCat Sirgue L. Pratt R.G., 2004 . Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies , Geophysics 69 1 , 231 – 248 . Google Scholar Crossref Search ADS WorldCat Tape C. Liu Q. Maggi A. Tromp J., 2009 . Adjoint tomography of the Southern California Crust , Science 325 5943 , 988 – 992 . Google Scholar Crossref Search ADS PubMed WorldCat Tarantola A. , 1984 . Inversion of seismic reflection data in the acoustic approximation , Geophysics 49 8 , 1259 – 1266 . Google Scholar Crossref Search ADS WorldCat Tarantola A. Valette B., 1982 . Inverse Problems = Quest for Information , Geophysics 50 159 – 170 . OpenURL Placeholder Text WorldCat Tkalčić H. Young M. Bodin T. Ngo S. Sambridge M., 2013 . The shuffling rotation of the Earth's inner core revealed by earthquake doublets , Nat. Geosci. 6 6 , 497 – 502 . Google Scholar Crossref Search ADS WorldCat Tsuji T. Kamei R. Pratt R.G., 2014 . Pore pressure distribution of a mega-splay fault system in the Nankai Trough subduction zone: insight into up-dip extent of the seismogenic zone , Earth planet. Sci. Lett. 396 C , 165 – 178 . Google Scholar Crossref Search ADS WorldCat Waldhauser F. Ellsworth W.L., 2000 . A double-difference earthquake location algorithm: method and application to the Northern Hayward Fault, California , Bull. seism. Soc. Am. 90 6 , 1353 – 1368 . Google Scholar Crossref Search ADS WorldCat Watanabe T. Shimizu S. Asakawa E. Matsuoka T., 2004 . Differential waveform tomography for time-lapse crosswell seismic data with application to gas hydrate production monitoring , SEG Technical Program Expanded Abstracts 2004 pp. 2323 – 2326 , Society of Exploration Geophysicists . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Williamson P.R. Worthington M.H., 1993 . Resolution limits in ray tomography due to wave behavior: numerical experiments , Geophysics 58 5 , 727 – 735 . Google Scholar Crossref Search ADS WorldCat Woodward M.J. , 1992 . Wave-equation tomography , Geophysics 57 1 , 15 – 26 . Google Scholar Crossref Search ADS WorldCat Wu R.-S. Toksoz M.N., 1987 . Diffraction tomography and multisource holography applied to seismic imaging , Geophysics 52 1 , 11 – 25 . Google Scholar Crossref Search ADS WorldCat Zhang H. Thurber C.H., 2003 . Double-difference tomography: the method and its application to the Hayward Fault, California , Bull. seism. Soc. Am. 93 5 , 1875 – 1889 . Google Scholar Crossref Search ADS WorldCat Zhang H. Thurber C.H., 2006 . Development and applications of double-difference seismic tomography , Pure appl. Geophys. 163 2/3 , 373 – 403 . Google Scholar Crossref Search ADS WorldCat Zhang Z. Huang L., 2013 . Double-difference elastic-waveform inversion with prior information for time-lapse monitoring , Geophysics 78 6 , R259 – R273 . Google Scholar Crossref Search ADS WorldCat Zhang F. Juhlin C. Ivandic M. Lüth S., 2013 . Application of seismic full waveform inversion to monitor CO2 injection: modelling and a real data example from the Ketzin site, Germany , Geophys. Prospect. 61 284 – 299 . Google Scholar Crossref Search ADS WorldCat Zhu H. Tromp J., 2013 . Mapping tectonic deformation in the crust and upper mantle beneath Europe and the North Atlantic Ocean , Science 341 6148 , 871 – 875 . Google Scholar Crossref Search ADS PubMed WorldCat Zoback M.D. , 2010 . Reservoir Geomechanics Cambridge Univ. Press . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC APPENDIX A: LAPLACE–FOURIER FULL WAVEFORM INVERSION Our implementation of acoustic Laplace–Fourier-domain full waveform inversion follows Kamei et al. (2013, 2014). We briefly review the theory here for time-invariant full waveform inversion, since time-lapse full waveform inversion is readily expandable as seen below. The readers are encouraged to see Kamei et al. (2013, 2014) for the details of the general Laplace–Fourier approach. The Laplace–Fourier transform of the time-domain wavefield u(x; t) is \begin{eqnarray} u(\mathbf {x};\omega ,s) = \int _{-\infty }^{\infty }u(\mathbf {x};t)\exp \left[-st\right]\exp \left[i\omega t\right] dt, \end{eqnarray} (A1) where ω is the real-valued angular frequency, and s is the real-valued Laplace constant (Shin & Cha 2009). Eq. (A1) suggests that the Laplace–Fourier wavefield can be interpreted as the Fourier transform of the original time-domain wavefield filtered by a time-damping function, exp [−st]. As using an appropriate value of s suppresses later arrivals, the Laplace–Fourier wavefield is originally employed to avoid time-wrapping of wavefields (Mallick & Frazer 1987; Sirgue & Pratt 2004), and then extended as an implementation of the multiscale approach (e.g. Kamei et al.2013). We recast eq. (A1) by using a complex frequency defined as \begin{equation} \Omega =\omega +is, \end{equation} (A2) such that \begin{equation} u(\mathbf {x};\Omega ) = \int _{-\infty }^{\infty }u(\mathbf {x};t)\exp \left[i\Omega t\right] dt. \end{equation} (A3) The time-invariant acoustic wave equation in the Laplace–Fourier domain at a complex frequency Ωl is \begin{eqnarray} \nabla \left(\frac{1}{\rho (\mathbf {x})}\nabla u(\mathbf {x},\mathbf {x}_{s_i};\Omega _l)\right)+\frac{\Omega _l^{2}}{\rho (\mathbf {x})c(\mathbf {x})^{2}}u(\mathbf {x},\mathbf {x}_{s_i};\Omega _l)=f(\mathbf {x}_{s_i};\Omega _l), \nonumber\\ \end{eqnarray} (A4) where u is the predicted pressure field, ρ is the density, c is the P-wave velocity, and  f  is the source term. We solve the wave equation by the finite-difference method following (Jo et al. 1996), and apply the absorbing boundary condition using the perfectly matching layer method (Roecker et al. 2010). We parametrize the subsurface by P-wave slowness, and estimate an optimum model by minimizing the objective function in eq. (3) such that \begin{eqnarray} E = \Vert \delta \mathbf {d}\Vert ^2 \end{eqnarray} (A5) \begin{eqnarray} \quad = \frac{1}{2}\sum _{i,j,l}\delta d(\mathbf {x}_{r_j},\mathbf {x}_{s_i};\Omega _l)\delta d(\mathbf {x}_{r_j},\mathbf {x}_{s_i};\Omega _l)^* \end{eqnarray} (A6) where * indicates the complex conjugate. δd is the data residual defined as \begin{equation} \delta d (\mathbf {x}_{r_j},\mathbf {x}_{s_i};\Omega _l)=u(\mathbf {x}_{r_j},\mathbf {x}_{s_i};\Omega _l)-d(\mathbf {x}_{r_j},\mathbf {x}_{s_i};\Omega _l), \end{equation} (A7) where |$u(\mathbf {x}_{r_j},\mathbf {x}_{s_i};\Omega _l)$| is the predicted waveforms obtained by solving eq. (A4), and |$d(\mathbf {x}_{r_j},\mathbf {x}_{s_i};\Omega _l)$| is the observed waveforms. We employ the nonlinear conjugate gradient method (Polak & Ribiére 1969) to optimize E. At each iteration, forward modelling is re-calculated with the updated model parameters and source wavelet, and the conjugate gradient is constructed from a set of the gradients of the objective function. The gradient ∇E is derived by the adjoint-state method (Tarantola 1984) following Pratt et al. (1998) and Kamei et al. (2013, 2014), and \begin{eqnarray} \nabla _{m(\mathbf {x})} E & = & \Re \left\lbrace \sum _{i,j,l} \frac{\partial u (\mathbf {x}_{r_j},\mathbf {x}_{s_i},\Omega _l) }{\partial m} \delta d^*(\mathbf {x}_{r_j},\mathbf {x}_{s_i},\Omega _l)\right\rbrace , \end{eqnarray} (A8) \begin{eqnarray} \qquad\quad= \Re \left\lbrace \sum _{i,j,l} -\frac{\Omega _l^2}{c(\mathbf {x}) }u(\mathbf {x},\mathbf {x}_{s_i},\Omega _l)v(\mathbf {x},\mathbf {x}_{r_j},\Omega _l)\right\rbrace , \end{eqnarray} (A9) where ℜ is the real part of a complex function, m(x) is the slowness at x, |$\frac{\partial u }{\partial m}$| is the Fréchet sensitivity kernel that can be obtained by solving the wave equation with a virtual source term |$-\frac{\Omega _l^2}{c(\mathbf {x})} u(\mathbf {x})$| , and v(x, xr; Ωl) is the backpropagated wavefield obtained by solving the wave equation with a source term |$\delta d^*(\mathbf {x}_{r_j})$|⁠. Because of computational efficiency, we implement gradient computation following eq. (A9). By using Green's function, we can denote the Fréchet sensitivity kernel as \begin{eqnarray} \frac{\partial u(\mathbf {x}_{r_j}, \mathbf {x}_{s_i},\Omega _l)}{\partial m(\mathbf {x})} & = & - \frac{\Omega _l^2}{c(\mathbf {x})} G(\mathbf {x},\mathbf {x}_{s_i},\Omega _l) f(\mathbf {x}_{s_i}, \Omega_l) G (\mathbf {x}_{r_j},\mathbf {x}, \Omega _l),\nonumber\\ \end{eqnarray} (A10) \begin{eqnarray} \qquad\qquad\qquad\,\,= K(\mathbf {x},\mathbf {x}_{r_j},\mathbf {x}_{s_i},\Omega _l) f(\mathbf {x}_{s_i}, \Omega_l) \end{eqnarray} (A11) where K is the source-free sensitivity kernel \begin{equation} K(\mathbf {x},\mathbf {x}_{r_j},\mathbf {x}_{s_i},\Omega _l) = -\frac{\Omega _l^2}{c(\mathbf {x}) }G(\mathbf {x},\mathbf {x}_{s_i},\Omega _l)G(\mathbf {x},\mathbf {x}_{r_j},\Omega _l), \end{equation} (A12) equivalent to eq. (40). Then we can write the gradient as \begin{eqnarray} \nabla _{m(\mathbf {x})} E = \Re \left\lbrace \sum _{i,j,l} K(\mathbf {x},\mathbf {x}_{r_j},\mathbf {x}_{s_i},\Omega _l)f(\mathbf {x}_{s_i},\Omega _l)\delta d^* (\mathbf {x}_{r_j},\mathbf {x}_{s_i},\Omega _l)\right\rbrace , \nonumber\\ \end{eqnarray} (A13) equivalent to eq. (41). The above derivation for the time-invariant medium can be readily expanded to the time-variant medium by defining the time-lapse wave equation \begin{eqnarray} {\nabla \left(\frac{1}{\rho ^{(k)}(\mathbf {x})}\nabla u^{(k)}(\mathbf {x},\mathbf {x}^{(k)}_{s_i},\Omega _l)\right)+\frac{\Omega _l^{2}}{\rho ^{(k)}(\mathbf {x})c^{(k)}(\mathbf {x})^{2}}u^{(k)}(\mathbf {x},\mathbf {x}^{(k)}_{s_i},\Omega _l)}\nonumber\\ =f^{(k)}(\mathbf {x}^{(k)}_{s_i},\Omega _l). \end{eqnarray} (A14) APPENDIX B: ERROR OF PARALLEL INVERSION To evaluate the error in the time-lapse model estimate |$\Delta \mathbf {m}_p^{(1,0)}$|⁠, we first define the error of the estimated models (⁠|$\mathbf {m}_p^{(0)}$| and |$\mathbf {m}_p^{(1)}$|⁠), |$\delta \mathbf {m}_p^{(0)}$| and |$\delta \mathbf {m}_p^{(1)}$| as \begin{eqnarray} \delta \mathbf {m}_p^{(0)} \left[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right] = \mathbf {m}^{(0)}_{\text{true}} - \mathbf {m}_p^{(0)}\left[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!, \end{eqnarray} (B1) \begin{eqnarray} \delta \mathbf {m}_p^{(1)} \left[\mathbf {d}^{(1)},\mathbf {m}_o,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right] = \mathbf {m}^{(1)}_{\text{true}} - \mathbf {m}_p^{(1)}\left[\mathbf {d}^{(1)},\mathbf {m}_o,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]\!. \end{eqnarray} (B2) Based on eq. (8), we can decompose the model errors into two terms, \begin{eqnarray} \delta \mathbf {m}_p^{(0)} \left[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right] = \delta \mathbf {m}^{bg}_p \left[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right] \nonumber\\ &&+\,\, \delta \Delta \mathbf {m}_p^{(0)} \left[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right] \!, \end{eqnarray} (B3) \begin{eqnarray} \delta \mathbf {m}_p^{(1)} \left[\mathbf {d}^{(1)},\mathbf {m}_o,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right] = \delta \mathbf {m}^{bg}_p \left[\mathbf {d}^{(1)},\mathbf {m}_o,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]\nonumber\\ &&+\,\, \delta \Delta \mathbf {m}_p^{(1)} \left[\mathbf {d}^{(1)},\mathbf {m}_o,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right] \!, \end{eqnarray} (B4) where |$\delta \mathbf {m}^{bg}_p$| is the error corresponding to the time-invariant (background) model, and |$\delta \Delta \mathbf {m}_p^{(0)}$| and |$\delta \Delta \mathbf {m}_p^{(1)}$| are time-variant components of the model at T(0) and T(1), respectively. As the two inversions are independent, |$\delta \mathbf {m}^{bg}_p[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}] \ne \delta \mathbf {m}^{bg}_p [\mathbf {d}^{(1)},\mathbf {m}_o,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}]$|⁠. Thus the inversion error in |$\Delta \mathbf {m}_p^{(1,0)}$|⁠, |$\delta \Delta \mathbf {m}_p^{(1,0)}$|⁠, is given by \begin{eqnarray} {\delta \Delta \mathbf {m}_{p}^{(1,0)} \left[\mathbf {d}^{(1)},\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]} \nonumber \\ && \quad= \Delta \mathbf {m}_{\text{true}}^{(1,0)} - \Delta \mathbf {m}_p^{(1,0)}\left[\mathbf {d}^{(1)},\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!, \end{eqnarray} (B5) \begin{eqnarray} \quad\,\,= \delta _p^{tl} + \delta _p^{bg}, \end{eqnarray} (B6) where \begin{eqnarray} \delta _p^{tl} & = & \delta \Delta \mathbf {m}_p^{(1)} \left[\mathbf {d}^{(1)},\mathbf {m}_o,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]- \delta \Delta \mathbf {m}_p^{(0)}\left[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!,\nonumber\\ \end{eqnarray} (B7) \begin{eqnarray} \delta _p^{bg} & = & \delta \mathbf {m}^{bg}_p \left[\mathbf {d}^{(1)},\mathbf {m}_o,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right] - \delta \mathbf {m}^{bg}_p\left[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!. \end{eqnarray} (B8) The above equations suggest that both the time-lapse model error |$\delta _p^{tl}$| and the background model error |$\delta _p^{bg}$| will be affected by repeatability of source, sensor parameters and data noise. As the repeatability increases, we expect that the two terms of the right hand of eq. (B7) (or eq. B8) become closer. The time-lapse model change error will be gradually reduced, but will not be zero, as the two independent inversions are conducted to fit different data sets (d(0) and d(1)). APPENDIX C: IMPLEMENTATION OF DOUBLE-DIFFERENCE INVERSION By following Watanabe et al. (2004), Zhang & Huang (2013) and Asnaashari et al. (2014), we rearrange eq. (28) using eq. (32) as \begin{eqnarray} \delta \Delta \mathbf {d}& = & \left(\mathbf {u}^{(1)}\left[\mathbf {m}_o+\Delta \mathbf {m}^{(1,0)},\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right]- \mathbf {u}^{(0)}\left[\mathbf {m}_o,\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right]\right)\nonumber\\ &&-\,\, \Delta \mathbf {d}^{(1,0)}\left[\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!, \end{eqnarray} (C1) \begin{eqnarray} \hphantom{\delta \Delta \mathbf {d}}& = & \mathbf {u}^{(1)}\left[\mathbf {m}_o+\Delta \mathbf {m}^{(1,0)},\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right] - \left(\mathbf {u}^{(0)}\left[\mathbf {m}_o,\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right]\right.\nonumber\\ &&\left.+\,\, \Delta \mathbf {d}^{(1,0)}\left[\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\right)\!, \end{eqnarray} (C2) \begin{eqnarray} \hphantom{\delta \Delta \mathbf {d}}& =& \mathbf {u}^{(1)}\left[\mathbf {m}_o+\Delta \mathbf {m}^{(1,0)},\mathbf {x}_r,\mathbf {x}_s\right]\nonumber\\ &&-\,\, \mathbf {d}_{d}\left[\mathbf {m}_o,\Delta \mathbf {d}^{(1,0)},\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!, \end{eqnarray} (C3) where dd is the ‘composite’ wavefield (Asnaashari et al. 2014): \begin{eqnarray} &&{\mathbf {d}_{d}\left[\mathbf {m}_o,\Delta \mathbf {d}^{(1,0)},\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]}\nonumber\\ &&{\quad= \mathbf {u}^{(0)}\left[\mathbf {m}_o,\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right] + \Delta \mathbf {d}^{(1,0)}\left[\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!.} \end{eqnarray} (C4) Eq. (C3) indicates that we can estimate mo + Δm(1, 0) by minimizing the difference between the composite and predicted waveforms. Note that the composite waveforms are still not linear with respect to the model change, and thus the inversion is nonlinear. We recast the objective function in eq. (31) as \begin{equation} E_{d} = \left\Vert \delta \mathbf {d}_{d}\right\Vert ^2, \end{equation} (C5) where \begin{eqnarray} \delta \mathbf {d}_{d} = \mathbf {u}^{(1)}\left[\mathbf {m}_o+\Delta \mathbf {m}^{(1,0)},\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s\right] \nonumber\\ &&-\,\, \mathbf {d}_{d}\left[\mathbf {m}_o,\Delta \mathbf {d}^{(1,0)},\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!. \end{eqnarray} (C6) Based on eqs (C4)–(C6), we implement the double-difference approach as follows: First, we calculate the waveform differences Δd in eq. (16), and synthesize the waveforms u(0) from the model mo. We compute the composite waveforms dd following eq. (C4) and then use them as ‘observed’ waveforms. We obtain |$\mathbf {m}_{d}[\Delta \mathbf {d}^{(1,0)},\mathbf {m}_o,\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}]$| by minimizing the L2 norm of the data residuals in eq. (C3). The velocity change |$\Delta \mathbf {m}_{d}^{(1,0)}$| is then calculated by \begin{eqnarray} &&{\Delta \mathbf {m}_{d}^{(1,0)}\left[\Delta \mathbf {d}^{(1,0)},\mathbf {m}_o,\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]}\nonumber\\ &&{\quad=\mathbf {m}_{d}\left[\Delta \mathbf {d}^{(1,0)},\mathbf {m}_o,\hat{\mathbf {x}}_r,\hat{\mathbf {x}}_s,\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)},\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right] - \mathbf {m}_o.} \end{eqnarray} (C7) As the starting model for the inversion of the composite waveforms is mo, |$\Delta \mathbf {m}_d^{(1,0)}$| is the total model update during the inversion. Thus the double-difference method does not require the subtraction of the two inversion results, and explicitly utilizes time-lapse differences in data. The resulting velocity change can be more accurate than the parallel inversion method, when we satisfy the two assumptions made previously; perfect data repeatability, and an accurate baseline velocity model. APPENDIX D: ERROR OF BOOTSTRAPPING INVERSION To evaluate error in the bootstrapping estimation, we define errors in the estimated models (⁠|$\mathbf {m}_b^{(0)}$| and |$\mathbf {m}_b^{(1)}$|⁠), |$\delta \mathbf {m}_b^{(0)}$| and |$\delta \mathbf {m}_b^{(1)}$| \begin{eqnarray} \delta \mathbf {m}_b^{(0)}\left[\mathbf {d}^{(0)}, \mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right] & = & \mathbf {m}_{\text{true}}^{(0)}-\mathbf {m}_b^{(0)} \left[\mathbf {d}^{(0)}, \mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!,\nonumber\\ \end{eqnarray} (D1) \begin{eqnarray} \delta \mathbf {m}_b^{(1)}\left[\mathbf {d}^{(1)}, \mathbf {m}^{(0)}_{b},\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]& = & \mathbf {m}_{\text{true}}^{(1)} - \mathbf {m}_b^{(1)}\left[\mathbf {d}^{(1)}, \mathbf {m}^{(0)}_{b},\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]\!.\nonumber\\ \end{eqnarray} (D2) We decompose the errors into two terms as performed in the analysis of parallel inversion \begin{eqnarray} \delta \mathbf {m}_b^{(0)} \left[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right] = \delta \mathbf {m}^{bg}_b\left[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\nonumber\\ &&+\,\, \delta \Delta \mathbf {m}_b^{(0)}\left[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!. \end{eqnarray} (D3) \begin{eqnarray} \delta \mathbf {m}_b^{(1)} \left[\mathbf {d}^{(1)},\mathbf {m}_b^{(0)},\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right] = \delta \mathbf {m}^{bg}_b\left[\mathbf {d}^{(1)},\mathbf {m}^{(0)},\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]\nonumber\\ &&+\,\, \delta \Delta \mathbf {m}_b^{(1)}\left[\mathbf {d}^{(1)},\mathbf {m}_b^{(0)},\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right]\!. \nonumber\\ \end{eqnarray} (D4) Then the inversion error in the time-lapse model change estimate is \begin{eqnarray} \delta \Delta \mathbf {m}_b^{(1,0)} = \Delta \mathbf {m}_{\text{true}} - \Delta \mathbf {m}_b^{(1,0)}\!, \end{eqnarray} (D5) \begin{eqnarray} \qquad\quad\,\,\,= \delta _{b}^{tl} + \delta _b^{bg}, \end{eqnarray} (D6) where \begin{eqnarray} \delta _b^{tl} & = &\delta \Delta \mathbf {m}_b^{(1)}\left[\mathbf {d}^{(1)},\mathbf {m}_b^{(0)},\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right] - \delta \Delta \mathbf {m}_b^{(0)}\left[\mathbf {d}^{(0)},\mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!,\nonumber\\ \end{eqnarray} (D7) \begin{eqnarray} \delta _b^{bg} & = & \delta \mathbf {m}_b \left[\mathbf {d}^{(1)},\mathbf {m}^{(0)},\mathbf {x}_r^{(1)},\mathbf {x}_s^{(1)}\right] -\delta \mathbf {m}^{bg}_b \left[\mathbf {d}^{(0)}, \mathbf {m}_o,\mathbf {x}_r^{(0)},\mathbf {x}_s^{(0)}\right]\!. \end{eqnarray} (D8) © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Full waveform inversion of repeating seismic events to estimate time-lapse velocity changes
JF - Geophysical Journal International
DO - 10.1093/gji/ggx057
DA - 2017-05-01
UR - https://www.deepdyve.com/lp/oxford-university-press/full-waveform-inversion-of-repeating-seismic-events-to-estimate-time-bnfCQsP5Sg
SP - 1239
VL - 209
IS - 2
DP - DeepDyve
ER -