Analytic and numerical solutions to the seismic wave equation in continuous media
aa r X i v : . [ phy s i c s . g e o - ph ] N ov rspa.royalsocietypublishing.org Research
Article submitted to journal
Subject Areas:
Geophysics, Applied Mathematics
Keywords:
Seismology, Wave Propagation,Spectral Method
Author for correspondence:
Stephen J Walterse-mail: [email protected]
Analytic and numericalsolutions to the seismic waveequation in continuous media
S. J. Walters , L. K. Forbes and A. M.Reading Mathematics and Physics, University of Tasmania
This paper presents two approaches to mathematicalmodelling of a synthetic seismic pulse, and acomparison between them. First, a new analyticalmodel is developed in two-dimensional Cartesiancoordinates. Combined with an initial condition ofsufficient symmetry, this provides a valuable check forthe validity of the numerical method that follows. Aparticular initial condition is found which allows fora new closed-form solution. A numerical scheme isthen presented which combines a spectral (Fourier)representation for displacement components andwave-speed parameters, a fourth order Runge-Kuttaintegration method, and an absorbing boundary layer.The resulting large system of differential equations issolved in parallel on suitable enhanced performancedesktop hardware in a new software implementation.This provides an alternative approach to forwardmodelling of waves within isotropic media whichis efficient, and tailored to rapid and flexibledevelopments in modelling seismic structure, forexample, shallow depth environmental applications.Visual comparisons of the analytic solution and thenumerical scheme are presented. © The Authors. Published by the Royal Society under the terms of theCreative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author andsource are credited. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A ..........................................................
1. Introduction
Wave propagation in elastic media, and the accurate capture of the resulting time series of velocityvariation with time for a given point, are of significant value in the subsequent determinationof hidden subterranean structure. Through these processes, seismic waves potentially yielduseful information on material properties and dimensions of the structures, and other bodies,comprising the buried parts of the Earth [1,2].Seismic wave propagation has attracted interest for many decades. An early focus was onthe propagation of surface waves, being the highest amplitude and therefore most destructivedisturbances resulting from shallow earthquakes. Lamb investigated the development of an initialimpulse to the surface of an elastic half space (often referred to as Lamb’s problem [3]). Duringthe twentieth century, this was further investigated, notably by Garvin [4]. Kausel [5] recast theproblem and solution into a relatively simple form. Analytic solutions to the seismic equationsin an infinite space have also been considered, although receiving less attention than Lamb’sproblem. Solutions have been found for simple initial conditions such as a step-function or a Diracdelta function [6,7]. Semi-analytic solutions have also been found in terms of Green’s functions [8].Using a numerical approach (see [9] for descriptions of the six most common methods),a number of advanced wave propagation and waveform modelling codes exist that weredeveloped for complex Earth models and global seismology applications [10–12]. However, asseismic methods gain increasing usage in environmental and other near-surface applications(e.g. [13–15]), the research imperative has emerged for an efficient waveform propagation codetailored to simple, readily adaptable, structures. Such code would run at high resolution onmodest hardware and hence be of wide, practical value where the application calls for an agileimplementation. Questions of interest include the optimisation and placement of seismic sensorsto resolve information about an underlying ice sheet and its interface with the bedrock beneath.In this contribution, following a summary of relevant theory, we present a new analyticsolution for the time-dependent propagation of seismic signals in an infinite 2D space. Use ofan infinite space allows freedom in the choice of initial conditions to suit the application ofinterest. A further advantage is the relative simplicity of solutions for the infinite space whichfacilitates testing of numerical codes for propagation within a medium. Choice of particular initialconditions leads to particular solutions, and we found one case which has a closed form solution.We also present a novel numerical implementation for media with smoothly varying wave-speedparameters. This numerical scheme is validated against the analytical solution, and shows thewave propagation and synthetic seismic waveforms for simple structures within a continuousisotropic medium.
2. Theoretical Background
In this section we briefly review the equations of motion for a linear, isotropic, elastic system. Webegin with Cauchy’s momentum equation (that is, Newton’s second law) in index notation: ρ ∂ u i ∂t = ∂ j τ ij + f i , (2.1)where ρ is the density of the medium, u i are the components of the displacement from equilibriumat each point, ∂ j ≡ ∂/∂x j is the derivative operator, τ ij is the stress tensor and f i are thecomponents of any forcing terms (e.g. gravity). In this paper we ignore such forcing terms. Also,we work exclusively in Cartesian coordinates, so all coordinates may be written as lower indices.We also adopt the Einstein summation convention, whereby a repeated index in a single termimplies summation over that index. Assuming a linear, isotropic stress-strain relation (Hooke’slaw), the stress tensor may be defined as τ ij = λδ ij e kk + 2 µe ij , (2.2) r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... where λ and µ are the Lamé parameters, δ ij is the Kronecker delta symbol, and Cauchy’s straintensor e ij is defined as e ij = ( ∂ i u j + ∂ j u i ) . (2.3)The generalized Hooke’s law (2.2) is given in [1], p.34, and we observe that the Lamé terms λ and µ may be functions of position within the medium. Substituting (2.3) into (2.2), we can write thestress tensor as τ ij = λδ ij ∂ k u k + µ ( ∂ i u j + ∂ j u i ) . (2.4)We may now substitute (2.4) into (2.1) to write out the displacement wave equation: ρ ∂ u i ∂t = ∂ j [ λδ ij ∂ k u k + µ ( ∂ i u j + ∂ j u i )] . In compact notation, this may be written as ρu i,tt = λ ,i u k,k + µ ,j ( u j,i + u i,j ) + ( λ + µ ) u j,ij + µu i,jj . (2.5)Here we have used the comma notation for derivatives, in which indices following the commaindicate differentiation with respect to the indicated coordinate (e.g. u i,jj ≡ ∂ u i /∂x j ) .In this paper, we will consider two-dimensional systems. Designating the two components ofdisplacement with upper indices for readability, u X and u Y , we write (2.5) as ρ ∂ u X ∂t = ∂λ∂x (cid:0) ∂u X ∂x + ∂u Y ∂y (cid:1) + 2 ∂µ∂x ∂u X ∂x + ∂µ∂y (cid:0) ∂u Y ∂x + ∂u X ∂y (cid:1) + ( λ + 2 µ ) ∂ u X ∂x + ( λ + µ ) ∂ u Y ∂x∂y + µ ∂ u X ∂y ρ ∂ u Y ∂t = ∂λ∂y (cid:0) ∂u X ∂x + ∂u Y ∂y (cid:1) + 2 ∂µ∂y ∂u Y ∂y + ∂µ∂x (cid:0) ∂u Y ∂x + ∂u X ∂y (cid:1) + ( λ + 2 µ ) ∂ u Y ∂y + ( λ + µ ) ∂ u X ∂x∂y + µ ∂ u Y ∂x . (2.6)
3. Analytical Model
In order to check the accuracy of any numerical model, and in particular the method which wepresent in Section 4, it is valuable to derive an exact non-trivial solution for the system. In order tofind an exact solution, we may consider a particular case of (2.6) for the parameters and boundaryconditions. As we are trying to find an exactly solvable system, we consider the case where theLamé parameters λ and µ and the density ρ are constant throughout the space. Equations 2.6 thensimplify to ρ ∂ u X ∂t = ( λ + 2 µ ) ∂ u X ∂x + ( λ + µ ) ∂ u Y ∂x∂y + µ ∂ u X ∂y ρ ∂ u Y ∂t = ( λ + 2 µ ) ∂ u Y ∂y + ( λ + µ ) ∂ u X ∂x∂y + µ ∂ u Y ∂x . (3.1)For the boundary conditions, we assume an infinite medium in which any initial disturbancenever reaches the boundary at infinity, so that u X , u Y → as x, y → ±∞ . Additionally, we specifyan initial disturbance to the displacement: u X ( x, y,
0) = f ( x, y ) and u Y ( x, y,
0) = g ( x, y ) for somespecified initial displacement functions f and g . The medium is set to be initially at rest: ∂u X /∂t = ∂u Y /∂t = 0 at t = 0 . r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... We now seek a solution using Fourier transforms. The transform and inverse transform for thedisplacement are ˆ U X ( ω , ω ; t ) = 12 π Z ∞−∞ Z ∞−∞ u X ( x, y, t ) e − iω x − iω y d x d yu X ( x, y, t ) = 12 π Z ∞−∞ Z ∞−∞ ˆ U X ( ω , ω ; t ) e iω x + iω y d ω d ω , and similar forms for u Y . The known initial conditions f and g also have Fourier Transforms ˆ F ( ω , ω ) = 12 π Z ∞−∞ Z ∞−∞ f ( x, y ) e − iω x − iω y d x d yf ( x, y ) = 12 π Z ∞−∞ Z ∞−∞ ˆ F ( ω , ω ) e iω x + iω y d ω d ω , and similar forms for g .Application of these Fourier Transforms to the partial differential equations (3.1) and initialconditions leads to a system of two ordinary differential equations in the Fourier space: ρ ( ˆ U X ) ′′ = − ( ω ( λ + 2 µ ) + ω µ ) ˆ U X − ( λ + µ ) ω ω ˆ U Y ρ ( ˆ U Y ) ′′ = − ( ω ( λ + 2 µ ) + ω µ ) ˆ U Y − ( λ + µ ) ω ω ˆ U X ˆ U X (0) = ˆ F , ˆ U Y (0) = ˆ G, d ˆ U X dt (0) = 0 , d ˆ U Y dt (0) = 0 , where the derivatives denote differentiation with respect to time t .Following some algebra in which we rearrange to eliminate the transformed variable ˆ U Y interms of ˆ U X , we obtain the fourth order ordinary differential equation for ˆ U X (3.2). This equationis constant coefficient and linear: ρ ( ˆ U X ) IV + ρ ( A + C )( ˆ U X ) ′′ + ( AC − B ) ˆ U X = 0 , (3.2)where A = ω ( λ + 2 µ ) + ω µB = ( λ + µ ) ω ω C = ω ( λ + 2 µ ) + ω µ. It is now straightforward to show that ˆ U X ( t ) = P exp (cid:20) it q ( ω + ω ) µ/ρ (cid:21) + P exp (cid:20) − it q ( ω + ω ) µ/ρ (cid:21) + Q exp (cid:20) it q ( ω + ω )( λ + 2 µ ) /ρ (cid:21) + Q exp (cid:20) − it q ( ω + ω )( λ + 2 µ ) /ρ (cid:21) ˆ U Y ( t ) = − P ω ω exp (cid:20) it q ( ω + ω ) µ/ρ (cid:21) − P ω ω exp (cid:20) − it q ( ω + ω ) µ/ρ (cid:21) + Q ω ω exp (cid:20) it q ( ω + ω )( λ + 2 µ ) /ρ (cid:21) + Q ω ω exp (cid:20) − it q ( ω + ω )( λ + 2 µ ) /ρ (cid:21) , r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... where the constants P , P , Q and Q are found from the Fourier transform of the initialconditions, ˆ F and ˆ G . Solving for these initial conditions gives the complete solution ˆ U X ( t ) = ω ( ˆ F ω − ˆ Gω ) ω + ω cos (cid:18) t r µρ ( ω + ω ) (cid:19) + ω ( ˆ F ω + ˆ Gω ) ω + ω cos (cid:18) t s λ + 2 µρ ( ω + ω ) (cid:19) ˆ U Y ( t ) = − ω ( ˆ F ω − ˆ Gω ) ω + ω cos (cid:18) t r µρ ( ω + ω ) (cid:19) + ω ( ˆ F ω + ˆ Gω ) ω + ω cos (cid:18) t s λ + 2 µρ ( ω + ω ) (cid:19) (3.3)in the Fourier-transformed space. (i) Initial Conditions It remains to determine the transformed quantities ˆ F and ˆ G in (3.3), through an appropriatechoice of initial conditions f ( x, y ) and g ( x, y ) in the physical space. Here we choose u X ( x, y, ≡ f ( x, y ) = F ∂∂x ( e − ( x + y ) /a )= − F xa e − ( x + y ) /a u Y ( x, y, ≡ g ( x, y ) = G ∂∂y ( e − ( x + y ) /a )= − G ya e − ( x + y ) /a , (3.4)for given constants F , G , a . This implies that the initial conditions in the Fourier space are ˆ F ( ω , ω ) = − F a π Z ∞−∞ Z ∞−∞ xe − ( x + y ) /a − iω x − iω y d x d y = iF a ω e a ( ω + ω )ˆ G ( ω , ω ) = − G a π Z ∞−∞ Z ∞−∞ ye − ( x + y ) /a − iω x − iω y d x d y = iG a ω e a ( ω + ω ) Inserting these initial conditions into (3.3) gives the following equations in Fourier space: ˆ U X ( t ) = ia ω ω ω + ω ) e a ( ω + ω ) ( F − G ) cos (cid:18) t r µρ ( ω + ω ) (cid:19) )+ ia ω F ω + G ω ω + ω e a ( ω + ω ) cos (cid:18) t s λ + 2 µρ ( ω + ω ) (cid:19) ˆ U Y ( t ) = ia ω ω ω + ω ) e a ( ω + ω ) ( G − F ) cos (cid:18) t r µρ ( ω + ω ) (cid:19) )+ ia ω F ω + G ω ω + ω e a ( ω + ω ) cos (cid:18) t s λ + 2 µρ ( ω + ω ) (cid:19) (3.5) r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... After some algebra, the solution in physical space is found from the inverse Fouriertransformation of (3.5) to be u X ( x, y, t ) = − a xr ( x F + y G ) B ( r, t ) − a ( F − G ) (cid:26) xy r A ( r, t )+ 2 x ( x − y ) r (cid:2) A ( r, t ) − B ( r, t ) (cid:3) − x ( x − y ) r (cid:2) A ( r, t ) − B ( r, t ) (cid:3)(cid:27) u Y ( x, y, t ) = − a yr ( x F + y G ) B ( r, t ) − a ( F − G ) (cid:26) x yr A ( r, t )+ 2 y ( y − x ) r (cid:2) A ( r, t ) − B ( r, t ) (cid:3) − y ( y − x ) r (cid:2) A ( r, t ) − B ( r, t ) (cid:3)(cid:27) , (3.6)where we have introduced the usual radial distance, r = p x + y . For readability, we have alsointroduced the following functions of r, t : A ( r, t ) = Z ∞ e − a k cos( kt √ α ) J ( kr ) d kA ( r, t ) = Z ∞ e − a k cos( kt √ α ) k J ( kr ) d kA ( r, t ) = Z ∞ e − a k cos( kt √ α ) kJ ( kr ) d kB ( r, t ) = Z ∞ e − a k cos( kt p β ) J ( kr ) d kB ( r, t ) = Z ∞ e − a k cos( kt p β ) k J ( kr ) d kB ( r, t ) = Z ∞ e − a k cos( kt p β ) kJ ( kr ) d k (3.7)and α = µ/ρ and β = ( λ + 2 µ ) /ρ are the squares of the azimuthal and radial wave speedsrespectively. The functions J and J are the zeroth-order and first-order Bessel functions of thefirst kind. We now have an exact solution which can be plotted over time, and will be used toinspect the accuracy of the numerical approach in section 4. This analytical solution requires thenumerical evaluation of integrals in the A and B functions. These can be evaluated extremelyrapidly and with any desired accuracy using numerical integration, such as the trapezoidalmethod or Gauss-Legendre quadrature. (ii) Closed-Form Solution Following some experimentation with various initial conditions, a particular initial displacementwas found to yield a closed-form solution. In addition to being useful in testing numericalschemes for accuracy, a closed-form solution is of theoretical and mathematical interest. Theseelasto-dynamical equations are of sufficient complexity that non-trivial closed-form solutionsmay be elusive. In the previous section, the chosen initial condition was the derivative of asmoothed Gaussian. This resulted in an analytic solution in terms of the integrals given in (3.7). r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... Here, we instead choose the following initial condition for the displacement: u X ( x, y, ≡ f ( x, y ) = ∂∂x ( 1( x + 1)( y + 1) )= − x ( x + 1) ( y + 1) u Y ( x, y, ≡ g ( x, y ) = ∂∂y ( 1( x + 1)( y + 1) )= − y ( x + 1)( y + 1) . (3.8)Following the Fourier transform of this initial condition, and solution of equations (3.3), theinverse Fourier transform yields the solutions u X ( x, y, t ) = 14 i Z π cos φ (cid:2)(cid:0) A + (cid:1) + (cid:0) A − (cid:1) (cid:3) d φu Y ( x, y, t ) = 14 i Z π sin φ (cid:2)(cid:0) A + (cid:1) + (cid:0) A − (cid:1) (cid:3) d φ, where A ± = | cos φ | + | sin φ | + ix cos φ + iy sin φ ± βt. (3.9)Again, β = ( λ + 2 µ ) /ρ is the square of the radial wave speed. Consideration of the four quadrantsis used to eliminate the absolute value signs, yielding the sum of eight integrals: u X ( x, y, t ) = 14 i X m = ± X n = ± X p = ± Z π/ m cos φ (cos φ + sin φ + i ( mx cos φ + py sin φ + nβt )) d φu Y ( x, y, t ) = 14 i X m = ± X n = ± X p = ± Z π/ m cos φ (cos φ + sin φ + i ( my cos φ − px sin φ + nβt )) d φ. It may be seen that these terms occur as four pairs of complex conjugates, so that the resultingdisplacement is entirely real. However, expanding in this way results in a cumbersome number ofterms, so the integrals are performed as they are. The solutions are included in an appendix. Thesesolutions are easily coded in Mathematica, and the results are shown alongside the correspondingnumerical solutions in Fig.3 below.
4. Numerical Approach
The modelling of the seismic equations (2.6) for a more general initial condition, or for morecomplicated functions for density and Lamé parameters requires an efficient and accuratenumerical method. The advantages and disadvantages of the most common methods aredescribed in [9]. We present here an alternative method, developed in the context of fluidflow studies [16]. This is a spectral method, whereby variables are represented as weightedsums of analytic functions. As with the pseudospectral method [9] the spatial derivatives aregenerated analytically from the spectral basis functions. This provides two advantages overa finite difference scheme in that derivatives are known and cached beforehand, which is anotherwise time consuming numerical procedure, and calculation of derivatives in this wayexhibits exponential convergence [17], pp. 45-46. However, unlike traditional pseudospectralschemes, our method is easily parallelisable, lending itself to efficient solution over large numbersof simple processing units, as found on modern graphics cards. Thus, we get the accuracyadvantages of the spectral method, combined with computational efficiency approaching thatof a finite difference scheme. The main disadvantage of a spectral method such as this, is thatit may perform poorly when dealing with discontinuities or very rapid changes in wave-speedparameters, unless those discontinuities are explicitly allowed for in the choice of basis functions. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... The numerical calculations are limited to a rectangular region, x < x < x , y < y < y . Withinthis region, we choose to represent the displacement as u X ( x, y, t ) = M X m =0 N X n =0 A mn ( t ) cos (cid:2) α m ( x − x ) (cid:3) cos (cid:2) α n ( y − y ) (cid:3) u Y ( x, y, t ) = M X m =0 N X n =0 B mn ( t ) cos (cid:2) α m ( x − x ) (cid:3) cos (cid:2) α n ( y − y ) (cid:3) , (4.1)where α m = mπ/ ( x − x ) and α n = nπ/ ( y − y ) . This particular choice of basis functionsallows the displacement to have any value at any point in the region and does not force anybi-lateral symmetry. However, the orthogonality condition of the Fourier series is preserved,allowing A mn and B mn to be obtained accurately and efficiently from u X and u Y . In orderto solve (2.6) numerically, the two second order equations are replaced with four first orderequations, resulting in two more sets of coefficients C mn , D mn for the velocities ∂u X /∂t, ∂u Y /∂t .The total system of first order equations is thendd t A mn B mn C mn D mn = C mn D mnγ mn v R x x R y y ¨ u X cos (cid:2) α m ( x − x ) (cid:3) cos (cid:2) α n ( y − y ) (cid:3) d y d x γ mn v R x x R y y ¨ u Y cos (cid:2) α m ( x − x ) (cid:3) cos (cid:2) α n ( y − y ) (cid:3) d y d x (4.2)To be clear, in these equations, we have taken the two second-order equations (2.6), representedthem as Fourier series, and are solving for the coefficents ( A mn ( t ) etc.) via these four first orderequations (4.2). This formulation allows for solution by explicit time-integration methods, such asthe Runge-Kutta scheme used in this current work. The symbol γ mn is due to the orthogonalitycondition of the Fourier series, and is if m and n are both zero, if exactly one of m or n is zero, otherwise. The volume of the region (or area in this 2D case) is v = ( x − x )( y − y ) . The valuesof ¨ u X and ¨ u Y at each time step are determined by the wave equation (2.6). The derivatives inequation (2.6) are calculated analytically from the Fourier representations (4.1). The trigonometricbasis functions in (4.1) are calculated only once and cached in memory for use throughout therunning of the algorithm.It is necessary in numerical calculations to avoid non-physical reflections, which may resultfrom the boundaries of the calculation region. To achieve this, a simple robust and efficientabsorbing boundary was developed. The region of interest is surrounded by a region in whichan increasing amount of damping is applied to the displacement, in proportion to the velocityof the disturbance at that point. In order to effect this, we initially added a damping term to thedisplacement equations (2.6). A damping field was implemented which increased exponentiallythroughout the absorbing region. This was applied to the displacement at each time stepaccording to: u Xd ( x, y ) = u X ( x, y ) − ˙ u X ( x, y ) D ( x, y ) u Yd ( x, y ) = u Y ( x, y ) − ˙ u Y ( x, y ) D ( x, y ) , where the subscript d refers to the (new) damped displacement, u X , u Y are the displacementcomponents, calculated from the forward integration of (2.6), ˙ u X and ˙ u Y are the velocitycomponents, and D is the damping field. Specifically, D is zero in the region of interest, but hasvalue D = p ( e p d − in the boundary layer, where d is the depth into the absorbing region.The optimal values of the damping parameters p , p were determined using the Nelder-Meaddownhill simplex method [18]. By damping in proportion to the velocity normal to the boundaries(i.e. ˙ u X and ˙ u Y ), we are effectively implementing a partial matching condition for waves ofvarying speeds. This will be useful in models where the wavespeed varies throughout the region.As an alternative method, we have also implemented the perfecly matched layer (PML) approachof [19]. In this method the second-order kinematic equations (2.6) are modified to r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... ρ (cid:2) ∂ u X ∂t + ( β x + β y ) ∂u X ∂t + β x β y u X (cid:3) = ∂∂x (cid:2) ( λ + 2 µ ) ∂u X ∂x + λ ∂u Y ∂y + ( β y − β x ) w (cid:3) + ∂∂y (cid:2) µ ( ∂u Y ∂x + ∂u X ∂y ) + ( β x − β y ) w (cid:3) ρ (cid:2) ∂ u Y ∂t + ( β x + β y ) ∂u Y ∂t + β x β y u Y (cid:3) = ∂∂y (cid:2) ( λ + 2 µ ) ∂u Y ∂y + λ ∂u X ∂x + ( β x − β y ) w (cid:3) + ∂∂x (cid:2) µ ( ∂u X ∂y + ∂u Y ∂x ) + ( β y − β x ) w (cid:3) . These equations introduce six new field variables: β x ( x ) and β y ( y ) are one-dimensionalfunctions, which are zero within the physical region, and increase throughout the damping regionas β x ( x ) = q ( x d − xx d − x ) q (4.3)and similar for β y ( y ) . The factor q and the index q are free parameters to be chosen. The w , w , w , w are functions of x, y, t which are initially set to zero, and evolved through timeaccording to the following equations: ∂w ∂t + β x w = ( λ + 2 µ ) ∂u X ∂x∂w ∂t + β y w = µ ∂u X ∂y∂w ∂t + β x w = µ ∂u Y ∂x∂w ∂t + β y w = ( λ + 2 µ ) ∂u Y ∂y . An additional "scaling coefficient" in each direction is employed in [19], but we found that itsuse led to instability and increased reflections. This problem was indicated in [20], along withthe recommendation, which we have followed here, not to use the scaling coefficient, unlessthe particular problem necessitated this due to certain instabilities. We thus have another two-parameter damping layer, the parameters in this case being q and q . Again, the downhillsimplex method was used to determine the optimum values for the two parameters, forvarious thicknesses of the damping boundary. Table 1 shows the greatest difference betweenthe numerical and analytic solutions at time t = 5 , for the two methods, along with the optimalvalues found ( p and p for the simple method, q and q for the PML). This is for the × simulation shown in Fig. 1 in section 5. Some of these 560 grid points are used for the absorbingboundary layer, from 14 to 56, as shown.It is clear from the results in Table. 1 that the PML boundary is far more effective at dampingreflections, producing a decrease in reflections by a factor of between 10 and 40 over the simplermethod. The additional cost for forward integration of the auxiliary variables resulted in nomore than twice the calculation time. Because of this, we have used the PML formulation forthe numerical calculations in this paper. However, we have included the simpler method asan alternative for situations where simplicity of implementation is valued over the additionaleffectiveness of the PML method.In all numerical simulations, the time step-size is reduced until any changes in the solutionare of an acceptable level (see [21] p.373-374 regarding stability of Runge-Kutta methods). In thefigures presented in this paper the model is no longer visibly changing with further reduction instep-size. Additionally, by comparing with analytical solutions, the step size was reduced untildifferences between analytical and numerical solutions are of the same order as the errors due tothe spurious boundary reflections. If further refinement is required, it is not difficult to implement r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... Table 1.
Effectiveness of the two types of damping layer as in Fig. 1. The total grid size is × , with the dampingboundary ranging from 14 to 56 grid points. The greatest absolute difference between the analytic and numerical solutionsis shown for time t = 5 , when the main body of the original impulsive wave has left the region. number of points 14 28 56Simple Damping p p q q . th − th Runge-Kutta method which adjusts the step-size automatically, such that the fourth-and fifth-order solutions differ by less than a specified tolerance.
5. Analytic-Numerical comparison
The numerical solution can now be compared against the analytical model as a measure of theaccuracy of the numerical approach. The first check will be with an initial condition which iscircularly symmetric. Using the initial condition in (3.4), we set F = G = 1 , and a = 0 . . Theanalytic solution is computed at time t = 1 , t = 2 , t = 3 , t = 4 (seconds), in an x − y grid of × regularly spaced points. A trapezoidal method was used to evaluate the integrals in (3.7)over k from to using integration points. These values for the truncation point of theintegral and the number of points were chosen by increasing them until there was no change inthe evaluation at double precision. The wave speeds V p = p β = s λ + 2 µρ ; V s = √ α = r µρ (5.1)were set at V p = 1 . km/s , V s = 0 . km/s and density ρ = 1000 kg/m throughout the space. Thenumerical solution was then calculated using the same parameters, and the same spatial grid. Theouter points on each edge of the space were used for the PML damping layer, implemented asdescribed in section 4. The classic fourth-order Runge-Kutta method was used with a timestepof 0.005 seconds. Fig. 1 shows a plot of the two methods, with the bottom row showing thedifference between the two. Only the undamped region is plotted, being the inner × points. It may be seen from the colour scales that the difference between the two methods islimited to approximately 0.1% of the initial amplitude.A second check of the numerical scheme is shown in Fig. 2 using the analytic solution witha different initial condition. In this case, we set initial condition parameters to be F = 1 , G = 0 .All other parameters remained the same as in Fig. 1. In this case the initial condition gives riseto both a compression wave and a slower moving shear wave. Again, the numerical solutionand analytic solution were run with identical initial condition, and the difference plotted. Thenumerical solution shows excellent agreement with the analytic model. It can be seen from thefirst frame of Fig. 2 that the amplitude of the pulse when first contacting the boundary has anamplitude in excess of , while reflected amplitudes are around the . level. As with Fig.1,the differences may be further reduced by increasing the number of grid points, spectral modes,integration steps and the absorbing boundary thickness. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... Figure 1.
A comparison of the expanding compressional wave at times (in seconds) from t = 1 (left) to t = 4 (right). Thetop row shows the analytic solution (3.6), the middle row shows the output from the numerical approach of section 4 andthe bottom row shows the difference between the two. The x and y axes are in kilometres and the z -axis shows themagnitude of the displacement based on a pulse with initial amplitudes F = G = 1 . In the numerical model, absorbingboundaries have been used to damp reflections. Figure 2.
A similar comparison of the numerical and analytical solutions as in Fig. 1. In this case the chosen initialcondition gives rise to both fast compression waves and slower moving shear waves. In the bottom row the dampedreflections from the numerical model are small but visible differences between the numerical and analytic solutions.
A final comparison is shown in Fig 3. This is the closed form solution, for initial displacement u X ( x, y,
0) = − x ( x + 1) ( y + 1) u Y ( x, y,
0) = − y ( x + 1)( y + 1) . r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... The closed-form solution was calculated for the specified times, from 0 to 4, in Mathematica.The numerical solution was computed using the numerical scheme described above, marchedforward in time in Fortran using parallel CUDA. The data for both was then plotted in M
ATLAB .The grid size is from -8 to 8 in both directions.
Figure 3.
A comparison of the numerical solution with the closed-form solution described in Section 3. This initial conditiongives rise to compression waves only. It may be seen qualitatively that the numerical solution matches the closed-formsolution well.
6. Variable Wavespeed Parameters
Having confidence that the numerical scheme is accurate in the case of constant Lamé parameters,we now move into the realm of varying parameters, for which we do not have any analyticsolutions. While we cannot compare directly in these cases, we are using exactly the samenumerical formulation as we used when comparing against the known constant parametersolutions shown in Figs. 1 and 2. We thus have at least some confidence that such a numericalscheme is a reasonable approach to inhomogeneous media. We first show a simulation of aninitially circular wavefront propagating through a medium with two different sets of wavespeeds.These wavespeeds change smoothly but rapidly across the boundary, as given below. Theupper right half of the space has wavespeeds v s = p µ /ρ = . km/s, v p = p ( λ + 2 µ ) /ρ =1 . km/s and the lower left half has v s = p µ /ρ = 1 km/s, v p = p ( λ + 2 µ ) /ρ = 2 . km/s .The equations really only have two independent parameters, so ρ and ρ were scaled out(set to 1) in the code. These scaled Lamé parameters were defined throughout the region as µ = µ + ( µ − µ ) / (1 + exp[ − x + y )]) , and similarly for λ . This change is shown in Fig.4where the value of µ changes from µ to µ over the interfacial region. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... Figure 4.
The scaled Lamé parameters change smoothly between the two regions in Fig.5. The change in µ is shown hereaccording to µ = µ + ( µ − µ ) / (1 + exp[ − D ]) , where D is the distance from the line at y = − x . This smoothchange results in a model containing refraction but not reflections (see Fig.5). The results of this simulation can be seen in Fig. 5, where the wavefront is seen to be travellingfaster after crossing the interface into the lower left region. A caustic has formed near the diagonalline due to the rapidly changing refractive index in the interfacial region. As with the simulationsin Fig. 1 and 2, a regular x − y grid of × points was used. We note that a steeper parametergradient requires a finer grid in order to avoid the Gibbs-type phenomena often associated withspectral methods. Figure 5.
A simulation of the propagation of an initially circular wave through a medium having two distinct sets of waveparameters, which change rapidly but smoothly across a diagonal line. The wave can be seen travelling more slowly inthe upper right half of the plane. Times are given in seconds, distances are given in kilometres. The upper row shows theabsolute magnitude of the displacement. The lower row shows the corresponding contour lines, along with a diagonal lineshowing the location of the "interface" and the star shows the centre of the initial disturbance.
A final simulation is shown in Fig. 6, where the density and wavespeeds increase rapidlywith depth. These parameters are loosely based on the seismic structure of the upper level ofan antarctic ice-sheet, in which the density in the packed snow near the top is quite low, butincreases rapidly over a vertical distance of approximately 400 metres [22,23]. To approximatesuch a situation, we have set the nominal upper density to ρ = 300 kg/m and wavespeeds to v p = 0 . km/s and v s = 0 . km/s . Lower nominal values are ρ = 900 kg/m with v p = 3 . km/s r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... and v s = 1 . km/s . These wavespeeds are related to the Lamé parameters by equation (5.1).These nominal values were used to create functions for density and wavespeeds increasingsmoothly with depth. The vertical profile for density was set to ρ ( y ) = ρ − ( ρ − ρ ) exp(20 y ) with − . < y < in kilometres. The wavespeed profiles were set in the same way. Absorbingboundary layers were implemented on left and right, and on the bottom edge, but not on the top.This allows the top boundary to act as a reflecting surface, although in this simple example, wehave not specifically implemented zero-normal-stress boundary conditions. The rapid change inwavespeed causes downward travelling signals generated near the surface to be refracted back upand thus to bounce along the underside of the surface at y = 0 . A row of sensors has been placed atthe top of the space, and the x and y velocities recorded for each sensor at each integration step.The formulation of the numerical scheme calculates both displacements and velocities at everytime step. In this case we show the velocities at the sensor points. This sensor data is shown in thelower panel, produced using the "wiggle" M ATLAB function [24]. This simulation was performedin a regularly spaced x − y grid × points.In the top frame, the initial, circular disturbance is shown. This disturbance expands, but thevariation in wavespeed causes the wavefront to distort, so that in the second frame, the front hasalready refracted back up to the top. The front then reflects off the surface at y = 0 , and proceedsin a downward direction, which again refracts back up to the top in the third frame and again inthe fourth. This series of reflections is seen in the series of wavefronts detected by the sensors inthe lower panel. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... Figure 6.
A simulation of the propagation of an initially circular wave through a medium having wavespeed increasingrapidly with depth is shown in the upper panel. The wavefront can be seen refracting back up to the surface. The topfigure shows the initial disturbance at t = 0 . Each succeeding figure is . seconds later than the previous. Distancesare given in kilometres. A row of sensors is placed at the top of the space. In the lower panel, the velocity data for each ofthe sensors is shown. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... Technical Details
All of the simulations were run on a computer employing an Intel I7-7700 CPU, 32 GB RAM and aQuadro GP100 GPU. The operating system was Ubuntu 18.04. The computations were performedusing Fortran, compiled using the PGI Fortran compiler [25], with most of the calculations beingcarried out by the GPU. The displacements u X and u Y were written to file and later readinto M ATLAB to produce the figures. All calculations were performed using double precisionarithmetic. For the × grid calculations in this paper, the analytic calculations took only 7seconds to run, this time being primarily for the calculations of the integrals (3.7). The numericalcalculations of Fig. 1 to 6 took 3-4 minutes. These consisted of 4 units of time, each unit using 280iterations through the four steps of the Runge-Kutta algorithm. A much smaller number of stepsmay be used, and doing so produces results which are not noticeably different, but this numberof steps was required to bring the difference between analytic and numerical models down to thelevel of the non-physical reflections (in our case, down to approximately 0.001). In all simulationsthe number of Fourier modes in each direction was one fifth of the number of grid points inthat direction. The majority of required computations are for the solution of equation (4.2). By re-writing the sums in these equations as pure Matrix products, they may be performed in parallelusing CUDA code for multiplying matrices efficiently on the GPU. This code was adapted fromthe sample provided in the PGI Fortran compiler user guide [25], and significantly improved forthe particular task described in this paper. These algorithmic changes reduced the runtime by atleast an order of magnitude.
7. Summary
For the case of homogeneous media, an exact solution to the elastic wave equation in twodimensions has been developed, with general initial conditions. For initial conditions of somesymmetry, along with infinite boundary conditions, an analytic function was found whichdescribes the displacement components for any chosen time. This provides an excellent test forthe accuracy of similar numerical methods.By combining powerful numerical techniques, a new scheme has been implemented forforward integration of the elastic wave equation. The spectral (Fourier) method allows forderivatives to be calculated immediately and analytically from the basis functions. The classicfourth order Runge-Kutta method enables quick and accurate integration. Evaluation ofmassively parallel systems of equations is well suited to computation on modern graphicshardware, which is implemented using the chosen Fortran compiler.Combining this numerical scheme with a suitable absorbing boundary layer allows forforward modelling of elastic waves. In the case of a homogeneous medium, the numerical methodhas been verified against the analytic solution. Waveforms from two inhomogeneous models areshown as illustrative examples.This numerical method and its implementation on GPU architecture is presented for use byresearchers where there is a need to investigate large numbers of relatively simple seismic models,tailored to environmental applications. The ability to produce results rapidly using desktopcomputing and the agility in use allows for experimentation over a range of model parameters.Data Accessibility.
The Fortran and M
ATLAB code for the simulations in this paper is available at https://github.com/StephenJWalters/math-seismic
Authors’ Contributions.
SJW wrote the code, developed the numerical model and prepared the initialdraft of the paper, LKF derived the analytic solution, AMR devised the overarching research program, andprovided the seismological applications context. All authors were involved in the revision and approval ofthe final manuscript.
Competing Interests.
We have no competing interests. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... Funding.
This research was supported under Australian Research Council’s Special Research Initiative forAntarctic Gateway Partnership (Project ID SR140300001), and Discovery Program (DP190100418).
Acknowledgements.
We are grateful to Professor Shaolin Liu and three anonymous referees, whoseconstructive comments have improved this paper considerably.
References
1. Aki K, & Richards PG. 2002
Quantitative seismology.
2. Kennett BL. 2001
The Seismic Wavefield: Volume 1, Introduction and Theoretical Development.
Cambridge University Press.3. Lamb H. 1904
On the propagation of tremors over the surface of an elastic solid.
Phil. Trans. R. Soc.Lond. A 203, 1-42.4. Garvin WW. 1956
Exact transient solution of the buried line source problem.
Proc. R. Soc. A234(1199), 528-541.5. Kausel E. 2013
Lamb’s problem at its simplest . Proc. R. Soc. A, 469(2149), doi:10.1098/rspa.2012.0462.6. Gosselin-Cliche B, & Giroux B. 2014
3D frequency-domain finite-difference viscoelastic-wavemodeling using weighted average 27-point operators with optimal coefficients.
Geophysics 79(3),T169-T188.7. Carcione JM. 1993
Seismic modeling in viscoelastic media.
Geophysics 58(1) 110-120.8. Diaz J, & Ezziani A. 2010
Analytical solution for waves propagation in heterogeneous acoustic/porousmedia. Part I: the 2D case.
Communications in Computational Physics, 7(1), 171.9. Igel H. 2017
Computational seismology: a practical introduction.
Oxford University Press.10. Fichtner A, Igel H, Bunge HP, & Kennett BL. 2009
Simulation and inversion of seismic wavepropagation on continental scales based on a spectral-element method.
JNAIAM 4(1-2), 11-22.11. Komatitsch D & Vilotte JP. 1998
The spectral element method: an efficient tool to simulate the seismicresponse of 2D and 3D geological structures.
Bulletin of the seismological society of America, 88(2),368-392.12. Maeda T, Takemura S & Furumura T. 2017
OpenSWPC: an open-source integrated parallelsimulation code for modeling seismic wave propagation in 3D heterogeneous viscoelastic media.
Earth,Planets and Space, 69(1), 102.13. Sens-Schönfelder C, & Wegler U. 2011
Passive image interferometry for monitoring crustal changeswith ambient seismic noise.
Comptes Rendus Geoscience, 343(8-9), 639-651.14. Tsai VC, Minchew B, Lamb MP & Ampuero JP. 2012
A physical model for seismic noise generationfrom sediment transport in rivers.
Geophysical Research Letters, 39(2).15. Paul Winberry J, Anandakrishnan S, Wiens DA, & Alley RB. 2013
Nucleation and seismic tremorassociated with the glacial earthquakes of Whillans Ice Stream, Antarctica.
Geophysical ResearchLetters, 40(2), 312-315.16. Walters SJ & Forbes LK. 2019
Fully 3d Rayleigh-Taylor instability in a Boussinesq fluid.
ANZIAMJ. 61(3), 286-304.17. Boyd JP. 2001
Chebyshev and Fourier spectral methods.
Dover, New York.18. Nelder JA & Mead R. 1965
A simplex method for function minimization . The Computer Journal,7(4), 308-313.19. Assi H, & Cobbold RS. 2017
Compact second-order time-domain perfectly matched layer formulationfor elastic wave propagation in two dimensions.
Mathematics and Mechanics of Solids, 22(1), 20-37.20. Assi H. 2016
Time-domain modeling of elastic and acoustic wave propagation in unbounded media,with application to metamaterials (Doctoral dissertation, University of Toronto, Canada).21. Atkinson KE. 1978
An introduction to numerical analysis
John Wiley & Sons.22. Reeh N, Fisher DA, Koerner RM, & Clausen HB. 2005
An empirical firn-densification modelcomprising ice lenses . Annals of Glaciology, 42, 101-106.23. Schlegel R, Diez A, Löwe H, Mayer C, Lambrecht A, Freitag J, Miller H, Hofstede C & EisenO. 2019
Comparison of elastic moduli from seismic diving-wave and ice-core microstructure analysis inAntarctic polar firn r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A ..........................................................
8. Appendix
After performing the integrals in equation (3.9), the displacements are written for readability inthe form u X ( x, y, t ) = X m = ± X n = ± X p = ± (cid:18) S + S + S + S + S D D D + S S D (cid:19) , where the auxiliary functions are S = m (cid:20) inβt (cid:16) − β t (cid:17) + 2 β t − β t + r − r + 4 (cid:21) S = nβ t (4 x − mpy ) S = 3 β t (cid:16) mx − my + pxy + impy − ix (cid:17) S = nβt [2 x − xy − x + mpy (4 x + y − − i (3 pxy + 5 mx + my )] S = − pxy − i ( r − x + mpy ) S = 3 nβt ( x − mi ) S = ArcTanh ( i − pyD ) + ArcTanh ( nβt − mx + pyD ) D = q r − β t − − i ( mx + py ) D = nβt + mx − iD = nβt + py − ir = p x + y Unsurprisingly, the solution for u Y is identical to that for u X but with ( x, y ) replaced by ( y, − x ) : u Y ( x, y, t ) = X m = ± X n = ± X p = ± (cid:18) T + T + T + T + T E E E + T T E (cid:19) , where T = m (cid:20) inβt (cid:16) − β t (cid:17) + 2 β t − β t + r − r + 4 (cid:21) T = nβ t (4 y + mpx ) T = 3 β t (cid:16) my − mx − pxy − impx − iy (cid:17) T = nβt [2 y − yx − y − mpx (4 y + x − − i (5 my + mx − pxy )] T = 8 pxy − i ( r − y − mpx ) T = 3 nβt ( y − mi ) T = ArcTanh ( i + pxE ) + ArcTanh ( nβt − my − pxE ) E = q r − β t − − i ( my − px ) E = nβt + my − iE = nβt − px − ir = p x + yy