Spectral-element simulations of elastic wave propagation in exploration and geotechnical applications

Abstract

We apply the spectral-element method (SEM), a high-order finite-element method (FEM) to simulate seismic wave propagation in complex media for exploration and geotechnical problems. The SEM accurately treats geometrical complexities through its flexible FEM mesh and accurately interpolates wavefields through high-order Lagrange polynomials. It has been a numerical solver used extensively in earthquake seismology. We demonstrate the applicability of SEM for selected 2D exploration and geotechnical velocity models with an open-source SEM software package SPECFEM2D. The first scenario involves a marine survey for a salt dome with the presence of major internal discontinuities, and the second example simulates seismic wave propagation for an open-pit mine with complex surface topography. Wavefield snapshots, synthetic seismograms, and peak particle velocity maps are presented to illustrate the promising use of SEM for industrial problems.

1 Introduction

Mapping the elastic and anelastic structures of the subsurface, for both global and industrial seismic applications, relies on the ability to solve the seismic wave equations accurately (e.g., Dahlen and Tromp 1998; Aki and Richards 2002). In particular, with the advancing of numerical methods and increasing computational power, full numerical calculations of synthetic seismograms have contributed significantly to the seismic inverse problem in seismological research since the early 80 (e.g., Tarantola 1984a, b; Chaljub et al. 2007; Virieux and Operto 2009; Liu and Gu 2012).

The finite-difference (FD) methods are the most widely used among all numerical approaches that model the propagation of seismic waves due to their relatively simple formulation. They employ the strong form of the differential equations of motion and essentially replace derivatives of velocity vectors and stress tensors by their finite difference equivalents between adjacent grid points (e.g., Madariaga 1976; Virieux 1988; Levander 1988; Komatitsch et al. 2005). Other classical gridding methods, such as spectral and pseudo-spectral methods, have also been applied to accurately solve seismic wave propagation in relatively smooth media (e.g., Tessmer et al. 1992; Carcione 1994; Komatitsch et al. 2005). However, these classical methods have inherent limitations for models that exhibit complex geometries, including undulated surface topography and major internal discontinuities, which although could be overcome, are numerically costly (e.g., Robertsson 1996; Ohminato and Chouet 1997; Komatitsch et al. 2005).

In comparison, the finite-element method (FEM), another category of numerical techniques in the space-time domain, honors the geometrical model complexity by sophisticated mesh design (Bathe 2008), and has been exploited to predict ground motions for complex media (e.g., Lysmer and Drake 1972; Bao et al. 1998; Bielak et al. 1999; Aagaard et al. 2001). The use of the unstructured meshing technique, which is based upon various types of elements (e.g., triangles, tetrahedra, and hexahedra), provides the FEM with the flexibility to honor strong surface topography and major internal discontinuities, as element boundaries coincide with these interfaces and discontinuities (Komatitsch et al. 2005; Tromp et al. 2008; Casarotti et al. 2008; Peter et al. 2011). Traditionally, low-degree polynomials are used in the FEM to represent shapes of elements as well as to express functions and integrals on the elements, which may be inadequate to capture large oscillations of wavefield across an element in seismic wave propagation problems. The use of higher-order elements in FEM can increase the accuracy of simulations although at the expense of significantly arising computation cost (Komatitsch et al. 2005).

The spectral-element method (SEM) is essential a higher-degree FEM, which combines the boundary composition flexibility of the FEM with the spectral accuracy of the pseudo-spectral method (Komatitsch and Tromp 1999; Komatitsch et al. 2004). Unlike the classic FEM, the use of high-degree polynomials defined over Gauss–Lobatto–Legendre (GLL) points to discretize functions on elements and as quadrature points for numerical integration, leads to a diagonal mass matrix, which can be easily inverted. It can be also efficiently implemented in parallel (e.g., Komatitsch and Tromp 1999; Komatitsch et al. 2005; Tromp et al. 2008) and makes it possible to solve problems of increasing scales. These favorable features have made SEM an important tool for seismic wave simulations in regional (e.g., Komatitsch and Tromp 1999; Komatitsch et al. 2004; Liu et al. 2004), as well as global (e.g., Komatitsch and Tromp 2002a, b; Chaljub and Valette 2004) earth models.

Recently, some efforts have also been made to employ SEM in industrial-scale problems (e.g., Zhu et al. 2009; Luo 2012; Gharti et al. 2012a, b). In this paper, we take advantage of an existing SEM software package SPECFEM2D (http://www.geodynamics.org) to study two cases of exploration and geotechnical problems, including a marine survey for a salt dome, and wave simulations for an open-pit mine scenario. We analyze wavefield snapshots, synthetic seismograms, and peak particle velocity (PPV) maps from these simulations to demonstrate the applicability of SEM in the forward modeling of wave propagation in complex media for industrial problems.

2 Brief introduction to the SEM

As the main focus of this paper is on the potential applications of SEM to specific types of industrial scenarios, we only briefly summarize the SEM technique below, with an emphasis on mesh design. More detailed comparisons to other numerical methods, such as FD and low-order FEM, can be found in Komatitsch and Tromp (1999), Komatitsch et al. (2004), Chaljub et al. (2007), and Peter et al. (2011).

2.1 Governing equations

Given an Earth model with volume $\mathbf{\Omega }$ and a free surface $\mathit{\partial }\mathbf{\Omega }$, the governing equations for elastodynamic problems can be generally written in the differential form

$$\mathit{\rho }{\mathit{\partial }}_t^2\mathit{s}=\mathrm{\nabla }·\mathit{T}+\mathit{f},$$

(1)

subject to the free-surface boundary condition

$$\widehat{\mathit{n}}·\mathit{T}=0\text{on}\mathit{\partial }\mathbf{\Omega },$$

(2)

wheres(x, t) denotes the displacement motion at position x, ρ(x) the distribution of mass density, and f(x, t) the external body force. The stress tensor T is related to the displacement gradient ∇s through the constitutive relation

$$\mathit{T}=\mathit{c}:\mathrm{\nabla }\mathit{s},$$

(3)

where c is a fourth-order elastic tensor. In addition, the displacement s and traction $\widehat{\mathit{n}}·\mathit{T}$ should be continuous across solid-solid interfaces. If both solid and fluid regions are present in the earth model, the SEM formulates the elastodynamic problem in terms of velocity potential in the fluid and displacement in solid regions, resulting in separate governing equations that are coupled at fluid-solid boundaries through continuity of traction and the normal component of velocity (Komatitsch et al. 2000a; Komatitsch and Tromp 2002a, b; Chaljub and Valette 2004; Tromp et al. 2008).

The integral or weak form of the elastodynamic equation, frequently used in both FEM and SEM, is obtained by dotting the differential equation (1)—the strong form for elastostatic problems—with an arbitrary test vector function w(x), and integrating over the entire model (Komatitsch and Tromp 1999; Komatitsch et al. 2005; Tromp et al. 2008)

$$\underset{\mathrm{\Omega }}{\int }\mathit{\rho }\mathit{w}·{\mathit{\partial }}_t^2\mathit{s}{(\text{d})}^3\mathit{x}=-\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }\mathit{w}:\mathit{T}{\text{d}}^3\mathit{x}+\underset{\mathrm{\Omega }}{\int }\mathit{f}·\mathit{w}{\text{d}}^3\mathit{x}.$$

(4)

Since Eq. (4) is valid for any test vector w, the weak and strong forms are mathematically equivalent. Note that the free-traction surface boundary condition (2) is satisfied naturally through integration by parts, so that the incorporation of surface topography is straightforward. Therefore, surface waves can be simulated more accurately by SEM than those methods based upon the strong formulation (Komatitsch and Tromp 1999). For regional problems, the differential equation (1) needs to be solved subject to absorbing boundary condition in which seismic waves travelling out of the volume $\mathbf{\Omega }$ are absorbed on fictitious boundaries and do not re-enter the domain (Komatitsch and Tromp 1999; Tromp et al. 2008). For the weak formulation, the absorbing boundary condition can be naturally imposed on Eq. (4) in terms of traction over the fictitious boundaries.

2.2 The spectral-element discretization

Similar to FEM, the first step of discretizing the governing equations is to subdivide the model volume $\mathbf{\Omega }$ into a number of non-overlapping elements (e.g., Komatitsch and Tromp 1999; Komatitsch et al. 2005; Tromp et al. 2008; Casarotti et al. 2008). Although a variety of elements such as tetrahedra, hexahedra, and prisms may be used in the FEM, the SEM is generally restricted to use meshes of unstructured hexahedral elements (deformed cube) in 3D case or quadrilateral elements (deformed square) in 2D case. Each hexahedral element can be mapped to a reference cube which is defined in terms of natural coordinates $\mathit{\xi }=\left(\mathit{\xi },\mathit{\eta },\mathit{\zeta }\right),-1\le \mathit{\xi },\mathit{\eta },\mathit{\zeta }\le 1$ using the invertible transformation.

$$\mathit{x}\left(\mathit{\xi }\right)=\sum _{a=1}^{n_a}{\mathbf{x}}_a{N}_a\left(\mathit{\xi }\right),$$

(5)

where ${\mathbf{x}}_a\left({\mathit{\xi }}_a,{\mathit{\eta }}_a,{\mathit{\zeta }}_a\right),a=1,\dots ,n_a$ denote n _a control points within a hexahedron and $N_a(\mathit{\xi })$ corresponding shape functions. The hexahedral shape functions N _a are triple products of one-dimensional Lagrange polynomials of degree 1 or 2. Generally, the n + 1 Lagrange polynomials of degree n are defined in terms of n + 1 control points $-1\le {\mathit{\xi }}_{\mathit{\alpha }}\le 1,\mathit{\alpha }=0,\dots ,n$, by

$${\ell }_{\mathit{\alpha }}=\frac{(\mathit{\xi }-{\mathit{\xi }}_0)\cdots (\mathit{\xi }-{\mathit{\xi }}_{\mathit{\alpha }-1})(\mathit{\xi }-{\mathit{\xi }}_{\mathit{\alpha }+1})\cdots (\mathit{\xi }-{\mathit{\xi }}_n)}{({\mathit{\xi }}_{\mathit{\alpha }}-{\mathit{\xi }}_0)\cdots ({\mathit{\xi }}_{\mathit{\alpha }}-{\mathit{\xi }}_{\mathit{\alpha }-1})({\mathit{\xi }}_{\mathit{\alpha }}-{\mathit{\xi }}_{\mathit{\alpha }+1})\cdots ({\mathit{\xi }}_{\mathit{\alpha }}-{\mathit{\xi }}_n)}.$$

(6)

To ensure the mapping from the element to the reference cube is unique and invertible, the transformation (5) should have nonzero Jacobian. Similarly, in 2D, each quadrilateral element is mapped to a reference square.

In the SEM, a function (or field) f is interpolated over a hexahedral element using the Lagrange polynomials of degree n defined over the Gauss-Lobatto-Legendre (GLL) points as (Komatitsch and Tromp 1999; Komatitsch et al. 2005; Tromp et al. 2008)

$$f\left(\mathit{x}\left(\mathit{\xi },\mathit{\eta },\mathit{\zeta }\right)\right)=\sum _{\mathit{\alpha },\mathit{\beta },\mathit{\gamma }=0}^n{f}^{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}{\ell }_{\mathit{\alpha }}\left(\mathit{\xi }\right){\ell }_{\mathit{\beta }}\left(\mathit{\eta }\right){\ell }_{\mathit{\gamma }}\left(\mathit{\zeta }\right),$$

(7)

where $f^{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}=f\left(\mathit{x}\left({\mathit{\xi }}_{\mathit{\alpha }},{\mathit{\eta }}_{\mathit{\beta }},{\mathit{\zeta }}_{\mathit{\gamma }}\right)\right)$ denotes the functional value at the GLL point $\mathbf{x}\left({\mathit{\xi }}_{\mathit{\alpha }},{\mathit{\eta }}_{\mathit{\beta }},{\mathit{\zeta }}_{\mathit{\gamma }}\right)$ and ℓ _α,β,γ the Lagrange polynomials. The n + 1 GLL points in each dimension are determined by the roots of the equation $\left(1-{\mathit{\xi }}^2\right)P_n^{\prime }\left(\mathit{\xi }\right)=0$, where $P_n^{\prime }$ denotes the first derivative of the Legendre polynomial of degree n. SEM simulations typically use polynomials of degree between 4 to 10 to represent functions on the elements, and degree 4 or 5 provides the best trade-off between accuracy and computational cost (Seriani and Priolo 1994; Tromp et al. 2008; Casarotti et al. 2008). In comparison, the classic FEM uses lower-degree polynomials (1 or 2 in general) as basis functions for both the representation of functions and the shape of elements.

To evaluate the integrals at the element level, the SEM utilizes the Gauss–Labotto–Legendre integration rule in conjunction with GLL interpolation points, instead of the Gauss quadrature frequently used in the classical FEM (Komatitsch and Tromp 1999; Komatitsch et al. 2005; Tromp et al. 2008; Casarotti et al. 2008). Although slightly less accurate, the use of the GLL quadrature significantly reduces the complexity and cost of the algorithm, and the mass matrix resulted from the evaluation of the left-hand term in the weak form (4) becomes diagonal (Komatitsch and Tromp 1999). Thus, the SEM can be efficiently run on parallel computers (Komatitsch and Tromp 1999).

2.3 Mesh design for the SEM

Generating a high-quality unstructured mesh appropriate for the model of interest, especially in 3D, is the first as well as one of the critical steps of both the FEM and SEM wave modeling processes (Casarotti et al. 2008). The fact that efficient 3D SEM simulations generally require unstructured all-hexahedral mesh significantly complicates the meshing process (Casarotti et al. 2008). As mentioned previously, an unstructured hexahedral mesh not only provides flexibility for the incorporation of topographic variations and major internal discontinuities (Komatitsch et al. 2005), but also leads to a diagonal mass matrix in conjunction with the use of the GLL integration points and quadrature. The effort of 2D/3D hexahedral mesh design is also justified by the fact that once a proper mesh is built for a region, wave simulations can be run repeated without any alteration to the mesh.

In order to ensure accuracy and computational stability, the unstructured mesh should be constructed so that the Jacobian of the transformation (5) varies smoothly throughout the model, and the distortion of elements is acceptable (Komatitsch et al. 2005; Casarotti et al. 2008). Generally, the equiangle skewness of the deformed hexahedral elements should be within 0.75 (0.8 is acceptable in some cases) to provide acceptable simulation accuracy (Komatitsch et al. 2005; Casarotti et al. 2008). To resolve a shortest period T₀ of interest, the choice of the element size $\mathrm{\Delta }h$ needs to be subject to the constraint

$$\mathrm{\Delta }h<\frac{n+1}{f(n)}{\mathit{\upsilon }}_{\min }{T}_0,$$

(8)

where υ_min denotes the minimum seismic wave speed in the model, and f(n) is an empirical function related to the degree n of polynomials used to interpolate functions on the elements (Casarotti et al. 2008). For example, if n is chosen as 4, the empirical function f(n) is equal to 5, ensuring that each wavelength λ contains at least 5 GLL points (Casarotti et al. 2008). The incremental time stepping duration is in turn constrained by the element size and the maximum seismic wave speed υ_max in the element by

$$\mathrm{\Delta }t<C\frac{\mathrm{\Delta }h}{{\mathit{\upsilon }}_{\max }},$$

(9)

where C is the Courant constant and can be chosen between 0.3 and 0.4 in the case of a deformed and heterogeneous mesh (Komatitsch et al. 2005; Casarotti et al. 2008). Moreover, the major discontinuities in wave speeds within the model need to be honored by the mesh to avoid numerical diffractions by a staircase sampling of the interfaces (e.g., Zahradník et al. 1993; Komatitsch et al. 2004; Casarotti et al. 2008).

3 SEM simulation examples

In this paper, a fully-developed 2D SEM community software package SPECFEM2D (downloaded from http://www.geodynamics.org) is applied to simulate seismic wave propagation in complex models that include both solid and fluid media, undulation of major discontinuities and strong lateral heterogeneities. In SPECFEM2D, polynomials of degree 4 associated with 5 GLL control points are used as basis functions to interpolate functions on each element. The SPECFEM2D package is an open-source MPI-based package that can be readily used in parallel. Another advantage of the SPECFEM2D package is that it supports not only the built-in mesher which is suitable for models with simple and regular geometries, but also mesh built by external quadrilateral/hexahedral meshing softwares such as CUBIT (http://cubit.sandia.gov), which provide more flexibility to the creation of high-quality mesh for complex media. All the quadrilateral meshes used in this paper are created with CUBIT. The SEM simulation examples presented are restricted to isotropic elastic and/or acoustic media. However, features such as full anisotropy, poroelasticity and attenuation are also implemented in the SPECFEM2D package and may be incorporated for future studies.

3.1 Simulation for exploration application

In this section, two simulations for a marine survey scenario are presented and analyzed based upon wavefield snapshots and seismograms, to illustrate the capability of the SEM, along with the use of unstructured quadrilateral mesh, in accommodating fluid layers (ocean), fluid–solid interface topography (bathymetry), and abrupt changes in the seismic velocity model. Both simulations utilize a similar geologic models representing a high-velocity salt dome beneath ocean sediments in a computation domain of 4,000 × 2,400 m (Fig. 1), and the only difference is that one includes a flat sea floor and the other one is simulated with an uneven sea floor (Fig. 1). In both models, The first layer (shown in red in Fig. 1) is a water layer with a thickness of approximately 600 m, and the sedimentary layer (shown in blue) extends to ∼2,400 m beneath the sea floor. A salt dome (shown in yellow) is placed in the sediments with significantly larger P- and S-velocity than the surrounding. The elastic properties of the different materials are listed in Table 1.

Fig. 1

Snapshots of velocity vectors at time t = 0.75 s on top of the unstructured quadrilateral meshes for the wave propagation simulations in Sect. 1. The size of modeling domains is 4,000 × 2,400 m for both models. The ocean layer with an average thickness of 600 m, the surrounding sediments beneath the ocean floor, and the salt dome are shown in red, blue, and yellow. The arrow indicates the position of the pressure source placed at 2 m below the ocean surface, diamonds indicate the position of the 101 receivers used to record seismograms as presented in Fig. 2, and the white triangle in the sediments indicates the position of the receiver used to record seismic wave propagation in the elastic regions as shown in Fig. 3. One can clearly observe the direct (a) and sea-floor reflected (b) P-waves in the ocean layer; the transmitted P (c) and P-to-S (d) converted waves through sea floor, as well as the reflected P (e) and S (f) waves off the top of the salt dome in the sediments. Phase (g) is the refracted P waves along the sea floor

Table 1 Elastic properties of models used in wave propagation simulations in Sect. 1

Taking into account the velocity variations, element sizes need to vary for different blocks to maintain on average a consistent number of grid (e.g., ≥5) points per S-wavelength for the simulations. The average element sizes for the fluid layer, sedimentary layer, and salt dome, are 12, 18, and 20 m. The CUBIT software automatically generates conforming unstructured quadrilateral meshes that accommodate the transition in element size between different materials (Fig. 1). Absorbing boundary conditions (Stacey 1988) are applied to the fictitious side and bottom boundaries of the model domains to prevent reflections. To investigate the impact of sea floor topography on seismic wave propagation, the same source mechanism and source-receiver geometry are used in both simulations. A Ricker pressure source (indicated by vertical arrows at [2,000, −2 m] in Fig. 1) with a dominant frequency of 20 Hz is injected 2 m below the ocean surface, to mimic pulses generated by airguns in marine surveys. A total number of 101 receivers are also evenly placed along a horizontal line 2 meters below the ocean surface with a source offset varying between −1,500 and 1,500 m to synthesize the seismic recordings of a string of hydrophones. Another receiver (at [2200, −800 m]) located in the surrounding sediment is used to record seismic waves motions in the elastic regions. Time steps are estimated for each block based on material wavespeed and mesh size according to Eq. (9), and the final time step employed by the simulation, $\mathrm{\Delta }t=1\times {10}^{-4}\mathrm{s},$ is chosen to be the smallest of all three estimates. A total duration of 1.5 s is simulated.

The source–receiver geometry as well as snapshots of the velocity vector field at time t = 0.75 s for both models are shown in Fig. 1. The corresponding horizontal- and vertical-component velocity seismogram profiles recorded by the horizontal receiver line are illustrated in Fig. 2. The two-component velocity seismograms recorded by the specific receiver located in the sediment for the models with flat and varying sea-floor topography are shown as the solid and dashed lines in Fig. 3. Several phases can be identified in all three figures, including (a) direct P waves, (b) reflected P waves from seabed, (c) transmitted P waves through seabed, (d) P-to-S converted waves across seabed, (e) reflected P waves of (c) off the top of salt dome, (e′) transmitted waves of (e) through seabed, (f) reflected S waves of (c) off the top of salt dome, (f′) transmitted waves of (f) converted to P waves in the ocean layer, and (g) refracted P waves along seabed. The presence of sea floor topography clearly affects the traveltime and amplitude of various phases, and is nicely captured by the SEM simulations. Both the snapshots and synthetic seismograms illustrate the ability of SEM method in simulating wave propagation in models with both fluid and solid regions as well as sharp velocity contrasts.

Fig. 2

Synthetic seismograms of the horizontal (top) and vertical (bottom) component velocity field recorded by the 101 receives denoted by the diamonds in Fig. 1 for the model with a flat ocean floor (a) and varying sea floor topography (b). Data traces are scaled by the maximum vertical velocity values in each scenario. The amplitudes of horizontal records are blown-up by a factor of 10 to better show the much weaker phases. The strong velocity contrast between the water and sedimentary layers results in the large amplitude of the reflected P wave (phase b). Phase (e′) and (f′) are the transmitted P waves in the ocean layer of phase (e) and (f), which are the P and S reflected waves off the top of salt dome

Fig. 3

Synthetic seismograms of velocity field recorded by the receiver in the sediment as denoted by the white triangle in Fig. 1. The dashed and solid lines correspond to the synthetics for model with and without sea floor topography. Both traces are scaled by a constant to allow for their proper display. The impact of the sea floor topography on seismic wave propagation is clearly demonstrated by the difference in arrival times of various phases, such as the transmitted P waves (phase c), the P-to-S converted wave at sea floor (phase d), the reflected P (phase e), and S (phase f) waves off the top of salt dome

3.2 Simulation for geotechnical application

In this section, we simulate wave propagation in three homogeneous 2D models with increasing free-surface complexity for an open-pit mine scenario (Fig. 4). Snapshots and synthetic seismograms, as well as the PPV maps for first breaks, are produced to illustrate the ability of the SEM in simulating seismic wave propagation in models that include interfaces with large topographic variations.

Fig. 4

(a) Geometry of the models described in Sect. 2. The size of all models is approximately equal to 2,216 × 1,226.98 m. The maximum widths and depths of both the bowl-shaped and staircase topographical regions are 1,616 and 426.78 m, respectively. A Ricker point force (denoted by the black star) with a dominant frequency of 100 Hz is excited vertically at (1108, −676.78 m). The recordings of 82 ground stations, located unevenly up to 1,010 m offset away from the source, are plotted in Fig. 6. (b) Bench geometry of the staircase model. (c) Unstructured quadrilateral mesh for the portion of the staircase model shown in (b)

The first model consists of a flat free-surface and is referred to as the flat-surface model; the second model is referred to as bowl-shape model because of the bowl-shaped surface depression; and the last one is inspired by section views of realistic open-pit mine models with staircase free-surface and referred to as the staircase model. In the staircase model, the bench width and height are 15 and 19.79 m, respectively, and the face angle is 56°. The geotechnical ramp width is 40 m and the depth of the pit is 426.78 m. The simulation domains for all models are approximately equal to 2,216 × 1,226.98 m, and the width of both the bowl-shape and the staircase topographical regions are 1,616 m.

To compare the effect of surface topography, the same seismic source mechanism and source-receiver layout, as well as elastic properties are used for all three simulations. The seismic source is given vertically at 676.98 m depth (i.e., 250 m below the bottom of the pit) as a Ricker point source with a dominant frequency of 100 Hz. In order to capture effects of topography on seismic wave propagation and energy radiation, 82 receivers are unevenly located along each ground surface up to a horizontal offset of 1,010 m to the left and right of the source. The homogeneous P- and S-wave velocities are 5,200 and 3,200 m/s, and mass density is 2,800 kg/m³. According to Eq. (8), to resolve waves with a dominant frequency of 100 Hz, the element size should be smaller than 32 m when 5 GLL points are used. Nonetheless, an average element size of 10 m is used to discretize the models so that each bench of the staircase model can be subdivided into two elements (Fig. 4c). The same time stepping $\mathrm{\Delta }t=5.0\times {10}^{-5}$ is used for all the simulations, and a total duration of 0.5 s is simulated. Stacey absorbing boundary condition are applied to all the simulations to prevent reflections from the artificial boundaries.

The displacement snapshot of seismic wave propagation in the three models are shown in Fig. 5 at t = 0.05 s, t = 0.15 s, and t = 0.25 s. For the homogeneous models, as wavefronts expand from the source, P and S waves clearly separate and show different radiation patterns (t = 0.15 s in Fig. 5a). Both waves are later reflected by the flat surface and travel downward (t = 0.25 s), including P-to-P, P-to-S, S-to-S, and S-to-P reflected and converted waves. In the model with bowl-shaped surface, P and S surface reflections occur much earlier (t = 0.15 s) due to the reduction of topography above the source. In comparison, in the staircase model, strong scattering effects is observed in addition due to abrupt variations of the ground surface.

Fig. 5

Snapshots of vertical displacement fields at t = 0.05, 0.15, and 0.25 s for simulations in Sect. 2. Grey color indicates the modeling domains. Red and blue colors represent positive and negative wavefield values

Figure 6 illustrates the vertical-component velocity seismograms recorded by the ground stations (shown in black), along with the derived PPVs for first breaks (in red). The scattered waves off the abrupt changes in surface topography in the staircase model are clearly shown in the synthetic seismograms (Fig. 6c), thanks to the power of SEM in accommodating strong free-surface topographic variations. Some artificial boundary reflected waves can be seen in the seismograms [denoted as (r) in Fig. 6], as Stacey condition works less effectively in this case probably due to the proximity of left and right fictitious boundaries to the surface regions with strong topography (i.e., edge of the model is only 98 m away from the nearest receiver), although amplitudes of these artificial reflections are much smaller than those of the direct P and S waves. The models used in these simulations are fully elastic with no attenuation. Wave amplitudes generally decrease with increasing distance from the source due to geometric spreading effect, as shown by the PPV plots, where smaller PPV value is generally observed for station further away from the source. However, due to the impact of large surface topographic undulation, lateral variations of PPVs can be observed, especially in the staircase simulation (Fig. 6c). These PPV values provide critical input to the analysis of slope stability for mining operations (Dowding 1985).

Fig. 6

Synthetic seismograms of the vertical velocity fields (black) recorded by the 82 receivers located along the ground surface and derived peak particle velocities (PPVs) for first breaks (red). The lateral variations of PPVs are due to the reflected and scattered waves caused by strong surface topography. The reflected waves (indicated by r) from the left and right absorbing boundaries can be observed, but their amplitudes are generally small compared with other phases

4 Conclusions

In this paper, we presented two simulation examples for exploration and geotechnical problems based upon the SEM. In Sect. 1, complex models that include the fluid layer (ocean), fluid–solid interface (sea floor) topography, and the asymmetric salt dome are easily accommodated by the SEM because of the use of the unstructured mesh generated by CUBIT as shown in Fig. 1. The snapshots of the velocity vector field and synthetic seismograms along with the distinct phases of the subsalt reflections illustrate that the SEM can accurately model seismic wave propagation in media with both fluid and solid regions as well as interface variations. The influence of the sea floor topography is demonstrated by differences in arrival times and amplitudes of various phases (see Fig. 3).

For geotechnical applications, Sect. 2 presents and analyzes SEM simulation results for an open-pit mine scenario, including highly irregular surfaces representing benches. With the simulations in the flat-surface, bowl-shape and staircase models, the capability of the SEM in handling strong surface topography is again proven by the displacement snapshots (Fig. 5) and synthetic seismograms (Fig. 6) recorded by ground stations.

Therefore, the SEM can be naturally applied to exploration and geotechnical applications as its weak formulation and flexible mesh design make it possible to deal with model complexities. It shows great potential in simulating seismic wave propagation more accurately and effectively for some industrial problems compared to other numerical methods, such as the FD method, the pseudo-spectral element method, and the traditional low-order FEM.