Microseismic wavefield modelling in anisotropic elastic media using integral equation method

In this paper, we present a frequency‐domain volume integral method to model the microseismic wavefield in heterogeneous anisotropic‐elastic media. The elastic wave equation is written as an integral equation of the Lippmann–Schwinger type, and the seismic source is represented as a general moment tensor. The actual medium is split into a background medium and a scattered medium. The background part of the displacement field is computed analytically, but the scattered part requires a numerical solution. The existing matrix‐based implementation of the integral equation is computationally inefficient to model the wavefield in three‐dimensional earth. An integral equation for the particle displacement is, hence, formulated in a matrix‐free manner through the application of the Fourier transform. The biconjugate gradient stabilized method is used to iteratively obtain the solution of this equation. The integral equation method is naturally target oriented, and it is not necessary to fully discretize the model. This is very helpful in the microseismic wavefield computation at receivers in the borehole in many cases; say, for example, we want to focus only on the fluid injection zone in the reservoir–overburden system and not on the whole subsurface region. Additionally, the integral equation system matrix has a low condition number. This provides us flexibility in the selection of the grid size, especially at low frequencies for given wave velocities. Considering all these factors, we apply the numerical scheme to three different models in order of increasing geological complexity. We obtain the elastic displacement fields corresponding to the different types of moment tensor sources, which prove the utility of this method in microseismic. The generated synthetic data are intended to be used in inversion for the microseismic source and model parameters.


INTRODUCTION
Activities such as fluid injection in the subsurface and hydraulic fracturing may lead to induced earthquakes and microseismicity.Microseismic wavefield modelling is a necessary step in characterization of microseismic events.In elastic media, we use the stiffness, density and microseismic source parameters to generate the waveform data at the receiver location.The source in microseismic events cannot be represented by a vectorial point source.Instead, moment tensors have been used to describe the source mechanisms of earthquakes (Aki & Richards, 2002;Jost & Herrmann, 1989).The simulation of seismic wave propagation in threedimensional (3D) elastic media with a general moment tensor source representation was first done by Graves (1996).He employed the staggered-grid finite difference method (Virieux, 1986) to solve the elastodynamic equation in the velocity-stress formulation.The work was carried out under the assumption that the earth model is isotropic.
The lithologies targeted for carrying out the fluid injection/extraction activities are usually stratified rocks like shale and sandstone.Along with the horizontal layering, a vertical fracture system is often observed in these rocks.This can be attributed to the triaxial stress state in the earth's crust.Seismic velocity anisotropy can be caused by factors, such as rock fabric, grain-scale microcracks, rock layering and aligned fractures at all scales, provided that the characteristic dimensions of these features are small relative to the seismic wavelength (Worthington, 2008).Seismic anisotropy can have a significant influence on the recorded wavefield, which ultimately leads to errors in localization and characterization of the microseismic source.Stierle et al. (2016) demonstrated that the retrieval of moment tensor and source mechanism critically depends on anisotropy using laboratory acoustic emission experiments.Therefore, the development of a full waveform modelling tool in general anisotropic media is a necessity for accurate inversion of microseismic data.
Microseismic full waveform modelling can be performed in the time domain or the frequency domain.Shi et al. (2018) and Lei et al. (2021) performed the time-domain seismic full-waveform modelling in anisotropic-elastic media with moment tensor source representation.To obtain a symmetric moment tensor source, the shear-stress increments were evenly distributed on the adjacent shear-stress grid points around the true moment tensor source location.In this research, we apply a volume integral method to model the microseismic wavefield in the frequency domain.The moment tensor source can be introduced at any grid point by using an equivalent body force term (Aki & Richards, 2002).While performing waveform inversion, the possibility to invert only a few frequencies in a sequential manner can save computational cost and also avoid getting trapped to local minima.Furthermore, it is easy to include attenuation effects.
The foundations of integral equations lie in the general principles of scattering theory.In the integral equation method, the actual medium is decomposed into a homogeneous reference or background medium and a contrast medium.Using Green's function for the background medium and numerical values of the physical properties in the contrast medium, we obtain the wavefield in the actual medium.The integral equation approach has been widely applied to solve the nonlinear scattering problem in electromagnetics (Jakobsen & Tveit, 2018).The same methodology can easily be extended to elastic wave scattering problems.Gibson and Ben-Menahem (1991) clearly demonstrated the application of elastic wave scattering theory to study wave propagations in fractured media.The Born approximation (Miles, 1960;Hudson & Heritage, 1982) to a fully elastic scattering obstacle is one of the frequently applied methods to compute the scattered field in elastodynamic wave propagation.It is used to calculate the wavefields of Rayleigh scatterings due to the perturbations to elastic constants considering only a single scattering of elastic waves.However, the iterative solver-based volume integral method accounts for multiple scatterings of elastic waves.The Born approximation also predicts incorrect travel times, whereas the full integral equation method accounts for the travel time difference between the actual medium and the background medium.
The volume integral method has some advantages and some disadvantages over the other modelling methods such as the finite difference method.The finite difference method is faster than the integral equation method and is geometrically flexible, such as it can use curved cells.However, it is only required to discretize the domain where the scattering potential is nonzero when using the volume integral method, whereas the finite difference method requires a full discretization of the model.There is no time stepping involved, and we do not encounter differentiation-related numerical artefacts in the simulations.The Toeplitz structure of the linear system matrix in the integral equation formulation (as a consequence of the translation-invariance of the homogeneous background medium) allows one to use the fast Fourier transform algorithm for efficient implementation of matrixvector or operator-function product.Malovichko et al. (2018) demonstrated that the integral equation system matrix is better conditioned than that of the finite-difference method.Additionally, integration is a smooth process and we can work over multiple well-separated domains (e.g., different reservoirs or a salt body and a reservoir within a sedimentary basin).It is still possible to take into account the interactions between these domains using the scattering-path operator (Jakobsen & Wu, 2018).Malovichko et al. (2017) discussed how the integral equation method naturally applies to the computation of the Fréchet derivatives, and for this reason it is very attractive for the solution of the inverse problem.One of the limitations of working with integral equation is that the full integral equation method is known for generating a dense system matrix.However, the integral equation can in principle be combined with a preconditioner that leads to a sparse linear system.Ying (2015) developed a sparsifying preconditioner for the integral equation representing the scalar wave equation.A preconditioner for elastic integral equation can also be designed by extending the idea of the contraction preconditioner originally developed for the integral equation method of electromagnetic modelling (Zhdanov & Fang, 1997).But these implementations are beyond the scope of the present study.
Within the framework of seismic wavefield modelling, a limited number of works have been concluded using the volume integral method, and most of those focused on acoustic waves only (Alles & van Dongen, 2011;Malovichko et al., 2017Malovichko et al., , 2018;;Jakobsen & Ursin, 2015;Xiang et al., 2021).An acoustic or elastic wave equation can be written as a corresponding integral equation of the Lippmann-Schwinger type, and different numerical techniques can be employed to solve this equation.Jakobsen and Ursin (2015) used the transition operator approach for seismic forward modelling in the acoustic approximation.Malovichko et al. (2018) presented frequency-domain acoustic 3D modelling using the integral equations method and proved that the method can accurately calculate the pressure field for models with sharp material boundaries.Tong and Chew (2009) developed a multilevel fast multipole algorithm as an accelerator for integral equation solvers for elastic wave scattering problems.Touhei (2011) proposed a fast method for the volume integral equation for elastic wave propagation in a half-space.The fast transform was implemented by decomposing the kernel of the transform into the ordinary Fourier and Laplace transforms.
Previous applications of the volume integral equation method to elastic wavefield modelling have ignored the effects of seismic anisotropy and have not presented synthetic seismograms for realistic seismic models.In the volume integral method, the heterogeneity and anisotropy parameters can be assigned directly; however, attaining the solution of the resulting integral equation is a challenging task.Jakobsen et al. (2020) discussed how to solve nonlinear inverse scattering problems in general anisotropic-elastic media with a variable density as well as stiffness tensor, but the numerical examples were restricted to media with constant density.They used a matrix-based solver for the two coupled integral equations to perform elastic wavefield modelling.Their methodology was inefficient in terms of computational cost and memory.
Our goal is to carry out microseismic wavefield modelling using the volume integral method.To simulate the wavefield in reservoir-scale 3D models, we focus on the development of a computationally efficient numerical solver for the integral equation.In the second section , we show how the integral equation can be formulated to incorporate the moment tensor source into the heterogeneous anisotropic-elastic media.We utilize the Lippmann-Schwinger equation for elastic waves (Jakobsen et al., 2020) and implement the moment tensor source (Stein & Wysession, 2003;Madariaga, 2015) in the equation using a body force term.In the third section , we demonstrate how to efficiently solve the integral equation in the particle displacement.We emphasize the Fourier transform-based iterative solver to perform the modelling work in a matrix-free manner.Symmetries in the off-diagonal components of the moment tensor, the stiffness tensor and the strain tensor allow for the use of the Voigt notation during the implementation.In the integral equation method, components of the strain field are also needed to obtain the particle displacement components.The finite difference method is used to numerically compute the strain.In the fourth section , we apply the method to a homogeneous 3D earth model.We observe the response of the wavefields to different types of moment tensor sources.We next test the accuracy of the integral equation method in computing the microseismic wavefield in the anisotropic-elastic media.This is accomplished by comparing the integral equation result with the finite difference time domain result for the same model parameters.Finally, we examine the microseismic wavefield in laterally inhomogeneous media to explore the applicability of the integral equation method for heterogeneous models.

INTEGRAL EQUATION FORMULATION
Waves generated in heterogeneous media can be reflected, refracted, transmitted and/or diffracted.The elastodynamic wave equation takes all these phenomena into account.The components of the particle displacement vector   () at point  due to a force source density with components   () in an anisotropic elastic medium with stiffness tensor component   () and mass density () satisfies the elastodynamic wave equation where  is the angular frequency.The value of each of the indices can be 1, 2 and 3, indicating the X, Y, and Z components, respectively.We assume that   () is a plane-wave solution to Equation (1).In the integral equation approach, the complex interactions of the wave phenomena are termed scatterings.In this method, we decompose the stiffness tensor components   () and the mass density () as where  (0)  and  (0) are the stiffness coefficients and mass density, respectively, of an arbitrary homogeneous reference/background medium and Δ  () and Δ() are the corresponding contrast terms.
From Equations ( 1 (4) The second and third terms on the right-hand side of Equation (4) represent the so-called contrast sources of the moment tensor and force vector type, respectively.By treating the contrast sources similar to ordinary sources and using Green's function representation (Jakobsen et al., 2020), we obtain the following integral equation: The integral equation in the particle displacement (Equation 5) demonstrates the linear superposition of the background and scattered fields.If the elastodynamic Green's functions  (0)   for the background model are known, it is only required to discretize the contrast model.The integral part has been computed over a scattered domain Ω, and  ′ is a position vector within this region.The superscript (0) denotes a physical quantity in the homogeneous background medium.The wavefield in the homogeneous background model due to the vectorial source with components   () is obtained by the following convolution integral: where Here, ( − ) is the three-dimensional (3D) Dirac-delta function and   is the Kronecker delta defined by   = 1, when  = , and   = 0 when  ≠ .In elasticity theory, Green's function gives the displacement generated by a point force in a certain direction (Snieder, 2002).Here, Green's tensor  (0)   (, ) is the displacement at location  in the  direction due to the th component of a unit force in the  direction at location .The details on Green's function and its spatial derivative can be found in the Appendix.
As the Green's function is inversely proportional to the distance between the source and the receiver, it is assumed that the wavefield reduces to zero at infinity.By performing a partial integration, and using the well-known symmetries of the elastic stiffness tensor, Equation ( 5) can be rewritten as (Červený, 2001) where the strain tensor is the symmetric gradient of the displacement vector and is given by and  (0)  (,  ′ ) is the first-order spatial derivative of Green's tensor for the background medium.It represents the th component of the elastic displacement due to the  component of the contrast stress tensor and is expressed as Equation ( 8) is known as the Lippmann-Schwinger equation for elastic wave scattering.In the derivation of this integral equation, we have made use of the fact that Green's function for the background medium only depends on the difference between the source and receiver coordinates since a homogeneous background medium is translation invariant.We also used the fact that any tensor can be decomposed into symmetric and anti-symmetric parts and that the contraction of a symmetric tensor with an anti-symmetric tensor is identical to zero (Pujol, 2003).
Since the third-rank tensor field  (0)  () by construction is symmetric with respect to an interchange of the last two indices, we can express Equation ( 8) in an abbreviated notation as Here, the small and capital indices take values from 1 to 3 and from 1 to 6, respectively, and Einstein's summation convention applies (see Auld, 1990).Finally, we note that the above equation can be expressed more compactly in the hybrid tensor-operator notation as In the above expression, we have kept the tensor indices and the abbreviated indices when expressing the integral equation in an operator form.The use of hybrid tensor-operator notation is convenient when discussing the algorithms for solving the integral equation iteratively.To simplify the presentation, we use the same notation for the integral operator   and its kernel function   (,  ′ ) defined in Equation ( 11).

Moment tensor source representation
A moment tensor can be used to represent the seismic source in the microseismic events related to the slip on a fault plane.The mathematical models of this type of source can be given by a 3 × 3 tensor   .Due to the symmetry of the offdiagonal elements, it can be written in the Voigt notation as a 6 × 1 vector   .The expression of the moment tensor source can be understood by analysing the kinematics behind the distributions of the transformation strain.The transformation strain describes an alteration of the stressfree configuration of a solid (Rice, 1998).Depending on the focal mechanism for microseismic events, the moment tensor source can be categorized as a single dipole, doublecouple (DC), explosion-type source, etc. (Stein & Wysession, 2003).Planar shear faulting produces a pure DC mechanism in isotropic media, but generally a non-double-couple (non-DC) mechanism in anisotropic media (Vavryčuk, 2005).The non-DC mechanism can be decomposed into an isotropic component and a compensated linear vector dipole component (Vavryčuk, 2015).The moment tensor component of DC microseismic sources in anisotropic formations have been eloquently derived by Grechka (2020).
The equivalent body force density (Aki & Richards, 2002) for the moment tensor source is where   is a moment tensor component,  represents the positioning of the microseismic event and  0 may be referred to as the rupture time or the timing of the microseismic event.Using Equation ( 13) in Equation ( 6), we obtain the background displacement field, It is a common practice to assume that the moment tensor component   is constant during the entire period of seismic activity.We can simulate the seismic radiation from a moment tensor source by using the spatial derivative of Green's function.We exploit the property of the Dirac-delta function; that is, if we integrate the product of a function and the derivative of the Dirac-delta function over the spatial domain, the result is the spatial derivative of that function computed at the source location.The particle displacement field at position  due to a moment tensor source at position  is now given by In an abbreviated notation by

EFFICIENT IMPLEMENTATION OF THE INTEGRAL EQUATION METHOD
We work in the frequency domain and want to demonstrate wave propagation in heterogeneous anisotropic media using the full integral equation solution.To develop a computationally efficient numerical scheme, we employ the following methodology.
In the integral equation formulation, there is a simultaneous computation of both the strain field and the displacement field.If we solve the two coupled integral equations for the particle displacement and strain (Jakobsen et al., 2020), the numerical computation is slow.This is because it requires the calculation of the second-order spatial derivative of the background Green's function.We compute the strain field using the finite difference operator on the elastic displacement vectors to calculate its derivative with respect to the position vector.Let the grid spacing in each of the directions be ℎ and the displacement components in the X-, Y-and Z-directions be  1 ,  2 and  3 , respectively.The strain field components at a grid node , ,  are given as 5 (, , ) = 1 2ℎ [  1 (, ,  + 1) −  1 (, ,  − 1) To compute the strain field near the model boundaries, depending on the position where we are computing the derivative, we use either the forward or the backward difference.We fixed a model having 64 grid points in each direction and used the same model parameters to calculate the strain field by using the two coupled integral equations and by this strategy.We found that the results are quite similar.The highest relative difference was 0.1.The major differences are observed at the source location and at a few other points in the near field.This is related to the singularity of Green's function, the frequency used in the simulation and the size of the grid.But in the far field, which is the main region of interest, we observe no striking differences.The difference in the computational time between the two approaches was 802 s.Though we do not intend to perform any further quantitative study here, we can say that this gap increases drastically as the size of the model increases.

Fast Fourier transform accelerated Krylov subspace method
Modelling seismic data through the implementation of a matrix-based method involves high computational costs and memory consumption.We utilize an iterative solver based on the fast Fourier transform (FFT; Stein & Wysession, 2003) to perform efficient seismic modelling, that is, in a matrix-free manner.The FFT iterative scheme is based on a recurrence relation and does not lead to the resolution of a linear system (Monchiet & Bonnet, 2012).The computational cost decreases from the order of  3 for the matrix-based direct solver method to the order of log for the matrix-free iterative method, where  is the number of grid points.The computational memory scales down from the order of  2 to the order of .
Equation ( 12) can be rearranged as A fast formulation of the displacement field can be obtained using the background Green's function and the spatial derivative of the background Green's function in the Fourier transformed coordinates (k-space): where  and  −1 are the three-dimensional fast Fourier and inverse fast Fourier transforms, respectively,  is a position vector in the real space and "•" denotes the elementwise multiplication operation.Although the formulation of the particle displacement has been obtained using the continuous form of FFT over the spatial domain, the discrete form of FFT is needed to perform the computational operation.The background medium should be homogeneous to apply the Fourier transform.Equation ( 21) can be written in the form  =  (0) .Here,  can be regarded as an operator acting on the elastic displacement component, and  (0) is the displacement field component in the homogeneous background medium.Instead of using a direct solver, we solve the system of equations in an iterative way using the conjugate-gradient (CG) method (Shewchuk, 1994;Nocedal & Wright, 2000); however, the CG is applicable to symmetric positive definite systems only.When the coefficient matrix is nonsymmetric and nonsingular, iterative methods such as the biconjugate gradient (BiCG), biconjugate gradient stabilized (BiCGSTAB) and conjugate gradient squared (CGS) are useful (Barrett et al., 1994).
T A B L E 1 Runtime for the matrix-free integral equation method at a single frequency (10 Hz) and using the homogeneous VTI model.and the other moment tensor components are zero.

Number of grid points in each direction
We apply the BiCGSTAB method for the numerical solution, as this method has a faster and smoother convergence than the other variants of the CG method (Van der Vorst, 1992).The computational cost per iteration for BiCGSTAB is similar to BiCG and CGS, but the method does not require the computation of the transpose of the coefficient matrix, unlike BiCG.Additionally, there is no issue of error in the solution like the ones observed with CGS.In CGS, local corrections to the current solution may be so large that cancellation effects occur.This may lead to a less accurate solution than suggested by the updated residual (see Van der Vorst, 1992).We do not want to use the generalized minimal residual method ( Saad, 2003) because of the increasing storage requirements per iteration step, despite the need to compute only matrix-vector products with the coefficient matrix.
The relative residual error (e) at each iteration step can be given as where ‖ ⋅ ‖ is Euclidean norm.We test the convergence of our scheme by plotting the relative residual error with the iteration numbers.Overall, this solver can be applied to a large enough model for which the matrix-based solver is unable to execute on (we consider a 16 GB memory processor for simulations).The approximate details on the runtime of simulation per frequency for different model sizes are given in Table 1.

Absorbing boundary condition
One of the challenges in simulating the seismic wavefield in any model is to avoid reflections at the boundaries.In the finite difference method, an artificial absorbing boundary can be implemented by either using the damping boundary condition or perfectly matched layers (Lei et al., 2021).Alles and van Dongen (2011) developed perfectly matched layers (PML) for frequency-domain integral equation acoustic scattering problems.A strong attenuation of scatter pressure fields is achieved in layers with a thickness of less than a wavelength by using a plane-wave function that has different solutions within and outside the PML domain.Osnabrugge et al. ( 2016) used a polynomial-based expression for the wavefield inside the absorbing layer for a similar type of scattering problem and found an expression for the scattering potential of the layer in terms of the solution for the polynomial.
The idea was to design a layer that has zero reflectivity for a normal incidence.
Green's function approach to calculate the seismic wavefield incorporates the exponential decay of the wavefield with the distance.Most of the absorbing boundary conditions given in the differential form can be transformed into an exponential decay in the integral equation.Strictly speaking, it is difficult to completely avoid any reflections from the boundaries; however, the amplitude of the reflected wave can be significantly minimized.The contrast between the actual medium and the background medium should be zero to have no reflections at all.In other words, the actual medium wavefield is the same as the background field.For the integral equation in the particle displacement, the absorbing boundary condition is introduced by formulating the scattered field equation such that it gradually vanishes as it approaches the boundaries.
Let Ω be the domain of the scattered region and Ω  be the domain of the boundary layer region for a model and   be a position within Ω  .Let   be the position vector for a point obtained by extending the position vector   to the external boundary of the model.Then, We can write the integral equation for the displacement as Then, the absorbing boundary condition can be applied through the following numerical expression: ] . (25) Here, ℎ is the grid spacing and   is the number of grid points in the  direction.
Choosing a background medium that is closer to the actual medium can sufficiently reduce the reflections from the boundaries.But this is restricted to homogeneous or slightly heterogeneous models.It is difficult to understand the effects of the boundary conditions in the frequency domain.Therefore, we will demonstrate its impact later in the numerical result section when we generate the time-domain data using the inverse Fourier transform.

NUMERICAL RESULTS AND DISCUSSION
We perform all the numerical simulations for a threedimensional (3D) Cartesian mesh.The grid is uniformly spaced in each direction.The volume of each cell is the product of the grid spacing in the X-, Y-and Z-directions.Each node characterizes a unique stiffness matrix and a density value.The physical properties do not vary laterally in the Y-direction.We generally show snapshots of the models on the XZ plane.The microseismic source is positioned at the centre of the model unless otherwise stated.To obtain an accurate scattered field, the grid size is chosen as one-fourth or smaller than the smallest wavelength of the S-waves in the actual medium.

Model 1: Homogeneous vertical axis of symmetry model
The homogeneous model represents the transversely isotropic model with a vertical axis of symmetry (VTI).The density value is 2500 kg/m 3 (constant within the model), and the values for the stiffness coefficients for the VTI medium (Jakobsen & Johansen, 2000 : table 3, sample at a depth of 3492 m) in units of GPa are 34.0 10.6 6.9 0 0 0 34.0 6.9 0 0 0 26.5 0 0 0 10.4 0 0 10.4 0 11.7 This VTI medium is the actual medium in the integral equation approach, and it is decomposed into an isotropic background medium and a contrast medium.For the isotropic background medium, the stiffness coefficients  11 and  44 are 26.5 and 10.4 GPa, respectively.The numerical values of the density at every grid point are the same in the background medium and the actual medium.The dimension of the model is 640 m in each direction, and the grid size is 10 m.

Displacement wavefields due to a moment tensor source
We want to demonstrate the microseismic wavefield in the isotropic background and actual VTI media.The absolute value of the elastic wavefield normalized with respect to the background medium wavefield on a decibel scale is given as The wavefront for a vector point source in an isotropic media is spherical.However, for the moment tensor source, the wavefront is lobe shaped.In the time-domain snapshot (Figure 1), the wavefront has a different polarity in opposite directions, which means that the pattern of compression and rarefaction advances in a different fashion around the source.Figures 1a and 1b show the X and Z components, respectively, of the time-domain wavefield.Figures 2-5 show the displacement field generated by a moment tensor source on the XZ-plane passing through the source at a particular frequency (20 Hz in this simulation).Figures 2 and 4 represent frequency-domain wavefields due to a single-dipole source and a double-dipole source, respectively, in isotropic media.Figures 3 and 5 represent frequency-domain wavefield due to a single-dipole source and a double-dipole source, respectively, in VTI media.Figure 6 shows the strain field generated by an explosion type of moment tensor source in the isotropic background medium.We label the strain field components in the Voigt notation.In Figure 7, we present the elastic wavefield in the VTI medium due to a double-couple source.The numerical expression for the source term is obtained following Grechka (2020).Figures 7a, 7b and 7c show the X, Y and Z components, respectively, of this wavefield.The convergence of the numerical scheme for this simulation is shown in Figure 8.

Comparison with the finite difference time domain method
To verify our methodology, we compare our results to the results obtained by using the fourth-order finite difference scheme in the time domain.The order of the differential operator, the discretization technique and the fashion in which the source and its adjacent nodes are treated can significantly affect the result of the finite difference method.The higher the order of the differential operator is the more accurate the wavefields obtained.However, the computational cost increases as more floating-point operations are needed (Sei & Symes, 1995).The finite difference solution to the elastic wave equation with the moment tensor source implementation can be obtained using the velocity-stress formulation on a staggered grid (Shi et al., 2018).The particle displacement can be obtained from the particle velocities by performing the integration over time.The elastic displacement field at the receiver position due to a moment tensor source, obtained by the integral equation method, is in the frequency domain.It is based on the evaluation of Green's tensor and its derivative.We need time domain data to compare this result with the wavefield of the displacement components obtained by the finite difference time domain (FDTD) method for the same model setup (Table 2).The source-receivers setup is shown in Figure 9.We set the values for the total time of F I G U R E 1 1 Comparison between the integral equation and the finite difference methods for the model parameters given in Table 2: (a) Synthetic seismograms (X-component of the particle displacement) generated using the integral equation method, (b) synthetic seismograms (X-component of the particle displacement) generated using the finite difference method, (c) difference between the seismograms obtained using two methods and (d) elastic displacement components at a common receiver.
T A B L E 3 Elastic stiffnesses (in units of GPa) for Model 2 -heterogeneous VTI model.

Stiffness coefficients
Top simulation and the time interval equal to the respective values used in the FDTD method.The chosen time step satisfies the Courant-Friedrichs-Lewy condition.First, we take a loop over positive and negative frequencies up to the Nyquist frequency and compute the wavefield corresponding to each frequency.Then, we perform the inverse Fourier transform to generate time-domain particle displacement data.We can see the effects of reflections from boundaries in the synthetic seismogram (Figure 10).Therefore, we need to apply the boundary condition introduced earlier in the implementation section to mitigate this effect (Figure 11a).Figures 11a and 11b show the wiggle trace generated using the integral equation and the finite difference methods, respectively.Figure 11c shows the difference between those two wiggle traces.The waveforms of the elastic displacement field obtained using the integral equation method and the FDTD method at a common receiver match well (Figure 11d).

Model 2: Heterogeneous vertical axis of symmetry model
The 3D model (Figure 12) has three layers.An isotropic layer is sandwiched between two VTI layers, the parameters of which are the same as the background medium in the integral equation.This model represents an underburden-reservoiroverburden rock system.The dimension of the model is 320 m in each direction, and the grid size is 2.5 m.The thickness of the top and bottom layers is 108 m.The middle layer is 104 m thick.The density value is 2500 kg/m 3 and is constant within the model.The values of the stiffness coefficients (Jakobsen & Johansen, 2000 : table 3) in the different layers are given in Table 3.

Displacement wavefields due to a moment tensor source
In Figure 13, we see the displacement field generated by an explosion-type moment tensor source on the XZ-plane passing through the source at a frequency of 50 Hz.We select the high frequency to demonstrate the scattering of waves from the layer interfaces.The strain field components  1 - 6 are shown in Figure 14.A significant amount of energy remains confined to the layer in which the source is present, that is, in the middle layer.The numerical convergence of the scheme is presented in Figure 15.

Comparison with the finite difference time domain method
We perform a comparison between the integral equation method and the finite difference time domain method to test the accuracy of the numerical scheme in heterogeneous media.The model parameters used in the comparison are given in Table 4.The source-receivers setup is shown in Figure 16.We obtain a satisfactory match between the synthetic seismograms obtained using these two methods (Figure 17).Figures 17a and 17b show the wiggle trace generated using the integral equation and the finite difference methods, respectively.Figure 17c shows the difference between those two wiggle traces.The waveforms of the elastic displacement field obtained using the integral equation method and the FDTD method at a common receiver also look identical (Figure 17d).The amplitude and phase information from the full waveform modelling can be used to demarcate the reflection and transmission events in the layered model.For this model, the strong amplitudes recorded at the receiver positions represent either the arrival of the direct waves or the arrival of the waves transmitted through layer interfaces, and the weak amplitudes are due to the arrival of the waves reflected from layer interfaces.For example, in receiver 10 (Figure 17a), the elastic displacement field at time 100 ms is due to the wave reflected at the interface between isotropic and upper VTI layers.The arrival prior to that is a direct wave.We can do a similar

Model 3: Hess model
We also apply our numerical scheme to the resampled Hess/SEG salt model (Figure 18).The grid size is 10 m.We use the two-dimensional dataset and extend the model in the Y-direction to 1 km to make it three dimensional.We want to observe the effect of the lateral variations in the physical properties of the medium on the wavefield.The density values vary within the medium.The model contains a fault zone.Along with the sedimentary layers, a salt body is also present.Hence, the geology of this model is much more complex than that of Models 1 and 2.

Wavefield of displacement components
We observe the wavefield on the XZ-plane passing through the microseismic source.The strain field components in the Voigt notation,  1 to  6 are shown in Figure 19.The seismic source is an explosion-type moment tensor.The X, Y and Z components of the elastic displacement written here as   ,   and   , respectively, are shown in Figure 20.The presence of the salt body and the fault can be inferred by observing the strain field components.The importance of the strain field calculation is that the strain components are needed during the inversion for the model parameters.
Although a large number of iterations are needed, the numerical scheme converges for this complex model (Figure 21).Furthermore, we tested our scheme by setting the contrast between the stiffness coefficient values at every grid point in the actual medium and the homogeneous background medium to be up to 10% and found that it still converges.4: (a) Synthetic seismograms (X-component of the particle displacement) generated using the integral equation method, (b) synthetic seismograms (X-component of the particle displacement) generated using the finite difference method, (c) difference between the seismograms obtained using two methods and (d) elastic displacement components at a common receiver.

CONCLUDING REMARKS
We have developed a scattering approach to frequency domain modelling of seismic waveform data in general anisotropic elastic media with variable mass density and elastic stiffness tensor.Our focus was on the modelling of microseismic waveform data due to a moment tensor source; however, our scattering approach can also be used for frequency domain full waveform inversion; that is, to construct the anisotropic elastic background model in a microseismic imaging experiment.Our scattering approach is based on the transformation of the elastodynamic wave equation into an equivalent volume integral equation of the Lippmann-Schwinger type that contains the boundary conditions.This volume integral equation involves decomposition of the elastic stiffness and mass density fields into a spatially invariant part related to an arbi-trary homogeneous reference medium and a spatially variable contrast model.It is only required to discretize that part of the model where the elastic stiffness and mass density fields differ from that of the isotropic elastic reference medium, making this method naturally target-oriented.
To solve this integral equation, we did not make use of any approximations that could limit the applicability of our scattering approach to media with small contrast.Instead, we solve the integral equation iteratively using the fast Fourier transform (FFT) accelerated Krylov subspace method; that is, a combination of FFT and the biconjugate gradient stabilized method.Previous attempts to solve this integral equation for anisotropic elastic media have typically assumed a fixed mass density and used a matrix representation of the relevant integral operators; where the computational cost and memory requirements scaled like  3 and  2 , respectively, where   is the number of grid points.In our FFT-accelerated iterative solution, however, the computation cost and memory requirements have been reduced to the order of  × log() and , respectively.Using the FFT in the matrix-vector multiplication comes with restrictions that the cells cannot be curved and the grid must be regular.In the context of full-waveform inversion, however, it is common to use a regular Cartesian grid.
Although we have now significantly reduced the computational cost of the integral equation method, it may still not be ready for large-scale three-dimensional (3D) applications.In one of our numerical examples of frequency wavefield simulation in a large 3D anisotropic elastic model with variable density, a relatively large number of iterations was required in order to reduce the relative residual error of the linear system down to 1e −6 .We think that the number of iterations can be significantly reduced if we employ a preconditioner to our linear system.Efficient preconditioning techniques like convergent Born series (Osnabrugge et al., 2016), renormalization and homotopy analysis exist for the integral equation method, but the most promising ones are limited to acoustic approximation or have been developed for electromagnetic scattering problems.Therefore, our plans for future research include a generalization of some of these preconditioners from the acoustic or electromagnetic domains to the anisotropic elastic case.Another suggestion for future work is to develop a matrix-free variant of the scattering domain decomposition method introduced in the transition operator paper of Jakobsen et al. (2020).A matrix-free domain decomposition strategy for the integral equation method was recently presented for the electromagnetic case (Saputera et al., 2023), suggesting that it is possible to proceed towards the anisotropic elastic problem.The integral equation method we have developed can easily be implemented on graphical processing units, which can lead to a significant improvement in efficiency.Also, the integral equation method allows us to treat different frequencies and sources in an embarrassingly parallel manner.Therefore, the CPU times observed in this study may be significantly reduced if we optimize our implementation.
One can in principle combine this integral equation method with ray theory and focus the modelling (or inversion) on a relatively small target.We have previously developed a method for target-oriented waveform inversion based on a combination of the integral equation method and ray-based methods (Huang et al., 2019) in acoustic approximation.The anisotropic elastic integral equation developed here may be particularly attractive for applications to time-lapse studies and near-surface characterization, in addition to microseismic imaging.Also, it may be attractive to studies of seismic anisotropy, since it can deal with arbitrary anisotropic media, although numerical experiments in this work involved transversely isotropic media with vertical axis of symmetry only.
Computational time (s) Nx = Ny = Nz = 16 108.32Nx = Ny = Nz = 32 612.16Nx = Ny = Nz = 64 1375.49Nx = Ny = Nz = 128 2181.70F I G U R E 1 Snapshot of the (a) X-component and (b) Z-component of the elastic displacement field in the isotropic background medium for the homogeneous model at time t = 0.12 s, when  11 =  33 = 1 √ 2 and the other moment tensor components are zero.F I G U R E 2 Displacement field at 20 Hz in the isotropic background medium for the homogeneous model when  33 = 1 Displacement field at 20 Hz in the homogeneous VTI medium when  33 = 1 √ 2

F
Displacement field at 20 Hz in the isotropic background medium for the homogeneous model when  11 =  33 = 1 Displacement field at 20 Hz in the homogeneous VTI medium when  11 =  33 = 1 √ 2 and the other moment tensor components are zero.F I G U R E 6 Strain field at 10 Hz in the homogeneous isotropic background medium due to an explosion-type moment tensor source.F I G U R E 7 3D view of the (a) X-component, (b) Y-component and (c) Z-component of elastic displacement field at 10 Hz in the homogeneous VTI model due to a double-couple source.

F
Plot of the relative residual error versus iteration number for the BiCGSTAB solver in the case of the homogeneous VTI model.
Abbreviations: FDTD, finite difference time domain; VTI, transversely isotropic media with vertical axis of symmetry.

F
Presentation for the source (*) and receivers (Δ) mutual placement in a 3D homogeneous VTI model.

F
Synthetic seismograms in the absence of absorbing boundaries.

F
Displacement field at 50 Hz for the heterogeneous VTI model due to an explosion-type moment tensor source.F I G U R E 1 4 Strain field at 50 Hz for the heterogeneous VTI model due to an explosion-type moment tensor source.
Abbreviation: VTI, transversely isotropic media with vertical axis of symmetry.

F
Comparison between the integral equation and the finite difference methods for the model parameters given in Table

F
Model 3: Hess model.F I G U R E 1 9 Strain field due to an explosion-type moment tensor source at 10 Hz frequency for the Hess model.

F
Displacement field due to an explosion-type moment tensor source at 10 Hz frequency for the Hess model.F I G U R E 2 1 Plot of the relative residual error versus iteration number for the BiCGSTAB solver in the case of the variable density Hess model.