SEM3D: A 3D High-Fidelity Numerical Earthquake Simulator for Broadband (0–10 Hz) Seismic Response Prediction at a Regional Scale

Abstract

In this paper, we present SEM3D: a 3D high-fidelity numerical earthquake simulator that is tailored to predict the seismic wave field of complex earthquake scenarios from the fault to the epicenter site. SEM3D solves the wave-propagation problem by means of the spectral element method (SEM). The presented demonstrative test case was a blind $M_W$6.0 earthquake scenario at the European experimental site located in the sedimentary basin of Argostoli on the island of Kefalonia (Western Greece). A well-constrained geological model, obtained via geophysical inversion studies, and seismological model, given the large database of seismic traces recorded by the newly installed ARGONET network, of the site were considered. The domain of interest covered a region of 44 km × 44 km × 63 km, with the smallest grid size of 130 m × 130 m × 35 m. This allowed us to simulate the ground shaking in its entirety, from the seismic source to the epicenter site within a 0–10 Hz frequency band. Owing to the pseudo-spectral nature of the numerical method and given the high polynomial order (i.e., degree nine), the model featured 1.35·10${}^{10}$ DOFs (degrees of freedom). The variability of the synthetic wave field generated within the basin is assessed herein, exploring different random realizations of the mean velocity structure and heterogeneous rupture path.

1. Introduction

The physics of an earthquake are naturally multi-scale. For instance, crustal earthquakes nucleate from seismogenic extended faults whose characteristic lengths range around ≈30 km for a magnitude of $M_W$6.0. On the other hand, the Earth’s crust presents layers that can range from a few kilometers to a few meters in thickness. Moreover, geomaterials show a small sub-meter scale heterogeneity. Seismic waves originate from active faults and travel toward the Earth’s surface through layers of heterogeneous, anisotropic rock and sedimentary or alluvial soil. The latter may display highly non-linear hysteretic behavior [1]. The recorded ground motion at the surface conserves the source footprint, along with the path and site effects of the earthquake, although the energy accumulated at fault discontinuities and released into the Earth’s crust is partially dispersed and absorbed while traveling toward the surface [2].

Due to this, the assessment of the seismic responses of critical structures and infrastructures during extreme earthquake events is evidently a difficult and delicate task to accomplish. Needless to say, the lack of direct and detailed information on the parameters at stake (including fault geometry, tectonic stress [3], slip patch [4], rupture path [5], geological configuration [6,7,8], multi-scale heterogeneities and anisotropies [9,10,11], non-linear site effects [12,13], the impact of surface topography [14], soil–structure interaction—SSI [15,16,17], etc.) increases the uncertainty of earthquake predictions with an unknown multivariate distribution of the seismic responses.

In recent years, the development of physics-based and engineering-oriented high performance computing (HPC) applications has widened the research horizons for seismology and earthquake engineering [18], allowing for the direct treatment of the previously mentioned uncertainty. High-fidelity numerical models can now be constructed [19,20,21,22,23], inputting the available information on the site/region of interest (including the extended source [24]) and exploring the uncertainty of the parameters and the panorama of possible scenarios [25,26]. In the past, one major drawback of the physics-based modeling approach resided in the difficulty of treating the high computational burdens that were required to routinely solve the wave propagation on such large domains (compared to the largest wavelength of interest [27]). At that point, it appeared necessary to build efficient multi-tool computational platforms for the assembly and calibration of the seismological model [28,29,30,31,32,33,34,35]. To this end, the main computational issues that need to be tackled are:

• To mesh the domain of interest [36], either by following the geological interfaces (honoring approach) or interpolating the mechanical properties on a structured mesh. The meshing scheme should follow the surface topography and the bathymetry (if present). The spatial refinement of the computational grid has to adapt to the minimum wavelength of interest [27], depending on the local value of wave speed. Ref. [24] points out that the surface geometry should be considered in numerical models for the accurate estimation of seismic forces. Topography modeling, however, creates new computational challenges;
• To represent the material rheology (i.e., elastic, viscoelastic, and non-linear hysteretic stress–strain relationships that are typical of geomaterials) [37,38];
• To describe the natural heterogeneity of the Earth’s crust and soil properties at different scales (i.e., regional geology, local basin-type structures, and the heterogeneity of granular materials) [10,39]. Particular attention should be paid to the remarkable amplification of surface waves due to basin-edge effects and shallow basin deposits [24];
• To introduce realistic rupture paths along the fault discontinuities through the distribution of moment tensor sources, i.e., the kinematic approach [40], or by simulating non-linear dynamic rupture along the fault plane [41]. However, the explicit inclusion of the extended source (planar, non-planar, or segmented) within the model is strongly advised [24];
• To describe the coupling between the wave-propagation kernel and the solvers for the structural dynamics, to simulate the SSI effects [16,17].

The focus of this paper is on two concurring aspects for improving the overall realism of earthquake predictions: (1) the adequate computational efforts required to increase the resolution of the wave-propagation kernel and (2) the inclusion of the heterogeneous fluctuations of the geological properties and source mechanism in order to accurately render the scattering effects. Regarding the first aspect, a new multi-tool numerical platform for generating 3D broadband source-to-site (BBS2S) earthquake scenarios is presented in Section 2. Figure 1 shows the scientific workflow of the BBS2S platform.

The core of the platform is represented by the wave-propagation solver, a software called SEM3D [42] tailored to efficiently solve the wave-propagation problem by means of the Spectral Element Method (SEM, see Section 2.1). SEM3D is complemented by a set of external libraries, either open-source or developed in-house, in order to construct the geological and seismological models. The computational improvements made in terms of BBS2S are the result of a synergistic effort of three French research teams at CentraleSupélec, Institut de Physique du Globe de Paris (IPGP), and the Commissariat à l’énergie atomique et aux énergies alternatives (CEA), mainly within framework of the French national project SINAPS@ (https://www.institut-seism.fr/projets/sinaps/, accessed on 24 January 2022).

Concerning the second main focus of this paper, our BBS2S is one of the few that integrates an efficient generator of random heterogeneous fluctuations of the mechanical properties of the Earth’s crust. As a matter of fact, smooth or discontinuous heterogeneity within geomaterials depends on how these effective grid parameters are evaluated. Kristek et al. [43] introduced an orthorhombic representation of a heterogeneous medium for the Finite Difference modeling. More recently,  Scalise et al. [11] performed large scale simulations of the Source Physics Experiment chemical explosions using 1-D and 3-D velocity models of the Yucca Flat basin in a 0–5 Hz frequency range. They added correlated stochastic velocity perturbations to the 3-D regional geological framework. In the same context, the heterogeneous fluctuations of the mechanical properties are introduced here by means of a HPC library, called randomField, developed within the BBS2S framework and linked to SEM3D (see Section 3). randomField represents a novel and competitive alternative to available tools to produce very large 3-D random fields without reducing the numerical efficiency of the wave-propagation kernel. Moreover, compared to [11], our BBS2S platform is hereafter tested in a 0–10 Hz frequency range and validated by drawing a digital twin of a European experimental site designed with the explicit purpose of testing physics-based numerical simulations. This site is represented by the Argostoli sedimentary basin, on Cephalonia Island (Western Greece) [44], which was the object of an extensive measurement campaign to monitor the seismic activity in the region [45,46], to constrain the geological basin-like model of the region [47], to assess the ground motion variability [48,49], and to test high-fidelity earthquake prediction by means of physics-based numerical simulation. All the geological and seismological information available at the site have been accurately modeled and verified in previous work [50,51], to validate the realism of the seismic response predicted at Argostoli site, following the validation matrix proposed by [52]—employing complex Intensity Measure (IM) metrics for site-specific ground-motion simulations—for both point-wise and extended source. The numerical model used in this paper is an improvement of the original one designed by [50], used to include the fluctuations of the mechanical properties within the sedimentary basin and assess the ground motion synthetic prediction at a high frequency (see Section 5). The outcome of the Argostoli earthquake simulation is discussed in Section 6, proving its realism compared to traditional ground-motion prediction equations.

The achievements presented in this paper represent a major step into the massive exploitation of high-fidelity and site-specific models to explore ground-shaking regional scenarios.

2. Wave-Propagation Numerical Solver

2.1. The Spectral Element Method in Seismology

The Spectral Element Method (SEM) is a powerful numerical method belonging to pseudo-spectral methods [53]. Based on the pioneering works in fluid dynamics performed by [54,55,56,57], the SEM was successfully applied to seismological applications to solve the 3-D wave equation in highly heterogeneous media [58,59], which reads:

$$\{\begin{array}{}\left(1\right)& {\nabla }_x·\sigma \left({\mathit{u}}_{}^{}\left({\mathit{x}}_{}^{};t_{}^{}\right)\right)+\mathit{b}\left({\mathit{x}}_{}^{};t\right)=\rho \left(\mathit{x}\right)\dot{\mathit{v}}\left(\mathit{x};t\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall \left(\mathit{x};t\right)\in \mathsf{\Omega }\times {\mathbb{I}}_t\hfill \left(2\right)& \rho \left(\mathit{x}\right)\mathit{v}\left(\mathit{x};t\right)=\rho \left(\mathit{x}\right)\dot{\mathit{u}}\left(\mathit{x};t\right)\hfill \end{array}$$

in an open domain $\mathsf{\Omega }\subseteq {\mathbb{R}}^3$ and a time interval ${\mathbb{I}}_t\subseteq {\mathbb{R}}^{+}$. ${\mathit{u}}_{}^{}\left({\mathit{x}}_{}^{};t_{}^{}\right)$ and ${\mathit{v}}_{}^{}\left({\mathit{x}}_{}^{};t_{}^{}\right)$ are the unknown displacement/velocity wave-field respectively; $\rho \left({\mathit{x}}_{}^{}\right)$ > 0 is the unit mass density, ${\mathit{b}}_{}^{}\left({\mathit{x}}_{}^{};t_{}^{}\right)$ is the body force density distribution (per unit mass) applied to the medium (e.g., the mass gravity), $\sigma \left({\mathit{x}}_{}^{};t_{}^{}\right)$ is Cauchy’s stress tensor (the small strain and small displacement assumptions are considered hereafter). $\dot{\square }$ represents the material time derivative $d\square /dt$. The relationship between and displacement field ${\mathit{u}}_{}^{}\left({\mathit{x}}_{}^{};t_{}^{}\right)$ is defined by heterogeneous linear elasticity: $\dot{\sigma }\left({\mathit{u}}_{}^{}\left({\mathit{x}}_{}^{};t_{}^{}\right)\right)=\mathbb{C}(\mathit{x}):{\mathit{v}}_{}^{}\left({\mathit{x}}_{}^{};t_{}^{}\right)$, with $\mathbb{C}(\mathit{x})$ being the elastic stiffness tensor.

The SEM belongs to the Finite Element Method (FEM) framework, with the unknown wave-field approximated by high-order Lagrangian polynomials, keeping low-order geometrical mapping on $N_e$ non-overlapping elements ${\overline{\mathsf{\Omega }}}_e$ (generally hexahedral elements), such that $\overline{\mathsf{\Omega }}={\cup }_{e=1,N_e}{\overline{\mathsf{\Omega }}}_e$ and the intersection between two distinct elements $e^{{}^{\prime }}$ and e ${\overline{\mathsf{\Omega }}}_{e^{{}^{\prime }}}\cap {\overline{\mathsf{\Omega }}}_e$ is an element’s vertex, edge, or face. In this context, the semi-discretized wave equation reads:

$$\{\begin{array}{}\left(3\right)& \mathbb{M}{\dot{\mathit{V}}}_G^h={\mathbb{F}}^{ext}-{\mathbb{F}}^{int}\left({\mathit{U}}_G^h;{\mathit{V}}_G^h\right)\hfill \left(4\right)& {\dot{\mathit{U}}}_G^h={\mathit{V}}_G^h\hfill \end{array}$$

with ${\mathit{U}}_G^h$, ${\mathit{V}}_G^h$ being the displacement and velocity vectors at the global degrees of freedom (DOF), respectively. $\mathbb{M}$ is the diagonal mass matrix (an interesting property, in terms of computational effort, inherited from the spectral discretization, and specifically from the orthogonal Lagrange polynomials employed). The vectors ${\mathbb{F}}^{ext}$ and ${\mathbb{F}}^{int}\left({\mathit{U}}_G^h;{\mathit{V}}_G^h\right)$ contain the external and internal forces, respectively [60]. The quadrature grid is constructed with ${(N+1)}^d$ integration points belonging to the Gauss–Lobatto–Legendre (GLL) set. The choice of GLL point positions, associated with Lagrangian shape functions, endows the SEM with so-called spectral precision [28,61]: the aliasing error exponentially vanishes with an increasing polynomial order. This aspect implies that with Lagrange polynomials of the fourth order, ≈5.7 GLL points per minimum S wave-length are sufficient to compute accurate wave-motion [27]. Another benefit resides in the fact that the vertices of the linear hexahedral element also belong to the GLL grid, leading to computing the exact stress field at the element interfaces.

The adopted time-marching scheme is a second-order accurate velocity–stress time-staggering Newmark scheme [62], exploiting the property of Lagrange polynomials being orthogonal on the GLL grid, thus providing a naturally diagonal mass matrix and a reduction in the non-zero stiffness coefficients to be computed.This explicit scheme requires a stability condition for numerical stability (the standard Courant–Friedrichs–Lewy condition, CFL), with the CFL number varying between 0.15 and 0.2 [63,64]. Accurate stability and convergence analyses are provided in [65,66,67]. Finally, efficient absorbing boundaries are provided by the Perfectly Matched Layers (PML) [68,69].

The SEM inherits from the FEM the ease of dealing with complex interfaces such as the basin edge. This represents a major advantage compared to the Finite Difference Method (FDM) for seismological applications [70] that requires more complicated numerical strategies to deal with those interfaces [7]. Moreover, the SEM formulation proposed by Festa and Vilotte [69] associates the mechanical properties to each GLL node, increasing the flexibility of spatial-interpolation schemes for highly heterogeneous media. The SEM can benefit from a Discontinuous Galerkin (DG) approach to effectively refine the mesh in the area of interest and mesh complex interfaces, such as sedimentary basins (see, for instance, [61,71]).

2.2. SEM3D

SEM3D is a high-performance numerical implementation of the SEM based on the original code called RegSEM (https://github.com/paulcup/RegSEM, accessed on 24 January 2022) [28,69], which was forked and further developed by CentraleSupélec, the IPGP, and the CEA since 2014. SEM3D represents a competitive alternative to other flagship codes employing the SEM, such as SPECFEM3D [53,59,70], EFISPEC3D [27], and SPEED [61].

SEM3D is a distributed-memory software (based on Message Passing Interface, MPI) written in Fortran 95 and C++, running on parallel MIMD (multiple instruction, multiple data) architectures. The computational domain $\mathsf{\Omega }$, duly discretized in $N_e$ hexahedral elements, is preliminary partitioned into ${\mathit{N}}_{MPI}$ subdomains (corresponding to the number of employed CPU cores) via a mesher that exploits the software library METIS 5.x (http://glaros.dtc.umn.edu/gkhome/views/metis, accessed on 24 January 2022), represented by a set of serial programs for partitioning graphs and finite-element meshes. The algorithms implemented are based on the multilevel recursive-bisection, multilevel k-way, and multi-constraint partitioning schemes (see Figure 2).

One original aspect is the optional SIMD (single instruction, multiple data) vectorization. The latter aspect is originally implemented via C++ macros in order to maximize the adequacy between data structures (member of a module) and architecture (vectorization possible or not). This is especially suitable for AVX512 vector machines. Other implementation schemes can be found in Sornet et al. [72]. Compared to RegSEM, I/O (Input/Output) operation is effectively performed using Hierarchical Data Format HDF5 (https://support.hdfgroup.org/HDF5, accessed on 24 January 2022) library. As a matter of fact, the HDF5 is associated with a programming interface that implements the abstract data and storage models. The library also implements a model of data transfer to efficiently move data from one stored representation to another stored representation. HDF5 implements the abstract model and maps the storage model to different storage mechanisms. HDF5 can effectively handle the data-aggregation problem that affects the I/O performances for simulation codes that compute for a large number of time steps, but each step only takes a small amount of time. As a matter of fact, data dumping at every step is known to cause I/O bottlenecking, which HDF5 libraries—properly tuned—can avoid. Therefore, this hierarchical data format effectively reduces the computational time required to handle partitioned mesh, to read material property files, and to write traces and snapshots. For instance,  Tang et al. [73] adopted HDF5 in order to efficiently tune parallel data compression and I/O for large-scale earthquake simulation, exploiting its data compression and chunking capabilities. In SEM3D, the I/O operations are shared by groups of 32 MPI processes and then separately written via HDF5.

Given the diagonal mass matrix provided by the quadrature rules, the SEM is naturally prone to parallelization schemes [74]: all local element matrices are computed independently, e.g., as in the elements-first algorithms [75]. In SEM3D, the time-marching Newmark’s scheme is solved on each DOF (i.e., at each GLL location on each element). The algorithm directly computes the product of the stiffness matrix by the displacement vector corresponding to ${\mathbb{F}}^{int}$ in Equation (Section 2.1). Therefore, the stiffness matrix is never assembled. This latter advantage, along with the diagonal mass matrix, leads to the asynchronous communication scheme employed by SEM3D to iterate at each time-step. This setup naturally accommodates for viscoelastic and non-linear constitutive relationships: in SEM3D, a generalized Zener model is implemented [76], as much as a simple non-linear rheological model [37]. The improvements made to the original RegSEM led to a considerable code speed-up, as Figure 3 shows. In [50], the author performed the verification and validation tests proposed by [6] for the Euroseistest Verification and Validation Project (E2VP) benchmark (the reference solution was uploaded on the SISMOWINE interactive seismological web interface used for numerical modeling benchmarking http://sismowine.org, accessed on 24 January 2022. HSP1a, Can2, and Can4 benchmarks were performed by [50]).

Poursartip et al. [24] compiled a table summarizing major large-scale time-domain simulations of seismic wave motion.

Table 1. Summary of major time-domain, large-scale simulations of seismic waves in the last 10 years (table after [24]). FDM: Finite Different Method, FEM: Finite Element Method, SEM: Spectral Element Method, DG: Discontinuous Galerkin. “-” indicates that the desired data are not available in the reference. This table is not exhaustive.

Ref.	Year	Method	Ressources	DOFs	Grid Size (m)	Size (km × km × km)	fmax (Hz)	Topography
[78]	2010	FDM	6 cores	-	25/125	-	2.5	no
	2010	SEM	32 cores	66,187,872	150	-	2.0	yes
	2010	SEM	63 cores	39,902,676	20–900	-	3.0	yes
	2010	DG	510 cores	-	200–5000	-	3.0	yes
[79]	2010	FEM	-	251,457,147	var	600 × 300 × 80	0.5	no
	2010	FDM	-	2.355 billions	200	500 × 250 × 50	0.5	no
	2010	FDM	-	5.419 billions	100	600 × 300 × 80	0.5	no
[80]	2010	SEM	192 GPU	131,000,256	-	chunk of earth	0.7	no
[81]	2010	FDM		1308 billions	40	810 × 405 × 85	2.0	no
[82]	2012	SEM	896 GPU	22 billions	24,000	Western Europe × 200	0.125	yes
[83]	2013	FEM	24,000 cores	15.9 billions	5.5 88	180 × 135 × 32	4.0	no
[84]	2014	DG	1,400,832 cores	96 billions	-	-	10	yes
[13]	2017	FEM	2,94,912 cores	10.7 billions	0.66	2 × 2 × 0.1	-	yes
[77]	2017	FDM	1,014,000 cores	23.4 trillions	8	320 × 312 × 40	18	yes
[85]	2021	FEM	1,179,648 cores	324 billions	0.125/64	256 × 205 × 100	-	yes
This paper	2022	SEM	4000 cores	13.5 billions	35 130	44 × 44 × 63	10	no

partially reports the mentioned table and it integrates the table compiled by [77] to show a panorama (non-exhaustive) of the most significant earthquake simulations performed in the last 10 years.

Table 1 highlights the fact that SEM3D represents a competitive alternative to modern broad-band synthetic earthquake simulators for engineering purposes (not considering pre-Exascale codes, such as SPECFEM3D [53,59,70] or ExaHyPE [86]): a rather large simulation (of ≈13.5 billions DOFs) could be run on a modern Petascale supercomputer, employing a rather small amount of computational resources.

A great effort has been made in order to grant SEM3D significant portability on several different supercomputers and grant it highly scalable properties. Figure 4 shows a comparison between the performances of SEM3D on two supercomputers: Fusion (this machine is owned by CentraleSupélec and Ecole Normale Supérieure de Paris Saclay, within the framework of Paris Saclay University (http://mesocentre.centralesupelec.fr/home/hardware/), accessed on 24 January 2022) (the Tier-2 [87] cluster is hosted at Mesocentre Moulon) and Occigen (this machine is hosted by French National Computing Center for Higher Education—CINES (https://www.cines.fr/en/supercomputing-2/), accessed on 24 January 2022) (a 3.5 pFlops Tier-1). The comparison is expressed in terms of CPU-time (mono-core equivalent simulation time) $t_{CPU}$ per second of real earthquake $t_{EQK}$ for different number of DOFs. The details of the computations (referenced as F1-F5 for Fusion and as O1-O2 for Occigen) are reported in Figure 4.

As shown in Figure 4, similar performances were obtained on both Fusion and Occigen, with an average $t_{CPU}/t_{EQK}$=1.772 10${}^3$ h/s. The most outstanding computation (i.e., O2, object of this article) was run on Occigen, where the maximum number of CPU cores per user is higher: the 3-D model consists of a box of 44 km × 44 km × 63 km, with the smallest grid-size of 130 m ×130 m × 35 m (at the surface) and a total of 1.3 × 10${}^{10}$ DOFs (polynomial order 9).

Concerning the weak scalability, Figure 5 shows SEM3D performances on a Tera 1000-1 supercomputer (https://www.top500.org/system/178790, accessed on 24 January 2022). In this graph, the number of elements treated per iteration is plotted for an increasing number of CPU cores.

The blue curve represents the theoretical linear scalability, i.e., considering the single chunk of 32 MPI processes (Tera 1000-1 is featured by bullx DLC 720 nodes, Xeon E5-2698v3 16C 2.3GHz, with Infiniband FDR). Green and red lines portray the best and average performances obtained over all the employed nodes. In each case, quasi-linear scalability curves are obtained, although the cluster working load may have an impact on the communication time.

2.3. Meshing Chunks of the Earth’s Crust: HexMesh

The code makes use of a library called HexMesh (https://github.com/jcamata/HexMesh.git, accessed on 24 January 2022), which implements an efficient, linear, 27-tree finite-element mesh-generation scheme based on the previous works of [88]) and is capable of generating large computational grids (i.e., ≈ 100 km) by extruding the provided Digital Elevation Model (DEM) and progressively coarsen it from top to bottom, so as to obtain a non-structured grid (see, for example, Figure 6). HexMesh handles coastlines and bathymetries.

3. Modeling Soil Heterogeneity by Means of Random Field

Modeling soil heterogeneity in large-scale domains implies a major computational effort: despite its characteristic dimension L (≈ 10${}^4$ m), the computational grid must resolve both the wavelength and the heterogeneity-correlation length ${\ell }_{C-\theta }$ (≈ 10${}^2$ m). Those large random fields can be effectively sampled over a coarse grid (with a step size relevant to the correlation length) and then interpolated onto the mesh of interest, provided by the GLL grid used in SEM3D.

To generate random samples, an open-source library randomField (https://github.com/cottereau/randomField.git, accessed on 24 January 2022) was developed by de Carvalho Paludo et al. [39] in the framework of BBS2S construction in order to be linked to SEM3D and feed it with large scale geological configurations with small random fluctuations. randomField is a part of SEM3D and exploits the spectral quadrature formula proposed by Shinozuka and Deodatis [89], Shinozuka and Deodatis [90], with a variation suited for isotropic media. This method generates the random field $G\left(\mathit{x}\right)$ as a sum of $N_{\varphi }=N_x{N}_y{N}_z$ cosines with random phases, namely:

$$G\left(\mathit{x}\right)=\sum _{n\le N_{\varphi }}\sqrt{2{\widehat{\mathcal{R}}}_{\theta }\left({\mathit{k}}_n^{}\right)|\Delta \mathit{k}|}{\zeta }_{n}cos\left({\mathit{k}}_n^{}·\mathit{x}+{\varphi }_n\right)$$

(5)

with ${\widehat{\mathcal{R}}}_{\theta }$ being the Fourier’s transform of the auto-covariance function ${\mathcal{R}}_{\theta }$. Following this approach, the wave-number domain spanned by k is discretized over a regular grid of size ${\mathit{N}}_{}^{}={\left[N_x,N_y,N_z\right]}^T$, indexed by ${\mathit{n}}_{}^{}={\left[n_x,n_y,n_z\right]}^T$. ${\widehat{\mathcal{R}}}_{\theta }\left({\mathit{k}}_n^{}\right)$ represents the power spectral density of the random field, whereas $|\Delta \mathit{k}|$ is the unit volume in the spectral domain. ${\varphi }_n$ are independent random phases with uniform probability density over $[0,2\pi ]$, whereas ${\zeta }_n$ are independent, random Gaussian amplitudes. By virtue of the central limit theorem, the generated random fields are Gaussian (for $N_{\varphi }$ asymptotically large). The Fast Fourier Transform is used (using the library FFTW (http://fftw.org/, accessed on 24 January 2022)) to bring the complexity of this generation method to $\mathcal{O}(N_{\varphi }log{N}_{\varphi })$) [89], but it requires a uniform grid. Since the GLL points are not uniformly distributed in space, the random field must then be interpolated to the $N_p$ GLL nodes.

randomField is written in Fortan95 and exploits HDF5 data structures for I/O operations. When dealing with large domains ($L\gg {\ell }_{C-\theta }$), it applies a generative strategy to overcome the scalability issue and efficiently reduce the generation time (see Figure 7). It generates realizations of $G\left(\mathit{x}\right)$ over the entire domain as superpositions of the $i\mathrm{th}$ smaller independent realizations $G_i\left(\mathit{x}\right)$ supported on overlapping subdomains ${\mathsf{\Omega }}_i$ of $\mathsf{\Omega }$ [39]:

$$G(\mathit{x})=\sum _{i\in I}\sqrt{{\psi }_i\left(\mathit{x}\right)}{G}_i(\mathit{x})$$

(6)

where the set of functions ${\psi }_i(\mathit{x})$ forms a partition of unity of $\mathsf{\Omega }$ (that is to say, ${\sum }_{i\in I}{\psi }_i(\mathit{x})=1$ for any $\mathit{x}\in \mathsf{\Omega }$), supported by the set of subdomains ${\mathsf{\Omega }}_i$.

The scheme decomposes the simulation domain into overlapping subdomains, each assigned to a single MPI process. Each processor generates its own statistically independent random field sample over its subdomain. After multiplying the local fields by the root of the partition of unity functions (see Equation (6)) and merging on the overlaps, the global fully correlated random field is obtained. Using this approach, the complexity becomes $\mathcal{O}(n_{p}log\left(n_p\right))$, where $n_p=N_{\varphi }/\mathrm{P}$ and $\mathrm{P}$ is the number of processors. Essentially, this means that the scheme is $\mathcal{O}\left(1\right)$ when we consider a constant number of GLL nodes per processor is considered. The overlapping Equation (6) involves an approximation that does not alter the average and variance of the resulting field $G(\mathit{x})$ [39]. The influence on the correlation structure depends on the overlap, relative to the correlation length, and has been studied by de Carvalho Paludo et al. [39].

Table 2. Weak scalability analysis for the random-field generation on Occigen cluster. Standard: wall-time required to generate a single random field over the entire domain. Localized: wall-time required to generate the random field with the localized approach mentioned in the text.

Cores	Nodes	Generation Time (s)
		Standard	Localized
16	$1\times {10}^6$	0.86	6.57
32	$2\times {10}^6$	1.55	13.87
256	$2\times {10}^7$	356.43	159.98
2048	$1\times {10}^8$	–	2037.48
4096	$3\times {10}^8$	–	3299.04

lists the results of the scalability test performed on randomField, with the standard and localized approaches, respectively (i.e., by generating a unique random field over the entire domain (standard approach) or by generating independent realizations over subdomains (localized approach). This scalability test was performed by running the generation code on an Occigen cluster (on 2.6 Ghz Intel Xeon Haswell dual-socket nodes with 12 cores each).

4. Modeling Extended Seismic Sources

Fault discontinuities are modeled by means of a kinematic approach: the slip across the fault surface $\Sigma$ (outward unit normal vector ${\mathit{n}}_{\Sigma }^{}$) is interpreted as a pre-existing displacement discontinuity $〚{\mathit{u}}_{}^{}〛\left({\mathit{x}}_{}^{};t_{}^{}\right)$ at a given surface. In practice, the kinematic assumption of extended faults consists in considering a distribution of double-coupled point-wise sources along the fault plane, with different time-varyng slip functions [91].

A very realistic kinematic source model, namely, Ruiz’s Integral Kinematic (RIK), was proposed by Ruiz et al. [92]. Specifically, the numerical simulations presented hereafter are based on the slip patches provided by the RIK implementation presented by Gallovič [40], called RIKsrf (https://github.com/fgallovic/RIKsrf, accessed on 24 January 2022). The rupture process on the fault plane (defined as a rectangular surface of along-strike length L and width W) is described as an assembly of numerous smaller events generated by circular cracks (multi-crack model). Crack’s radius distribution follows a truncated fractal distribution bounded by the lower and upper limits of crack size $\{R_{min};R_{max}\}$, with $R_{min}$ generally set to half the discretization step $\Delta x/2$ and $R_{max}=\alpha min\left(W,L\right)$, with $\Delta x$ being the minimum distance between the centers of two cracks and $\alpha <0.5$.

The distribution of the crack’s centers $C\left({\mathit{x}}_{\Sigma }^{}\right)$ (with C being a random field representing the position of a crack’s center) is randomly defined by sampling its probability density function $p_C$. The latter can be chosen arbitrarily (Gaussian or uniform, for instance): in this study, the probability distribution is conditioned by the slip distribution patch obtained by source-inversion techniques. Each sub-event is characterized by a maximum slip function $\Delta u_c\left(r\right)$ defined along the Eshelby’s crack surface [93].

The nucleation point within each crack is randomly chosen to grant better representation of the high-frequency directivity [92]. in this way, a sub-event does not activate when the rupture front perturbs the source, but it starts radiating only when the rupture front encounters the nucleation point within the nucleation region. The latter is a scale-dependent region, defined by the distance $d\left(R\right)$:

$$d\left(R\right)\equiv \{\begin{array}{}\left(7\right)& 2hR,\hspace{1em}R<R\hfill \left(8\right)& 2h{R}_c,\hspace{1em}R\ge R\hfill \end{array}$$

with $R_c$ being a predefined critical radius and h being the parameter controlling the extension of this region. If $h=0$, the nucleation point is located at the intersection between the rupture front and the crack’s edge. For $h>0$, the nucleation point is randomly picked within the defined nucleation region.

When the nucleation starts, the rupture front crosses the crack, which radially radiates within its surface, until the $\Delta u_c^{max}\left(r\right)$ is reached. The time evolution of the slip at each crack follows a boxcar source-time function with a scale-dependent rise time, with $\tau \left(R\right)$ defined as:

$$\tau \left(R\right)\equiv \{\begin{array}{}\left(9\right)& \frac{2R}{V_r},\hspace{1em}R<R_p\hfill \left(10\right)& \frac{2R_p}{V_r},\hspace{1em}R\ge R_p\hfill \end{array}$$

where $V_r$ is the constant rupture velocity and $R_p$ is the sub-event size that has the maximum rise-time ${\tau }_{max}=R_p/V_r$. Further details on this multi-crack model are provided in Ruiz et al. [92].

In order to integrate the RIK extended fault model into the BBS2S platform, the point-wise seismic moment rate time-histories generated by RIKsrf are stored in HDF5 files, which are then read by SEM3D at run time in order to assemble the external force vector ${\mathbb{F}}^{ext}$. Each of the point-wise sources produced by RIKsrf generates a couple of forces that are integrated over the mesh element ${\mathsf{\Omega }}_e$, which belong to [94]. The advantage of the described strategy to include an extended source into the numerical model (as suggested by Poursartip et al. [24] as best practice) is twofold: (1) an arbitrary slip patch provided, for instance, by source wave-form inversion, can serve as a background model to which we can add high-frequency content via heterogeneous crack distributions; (2) the mesh is not required to comply with the distribution of point sources constituting the kinematic extended fault, granting a certain modularity to the whole platform (although numerical accuracy relies on the level of discretization of both the source plane and the portion of the mesh accommodating it).

5. Source-to-Site Earthquake Simulation at the Argostoli Site

In order to show an application of the features described above, the case study of the Argostoli site (Greece) was performed to prove the predictive capability of high-fidelity earthquake simulations based on all the geological and seismological information available at the site of interest. The Argostoli-Koutavos sedimentary basin is located on the Kefalonia Island in Greece, at the bottom of a lagoon filled with quaternary and Pliocene detritic deposits [47]. This area represents the seismically most active part of Europe. High seismic activity has been recorded in its proximity, caused by the Kefalonia Transform Fault [45], which is a connection boundary close to the island of two subduction troughs: the Hellenic Arc in the south and the Adriatic Fault Zone [95]. However, fault modeling in this region is a major challenge because of the structural complexity and poor coverage by seismological, GPS and InSAR data [96]. The site was originally chosen for its basin-like geological configuration, which makes it attractive to study 3-D site effects and surface waves. An extremely interesting amount of data regarding the area is fully documented in the NERA (Network of European Research Infrastructures for Earthquake Risk Assessment and Mitigation) project deliverable D11.6 Comparison between data and numerical models (https://www.orfeus-eu.org/other/projects/nera/, accessed on 24 January 2022). Based on this state of the art,  Touhami [50] constructed the numerical model used hereafter. Two possible realizations of the earthquake scenario are hereafter compared, run on two different supercomputers (Fusion and Occigen).

5.1. Geological Characteristics

Two velocity models were combined in this work, namely: (1) the regional crustal velocity model implied by Haslinger et al. [97] via a wave-form inversion procedure, refined at the surface by (2) the layered geology profile proposed by Hollender et al. [98] (see Figure 8 and

Table 3. Mechanical properties of the geological profile proposed by [98] for the first 5 km down. $V_P$ is the pressure-wave velocity, $V_S$ is the shear-wave velocity, $Q_P$ and $Q_S$ are the quality factors for P- and S-waves, respectively.

${\mathit{z}}_{\mathbf{top}}$ (m)	${\mathit{V}}_{\mathit{P}}$ (m/s)	${\mathit{V}}_{\mathit{S}}$ (m/s)	${\mathit{Q}}_{\mathit{P}}$ (1)	${\mathit{Q}}_{\mathit{S}}$ (1)
0	2000	600	300	150
300	2400	1000	300	150
400	4600	2700	300	150
1000	6000	3200	300	150
2000	6200	3200	300	150

). The basin-like structure was obtained via in situ geophysical measurements by Cushing et al. [47], Cushing et al. [99]. The mechanical properties within the basin were taken from Hollender et al. [98].

Figure 9a is a top view of the area surrounding the Argostoli basin, depicting the contour of shear-wave velocity $V_S$, as per Cushing et al. [99]. For the sake of clarity, Figure 9b reports a 3-D view of two cross-sections cutting the Argostoli basin across its longitudinal and transversal extensions. The 3-D basin-like structure can be effectively noticed. Table 3 reports the geological characteristics provided by Cushing et al. [47] for the soft soil sediments within the Argostoli basin depicted in Figure 9.

Figure 10 shows a clip of the mesh generated for the case (44 km × 44 km × 63 km with a minimum element size of 130 m × 130 m × 35 m), featured by $\approx 4.5·{10}^6$ elements and ≈9.8–13.5$·{10}^9$ (10 × 10 × 10 GLL × 3 per element). No surface topography was considered at this stage (see Figure 10a) in order to concentrate on the influence of the small-scale heterogeneities on the seismic site response. Further considerations of the effects of topography on the Argostoli site response simulation can be found in [51].

The mesh is designed to accommodate the whole geological profile provided by Haslinger et al. [97] so as to include all the information available and exploit it for different earthquake scenarios. In any case, a progressive refinement is applied towards the free surface, and the extra layers at depth do not affect the overall computational burden. Moreover, those deep layers help in attenuating the wave-field radiated towards the PML layers enveloping the model.

Touhami et al. [51] performed the verification of the numerical model to compare the recorded and synthetic rock-to-site amplification functions, as well as the seismic response to an aftershock belonging to the sequence described in [100].  Touhami [50] performed a stress test of the seismic response of the basin for small point-wise aftershocks and obtained a satisfactory fit of the two first resonant frequencies measured in situ by Perron et al. [101]. This works represent the necessary verification task recommended by Bradley [52] and following the methodological approach presented in [102,103] for the Mygdonian basin E2VP test case (Greece). Moreover,  Touhami [50] performed three extended fault seismic scenarios of the 2014 $M_W$6.0 Cephalonia earthquake, following the refined fault models proposed by Saltogianni et al. [96]. One of those scenarios, corresponding to a strike-slip fault model, was selected for this study, enlarging the frequency content of the numerical analysis up to 10 Hz.

One of the novelties of this paper consists of adding the fluctuation of the soil properties within the sedimentary basin by employing randomField coupled with SEM3D so as to obtain a spatially variable non-stationary random field for the unit mass density $\rho \left({\mathit{x}}_{}^{}\right)$ for the light-gray area shown in Figure 10a and therefore a non-stationary elastic shear modulus $\mu \left({\mathit{x}}_{}^{}\right)=\rho \left({\mathit{x}}_{}^{}\right)c_S^2\left({\mathit{x}}_{}^{}\right)$. $c_S\left({\mathit{x}}_{}^{}\right)$ is the average shear wave velocity described above. $\rho \left({\mathit{x}}_{}^{}\right)$ is generated via randomField and it is characterized by an anisotropic Von Karman Correlation structure (${\ell }_{C-x}={\ell }_{C-y}$ = 90 m, ${\ell }_{C-z}$ = 40 m) over a domain of ≈7.3 km × 6.5 km × 0.1 km. Its standard coefficient of variation CoV = ${\sigma }_{\rho }/{\mu }_{\rho }$ is equal to 20%, with a first-order marginal that is assumed to follow a log-normal pdf distribution.

5.2. Seismic Source

A target $M_W$6.0 earthquake was simulated using the implemented RIK kinematic source model described in Section 4. The fault segment has a 15 km × 10 km extension and its orientation is defined by the strike, dip, rake triplet ${\varphi }_S$ = 38${}^{\circ }$, $\delta$ = 63${}^{\circ }$, $\lambda$ = 172${}^{\circ }$. The hypocenter (represented by a red star in Figure 11) is placed at a depth of 10 km. This fault configuration corresponds to the Cephalonia segment, one of the branches identified along the Cephalonia Transform Fault [45]. Moreover, the adopted slip patch corresponds to a plausible source model for a target earthquake of $M_W$6.0 (compatible with [100]). Other plausible fault models were provided by [96]. The slip patch is discretized on a grid of 150 × 100 points. Figure 11 shows the slip distribution obtained with the RIK integral model.

6. Results and Discussion

In this section, the seismic response of the Argostoli basin is studied by comparing two simulations with two different levels of accuracy (i.e., the same mesh grid and material parameters but with different numbers of GLL integration points), namely: (1) the simulation labeled as F5 and (2) the simulation labeled as O2 in Figure 4b. Those simulations are respectively carried out on Fusion and Occigen supercomputers (see Section 2.2) and use 7 × 7 × 7 and 10 × 10 × 10 GLL points, respectively (p-refinement). In terms of accuracy, the two simulations have a theoretical $f_{max}$ = 7 Hz and $f_{max}$ = 10 Hz within the softest sedimentary layer (the sedimentary basin, with $V_S$ = 650 m/s).

For the sake of comparison, three monitoring points were selected (see Figure 9): P3 and P6, located within the basin (P3 at -83 m G.L. and P6 at 0 m G.L.) and P262, located at −200 m G.L. outside the basin.

O2 simulation is selected to test the effect of crustal heterogeneity, simulated via randomField (see Section 3). A heterogeneous density $\rho \left(x\right)$ field was generated, featuring a von Karman auto-correlation function with correlations lengths $\left({\ell }_{C-x},{\ell }_{C-y},{\ell }_{C-z}\right)$ = (90.0 m, 90.0 m, 40.0 m) and a Coefficient of Variation CoV=0.2. Three monitoring points were selected: P3 (located at −83 m G.L., within the basin), P6 (located at the surface), and P262 (located at −200 m G.L., below the heterogeneous layer). Compared to [11], the ratio between the domain domain sizes $L_x,L_y,L_z$ and the correlation lengths was $\left({\ell }_{C-x},{\ell }_{C-y},{\ell }_{C-z}\right)$ .

From a regional point of view, both numerical simulations seem to agree with the Ground Motion Prediction Equation (GMPE) proposed by [104], as Figure 12 shows.

Both the $M_W$6.0 simulations comply with the GMPE proposed by Akkar et al. [104]: most of the stations considered in Figure 12 showed a Pseudo Spectral Acceleration ($PSA$) value at T = 0.1 s close to the median value m and within $m\pm \sigma$ (standard deviation) uncertainty margins. All the considered stations are located in very near field, i.e., around a hypocentral distance $R_{hyb}\approx$10–20 km. A few stations, such as P3 (see Figure 9), underestimate the GMPE prediction, because they are located underground, below the heterogeneous basin interface. This explains the fact that the $PSA$ predicted by F5 and O2 simulations at those locations are the same. Such control points were introduced, as suggested by Berge-Thierry et al. [105], in order to check that the regional ground motion model was not affected by the change in mechanical properties within the basin and that the incident wave-field remained unchanged from one simulation to another. The large dispersion in Figure 12 may be essentially linked to two major aspects: (1) the directivity effects engendered by this particular realization of the source scenario and (2) the arbitrary slip patch adopted for the blind prediction. The near-field earthquake ground motion prediction is extremely sensitive to local details, both in terms of geology and source. Figure 13 shows the Peak Ground Acceleration (PGA) contour plot across the Argostoli sedimentary basin, obtained by post-processing the outcome of the F5 simulation. The extended fault is located northeastern-ward. The major amplification is shown within the basin edges, within soft sediments. The average horizontal PGA is approximately five times larger than the vertical one.

Figure 14 and Figure 15 picture the comparison between homogeneous and heterogeneous crustal properties (see Section 5.1) in terms of time-histories extracted from the SEM3D outcome for model O2. Those time-histories are computed by SEM3D at the receiver’s locations, using the SEM shape function to interpolate the solution at each DOF.

As expected, the comparison between the results of the two models (same source and large scale geology) provides fairly similar synthetics at P262 (Figure 14) located below the 200 m heterogeneous soil deposit. Major discrepancies arise at the surface, though (at P6, Figure 15): the homogeneous model predicts larger amplitudes by an approximate factor of 2. This was also reported by, e.g.,  Jussila et al. [106], and it has an important consequence. Namely, the homogeneous model generally gives an upper bound and conservative prediction. From the practical standpoint, this aspect may guide the modeling of the small-scale heterogeneity. The corresponding Fourier’s spectra show the lower amplitude of the propagated ground shaking when heterogeneity is considered, and even a predominant frequency shift for EW direction (due to the effect of multiple scatterings).

A major difference between the two simulations is noticed in the low-frequency part of the propagated synthetic ground motion at the surface. Figure 16—where the signals are band-pass filtered between 0.1 and 5 Hz)—shows out this discrepancy quite clearly for monitor point P3.

The outcome of those two simulations proves the influence of the chosen polynomial order on the low-frequency part of the ground motion for complex 3-D test cases and complex material properties. The high amplification on the NS component is presumably due to directivity effects given the assumed focal mechanism. A higher-order polynomial interpolation better renders the complex source (a finer discretization renders the larger wavelengths generated by the source mechanism more accurately) and path effects (geological interfaces are better interpolated onto a finer grid, modifying the local medium configuration and, indirectly, the propagated wave field and scattering at high-frequency), as it is noticeable in the larger duration of the S-wave component of ground motion. This spatial-aliasing effect reduces the surface maximum amplification at P3 by an approximate factor of 1.5 on the horizontal components and 2.5 on the vertical one. Nevertheless, better accuracy and dispersion analyses are required for this complex 3-D highly heterogeneous wave-propagation problem. According to the theoretical results of Oliveira [107], the a priori error is estimated to be around 10${}^{-9}$, given the polynomial order 9 (Occigen run) and a time step $\Delta t$ = 10${}^{-4}$. This error can preferably be reduced by h-refinement [107].

Figure 17 displays both the evolution of the particle velocity contour plots nearby the basin and the slip rate in the extended source for four time steps. The contours are extracted by post-processing the results of the O2 simulation. The heterogeneous nature of the slip patch evolving in time (the blurred rupture front is due to the simultaneous rupture evolution of different sub-asperities from random nucleation points) can be appreciated, as well as the ground-shaking amplification when the incident wave-field enters the basin-like structure.

7. Conclusions

In this paper, we present an HPC platform for source-to-site wave propagation and broad-band regional-scale scenarios. The performances of the code are highlighted and tested on the synthetic simulation of the seismic response of the experimental site of Argostoli (Greece). The earthquake simulator is quasi-physics-based, since the source mechanism is not reproduced by dynamic rupture, but by following the classical kinematic description [50]. The present contribution aims at showing the high-frequency predictive power of the crafted numerical model for future plausible earthquake scenarios. The results of two simulations are compared, differing in terms of numerical accuracy (7 Hz and 10 Hz). As an improvement of the study performed by Touhami [50], Touhami et al. [51], the introduction of heterogeneous slip patches and heterogeneous fluctuations of the mean geological properties entails the generation of realistic coda-waves due to multiple scatterings. The low-frequency portion of the synthetic earthquake response undergoes slight modification when increasing the polynomial order of the underlying shape functions, mostly due to different spatial interpolations of the source points within each element. This aspect suggests great care when transitioning to the higher-frequency range and when introducing heterogeneous material properties. One realization of the seismic scenario (both in terms of source characteristics and random fluctuations of the material properties) cannot be considered exhaustive: however, this paper shows the feasibility of this strategy exploiting physics-based earthquake simulators to span the scenario and uncertainty space in an almost blind predictive framework. SEM3D can be effectively used to generate hybrid broad-band synthetic ground motion and to estimate the intra-event and inter-event residuals of ground-motion prediction equations, in areas of low–moderate seismicity, such as France, where a poor database can prevent the definition of representative ground-motion prediction equations.