33D virtual seismology
Joeri Brackenhoff, Jan Thorbecke and Kees Wapenaar
Department of Geoscience and Engineering, Delft University of Technology,P.O. Box 5048, 2600 GA Delft, The Netherlands
Abstract
We create virtual sources and receivers in a 3D subsurface using the single-sided homogeneous Green’sfunction representation. We employ Green’s functions and focusing functions that are obtained us-ing reflection data at the Earth’s surface, a macro velocity model and the Marchenko method. Thehomogeneous Green’s function is a Green’s function superposed with its time-reversal. Unlike the clas-sical homogeneous Green’s function representation, this approach requires no receivers on an enclosingboundary, however, it does require the source signal to be symmetric in time. We demonstrate that thesingle-sided representation is an improvement over the classical representation by applying the represen-tations to numerical data that are modeled in a complex 3D model. We retrieve responses to virtualpoint sources with an isotropic and with a double-couple radiation pattern and compare the results to adirectly modeled reference result. We also demonstrate the application of the single-sided representationfor retrieving the response to a virtual rupture that consists of a superposition of point sources. Thisis achieved by obtaining the homogeneous Green’s function for each source separately, before they aretransformed to the causal Green’s function, time-shifted and superposed. The single-sided representationis also used to monitor the complete wavefield that is caused by a numerically modeled rupture. How-ever, the source signal of an actual rupture is not symmetric in time and the single-sided representationcan therefore only be used to obtain the causal Green’s function. This approach leaves artifacts in thefinal result, however, these artifacts are limited in space and time.
PACS numbers: a r X i v : . [ phy s i c s . g e o - ph ] S e p . INTRODUCTION Over the past few decades, the amount of induced seismicity has increased and is occurringat locations around the world [10]. While the effects of induced seismicity are often harmful, themeasurements of these events can be used to gain more insight into the mechanics of earthquakerupture [8]. For example, the measurements can be used in an inversion process to obtain the seis-mic moment tensor, which describes the source mechanism of a seismic event [1]. The knowledgeof the moment tensor as well as the location of the source can help to determine what caused theinduced seismicity. These inversions often rely on an accurate velocity model of the subsurfaceto obtain the required wavefields [28], because errors in the velocity model can cause mistakes inthe inversion result [20].A recent development for obtaining accurate wavefields in the subsurface is the homogeneousGreen’s function retrieval method. A homogeneous Green’s function is a Green’s function super-posed with its time-reversal. The classical representation of the homogeneous Green’s functionwas derived in 1970 [17] and was later used for inverse source problems [18] and further researchedfor inverse scattering methods [15], seismic imaging [7] and seismic holography [12]. The classicalrepresentation of the homogeneous Green’s function involves an integral over a closed boundary.In practical situations, data are usually available only on an open boundary. Methods like seis-mic imaging and holography still work well for this situation as long as only primary waves areconsidered. However, internal multiples are incorrectly handled and lead to artifacts when theclassical representation is approximated by an integral along an open boundary.Instead of the classical representation of the homogeneous Green’s function, a single-sidedrepresentation can be used, which is designed to work with an open boundary, typically thesurface of the Earth [26]. This single-sided representation is designed to correctly handle theinternal multiples by employing so-called focusing functions. These focusing functions can beobtained through the use of the Marchenko method, which employs reflection data at the surfaceof the Earth [13]. The single-sided representation has been successfully applied to field data [6].While many applications of the Marchenko method have been performed on 2D data, recentlymore applications on 3D data have been achieved [16, 22]. Especially in areas where there arestrong out-of-plane effects, the 2D approximation on 3D data can cause errors in the result [14].To properly take into account the effects of wave propagation and scattering in 3D, the single-sided retrieval scheme for the homogeneous Green’s function needs to be employed together witha 3D version of the Marchenko method.In this paper, we present the retrieval of the 3D homogeneous Green’s function. We first2eview the classical and single-sided homogeneous Green’s function retrieval schemes and applythe schemes to single source-receiver pairs. We use a 3D Opensource Marchenko method [4] ona synthetic reflection response, that was modeled using a subset of the Overthrust model [3],to create the required Green’s functions and focusing functions for the retrieval schemes. Wedemonstrate the method for point sources that have an isotropic radiation pattern and comparethe retrieved Green’s functions to directly modeled data. Furthermore, we also retrieve snapshotsof wavefields at virtual receivers in 3D to observe the propagation of the wavefield through themedium over time. Aside from considering an isotropic radiation pattern, we also consider thenon-isotropic double-couple radiation pattern, which describes the seismic response to a pure shearfault [1]. Furthermore, we consider the retrieval of a response caused by a rupture in the subsurfaceby employing a series of superposed point sources with varying amplitudes and activation timesand a double-couple radiation pattern, similar to previous research on 2D data [5], but extendedto a 3D medium. For this latter situation we use two different approaches. One is a one-stepprocess, where we assume that we measure the response from the rupture directly, so that we canmonitor the wavefield as it propagates through the subsurface. Hence, in this one-step process wecreate virtual receivers to monitor the response to a real source. The other is a two-step process,where we use the Marchenko method to obtain the homogeneous Green’s function for each virtualsource point separately, and superpose them after each homogeneous Green’s function has beenobtained. Hence, in this two-step process we create virtual receivers and virtual sources. This isa way to forecast the wavefield that would be caused by the rupture, given the properties of therupture and reflection data at the surface. We illustrate the methods with numerical examples.When we speak, for the sake of argument, of measurements of the response to a real source, inthe examples these measurements are simulated by numerical modeling.
II. 3D VIRTUAL SEISMOLOGYA. Wavefields
We consider a Green’s function, G = G ( x , x A , t ), which describes the response of a mediumat time t and position x = ( x , x , x ) . due to an impulsive point source at x A , using a Cartesiancoordinate system. In the coordinate system that we use, the third principal direction pointsdownwards. The Green’s function is the solution to the following acoustic wave equation: ∂ i ( ρ − ∂ i G ) − κ∂ t G = − δ ( x − x A ) ∂ t δ ( t ) , (1)3here ρ = ρ ( x ) is the density of the medium in kg m − , κ = κ ( x ) is the compressibility in kg − m s , ∂ i = ∂/∂x i is the component of a vector consisting of the partial differential operatorsin the three principal directions of the coordinate system, ∂ t = ∂/∂t is the temporal partialdifferential operator and δ ( · ) is a Dirac delta function. In case of repeating subscripts, Einstein’ssummation convention applies. The Green’s function is causal; i.e. G ( x , x A , t ) = 0 for t < G ( x , x A , t ) = G ( x A , x , t ). Because the wave equation for the Green’s function contains a temporalderivative in the source term, the source is defined as a volume injection rate source.We also consider the homogeneous Green’s function G h = G h ( x , x A , t ), which is defined as G h ( x , x A , t ) = G ( x , x A , t ) + G ( x , x A , − t ) , (2)where G ( x , x A , − t ) is the time-reversed Green’s function, which is acausal; i.e. G ( x , x A , − t ) = 0for t >
0, hence, it propagates towards the source. By combining Equations (1) and (2), weobtain the acoustic homogeneous wave equation ∂ i ( ρ − ∂ i G h ) − κ∂ t G h = 0 , (3)where the right hand side vanishes, because the source term on the right hand side of Equation(1) contains a temporal derivative, hence, the wave equation for the time reversal of the Green’sfunction causes the source term to obtain the opposite sign.In this paper, we will make use of the frequency domain versions of the Green’s function andother quantities. A time-dependent function u ( x , t ) is related to the frequency-dependent function u ( x , ω ) by u ( x , ω ) = (cid:90) ∞−∞ u ( x , t ) e iωt d t, (4a) u ( x , t ) = 12 π (cid:90) ∞−∞ u ( x , ω ) e − iωt d ω, (4b)where ω is the angular frequency in rad s − and i is the imaginary unit. Using Equation (4a),the Green’s function can be transformed to the frequency domain, for the sake of efficientlyperforming certain operations, such as convolution. The data that are considered in this paperare band-limited and therefore we define a pressure wavefield p ( x , x A , t ), which is related to theGreen’s function by p ( x , x A , t ) = (cid:90) ∞−∞ G ( x , x A , t − t (cid:48) ) s ( t (cid:48) )d t (cid:48) , (5a) p ( x , x A , ω ) = G ( x , x A , ω ) s ( ω ) , (5b)4here s ( t ) and s ( ω ) are the time domain and frequency domain versions of a source signal. Wealso define a homogeneous pressure wavefield, similar to Equation (2), p h ( x , x A , t ) = (cid:90) ∞−∞ G h ( x , x A , t − t (cid:48) ) s ( t (cid:48) )d t (cid:48) , (6a) p h ( x , x A , ω ) = G h ( x , x A , ω ) s ( ω ) = { G ( x , x A , ω ) + G ∗ ( x , x A , ω ) } s ( ω )= 2 (cid:60){ G ( x , x A , ω ) } s ( ω ) , (6b)where (cid:60) indicates the real part of a complex function and * indicates complex conjugation. Notethat in Equation (6), we have defined that homogeneous wavefield as the convolution of thesource wavelet with the homogeneous Green’s function; i.e. the Green’s function is superposedwith its time-reversal before the convolution. If the Green’s function is convolved with a waveletbefore the superposition is applied, the time-reversal will affect the source wavelet as well. Onlyif s ( t ) = s ( − t ) and hence s ( ω ) = s ∗ ( ω ) can the convolution be applied before the superposition.In other words, { G ( x , x A , ω ) + G ∗ ( x , x A , ω ) } s ( ω ) = G ( x , x A , ω ) s ( ω ) + G ∗ ( x , x A , ω ) s ∗ ( ω ) , (7a) p h ( x , x A , ω ) = p ( x , x A , ω ) + p ∗ ( x , x A , ω ) , (7b)only holds if the source spectrum s ( ω ) is purely real-valued. B. Homogeneous Green’s function retrieval x x x Homogeneous half-space G ( x , x (1) B , ω ) f +1 ( x , x A , ω ) −{ f − ( x , x A , ω ) } ∗ G ( x A , x (1) B , ω ) D ∂ D ∂ D A • x (1) B • x A x (a) Homogeneous half-space D ∂ D ∂ D A G ( x , x (2) B , ω ) f +1 ( x , x A , ω ) −{ f − ( x , x A , ω ) } ∗ G ( x A , x (2) B , ω ) +2 i ={ f ( x (2) B , x A , ω ) } • x (2) B • x A x (b) FIG. 1: Schematic representation of the single-sided Green’s function retrieval scheme from Equations(10) and (11). (a) Retrieval of G ( x A , x B , ω ) when the receiver location is located above the source and noartifacts are present and (b) retrieval of G ( x A , x B , ω ) when the virtual receiver location is located belowthe virtual source and artifacts are present in the form of 2 i (cid:61){ f ( x B , x A , ω ) } . The figure is adaptedfrom a previous publication [5]. G h ( x A , x B , ω ) = (cid:73) ∂ D − iωρ ( x ) (cid:110) ∂ i G ∗ ( x A , x , ω ) G ( x , x B , ω ) − G ∗ ( x A , x , ω ) ∂ i G ( x , x B , ω ) (cid:111) n i d x , (8)where n i is the i th component of the outward pointing normal vector on ∂ D . In Equation (8), G ( x , x B , ω ) describes the response to a source at x B , inside the medium in D , at location x on a boundary ∂ D , which encloses the medium. G ∗ ( x A , x , ω ) back propagates this responsefrom location x at the boundary to receiver location x A inside D . This creates the response G h ( x A , x B , ω ), with a source at location x B and a receiver at location x A . The main practicaldisadvantage of this approach is that a closed boundary around the medium is required, which isusually not feasible for seismological applications. More realistically, the boundary will be openand situated on a single side of the medium, which is often the surface of the Earth. In this case,the representation is approximated as G h ( x A , x B , ω ) ≈ (cid:90) ∂ D iωρ (cid:110) ∂ G ∗ ( x A , x , ω ) G ( x , x B , ω ) (cid:111) d x , (9)where ρ is the density at a horizontal single open boundary ∂ D and we used n = (0 , , − ∂ D is homogeneous. Applying the representationin this way introduces significant artifacts in the homogeneous Green’s function [6].In more recent years, the homogeneous Green’s function representation has been adjusted totake into account the single-sided open boundary [26]. The scheme that is used in this paper istaken from previous research [5, Equations (10) and (11)], G ( x A , x B , ω ) + χ ( x B )2 i (cid:61){ f ( x B , x A , ω ) } = (cid:90) ∂ D iωρ G ( x , x B , ω ) ∂ (cid:16) f +1 ( x , x A , ω ) − { f − ( x , x A , ω ) } ∗ (cid:17) d x , (10)where f ( x B , x A , ω ) = f +1 ( x B , x A , ω ) + f − ( x B , x A , ω ), (cid:61) denotes the imaginary part of a complexfunction and χ ( x B ) is a characteristic function that is defined as χ ( x B ) = , for x B in D , , for x B on ∂ D = ∂ D ∪ ∂ D A , , for x B outside D ∪ ∂ D , (11)6here ∂ D A is a horizontal open boundary inside the subsurface of the Earth at the same depthas x A . The medium in D is assumed to be lossless and evanescent waves are ignored. Note, thatin Equation (10), we retrieve the causal Green’s function instead of the homogeneous Green’sfunction. In this representation, no time-reversed Green’s function is employed, but rather thedecomposed focusing functions f +1 ( x , x A , ω ) and f − ( x , x A , ω ) are used, where the superscripts +and − indicate a downgoing and upgoing wavefield, respectively. These focusing functions aredesigned to focus from a single-sided open boundary ∂ D to a location x A inside the subsurfaceof the Earth, generally referred to as the focal location, without artifacts caused by multiplescattering in the overburden. The downgoing focusing function is defined as the inverse of thetransmission response of a medium that is truncated below the focal location [21, 27]. In Equation(10), the focusing functions f +1 ( x , x A , ω ) and f − ( x , x A , ω ) operate in a similar way as the time-reversed Green’s function G ∗ ( x A , x , ω ) does in Equation (9), backpropagating the response fromthe boundary ∂ D to location x A . The main difference is that unlike Equation (9), Equation (10)is specifically designed for application to the open boundary.The representation in Equation (10) does have an issue on the left hand side of the equation inthe form of the term χ ( x B )2 i (cid:61){ f ( x B , x A , ω ) } . Depending on the relative locations of the receiver x A and the source x B , as formulated by Equation (11), artifacts related to the focusing functionbetween the two locations are introduced in the obtained Green’s function. This is schematicallyshown in Figure 1, where in (a) the receiver is located above the source and the Green’s functionis retrieved without artifacts. When the virtual source is located at the same depth level orabove the virtual receiver, artifacts are present in the retrieved Green’s function. By combiningEquations (10) and (2), we obtain the single-sided retrieval scheme for the homogeneous Green’sfunction [25, Equation (33)]: G h ( x A , x B , ω ) = 4 (cid:60) (cid:90) ∂ D iωρ G ( x , x B , ω ) ∂ (cid:16) f +1 ( x , x A , ω ) − { f − ( x , x A , ω ) } ∗ (cid:17) d x . (12)Equation (12) expresses the retrieval of the homogeneous Green’s function between two locationsin the subsurface using a single-sided boundary, without any artifacts from the focusing function. C. Implementation of Green’s function retrieval
We will demonstrate the results of the retrieval schemes in Equations (9), (10) and (12) withnumerical examples. In order to obtain the required Green’s functions and focusing function,we employ the 3D Marchenko method on acoustic reflection data [27]. The scheme allows one7 (m) x (m) x ( m ) V e l o c i t y ( m s − ) D e n s i t y ( k g m − ) x (m) x (m) T i m e ( m s ) A m p lit ud e (-) A m p lit ud e (-) RickerConvolution 0.00.51.0 | A m p lit ud e | (-) | A m p lit ud e | (-) RickerConvolution (a) (b)(c) (d)(e) (f)
FIG. 2: (a) Velocity and density model of the subsurface, based on the Overthrust model by [3]. (b)Common-source record of a source, located at (0,0,0)m, recorded at the surface of the model in (a).Wavelet with a flat spectrum in (c) the time domain and (d) the frequency domain. Ricker wavelets in(e) the time domain and (f) the frequency domain. The common-source record in (b) is modeled usingthe black wavelet in (c) and (d). The gray dashed wavelet in (e) and (f) is the result of the temporalconvolution of the black wavelet in (c) and (d) with the black wavelet in (e) and (f). to retrieve the Green’s function and focusing function between receivers at the surface of theEarth and a focal location in the subsurface of the Earth. To obtain these functions, a reflectionresponse without surface related multiples at the surface of the Earth is required, as well as anestimation of the direct arrival from the surface of the Earth to the focal location. Usually, thetime-reversed direct arrival of the Green’s function from the focal location to the surface is usedfor this, even though this introduces errors proportional to the transmission losses into the final8esult [24].In this paper, we make use of an opensource 3D implementation of the Marchenko method [4].To obtain the reflection data, we use a 3D finite-difference modeling code [23] together with asubset of the 3D Overthrust model [3], which is shown in Figure 2(a). To ensure strong reflections,the same model is used for the density and velocity values. To model the data, a fixed-spreadacquisition is utilized, where a source is modeled at every receiver location. The source/receiverlocations vary from -2250 to 2250m in the inline ( x ) direction, with a spacing of 25m, whilethe locations in the crossline ( x ) direction vary from -1250 to 1250m, with a spacing of 50m.We define a common-source record as the reflection response to a fixed source, observed by allreceivers. The recording length of each common-source record is 4.0s with a temporal samplingof 4ms. An example of a common-source record is shown in Figure 2(b). The data are modeledusing a wavelet with a flat spectrum, shown as the black wavelet in Figure 2(c) and (d). Examplesof Green’s functions and a focusing function obtained from these data can be found in Figure 3.Once we obtain the required Green’s functions and focusing functions, we use them in the var-ious retrieval schemes. The schemes so far have been formulated assuming we have the analyticalrepresentations instead of the numerical data that we are actually dealing with. In reality, wewill not have impulse responses, but rather wavefields with a band-limited signal as defined byEquations (5a) and (5b). To account for this, we use numerical approximations of the schemesand make use of pressure wavefields with a band-limited source signature. We rewrite Equations(9), (10) and (12) as p h ( x A , x B , ω ) ≈ n R (cid:88) j iωρ (cid:110) ∂ G ∗ ( x A , x j , ω ) p ( x j , x B , ω ) (cid:111) ∆ x j , (13) p ( x A , x B , ω ) + χ ( x B )2 is ( ω ) (cid:61){ f ( x B , x A , ω ) } = n R (cid:88) j iωρ p ( x j , x B , ω ) ∂ (cid:16) f +1 ( x j , x A , ω ) − { f − ( x j , x A , ω ) } ∗ (cid:17) ∆ x j , (14) p h ( x A , x B , ω ) = 4 (cid:60) n R (cid:88) j iωρ p ( x j , x B , ω ) ∂ (cid:16) f +1 ( x j , x A , ω ) − { f − ( x j , x A , ω ) } ∗ (cid:17) ∆ x j , (15)where x j is the location of the j th receiver at the surface of the Earth, n R is the amount ofreceivers and ∆ x j indicates the receiver sampling distance. While ∆ x j can be unique for eachreceiver position, in our fixed spread acquisition the value is the same for all receivers, namely∆ x j = ∆ x ∆ x = 25 . · . G h ( x A , x B , ω ), G ( x A , x B , ω ) and G ( x , x B , ω ) by p h ( x A , x B , ω ), p ( x A , x B , ω ) and p ( x , x B , ω ), respectively, while some of the other quantities are still denoted by their originalsymbol. In the application of Equations (13)-(15), we assume that p ( x , x B , ω ) is obtained either9 (m) x (m) T i m e ( m s ) x (m) x (m) T i m e ( m s ) x (m) x (m) T i m e ( m s ) x (m) x (m) T i m e ( m s ) (a) (b)(c) (d) FIG. 3: Examples of wavefields obtained using the Marchenko method. (a) Green’s function G ( x A , x j , ω )and (b) focusing function f +1 ( x j , x A , ω ) −{ f − ( x j , x A , ω ) } ∗ convolved with the black flat spectrum waveletfrom Figure 2(c), with x A = ( − , , p ( x j , x B , ω ), i.e., a Green’s functionconvolved with the dashed gray Ricker wavelet from Figure 2(e), with (c) an isotropic source and (d)pressure wavefield D θB { p ( x j , x B , ω ) } with a double-couple source, both with x B = (500 , − , through the use of the Marchenko method (in the two-step method) or by a direct measurement(in the one-step method), while G ( x A , x , ω ), f +1 ( x , x A , ω ) and f − ( x , x A , ω ) are always obtainedthrough the use of the Marchenko method. Therefore, we can control the source spectrum of thedata that are used to generate the virtual receiver data. We ensure that G ( x A , x , ω ), f +1 ( x , x A , ω )and f − ( x , x A , ω ) have a source signature with a flat spectrum of amplitude 1.0 for a certainfrequency range, so that the convolution with a unique source signature in that frequency rangewill produce the response to the latter source signal. This is schematically shown in Figures2(c)-(f). The reflection data in Figure 2(b) are convolved with the black wavelet from Figures2(c) and (d), which has a flat spectrum in a frequency range of 5Hz to 25Hz, and the estimation10f the direct arrival is modeled with the same wavelet. The result of the convolution of thesewavefields is a wavefield that contains a wavelet very similar to the original flat spectrum wavelet.If a wavefield with this source signature is convolved with a wavefield that contains a differentsource spectrum, for example, the Ricker wavelet that is shown as the black wavelet in Figures2(e) and (f), the resulting wavefield will contain a wavelet that is almost identical to the originalRicker wavelet, as is shown by the gray dashed line in Figures 2(e) and (f). There is some slightattenuation of the lowest and highest frequencies, as can be seen in Figure 2(f), however, thishas little effect on the wavelet shape in Figure 2(e). The versions of p h ( x A , x B , ω ), p ( x A , x B , ω )and p ( x , x B , ω ) that are used in Equations (13)-(15) all include the original Ricker wavelet fromFigure 2(e); an example of such a pressure wavefield can be found in Figure 3(c), with its sourceat location x B = (500 , − , x A = ( − , , p ( x A , x B , ω ) is purely real valued,which holds for the source spectrum of the zero-phase Ricker wavelet.To demonstrate the validity of our implementation, we show the result of the retrieval schemesin Figure 4, using the two-step method. Each column corresponds to a different pair of virtualsource and virtual receiver positions, while each row corresponds to a different type of retrievalmethod. The first column has a virtual receiver located above the virtual source and the positionsonly differ in depth. In the second column the virtual receiver is located below the virtual sourceand the locations differ in both the inline direction and depth. For the third column the virtualreceiver is located below the virtual source and the locations differ in all three principal directions.The required Green’s function and focusing function are obtained using the Marchenko methodand a first arrival that was obtained by modeling in the exact medium. We invert the first arrivalinstead of only time-reversing it, to avoid the transmission losses. While this is not a realisticscenario, as for field data we would not be able to use the exact model, we wish to demonstratethat the method is, at least in theory, capable of obtaining the exact amplitudes. The source hasan isotropic radiation pattern. For each panel, the result that is obtained through the use of aretrieval scheme is plotted in dashed gray, while a directly modeled reference solution is shownin solid black.The homogeneous wavefield that is obtained using Equation (13) is shown in the first row ofFigure 4. Both the Green’s function for the virtual source and the virtual receiver were obtainedusing the Marchenko method. For all location pairs, the results are poor. While the order of11 P r e ss u r e ( P a ) P r e ss u r e ( P a ) P r e ss u r e ( P a ) | P r e ss u r e | (-) (a) x B = (0 , , x A = (0 , , p h retrieval (b) x B = (500 , , x A = ( − , , p h retrieval (c) x B = (500 , − , x A = ( − , , p h retrieval (d) x B = (0 , , x A = (0 , , p retrieval (e) x B = (500 , , x A = ( − , , p retrieval (f) x B = (500 , − , x A = ( − , , p retrieval (g) x B = (0 , , x A = (0 , , p h retrieval FD (h) x B = (500 , , x A = ( − , , p h retrieval FD (i) x B = (500 , − , x A = ( − , , p h retrieval FD (j) x B = (0 , , x A = (0 , , p h retrieval Eikonal (k) x B = (500 , , x A = ( − , , p h retrieval Eikonal (l) x B = (500 , − , x A = ( − , , p h retrieval Eikonal FIG. 4: Green’s funcions of pairs of virtual sources and virtual receivers for different locations anddifferent types of retrieval scheme. The solid black lines are the exact (directly modeled) Green’s functionsand the dashed gray lines are the retrieved functions. Each column corresponds to a different pair oflocations. The first row corresponds to the classical retrieval scheme of Equation (13), the second rowto the Marchenko retrieval scheme of Equation (14) and the third row to the homogeneous Marchenkoretrieval scheme of Equation (15). For all of these rows, the first arrival required for the Marchenkomethod is obtained using finite-difference modeling. For the fourth row the same retrieval method isused as in the third row, except the first arrival is obtained using an Eikonal solver, instead of finite-difference modeling. The traces in the final row are all normalized. All traces contain the Ricker waveletfrom Figure 2(e). magnitude of the retrieved wavefield is similar to that of the direct modeling, the exact amplitudeshave a strong mismatch and there are artifacts present at all times. The exceptions are the firstarrival in Figure 4(b) and the early coda in Figure 4(c). Aside from these events however, allother events are wrong and there are still significant artifacts for both examples.The second row shows the pressure wavefield that is obtained using the open boundary retrievalscheme from Equation (14). When the source is located below the virtual receiver, as is the casein Figure 4(d), the result shows a good match to the reference solution in both amplitude andarrival time. Because the first arrival is isolated, we apply a muting window to the data before thefirst arrival to remove numerical artifacts. When the source is located above the receiver, however,the result degrades in quality. There are strong artifacts present in the result at times before thefirst arrival and the first arrival has the wrong polarity and amplitude. As these artifacts are12f a similar magnitude as the first arrival, we cannot apply the muting window. The retrievedcoda in these two latter cases is still accurate, with some slight mismatch in the amplitude of theevents. This is caused by the different lateral positions of the source and receiver. The apertureand sampling of the data in both the inline and crossline direction are limited, so the exact eventsbecome harder to obtain. The overall coda shows a good match to the reference.To improve the results of the retrieval, the representation from Equation (15) is used to retrievethe homogeneous Green’s function, shown in the third row of Figure 4. The results in Figure 4(d)and (g) are identical, which corresponds to the condition in Equation (11). The improvementis apparent when the source is located above the receiver as is the case in Figure 4(h) and (i).Compared to the results in Figure 4(e) and (f), the unwanted artifacts are removed and the firstarrival is retrieved properly. Here, we once again apply the muting window before the first arrival.The amplitude mismatch in the coda is still present, indicating that this is a limitation causedby the aperture of the recording array and not of the type of retrieval method.Finally, in the bottom row of Figure 4, we apply Equation (15) again, however, this time thefirst arrivals used in the Marchenko method are obtained using an Eikonal solver in a smoothedversion of the velocity model. This is to simulate a more realistic situation, where accurate modelinformation would not be available. Because the exact amplitude of the first arrival cannot beobtained in this case, the retrieved homogeneous Green’s function is normalized and compared toa normalized version of the reference solution. This is intended to show that even when the exactamplitude cannot be obtained, the relative amplitude can be properly obtained. The matchesfor all three source-receiver pairs are good, but of a lesser quality than when the finite-differencemodeling is employed. Due to the complexity of the model, as well as the smoothing, the Eikonalsolver can encounter issues with obtaining the correct arrival times. Furthermore, we only usean estimation of the amplitude distribution along the wavefront, which also will not properlyrepresent the true effect that the subsurface would have on the amplitude. However, the resultsstill support that use of an Eikonal solver for 3D media can yield useful results.
D. Visualization of the 3D results
While the traces in Figure 4 demonstrate the validity of our approach, they are limited inscope. To further test our approach in 3D, we obtain the results for not just a single source-receiver pair. Instead we retrieve a large amount of focusing functions and use these in Equation(15) to visualize the retrieved Green’s functions evolving in time through the 3D medium. Toobtain the results, we use the approach employing the Eikonal solver, similarly to how we obtained13 (m) x (m) x ( m ) (a) x (m) x (m) x ( m ) (b) x (m) x (m) x ( m ) (c) x (m) x (m) x ( m ) (d) x (m) x (m) x ( m ) (e) x (m) x (m) x ( m ) (f) x (m) x (m) x ( m ) (g) x (m) x (m) x ( m ) (h) t=0ms t=0mst=300ms t=300mst=600ms t=600mst=900ms t=900ms FIG. 5:
3D snapshots of the homogeneous Green’s function retrieval of an isotropic virtual source in the Overthrust model at (a)0ms, (b) 300ms, (c) 600 ms and (d) 900ms. For comparison, (e), (f), (g) and (h) are snapshots of a directly modeled homogeneouspressure wavefield at 0ms, 300ms, 600ms and 900ms, respectively. The source is located at x B = (0 , , p ( x j , x B , t ), we obtain a single pressure wavefield due to a source withan isotropic radiation pattern at x B = (0 , , III. MOMENT TENSOR MONITORINGA. Non-isotropic point source
In reality, an event in the subsurface is seldom generated by an isotropic point source. Insteadthe source wavefield is often caused by faulting, the mechanism of which can be described by amoment tensor [1], which causes the amplitude along the wavefront to vary. The double-couplesource mechanism is often used, which is a moment tensor that describes a pure shear fault, byits strike, rake and dip [11]. In previous work, the double-couple source mechanism was combinedwith the Marchenko method to obtain the virtual response of a double-couple point source in thesubsurface, as well as that of a rupture plane [6]. Here, we wish to demonstrate that similar resultscan be achieved in 3D. We repeat the examples of the isotropic point source, using Equations(13)-(15), however, we replace the isotropic source at x B in p ( x j , x B , ω ) by a double-couple sourcegenerated by a moment tensor. We use an operator D θB {·} , which transforms the radiation patternof the source at x B from an isotropic radiation pattern to a double-couple radiation pattern. Itis defined as D θB {·} = ( θ (cid:107) i + θ ⊥ i ) ∂ i,B , (16)where ∂ i,B is a component of the vector containing the partial derivatives acting on the monopolesignal originating from source location x B , which alters the radiation pattern, θ (cid:107) i is a componentof a vector that orients one couple of the signal parallel to the fault plane and θ ⊥ i is a componentof a vector that orients the other couple perpendicular to the fault plane. Because we are dealingwith acoustic reflection data, we only model the P-waves of the double-couple source, select thefirst arrival and use it in the Marchenko method to obtain the desired virtual double-coupleresponse D θB { p ( x j , x B , ω ) } . An example of such a wavefield can be found in Figure 3(d), whichhas a source at the same position as the pressure wavefield with an isotropic source in (c). Thiswavefield is then used in Equation (13) or (15) to obtain D θB { p h ( x A , x B , ω ) } or in Equation (14)to obtain D θB { p ( x A , x B , ω ) } . Previous research has suggested that the double-couple source isnot always a sufficient description of an earthquake source [9]. Our wavefield retrieval methodis valid for any type of moment tensor, however, for the sake of simplicity, we stick with thedouble-couple representation. 16 P r e ss u r e ( P a ) P r e ss u r e ( P a ) P r e ss u r e ( P a ) | P r e ss u r e | (-) (a) x B = (0 , , x A = (0 , , p h retrieval (b) x B = (500 , , x A = ( − , , p h retrieval (c) x B = (500 , − , x A = ( − , , p h retrieval (d) x B = (0 , , x A = (0 , , p retrieval (e) x B = (500 , , x A = ( − , , p retrieval (f) x B = (500 , − , x A = ( − , , p retrieval (g) x B = (0 , , x A = (0 , , p h retrieval FD (h) x B = (500 , , x A = ( − , , p h retrieval FD (i) x B = (500 , − , x A = ( − , , p h retrieval FD (j) x B = (0 , , x A = (0 , , p h retrieval Eikonal (k) x B = (500 , , x A = ( − , , p h retrieval Eikonal (l) x B = (500 , − , x A = ( − , , p h retrieval Eikonal FIG. 6: As Figure 4, but for a virtual double-couple source with a strike, rake and dip of 19, 68 and 25degrees, respectively.
For our example, we use a double-couple source with a strike, rake and dip of, 19, 68 and 25degrees, respectively, and obtain the response between the virtual source-receiver pairs, similarto the examples in Figure 4, using the same color scheme as in that figure. The results are shownin Figure 6, at the same locations as the source-receiver pairs for the isotropic results. Each rowshows a different retrieval method, while each column shows a different source-receiver pair. Thefirst column shows very comparable results to those of the first column of Figure 4. The amplitudeand shape of the events are different, caused by the different source mechanism. Once again, theclassical retrieval using Equation (13) shows poor results, while the results of the Green’s functionretrieval using Equation (14) and the homogeneous Green’s function retrieval using Equation (15)are identical, because the virtual receiver is located above the virtual source, and the match tothe direct modeling is good.If the second and third column are considered, the retrieval using Equation (14) shows strongerrors around the first arrival, while the coda is retrieved properly. When the homogeneousGreen’s function retrieval from Equation (15) is used, the result improves around the first arrival,as can be seen in the bottom two rows of Figure 6. There are still some errors in the early coda.The results in the final column show that parts of the coda are also not properly obtained. Dueto the more complex source signature, it becomes harder for the method to resolve all events,especially with the limited aperture. The Eikonal solver once again can only obtain the relativeamplitudes, and has similar issues with retrieving the proper events in the coda. Overall theresults are encouraging. 17o further investigate the effects of using a double-couple source mechanism, we retrieve thewavefield for the same three 3D slices as we did for Figure 5. The result is shown in Figure 7, wherethe left column shows the result for the retrieval using the Eikonal solver and Equation (15), whilea direct modeling is shown in the right column for comparison. Note that the direct modelingcontains a source artifact caused by the modeling of only the P-waves. When the results of theretrieval and the direct modeling are compared, most of the nearly vertically traveling eventsare properly retrieved, not only in arrival time, but also in polarity. For events traveling nearlyhorizontally, the retrieval is once again poorer. Overall, the result using the double-couple sourcehave a similar quality as the results for the isotropic source, which demonstrates that in 3D,the double-couple source can be successfully integrated into the homogeneous Green’s functionretrieval.
B. Rupture
In previous sections, we have only considered point sources, however, in the field, an earthquakeis seldom a single event, rather, it consists of a cluster of several events that are activated overan area for a period of time [19]. Hence, the total wavefield of an earthquake is not the resultof a single instantaneous source, instead, it consists of a superposition of wavefields caused bydifferent sources that are activated at different times. To approximate this kind of wavefield, wedefine a total wavefield P ( x A , t ) that consists of a superposition of wavefields that are caused bydouble-couple point sources. The superposition can be expressed as P ( x A , t ) = n S (cid:88) k D θ, ( k ) B { p ( x A , x ( k ) B , t ) } = n S (cid:88) k (cid:90) ∞−∞ D θ, ( k ) B { G ( x A , x ( k ) B , t − t (cid:48) ) } s ( k ) ( t (cid:48) )d t (cid:48) , (17a) P ( x A , ω ) = n S (cid:88) k D θ, ( k ) B { p ( x A , x ( k ) B , ω ) } = n S (cid:88) k D θ, ( k ) B { G ( x A , x ( k ) B , ω ) } s ( k ) ( ω ) , (17b)where x ( k ) B indicates the location of the k th source of a total of n S sources, D θ, ( k ) B {·} is the double-couple operator for each location and s ( k ) ( t ) is the corresponding source signal for each locationthat contains all the information for the source strength, activation time and duration. Becauseof the different activation times, the source spectrum of P ( x A , t ) is no longer purely real-valuedand can therefore not be used in Equation (15). However, using it in Equation (14) is still valid,18 (m) x (m) x ( m ) (a) x (m) x (m) x ( m ) (b) x (m) x (m) x ( m ) (c) x (m) x (m) x ( m ) (d) x (m) x (m) x ( m ) (e) x (m) x (m) x ( m ) (f) x (m) x (m) x ( m ) (g) x (m) x (m) x ( m ) (h) t=0ms t=0mst=300ms t=300mst=600ms t=600mst=900ms t=900ms FIG. 7: As Figure 5, but for a virtual double-couple source with a strike, rake and dip of 19, 68 and 25degrees, respectively. Note that the direct modeling contains a source artifacts caused by only modelingthe P-waves.
19s no time-reversal is applied. We rewrite Equation (14) for this purpose as P ( x A , ω ) + n S (cid:88) k D θ, ( k ) B { χ ( x ( k ) B )2 is ( k ) ( ω ) (cid:61){ f ( x ( k ) B , x A , ω ) }} = n R (cid:88) j iωρ P ( x j , ω ) ∂ (cid:16) f +1 ( x j , x A , ω ) − { f − ( x j , x A , ω ) } ∗ (cid:17) ∆ x j = n R (cid:88) j iωρ n S (cid:88) k D θ, ( k ) B { p ( x j , x ( k ) B , ω ) } ∂ (cid:16) f +1 ( x j , x A , ω ) − { f − ( x j , x A , ω ) } ∗ (cid:17) ∆ x j . (18)In Equation (18), we can retrieve P ( x A , ω ), however, we will also obtain the focusing functionartifacts that are related to each source position, below the source depth. As there are multiplesources, that can have different depths, the artifacts related to one source can interfere withthe part of the signal that originates from a deeper part of the medium. Consequently, onlyabove the shallowest source depth can we expect to obtain the correct wavefield at all times. Fordeeper parts of the medium, we expect to retrieve artifacts before and around the first arrivaltime and the correct coda at later times, similar to the results that were shown in Figures 4 and6. Obtaining the wavefield in this way is a one-step process, where we first measure the totalwavefield of an actual rupture and use it in combination with the focusing functions, obtainedthrough the Marchenko method, to monitor the subsurface with virtual receivers.On the other hand, to obtain the response to a virtual rupture, we can use a two-step process toretrieve P ( x A , t ). This has a significant advantage over the one-step process. Instead of measuringthe resulting wavefield of the superposed sources, we use the Marchenko method to retrieve theindividual wavefields D θ, ( k ) B { p h ( x A , x ( k ) B , ω ) } related to each source position. In this case we do notmeasure the total wavefield, but predict it by using the Marchenko method to obtain the sourcewavefield D θ, ( k ) B { p ( x j , x ( k ) B , ω ) } before using it in Equation (15). Because of this, we can ensurethat the source spectrum of D θ, ( k ) B { p ( x j , x ( k ) B , ω ) } for each individual virtual source is purely real-valued, before we apply the time-reversal. The wavefields that we retrieve in this way are free ofthe artifacts related to the focusing function and can be combined to form P ( x A , t ): P ( x A , t ) = n S (cid:88) k H ( t − t ( k ) ) D θ, ( k ) B { p h ( x A , x ( k ) B , t − t ( k ) ) } , (19)where H is the Heaviside function and t ( k ) is the activation time of the source. In Equation (19),we shift the signals in time by t ( k ) before superposition is applied. Because these wavefields aretime-shifted and homogeneous, i.e. they contain time-shifted versions of D θ, ( k ) B { p ( x A , x ( k ) B , t ) } and D θ, ( k ) B { p ( x A , x ( k ) B , − t ) } , the acausal part of one wavefield may interfere with the causal part ofanother wavefield. The Heaviside function is applied to remove all acausal parts of the wavefieldsto avoid such an issue. While this approach cannot be used for the monitoring of wavefields20easured in the field that are caused by sources that are active over a period of time, the approachcan be used to forecast the total wavefield of a virtual rupture, given a specific distribution ofsources.To demonstrate the monitoring and forecasting of the total wavefield, we consider a ruptureplane in the Overthrust model that consists of a cluster of 61 point sources with a double-couple radiation pattern and that are activated at different points in time. Instead of retrievingwavefields that contain the zero-phase wavelet, like we have done in the previous examples, weretrieve wavefields that contain a unique causal wavelet for each source position, as wavefields inthe real subsurface will be causal and not zero-phase. We choose the Berlage wavelet, which isdefined as [2]: W ( t ) = AH ( t ) t n e − αt cos(2 πf t + φ ) , (20)where A is the amplitude of the wavelet. The time exponent n , exponential decay factor α ,initial phase angle φ and peak frequency f control the shape of the wavelet. To ensure that thewavelet has an amplitude equal to zero at t = 0, we use an initial phase angle of -90 degrees. Forthe peak frequency, we use the same peak frequency as we used for the Ricker wavelet in Figure2(f). However, for the amplitude, time exponent and exponential decay factor, we take randomvalues, to simulate a heterogeneous region along the rupture plane. The schematic overview forthe rupture simulation can be found in Figure 8. The sources are located along a fault in themodel, where each source has a strike and rake of 90 and 0 degrees, respectively, and is locatedat a fixed crossline position of 0m. The dip of the source is dictated by the fault orientation ateach source location. Figure 8(a) contains the locations of the sources, while Figure 8(b) showsthe activation time and random amplitude and Figure 8(c) shows the random time exponentand exponential decay factor that are used for the Berlage wavelets. The activation time for thesources is linear, with a time delay of 24ms between the activation of subsequent sources, exceptfor the positions where the depth of the source changes. In these cases the time delay is increasedto 32ms to account for the increase in step size. In this way, we simulate a rupture activatingand propagating along the rupture plane with a velocity of 520m s − .The results of both the one-step and the two-step process can be found in Figure 9. Theleft column of this figure shows the result for the one-step process for monitoring a signal, usingEquation (18), and the right column shows the result for the two-step process for forecasting asignal, using Equation (19). For the monitoring process, we convolve the Berlage wavelets withthe source wavefield before we employ the causal Green’s function retrieval. For the forecastingprocess, we use a flat spectrum wavelet to obtain the individual homogeneous Green’s functions.Because we are creating everything from the data, this is a valid approach. After we obtain these21
00 200 0 200 400 x (m)240026002800 x ( m ) A c ti v a ti on ti m e ( s ) A m p lit ud e (-)
400 200 0 200 400 x (m)1.52.02.5 n (-) (-) (a)(b)(c) FIG. 8: (a) Locations of individual double-couple sources with a strike and rake of 90 and 0 degrees,respectively. The dip is oriented along the fault direction for each location. The slice is located along aconstant crossline position of 0m. (b) Activation time and amplitude in black and gray, respectively, and(c) time exponent n and exponential decay factor α in black and gray, respectively, for computing theBerlage wavelets using Equation (20). The horizontal positions of the sources in (a), (b) and (c) match. IV. CONCLUSIONS
We have shown that the Marchenko method can be applied to 3D reflection data at thesurface of the Earth to obtain the responses for virtual source-receiver pairs in the subsurface ina data-driven way. We did this by considering the 3D single-sided representation for obtainingthe homogeneous Green’s function in the subsurface. For the input reflection data, we modeleda reflection dataset in a subsection of the 3D Overthrust model. We compared the single-sidedapproach to the classical representation and showed for a number of selected source-receiver pairsthat the single-sided representation can obtain accurate results, whereas the result of the classicalrepresentation contains artifacts. The single-sided approach can be applied in two ways, one wherethe causal Green’s function and one where the homogeneous Green’s function is retrieved. Theretrieval of the causal Green’s function is accurate between the surface and the source location.When the virtual receiver is located below the source, artifacts are created that are related to thefocusing function from the source to the virtual receiver. These artifacts are present before thecoda, and therefore affect the first arrival, however, the coda is obtained accurately. When the23 (m) x (m) x ( m ) (a) x (m) x (m) x ( m ) (b) x (m) x (m) x ( m ) (c) x (m) x (m) x ( m ) (d) x (m) x (m) x ( m ) (e) x (m) x (m) x ( m ) (f) x (m) x (m) x ( m ) (g) x (m) x (m) x ( m ) (h) t=0ms t=0mst=640ms t=640mst=1280ms t=1280mst=1920ms t=1920ms FIG. 9:
3D snapshots of the Green’s function retrieval for a wavefield caused by a rupture in the Overthrust model at (a) 0ms,(b) 640ms, (c) 1280 ms and (d) 1920 ms, using Equation (18), and 3D snapshots of superposed and time-shifted wavefields in theOverthrust model, obtained using homogeneous Green’s function retrieval using Equation (19), at (e) 0ms, (f) 640ms, (g) 1280 msand (h) 1920 ms. All wavefields have an overlay of a cross section of the Overthrust model to indicate the locations where we expectscattering to take place. Details about the locations, activation times and the wavelets of each source can be found in Figure 8. unding
This work has received funding from the European Unions Horizon 2020 research and innova-tion program: European Research Council (grant agreement no. 742703).
Copyright notice
This works has been submitted to IEEE Transactions on Geoscience and Remote Sensing.c (cid:13)
1] Aki, K. and Richards, P. G. (2002).
Quantitative Seismology . W.H. Freeman and Company, SanFransisco.[2] Aldridge, D. F. (1990). The Berlage wavelet.
Geophysics , 55(11):1508–1511.[3] Aminzadeh, F., Brac, J., and Kunz, T. (1997). SEG/EAGE 3-D Salt and Overthrust models.SEG/EAGE 3-D modeling series, no. 1: Distribution CD of Salt and Overthrust models.
SEG bookseries .[4] Brackenhoff, J., Thorbecke, J., Koehne, V., Barrera, D., and Wapenaar, K. (2020). Implementationof the 3D Marchenko method. arXiv preprint arXiv:2004.00896 .[5] Brackenhoff, J., Thorbecke, J., and Wapenaar, K. (2019a). Monitoring induced distributed double-couple sources using Marchenko-based virtual receivers.
Solid Earth , 10(1):1301–1319.[6] Brackenhoff, J., Thorbecke, J., and Wapenaar, K. (2019b). Virtual sources and receivers in the realEarth: Considerations for practical applications.
Journal of Geophysical Research: Solid Earth ,124(11):11802–11821.[7] Esmersoy, C. and Oristaglio, M. (1988). Reverse-time wave-field extrapolation, imaging, and inver-sion.
Geophysics , 53(7):920–931.[8] Galis, M., Ampuero, J. P., Mai, P. M., and Cappa, F. (2017). Induced seismicity provides insightinto why earthquake ruptures stop.
Science advances , 3(12):eaap7528.[9] Julian, B. R., Miller, A. D., and Foulger, G. (1998). Non-double-couple earthquakes 1. Theory.
Reviews of Geophysics , 36(4):525–549.[10] Keranen, K. M. and Weingarten, M. (2018). Induced seismicity.
Annual Review of Earth andPlanetary Sciences , 46(1):149–174.[11] Li, D., Helmberger, D., Clayton, R. W., and Sun, D. (2014). Global synthetic seismograms using a2-D finite-difference method.
Geophysical Journal International , 197(2):1166–1183.[12] Lindsey, C. and Braun, D. (2004). Principles of seismic holography for diagnostics of the shallowsubphotosphere.
The Astrophysical Journal Supplement Series , 155(1):209.[13] Lomas, A. and Curtis, A. (2019). An introduction to Marchenko methods for imaging.
Geophysics ,84(2):F35–F45.[14] Lomas, A. and Curtis, A. (2020). Marchenko methods in a 3-D world.
Geophysical Journal Inter-national , 220(1):296–307.[15] Oristaglio, M. L. (1989). An inverse scattering formula that uses all the data.
Inverse Problems ,5(6):1097–1105.
16] Pereira, R., Ramzy, M., Griscenco, P., Huard, B., Huang, H., Cypriano, L., and Khalil, A. (2019).Internal multiple attenuation for OBN data with overburden/target separation. In
SEG annualmeeting, 2019 , pages 4520–4524. Society of Exploration Geophysicists.[17] Porter, R. P. (1970). Diffraction-limited, scalar image formation with holograms of arbitrary shape.
Journal of the Optical Society of America , 60(8):1051–1059.[18] Porter, R. P. and Devaney, A. J. (1982). Holography and the inverse source problem.
Journal ofthe Optical Society of America , 72(3):327–330.[19] Schorlemmer, D., Wiemer, S., and Wyss, M. (2005). Variations in earthquake-size distributionacross different stress regimes.
Nature , 437(7058):539–542.[20] ˇS´ılen`y, J. (2009). Resolution of non-double-couple mechanisms: Simulation of hypocenter mis-location and velocity structure mismodeling.
Bulletin of the Seismological Society of America ,99(4):2265–2272.[21] Slob, E., Wapenaar, K., Broggini, F., and Snieder, R. (2014). Seismic reflector imaging usinginternal multiples with Marchenko-type equations.
Geophysics , 79(2):S63–S76.[22] Staring, M. and Wapenaar, K. (2019). Interbed demultiple using Marchenko redatuming on 3Dfield data of the Santos Basin. In . Brazilian Geophysical Society, SBGf.[23] Thorbecke, J. and Brackenhoff, J. (2019). OpenSource: Code for geophysical 3D/2D Finite Differ-ence modelling, Marchenko algorithms, 2D/3D x-w migration and utilities.[24] Thorbecke, J., Slob, E., Brackenhoff, J., van der Neut, J., and Wapenaar, K. (2017). Implementationof the Marchenko method.
Geophysics , 82(6):WB29–WB45.[25] Wapenaar, K., Brackenhoff, J., and Thorbecke, J. (2019). Green’s theorem in seismic imaging acrossthe scales.
Solid Earth , 10(2):517–536.[26] Wapenaar, K., Thorbecke, J., and van der Neut, J. (2016). A single-sided homogeneous Green’sfunction representation for holographic imaging, inverse scattering, time-reversal acoustics and in-terferometric Green’s function retrieval.
Geophysical Journal International , 205(1):531–535.[27] Wapenaar, K., Thorbecke, J., van Der Neut, J., Broggini, F., Slob, E., and Snieder, R. (2014).Marchenko imaging.
Geophysics , 79(3):WA39–WA57.[28] Willacy, C., van Dedem, E., Minisini, S., Li, J., Blokland, J. W., Das, I., and Droujinine, A. (2018).Application of full-waveform event location and moment-tensor inversion for Groningen inducedseismicity.
The Leading Edge , 37(2):92–99., 37(2):92–99.