{'fraction_non_alphanumeric': 0.036212691134901734, 'fraction_numerical': 0.012130249949920987, 'mean_word_length': 4.584141188168034, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Upwind summation by parts finite difference methods for large scale elastic wave simulations in 3D complex geometries

High-order accurate summation-by-parts (SBP) finite difference (FD) methods constitute efficient numerical methods for simulating large-scale hyperbolic wave propagation problems. Traditional SBP FD operators that approximate first-order spatial derivatives with central-difference stencils often have spurious unresolved numerical wave-modes in their computed solutions. Recently derived high order accurate upwind SBP operators based upwind FD stencils have the potential to suppress these poisonous spurious wave-modes on marginally resolved computational grids. In this paper, we demonstrate that not all high order upwind SBP FD operators are applicable. Numerical dispersion relation analysis shows that odd-order upwind SBP FD operators also support spurious unresolved high-frequencies on marginally resolved meshes. Meanwhile, even-order upwind SBP FD operators (of order 2, 4, 6) do not support spurious unresolved high frequency wave modes and also have better numerical dispersion properties. We discretise the three space dimensional (3D) elastic wave equation on boundary-conforming curvilinear meshes. Using the energy method we prove that the semi-discrete approximation is stable and energy-conserving. We derive a priori error estimate and prove the convergence of the numerical error. Numerical experiments for the 3D elastic wave equation in complex geometries corroborate the theoretical analysis. Numerical simulations of the 3D elastic wave equation in heterogeneous media with complex non-planar free surface topography are given, including numerical simulations of community developed seismological benchmark problems. Computational results show that even-order upwind SBP FD operators are more efficient, robust and less prone to numerical dispersion errors on marginally resolved meshes when compared to the odd-order upwind and traditional SBP FD operators.


Introduction
High fidelity numerical simulations of seismic (elastic) waves are prevalent in many applications such as earthquake engineering, natural minerals and energy resources exploration, strong-ground motion analysis, and underground fluid injection monitoring. Seismic waves emanating from geophysical events can propagate hundreds to thousands of kilometres interacting with geological structure and complicated free-surface topography. Exploration seismology and natural earthquake hazard mitigation increasingly rely on multi-scale and high frequency (0-20 Hz) simulations. Often, surface and interface waves [17] are the largest amplitude wave modes and are arguably the most important wave modes in the medium. Therefore accurate and efficient numerical simulation of seismic surface and interface waves, and scattering of high-frequency wave modes by complex non-planar free-surface topography, are critical for assessing and quantifying seismic risks and hazards [8].
In this study, we derive efficient high order accurate numerical methods for large scale numerical simulations of seismic waves in complex geometries. We consider the elastic wave equation in the first-order form, where the unknowns are the particle velocity vector and stress fields. The main classes of numerical methods for the solution of time-dependent PDEs are the spectral method, finite element (FE) method, finite difference (FD) method, finite volume (FV) method and the discontinuous Galerkin (DG) finite element method. They all have different strengths and weaknesses. Computational efficiency has continued to make the use of FD methods on structured grids attractive. However, the presence of boundary conditions and complex non-planar free-surface topography make the design of stable and accurate FD methods challenging.
For computational seismology, staggered grids FD methods on Cartesian meshes are the industry standard because they are efficient and have optimal numerical dispersion properties. However, the design of high order accurate and stable staggered FD methods for the elastic wave equation in complex geometries is a challenge. Although some progress is being made in this direction [15].
For well-posed initial boundary value problems (IBVP), the summation-by-parts (SBP) FD [20] with the simultaneous approximation term (SAT) [13,11,14] technique for implementing boundary conditions enables the development of stable numerical approximations on smooth geometries. The methods can be extended to complex geometries using curvilinear coordinate transforms and multi-block schemes [4,21].
Traditional SBP FD operators are based on central finite difference formula with special one-sided boundary closures, designed such that the operator preserves the integration by parts principle [20]. Often, traditional SBP FD operators which approximate the spatial derivative suffer from spurious unresolved wave-modes in their numerical solutions. For marginally resolved solutions, these spurious wave-modes have the potential to destroy the accuracy of numerical solutions for a first-order hyperbolic partial differential equation, such as the elastic wave equation.
To ensure the accuracy of numerical solutions of elastic wave equations in complex geometries, we discretise the 3D elastic wave equation with a pair of non-central (upwind) FD stencils [12], on boundaryconforming curvilinear meshes. The main benefit for these operators [12,3] is that they have the potential to suppress poisonous spurious oscillations from unresolved wave-modes, which can destroy the accuracy of numerical simulations. However, these operators are asymmetric and dissipative, can potentially destroy symmetries that exist in the continuum problem. For example, the linear elastic wave equation with freesurface boundary conditions preserves the mechanical energy for all time. It is imperative that a stable numerical approximation must preserve these symmetries, by mimicking the corresponding continuous energy estimate at the discrete level. Otherwise numerical simulations on marginally resolved meshes can be polluted by numerical artefacts, in particular for long time simulations. A good case in point is that important information in the medium propagated by scattered high frequency surface waves could be corrupted through numerical dissipation.
Our first objective is to carefully combine the upwind SBP operator pairs which have good dispersion relation properties so that we preserve the discrete anti-symmetric property and invariants of the underlying IBVP. We then achieve energy stability by imposing boundary conditions weakly with penalty terms, in a manner that leads to bounded discrete energy estimate. In particular, if the IBVP is energy conserving, the numerical approximation is also energy conserving. We derive a priori error estimates and prove the convergence of the numerical error.
The second objective of this study is to reveal the numerical dispersion properties of the discrete upwind SBP FD operators [12]. By the analysis of the numerical dispersion relation we show that not all the upwind SBP FD operators derived in [12] are capable of suppressing poisonous spurious unresolved wave modes for the elastic wave equation. In particular, we show that odd-order upwind SBP FD operators also support spurious unresolved high frequency wave modes on marginally resolved meshes. Even-order upwind SBP FD operators (of order 2, 4, 6) do not support spurious unresolved wave modes and have better numerical dispersion properties than traditional SBP FD operators and odd-order upwind SBP FD operators. The numerical dispersion relation reaches near optimal properties for the 6th order accurate upwind SBP FD operator. Beyond the 6th order accurate upwind SBP FD operator, higher order accuracy does not improve the dispersion properties of the upwind SBP operators. Numerical experiments for the 3D elastic wave equation in complex geometries verify the analysis.
Numerical solutions are integrated in time using the fourth-order accurate low-storage Runge-Kutta method [1]. The numerical method is implemented in WaveQLab [4], a petascale elastic wave solver. Simulations of elastic waves in heterogeneous media with free surface topography are presented, including the numerical simulation of community developed seismological benchmark problems. Our results show that even-order upwind SBP FD operators are more robust and less prone to numerical dispersion on marginally resolved meshes when compared to traditional SBP operators [20,4]. Consequently, more accurate numerical solutions on marginally refined meshes can be achieved with less computational effort by carefully using certain (6th order acurate) upwind SBP operators in a way that mimic the continuous problem.
We have used the perfectly matched layer (PML) [5] to enable efficient domain truncation and prevent artificial numerical reflections, from the computational boundaries, from contaminating the numerical sim-ulations. The effective absorption properties of the PML allows us to sufficiently limit the computational domain with only a few grid points around the computational boundaries where the PML is active. Thus saving as much as %97.9017 of the required computational resources for the 3D seismological benchmark problem [19]. Because of the asymmetric properties of the PML and the upwind SBP operators, a stable implementation of the PML for the 3D IBVP elastic wave equation using the upwind SBP operators is a non-trivial task. The details of the numerical treatment of the PML using upwind SBP operators will be reported in a forthcoming paper.
For large scale numerical simulations, it is important that the parallel numerical software is efficient and scalable. We have performed strong scaling tests, demonstrating nearly perfect scaling and verifying the efficiency of our parallel implementation of high order upwind SBP operators for large scale elastic wave simulations in 3D complex geometries.
The structure of the paper is as follows. In the next section, we present the elastic wave equation in general curvilinear coordinates and prove anti-symmetric properties that must be preserved. In section 3, we introduce general linear well-posed boundary conditions and derive continuous energy estimates. In section 4, we discretise in space, introduce upwind SBP operators, and approximate the elastic wave equation in space. Numerical boundary conditions are derived in section 5 as well as the derivation of semi-discrete energy estimates, proving stability. In section 6, we provide an error analysis for our discretised schemes. We also provide a dispersion relation analysis, showing that even-order upwind SBP operators have dispersion relations which better mimic the continuous problem when compared to their odd-order and traditional counterparts. In section 7, we present numerical simulations verifying the accuracy and demonstrating the efficacy of the numerical method in complex geometry, with geologically constrained non-planar free-surface topography. We also provide scaling tests that indicate our code achieves near-perfect strong scaling. In section 8, we summarise our findings and speculate on the directions for future work.

Preliminary
In this section we introduce the elastic wave equation in general curvilinear coordinates and derive invariants and anti-symmetric properties that must be preserved.

Physical model
Let Ω be a connected compact subset of R 3 with a piecewise-smooth boundary Γ, t > 0 be the time variable, σ := (σ xx , σ yy , σ zz σ xy , σ xz , σ yz ) T be the vector of stresses, and v := (v x , v y , v z ) T be the particle velocities. Throughout, we assume that the velocities and stresses are bounded functions of space and time.
The first order time-dependent elastic wave equations in a source free, heterogeneous medium are where S = S T > 0 is the compliance matrix given in Equation (3) and ρ : Ω → R + is the density of the medium. The first three equations in (1) describe the conservation of momentum, and the latter six encode time derivative of Hookes law in three space dimensions.
It is convenient to work with Equation (1) in its general basis form. Let e := {e x , e y , e z } be a basis for R 3 , then we may recast this PDE system into its conservative and non-conservative components by introducing the anti-symmetric form [7] where Q := (v, σ) T , with µ, λ ∈ R are independent Lamé parameters which describe an isotropic medium and e ξx σ xx + e ξy σ xy + e ξz σ xz e ξx σ xy + e ξy σ yy + e ξz σ yz e ξx σ xz + e ξy σ yz + e ξz σ zz with F := (F x , F y , F z ) T and the basis vectors e ξ = (e ξx , e ξy , e ξz ) T . For instance in the standard Cartesian co-ordinate systems we have which recovers Equation 1.
Lemma 2.1. Consider the anti-symmetric form given in Equation (2). For any basis e that spans Ω we have ∂Q ∂ξ Proof. Expanding the matrix multiplication and simplifying yields the result.
We will find Lemma 2.1 useful in proving stability in our discretisation scheme in complex geometries. It will often be the case that our basis will be defined locally, through finding a set of functions whose partial derivatives evaluated at each point in Ω span R 3 . For instance, consider the functions proj ξ : Ω → R for ξ ∈ {x, y, z} given through These functions form the (local) basis vectors so e ξ are the contravariant basis defined in (5).

Curvilinear coordinates
Assume Ω ⊂ R 3 to be sufficiently smooth such that it can be mapped to the unit cube Ω := [0, 1] 3 . See also Figure 1. If Ω is piece-wise smooth we can partition it into locally smooth sub-blocks and map each sub-block to the unit cube. For simplicity, we will only consider one sub-block here. Let Φ : Ω → Ω be a diffeomorphism and adopt the notation Φ(x, y, z) := (q(x, y, z), r(x, y, z), s(x, y, z)), Assume that Φ(Ω) = Ω, and furthermore the boundary interaction for a smooth function X : R 2 → R. For example X(y, z) could describe a complex free-surface topography. The Jacobian determinant for Φ −1 can be written as Here ξ η is the partial derivative ∂ξ/∂η for ξ, η ∈ {x, y, z, q, r, s}. Similarly, The spatial derivatives in the transformed coordinates have: the conservative form and the non-conservative form Although the conservative (8) and non-conservative (9) transformations of the spatial derivatives are mathematically equivalent, when discretised they give different approximations. Specifically, in the discrete setting, the conservative form (8) approximated with an SBP operator preserves the divergence theorem.

Remark 1. The coordinate transformation in Equation
for all ξ ∈ {q, r, s}. This will be crucial in deriving high order accurate, structure preserving and provably energy stable scheme for the elastic wave equation in complex geometries.

Boundary Conditions
In this section, we formulate linear well-posed boundary conditions in complex geometries. As shown in Figure 1, define the faces of the boundary Γ as for ξ ∈ {q, r, s} and i ∈ {0, 1}. Each of the F ξ,i are Lebesgue-almost disjoint and form the boundary of Ω, that is Similarly the faces of the computational boundary Γ are given through and the boundary Γ of Ω is made from the level sets of the functions q, r, s, For a point (x 0 , y 0 , z 0 ) ∈ F ξ,i , the unit normal vector to the surface F ξ,i is given by where ξ x , ξ y , ξ z are the partial derivatives of ξ with respect to x, y, z respectively. On the boundary surface, we extract the particle velocity vector and the traction vector, and the local where det(R) = 0 and R −1 = R T . Next, rotate the particle velocity and traction vectors into the local orthonormal basis, l , m and n, Plane p-waves and plane s-waves propagating along the normal vector n on the boundary are given by At the boundary faces F ξ,i defined in (14) we consider the linear boundary conditions, Here γ η are real parameters with 0 ≤ |γ η | ≤ 1. The boundary condition (20), can describe several physical situations. We have a free-surface boundary condition if γ η = 1, an absorbing boundary condition if γ η = 0 and a clamped boundary condition if γ η = −1. Note that the boundary condition (20) satisfy the inequalities We introduce the boundary terms which are surface integrals encoding the work done by the traction force on the boundary Finally, the identity (21) completes the proof of the lemma.
In the continuous setting, we can show that our PDE has finite energy controlled by the boundary terms, BTs. To begin, for real functions we introduce the L 2 inner product, and the corresponding energy-norm The weighted L 2 -norm Q (·, ·, ·, t) 2 P is the mechanical energy, which is the sum of the kinetic energy and the strain energy.
Expanding the right hand side and applying integration by parts yields which from Lemma 2.1 and 3.1 gives the result.
In the next section, we will derive a numerical approximation of the transformed elastic wave equation (11) in curvilinear coordinates subject to the boundary conditions (20). We will approximate spatial derivatives using upwind SBP operators and imposed boundary conditions using penalties. To guarantee numerical stability we will prove numerical results analogous to Theorem 3.2.

Discretisation
In this section, stable discrete numerical approximation of the transformed elastic wave equation (11) in curvilinear coordinates subject to the boundary conditions (20) is derived. We use upwind SBP operators [12] to discretise the spatial derivatives, and keep the time variable continuous. The upwind SBP operators come in pairs, the forward difference operator D + and the backward difference operator D − . Our ultimate goal is to carefully combine this pair of upwind SBP operators to preserve, in the discrete setting, the antisymmetric property given in Lemma 2.1, and derive a conservative scheme. Boundary conditions are then enforced weakly using penalties in a manner that leads to energy stability.

Discrete spatial derivative
We will use a reference mesh that is uniform across each of the axes to discretise the reference computational cube Ω = [0, 1] 3 . For each ξ ∈ {q, r, s}, consider the uniform discretisation of the unit interval ξ ∈ [0, 1] where n ξ + 1 is the number of grid-points on the ξ-axis.
We will use upwind SBP operators introduced in [12] to approximate the spatial derivatives, ∂ ∂ξ . By combining the operators D − , D + we can respect the integration by parts formula: This will be critical in deriving a stable and conservative numerical approximation of the elastic wave equation in complex geometries.
j > 0 for all j ∈ {0, 1, . . . , n ξ }. We mimic the integration by parts property through finding the dual pairing of linear operators D +ξ , D −ξ : R n ξ +1 → R n ξ +1 so that is negative semi-definite to introduce efficient numerical suppression of unresolved high frequency wave modes, for more details see [12,3].
We make the discussion more formal.
Similar to traditional SBP operators, upwind SBP operators have higher accuracy in the interior, away from the boundaries. The accuracy of the operators is lowered close to the boundaries where special boundary closures are used. Upwind SBP operators with even-order (2p)-th accuracy in the interior are closed with p-th order accurate stencils close to boundaries. Upwind SBP operators with odd-order (2p + 1)-th accuracy in the interior are closed with p-th order accurate stencils close to boundaries. These operators can yield (p + 1)-th global order of accuracy.
The 1D SBP operators can be extended to higher space dimensions using tensor products ⊗. Let f (q, r, s) denote a 3D scalar function, and f ijk := f (q i , r j , s k ) denote the corresponding 3D grid function. The 3D scalar grid function f ijk is rearranged row-wise as a vector f of length n q n r n s . For ξ ∈ {q, r, s} and η ∈ {−, +} define: where I n ξ is the identity matrix of size n ξ × n ξ , and we take χ k=ξ := χ {ξ} (k) and χ k =ξ := 1 − χ k=ξ . So D ξη approximates the partial derivative operator in the ξ direction. An inner product on where the corresponding discrete energy-norm is given by Further, we have the multi-dimensional SBP property ξ∈{q,r,s} where S ξ (f , g) in the right hand side is the surface cubature, defined by Note that ξ 0 = 0 and ξ n ξ = 1, for all ξ ∈ {q, r, s}.

Numerical approximation in space
Consider the transformed elastic wave equation (11) in curvilinear coordinates, and approximate the spatial operators using the upwind SBP operators. Note that every 3D scalar grid function is rearranged row-wise as a vector of length n q n r n s . Therefore the unknown vector field Q is a vector of length 9n q n r n s .
Remark 2. The backward difference operator D − is used to approximate the spatial derivative for the conservative flux term, whilst the forward difference operator D + is an approximant for the non-conservative product term. This combination of upwind operators and the specific choice of the anti-symmetric form (11) is critical to deriving a conservative and energy stable scheme for the elastic wave equation in complex geometries.
Note that we have not imposed boundary conditions yet. Numerical boundary treatment will be discussed in the next subsection.
We will now prove the discrete equivalence of Lemma 2.1.

Lemma 4.2. Consider the semi-discrete approximation given in Equation 40
. We have the discrete antisymmetric form Proof. As before, expanding the matrix multiplication and simplifying yields the result.
Further, for a 3D scalar field f ijk = f (x i , y j , z k ) we also introduce the surface cubature and Therefore we have where the surface cubature S ξ is defined in (37)-(39). Here S ξ approximates integrals on faces of Ω, whilst Expanding the right hand side and applying the multi-dimensional SBP property (36) yields ξ∈{q,r,s} which from Lemma 4.2 gives the result.

Boundary treatment
We will now implement the boundary conditions (20) weakly using penalties. The idea is to impose the boundary conditions as Simultaneous Approximating Terms (SAT) in (40) with appropriate penalty parameters such that the numerical boundary terms do not permit energy growth. We will choose penalty parameters such that a discrete energy estimate is derived.
The semi-discrete approximation with weak enforcement of boundary conditions is where SAT ξ,i are penalty terms added to the discrete equation (40) at the boundaries to enforce the boundary conditions (20). The SAT penalty terms are not unique, they are designed such that the boundary procedure is consistent and the discrete approximation is energy stable. We will consider first the case of a free-surface boundary condition, and proceed later to the general case.

SAT term for the free-surface boundary condition
We consider specifically the free-surface boundary condition at all boundary surfaces, F ξ,0 , F ξ,1 for all ξ ∈ {q, r, s}. With the free-surface boundary condition, at F ξ,0 , F ξ,1 , the traction vector vanishes H s = I 9 ⊗ I nq ⊗ I nr ⊗ H s , e 0q = (I 9 ⊗ e 0q e T 0q ⊗ I nr ⊗ I ns ), e nq = (I 9 ⊗ e nq e T nq ⊗ I nr ⊗ I ns ), e 0r = (I 9 ⊗ I nq ⊗ e 0r e T 0r ⊗ I ns ), e nr = (I 9 ⊗ I nq ⊗ e nr e T nr ⊗ I ns ), e 0s = (I 9 ⊗ I nq ⊗ I nr ⊗ e 0s e T 0s ), e ns = (I 9 ⊗ I nq ⊗ I nr ⊗ e ns e T ns ), Here I 9 and I n ξ are identity matrices of size 9 × 9 and n ξ × n ξ , respectively, and e 0 ξ , e n ξ are boundary projection operators.
We state the first main result where ξ 0 = 0 and ξ n ξ = 1. We therefore have which completes the proof.

SAT terms for general boundary conditions
We will now construct the SAT terms for the general boundary condition (20). Similar to the DG framework [6,18,7], a weak boundary procedure can be derived by constructing boundary data, v η , T η , which are the solution of a Riemann-like problem constrained to satisfy the boundary condition (20) exactly.
SAT penalty terms are constructed by penalizing data, that is v η , T η , against the in-going waves only. The construction of boundary data, v η , T η , can be found in Appendix Appendix A. We also refer the reader to [6,18] for more detailed discussions.
Introduce the penalty terms The penalty terms are computed in the transformed coordinates l, m, n. We will now rotate them to the physical coordinates x, y, z, we have Note that We introduce the SAT vector that matches the eigen-structure of the elastic wave equation Here, n = (n x , n y , n z ) T is the unit normal vector on the boundary defined in (16). Note that Finally, the SAT terms for the general boundary conditions are defined as follows Introduce the fluctuation term and discrete boundary surface terms I v T T . Note that where the surface cubature S ξ is defined in (37)-(39). Note also that by (A.6) the boundary term is never positive, I v T T ≤ 0 for all |γ η | ≤ 1, and by (56) the fluctuation term is never positive, F luc (G, Z) ≤ 0.
We state the second main result Theorem 5.2. Consider the semi-discrete approximation (47) of the elastic wave equation with the SATterms SAT ξ,i defined in (55). We have with By using the identity (52) we have The proof is complete.
The fluctuation term F luc (G, Z) ≤ 0 adds a little numerical dissipation on the boundary. However, in the limit of mesh refinement the fluctuation term vanishes, that is F luc (G, Z) → 0 + as h → 0 + , and we have

Numerical error analysis
In this section we will analyse the numerical errors for the semi-discrete approximation (47). We will derive a priori error estimate and prove convergence of the error as h → 0 + . Next, we will discuss numerical dispersion errors that are peculiar to wave propagation problems, and which are most prominent at high frequencies.

A priori error estimate
Let Q denote the exact solution of the IBVP, and Q(q i , r j , s k , t) denote the restriction of the exact solution on the grid (q i , r j , s k ) at time t. We introduce the numerical error on the grid The error E ∈ R 9nqnrns satisfies the error equation where T ∈ R 9nqnrns is the truncation error of the SBP FD operator. Note that the truncation error T is a 3D vector field. However, it has a structure which is similar in all spatial directions. In particular, for grid points ξ j = jh ξ in the spatial direction ξ ∈ {q, r, s} the truncation error is of the form where T ijk = T q,i + T r,j + T s,k . Here, β j , are mesh independent constants, h ξ > 0 is the grid spacing, γ ∈ {1, 2, · · · } is the order of accuracy of the SBP FD stencils close to the boundary and ν ∈ {1, 2, · · · } is the order of accuracy of the SBP FD stencils in the interior, away from the boundaries. For traditional SBP operators based on central difference stencils the interior accuracy is always even, and we have (γ, ν) = (p, 2p), for p ∈ N. For upwind SBP operators the interior order of accuracy can be odd or even. As discussed in section 4, see also [12], upwind SBP operators with even-order (2p)-th accuracy in the interior are closed with p-th order accurate stencils close to boundaries, and we also have (γ, ν) = (p, 2p). Upwind SBP operators with odd-order (2p + 1)-th accuracy in the interior are closed with p-th order accurate stencils close to boundaries, giving (γ, ν) = (p, 2p + 1). The traditional SBP operators and upwind operators can yield (p + 1)-th global order of accuracy, for smooth solutions.  Note also that Theorem 6.1 hold for traditional and upwind SBP operators. Thus traditional and evenorder upwind SBP operators with (γ, ν) = (p, 2p) will have the same asymptotic error E = O(h p+1 ). However, as we will see from numerical experiments performed in the next section, on a marginally resolved mesh the numerical errors could completely differ, where for example the upwind SBP operators yield optimal numerical errors than the traditional SBP operator. This can be explained by analysing the numerical dispersion properties of the operators.

Numerical dispersion relation analysis
In order to understand the numerical dispersion properties of the SBP operators we consider the 1D shear plane wave propagating along the x-axis with the shear wave speed c s = µ/ρ > 0, where µ > 0 is the shear modulus and ρ > 0 is the density of the medium. Here, Q 0 = (v 0 , σ 0 ) T is the constant polarisation vector, ω ∈ R is the temporal frequency and k ∈ R is spatial wave number. The plane wave Q(x, t) solves the 1D elastic wave equation subject to the solvability condition called the dispersion relation ω = c s k.
We introduce the phase velocity V p and the group velocity V g defined by For the simple 1D model (63) the dispersion relation (64) is linear, the phase and group velocities are constant c s > 0, and are independent of the frequency and the wave number. Now let us consider the semi-discrete counterpart where D + and D − are the upwind SBP operators. We consider the interior stencils only, thus with the consistency requirements Here, h > 0 is the uniform grid spacing, α j , β j are the non-dimensional constant coefficients defining the upwind finite difference stencils. Note that β j = −α −j for j = −q, −(q − 1), · · · r. The Upwind SBP finite difference coefficients α j for the interior stencils with order of accuracy 2, 3, 4, 5, 6, 7, 8, 9 are given in Table   1.  For the traditional SBP finite difference operator we have where γ j are the non-dimensional constant coefficients of the finite difference operator, with γ −j = −γ j , γ 0 = 0, and satisfying the consistency requirements q j=−q γ j = 0, q j=1 2jγ j = 1.
The traditional SBP finite difference coefficients γ j for the interior stencils with order of accuracy 2, 4, 6, 8 are given in Table 2.   In general, a mode ω( k 1 ) is spurious if ∂ ω( k)/∂ k k= k1 < 0. The negative group velocity will propagate energy in the opposite (wrong) direction. This also implies that there are ( k 0 , k 1 ) with 0 ≤ k 0 < k 1 such that ω( k 0 ) = ω( k 1 ). Note that ω( k 0 ) would correspond to a resolved wave mode, and ω( k 1 ) is an unresolved spurious wave mode which can be poisonous to numerical simulations.
For simplicity we set the shear wave speed to c s = 1, and because of symmetry we plot the dispersion relations in the positive sub-interval k ∈ [0, π]. The numerical dispersion relations are displayed in Figure 2 for upwind operators with interior order of accuracy 2, 3, 4, 5, 6, 7, 8, 9 and for traditional centered difference SBP operators with interior order of accuracy 2, 4, 6, 8.
We summarise the significant observations below: 1. In general upwind SBP operators have better dispersion properties than traditional SBP operators based on centered finite difference stencils.
2. The properties of the numerical dispersion relation may improve with increasing accuracy.
3. The numerical dispersion relation reaches near optimal properties for 6th order accurate upwind SBP operator.
4. Beyond the 6th order accurate upwind SBP operator, higher order accuracy does not improve the dispersion properties of the upwind SBP operators.
5. even-order upwind SBP operators of order 2, 4, 6 do support not spurious unresolved wave mode, since ∂ ω( k)/∂ k 0 for all k. All odd-order upwind SBP operator including upwind SBP operator of order 8 support spurious unresolved high frequency wave modes, ∂ ω( k)/∂ k < 0 for some k.

For all traditional SBP operators almost half of the spectrum is dominated by spurious unresolved
wave modes with ∂ ω( k)/∂ k < 0.

Numerical simulations in 3D
In this section, we present numerical simulations in 3D. The experiments are designed to evaluate accuracy and demonstrate the efficiency of the upwind SBP operators in resolving scattered high frequency waves from complex geometries on marginally resolved meshes. We will also demonstrate parallel efficiency and perfect scaling of our parallel implementation and numerical simulation on Gadi 1 , Australia's newest supercomputer.
We will use WaveQLab3D's default implementation [4] of traditional SBP operator of interior order 6 and upwind operators [12] of interior order 2, 3, 4, 5, 6, 7, 8, 9. Then we will compare the numerical accuracy of traditional SBP operator of interior order 6 and upwind operators of interior order 6. These SBP operators are closed with 3rd order accurate stencils have 4th order global accuracy. The solutions are integrated in time using 4-th order acurate low-storage Runge-Kutta time stepping scheme [1].
To verify accuracy, we compute numerical solution of the 3D benchmark problem LOH1 [9,10], which has a semi-analytic-solution and compare results. We then show the potential of the upwind scheme by simulating a large scale 3D wave propagation problem with complex geometry at different grid resolutions.
For the 6-th order accurate operators, the upwind scheme provides comparable results on a marginally resolved grid to the traditional SBP operator on a finer grid, thereby improving computational efficiency for 3D numerical seismic wave simulations.
In all simulations, we have used the PML [5] to prevent artificial numerical reflections, from the computational boundaries, from contaminating the numerical simulations. A stable implementation of the PML for the 3D elastic wave equation IBVP using the upwind SBP operators is a non-trivial task. The details of the numerical treatment of the PML using upwind SBP operators will be reported in a forthcoming paper.

Layer over homogeneous half-space (LOH1) 3D benchmark problem
To verify and assess the numerical accuracy of our upwind method, we choose the 3D seismological benchmark problem, Layer Over Homogeneous Half-space (LOH1) [19,9,10] benchmark problem, a Seismic wave Propagation and Imaging in Complex media (SPICE) validation code. The LOH1 benchmark has a planar free surface and an internal interface between a thin low velocity (soft) upper-layer and high velocity (hard) lower crust, see Figure 3. The material properties for the soft upper-layer and hard lower-half-space  where M 0 = 10 18 Nm is the moment magnitude. Note that the moment tensor source (68) is spatially singular. Our numerical implementation approximates the singular source to high order accuracy as in [16].
In the z and y directions, the domain of the problem is unbounded. In the positive x direction (intowards the Earth), the domain is also unbounded with the Earth's surface x = 0 having the free surface, traction-free, boundary condition (T x , T y , T z ) = 0.
The SPICE code validation project [19] has suggested to use large enough computational model, namely  To evaluate the results of our schemes, simulated data is compared at the receiver, labelled station 6 and 9 in the LOH1 documentation [19], located on the face F q,0 at the locations (0, 7.348, 7.348) and  In Figures 5-6 we compare the analytical solution with specifically the numerical seismograms for the 6th order upwind SBP operator and the 6th order traditional SBP operator, at two mesh resolutions. In Figure 5, it can be seen that both the upwind and traditional schemes have relatively comparable accuracy for the benchmark at 100 m grid resolution. Both of these schemes approximate the analytical solution, with certain peaks in the data containing discrepancies between the numerical and analytical solutions, such as the peak at ≈ 6s in the V x data at station 6. We note however that the upwind SBP numerical seismogram matches the analytical solution better than the traditional SBP numerical seismogram in most parts of the solution.
Next we refine the mesh by doubling the number of grid points in each direction. That is, in the thin low

Large-scale numerical simulations in a 3D complex geometry
We will now present numerical simulations in a complex geometry, with a geologically constrained complex non-planar free-surface topography. Zugspitze is the tallest mountain in Germany, lying in the Wetterstein mountain range. The topography of this region is complex, with large variations in altitude across the Earths surface. We extracted the topography data of the Zugspitze region from the high resolution Alpine topography [2]. Accurate and stable simulation of seismic wave propagation in this region is a computationally expensive task, the main reason being the high frequency wave modes generated by scattering from the complex non-planar topography. Effective numerical simulations will require an efficient HPC and parallel numerical elastic wave solver that scales perfectly with increasing supercomputing resources. We will present first some scaling tests to demonstrate efficient parallel implementation of the upwind SBP FD operators for the numerical simulations of 3D elastic waves in complex geometries, and will proceed later to the numerical simulations of seismic waves in the Zugspitze region.
The Zugspitze setup. The modelling domain is Ω = y,z∈ [−5,85] [ X(x, y), 80]×[−5, 85] 2 with the x-co-ordinate being positive in-towards the Earth, like our previous example, and X(y, z) parameterising the Earth's surface. Our primary objective is to understand how different upwind SBP FD operators and the traditional SBP FD operator will resolve high frequency scattered wave modes generated by the complex free-surface topography without introducing spurious wave modes in the solution. To isolate scattering from complex geometries, the material parameters of the region is assumed to be constant and given by ρ = 2700 kg/m 3 , c p = 6000 m/s and c s = 3464 m/s. Therefore scattered wave fields seen in our simulations are primarily the effects of complex non-planar topography of the Zugspitze region. The location of the moment tensor point source is at (10,10,10), that is 10 km at depth. As per the LOH1 experiment, we have free surface boundary conditions at the complex topography X(y, z) and use the PML to prevent artificial reflections from the computational boundaries from contaminating the solution.

Scaling tests
For large scale numerical simulations of PDEs, it is imperative that the parallel numerical software is efficient and scalable. Here, we will perform strong scaling tests to verify the efficiency of our parallel implementation of high order upwind SBP operators for large-scale elastic wave simulations in 3D geometrically complex elastic solids. We will demonstrate strong perfect scaling with increasing supercomputing resources.
To conduct these scaling tests we consider the Zugspitze setup as described above with h = 200 m uniform resolution of the complex topography and the model volume. At this mesh resolution, the discritisation generates about 1 billion DoF for the evolving unknown vector field, and about 1.45 billion DoF needed to store the mesh, material parameters, the Jacobian and metric parameters. We run the simulation for 1 s using the same time step for all order of accuracy. In these runs, the I/O is disabled as to only evaluate the performance of the numerical solver and the time-to-solution.
The HPC implementation of WaveQLab is parallelised with MPI 2 . The time-to-solution can be split into the CPU-time used for floating point operations to evaluate the FD stencil and the MPI-time used to communicate ghost nodes shared between adjacent processors. The CPU-time will be proportional to the width of the finite difference stencil, (r + q + 1) for upwind operator and (2q + 1) for the traditional operator, and the MPI-time which is proportional to the number of ghost q, where r and q are listed in Tables 1-2   We will now formulate a simple theoretical model for the scaling tests and speed-up. Let T N denote the wall-clock time it takes to reach the final time for N number of compute nodes, and S N denotes the achieved speed-up. Ideally we have T N := T 1 N ⇐⇒ log 10 (T N ) = − log 10 (N ) + log 10 (T 1 ), where T 1 is the wall-clock time it takes to reach the final time on single node. Thus, the experimental results are modelled by where b 1 = log 10 (T 1 ), and the model parameters m 1 and m 2 are determined experimentally. Note that in the theoretically ideal case we have m 1 = −1 and m 2 = 1. The experimentally determined constant parameters m 1 and m 2 are given in Table 3 for each of the upwind operators.    We present the computational results in Figures 9, 10, 11, 12. Figure      on the Earths surface. Note that the numerical solutions converge to the reference solution as the order of accuracy increases. Again, for the marginally resolved mesh the accuracy of the solution becomes nearly optimal for the 6th order accurate upwind SBP operator. Beyond the 6th order accurate upwind SBP operator, the higher order accurate SBP operators do not necessarily improve the accuracy of numerical solution for the upwind SBP operators without introducing spurious oscillations. We also identify that the odd-order operators have more spurious high frequency modes in their numerical solutions than their evenorder operator counterparts. As above, this is consistent with the numerical dispersion relation analysis of Section 6.2 and numerical dispersion plots displayed in Figure 2. Much of the spurious oscillation in the traditional scheme has diminished with increasing mesh resolution, that is for h = 100 m grid spacing. However, doubling the mesh resolution by using h = 100 m grid spacing, increases the computational resources for the 3D problem by several orders of magnitude.   to see the coda waves. All odd-order (3,5,7,9) accurate SBP FD operators generate high frequency spurious modes.

Summary and outlook
High-order accurate SBP FD methods are attractive for efficient large-scale numerical simulation of hyperbolic wave propagation problems. Traditional SBP FD operators that approximate first-order spatial derivatives with central-difference stencils often have spurious unresolved wave-modes in their numerical solutions. On marginally resolved computational grids, these spurious wave-modes have the potential to destroy the accuracy of numerical solutions for a first-order hyperbolic partial differential equation, such as the elastic wave equation. Recently derived high order accurate upwind SBP operators [12] based non-central (upwind) FD stencils have the potential to suppress these poisonous spurious wave-modes on marginally resolved computational grids.
We have demonstrated in this paper that not all high order upwind SBP FD operators are applicable.
Numerical dispersion relation analysis shows that odd-order upwind SBP FD operators also support spurious unresolved high frequency wave modes on marginally resolved meshes. Even-order upwind SBP FD operators (of order 2, 4, 6) do not support spurious unresolved wave modes and have better numerical dispersion properties than traditional SBP FD operators and odd-order upwind SBP FD operators.
To ensure the accuracy of numerical solutions of the three space dimensional (3D) elastic wave equation in complex geometries, we discretise the 3D elastic wave equation with a pair of upwind SBP operators, on boundary-conforming curvilinear meshes. Using the energy method we prove that the numerical method is stable, and energy conserving. We derive a priori error estimates and prove the convergence of the numerical error.
We presented numerical simulations of the 3D elastic wave equation in heterogeneous media with complex non-planar free surface topography, including numerical simulations of community developed seismological benchmark problems. Our results show that even-order upwind SBP FD operators of order are more robust and less prone to numerical dispersion errors on marginally resolved meshes when compared to odd-order upwind SBP FD operators and traditional SBP FD operators, thereby increasing efficiency.
We performed strong scaling tests, demonstrating nearly perfect strong scaling and verifying the efficiency of our parallel implementation of high order upwind SBP operators for large scale elastic wave simulations in 3D complex geometries. We believe that the method and the software will increase the efficiency of numerical simulations of elastic waves in many applications such as earthquake engineering, natural minerals and energy resources exploration, as well in strong-ground motion analysis and underground fluid injection monitoring.
Our preliminary 3D dynamic earthquake rupture simulations show promise in increasing efficiency through using upwind SBP operators with good dispersion relation properties for simulating nonlinear friction laws in elastic solids, and earthquake source modelling.
We have used the PML to enable efficient domain truncation and prevent artificial numerical reflections, from the computational boundaries, from contaminating the numerical simulations. Because of the asymmetric properties of the PML and the upwind SBP operators, a stable implementation of the PML for the 3D IBVP elastic wave equation using the upwind SBP operators is a non-trivial task. The details of the numerical treatment of the PML using upwind SBP operators will be reported in a forthcoming paper.