The spectral-infinite-element method combines the highly accurate spectral-element method with the mapped-infinite element method to solve the unbounded problems efficiently and accurately. Spectral elements represent the domain of interest, and a single layer of infinite elements captures outer space. To solve the weak form of the unbounded equation, we use Gauss-Legendre-Lobatto quadrature in spectral elements inside the domain of interest. Outside the domain, we use Gauss-Radau quadrature in the infinite direction, and Gauss-Legendre-Lobatto quadrature in the other directions.

The software package SPECFEM-X based on the SIEM aims to solve most (global) quasistatic problems.

Although earthquake-induced gravity perturbations are frequently observed, numerical modeling of this phenomenon has remained a challenge. Due to the lack of reliable and versatile numerical tools, induced-gravity data have not been fully exploited to constrain earthquake source parameters. From a numerical perspective, the main challenge stems from the unbounded Poisson/Laplace equation that governs gravity perturbations. Additionally, the Poisson/Laplace equation must be coupled with the equation of conservation of linear momentum that governs particle displacement in the solid. Most existing methods either solve the coupled equations in a fully spherical harmonic representation, which requires models to be (nearly) spherically symmetric, or they solve the Poisson/Laplace equation in the spherical harmonics domain and the momentum equation in a discretized domain, a strategy that compromises accuracy and efficiency. We present a spectral-infinite-element approach that combines the highly accurate and efficient spectral-element method with a mapped-infinite-element method capable of mimicking an infinite domain without adding significant memory or computational costs. We solve the complete coupled momentum-gravitational equations in a fully discretized domain, enabling us to accommodate complex realistic models without compromising accuracy or efficiency. We present several coseismic and post-earthquake examples and benchmark the coseismic examples against the Okubo analytical solutions. Finally, we consider gravity perturbations induced by the 1994 Northridge earthquake in a 3D model of Southern California. The examples show that our method is very accurate and efficient, and that it is stable for post-earthquake simulations.

Most earthquake location methods require phase identification and arrival-time measurements. These methods are generally fast and efficient but not always applicable to microearthquake data with low signal-to-noise ratios because the phase identification might be very difficult. The migration-based source location methods, which do not require an explicit phase identification, are often more suitable for such noisy data. Whereas some existing migration-based methods are computationally intensive, others are limited to a certain type of data or make use of only a particular phase of the signal. We have developed a migration-based source location method especially applicable to data with relatively low signal-to-noise ratios. We projected seismograms on to the ray coordinate system for each potential source-receiver configuration and subsequently computed their envelopes. The envelopes were time shifted according to synthetic P- and S-wave arrival times (computed using an eikonal solver) and stacked for a predefined time window centered on the arrival time of the corresponding phase. This was done for each component and phase individually, and the squared sum of the stacks was defined as the objective function. We applied a robust global optimization routine called differential evolution to maximize the objective function and thereby locate the seismic event. Our source location method provides a complete algorithm with only a few control parameters, making it suitable for automatic processing. We applied this method to single and multicomponent data using P and/or S phases. We conducted controlled tests using synthetic seismograms contaminated with a minimum of 30% white noise. The synthetic data were computed for a complex and heterogeneous model of the Pyhäsalmi ore mine in Finland. We also successfully applied the method to real seismic data recorded with the in-mine seismic network of the Pyhäsalmi mine.

Accurate and efficient simulations of coseismic and post-earthquake deformation are important for proper inferences of earthquake source parameters and subsurface structure. These simulations are often performed using a truncated halfspace model with approximate boundary conditions. The use of such boundary conditions introduces inaccuracies unless a sufficiently large model is used, which greatly increases the computational cost. To solve this problem, we develop a new approach by combining the spectral-element method with the mapped infinite-element method. In this approach, we still use a truncated model domain, but add a single outer layer of infinite elements. While the spectral elements capture the domain, the infinite elements capture the far-field boundary conditions. The additional computational cost due to the extra layer of infinite elements is insignificant. Numerical integration is performed via Gauss-Legendre-Lobatto and Gauss-Radau quadrature in the spectral and infinite elements, respectively. We implement an equivalent moment-density tensor approach and a split-node approach for the earthquake source, and discuss the advantages of each method. For post-earthquake deformation, we implement a general Maxwell rheology using a second-order accurate and unconditionally stable recurrence algorithm. We benchmark our results with the Okada analytical solutions for coseismic deformation, and with the Savage & Prescott analytical solution and the PyLith finite-element code for post-earthquake deformation.

Gravity anomalies induced by density heterogeneities are governed by Poisson’s equation. Most existing methods for modelling such anomalies rely on its integral solution. In this approach, for each observation point, an integral over the entire density distribution needs to be carried out, and the computational cost is proportional to the number of observation points. Frequently, such methods are sensitive to high density contrasts due to inaccurate resolution of the volume integral. We introduce a new approach which directly solves a discretized form of the Poisson/Laplace equation. The main challenge in our approach involves the unbounded nature of the problem, because the potential exists in all of space. To circumvent this challenge, we combine a mapped infinite-element approach with a spectral-element method. Spectral elements represent the domain of interest, and a single layer of infinite elements captures outer space. To solve the weak form of the Poisson/Laplace equation, we use Gauss–Legendre–Lobatto (GLL) quadrature in spectral elements inside the domain of interest. Outside the domain, we use Gauss–Radau quadrature in the infinite direction, and GLL quadrature in the other directions. We illustrate the efficiency and accuracy of our method by comparing calculated gravity anomalies for various density heterogeneities with corresponding analytical solutions. Finally, we consider a complex 3-D model of an ore mine, which consists of both positive and negative density anomalies.

Magnetic anomalies induced by a magnetization distribution are governed by Poisson’s equation. Analogous to gravity anomalies, magnetic anomalies are frequently computed using a direct integral approach, involving an integral over the entire magnetization distribution for each observation point. Alternatively, magnetic anomalies can be computed using the so-called Poisson’s relation, in which the magnetic scalar potential is approximated by the derivative of the gravitational potential based on the analogy between the gravity and magnetic fields. However, this analogy is generally limited to simple isolated objects, and computing the magnetic field involves the second derivative of the gravitational potential. Both methods suffer from an often inaccurate volume integral over complex objects, and the computational cost scales with the number of observation points. We implement a spectral-infinite-element method to directly solve a discretized form of the Poisson/Laplace equation for magnetic anomalies. The spectral-infinite-element method combines the highly accurate spectral-element method with the mapped-infinite element method, which reproduces an unbounded domain accurately and efficiently. This combination is made possible by coupling Gauss-Legendre-Lobatto quadrature in spectral-elements with Gauss-Radau quadrature in infinite elements along the infinite directions. We illustrate the efficiency and accuracy of our method by comparing calculated magnetic anomalies for various magnetization heterogeneities with corresponding analytical and commonly used computational solutions. Finally, we consider a practical example involving a complex 3D model of an ore mine.

We implement a 3D spectral-element method for multistage excavation problems. To simulate excavation in elastoplastic soils, we employ a Mohr–Coulomb yield criterion using an initial strain method. We parallelize the software based on non-overlapping domain decomposition using MPI. We verify the uniqueness principle for multistage excavation in linear elastic materials. We validate our serial and parallel programs, and illustrate several examples of multistage excavation in elastoplastic materials. Finally, we apply our software to a model of the Pyhäsalmi ore mine in Finland. Strong-scaling performance tests involving multistage excavation show that the parallel program performs reasonably well for large-scale problems.

We implement a spectral‐element method for 3D time‐independent elastoplastic problems in geomechanics. As a first application, we use the method for slope stability analysis ranging from small to large scales. The implementation employs an element‐by‐element preconditioned conjugate‐gradient solver for efficient storage. The program accommodates material heterogeneity and complex topography. Either simple or complex water table profiles may be used to assess effects of hydrostatic pressure. Both surface loading and pseudostatic seismic loading are implemented. For elastoplastic behavior of slopes to be simulated, a Mohr–Coulomb yield criterion is employed using an initial strain method (i.e., a viscoplastic algorithm). For large‐scale problems, the software is parallelized on the basis of domain decomposition using Message Passing Interface. Strong‐scaling measurements demonstrate that the parallelized software performs efficiently. We validate our spectral‐element results against several other methods and apply the technique to simulate failure of an earthen embankment and a mountain slope.