TY - JOUR
AU1 - Yang,, Jidong
AU2 - Zhu,, Hejun
AB - SUMMARY Locating and monitoring passive seismic sources provides us important information for studying subsurface rock deformation, injected fluid migration, regional stress conditions as well as fault rupture mechanism. In this paper, we present a novel passive-source monitoring approach using vector-based elastic time-reversal imaging. By solving the elastic wave equation using observed multicomponent records as boundary conditions, we first compute back-propagated elastic wavefields in the subsurface. Then, we separate the extrapolated wavefields into compressional (P-wave) and shear (S-wave) modes using the vector Helmholtz decomposition. A zero-lag cross-correlation imaging condition is applied to the separated pure-mode vector wavefields to produce passive-source images. We compare imaging results using three implementations, that is dot-product, energy and power. Numerical experiments demonstrate that the power imaging condition gives us the highest resolution and is less sensitive to the presence of random noises. To capture the propagation of microseismic fracture and earthquake rupture, we modify the traditional zero-lag cross-correlation imaging condition by summing the multiplication of the separated P and Swavefields within local time windows, which enables us to capture the temporal and spatial evolution of earthquake rupture. 2-D and 3-D numerical examples demonstrate that the proposed method is capable of accurately locating point sources, as well as delineating dynamic propagation of hydraulic fracture and earthquake rupture. Numerical solutions, Computational seismology, Earthquake source observations 1 INTRODUCTION Locating and real-time monitoring passive sources are important for understanding subsurface structural conditions, as well as corresponding stress and stain states at a variety of scales. For instance, at a reservoir scale, hydraulic fracturing becomes a widely used technology to enhance oil recovery in unconventional gas and oil fields, which can induce a lot of non-neglectable microseismic events (McGarr 2014; McGarr et al. 2015; Weingarten et al. 2015). Analysing these microseismic data is helpful for delineating local faults and fracture trajectories, injected and in situ fluid migration, as well as rock deformation mechanisms. These information can be further used to predict new or reactivating failure regions (Baig & Urbancic 2010; Dai et al. 2016; Jonas et al. 2016; Ma et al. 2016). At crustal and global scales, identification and characterization of earthquake location and size can be used to guide effective earthquake response efforts, and estimate the potential impacts of subsequent devastating tsunamis (Kiser & Ishii 2017; Lay 2018). In addition, studying earthquake rupture processes enables us to better understand fault failure mechanisms. Therefore, it is important to develop accurate and efficient passive-source locating and monitoring approaches for investigating microseismicities and earthquakes. Despite significant advances in theoretical and computational seismology, most imaging methods for microseismic events and earthquakes still rely on picking and inverting P- and S-wave arrival times. This phase-picking-based approach is effective for simple events, that traveltimes can be separated and identified, and signal-to-noise ratio (SNR) of data is high. However, if multiple waves arrive at the same station and at the same time, it is difficult to identify the onset of individual events, especially with strong ambient noises (Kao & Shan 2004; Bose et al. 2009; Kummerow 2010; Maxwell et al. 2010; Hayles et al. 2012). Strong dependence of such locating methods on the picking accuracy of traveltimes also leads to a low degree of automation (Rutledge & Phillips 2003). Using traveltime differences between pairs of events or stations, the precision of relative locations can be improved for a cluster of seismic events (Zhou 1994; Waldhauser & Ellsworth 2000; Zhang & Thurber 2003). But this method is still subject to the same problem described above, that is the robustness of traveltime picking. An alternative strategy for imaging earthquake sources is the waveform-based time-reversal imaging method. Similar to reverse-time migration (RTM) used in reflection seismology (Claerbout 1971; Baysal et al. 1983; McMechan 1983), in this technology, a 4-D wavefield volume is first computed by solving the acoustic or elastic wave equation, using observed passive-source data as boundary conditions. Then, a designed imaging condition, such as the zero-lag cross-correlation, can be utilized to collapse the time axis of the 4-D cube, and produce an image for the passive source distribution. In comparison with the traveltime-based method, the time-reversal imaging approach has the following advantages. First, without picking arrival times for particular phases, it enables us to avoid the locating error caused by inaccurate phase onset identification, and also has a high level of automation (Douma & Snieder 2015). Second, with a relatively accurate background velocity model, the back-propagated wavefields are focused at the true source locations. Coherent staking of wavefields from different stations can suppress ambient noises and produce a high-SNR image (Nakata & Beroza 2016). Third, the time-reversal extrapolation enables us to reconstruct traveling wavefields in the subsurface. With specially designed time-lapse imaging strategies, we are able to produce a series of dynamic images for monitoring cracks and fault rupture (Ishii et al. 2005; Walker et al. 2005; White 2013). Fourth, instead of only relying on the traveltime information, the amplitudes and polarities of passive source images can be used to estimate earthquake magnitude, faulting mechanisms, as well as subsurface strain and stress states (Baig & Urbancic 2010; Kiser & Ishii 2017). The basic framework of using acoustic and elastic reverse-time imaging for transmission data to locate passive sources can be found in Artman et al. (2010). In an acoustic case, the passive source image can be obtained by summing the autocorrelation of the extrapolated wavefield over the time axis. By modifying the autocorrelation procedure with a multidimensional cross-correlation imaging condition, Nakata & Beroza (2016) developed a geometric-mean imaging condition to enhance the SNR of imaging results. In an elastic case, Artman et al. (2010) proposed to first separate P and Swaves, and then evaluated the elastic imaging condition. Later, Shang et al. (2012) applied this workflow to migrate teleseismic data to image crust and mantle structures. In Artman et al. (2010)’s method, the elastic wavefield decomposition is implemented by computing scalar and vector potential wavefields using divergence and curl operators (Yan & Sava 2008). As analysed by Zhu (2017) in elastic RTM for reflection data, there are several issues when applying this approach in practice. First, the separated scalar and vector potential wavefields do not have the same phases, amplitudes and physical units as the input extrapolated wavefields, leading to inaccurate migration results (Sun & McMechan 2001). Second, there are three-component PS images for 3-D problems, and they are difficult to be combined for structure interpretation. Third, using this type of wavefield decomposition and the conventional imaging condition leads to the polarity reversal problem for PS images. Without appropriate angular dependent corrections, summation of PS images over different sources usually yields nonconstructive stacking results, which degrades the image quality of final solutions (Du et al. 2012; Duan & Sava 2015). In the passive source imaging, this polarity reversal issue leads to asymmetric imaging results even for a single point source (Artman et al. 2010), which might be erroneously interpreted as the occurrence of multiple microseismic events. In reflection seismology, the state-of-the-art elastic wavefield decomposition is to produce pure-mode vector wavefields, which have the same phases and amplitudes as the input coupled elastic wavefields. This can be implemented by either using the vector Helmholtz decomposition (Zhang & McMechan 2010; Zhu 2017; Yang et al. 2017), or introducing auxiliary P- and S-wave equations (Zhang et al. 2007; Xiao & Leaney 2010; Wang & McMechan 2015). With the separated vector wavefields, a dot-product (Wang et al. 2016; Zhu 2017) or a modified dot-product (Yang et al. 2018a, b) imaging condition can be used to produce PP and PS reflectivity images. This vector-based elastic RTM workflow produces only one PS image for both 2-D and 3-D problems, and also enables us to avoid the polarity reversal issue. In this paper, we introduce the concept of vector-based elastic time-reversal imaging for locating and monitoring passive sources. Similar to reflection RTM, the backward propagated wavefields are first computed by solving the elastodynamic wave equation, using observed multicomponent records as boundary conditions. Then, P- and S-mode wavefields are separated using the vector Helmholtz decomposition (Zhu 2017) or an efficient phase- and amplitude-correction method (Yang et al. 2018b). Finally, a zero-lag cross-correlation imaging condition with three implementations, that is, dot-product, energy and power, is used to produce an image for passive-source distributions. Numerical tests demonstrate that the power-based imaging condition has the highest spatial resolution and is less sensitive to the presence of noises. To capture the dynamic behaviours of crack and rupture propagation, we modify the traditional zero-lag cross-correlation imaging condition by summing the separated P and Swavefields within local time windows, which produces a series of time-lapse images. Numerical experiments demonstrate that the proposed method enables us to accurately locate passive point sources, and characterize the dynamic propagation of hydraulic fractures and earthquake ruptures. The outline of this paper is as follows. First, we review the elastic wavefields decomposition strategy using the vector Helmholtz decomposition. Then, we present the method of vector-based elastic time-reversal seismic-source imaging. Finally, several 2-D and 3-D synthetic examples are used to test the feasibility and adaptability of the proposed method for locating and monitoring point sources, as well as fracture and rupture propagations. 2 METHODOLOGY 2.1 Review of elastic wavefield separation using the vector Helmholtz decomposition In isotropic elastic media, the Helmholtz decomposition theory (Morse & Feshbach 1946; Aki & Richards 1980) states that the vector displacement or particle velocity wavefields |$\mathbf {u}$| can be decomposed as \begin{eqnarray*} \mathbf {u} = \mathbf {u}^{P}+\mathbf {u}^{S}=\nabla \phi +\nabla \times \mathbf {\Psi }, \end{eqnarray*} (1) where |$\mathbf {u}^P$| and |$\mathbf {u}^{S}$| are the P and Swavefields, ϕ and |$\mathbf {\Psi }$| are their scalar and vector potentials, respectively. By introducing an auxiliary wavefield |$\mathbf {w}$| that satisfies the following Poisson’s equation \begin{eqnarray*} \nabla ^{2}\mathbf {w} = \mathbf {u}, \end{eqnarray*} (2) and using the identical relation for the Laplacian operator \begin{eqnarray*} \nabla ^{2}\mathbf {w} = \nabla \left(\nabla \cdot \mathbf {w}\right)-\nabla \times \nabla \times \mathbf {w}, \end{eqnarray*} (3) Zhu (2017) proposed the following wavefield separation formulae \begin{eqnarray*} \mathbf {u}^{P} = \nabla \left(\nabla \cdot \mathbf {w}\right), \quad \mathbf {u}^{S} = -\nabla \times \nabla \times \mathbf {w}, \end{eqnarray*} (4) where ∇, ∇ · and ∇ × are the gradient, divergence and curl operators, respectively. By using Fourier-based methods to solve eq. (2), this wavefield separation approach is equivalent to solving a linear system in the wavenumber domain as proposed by Zhang & McMechan (2010). Different from the divergence and curl operators proposed by Yan & Sava (2008), Zhu (2017)’s and Zhang & McMechan (2010)’s methods produce vector P and Swavefields, which have the same phases and amplitudes as the input coupled elastic wavefields. In terms of computational costs, these two vector wavefield separation methods require either a pair of forward and inverse Fourier transform or solving a vector Poisson’s equation at each time step during the wavefield extrapolation. It is still challenging for large-scale 3-D problems under current computational capability (Yang et al. 2017). Using smoothed migration velocity models in elastic RTM, Yang et al. (2018b) simplify Zhu (2017)’s method by correcting the phases and amplitudes of the source and receiver wavefields. This strategy enables us to avoid solving the Poisson’s equation, and therefore significantly reduce computational costs. Figs 1–3 illustrate a 2-D wavefield separation experiment using different methods for a homogeneous medium with P-wave velocity of 3 km s−1 and S-wave velocity of 1.5  km s−1. The multicomponent wavefields (Fig. 1) are computed by solving the elastic wave equation using the staggered-grid finite-difference scheme. A double-couple source is deployed in the centre of this model. The separation results using the divergence and curl operators are shown in Fig. 2. As expected, the outputs are scalar wavefields for this 2-D problem, and their phases and amplitudes are changed due to the usages of the divergence and curl operators. In addition, it is notable that the spatial radiation patterns of the separated wavefields are also changed. For instance, the Pwavefield (Fig. 2a) is neither the same as the input elastic wavefield, nor as that in the acoustic medium (homogeneous radiation). For the S wavefield (Fig. 2b), the quadrant number of the polarity variation reduces from six to four. In contrast, eq. (4) allows us to produce the decomposed wavefields with the same phases, amplitudes and radiation patterns as the input elastic wavefields (Fig. 3). The relative errors between the input and reconstructed wavefields are 3.6 × 10−6 and 3.5 × 10−6 for horizontal and vertical components, respectively. The reconstructed wavefields are calculated by summing the separated P and Swavefields along the horizontal and vertical directions. These small errors indicate that the elastic wavefield decomposition using eq. (4) has a good accuracy. Figure 1. View largeDownload slide Multicomponent wavefields in a homogeneous medium with vp = 3 km s−1 and vs = 1.5 km s−1. Panels (a) and (b) are the horizontal and vertical components, respectively. An double-couple source is located in the middle of this model. Figure 1. View largeDownload slide Multicomponent wavefields in a homogeneous medium with vp = 3 km s−1 and vs = 1.5 km s−1. Panels (a) and (b) are the horizontal and vertical components, respectively. An double-couple source is located in the middle of this model. Figure 2. View largeDownload slide Wavefield separation experiment for the snapshots in Fig. 1 using the divergence and curl operators. Panel (a) and (b) are the P- and S-mode wavefields, respectively. Note that the divergence and curl operations change the phases, amplitudes and radiation patterns of the original elastic wavefields. Figure 2. View largeDownload slide Wavefield separation experiment for the snapshots in Fig. 1 using the divergence and curl operators. Panel (a) and (b) are the P- and S-mode wavefields, respectively. Note that the divergence and curl operations change the phases, amplitudes and radiation patterns of the original elastic wavefields. Figure 3. View largeDownload slide Wavefield separation experiment for the snapshots in Fig. 1 using the vector Helmholtz decomposition. Panels (a) and (b) are the horizontal and vertical components of separated P-mode wavefields. Panels (c) and (d) are the horizontal and vertical components of separated S-mode wavefields. Panel (e) and (f) are the horizontal and vertical differences between the original and reconstructed wavefields. Note that the separated wavefields have the same phase, amplitude and radiation pattern as the input elastic wavefields. Figure 3. View largeDownload slide Wavefield separation experiment for the snapshots in Fig. 1 using the vector Helmholtz decomposition. Panels (a) and (b) are the horizontal and vertical components of separated P-mode wavefields. Panels (c) and (d) are the horizontal and vertical components of separated S-mode wavefields. Panel (e) and (f) are the horizontal and vertical differences between the original and reconstructed wavefields. Note that the separated wavefields have the same phase, amplitude and radiation pattern as the input elastic wavefields. 2.2 Vector-based elastic time-reversal source imaging One difference between the passive-source survey and the reflection survey is that the former acquisition array is commonly irregular and sparsely distributed. If only using one wave-mode (e.g. P waves in acoustic RTM), the imperfect acquisition can introduce strong swinging (aliasing) artefacts in passive source images (Douma & Snieder 2015; Nakata & Beroza 2016). If a given passive-source generates both compressional and shear waves, elastic time-reversal back-propagation makes P and Swaves concurrently focusing at the true source location, and defocussing at other locations. Therefore, cross-correlating the separated P and Swavefields helps us to reduce the swinging artifacts, and produce a focused image (Artman et al. 2010). In this section, we introduce the vector-based elastic RTM (Wang & McMechan 2015; Wang et al. 2016; Zhu 2017) into the passive source imaging, and propose several optimal imaging conditions according to the characteristics of transmission waves. Similar as the reflection RTM workflow, the first step of elastic time-reversal source imaging is to compute the receiver wavefield using multicomponent records as the boundary condition. Taking the first-order velocity–stress wave equation as an example, the back-propagation can be expressed as \begin{eqnarray*} \begin{split} \rho (\mathbf {x})\frac{\partial {\mathbf {v}(\mathbf {x},t)}}{\partial {t}} = \mathbf {L}\boldsymbol {\sigma }(\mathbf {x},t), \quad \frac{\partial \boldsymbol {\sigma }(\mathbf {x},t)}{\partial {t}} = \mathbf {C}\mathbf {L}^{T}\mathbf {v}(\mathbf {x},t), \quad \text{with} \quad \mathbf {v}(\mathbf {x}_r,t) = \mathbf {d}(\mathbf {x}_r,t_{\rm {max}}-t), \end{split} \end{eqnarray*} (5) where |$\mathbf {x}=(x,y,z)$| is the location of an image point. |$\mathbf {x}_r$| is the receiver location. t is the time variable, and tmax is the record duration. Superscript T denotes the transpose. |$\mathbf {v}=[v_x,v_y,v_z]^{T}$| is the particle velocity wavefields. |$\boldsymbol {\sigma }=[\sigma _{xx},\sigma _{yy},\sigma _{zz},\sigma _{xy},\sigma _{xz},\sigma _{yz}]^{T}$| is the stress wavefields. |$\mathbf {d}=[{d}_x,{d}_y,{d}_z]^{T}$| is the observed multicomponent records. |$\mathbf {L}$| is a spatial partial-derivative matrix, and has the following expression: \begin{eqnarray*} \mathbf {L} = {\left[\begin{array}{*6c}\displaystyle\frac{\partial }{\partial {x}} & \quad 0 & \quad 0 & \quad \displaystyle\frac{\partial }{\partial {y}} & \quad \displaystyle\frac{\partial }{\partial {z}} & \quad 0 \\ \\ 0 & \quad \displaystyle\frac{\partial }{\partial {y}} & \quad 0 & \quad \displaystyle\frac{\partial }{\partial {x}} & \quad 0 & \quad \displaystyle\frac{\partial }{\partial {z}} \\ \\ 0 \quad & \quad 0 & \quad \displaystyle\frac{\partial }{\partial {z}} & \quad 0 & \quad \displaystyle\frac{\partial }{\partial {x}} & \displaystyle \frac{\partial }{\partial {y}} \end{array}\right]}. \end{eqnarray*} (6)|$\mathbf {C}$| is the model parameter matrix. For isotropic elastic media, it takes the form \begin{eqnarray*} \mathbf {C} = {\left[\begin{array}{*6c}\rho (\mathbf {x})v_p^2(\mathbf {x})\quad & \quad \rho (\mathbf {x})[v_p^2(\mathbf {x})-2v_s^2(\mathbf {x})] & \quad \rho (\mathbf {x})[v_p^2(\mathbf {x})-2v_s^2(\mathbf {x})] & \quad \rho (\mathbf {x})v_s^2(\mathbf {x}) & \quad 0 & \quad 0 \\ \\ \rho (\mathbf {x})[v_p^2(\mathbf {x})-2v_s^2(\mathbf {x})] & \quad \rho (\mathbf {x})v_p^2(\mathbf {x}) & \quad \rho (\mathbf {x})[v_p^2(\mathbf {x})-2v_s^2(\mathbf {x})]& \quad 0 & \quad \rho (\mathbf {x})v_s^2(\mathbf {x}) & \quad 0 \\ \\ \rho (\mathbf {x})[v_p^2(\mathbf {x})-2v_s^2(\mathbf {x})] & \quad \rho (\mathbf {x})[v_p^2(\mathbf {x})-2v_s^2(\mathbf {x})] & \quad \rho (\mathbf {x})v_p^2(\mathbf {x}) & \quad 0 & \quad 0 & \quad \rho (\mathbf {x})v_s^2(\mathbf {x}) \end{array}\right]}, \end{eqnarray*} (7) where |$v_p(\mathbf {x})$| is the P-wave velocity, |$v_s(\mathbf {x})$| is the S-wave velocity, and |$\rho (\mathbf {x})$| is the density. eq. (5) can be numerically solved using finite-difference (Virieux 1986; Graves 1996), finite-element (Marfurt 1984; Komatitsch et al. 2010), or spectra-element (Komatitsch & Tromp 2002) methods. The second step is to separate the back-propagated wavefields into P and Swaves using the vector Helmholtz decomposition introduced in the previous section, that is, \begin{eqnarray*} \begin{split} \mathbf {v}^{P}(\mathbf {x},t) = \nabla \left(\nabla \cdot \mathbf {w}(\mathbf {x},t)\right), \quad \mathbf {v}^{S}(\mathbf {x},t) = -\nabla \times \nabla \times \mathbf {w}(\mathbf {x},t), \quad \end{split} \end{eqnarray*} (8) where |$\mathbf {w}(\mathbf {x},t)$| satisfies \begin{eqnarray*} \nabla ^2\mathbf {w}(\mathbf {x},t) = \mathbf {v}(\mathbf {x},t). \end{eqnarray*} (9) The vector Poisson’s equation can be solved using numerical algorithms, such as Fourier transform, LU factorization and multigrid methods (Yang et al. 2017; Zhu 2017). For 2-D problems, the methods based on Fourier transform and LU factorization have good accuracy and acceptable computational costs. For 3-D problems, the multigrid method is the most applicable approach in terms of both computational precision and efficiency. Without a strict requirement for the imaging accuracy, the phase- and amplitude-correction strategy proposed by Yang et al. (2018b) can be used, which gives us similar computational costs as Yan & Sava (2008)’s method. The final step is to extract the source information from the 4-D P- and S-wavefield volumes. Since the P and S wavefields are all vectorial, we use a dot-product imaging condition (Wang et al. 2016; Zhu 2017) to collapse the time axis, and produce a PS image for subsurface passive sources. This can be expressed as \begin{eqnarray*} I_{1}(\mathbf {x}) = \int _{0}^{t_{\rm {max}}}\mathbf {v}^{P}(\mathbf {x},t)\cdot \mathbf {v}^{S}(\mathbf {x},t){\rm d}t, \end{eqnarray*} (10) where · denotes the dot-product, |$\mathbf {v}^{P}(\mathbf {x},t)$| and |$\mathbf {v}^{S}(\mathbf {x},t)$| are the separated P and Swavefields. As analysed by Artman et al. (2010), correlating the P- and S-wave energy density functions helps to suppress weak scattering artifacts and produce a high SNR passive source image. Here, we propose two additional imaging conditions for the separated vector wavefields: \begin{eqnarray*} \begin{split} I_{2}(\mathbf {x}) &= \int _{0}^{t_{\rm {max}}}\left|\mathbf {v}^{P}(\mathbf {x},t)\right|^2\left|\mathbf {v}^{S}(\mathbf {x},t)\right|^2{\rm d}t, \\ \end{split} \end{eqnarray*} (11) and \begin{eqnarray*} \begin{split} I_{3}(\mathbf {x}) &= \int _{0}^{t_{\rm {max}}}\left|\mathbf {v}^{P}(\mathbf {x},t)\cdot \mathbf {v}^{S}(\mathbf {x},t)\right|^{p}{\rm d}t, \end{split} \end{eqnarray*} (12) where |$|\mathbf {v}^{P}(\mathbf {x},t)|^2$| is the energy of the Pwavefield, which can be computed as |$|\mathbf {v}^{P}(\mathbf {x},t)|^2={\sum \limits _{x_i}[{v}^{P}_{x_i}(\mathbf {x},t)]^2}$|⁠, xi = x, y, z, |$|\mathbf {v}^{S}(\mathbf {x},t)|^2$| has the same definition for the S wavefield. Superscript p is the power parameter, and we choose p = 2 for the following numerical examples. Note that, |$I_2(\mathbf {x})$| is computed by summing the multiplication of P- and S-wave energies, and we refer to eq. (11) as the energy imaging condition. Similarly, we refer to eq. (12) as the power imaging condition. Figs 4 and 5 illustrate an experiment for imaging a double-couple point source in a homogeneous medium. P- and S-wave velocities are 3.0  and 1.5  km s−1, respectively. The double-couple source is located at the depth of 2.5  km and the distance of 2  km. 401 receivers are evenly deployed at the surface with a spacing of 10 m. For comparison, we also compute the imaging results (Figs 4a–c) using the divergence and curl operators (Yan & Sava 2008; Artman et al. 2010). As noted in Artman et al. (2010), PS image has a better resolution than PP and SS images, and is robust for different types of sources. However, the divergence and curl operations change the phases and amplitudes of the decomposed wavefields (Fig. 2), and their resulting image is asymmetric (Fig. 4a). This result might be erroneously interpreted as the occurrence of four small events. By cross-correlating P- and S-wave energies or taking the squares of P- and S-wavefield multiplication, the scalar energy and power imaging conditions (Figs 4b and c) enable us to suppress certain imaging artifacts, and make the point-source image more focused. But they are still asymmetric. In contrast, the vector-based dot-product image condition produces a focused image with one dominant peak at the correct source location (Figs 4d, and 5a, d). Corresponding energy and power imaging conditions achieve a better spatial resolution for the event location (Figs 4e, f and 5b, c, e, f). Figure 4. View largeDownload slide Passive source imaging experiments for a double-couple point source in a homogeneous medium with vp = 3  km s−1 and vs = 1.5 km. The model size is 4 km × 4  km, and only the area around the source is plotted. The point source is located at the depth of 2.5  km and the distance of 2  km. The top row denoted by 'Scalar' is the imaging results with the separated wavefields using the divergence and curl operators, and the bottom row denoted by 'Vector' is the results using the vector Helmholtz decomposition. The three columns correspond to waveform, energy and power imaging conditions, respectively. Figure 4. View largeDownload slide Passive source imaging experiments for a double-couple point source in a homogeneous medium with vp = 3  km s−1 and vs = 1.5 km. The model size is 4 km × 4  km, and only the area around the source is plotted. The point source is located at the depth of 2.5  km and the distance of 2  km. The top row denoted by 'Scalar' is the imaging results with the separated wavefields using the divergence and curl operators, and the bottom row denoted by 'Vector' is the results using the vector Helmholtz decomposition. The three columns correspond to waveform, energy and power imaging conditions, respectively. Figure 5. View largeDownload slide Detailed views of the passive source images in Fig. 4. Panels (a)–(c) are the horizontal results at the depth of 2.5  km. Panels (d)–(f) are the vertical results at the distance of 2.0  km. Black solid lines are the true source location, blue dashed lines are the results from the scalar imaging condition and red dashed lines are the results from the vector imaging condition. Figure 5. View largeDownload slide Detailed views of the passive source images in Fig. 4. Panels (a)–(c) are the horizontal results at the depth of 2.5  km. Panels (d)–(f) are the vertical results at the distance of 2.0  km. Black solid lines are the true source location, blue dashed lines are the results from the scalar imaging condition and red dashed lines are the results from the vector imaging condition. The asymmetry in the PS image from Artman et al. (2010)’s method can be explained using the polarity reversal phenomena that occurs in reflection elastic RTM (Yan & Sava 2008). Since the curl operator produces opposite polarity directions for Swavefields that propagate along opposite directions, the migrated PS-reflector polarity is reversed (Du et al. 2012). By separating elastic wavefields using the Helmholtz decomposition, and applying the dot-product imaging condition, the vector-based elastic RTM avoids the polarity reversal problem, and can produce constructively stacked PS images without particular polarity correction (Zhu 2017). This advantage also helps the vector-based time-reversal source imaging to produce a symmetric image for a single point source (Figs 4d–f). Eqs (10)–(12) collapse the time axis of the backward propagated wavefields by summing the entire record duration. When there are multiple events excited at different times, this method produces an average imaging result. But for some natural phenomenon, such as hydraulic fractures and earthquake ruptures, monitoring their dynamic propagation is important for guiding on-site fracturing operation and studying fault failure mechanism. To capturing these dynamic procedures, we modify the dot-product imaging condition in eq. (10) as \begin{eqnarray*} I_{1}(\mathbf {x},\tau ) = \int _{\tau _1(\tau )}^{\tau _2(\tau )}\mathbf {v}^{P}(\mathbf {x},t)\cdot \mathbf {v}^{S}(\mathbf {x},t){\rm d}t, \end{eqnarray*} (13) where τ is the time variable for the dynamic imaging results, and [τ1(τ), τ2(τ)] is local time window. Similarly, the dynamic energy and power imaging conditions can be formulated as \begin{eqnarray*} \begin{split} I_{2}(\mathbf {x},\tau ) &= \int _{\tau _1(\tau )}^{\tau _2(\tau )}\left|\mathbf {v}^{P}(\mathbf {x},t)\right|^2\left|\mathbf {v}^{S}(\mathbf {x},t)\right|^2{\rm d}t, \\ I_{3}(\mathbf {x},\tau ) &= \int _{\tau _1(\tau )}^{\tau _2(\tau )}\left|\mathbf {v}^{P}(\mathbf {x},t)\cdot \mathbf {v}^{S}(\mathbf {x},t)\right|^{p}{\rm d}t. \end{split} \end{eqnarray*} (14) By setting different values for τ1(τ) and τ2(τ), eqs (13) and (14) can produce different dynamic imaging results. For example, setting τ1(τ) = 0 and τ2(τ) = τ produces a series of time-lapse migration results. This setting means that as receivers continuously record surface motions, the recorded seismograms with different end times τ are used for time-reversal source imaging. This approach requires multiple wavefield extrapolations and is applicable for real-time monitoring problems, such as hydraulic fracturing monitoring. Another choice for the cross-correlation window is τ1(τ) = τ − Δτ and τ2(τ) = τ + Δτ. Δτ is a short time window, which can be set as the half period of recorded seismograms. This approach can be implemented using one wavefield back-propagation, and each resulting image illustrates subsurface source activities at a given excitation-time isochron τ. This imaging condition can be used to locate earthquake hypocentres, and monitor corresponding rupture propagation for moderate- and large-sized earthquakes. 3 NUMERICAL EXAMPLES In this section, three synthetic examples are used to illustrate the performance of the proposed passive-source locating and monitoring method. In the first example, we use the Marmousi model to test the feasibility of the proposed method for locating a point source, and its sensitivity to noisy and irregularly acquired data. In the second example, we design static and dynamic source signatures for the 2-D SEAM model to mimic a realistic hydraulic fracturing procedure. Then, we use corresponding seismic data to test the performance of the proposed method for monitoring crack propagation. In the final example, we first compute a 3-D data set using a kinematic finite-fault rupture model as the earthquake source, and use this dataset to test the feasibility of the proposed method for imaging earthquake rupture propagation. 3.1 A point source The P- and S-wave velocities of the Marmousi model are shown in Fig. 6. Corresponding migration velocities are built by smoothing the true models using a 180 m × 180 m Gaussian filter. We discretize this model using a grid of 251 × 767 with a spatial spacing of 12 m. A vertical single-force is deployed in the interior of the bottom anticline reservoir (black dot in Fig. 6 a at the depth of 2.01 km and horizontal position of 5.68 km). In this example, a surface survey (black triangles in Fig. 6a) and a well survey (red squares in Fig. 6a) are used. Synthetic passive-source seismograms are computed by solving the elastic wave equation using the staggered-grid finite-difference scheme with eighth-order accuracy in space and second-order accuracy in time. Figure 6. View largeDownload slide The Marmousi model. (a) is the P-wave velocity. (b) is the S-wave velocity. The black dot denotes the passive point source. Black triangles denote the receivers in the surface survey. Red squares denote the receivers in the well survey. Figure 6. View largeDownload slide The Marmousi model. (a) is the P-wave velocity. (b) is the S-wave velocity. The black dot denotes the passive point source. Black triangles denote the receivers in the surface survey. Red squares denote the receivers in the well survey. For the surface survey, 383 receivers are evenly distributed at the surface. The time sampling space is 1 ms and the record duration is 6 s. Fig. 7 compares the scalar- and vector-based passive-source imaging results. Since we have a good receiver coverage at the surface, both imaging methods successfully resolve the point source. Due to the wave-mode conversions at stratum interfaces during elastic wave propagation, Figs 7(a) and (b) not only locate the point source but also image certain geological structures. This property can be used to image Earth’s discontinuities with large velocity contrasts, such as the Moho and subduction zone boundaries (Shang et al. 2012, 2017). The energy and power imaging conditions allow us to suppress these weak structure features, and produce higher resolution for the point source location (Figs 7c–f). As analysed in the previous section, because the divergence and curl operators change the phases and amplitudes of the extrapolated wavefields, the scalar imaging condition produces asymmetric results (Figs 7a, c and e). Their two lobes might lead to the interpretation of two microseismic events. Even worse, the true point source location (5.68  km, 2.01  km) falls into the gaps between these two lobes. This means that although the PS image obtained using Artman et al. (2010)’s method can be used to find passive sources, their resulting locations are inaccurate. In contrast, the vector Helmholtz decomposition produces the same phases and amplitudes as the coupled elastic wavefields, which enables the dot-product, energy and power imaging conditions to produce symmetric imaging results (Figs 7b, d and f). Their dominant peaks correspond to the true point source location. Figure 7. View largeDownload slide Passive-source imaging results of the Marmousi model for the surface survey. The first column shows the scalar-based elastic RTM results. The second column shows the vector-based elastic RTM results. Green stars denote the source location. The insets at the bottom left corners are the zoomed views for the area around the source location. Figure 7. View largeDownload slide Passive-source imaging results of the Marmousi model for the surface survey. The first column shows the scalar-based elastic RTM results. The second column shows the vector-based elastic RTM results. Green stars denote the source location. The insets at the bottom left corners are the zoomed views for the area around the source location. For the well survey, 125 receivers are evenly deployed along the well path (red squares in Fig. 6a). The time sample interval is 1 ms and the record duration is 4 s. Corresponding scalar and vector elastic RTM results are shown in Fig. 8. Since the illumination of the incident angles is insufficient, there are some artifacts to the left of the point source image. Similar to the surface survey, scalar-based RTMs produce asymmetric results with multiple lobes, and neither lobes pinpoint the true point source location (Figs 8a–c). On the contrary, vector-based elastic RTMs produce symmetric images with a dominant peak at the right source position (Figs 8d–f). The dot-product imaging condition produces a few structural reflectors (Fig. 8d). The energy imaging condition reduces these structural features but yields a smearing source image. The result derived from the power imaging condition has the highest spatial resolution, and can be easily interpreted as one microseismic event (Fig. 8f). Figure 8. View largeDownload slide Passive-source imaging results of the Marmousi model with the well survey. The first row shows the scalar-based elastic RTM results. The second row shows the vector-based elastic RTM results. Green stars denote the source location. The insets at the upper right corners are the zoomed views for the area around the source location. Figure 8. View largeDownload slide Passive-source imaging results of the Marmousi model with the well survey. The first row shows the scalar-based elastic RTM results. The second row shows the vector-based elastic RTM results. Green stars denote the source location. The insets at the upper right corners are the zoomed views for the area around the source location. To test the sensitivity of the proposed elastic RTM approach for the presence of random noises and sparsely acquired data, we generate two additional datasets. One is built by adding Gaussian white noises to the surface survey gathers (Figs 9a and b), whose SNR is 1. The other is obtained by randomly decimating |$90 \, {per \,cent}$| of the noisy gathers (Figs 9c and d). Their corresponding passive-source imaging results are displayed in Fig. 10. It is observable that Gaussian noises introduce strong artifacts at the near surface. This is because that strong converted waves, such as P-to-S or S-to-P, are generated when inserting the noisy multicomponent records into the boundaries of the extrapolated elastic wavefields. Similar wave-mode conversions at the stratum interfaces also produce deep structural information (Figs 10a, c and e). However, the point source is still imaged clearly with stronger energies than migration noises. With irregularly decimated data, aliasing artifacts occur at shallow areas, especially in the energy and power imaging condition results (red smearing at depths shallower than 1.5  km in Figs 10b, d, and f). Coherent stacking from different traces enhances the imaging SNR for the point source, and produces a clear energy peak (Figs 10 d and f). Figure 9. View largeDownload slide Noisy and decimated common-shot gathers of the Marmousi model with the surface survey. Panels (a) and (b) are the horizontal and vertical components of the noisy data. Panels (c) and (d) are the horizontal and vertical components decimated noisy data. The signal-to-noise ratio is 1 for all panels and |$10 \, {per \,cent}$| data is left in (c) and (d). Note that no continuous seismic events can be visually picked in panels (c) and (d). Figure 9. View largeDownload slide Noisy and decimated common-shot gathers of the Marmousi model with the surface survey. Panels (a) and (b) are the horizontal and vertical components of the noisy data. Panels (c) and (d) are the horizontal and vertical components decimated noisy data. The signal-to-noise ratio is 1 for all panels and |$10 \, {per \,cent}$| data is left in (c) and (d). Note that no continuous seismic events can be visually picked in panels (c) and (d). Figure 10. View largeDownload slide Passive-source imaging results of the Marmousi model for the surface noisy and decimated data. The first column shows results for the noisy data (Figs 9 a and b). The second column shows results for the noisy decimated data (Figs 9 c and d). Green stars denote the source location. The insets at the bottom left corners are the zoomed views for the area around the source location. Figure 10. View largeDownload slide Passive-source imaging results of the Marmousi model for the surface noisy and decimated data. The first column shows results for the noisy data (Figs 9 a and b). The second column shows results for the noisy decimated data (Figs 9 c and d). Green stars denote the source location. The insets at the bottom left corners are the zoomed views for the area around the source location. 3.2 Hydraulic crack propagation The simple portion of the 2-D SEAM model is used in this example. Its P- , S-wave velocities and density models are shown in Figs 11(a)–(c). Migration velocities (Fig. 11d) are built by applying a 200 m × 200 m Gaussian filter to the true velocity models. Density is set as a constant of 2.0 g cm–3 in migration. This model is discretized with a grid of 600 × 700 and a spacing of 20 m. 350 receivers are evenly deployed at the surface to record passive-source seismograms. The time sampling interval is 1 ms and the record duration is 10 s. We design two passive source signatures to simulate hydraulic fracture: one is static, and the other is dynamic. The static source is shown in Fig. 12(a), in which case, all microseismic sources along the crack path are assumed to be simultaneously excited. The source time function is a Ricker wavelet with the peak frequency of 10 Hz (Fig. 12b). The dynamic source is shown in Fig. 14, in which case, hydraulic fractures propagate from right to left with a fracture velocity of 5  km s−1. Figure 11. View largeDownload slide The 2-D SEAM model: panel (a) is the P-wave velocity; panel (b) is the S-wave velocity; panel (c) is the density; panel (d) is the migration velocity, which is built by smoothing the true model using a Gaussian filter of 200 m × 200 m. Triangles in panel (a) denote the receivers, and the dendritic curves denote the fracturing propagation path. Figure 11. View largeDownload slide The 2-D SEAM model: panel (a) is the P-wave velocity; panel (b) is the S-wave velocity; panel (c) is the density; panel (d) is the migration velocity, which is built by smoothing the true model using a Gaussian filter of 200 m × 200 m. Triangles in panel (a) denote the receivers, and the dendritic curves denote the fracturing propagation path. Figure 12. View largeDownload slide A static crack signature (a) and the corresponding source time function for the 2-D SEAM model. Figure 12. View largeDownload slide A static crack signature (a) and the corresponding source time function for the 2-D SEAM model. Elastic time-reversal imaging results for the static fracture sources are shown in Fig. 13. Since there are many stratum interfaces with large velocity and density contrasts in the SEAM model, the dot-product imaging condition not only produces clear hydraulic fracturing paths, but also images certain structures above the passive sources (Figs 13a and d). The energy imaging condition removes most of these structural features, but its resulting image has a degraded resolution and just captures several dominant fracture segments (Figs 13b and d). The power imaging condition produces the best passive-source image (Figs 13c and d). Most structural features are suppressed and the spatial resolution is enhanced in comparison to the dot-product imaging condition result. Figure 13. View largeDownload slide Static passive-source imaging results for the 2-D SEAM model: panel (a) is the result using the dot-product imaging condition; panel (b) is the result using the energy imaging condition; panel (c) is the result using the power imaging condition; panel (d) shows the comparisons of these three results with the true source signature. Figure 13. View largeDownload slide Static passive-source imaging results for the 2-D SEAM model: panel (a) is the result using the dot-product imaging condition; panel (b) is the result using the energy imaging condition; panel (c) is the result using the power imaging condition; panel (d) shows the comparisons of these three results with the true source signature. Figure 14. View largeDownload slide A dynamic fracture signature for the 2-D SEAM model. Panels (a), (b) and (c) are the source distribution at 0.3, 0.4 and 0.6 s. Panel (d) is the source time function at four different locations. The fracture is excited at the right part and propagates toward left. Figure 14. View largeDownload slide A dynamic fracture signature for the 2-D SEAM model. Panels (a), (b) and (c) are the source distribution at 0.3, 0.4 and 0.6 s. Panel (d) is the source time function at four different locations. The fracture is excited at the right part and propagates toward left. In the dynamic source case, multicomponent records within four time windows of [0, 2.65 s], [0, 2.9 s], [0, 3.4 s] and [0, 4.9 s] are used during the back-propagation, and their imaging results are shown in Fig. 15. Initial fracturing are accurately resolved in the elastic time-reversal source image using early recorded seismic data (Figs 15a and b). With the increasing recording duration, later fractures are gradually imaged (Figs 15c–f). After all microseismic events are recorded, the complete fracturing paths are illuminated (Figs 15g and h). As the previous static case, the image computed using the dot-product imaging condition contains a few structural reflectors, while the power imaging condition produces a clear and high-resolution image for the fracturing paths. This test suggests that the proposed elastic RTM can be used to locate subsurface microseismic sources, as well as monitoring their dynamic propagation. Figure 15. View largeDownload slide Dynamic passive-source imaging results for the 2-D SEAM model. The first column shows the results using the dot-product imaging condition. The second column shows the results using the power imaging condition. Four rows correspond to four different back-propagation time windows. The time marked on each panel is τ2(τ) in eq. (14), and τ1(τ) is set to 0. Figure 15. View largeDownload slide Dynamic passive-source imaging results for the 2-D SEAM model. The first column shows the results using the dot-product imaging condition. The second column shows the results using the power imaging condition. Four rows correspond to four different back-propagation time windows. The time marked on each panel is τ2(τ) in eq. (14), and τ1(τ) is set to 0. 3.3 Earthquake kinematic rupture imaging In this example, we apply the proposed elastic time-reversal imaging to monitor earthquake rupture propagation. The 3-D SEG/EAGE overthrust model (Fig. 16) is used as the background velocity in the finite-fault modeling. S-wave velocity model is built by scaling P-wave velocity with a factor of |$1/\sqrt{3}$|⁠. Density is a constant of 2.0 g cm–3. We discretize this model with a grid of 187 × 401 × 401 and a spatial spacing of 25 m. A vertical fault plane located at the cross line position of 5  km is represented by 11 × 28 cells, and each cell is considered as a subfault. The CMT parameters (Figs 17a) for these subfaults are obtained from the finite-fault inversion result for 2007 Mw 8.39 Sumatra Earthquake Ji et al. (2002). The rupture starts from the centre of the fault plane with a large slip, and propagates along all directions with a rupture speed of 1.2  km s−1 (yellow isochron in Fig. 17). Multicomponent seismograms are computed using a 3-D staggered-grid finite-difference scheme. 200 × 200 stations are regularly deployed at the Earth’s surface. Figure 16. View largeDownload slide P-wave velocity for the 3-D SEG/EAGE Overthrust model. S-wave velocity model is built by applying a scale factor of |$1/\sqrt{3}$| to P-wave velocity. Density is a constant of 2.0 g cm–3. Figure 16. View largeDownload slide P-wave velocity for the 3-D SEG/EAGE Overthrust model. S-wave velocity model is built by applying a scale factor of |$1/\sqrt{3}$| to P-wave velocity. Density is a constant of 2.0 g cm–3. Figure 17. View largeDownload slide A finite-fault slip distributions (a) and its time-reversal imaging results (b) for a vertical fault plane in the Overthrust model. Dip is a constant of 12°. The contour curves denote the initial time isochron of the rupture propagation. Detailed CMT parameters can be found at http://www.geol.ucsb.edu/faculty/ji/big_earthquakes/2007/09/result_sumatra/static_out. Figure 17. View largeDownload slide A finite-fault slip distributions (a) and its time-reversal imaging results (b) for a vertical fault plane in the Overthrust model. Dip is a constant of 12°. The contour curves denote the initial time isochron of the rupture propagation. Detailed CMT parameters can be found at http://www.geol.ucsb.edu/faculty/ji/big_earthquakes/2007/09/result_sumatra/static_out. We set τ1(τ) = τ − 0.2 s and τ2(τ) + 0.2 s for the imaging conditions in eqs (13) and (14), and obtain a series of dynamic images (Figs 18 and 19). From imaging profiles, we see that the rupture is activated at 0 s, and at the horizontal location of 5  km and the vertical location of 2.7 km (top right panel in Fig. 18). In addition, its amplitude is larger in comparison with later imaging results, indicating the initial rupture has large slip displacements (top left-hand panel in Fig. 18). As time increasing, the rupture gradually propagates in different directions. Small energy peaks, for example, at 0.8 and 1.2 s, indicate the initial dominant fault beaks into several subfaults during rupture propagation. The locating times for these subfaults are consistent with the isochron used in the finite-fault modelling (left-hand column in Figs 18 and 19). It is notable that these two small events at the upper left (2.2  km vertically and 3.1  km horizontally around 1.5 s isochron in Fig. 17a) and the bottom right (3.6  km vertically and 7.8  km horizontally around 2.5 s isochron) corners can also be accurately resolved in the panels at 1.6 and 2.4 s (right-hand column in Fig. 19). These results indicate that the proposed elastic time-reversal imaging method is capable of accurately locating and monitoring the kinematic rupture propagation. Fig. 17(b) shows the cumulative depth profile of the fault plane computed using eq. (12). Both the dominant rupture in the middle of the fault plane, and the secondary ruptures in the upper left and bottom right corners are imaged. Their spatial patterns are consistent with the input finite fault slip distribution (Fig. 17a). But due to unbalanced energy distribution at different depths, it is difficult for our method to locate small slips at great depth, such as the the event at depth of 4  km around 1.2 s isochron (Fig. 18), and the event at the depth of 3.5  km and the distance of 7  km around 2.0 s isochron. This is because the proposed time-reversal imaging method is an adjoint instead of the inverse operator. To produce more accurate imaging results, especially for the deep small slips, the inversion-based methods should be used (Melgar et al. 2015; Zhang et al. 2015). Figure 18. View largeDownload slide Comparisons of slip distribution (left-hand column) and dynamic passive-source imaging results (right-hand column) for the Overthrust model from 0 to 1.2 s. The power imaging condition is used in this experiment. The time marked in each panel is the rupture propagation time. Figure 18. View largeDownload slide Comparisons of slip distribution (left-hand column) and dynamic passive-source imaging results (right-hand column) for the Overthrust model from 0 to 1.2 s. The power imaging condition is used in this experiment. The time marked in each panel is the rupture propagation time. Figure 19. View largeDownload slide Comparisons of slip distribution (left-hand column) and dynamic passive-source imaging results (right-hand column) for the Overthrust model from 01.6 to 2.8 s. The power imaging condition is used in this experiment. The time marked in each panel is the rupture propagation time. Figure 19. View largeDownload slide Comparisons of slip distribution (left-hand column) and dynamic passive-source imaging results (right-hand column) for the Overthrust model from 01.6 to 2.8 s. The power imaging condition is used in this experiment. The time marked in each panel is the rupture propagation time. 4 DISCUSSION A vector-based elastic time-reversal seismic-source imaging method is developed for transmission multicomponent seismic data. In this procedure, subsurface elastic wavefields are first constructed by solving the elastic wave equation using multicomponent seismograms as boundary conditions. Then, the vector Helmholtz decomposition is used to separate pure-mode compressional (P) and shear (S) wavefields, which have the same amplitudes, phases and physical units as the input coupled elastic wavefields. Finally, the zero-lag cross-correlation imaging condition is applied to the separated P and Swavefields for the entire record duration to produce a complete passive-source image, or within local time windows to produce a series of dynamic images. In real applications, the source time function is unknown for hydraulic fracture and earthquake rupture. The imaging time-window intervals can be set to the average period of the peak frequency of observed seismic events. The stacking within local time windows helps to suppress imaging artifacts and improve imaging quality. Three calculations, that is, dot-product, energy and power, are used to implement the cross-correlation image condition. Numerical examples reveal that the proposed transmission elastic RTM method can be used to locate point sources, as well as monitoring dynamic hydraulic fracture and earthquake rupture propagation. In the Marmousi example, we analyse the sensitivity of the proposed method to the presence of noisy and sparsely acquired data. Strong noises in the records introduce many migration artifacts in the passive-source image, especially at near surface. Irregularly acquired seismograms also cause strong aliasing effects at shallow depths. If the passive sources are located at deep areas, coherent stacking helps to reduce aliasing artefacts and produce clear source images. If the passive sources are located at near surface, high SNR and regularly sampled data are important for producing high-quality source images. However, field data commonly has low SNR and is sparsely acquired at surface. In these cases, proper denoising and regularization techniques should be used in the pre-processing step. For instance, the sparsity-promoted denoising and interpolation approaches, such as in the Fourier domain Xu et al. (2010) and curvelet domain (Shang et al. 2017), can be used to enhance the effective transmission events. For global earthquake data, using multiple networks and seismic phases can help to improve vertical and horizontal image resolutions (Kiser & Ishii 2017). In all RTM examples, we use smoothed models as background velocities. When migration velocities are inaccurate, the P and S wavefields will not be concurrently focused at the right source locations. Fig. 20 illustrate the point-source image for the Marmousi surface survey using inaccurate velocity. It is notable that inaccurate velocity models give rise to erroneous source location and degraded resolution. Accurate background velocity models are important to produce high-resolution, focused passive-source images. At a reservoir scale, the background velocity models can be computed using ray-based tomography (Dapeng 2008; Woodward et al. 2008), or wave-equation-based waveform inversion (Tromp et al. 2005; Virieux & Operto 2009). Since hydraulic microseismic sources are commonly distributed sparsely and irregularly, it is difficult to directly use microseismic data to simultaneously invert for velocity models and source signatures. One way to mitigate this problem is to use both active- and passive-source data to constrain the background velocities. This requires further investigations. At crustal and global scales, 1-D reference velocity model, such as IASP91 (Kennett & Engdahl 1991) or more accurate 3-D tomographic velocity models, can be used as the background velocities for wavefield extrapolation and passive-source imaging. In addition, since the proposed method enables us to delineate dynamic crack and rupture propagation, it can also be used for monitoring volcano (Sens-Schönfelder & Wegler 2006) and glacier activities (Köhler et al. 2015). Figure 20. View largeDownload slide Elastic time-reverse imaging results of Marmousi model using inaccurate migration velocity models. The migration velocity models for the results in panels (a) and (b) are constructed by scaling the smoothed model with factors of 0.95 and 0.9, respectively. The power imaging condition is used in this experiment. Figure 20. View largeDownload slide Elastic time-reverse imaging results of Marmousi model using inaccurate migration velocity models. The migration velocity models for the results in panels (a) and (b) are constructed by scaling the smoothed model with factors of 0.95 and 0.9, respectively. The power imaging condition is used in this experiment. 5 CONCLUSION In this paper, we introduce the vector Helmholtz decomposition into elastic time-reversal imaging and develop an efficient and accurate passive-source locating and monitoring method. Three implementations, that is, dot-product, energy and power, are used in the cross-correlation imaging condition for the separated vector P and Swavefields, in which the power imaging condition has the best imaging resolution. These new imaging conditions not only avoid the asymmetric problem of conventional divergence-curl-based time-reversal imaging method for a single point source, but also can locate the true source location with a dominant energy peak. By cross-correlating the separated wavefields within local time windows, the proposed method can be used to locate and monitor the dynamic behaviour of subsurface passive sources. Numerical examples demonstrate its feasibility and adaptability for simple point source as well as complex hydraulic fracture and earthquake rupture. ACKNOWLEDGEMENTS This research is supported by the Geosciences Department Ph.D. Fellowship. We are grateful to the computational resources provided by the Texas Advanced Computing Center. We appreciate comments from editor H. Yao and two anonymous reviewers that significantly improved the manuscript. This paper is contribution no. 1338 from the Department of Geosciences at the University of Texas at Dallas. REFERENCES Aki K. , Richards P. , 1980 . Quantitative Seismology: Theory and Methods , Vol. 1 . WH Freeman and Company . Artman B. , Podladtchikov I. , Witten B. , 2010 . Source location using timePLBIBITALIC–reverse imaging , Geophys. Prospect. , 58 ( 5 ), 861 – 873 . https://doi.org/10.1111/j.1365-2478.2010.00911.x Google Scholar Crossref Search ADS Baig A. , Urbancic T. , 2010 . Microseismic moment tensors:a path to understanding frac growth , Leading Edge , 29 ( 3 ), 320 – 324 . https://doi.org/10.1190/1.3353729 Google Scholar Crossref Search ADS Baysal E. , Kosloff D.D. , Sherwood J.W.C. , 1983 . Reverse time migration , Geophysics , 48 ( 11 ), 1514 – 1524 . https://doi.org/10.1190/1.1441434 Google Scholar Crossref Search ADS Bose S. , Valero H.P. , Liu Q. , Shenoy R.G. , Ounadjela A. , 2009 . An automatic procedure to detect microseismic events embedded in high noise , SEG Technical Program Expanded Abstracts , 1537 – 1541 . Claerbout J.F. , 1971 . Toward a unified theory of reflector mapping , Geophysics , 36 ( 3 ), 467 – 481 . https://doi.org/10.1190/1.1440185 Google Scholar Crossref Search ADS Dai F. , Li B. , Xu N. , Fan Y. , Zhang C. , 2016 . Deformation forecasting and stability analysis of large-scale underground powerhouse caverns from microseismic monitoring , Int. J. Rock Mech. Mining Sci. , 86 , 269 – 281 . https://doi.org/10.1016/j.ijrmms.2016.05.001 Google Scholar Crossref Search ADS Dapeng Z. , 2008 . New advances of seismic tomography and its applications to subduction zones and earthquake fault zones: a review , Island Arc , 10 ( 1 ), 68 – 84 . Douma J. , Snieder R. , 2015 . Focusing of elastic waves for microseismic imaging , Geophys. J. Int. , 200 ( 1 ), 390 – 401 . https://doi.org/10.1093/gji/ggu398 Google Scholar Crossref Search ADS Du Q. , Zhu Y. , Ba J. , 2012 . Polarity reversal correction for elastic reverse time migration , Geophysics , 77 ( 2 ), S31 – S41 . Google Scholar Crossref Search ADS Duan Y. , Sava P. , 2015 . Scalar imaging condition for elastic reverse time migration , Geophysics , 80 ( 4 ), S127 – S136 . https://doi.org/10.1190/geo2014-0453.1 Google Scholar Crossref Search ADS Graves R.W. , 1996 . Simulating seismic wave propagation in 3D elastic media using staggeredPLBIBITALIC–grid finite differences , Bull. seism. Soc. Am. , 86 ( 4 ), 1091 . Hayles K. , Horine R.L. , Checkles S. , Blangy J.P. , 2012 . Comparison of microseismic results from the bakken formation processed by three different companies: Integration with surface seismic and pumping data , SEG Technical Program Expanded Abstracts 2011 , 1468 – 1472 . Ishii M. , Shearer P.M. , Houston H. , Vidale J.E. , 2005 . Extent, duration and speed of the 2004 Sumatra-Andaman earthquake imaged by the Hi-Net array , Nature , 435 ( 7044 ), 933 . https://doi.org/10.1038/nature03675 Google Scholar Crossref Search ADS PubMed Ji C. , Wald D.J. , Helmberger D.V. , 2002 . Source description of the 1999 Hector Mine, California, Earthquake, part I: Wavelet domain inversion theory and resolution analysis , Bull. seism. Soc. Am. , 92 ( 4 ), 1192 . https://doi.org/10.1785/0120000916 Google Scholar Crossref Search ADS Jonas F. , Jrn K.A.S.S. , Markus H. , Hiroshi A. , 2016 . Rupture directivity of fluid-induced microseismic events: observations from an enhanced geothermal system , J. geophys. Res.: Solid Earth , 121 ( 11 ), 8034 – 8047 . https://doi.org/10.1002/2016JB013078 Google Scholar Crossref Search ADS Kao H. , Shan S.-J. , 2004 . The source-scanning algorithm: mapping the distribution of seismic sources in time and space , Geophys. J. Int. , 157 ( 2 ), 589 – 594 . https://doi.org/10.1111/j.1365-246X.2004.02276.x Google Scholar Crossref Search ADS Kennett B.L.N. , Engdahl E.R. , 1991 . Traveltimes for global earthquake location and phase identification , Geophys. J. Int. , 105 ( 2 ), 429 – 465 . https://doi.org/10.1111/j.1365-246X.1991.tb06724.x Google Scholar Crossref Search ADS Kiser E. , Ishii M. , 2017 . Back-projection imaging of earthquakes , Ann. Rev. Earth Planet. Sci. , 45 ( 1 ), 271 – 299 . https://doi.org/10.1146/annurev-earth-063016-015801 Google Scholar Crossref Search ADS Köhler A. , Nuth C. , Schweitzer J. , Weidle C. , Gibbons S.J. , 2015 . Regional passive seismic monitoring reveals dynamic glacier activity on Spitsbergen, Svalbard , Polar Res. , 34 ( 1 ), 26178 . https://doi.org/10.3402/polar.v34.26178 Google Scholar Crossref Search ADS Komatitsch D. , Tromp J. , 2002 . Spectral-element simulations of global seismic wave propagation-ii. three-dimensional models, oceans, rotation and self-gravitation , Geophys. J. Int. , 150 ( 1 ), 303 – 318 . https://doi.org/10.1046/j.1365-246X.2002.01716.x Google Scholar Crossref Search ADS Komatitsch D. , Erlebacher G. , Gddeke D. , Micha D. , 2010 . High-order finite-element seismic wave propagation modeling with mpi on a large gpu cluster , J. Comput. Phys. , 229 ( 20 ), 7692 – 7714 . https://doi.org/10.1016/j.jcp.2010.06.024 Google Scholar Crossref Search ADS Kummerow J. , 2010 . Using the value of the crosscorrelation coefficient to locate microseismic events , Geophysics , 75 ( 4 ), MA47 . https://doi.org/10.1190/1.3463713 Google Scholar Crossref Search ADS Lay T. , 2018 . A review of the rupture characteristics of the 2011 Tohoku-oki Mw 9.1 earthquake , Tectonophysics , 733 , 4 – 36 . https://doi.org/10.1016/j.tecto.2017.09.022 Google Scholar Crossref Search ADS Ma C. , Jiang Y. , Li T. , Chen G. , 2016 . Microseismic characterization of brittle fracture mechanism in highly stressed surrounding rock mass , ASEG Extended Abstracts , 2016 ( 1 ), 1 – 5 . https://doi.org/10.1071/ASEG2016ab203 Google Scholar Crossref Search ADS Marfurt K.J. , 1984 . Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations , Geophysics , 49 ( 5 ), 533 – 549 . https://doi.org/10.1190/1.1441689 Google Scholar Crossref Search ADS Maxwell S.C. , Rutledge J. , Jones R. , Fehler M. , 2010 . Petroleum reservoir characterization using downhole microseismic monitoring , Geophysics , 75 ( 5 ), 75A129 – 75A137 . https://doi.org/10.1190/1.3477966 Google Scholar Crossref Search ADS McGarr A. , 2014 . Maximum magnitude earthquakes induced by fluid injection , J. Geophys. Res.: Solid Earth , 119 ( 2 ), 1008 – 1019 . https://doi.org/10.1002/2013JB010597 Google Scholar Crossref Search ADS McGarr A. et al. , 2015 . Coping with earthquakes induced by fluid injection , Science , 347 ( 6224 ), 830 – 831 . https://doi.org/10.1126/science.aaa0494 Google Scholar Crossref Search ADS PubMed McMechan G.A. , 1983 . Migration by extrapolation of time-dependent boundary values , Geophys. Prospect. , 31 ( 3 ), 413 – 420 . https://doi.org/10.1111/j.1365-2478.1983.tb01060.x Google Scholar Crossref Search ADS Melgar D. , Geng J. , Crowell B.W. , Haase J.S. , Bock Y. , Hammond W.C. , Allen R.M. , 2015 . Seismogeodesy of the 2014 Mw6.1 Napa earthquake, California: rapid response and modeling of fast rupture on a dipping strike-slip fault , J. geophys. Res.: Solid Earth , 120 ( 7 ), 5013 – 5033 . https://doi.org/10.1002/2015JB011921 Google Scholar Crossref Search ADS Morse P.M. , Feshbach H. , 1946 . Methods of Theoretical Physics . Technology Press . Nakata N. , Beroza G.C. , 2016 . Reverse time migration for microseismic sources using the geometric mean as an imaging condition , Geophysics , 81 ( 2 ), KS51 – KS60 . https://doi.org/10.1190/geo2015-0278.1 Google Scholar Crossref Search ADS Rutledge J.T. , Phillips W.S. , 2003 . Hydraulic stimulation of natural fractures as revealed by induced microearthquakes, carthage cotton valley gas field, east texas , Geophysics , 68 ( 2 ), 441 – 452 . https://doi.org/10.1190/1.1567214 Google Scholar Crossref Search ADS Sens-Schönfelder C. , Wegler U. , 2006 . Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia , Geophys. Res. Lett. , 33 ( 21 ), L21302 . https://doi.org/10.1029/2006GL027797 Google Scholar Crossref Search ADS Shang X. , de Hoop M.V. , Hilst R.D. , 2012 . Beyond receiver functions: Passive source reverse time migration and inverse scattering of converted waves , Geophys. Res. Lett. , 39 ( 15 ), L15308 . https://doi.org/10.1029/2012GL052289 Google Scholar Crossref Search ADS Shang X. , de Hoop M.V. , van der Hilst R.D. , 2017 . Common conversion point stacking of receiver functions versus passive-source reverse time migration and wavefield regularization , Geophys. J. Int. , 209 ( 2 ), 923 – 934 . https://doi.org/10.1093/gji/ggx069 Google Scholar Crossref Search ADS Sun R. , McMechan G.A. , 2001 . Scalar reverse-time depth migration of prestack elastic seismic data , Geophysics , 66 ( 5 ), 1519 – 1527 . https://doi.org/10.1190/1.1487098 Google Scholar Crossref Search ADS Tromp J. , Tape C. , Liu Q. , 2005 . Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels , Geophys. J. Int. , 160 ( 1 ), 195 – 216 . https://doi.org/10.1111/j.1365-246X.2004.02453.x Google Scholar Crossref Search ADS Virieux J. , 1986 . P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method , Geophysics , 51 ( 4 ), 889 – 901 . https://doi.org/10.1190/1.1442147 Google Scholar Crossref Search ADS Virieux J. , Operto S. , 2009 . An overview of full-waveform inversion in exploration geophysics , Geophysics , 74 ( 6 ), WCC1 – WCC26 . https://doi.org/10.1190/1.3238367 Google Scholar Crossref Search ADS 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 . https://doi.org/10.1785/0120000006 Google Scholar Crossref Search ADS Walker K.T. , Ishii M. , Shearer P.M. , 2005 . Rupture details of the 28 March 2005 Sumatra Mw 8.6 earthquake imaged with teleseismic P waves , Geophys. Res. Lett. , 32 ( 24 ), L24303 . https://doi.org/10.1029/2005GL024395 Google Scholar Crossref Search ADS Wang C. , Cheng J. , Arntsen B. , 2016 . Scalar and vector imaging based on wave mode decoupling for elastic reverse time migration in isotropic and transversely isotropic media , Geophysics , 81 , S383 – S398 . Google Scholar Crossref Search ADS Wang W. , McMechan G.A. , 2015 . Vector-based elastic reverse time migration , Geophysics , 80 ( 6 ), S245 – S258 . Google Scholar Crossref Search ADS Weingarten M. , Ge S. , Godt J.W. , Bekins B.A. , Rubinstein J.L. , 2015 . High-rate injection is associated with the increase in u.s. mid-continent seismicity , Science , 348 ( 6241 ), 1336 – 1340 . https://doi.org/10.1126/science.aab1345 Google Scholar Crossref Search ADS PubMed White D. , 2013 . Seismic characterization and time-lapse imaging during seven years of co2 flood in the weyburn field, saskatchewan, canada , Int. J. Greenhouse Gas Control , 16 , S78 – S94 . Google Scholar Crossref Search ADS Woodward M.J. , Nichols D. , Zdraveva O. , Whitfield P. , Johns T. , 2008 . A decade of tomography , Geophysics , 73 ( 5 ), VE5 – VE11 . https://doi.org/10.1190/1.2969907 Google Scholar Crossref Search ADS Xiao X. , Leaney W.S. , 2010 . Local vertical seismic profiling (VSP) elastic reverse-time migration and migration resolution: Salt-flank imaging with transmitted P-to-S waves , Geophysics , 75 , S35 – S49 . https://doi.org/10.1190/1.3309460 Google Scholar Crossref Search ADS Xu S. , Zhang Y. , Lambar G. , 2010 . Antileakage fourier transform for seismic data regularization in higher dimensions , Geophysics , 75 ( 6 ), WB113 – WB120 . Google Scholar Crossref Search ADS Yan J. , Sava P. , 2008 . Isotropic angle-domain elastic reverse-time migration , Geophysics , 73 ( 6 ), S229 – S239 . Google Scholar Crossref Search ADS Yang J. , Zhang S. , Zhu H. , 2017 . Isotropic elastic wavefields decomposition using fast poisson solvers , SEG Technical Program Expanded Abstracts , 4716 – 4720 . Yang J. , Zhu H. , Huang J. , Li Z. , 2018a . 2D isotropic elastic Gaussian beam migration for common-shot multicomponent records , Geophysics , 83 ( 2 ), S127 – S140 . Google Scholar Crossref Search ADS Yang J. , Zhu H. , Wang W. , Zhao Y. , Zhang H. , 2018b . Isotropic elastic reverse time migration using the phase- and amplitude-corrected vector P- and S-wavefields , Geophysics , 83 ( 6 ), S489 – S503 . Google Scholar Crossref Search ADS 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 . https://doi.org/10.1785/0120020190 Google Scholar Crossref Search ADS Zhang J. , Tian Z. , Wang C. , 2007 . P- and S-wave-separated elastic wave-equation numerical modeling using 2d staggered grid , SEG Technical Program Expanded Abstracts , 2104 – 2109 . Zhang Q. , McMechan G.A. , 2010 . 2D and 3D elastic wavefield vector decomposition in the wavenumber domain for VTI media , Geophysics , 75 ( 3 ), D13 – D26 . Google Scholar Crossref Search ADS Zhang Y. , Wang R. , Chen Y.-T. , 2015 . Stability of rapid finite-fault inversion for the 2014 Mw6.1 South Napa earthquake , Geophys. Res. Lett. , 42 ( 23 ), 10,263 – 10,272 . Google Scholar Crossref Search ADS Zhou H. , 1994 . Rapid three-dimensional hypocentral determination using a master station method , J. Geophys. Res.: Solid Earth , 99 ( B8 ), 15 439 – 15 455 . https://doi.org/10.1029/94JB00934 Google Scholar Crossref Search ADS Zhu H. , 2017 . Elastic wavefield separation based on the Helmholtz decomposition , Geophysics , 82 , S173 – S183 . Google Scholar Crossref Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Locating and monitoring microseismicity, hydraulic fracture and earthquake rupture using elastic time-reversal imaging
JF - Geophysical Journal International
DO - 10.1093/gji/ggy460
DA - 2019-01-01
UR - https://www.deepdyve.com/lp/oxford-university-press/locating-and-monitoring-microseismicity-hydraulic-fracture-and-y70r7SS1UJ
SP - 726
VL - 216
IS - 1
DP - DeepDyve
ER -