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

Sara Touhami, Filippo Gatti, Fernando Lopez-Caballero, Régis Cottereau,
Lúcio de Abreu Corrêa, Ludovic Aubry, Didier Clouteau

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sara Touhami, Filippo Gatti, Fernando Lopez-Caballero, Régis Cottereau, Lúcio de Abreu Corrêa, et al.. SEM3D: A 3D High-Fidelity Numerical Earthquake Simulator for Broadband (0-10 Hz) Seismic Response Prediction at a Regional Scale. Geosciences, 2022, 12 (112), 10.3390/geosciences12030112. hal-03593533

HAL Id: hal-03593533

https://hal.science/hal-03593533

Submitted on 2 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

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

Sara Touhami <fnref n="1" />,† ID, Filippo Gatti <fnref n="1" />,*,† ID, Fernando Lopez-Caballero <fnref n="1" /> ID, Régis Cottereau <fnref n="2" /> ID,
Lúcio de Abreu Corrêa <fnref n="3" /> ID, Ludovic Aubry <fnref n="4" /> and Didier Clouteau <fnref n="1" /> ID

<figure><img src="image_3.png" /></figure>

### check for updates

**Citation:** Touhami, S.; Gatti, F.; Lopez-Caballero, F.; Cottereau, R.; de Abreu Corrêa, L.; Aubry, L.; Clouteau, D. SEM3D: A 3D High-Fidelity Numerical Earthquake Simulator for Broadband (0–10 Hz) Seismic Response Prediction at a Regional Scale. *Geosciences* **2022**, *12*, 112. https://doi.org/10.3390/geosciences12030112

Academic Editors: Maria D'Amico, Chiara Smerzini and Jesus Martinez-Frias

Received: 25 January 2022

Accepted: 23 February 2022

Published: 2 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**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¹⁰ 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.

**Keywords:** computational earthquake engineering; strong ground motion prediction; site effects; HPC numerical simulation

## 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, ge materials 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].

<footnote n="1">¹ Laboratoire de Mécanique Paris-Saclay, Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, 8-10 rue Joliot Curie, 91190 Gif sur Yvette, France; sara.touhami@centralesupelec.fr (S.T.); fernando.lopez-caballero@centralesupelec.fr (F.L.-C.); didier.clouteau@centralesupelec.fr (D.C.)</footnote>

<footnote n="2">² UMR 7031 Laboratoire de Mécanique et Acoustique, Aix-Marseille University, CNRS, Centrale Marseille, impasse Nikola Tesla, 13453 Marseille, France; cottereau@lma.cnrs-mrs.fr</footnote>

<footnote n="3">³ Post-Graduation Program in Mining, Metallurgy and Materials Engineering (PPGEM), Federal University of Rio Grande do Sul (UFRGS), Porto Alegre 90040-060, Brazil; lcorrea@demet.ufrgs.br</footnote>

<footnote n="4">⁴ CEA, DAM, DIF, F-91297 Arpajon, France; ludovic.aubry@cea.fr</footnote>

<footnote>* Correspondence: filippo.gatti@centralesupelec.fr</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
Due to this, the assessment of the seismic responses of critical structures and infras-
tructures 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], geolog-
ical configuration [6–8], multi-scale heterogeneities and anisotropies [9–11], non-linear site
effects [12,13], the impact of surface topography [14], soil-structure interaction—SSI [15–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 seismol-
ogy and earthquake engineering [18], allowing for the direct treatment of the previously
mentioned uncertainty. High-fidelity numerical models can now be constructed [19–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 pos-
sible 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–35]. To this end, the main computational issues that need to be tackled are:

1. 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;

2. To represent the material rheology (i.e., elastic, viscoelastic, and non-linear hysteretic stress–strain relationships that are typical of geomaterials) [37,38];

3. 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];

4. 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];

5. 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.
<figure><img src="image_5.png" /></figure>

**Figure 1.** BBS2S flowchart. The wave-propagation kernel is represented by the HPC code SEM3D. The latter is fed by a hexahedral mesh, created with the software HexMesh, which extrudes the Digital Elevation Model (DEM), creating the connectivity matrix. The mesh is decomposed in subdomains by employing the open-source METIS library. The available geological information (interfaces, layers) is introduced and small heterogeneous fluctuations are added by using the library randomField, developed within the BBS2S project. The seismological information, represented by the available fault models (focal mechanism and slip inversion) are used to generate heterogeneous and broad-band slip distributions on a fault plane using the software RIKsrf. Topographical, geological, and seismological data are introduced into SEM3D via HDF5 library. The outcome is represented by regional snapshots and point-wise time-histories.

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–57], the SEM was successfully applied to seismological applications to solve the 3-D wave equation in highly heterogeneous media [58,59], which reads:

$$ \left\{ \begin{array}{l} \nabla_x \sigma(\mathbf{u}(x;t)) + \mathbf{b}(x;t) = \rho(x) \dot{\mathbf{v}}(x;t) \\ \rho(x) \mathbf{v}(x;t) = \rho(x) \dot{\mathbf{u}}(x;t) \end{array} \right. \quad \forall (x;t) \in \Omega \times \mathbb{I}_t \quad (1) $$

$$ \left\{ \begin{array}{l} \dot{\mathbf{u}}_G^h = \mathbf{v}_G^h \\ \dot{\mathbf{v}}_G^h = \mathbf{u}_G^h \end{array} \right. \quad (2) $$

in an open domain $\Omega \subseteq \mathbb{R}^3$ and a time interval $\mathbb{I}_t \subseteq \mathbb{R}^+$. $\mathbf{u}(x;t)$ and $\mathbf{v}(x;t)$ are the unknown displacement/velocity wave-field respectively; $\rho(x) > 0$ is the unit mass density, $\mathbf{b}(x;t)$ is the body force density distribution (per unit mass) applied to the medium (e.g., the mass gravity), $\sigma(x;t)$ is Cauchy’s stress tensor (the small strain and small displacement assumptions are considered hereafter). $\square$ represents the material time derivative $d\square/dt$. The relationship between $\sigma$ and displacement field $\mathbf{u}(x;t)$ is defined by heterogeneous linear elasticity: $\dot{\sigma}(\mathbf{u}(x;t)) = \mathbb{C}(x) : \mathbf{v}(x;t)$, with $\mathbb{C}(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 $\bar{\Omega}_e$ (generally hexahedral elements), such that $\bar{\Omega} = \cup_{e=1,N_e} \bar{\Omega}_e$ and the intersection between two distinct elements $e'$ and $e$ $\bar{\Omega}_{e'} \cap \bar{\Omega}_e$ is an element’s vertex, edge, or face. In this context, the semi-discretized wave equation reads:

$$ \left\{ \begin{array}{l} \mathbb{M} \dot{\mathbf{V}}_G^h = \mathbb{F}^{\text{ext}} - \mathbb{F}^{\text{int}} (\mathbf{u}_G^h; \mathbf{V}_G^h) \\ \dot{\mathbf{u}}_G^h = \mathbf{V}_G^h \end{array} \right. \quad (3) $$

$$ \dot{\mathbf{u}}_G^h = \mathbf{V}_G^h \quad (4) $$

with $\mathbf{u}_G^h, \mathbf{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}^{\text{ext}}$ and $\mathbb{F}^{\text{int}}(\mathbf{u}_G^h; \mathbf{V}_G^h)$ 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-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 $\Omega$, duly discretized in $N_e$ hexahedral elements, is preliminary partitioned into $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).

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Example of the mesh partitioning performed by means of METIS software.</figcaption></figure>
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 indepen-
dently, 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 correspond-
ing to *F*int in Equation (4). 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]).

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Improvement of SEM3D compared to initial RegSEM version. The values of CPU time refer to the equivalent-mono-core time required to solve a wave-propagation problem in a 600 m × 600 m × 600 m cube, with PML on all sides and 5×5×GLL per element of 50 m × 50 m × 50 m.</figcaption></figure>

Poursartip et al. [24] compiled a table summarizing major large-scale time-domain
simulations of seismic wave motion. Table 1 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><caption>**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.</caption><thead><tr><th>Ref.</th><th>Year</th><th>Method</th><th>Ressources</th><th>DOFs</th><th>Grid Size (m)</th><th>Size (km×km×km)</th><th>f<sub>max</sub> (Hz)</th><th>Topography</th></tr></thead><tbody><tr><td rowspan="4">[78]</td><td>2010</td><td>FDM</td><td>6 cores</td><td>-</td><td>25/125</td><td>-</td><td>2.5</td><td>no</td></tr><tr><td>2010</td><td>SEM</td><td>32 cores</td><td>66,187,872</td><td>150</td><td>-</td><td>2.0</td><td>yes</td></tr><tr><td>2010</td><td>SEM</td><td>63 cores</td><td>39,902,676</td><td>20-900</td><td>-</td><td>3.0</td><td>yes</td></tr><tr><td>2010</td><td>DG</td><td>510 cores</td><td>-</td><td>200-5000</td><td>-</td><td>3.0</td><td>yes</td></tr><tr><td rowspan="3">[79]</td><td>2010</td><td>FEM</td><td>-</td><td>251,457,147</td><td>var</td><td>600×300×80</td><td>0.5</td><td>no</td></tr><tr><td>2010</td><td>FDM</td><td>-</td><td>2.355 billions</td><td>200</td><td>500×250×50</td><td>0.5</td><td>no</td></tr><tr><td>2010</td><td>FDM</td><td>-</td><td>5.419 billions</td><td>100</td><td>600×300×80</td><td>0.5</td><td>no</td></tr><tr><td>[80]</td><td>2010</td><td>SEM</td><td>192 GPU</td><td>131,000,256</td><td>-</td><td>chunk of earth</td><td>0.7</td><td>no</td></tr><tr><td>[81]</td><td>2010</td><td>FDM</td><td>-</td><td>1308 billions</td><td>40</td><td>810×405×85</td><td>2.0</td><td>no</td></tr><tr><td>[82]</td><td>2012</td><td>SEM</td><td>896 GPU</td><td>22 billions</td><td>24,000</td><td>Western Europe × 200</td><td>0.125</td><td>yes</td></tr><tr><td>[83]</td><td>2013</td><td>FEM</td><td>24,000 cores</td><td>15.9 billions</td><td>5.5 88</td><td>180×135×32</td><td>4.0</td><td>no</td></tr><tr><td>[84]</td><td>2014</td><td>DG</td><td>1,400,832 cores</td><td>96 billions</td><td>-</td><td>-</td><td>10</td><td>yes</td></tr><tr><td>[13]</td><td>2017</td><td>FEM</td><td>2,94,912 cores</td><td>10.7 billions</td><td>0.66</td><td>2 × 2 × 0.1</td><td>-</td><td>yes</td></tr><tr><td>[77]</td><td>2017</td><td>FDM</td><td>1,014,000 cores</td><td>23.4 trillions</td><td>8</td><td>320×312×40</td><td>18</td><td>yes</td></tr><tr><td>[85]</td><td>2021</td><td>FEM</td><td>1,179,648 cores</td><td>324 billions</td><td>0.125/64</td><td>256 × 205 × 100</td><td>-</td><td>yes</td></tr><tr><td>This paper</td><td>2022</td><td>SEM</td><td>4000 cores</td><td>13.5 billions</td><td>35 130</td><td>44 × 44 × 63</td><td>10</td><td>no</td></tr></tbody></table>

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 (monocore 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 \times 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.
<figure><img src="image_8.png" /></figure>

<table><thead><tr><th></th><th>N<sub>EL</sub></th><th>N<sub>GLL</sub></th><th>N<sub>MPI</sub></th></tr></thead><tbody><tr><td>F1</td><td>2×10<sup>6</sup></td><td>5</td><td>648</td></tr><tr><td>F2</td><td>-</td><td>7</td><td>648</td></tr><tr><td>F3</td><td>-</td><td>10</td><td>648</td></tr><tr><td>F4</td><td>18×10<sup>6</sup></td><td>5</td><td>648</td></tr><tr><td>F5</td><td>4.5×10<sup>6</sup></td><td>7</td><td>648</td></tr><tr><td>O1</td><td>-</td><td>5</td><td>4800</td></tr><tr><td>O2</td><td>-</td><td>10</td><td>4000</td></tr></tbody></table>

**Figure 4.** *N*<sub>E</sub>: Hexahedral element number; *N*<sub>GLL</sub>: number of GLL points per element); *N*<sub>MPI</sub>: number of CPU cores. F1-F4 and O1 refer to the simulation of the seismic response of the Kashiwazaki-Kariwa Nuclear Power Plant (Japan) during the 2007 Niigata *M*<sub>W</sub>6.6 earthquake. F5 and O2 refer instead to the seismic response of the Argostoli site, object of this paper. (a) CPU-time per real earthquake simulation *t*<sub>CPU</sub>/*t*<sub>EQK</sub>, for different DOF numbers. (b) Synoptic summary of SEM3D performances.

<figure><img src="image_9.png" /></figure>

**Figure 5.** Weak scalability curves obtained for SEM3D on Tera 1000-1 (Europe's first petaflop computer, hosted at CEA, France), expressed as time rate of mesh elements *N*<sub>el</sub> processed per chunk of 32 MPI core each. In blue, the linear scalability curve represents the ideal *optimum* result (perfect linear scalability). Green and red lines portray the *best* and *average* performances obtained instead.

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.

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Example of Digital Elevation Model including bathymetry line of the Niigata region (Japan) and hexahedral mesh obtained by `Hexmesh` code. (a) Digital Elevation Model (DEM) of the Niigata area (Japan); (b) Mesh of the Niigata area, with detail on the transition elements (reprinted from [35]).</figcaption></figure>

# 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,90], with a variation suited for isotropic media. This method
generates the random field $G(x)$ as a sum of $N_\phi = N_x N_y N_z$ cosines with random phases, namely:

$$G(x) = \sum_{n \le N_\phi} \sqrt{2\hat{R}_\theta(k_n) |\Delta k| \zeta_n \cos(k_n \cdot x + \phi_n)} \quad (5)$$

with $\hat{R}_\theta$ being the Fourier's transform of the auto-covariance function $R_\theta$. Following this approach, the wave-number domain spanned by $k$ is discretized over a regular grid of size $N = [N_x, N_y, N_z]^T$, indexed by $n = [n_x, n_y, n_z]^T$. $\hat{R}_\theta(k_n)$ represents the power spectral density of the random field, whereas $|\Delta k|$ is the unit volume in the spectral domain. $\phi_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_\phi$ 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 $O(N_\phi \log N_\phi)$ [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 Fortran95 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(x)$ over the entire domain as superpositions of the $i^{th}$ smaller independent realizations $G_i(x)$ supported on overlapping subdomains $\Omega_i$ of $\Omega$ [39]:

$$G(x) = \sum_{i \in I} \sqrt{\psi_i(x)} G_i(x) \quad (6)$$

where the set of functions $\psi_i(x)$ forms a partition of unity of $\Omega$ (that is to say, $\sum_{i \in I} \psi_i(x) = 1$ for any $x \in \Omega$), supported by the set of subdomains $\Omega_i$.

<figure><img src="image_11.png" /><figcaption>Figure 7. Generation of four independent realizations $G_i(x)$ on local supports $\Omega_i$. The green band represents the overlapping region, i.e., the support of the weighting functions $\phi_i$.</figcaption></figure>

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 $O(n_p \log(n_p))$, where $n_p = N_\phi/P$ and P is the number of processors. Essentially, this means that the scheme is $O(1)$ 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(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 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).

<table><caption>**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.</caption><thead><tr><th rowspan="2">Cores</th><th rowspan="2">Nodes</th><th colspan="2">Generation Time (s)</th></tr><tr><th>Standard</th><th>Localized</th></tr></thead><tbody><tr><td>16</td><td>1 × 10<sup>6</sup></td><td>0.86</td><td>6.57</td></tr><tr><td>32</td><td>2 × 10<sup>6</sup></td><td>1.55</td><td>13.87</td></tr><tr><td>256</td><td>2 × 10<sup>7</sup></td><td>356.43</td><td>159.98</td></tr><tr><td>2048</td><td>1 × 10<sup>8</sup></td><td>–</td><td>2037.48</td></tr><tr><td>4096</td><td>3 × 10<sup>8</sup></td><td>–</td><td>3299.04</td></tr></tbody></table>

## 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 $n_\Sigma$) is interpreted as a pre-existing displacement discontinuity $[u](x; t)$ 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-varying 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(W, L)$, 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(x_\Sigma)$ (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(r)$ 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(R)$:

$$d(R) = \begin{cases} 2hR, & R < R_c \\ 2hR_c, & R \ge R_c \end{cases} \qquad (7)$$

$$d(R) = \begin{cases} 2hR_c, & R \ge R_c \\ 2hR, & R < R_c \end{cases} \qquad (8)$$

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}(r)$ 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(R)$ defined as:

$$ \tau(R) = \begin{cases} \frac{2R}{V_r}, & R < R_p \\ \frac{2R_p}{V_r}, & R \ge R_p \end{cases} \qquad (9) $$

$$ \tau(R) = \begin{cases} \frac{2R}{V_r}, & R < R_p \\ \frac{2R_p}{V_r}, & R \ge R_p \end{cases} \qquad (10) $$

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 $\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). The basin-like structure was obtained via in situ geophysical measurements by Cushing et al. [47,99]. The mechanical properties within the basin were taken from Hollender et al. [98].
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Regional crustal models proposed by Haslinger et al. [97] and refined by Hollender et al. [98]. Courtesy of Touhami [50].</figcaption></figure>

<table><caption>**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.</caption><thead><tr><th>$z_{\text{top}}$ (m)</th><th>$V_p$ (m/s)</th><th>$V_s$ (m/s)</th><th>$Q_p$ (1)</th><th>$Q_s$ (1)</th></tr></thead><tbody><tr><td>0</td><td>2000</td><td>600</td><td>300</td><td>150</td></tr><tr><td>300</td><td>2400</td><td>1000</td><td>300</td><td>150</td></tr><tr><td>400</td><td>4600</td><td>2700</td><td>300</td><td>150</td></tr><tr><td>1000</td><td>6000</td><td>3200</td><td>300</td><td>150</td></tr><tr><td>2000</td><td>6200</td><td>3200</td><td>300</td><td>150</td></tr></tbody></table>

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><img src="image_13.png" /><figcaption>**Figure 9.** *V*<sub>S</sub> structure of the numerical model of the Argostoli basin. Clips B-B' and C-C' show the location of the principal monitoring points considered (P3, P6, P7, and P262).</figcaption></figure>

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 ≈ 4.5·10<sup>6</sup> elements and ≈9.8–13.5·10<sup>9</sup> (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.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Details of the computational domain surrounding the Argostoli basin. (a) Top view of the domain of interest, including PML layers (flat surface topography). (b) Cut AA' with associated element sizes.</figcaption></figure>

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,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 $ρ(x)$ for the light-gray area shown in Figure 10a and therefore a non-stationary elastic shear modulus $μ(x) = ρ(x)c_S^2(x)$. $c_S(x)$ is the average shear wave velocity described above. $ρ(x)$ is generated via randomField and it is characterized by an anisotropic Von Karman Correlation structure ($ℓ_{C-x} = ℓ_{C-y} = 90$ m, $ℓ_{C-z} = 40$ m) over a domain of ≈7.3 km × 6.5 km × 0.1 km. Its standard coefficient of variation $CoV = σ_ρ/μ_ρ$ 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 $φ_S = 38°$, $δ = 63°$, $λ = 172°$. 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_{W6.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.

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Maximum slip contour across the assumed fault plane. The slip patch was obtained via RIKsurf.</figcaption></figure>

**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 $ρ(x)$ field was generated, featuring a von Karman auto-correlation function with correlations lengths ($ℓ_{C-x}$, $ℓ_{C-y}$, $ℓ_{C-z}$) = (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 ($ℓ_{C-x}$, $ℓ_{C-y}$, $ℓ_{C-z}$).

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.
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Comparison between SEM3D numerical simulations F5 (homogeneous) and O2 (heterogeneous) and the GMPE proposed by [104]. Pseudo Spectral Acceleration (PSA) at natural period *T* = 0.1 s is plotted against hypocentral distance *R*<sub>hyp</sub>. *m*: median ln(PSA); *m* ± *σ*: median ln(PSA) plus/minus standard deviation *σ*. Three horizontal grids of points, located at 0 m G.L., -50 m G.L., and -150 m G.L. were considered. Records were low-pass filtered at 10 Hz before computing PSA.</figcaption></figure>

Both the *M*<sub>W</sub>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* ± *σ* (standard deviation) uncertainty margins. All the considered stations are located in very near field, i.e., around a hypocentral distance *R*<sub>lyb</sub> ≈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.

Figures 14 and 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.
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Horizontal and vertical PGA contour plot from F5 simulation in the surroundings of the Argostoli sedimentary basin. (a) Geometric horizontal average peak ground acceleration $a_{max}$. (b) Vertical peak ground acceleration $a_{max}$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Synthetic acceleration time-histories (a–c) and Fourier's spectra (d–f) filtered at 0–10 Hz, at P262 (–200 m G.L.) obtained in O2 simulation. Comparison between homogeneous (blue) and heterogeneous (red) models.</figcaption></figure>

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.
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Synthetic acceleration time-histories (a–c) and Fourier’s spectra (d–f) filtered at 0–10 Hz at P6 (0 m G.L.) obtained in O2 simulation. Comparison between homogeneous (blue) and heterogeneous (red) models.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Synthetic acceleration time-histories (a–c) and Fourier’s spectra (d–f) filtered at 0–5 Hz at P3. Comparison between O2 (Occigen, blue) and F5 (Fusion, red).</figcaption></figure>

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.

<figure><img src="image_21.png" /><figcaption>**Figure 17.** Contour plots at different time steps, of the particle velocity magnitude (in m/s) in a chunk of the Earth surrounding the sedimentary basin and of the slip rate on the fault plane (in m/s). (a) Snapshots of the particle velocity magnitude (in the crust) and slip rate magnitude (on the fault plane) at *t* = 3 s. (b) Snapshots of the particle velocity magnitude (in the crust) and slip rate magnitude (on the fault plane) at *t* = 4 s. (c) Snapshots of the particle velocity magnitude (in the crust) and slip rate magnitude (on the fault plane) at *t* = 5 s. (d) Snapshots of the particle velocity magnitude (in the crust) and slip rate magnitude (on the fault plane) at *t* = 6 s.</figcaption></figure>

**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.

**Author Contributions:** Conceptualization, F.G. and S.T.; methodology, S.T.; software, F.G., L.d.A.C., R.C., and L.A.; validation, S.T., F.G., and F.L.-C.; formal analysis, S.T. and F.G.; investigation, S.T.; resources, F.G.; data curation, S.T.; writing—original draft preparation, F.G.; writing—review and editing, S.T., F.G., F.L.-C., and D.C.; visualization, F.G. and S.T.; supervision, D.C.; project administration, F.L.-C.; funding acquisition, D.C. and F.L.-C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the French state funding managed by the National Research Agency under the program RNSR Future Investments bearing reference no. ANR-11-RSNR-0022-04 within the SINAPS@ project. The research reported in this paper has been supported in part by the SEISM Paris Saclay Research Institute. This work was granted access to the HPC resources of CINES under the allocations 2018-A0040410444, 2019-A0060410444, 2020-A0070411083, and 2020-A0080410444 made by GENCI (Grand équipement national de calcul intensif). Computations were also performed using HPC resources of the Moulon Mésocentre, hosted by CentraleSupélec and ENS Paris-Saclay.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data generated during the study are available from the corresponding author upon reasonable request.

**Acknowledgments:** The present work is placed within the framework of the French research project SINAPS@ (Séismes Installations Nucléaires, Améliorer et Préserver la Sûreté). This project has been initiated by the members of the institute SEISM after the nuclear disaster of Fukushima (Japan, 2011). The present work is part of the work package Non-Linear Sites Effects and Soil-Structure Interaction [49].

**Conflicts of Interest:** The authors declare no conflicts of interest.

**Abbreviations**

The following abbreviations are used in this manuscript:
<table><tr><td>3-D</td><td>Three-dimensional</td></tr><tr><td>BBS2S</td><td>Broad-Band Source-to-Site</td></tr><tr><td>CEA</td><td>Commissariat à l'énergie atomique et aux énergies alternatives</td></tr><tr><td>CFL</td><td>Courant-Friedrichs-Lewy</td></tr><tr><td>CPU</td><td>Central Process Unit</td></tr><tr><td>CoV</td><td>Coefficient of Variation</td></tr><tr><td>DEM</td><td>Digital Elevation Model</td></tr><tr><td>E2VP</td><td>EUROSEISTEST Verification and Validation Project</td></tr><tr><td>DG</td><td>Discontinuous Galerkin</td></tr><tr><td>FDM</td><td>Finite-Difference Method</td></tr><tr><td>FEM</td><td>Finite-Element Method</td></tr><tr><td>G.L.</td><td>Ground Level</td></tr><tr><td>GLL</td><td>Gauss-Lobatto-Legendre</td></tr><tr><td>GMPE</td><td>Ground Motion Prediction Equation</td></tr><tr><td>GPS</td><td>Global Positioning System</td></tr><tr><td>GPU</td><td>Graphic Process Unit</td></tr><tr><td>I/O</td><td>Input/Output</td></tr><tr><td>INSAR</td><td>Interferometric synthetic aperture radar</td></tr><tr><td>IPGP</td><td>Institut de Physique du Globe de Paris</td></tr><tr><td>MPI</td><td>Message Passing Interface</td></tr><tr><td>NERA</td><td>Network of European Research Infrastructures for Earthquake Risk Assessment and Mitigation</td></tr><tr><td>PGA</td><td>Peak Ground Acceleration</td></tr><tr><td>PML</td><td>Perfectly Matched Layer</td></tr><tr><td>PSA</td><td>Pseudo Spectral Acceleration</td></tr><tr><td>RIK</td><td>Ruiz Integral Kinematic</td></tr><tr><td>SEM</td><td>Spectral Element Method</td></tr></table>

## References

1. Aki, K.; Richards, P.G. *Quantitative Seismology*; Books in Geology; W.H Freeman and Company: San Francisco, CA, USA, 2002; Volume I-II.

2. Aki, K.; Chouet, B. Origin of coda waves: Source, attenuation, and scattering effects. *J. Geophys. Res.* **1975**, *80*, 3322–3342. [CrossRef]

3. Nakano, K.; Matsushima, S.; Kawase, H. Statistical Properties of Strong Ground Motions from the Generalized Spectral Inversion of Data Observed by K-NET, KiK-net, and the JMA Shindokei Network in Japan. *Bull. Seismol. Soc. Am.* **2015**, *105*, 2662–2680. [CrossRef]

4. Mai, P.M.; Beroza, G.C. A spatial random field model to characterize complexity in earthquake slip. *J. Geophys. Res. Solid Earth* **2002**, *107*, ESE 10-1–ESE 10-21. [CrossRef]

5. Aochi, H. Input Ground Motion Calculation Based on Dynamic Rupture Modeling on Segmented Faults. In Proceedings of the International Workshop on Best Practices in Physics-based Fault Rupture Models for Seismic Hazard Assessment of Nuclear Installations, Vienna, Austria, 18–20 November 2015.

6. Chaljub, E.; Maufroy, E.; Moczo, P.; Kristek, J.; Hollender, F.; Bard, P.Y.; Priolo, E.; Klin, P.; De Martin, F.; Zhang, Z.; et al. 3-D numerical simulations of earthquake ground motion in sedimentary basins: Testing accuracy through stringent models. *Geophys. J. Int.* **2015**, *201*, 90–111. [CrossRef]

7. Moczo, P.; Kristek, J.; Bard, P.Y.; Stripajová, S.; Hollender, F.; Chovanová, Z.; Kristeková, M.; Sicilia, D. Key structural parameters affecting earthquake ground motion in 2D and 3D sedimentary structures. *Bull. Earthq. Eng.* **2018**, *16*, 2421–2450. [CrossRef]

8. De Martin, F.; Chaljub, E.; Thierry, P.; Sochala, P.; Dupros, F.; Maufroy, E.; Hadri, B.; Benaichouche, A.; Hollender, F. Influential parameters on 3-D synthetic ground motions in a sedimentary basin derived from global sensitivity analysis. *Geophys. J. Int.* **2021**, *227*, 1795–1817. [CrossRef]

9. Faccioli, E.; Tagliani, A.; Paolucci, R. Effects of the Wave-Propagation in Random Earth's Media on the Seismic Radiation Spectra. In Proceedings of the 4th International Conference on Soil Dynamics and Earthquake Engineering, Mexico City, Mexico, 23–26 October 1989; Cakmak, A.S., Herrera, I., Eds.; Computational Mechanics Publications: Berlin, Germany, 1989; pp. 61–75.

10. Moczo, P.; Kristek, J.; Vavryčuk, V.; Archuleta, R.J.; Halada, L. 3D heterogeneous staggered-grid finite-difference modeling of seismic motion with volume harmonic and arithmetic averaging of elastic moduli and densities. *Bull. Seismol. Soc. Am.* **2002**, *92*, 3042–3066. [CrossRef]

11. Scalise, M.; Pitarka, A.; Louie, J.N.; Smith, K.D. Effect of Random 3D Correlated Velocity Perturbations on Numerical Modeling of Ground Motion from the Source Physics Experiment. *Bull. Seismol. Soc. Am.* **2021**, *111*, 139–156. [CrossRef]
12. Ichimura, T.; Fujita, K.; Quinay, P.E.B.; Maddegedara, L.; Hori, M.; Tanaka, S.; Shizawa, Y.; Kobayashi, H.; Minami, K. Implicit nonlinear wave simulation with 1.08 T DOF and 0.270 T unstructured finite elements to enhance comprehensive earthquake simulation. In Proceedings of the High Performance Computing, Networking, Storage and Analysis, Austin, TX, USA, 15–20 November 2015; pp. 1–12. [CrossRef]

13. Ichimura, T.; Fujita, K.; Tanaka, S.; Hori, M.; Lalith, M.; Shizawa, Y.; Kobayashi, H. Physics-Based Urban Earthquake Simulation Enhanced by 10.7 BlnDOF x 30 K Time-Step Unstructured FE Non-Linear Seismic Wave Simulation. In Proceedings of the SC '14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, New Orleans, LA, USA, 16–21 November 2014; pp. 15–26. [CrossRef]

14. Restrepo, D.; Bielak, J.; Serrano, R.; Gomez, J.; Jaramillo, J. Effects of realistic topography on the ground motion of the Colombian Andes—A case study at the Aburrá Valley, Antioquia. *Geophys. J. Int.* **2016**, *204*, 1801–1816. [CrossRef]

15. Ichimura, T.; Fujita, K.; Yoshiyuki, A.; Quinay, P.E.; Hori, M.; Sakanoue, T. Performance Enhancement of Three-Dimensional Soil Structure Model via Optimization for Estimating Seismic Behavior of Buried Pipelines. *J. Earthq. Tsunami* **2017**, *11*, 1750019. [CrossRef]

16. Zuchowski, L.; Brun, M.; De Martin, F. Co-simulation coupling spectral/finite elements for 3D soil/structure interaction problems. *C. R. Méc.* **2018**, *346*, 408–422. [CrossRef]

17. Brun, M.; De Martin, F.; Richart, N. Hybrid asynchronous SEM/FEM co-simulation for seismic nonlinear analysis of concrete gravity dams. *Comput. Struct.* **2021**, *245*, 106459. [CrossRef]

18. Paolucci, R.; Mazzieri, I.; Smerzini, C.; Stupazzini, M. Physics-Based Earthquake Ground Shaking Scenarios in Large Urban Areas. In Perspectives on European Earthquake Engineering and Seismology; Ansal, A., Ed.; Springer: Cham, Switzerland, 2014; Volume 34, pp. 331–359.

19. Carrington, L.; Komatitsch, D.; Laurenzano, M.; Tikir, M.; Michéa, D.; Le Goff, N.; Snavely, A.; Tromp, J. High-frequency simulations of global seismic wave propagation using SPECFEM3D_GLOBE on 62 thousand processor cores. In Proceedings of the SC'08 ACM/IEEE conference on Supercomputing, Austin, TX, USA, 15–21 November 2008; pp. 1–11. [CrossRef]

20. Dupros, F.; De Martin, F.; Foerster, E.; Komatitsch, D.; Roman, J. High-performance finite-element simulations of seismic wave propagation in three-dimensional nonlinear inelastic geological media. *Parallel Comput.* **2010**, *36*, 308–325. [CrossRef]

21. Paolucci, R.; Mazzieri, I.; Smerzini, C. Anatomy of strong ground motion: Near-source records and 3D physics-based numerical simulations of the Mw 6.0 May 29 2012 Po Plain earthquake, Italy. *Geophys. J. Int.* **2015**, *203*, 2001–2020. [CrossRef]

22. Smerzini, C.; Pitilakis, K.; Hashemi, K. Evaluation of earthquake ground motion and site effects in the Thessaloniki urban area by 3D finite-fault numerical simulations. *Bull. Earthq. Eng.* **2017**, *15*, 787–812. [CrossRef]

23. Gatti, F.; Lopez-Caballero, F.; Clouteau, D.; Paolucci, R. On the effect of the 3-D regional geology on the seismic design of critical structures: The case of the Kashiwazaki-Kariwa Nuclear Power Plant. *Geophys. J. Int.* **2018**, *213*, 1073–1092. [CrossRef]

24. Poursartip, B.; Fathi, A.; Tassoulas, J.L. Large-scale simulation of seismic wave motion: A review. *Soil Dyn. Earthq. Eng.* **2020**, *129*, 105909. [CrossRef]

25. Paolucci, R.; Infantino, M.; Mazzieri, I.; Özcebe, A.G.; Smerzini, C.; Stupazzini, M. 3D Physics-Based Numerical Simulations: Advantages and Current Limitations of a New Frontier to Earthquake Ground Motion Prediction. The Istanbul Case Study. In Proceedings of the Recent Advances in Earthquake Engineering in Europe: 16th European Conference on Earthquake Engineering, Thessaloniki, Greece, 18–21 June 2018; Volume 46, p. 203.

26. Dujardin, A.; Causse, M.; Berge-Thierry, C.; Hollender, F. Radiation Patterns Control the Near-Source Ground-Motion Saturation Effect. *Bull. Seismol. Soc. Am.* **2018**, *108*, 3398. [CrossRef]

27. De Martin, F. Verification of a Spectral-Element Method Code for the Southern California Earthquake Center LOH.3 Viscoelastic Case. *Bull. Seismol. Soc. Am.* **2011**, *101*, 2855–2865. [CrossRef]

28. Cupillard, P.; Delavaud, E.; Burgos, G.; Festa, G.; Vilotte, J.P.; Capdeville, Y.; Montagner, J.P. RegSEM: A versatile code based on the spectral element method to compute seismic wave propagation at the regional scale. *Geophys. J. Int.* **2012**, *188*, 1203–1220. [CrossRef]

29. Quinay, P.E.B.; Ichimura, T.; Hori, M.; Nishida, A.; Yoshimura, S. Seismic Structural Response Estimates of a Fault-Structure System Model with Fine Resolution Using Multiscale Analysis with Parallel Simulation of Seismic-Wave Propagation. *Bull. Seismol. Soc. Am.* **2013**, *103*, 2094–2110. [CrossRef]

30. He, C.H.; Wang, J.T.; Zhang, C.H. Nonlinear Spectral-Element Method for 3D Seismic-Wave Propagation. *Bull. Seismol. Soc. Am.* **2016**, *106*, 1074–1087. [CrossRef]

31. Ichimura, T.; Agata, R.; Hori, T.; Hirahara, K.; Hashimoto, C.; Hori, M.; Fukahata, Y. An elastic/viscoelastic finite element analysis method for crustal deformation using a 3-D island-scale high-fidelity model. *Geophys. J. Int.* **2016**, *206*, 114–129. [CrossRef]

32. Tsuboi, S.; Ando, K.; Miyoshi, T.; Peter, D.; Komatitsch, D.; Tromp, J. A 1.8 trillion degrees-of-freedom, 1.24 petaflops global seismic wave simulation on the K computer. *Int. J. High Perform. Comput. Appl.* **2016**, *30*, 411–422. [CrossRef]

33. Baltaji, O.; Numanoglu, O.A.; Veeraraghavan, S.; Hashash, Y.M.; Coleman, J.L.; Bolisetti, C. Non-linear Time Domain Site Response and Soil Structure Analyses for Nuclear Facilities using MASTODON. In Proceedings of the 24th Conference on Structural Mechanics in Reactor Technology, Busan, Korea, 20–25 August 2017.
56. Maday, Y.; Patera, A.T.; Ronquist, E.M. *A Well-Posed Optimal Spectral Element Approximation for the Stokes Problem*; Technical Report, ICASE Report, no. 87-48.; NASA Contractor Report, NASA CR-178343; NASA: Washington, DC, USA, 1987.

57. Maday, Y.; Patera, A.; Rønquist, E. Optimal Legendre spectral element methods for the multi-dimensional Stokes problem. *SIAM J. Num. Anal.* **1989**.

58. Seriani, G. 3-D large-scale wave propagation modeling by spectral element method on Cray T3E multiprocessor. *Comput. Methods Appl. Mech. Eng.* **1998**, *164*, 235–247. [CrossRef]

59. Komatitsch, D.; Tromp, J. Introduction to the spectral element method for three-dimensional seismic wave propagation. *Geophys. J. Int.* **1999**, *139*, 806–822. [CrossRef]

60. Delavaud, E. *Simulation Numérique de la Propagation D'ondes en Milieu Géologique Complexe: Application à L'évaluation de la Réponse Sismique du Bassin de Caracas (Venezuela)*; Institut de Physique du Globe de Paris: Paris, France, 2007.

61. Mazzieri, I.; Stupazzini, M.; Guidotti, R.; Smerzini, C. SPEED: SPpectral Elements in Elastodynamics with Discontinuous Galerkin: A non-conforming approach for 3D multi-scale problems. *Int. J. Numer. Methods Eng.* **2013**, *95*, 991–1010. [CrossRef]

62. Simo, J.; Tarnow, N.; Wong, K. Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics. *Comput. Methods Appl. Mech. Eng.* **1992**, *100*, 63–116. [CrossRef]

63. Courant, R.; Friedrichs, K.; Lewy, H. On the partial difference equations of mathematical physics. *IBM J. Res. Dev.* **1967**, *11*, 215–234. [CrossRef]

64. Cohen, G.C. *High-Order Numerical Methods for Transient Wave Equations*; Springer: Berlin, Germany, 2010; Volume 63, p. 49. [CrossRef]

65. Sevilla, R.; Cottereau, R. Influence of periodically fluctuating material parameters on the stability of explicit high-order spectral element methods. *J. Comput. Phys.* **2018**, *373*, 304–323. [CrossRef]

66. Cottereau, R.; Sevilla, R. Stability of an explicit high-order spectral element method for acoustics in heterogeneous media based on local element stability criteria. *Int. J. Numer. Methods Eng.* **2018**, *116*, 223–245. [CrossRef]

67. Seriani, G.; Oliveira, S.P. Dispersion analysis of spectral element methods for elastic wave propagation. *Wave Motion* **2008**, *45*, 729–744. [CrossRef]

68. Berenger, J.P. A perfectly matched layer for the absorption of electromagnetic waves. *J. Comput. Phys.* **1994**, *114*, 185–200. [CrossRef]

69. Festa, G.; Vilotte, J.P. The Newmark scheme as velocity-stress time-staggering: An efficient PML implementation for spectral element simulations of elastodynamics. *Geophys. J. Int.* **2005**, *161*, 789–812. [CrossRef]

70. Chaljub, E.; Komatitsch, D.; Vilotte, J.P.; Capdeville, Y.; Valette, B.; Festa, G. Spectral-element analysis in seismology. *Adv. Geophys.* **2007**, *48*, 365–419.

71. Chabot, S.; Glinsky, N.; Diego Mercerat, E.; Bonilla, L.F. A High-Order Discontinuous Galerkin Method for Coupled Wave Propagation in 1D Elastoplastic Heterogeneous Media. *J. Theor. Comput. Acoust.* **2018**, *26*, 1850043. [CrossRef]

72. Sornet, G.; Jubertie, S.; Dupros, F.; De Martin, F.; Limet, S. Performance Analysis of SIMD Vectorization of High-Order Finite-Element Kernels. In Proceedings of the 2018 International Conference on High Performance Computing Simulation (HPCS), Orleans, France, 16–20 July 2018; pp. 423–430. [CrossRef]

73. Tang, H.; Byna, S.; Petersson, N.A.; McCallen, D. Tuning Parallel Data Compression and I/O for Large-scale Earthquake Simulation. In Proceedings of the 2021 IEEE International Conference on Big Data (Big Data), Orlando, FL, USA, 15–18 December 2021; pp. 2992–2997. [CrossRef]

74. Göddeke, D.; Komatitsch, D.; Möller, M. Finite and Spectral Element Methods on Unstructured Grids for Flow and Wave Propagation Problems. In *Numerical Computations with GPUs*; Springer International Publishing: Berlin, Germany, 2014; pp. 183–206. [CrossRef]

75. Markall, G.R.; Slemmer, A.; Ham, D.A.; Kelly, P.H.J.; Cantwell, C.D.; Sherwin, S.J. Finite element assembly strategies on multi-core and many-core architectures. *Int. J. Numer. Methods Fluids* **2013**, *71*, 80–97. [CrossRef]

76. Lombard, B.; Piraux, J. Numerical modeling of transient two-dimensional viscoelastic waves. *J. Comput. Phys.* **2011**, *230*, 6099–6114. [CrossRef]

77. Fu, H.; He, C.; Chen, B.; Yin, Z.; Zhang, Z.; Zhang, W.; Zhang, T.; Xue, W.; Liu, W.; Yin, W.; et al. 18.9Pflopss Nonlinear Earthquake Simulation on Sunway TaihuLight: Enabling Depiction of 18-Hz and 8-meter Scenarios. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, New York, NY, USA, 12–17 November 2017; pp. 2:1–2:12. [CrossRef]

78. Chaljub, E.; Moczo, P.; Tsuno, S.; Bard, P.Y.; Kristek, J.; Käser, M.; Stupazzini, M.; Kristekova, M. Quantitative comparison of four numerical predictions of 3D ground motion in the Grenoble Valley, France. *Bull. Seismol. Soc. Am.* **2010**, *100*, 1427–1455. [CrossRef]

79. Bielak, J.; Graves, R.W.; Olsen, K.B.; Taborda, R.; Ramirez-Guzman, L.; Day, S.M.; Ely, G.; Roten, D.; Jordan, T.; Maechling, P.; et al. The ShakeOut earthquake scenario: Verification of three simulation sets. *Geophys. J. Int.* **2010**, *180*, 375–404. [CrossRef]

80. Komatitsch, D.; Göddeke, D.; Erlebacher, G.; Michéa, D. Modeling the propagation of elastic waves using spectral elements on a cluster of 192 GPUs. *Comput. Sci. Res. Dev.* **2010**, *25*, 75–82. [CrossRef]
81. Cui, Y.; Olsen, K.B.; Jordan, T.H.; Lee, K.; Zhou, J.; Small, P.; Roten, D.; Ely, G.; Panda, D.K.; Chourasia, A.; et al. Scalable earthquake simulation on petascale supercomputers. In Proceedings of the SC '10: Proceedings of the 2010 ACM/IEEE International Conference for High Performance Computing, Networking, Storage and Analysis, New Orleans, LA, USA, 13–19 November 2010; pp. 1–20.

82. Rietmann, M.; Messmer, P.; Nissen-Meyer, T.; Peter, D.; Basini, P.; Komatitsch, D.; Schenk, O.; Tromp, J.; Boschi, L.; Giardini, D. Forward and adjoint simulations of seismic wave propagation on emerging large-scale GPU architectures. In Proceedings of the SC '12: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, Salt Lake City, UT, USA, 10–16 November 2012; pp. 1–11. [CrossRef]

83. Taborda, R.; Bielak, J. Ground-Motion Simulation and Validation of the 2008 Chino Hills, California, Earthquake. *Bull. Seismol. Soc. Am.* **2013**, *103*, 131–156. [CrossRef]

84. Heinecke, A.; Breuer, A.; Rettenberger, S.; Bader, M.; Gabriel, A.; Pelties, C.; Bode, A.; Barth, W.; Liao, X.; Vaidyanathan, K.; et al. Petascale High Order Dynamic Rupture Earthquake Simulations on Heterogeneous Supercomputers. In Proceedings of the SC '14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, New Orleans, LA, USA, 16–21 November 2014; pp. 3–14. [CrossRef]

85. Ichimura, T.; Fujita, K.; Koyama, K.; Kusakabe, R.; Kikuchi, Y.; Hori, T.; Hori, M.; Maddegedara, L.; Ohi, N.; Nishiki, T.; et al. 152K-Computer-Node Parallel Scalable Implicit Solver for Dynamic Nonlinear Earthquake Simulation. In Proceedings of the International Conference on High Performance Computing in Asia-Pacific Region, New York, NY, USA, 12–14 January 2022; pp. 18–29. [CrossRef]

86. Duru, K.; Rannabauer, L.; Gabriel, A.A.; Ling, O.K.A.; Igel, H.; Bader, M. A stable discontinuous Galerkin method for linear elastodynamics in 3D geometrically complex elastic solids using physics based numerical fluxes. *Comput. Methods Appl. Mech. Eng.* **2022**, *389*, 114386. [CrossRef]

87. Bachem, A. *PRACE*—*Strategy for HPC in Europe*; Technical Report; Forschungszentrum Jülich: Jülich, Germany, 2008.

88. Camata, J.J.; Coutinho, A.L.G.A. Parallel implementation and performance analysis of a linear octree finite element mesh generation scheme. *Concurr. Comput. Pract. Exp.* **2013**, *25*, 826–842. [CrossRef]

89. Shinozuka, M.; Deodatis, G. Simulation of stochastic processes by spectral representation. *Appl. Mech. Rev.* **1991**, *44*, 191–204. [CrossRef]

90. Shinozuka, M.; Deodatis, G. Simulation of multi-dimensional Gaussian stochastic fields by spectral representation. *Appl. Mech. Rev.* **1996**, *49*, 29–53. [CrossRef]

91. Crempien, J.G.F.; Archuleta, R.J. UCSB Method for Simulation of Broadband Ground Motion from Kinematic Earthquake Sources. *Seismol. Res. Lett.* **2015**, *86*, 61–67. [CrossRef]

92. Ruiz, J.A.; Baumont, D.; Bernard, P.; Berge-Thierry, C. Modelling directivity of strong ground motion with a fractal, k-2, kinematic source model. *Geophys. J. Int.* **2011**, *186*, 226. [CrossRef]

93. Eshelby, J.D. The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems. *Proc. R. Soc. Lond. Math. Phys. Eng. Sci.* **1957**, *241*, 376–396. [CrossRef]

94. Faccioli, E.; Maggio, F.; Paolucci, R.; Quarteroni, A. 2D and 3D elastic wave propagation by a pseudo-spectral domain decomposition method. *J. Seismol.* **1997**, *1*, 237–251. [CrossRef]

95. Louvari, E.; Kiratzi, A.; Papazachos, B. The Cephalonia Transform Fault and its extension to western Lefkada Island (Greece). *Tectonophysics* **1999**, *308*, 223–236. [CrossRef]

96. Saltogianni, V.; Moschas, F.; Stiros, S. The 2014 Cephalonia Earthquakes: Finite Fault Modeling, Fault Segmentation, Shear and Thrusting at the NW Aegean Arc (Greece). *Pure Appl. Geophys.* **2018**, *175*, 4145–4164. [CrossRef]

97. Haslinger, F.; Kissling, E.; Ansorge, J.; Hatzfeld, D.; Papadimitriou, E.; Karakostas, V.; Makropoulos, K.; Kahle, H.G.; Peter, Y. 3D crustal structure from local earthquake tomography around the Gulf of Arta (Ionian region, NW Greece). *Tectonophysics* **1999**, *304*, 201–218. [CrossRef]

98. Hollender, F.; Perron, V.; Imtiaz, A.; Svay, A.; Mariscal, A.; Bard, P.; Cottereau, R.; Lopez-Caballero, F.; Cushing, M.; Theodoulidis, N.; et al. *A Deux Pas du Repaire du Cyclope D'Ulysse: La Campagne Post-Sismique et le Démarrage du Réseau Accélérométrique SINAPS@ sur l'île de Céphalonie*; Colloque AFPS: Paris, France, 2015.

99. Cushing, E.M.; Hollender, F.; Moiriat, D.; Guyonnet-Benaize, C.; Theodoulidis, N.; Pons-Branchu, E.; Sépulcre, S.; Bard, P.Y.; Cornou, C.; Dechamp, A.; et al. Building a three dimensional model of the active Plio-Quaternary basin of Argostoli (Cephalonia Island, Greece): An integrated geophysical and geological approach. *Eng. Geol.* **2020**, *265*, 105441. [CrossRef]

100. Sokos, E.; Zahradnik, J.; Gallovic, F.; Serpetsidaki, A.; Plicka, V.; Kiratzi, A. Asperity break after 12 years: The Mw6.4 2015 Lefkada (Greece) earthquake. *Geophys. Res. Lett.* **2016**, *43*, 6137–6145. [CrossRef]

101. Perron, V.; Hollender, F.; Froment, B.; Gelis, C.; Cushing, E.M.; Bard, P.Y.; Cultrera, G. Can broad-band earthquake site responses be predicted by the ambient noise spectral ratio? Insight from observations at two sedimentary basins. *Geophys. J. Int.* **2018**, *215*, 1442–1454. [CrossRef]

102. Maufroy, E.; Chaljub, E.; Hollender, F.; Kristek, J.; Moczo, P.; Klin, P.; Priolo, E.; Iwaki, A.; Iwata, T.; Etienne, V.; et al. Earthquake Ground Motion in the Mygdonian Basin, Greece: The E2VP Verification and Validation of 3D Numerical Simulation up to 4 Hz. *Bull. Seismol. Soc. Am.* **2015**, *105*, 1398–1418. [CrossRef]
103. Maufroy, E.; Chaljub, E.; Theodoulidis, N.P.; Roumelioti, Z.; Hollender, F.; Bard, P.Y.; de Martin, F.; Guyonnet-Benaize, C.; Margerin, L. Source-Related Variability of Site Response in the Mygdonian Basin (Greece) from Accelerometric Recordings and 3D Numerical Simulations. *Bull. Seismol. Soc. Am.* **2017**, *107*, 787–808. [CrossRef]

104. Akkar, S.; Sandıkkaya, M.A.; Bommer, J.J. Empirical ground-motion models for point- and extended-source crustal earthquake scenarios in Europe and the Middle East. *Bull. Earthq. Eng.* **2014**, *12*, 359–387. [CrossRef]

105. Berge-Thierry, C.; Voldoire, F.; Ragueneau, F.; Lopez-Caballero, F.; Le Maoult, A. Main Achievements of the Multidisciplinary SINAPS Research Project: Towards an Integrated Approach to Perform Seismic Safety Analysis of Nuclear Facilities. *Pure Appl. Geophys.* **2019**, *177*, 2299–2351. [CrossRef]

106. Jussila, V.; Fälth, B.; Mäntyniemi, P.; Voss, P.H.; Lund, B.; Fülöp, L. Application of a Hybrid Modeling Method for Generating Synthetic Ground Motions in Fennoscandia, Northern Europe. *Bull. Seismol. Soc. Am.* **2021**, *111*, 2507–2526. [CrossRef]

107. Oliveira, S.P. Error Analysis of Chebyshev Spectral Element Methods for the Acoustic Wave Equation in Heterogeneous Media. *J. Theor. Comput. Acoust.* **2018**, *26*, 1850035. [CrossRef]