# Synthetic ground motions in heterogeneous geologies from various sources: the HEMEW S -3D database

Fanny Lehmann, Filippo Gatti, Michaël Bertin, Didier Clouteau

► To cite this version:

Fanny Lehmann, Filippo Gatti, Michaël Bertin, Didier Clouteau. Synthetic ground motions in heterogeneous geologies from various sources: the HEMEW S -3D database. Earth System Science Data, 2024, 16 (9), pp.3949-3972. 10.5194/essd-16-3949-2024. hal-04699436

HAL Id: hal-04699436

https://hal.science/hal-04699436v1

Submitted on 31 Mar 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Synthetic ground motions in heterogeneous geologies from various sources: the HEMEW<sup>S</sup>-3D database

Fanny Lehmann<sup>1,2</sup>, Filippo Gatti<sup>2</sup>, Michaël Bertin<sup>1</sup>, and Didier Clouteau<sup>2</sup>

<sup>1</sup>CEA, DAM, DIF, 91297 Arpajon, France

<sup>2</sup>LMPS – Laboratoire de Mécanique Paris-Saclay, Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, Gif-sur-Yvette, France

**Correspondence:** Fanny Lehmann (fanny.lehmann@centralesupelec.fr)

Received: 14 November 2023 – Discussion started: 9 January 2024

Revised: 12 July 2024 – Accepted: 18 July 2024 – Published: 3 September 2024

**Abstract.** The ever-improving performances of physics-based simulations and the rapid developments of deep learning are offering new perspectives to study earthquake-induced ground motion. Due to the large amount of data required to train deep neural networks, applications have so far been limited to recorded data or two-dimensional (2D) simulations. To bridge the gap between deep learning and high-fidelity numerical simulations, this work introduces a new database of physics-based earthquake simulations.

The Heterogeneous Materials and Elastic Waves with Source variability in 3D (HEMEW<sup>S</sup>-3D) database comprises 30 000 simulations of elastic wave propagation in 3D geological domains. Each domain is parametrized by a different geological model built from a random arrangement of layers augmented by random fields that represent heterogeneities. Elastic waves originate from a randomly located pointwise source parametrized by a random moment tensor. For each simulation, ground motion is synthesized at the surface by a grid of virtual sensors. The high frequency of waveforms ($f_{max}$ = 5 Hz) allows for extensive analyses of surface ground motion.

Existing and foreseen applications range from statistical analyses of the ground motion variability and machine learning methods on geological models to deep-learning-based predictions of ground motion that depend on 3D heterogeneous geologies and source properties. Data are available at https://doi.org/10.57745/LAI6YU (Lehmann, 2023).

## 1 Introduction

Deep learning has a long tradition in seismology thanks to large networks of sensors recording earthquakes worldwide. Applications are extremely diverse in terms of methods, data, and scientific goals (see, for example, Mousavi and Beroza, 2023, for a review). Detecting earthquakes and distinguishing them from other events such as explosions, quarry blasts, or seismic noise are the most common applications of deep learning in seismology (Mousavi and Beroza, 2023). A wide variety of methods is also devoted to characterizing earthquakes from ground motion recordings – for instance, to estimate source mechanisms, earthquake location, and magnitude. The rapid improvements in deep learning in the last few years have even enabled its use in operational frame-

works, thereby providing real-time predictions of earthquake parameters (Zhu et al., 2022).

However, all those methods rely on databases of seismic waveforms. While there exist several curated databases of recorded ground motion, they are sparse in regions with low to moderate seismicity or poor instrumental coverage (Bahrampouri et al., 2021; Michelini et al., 2021; Mousavi et al., 2019). In those cases, numerical simulations are a great opportunity to complement existing databases. Simulations rely on computational schemes to solve the wave propagation equations from the earthquake’s source to the Earth’s surface and provide synthetic waveforms at any spatial point of the simulation domain. Results of three-dimensional (3D) physics-based simulations have been compiled for several past earthquakes in, for example, the BB-SPEEDset dataset (Paolucci et al., 2021) and the Southern California Earth-
quake Center (SCEC) Broadband Platform (Maechling et al., 2015), but the number of simulations is not appropriate for machine learning approaches.

In fact, physics-based simulations show several limitations. Firstly, they require a detailed description of the ground properties that define the physical behaviour of the waves propagating in the Earth. Especially, ground properties should be given in the form of 3D geological models since 3D features have crucial effects that are not accounted for in 2D settings (e.g. sedimentary basins leading to site effects) (Moczo et al., 2018; Smerzini et al., 2011; Zhu et al., 2020). Since extensive geophysical investigations are needed to obtain 3D geological models, they are rare, and when they exist, they are still limited by epistemic uncertainties. Therefore, when trying to reproduce an earthquake with physics-based numerical simulations, uncertainties can be represented by random heterogeneities added to the reference model to introduce variability (Chaljub et al., 2021; Lehmann et al., 2022).

Quantifying the effects of 3D geological features is made more difficult by the second limitation of physics-based simulations, which is their high computational cost, especially when dealing with high frequencies and large spatial domains. Despite relying on high-performance computing (HPC) frameworks, seismic waves propagation simulations can reach tens to hundreds of thousands of equivalent core hours (Fu et al., 2017; Heinecke et al., 2014; Poursar tip et al., 2020). Since computational costs prevent statistical studies on synthetic waveforms due to a limited number of simulations per geological model, deep learning represents a promising alternative to obtain waveforms.

When predicting the surface ground motion generated by an earthquake, it is important to obtain time series that describe the temporal evolution of shaking and not only scalar features (such as peak ground acceleration and cumulative absolute velocity) that give useful but limited information. Physics-informed neural networks (PINNs; Raissi et al., 2019) successfully solved the wave equation (Ding et al., 2023; Karimpouli and Tahmasebi, 2020; Moseley et al., 2020; Rasht-Behesht et al., 2022; Ren et al., 2024; Song et al., 2023; Wu et al., 2023). However, applications are mainly limited to 2D domains, and models cannot extrapolate to different geological configuration to the one used in the training phase. Alternatively, generative methods have been used to enhance existing numerical simulations by increasing their spatial resolution (e.g. Gadylishin et al., 2021) or their frequency content (e.g. Gatti and Clouteau, 2020).

The recent emergence of scientific machine learning (SciML) is offering a new paradigm for the prediction of physics-based ground motion parametrized by 3D ground properties and source parameters, with intrinsic ability to generalize to various resolutions and geological configurations. SciML has led to significant scientific developments in communities with large, reliable, and freely available databases. For instance, in numerical weather prediction,

Bonev et al. (2023) and Pathak et al. (2022) took advantage of the ERA5 dataset provided by the European Centre for Medium-Range Weather Forecasts (Hersbach et al., 2020). In seismology, Mousavi and Beroza (2023) pointed out that “the limitations on training data and generalization are the main challenges in solving inverse and forward problems using supervised [deep neural networks].”

In this work, we describe the first open database of seismic simulations associated with 3D heterogeneous geological models. The HEMEWS³-3D (HEterogeneous Materials and Elastic Waves with Source variability in 3D) database contains 30 000 high-fidelity simulations in 3D domains of 9.6 × 9.6 × 9.6 km³ in size. This represents a challenging computational task, accounting for 9 × 10⁵ core hours and 4.4 MWh in total. Ground motion was synthesized at the surface of the simulation domain for 8 s on a grid of 32 × 32 virtual sensors. Data are available at https://doi.org/10.57745/LAI6YU (Lehmann, 2023).

In the following text, Sect. 2 provides an overview of existing datasets in related fields; Sect. 3 describes the geological models, sources, and surface wavefields in the database; Sect. 4 analyses physical characteristics; Sect. 5 illustrates applications; and Sect. 6 discusses limitations and perspectives.

## 2 Related work

Datasets of recorded ground motion have enabled major deep learning applications in seismology, but they have several limitations in data-scarce regions. In this section, we focus on datasets with 2D or 3D data used in geophysics and seismology with SciML applications. Due to the mathematical similarities between wave propagation and fluid flow (both are governed by hyperbolic equations), related studies are reviewed beyond the field of seismology. This highlights the challenges of high-fidelity numerical simulations for deep learning applications.

### 2.1 Datasets of ground motion simulations

Ground motion simulations of past earthquakes have been collected in databases for model verification, characterization of complex near-field conditions, and machine learning purposes. The BB-SPEEDset provides 3D simulations of 16 earthquake scenarios in various regions of the globe (Paolucci et al., 2021). The SCEC Broadband Platform also provides simulations of 17 past earthquakes with different source models (Maechling et al., 2015). Other studies focus on specific regions, such as the Hayward Fault (Petrone et al., 2021) and Türkiye (Altindal and Askan, 2022). The limited number of scenarios makes the abovementioned databases less suitable for machine learning purposes, where dataset variability and size are crucial factors. In addition, although they provide some variability on the source, those datasets do
<table><caption>**Table 1.** Summary of datasets providing geological models and seismic wavefields. Domain: size of the physical domain with the number of grid points given in parenthesis, in the order width, depth for 2D datasets and width, length, depth for 3D datasets. The dimension of seismic wavefields is given in the format receivers along width, time steps for 2D datasets and receivers along width, receivers along length, time steps for 3D datasets.</caption><thead><tr><th rowspan="2">Dataset</th><th rowspan="2">Train/test</th><th colspan="3">Geological models</th><th colspan="3">Seismic wavefields</th></tr><tr><th>Domain</th><th>Values of V<sub>s</sub></th><th>Construction</th><th>Dimensions</th><th>Equation</th><th>Source</th></tr></thead><tbody><tr><td>Noddverse (Jessell et al., 2022)</td><td>1 000 000/-</td><td>4 × 4 × 4 km<sup>3</sup> (200 × 200 × 200)</td><td>categorical</td><td>succession of geological events</td><td>not provided</td><td>not provided</td><td>not provided</td></tr><tr><td>OpenFWI (Deng et al., 2022)</td><td>408 000/62 000</td><td>0.7 × 0.7 km<sup>2</sup> (70 × 70)</td><td>1500; 4500 m s<sup>-1</sup></td><td>mathematical and recorded images and geological faults</td><td>70 × 1000</td><td>2D acoustic</td><td>5 fixed sources at the surface</td></tr><tr><td>OpenFWI Kimberlina CO<sub>2</sub> (Deng et al., 2022)</td><td>15 000/4000</td><td>4 × 1.4 km<sup>2</sup> (401 × 141)</td><td>1200; 3600 m s<sup>-1</sup></td><td>from real data</td><td>101 × 1251</td><td>2D acoustic</td><td>9 fixed sources at the surface</td></tr><tr><td>OpenFWI Kimberlina 3D (Deng et al., 2022)</td><td>1600/163</td><td>4 × 4 × 3.5 km<sup>3</sup> (400 × 400 × 350)</td><td>not provided</td><td>from real data</td><td>40 × 40 × 5001</td><td>3D acoustic</td><td>25 fixed sources at the surface</td></tr><tr><td>EFWI (Feng et al., 2023)</td><td>144 000/24 000</td><td>0.35 × 0.35 km<sup>2</sup> (70 × 70)</td><td>612; 3000 m s<sup>-1</sup></td><td>mathematical and from geological faults</td><td>70 × 1000</td><td>2D elastic</td><td>5 fixed sources at the surface</td></tr><tr><td>HEMEWS<sup>S</sup>-3D (this study)</td><td>27 000/3000</td><td>9.6 × 9.6 × 9.6 km<sup>3</sup> (32 × 32 × 32)</td><td>1071; 4500 m s<sup>-1</sup></td><td>random fields added to horizontal layers</td><td>32 × 32 × 800</td><td>3D elastic</td><td>1 source with a random location and a random orientation</td></tr></tbody></table>
not consider varying geological models but focus on regional validated geological models.

## 2.2 3D datasets

Due to the high computational costs of solving 3D partial differential equations (PDEs), only very few 3D datasets are available. CO₂ underground storage has been explored with SciML based on 3D numerical simulations (Grady et al., 2023; Wen et al., 2023; Witte et al., 2023). To support the study of Witte et al. (2023), Annon (2022) provided 4000 simulation results for 3D CO₂ flow through geological models based on the Sleipner dataset complemented by random fields (Equinor, 2020). The Kimberlina dataset also contains 6000 CO₂ leakage rate simulations (Mansoor et al., 2020). However, the geological models in both of those databases are all variants of the geological model that was carefully estimated for a given region, thereby limiting the reproducibility in other areas.

## 2.3 Geophysical datasets

A few datasets of realistic geological units have been developed, such as the Noddyverse dataset of 3D geological models (Jessell et al., 2022). In this dataset, geological models result from the deformation of horizontal layers by successive geological events (folds, faults, unconformities, dykes, plugs, shear zones, and tilts), but no associated ground motion is provided. Along the same line of geological deformation, the OpenFWI database combines geological models with associated waveforms and targets 2D geophysical inversion as the main application (Deng et al., 2022). OpenFWI contains geological models made of horizontal and non-horizontal layers with various folds. It also includes real geological models from field survey areas and models of CO₂ geological storage. To generate the wavefields, the acoustic wave equation is solved in the 2D domains. Waves originate from a line of sources at the surface and wavefields are acquired on a line of receivers at depth. The EFWI database is an extension of OpenFWI to the elastic wave equation, providing two-component ground motion time series (Feng et al., 2023).

Several other studies have computed simulation outputs for the acoustic or elastic wave equation, but their data are not public (e.g. Liu et al., 2021; Ovadia et al., 2023; Zhang et al., 2023). Table 1 summarizes the characteristics of the public datasets and shows that no database provides solutions of the elastic wave equation in 3D domains. Our HEMEW<sup>S</sup>-3D database intends to fill this gap.

# 3 Dataset creation

## 3.1 The elastic wave equation

Elastodynamics describes reversible wave propagation phenomena in solid and fluid domains. In solid mechanics, the solution is represented by a displacement field, $u \in \mathbb{R}^3$, that propagates in a 3D Euclidean space. We consider a truncated propagation domain, $\Omega = [0, L]^3$, with absorbing boundary conditions all around, except the traction-free top surface, and a solution, $u: \Omega \times [0, T] \to \mathbb{R}^3$. The domain length is fixed to $L = 9.6$ km, and the total time is $T = 8$ s. In its most general form, the elastic wave equation is written as

$$ \rho \frac{\partial^2 u}{\partial t^2} = \nabla \lambda (\nabla \cdot u) + \nabla \mu [\nabla u + (\nabla u)^T] + (\lambda + 2\mu) \nabla (\nabla \cdot u) - \mu \nabla \times \nabla \times u + f, \quad (1) $$

where $\rho: \Omega \to \mathbb{R}$ is the material unit mass density; $\lambda: \Omega \to \mathbb{R}$ and $\mu: \Omega \to \mathbb{R}$ are the Lamé parameters, characterizing the thermodynamically reversible mechanical behaviour of the material; and $f$ is the body force distribution. In geomechanics, properties $\rho$, $\lambda$, and $\mu$ are rarely independently characterized due to a lack of measurements. Therefore, it is legitimate to assume that there is a single informative variable from which all parameters can be deduced. In this work, the velocity of shear waves, $v_S$, is the informative variable. Equation (1) can then be rewritten in the following general form:

$$ \mathcal{L}(v_S, u) = f. \quad (2) $$

## 3.2 Earthquake source

In our database, the forcing term $f(x, t) = \text{div } m(x) \cdot s(t)$ is the divergence of a moment tensor density, $m$, localized at a pointwise location. $m$ encodes the source radiation patterns as a double couple representing a pointwise kinematic discontinuity in the media.

The source position is represented by the coordinates $(x_s, y_s, z_s) \in \Omega$, which is not too close from the boundaries to avoid numerical issues due to absorbing boundary conditions. The position is chosen from a product of Latin hypercube sampling (LHS), with

$$ \begin{aligned} x_s &\in [1.2, 8.4]\,\text{km}, \\ y_s &\in [1.2, 8.4]\,\text{km}, \\ z_s &\in [-9.0, -0.6]\,\text{km}. \end{aligned} \quad (3) $$

In addition to the source position, the source is parametrized by the symmetric $3 \times 3$ moment tensor. An equivalent formulation is obtained with the three angles, strike, dip, and rake (Aki and Richards, 1980). With this representation, the angles are samples from a LHS with a strike between 0 and 360°, dip between 0 and 90°, and rake between 0 and 360°.
The source amplitude corresponds to a seismic moment, $M_0 = 2.47 \times 10^{16}$ Nm, and the source time evolution is a spice bench given by $s(t) = 1 - (1 + \frac{t}{\tau}) e^{-\frac{t}{\tau}}$, with $\tau = 0.1$ s (Fig. A1). Due to the linearity of the elastic wave equation (Eq. 1), it is important to notice that the choice of the source time function in the HEMEWS-3D database does not constrain the variability in resulting ground motions.

First, any magnitude can be obtained by applying a scalar factor to the ground motion wavefields. Second, the response to any source time function can be computed from the Green function, $G(x, t)$, which is the fundamental solution of the elastic wave equation, Eq. (1), when the source is an impulse point force located at $x_s$ and occurring at $t = t_0$. For the reference source time function, $s(t)$, the solution, $u(x, t)$, provided in HEMEWS-3D can be written as the convolution of the Green function with the source time function, i.e. $u(x, t) = G(x, t) * s(t)$.

Computing the solution, $u_1$, for a new source time function, $s_1$ (provided that the moment tensor density $m$ and the geological parameter $V_S$ remain the same), is straightforward following the steps below:

1. compute the Fourier transform of the reference source time function, $\hat{s} := \mathcal{F}(s)$, and the solution, $\hat{u} := \mathcal{F}(u)$;

2. derive the Green function in the frequency domain, $\hat{G} = \frac{\hat{u}}{\hat{s}}$;

3. compute the Fourier transform of the new source time function, $\hat{s}_1$;

4. compute the new solution in the frequency domain, $\hat{u}_1 = \hat{G} * \hat{s}_1$;

5. deduce the new solution in the temporal domain, $u_1 = \mathcal{F}^{-1}(\hat{u}_1)$.

From these remarks, one should remember that ground motion wavefields in the HEMEWS-3D database originate from pointwise sources with a different location, $x_s \in \mathbb{R}^3$, and orientation, $\theta_s \in \mathbb{R}^3$, but the same source time function.

## 3.3 Heterogeneous geological models

The HEMEWS-3D database contains samples $\{V_S^{(i)}, x_s^{(i)}, \theta_s^{(i)}, \hat{u}^{(i)}\}_i$ that satisfy Eq. (2) ($\hat{u}$ denotes the velocity field obtained as the time derivative of the displacement field $u$). The 3D geological model set, $V_S(x)$, is a set of non-stationary random fields defined as a mean stair function (horizontal homogeneous layers) to which fluctuations are added, as illustrated in Fig. 1.

### 3.3.1 Homogeneous models

A homogeneous layer that is 1.8 km thick is imposed at the bottom of each geological model, with a $V_S$ value of $V_{S,min} = 4500$ m s$^{-1}$. The minimum S-wave velocity is

$V_{S,min} = 1071$ m s$^{-1}$. Above the bottom layer, the number of horizontal layers and their thickness are randomly chosen for each sample $V_S^{(i)}$, with the sole constraint to fill the total depth with two to seven layers. In particular, this means that velocity values are not sorted by depth. This choice is discussed in depth in Sect. 6. Then, the mean layerwise value is drawn from the uniform distribution, $\mathcal{U}([\mu_1, \mu_2])$. All layerwise values are chosen independently. Values of $\mu_1 = V_{S,min}/0.6 = 1785$ m s$^{-1}$ and $\mu_2 = V_{S,max}/1.4 = 3214$ m s$^{-1}$ were determined to ensure that most values remain bounded within $[V_{S,min}, V_{S,max}]$ after the addition of random fields in each layer (see Table 2 for a summary of the parameters).

To recover the other geological properties, the ratio of P- to S-wave velocity was fixed to $V_P/V_S = 1.7$. The density, $\rho$, is computed as a function of the P-wave velocity (Molinari and Morelli, 2011):

$$ \rho = 1.6612V_P - 0.4721V_P^2 + 0.0671V_P^3 - 0.0043V_P^4 + 0.000106V_P^5. \quad (4) $$

Attenuation factors for P waves ($Q_P$) and S waves ($Q_S$) are computed as follows:

$$ Q_P = \max \left( \frac{V_P}{20}, \frac{V_S}{5} \right), \quad Q_S = \frac{V_S}{10}. \quad (5) $$

### 3.3.2 Addition of heterogeneities

The layers' thicknesses and mean values describe the general structure of the propagation domain and they correspond to the prior physical information usually available. However, geomaterials of the Earth's crust contain a lot of variability, especially along the horizontal directions. This heterogeneity can be represented by random fields, characterized by their correlation length and coefficient of variation. Following previous studies on geological heterogeneity (e.g. Hartzell et al., 2010; Imperatori and Mai, 2013; Khazaie et al., 2016; Scalise et al., 2021; Thompson et al., 2007), we drew random fields with a von Kármán correlation kernel and a Hurst exponent of 0.1 (Chernov, 1960) (marginal distributions are log-normal to preserve positive values).

In order to provide a sufficient dataset variability, the choice of correlation lengths and coefficients of variation is tricky yet crucial (Colvez, 2021). The correlation length gives an idea of the distance above which two points, $x_A$ and $x_B$, have independent geological properties, $V_S(x_A)$ and $V_S(x_B)$. We chose correlation lengths randomly from the set $\{1.5, 3, 4.5, 6\}$ km to mix samples with small- and large-scale heterogeneity. In addition, large coefficients of variation were chosen to provide high geological contrasts, following the folded normal distribution, $|\mathcal{N}(0.2, 0.1)|$, with a mean of 0.2 and coefficient of variation of 0.1. Coefficients of variation of around 20% are common at the surface (Arroucau, 2020), while it is known that values of up to 40% can be found locally (El Haber et al., 2021).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Geological models are built by adding heterogeneities to randomly chosen horizontal layers. Elastic waves are then propagated from a source with a random position and random orientation to the surface, where velocity wavefields are synthesized.</figcaption></figure>

<table><caption>**Table 2.** Statistical distribution of each parameter describing the geological models. Mean $V_S$ values, coefficients of variation, and correlation lengths are chosen independently in each layer. Since the bottom layer has a constant thickness of 1.8 km, it is not included in these parameters.</caption><thead><tr><th>Parameter</th><th>Statistical distribution</th></tr></thead><tbody><tr><td>Number of heterogeneous layers ($N_ℓ$)</td><td>$U(\{1, 2, 3, 4, 5, 6\})$</td></tr><tr><td>Layer thickness ($h_1, \dots, h_{N_ℓ}$)</td><td>$U(\{(h_1, \dots, h_{N_ℓ}) > 0 | h_1 + \dots + h_{N_ℓ} = 7.8\})$</td></tr><tr><td>Mean $V_S$ value per layer</td><td>$U([1785, 3214])$</td></tr><tr><td>Layerwise coefficient of variation</td><td>$|N(0.2, 0.1)|$</td></tr><tr><td>Layerwise correlation length along x</td><td>$U(\{1.5, 3, 4.5, 6 \text{ km}\})$</td></tr><tr><td>Layerwise correlation length along y</td><td>$U(\{1.5, 3, 4.5, 6 \text{ km}\})$</td></tr><tr><td>Layerwise correlation length along z</td><td>$U(\{1.5, 3, 4.5, 6 \text{ km}\})$</td></tr></tbody></table>

The 3D random field computation is made highly efficient by the use of the spectral representation (Shinozuka and Deodatis, 1991; de Carvalho Paludo et al., 2019). With this formulation, a centred Gaussian random field, $V_S$, determined by its autocovariance function, $\mathcal{R}$, can be decomposed as a sum of independent identically distributed random variables $(V_{S,n})_{-N \le n \le N}$, with a uniform distribution over $[0, 2\pi]$.

$$V_S(x) = \sum_{n=-N}^{N} \sqrt{2\hat{\mathcal{R}}(n\Delta k)} \cos(n\Delta k \cdot x + V_{S,n}),$$

where $\hat{\mathcal{R}}$ is the Fourier transform of the autocovariance function $\mathcal{R}$ and $\Delta k$ is the unit volume in $\mathbb{R}^3$.

Finally, $V_S$ values are clipped between $V_{S,min} = 1071 \text{ m s}^{-1}$ and $V_{S,max} = 4500 \text{ m s}^{-1}$. These bounds correspond to the velocity of shear waves in hard sediments and at the bottom of the continental crust (Molinari and Morelli, 2011).

It should be noted that all layers have distinct coefficients of variation and correlation lengths, meaning that different random fields are drawn inside each layer. Also, random fields are drawn only once for each set of parameters.

### 3.3.3 Representation in the database

Geological realization, $V_S^{(i)}$, is discretized over a grid of 32 × 32 × 32 elements (corresponding to x, y, and z axes). The geological models are provided as `.npy` arrays. The total size of the geological dataset is 3.9 GB, split in 15 files of 2000 geological models for easier data management. Additionally, metadata give parameters of each layer: the mean $V_S$ value; the thickness; the coefficient of variation; and the correlation lengths along x, y, and z.

## 3.4 Solutions of the wave equation

The elastic wave equation was solved in each domain $i$ by means of the open-source code SEM3D (https://github.com/sem3d/SEM, CEA et al., 2017) (Touhami et al., 2022) based on the spectral element method (Faccioli et al., 1997; Komatitsch and Tromp, 1999). The dimension of the simulation mesh is prescribed by the maximum frequency, $f_{max}$, one aims at exactly resolving. In this study, $f_{max}$ was fixed at 5 Hz, which is relatively high for this type of simulation. Many simulations have been conducted so far with an
accuracy of up to 1 or 2 Hz (Rekoske et al., 2023; Rosti et al., 2023), while high-fidelity simulations for local realistic earthquake scenarios extend to up to 10 Hz (Castro-Cruz et al., 2021; De Martin et al., 2021; Heinecke et al., 2014) (and exceptionally to up to 18 Hz, such as in Fu et al., 2017). Then, the smallest wavelength, $\lambda_{min} = V_{S,min}/f_{max}$, must be described on the mesh by at least five quadrature points. With seven Gauss-Lobatto-Legendre quadrature points per mesh element, this leads to elements of size $h = \frac{7}{5} \cdot \frac{V_{S,min}}{f_{max}} = 300$ m. This explains how 32 elements in each direction amount to a domain size of $L = 9.6$ km. The time-marching scheme is a leap-frog second-order accurate explicit scheme solved for velocity fields.

To maintain reasonable computational loads and reflect realistic scenarios, velocity fields were recorded only at the surface of the propagation domain. A regular grid of 32 × 32 sensors was placed between 150 and 9450 m in both horizontal directions (implying a distance of 300 m between two neighbouring sensors). At each monitoring point, the three-component velocity field is synthesized with a 100 Hz sampling frequency for between 0 and 8 s. Although the sampling frequency is higher than the Nyquist frequency (i.e. $2 \times f_{max} = 10$ Hz), the value of 100 Hz was chosen to match the temporal resolution of recorded time series in several publicly accessible datasets (e.g. STEAD, Mousavi et al., 2019; INSTANCE, Michelini et al., 2021). The sampling frequency is sufficient to allow for an accurate computation of peak ground velocity (PGV), derive the acceleration time series with finite differences, and compute the peak ground acceleration (PGA). Figure 2 illustrates velocity waveforms at eight virtual sensors.

### 3.4.1 Representation in the database

Velocity fields are provided as .h5 files. Each file contains three keys: uE, uN, and uZ that correspond to the three components of ground motion (east-west, north-south, and vertical). Each velocity field is of shape 32 × 32 × 800, where the first index corresponds to the y axis, the second index to the x axis, and the third index to the temporal axis.

Files are gathered in .zip archives containing 100 simulation results. The 300 .zip files amount to 263.4 GB. They are downloadable individually (0.87 GB per file).

## 4 Dataset analysis

### 4.1 Descriptive statistics of the temporal evolution

Since most of the geological parameters are chosen uniformly randomly (Table 2), the geological dataset is well balanced: geological models with one to six layers are equipartitioned, and all random field parameters have approximately the same frequency. Mean $V_S$ values range from 1756 to 3145 m s⁻¹.

The first-wave arrival time is a crucial parameter for earthquake early warning, and seismic phase picking is a common task with deep learning models. Arrival time depends on the distance between the earthquake source and the monitoring sensor as well as the geological properties on the propagation path. Wave arrival times are usually determined from recordings, either manually by experts or using machine learning methods. However, it is possible to compute almost exact arrival times from synthetic velocity fields since ground motion is almost zero before the first wave arrival. Therefore, we obtained wave arrival times for P-waves as the earliest time step where the amplitude exceeds 0.1 % of the maximum amplitude. Due to the source depth variability and the different wave velocities in the geological models, first wave arrival times significantly vary among samples and among sensors. Figure 3a shows that 10 % of velocity time series are initiated before 0.65 s, while 10 % of time series are still null after 2.17 s.

As expected, the P-wave arrival time is strongly correlated with the hypocentral distance (Fig. 4a) since shorter hypocentral distances are associated with shorter propagation paths. The mean velocity on the propagation path also influences the first wave arrival time, but variability is higher. Figure 4b indeed shows that the P-wave arrival time is negatively correlated with the mean S-wave velocity in the whole domain. It confirms that waves propagate more slowly when the mean velocity is lower. The mean velocity gives an approximation of the velocity values encountered by the waves along the propagation path. In particular, the mean velocity does not depend on the sensor position with this approximation.

The temporal evolution of ground motion can also be characterized by its relative significant duration (RSD). It corresponds to the duration of the signal between 5 % and 95 % of the Arias intensity ($I_A$) (Arias, 1970).

$$I_A = \frac{\pi}{2g} \int_{0}^{T} a^2(t)dt, \qquad (6)$$

where $a(t)$ is the acceleration field and $T$ is the total duration of the signal. Figure 3b shows that RSD covers a large variation range, from 0.17 to 7.60 s. This variability is illustrated in Fig. 2, where the grey areas represent the RSD. One can especially notice that samples with a strong pulse have a short RSD. Indeed, most of the energy is concentrated around the pulse.

Quantitatively, the HEMEWS-3D database contains short ground motions since 10 % of time series have a RSD lower than 0.24 s as well as longer ground motions, where 10 % of time series have a RSD longer than 3.42 s. The median RSD is 1.06 s. These low RSD values are related to the absence of high-frequency components in the coda and the dominance of high pulse-like time series in cases with shallow sources and low heterogeneity contrasts. Combined with the short P-
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Three-component velocity waveforms synthesized at eight virtual sensors on a line parallel to the y axis at x = 8.85 km. The shaded area extends from 5% to 95% of the Arias intensity; hence, its length equals the relative significant duration (RSD). The corresponding geological model is shown at the top. The source is located at (2.04, 3.64, -2.17) km. The associated frequency spectra are given in Fig. A3.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Distributions of the temporal features of velocity time series at each monitoring sensor and for 30 000 samples. (a) The first P-wave arrival time is computed on the vertical component. (b) The relative significant duration (RSD) is shown for the east-west component, and the results are very similar for the two other components.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4. For each sample and each sensor, the P-wave arrival time is shown against (a) the hypocentral distance and (b) the mean S-wave velocity in the 3D domain.</figcaption></figure>

wave arrival times, RSD values also justify the fact that the 8 s window contains a significant part of ground motion.

## 4.2 Descriptive statistics related to energy content

The PGV is computed as the maximum absolute value over all time steps separately for each component. The PGV is slightly lower on the vertical component than the two horizontal components (Fig. 5). It is very similar between the east-west and north-south components, which is expected since the HEMEWS-3D database is statistically invariant per horizontal rotation. Figure 5 shows that the PGV extends over 3 orders of magnitude, with the first percentile being equal to 0.89 cm s⁻¹, while the 99th percentile is equal to 129.3 cm s⁻¹. The median PGV is 8.9 cm s⁻¹. The 1st-to-99th-percentile interval is in line with ground motion observed within a few kilometres of moderate-magnitude earthquakes (e.g. Convertito et al., 2022).

When the propagation path is longer, seismic waves encounter more geological heterogeneities. They create a dispersion and diffraction of waves that spread the energy signal over time. Larger hypocentral distances are associated with longer propagation paths. Figure 6a then shows that the PGV is negatively correlated with the hypocentral distance.

It is also known that the seismic energy, $E_s$, generated by a fault rupture is

$$E_s = \frac{M_0 \Delta \sigma}{2\mu}, \qquad (7)$$

where $M_0$ is the seismic moment, $\Delta\sigma$ is the stress drop, and $\mu$ is the shear modulus at the fault location. Knowing that the shear wave velocity is written as $V_S = \sqrt{\mu/\rho}$, Eq. (7) indicates that the seismic energy is inversely proportional to $V_S^2$. Furthermore, Fig. 6b confirms that the PGV is negatively correlated with the velocity of S waves at the source location.

## 4.3 Distribution of pseudo-spectral acceleration (PSA)

The pseudo-spectral acceleration (PSA) is a commonly used metric to estimate structural response. It evaluates the maximal acceleration of a 1-degree-of-freedom oscillator (with a 5% damping), with a natural period, $T$. At $T = 0.2$ s, the PSA value in the HEMEWS-3D database is between $2.3 \times 10^{-3}$ and 81.2 g. It decreases to the interval between $5.8 \times 10^{-4}$ g and 6.1 g at $T = 1.0$ s. Figure 7 additionally shows that there exists a negative correlation between the PSA and the hypocentral distance. The distance-dependent PSA values can be compared with existing ground motion models (GMMs), as shown in Fig. 8.

GMMs provide an analytical formula to compute intensity measures, such as PSA and PGV, based on regression analyses. They are mainly derived from databases of recorded earthquakes, although numerical simulations can also be used. The PSA estimated from the HEMEWS-3D database is compared with four GMMs (all taken with a moment magnitude of 4.9 that corresponds to the seismic moment of the HEMEWS-3D database):

1. GMM from Atkinson and Boore (2006), computed with $V_{S,30} = 2000 \text{ m s}^{-1}$;

2. GMM from Atkinson (2015), computed with a depth of 1 km (solid line in Fig. 8) and a depth of 7.5 km (dashed line in Fig. 8);

3. GMM from Chiou and Youngs (2014), computed with $V_{S,30} = 2000 \text{ m s}^{-1}$, a depth of 1 km, and different dip and rake angles (these last two factors having little influence on the PSA);

4. GMM from Shahjouei and Pezeshk (2016).

Figure 8 shows that the mean PSA computed from the HEMEWS-3D database is in good agreement with all GMMs and that the standard deviation is smaller than GMMs. The Atkinson (2015) GMM illustrates the influence of
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The peak ground velocity (PGV) is computed as the maximum absolute value over all time steps separately on each component. There is one value for each of the 32 × 32 sensors and each of the 30 000 samples.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** For each sample and each sensor, the PGV is shown against (a) the hypocentral distance and (b) the S-wave velocity at the source location. The PGV is computed on the east-west component, results are very similar for the two other components.</figcaption></figure>

depth, especially for the smallest hypocentral distances. The HEMEWS$^S$-3D PSA values fit better with the GMM PSA for shallow sources (solid purple line in Fig. 8) than deeper sources (dashed purple line).

## 4.4 Dimensionality

In supervised deep learning, it is always challenging to determine whether the size of the database (i.e. the number of samples) is sufficient to represent its variability. This question relates to the definition of the intrinsic dimension of the dataset, which indicates the number of hidden variables that should be necessary to represent the main features of the samples. In the following sections, we provide insights into this question, with the intrinsic dimension based on the principal component analysis (Sect. 4.4.1), the correlation dimension (Sect. 4.4.2), the maximum likelihood estimate (Sect. 4.4.3), and the structural similarity index (Sect. 4.5).

### 4.4.1 Principal component analysis (PCA)

The principal component analysis (PCA) decomposes data in principal components that correspond to the directions where data vary the most. For different sizes of datasets, we compute the number of principal components required to retain 95% of variance and define this number as the intrinsic dimension of data. The 3D geological models and the 3D ground motion wavefields are transformed into 1D vectors to perform the PCA. To reduce the memory requirements, ground motions are analysed only on the east-west component. Geological models are represented by 32 × 32 × 32 = 32 768 points, and ground motions contain 16 × 16 × 320 = 81 920 points (16 sensors in directions x and y and 320 time steps between 0 and 6.4 s). To ease the computation on the large sample covariance matrix, an incremental PCA algorithm was used (Ross et al., 2008).

Table 3 and Fig. B1 show that more than 1000 principal components are needed to reconstruct the geological models with high accuracy, whereas the intrinsic dimension of ground motion wavefields is around 4900. It is a reasonable
<figure><img src="image_9.png" /><figcaption>Figure 7. For each sample and each sensor, the pseudo-spectral acceleration is shown against the hypocentral distance at period T = 0.2 s (a) and T = 1.0 s (b). The PSA is computed on the east-west component, and the results are very similar for the two other components.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. Horizontal PSA at period T = 0.2 s as a function of hypocentral distance (a) and epicentral distance (b) for GMMs by Atkinson (2015), with a depth of 1 km (solid purple line, a); Atkinson (2015), with a depth of 7.5 km (dashed purple line, a); Atkinson and Boore (2006) (orange, b); Chiou and Youngs (2014) (blue, b); Shahjouei and Pezeshk (2016) (green, b); and our HEMEW$^{S}$-3D database (black). Solid lines correspond to the mean PSA and shaded areas to 1 standard deviation. Horizontal PSA is computed as the geometrical mean of east-west and north-south components.</figcaption></figure>

fact that the wavefield intrinsic dimension is larger than the geological dimension since wavefield variability is created by geological variations and the source position. However, due to its linearity, the PCA requires a large number of components to accurately represent complex patterns. Therefore, it may overestimate the intrinsic data dimension.

### 4.4.2 Correlation dimension

An alternative dimensionality measure was introduced by Grassberger and Procaccia (1983) as the correlation dimension, which characterizes the distance between pairs of samples. For a dataset of N samples {$V_S^{(i)}$}$_{1 \le i \le N}$ and a given radius r, the correlation dimension ($C_N(r)$) is defined as the ratio of sample pairs ($V_S^{(i)}, V_S^{(j)}$)$_{i \ne j}$ that are at a distance lower than r:

$$C_N(r) = \frac{2}{N(N-1)} \sum_{i=1}^{N} \sum_{j=i+1}^{N} 1 (\|V_S^{(i)} - V_S^{(j)}\| \le r). \quad (8)$$

Figure B2 indicates a correlation dimension of 8 for the geological dataset, which is significantly lower than the PCA dimension. In fact, it is known that the correlation dimension may underestimate the intrinsic dimension, especially “when data are scattered” (Qiu et al., 2023), which is likely to be the case in high-dimensional spaces. However, the correlation dimension of the ground motion wavefields is debatable as it drops to 2. Figure B3c shows that the log-log representation of $C_N(r)$ does not produce an obvious linear part, which makes it difficult to identify the correlation dimension.
<table><caption>**Table 3.** Database intrinsic dimension estimated by PCA, correlation dimension, and maximum likelihood estimator (MLE) for the geological database and the velocity fields database, depending on the number of data samples.</caption><thead><tr><th rowspan="2">No. of samples (×10<sup>3</sup>)</th><th colspan="5">Geological database</th><th colspan="5">Velocity field database</th></tr><tr><th>2</th><th>6</th><th>10</th><th>20</th><th>30</th><th>2</th><th>6</th><th>10</th><th>20</th><th>30</th></tr></thead><tbody><tr><td>PCA</td><td>491</td><td>766</td><td>880</td><td>1006</td><td>1094</td><td>853</td><td>2057</td><td>2853</td><td>4073</td><td>4875</td></tr><tr><td>Correlation dimension</td><td>8.8</td><td>8.3</td><td>8.4</td><td>8.3</td><td>8.2</td><td>2.2</td><td>2.3</td><td>2.3</td><td>2.3</td><td>2.3</td></tr><tr><td>MLE</td><td>17.9</td><td>23.4</td><td>26.7</td><td>31.5</td><td>33.9</td><td>107.2</td><td>129.0</td><td>116.6</td><td>113.8</td><td>110.9</td></tr></tbody></table>

### 4.4.3 Maximum likelihood estimator (MLE) intrinsic dimension

Levina and Bickel (2004) proposed another approach based on the maximum likelihood estimator (MLE) of the distance to the closest neighbours. Figure B4 shows the evolution of the intrinsic dimension as a function of the number of samples for geological models and velocity wavefields. The intrinsic dimension of geological models is 34, while the dimension of velocity wavefields is higher (around 110). Although this method may still underestimate data with high intrinsic dimensionality (Qiu et al., 2023), it provides higher estimates than the correlation dimension.

It can also be noted that the intrinsic dimension increases with the number of samples, as was observed for the PCA. This may reflect a flaw in the intrinsic dimension's definition, or it may indicate that despite being already large, our database of 30 000 samples does not capture all the variables.

### 4.5 Structural similarity

The correlation dimension is computed from the Euclidean distance between pairs of geological models. However, pointwise metrics do not necessarily best represent similarities between geological models, and alternative metrics such as the structural similarity index measure (SSIM) have been introduced for this purpose (Wang et al., 2004). This index theoretically ranges from 0 to 1, with 0 indicating no similarity and 1 indicating perfectly similar geological models (although values between -1 and 0 can be obtained numerically from the covariance computation). The SSIM of two geological models A and B is defined as

$$SSIM(A, B) = \frac{(2\mu_A\mu_B + C_1)(2\sigma_{AB} + C_2)}{(\mu_A^2 + \mu_B^2 + C_1)(\sigma_A^2 + \sigma_B^2 + C_2)}, \quad (9)$$

where $\mu_A$ and $\mu_B$ are the means of A and B, $\sigma_A$ and $\sigma_B$ are the unbiased estimators of the variance of A and B, $\sigma_{AB}$ is the unbiased estimator of the covariance of A and B, and $C_1$ and $C_2$ are constants determined from the range of A and B values.

Figure B5a and b illustrate two pairs of geological models with the same SSIM (of 0.6), meaning there is a rather high similarity. The first geologies have similar mean values but different heterogeneities, resulting in a low Euclidean

<figure><img src="image_11.png" /><figcaption>**Figure 9.** The structural similarity index measure (SSIM) quantifies the visual resemblance between images, in a way that should mimic human perception. For each SSIM value *s* on the *x* axis, the percentage of geological pairs being more similar than *s* is reported on the *y* axis.</figcaption></figure>

distance (Fig. B5a), while the second geologies have different mean values, leading to a higher Euclidean distance (Fig. B5b).

To give insights into the sparsity of the geological database, Fig. 9 shows that only 1.4% of geological pairs have a SSIM greater than 0.2. This means that geological models are generally very distinct from each other in the HEMEWS-3D database. For comparison, the 2D OpenFWI dataset leads to significantly higher SSIM, with 31% of geologies having a SSIM larger than 0.2 (3000 models were chosen from each of the 10 OpenFWI families; Deng et al., 2022).

## 5 Applications

### 5.1 Dimensionality reduction in geological models

Dimensionality analyses have shown that at least 1000 principal components are necessary to represent geological models with enough accuracy, as measured by the reconstructed variance. This means that the PCA provides a basis for 3D models to decompose a wide diversity of geological models. One can consider geological models that are very differ-
ent from the random fields contained in the HEMEWS-3D database – for instance, embedding a basin shape.

Figure 10 shows that 1000 principal components allow for a good reconstruction of the basin shape with correct velocity values inside and outside the basin. Edges are slightly blurred, which is expected since sharp contrasts correspond to high spatial frequencies that require many principal components. This example illustrates the generalization ability of the HEMEWS-3D database from a geometrical point of view. To match the design of the HEMEWS-3D database, the velocity values are chosen within the same bounds. If one were to consider real sedimentary basins, they should apply rescaling to target lower velocity values.

The influence of the PCA reconstruction on the generated velocity wavefields was investigated in more details in Lehmann et al. (2022). It was shown that wavefields created by the propagation of seismic waves inside the reconstructed geological models and the reference model are very similar. When the initial geological model has strong heterogeneities, heterogeneities tend to be blurred in the PCA reconstruction, which reduces the dispersion of seismic waves. As a consequence, velocity wavefields generated inside the reconstructed geological model have slightly larger amplitudes.

## 5.2 Velocity fields predictions

Since the HEMEWS-3D database associates geological models and sources with their corresponding velocity wavefields, it can serve to predict the latter using the former. Neural operators are one class of SciML models that have shown great success in the prediction of parametric PDEs. One can mention in particular the Multiple-Input Fourier Neural Operator (MIFNO; Lehmann et al., 2024) that uses the fast Fourier transform to learn the frequential representation of the elastic wave equation and a dedicated handling of the source term (Fig. C1).

For each geological model and source in the HEMEWS-3D database, the MIFNO predicts the velocity field at each surface point. Figure 11 illustrates that the MIFNO gives accurate predictions for samples with different ground motions. This shows that the variability and size of the HEMEWS-3D database are appropriate to train complex SciML models.

## 5.3 Other potential applications

Thanks to the large number of simulations, one can also envision studying the variability in ground motion to capture its statistical distribution. In particular, one can investigate the best sampling that minimizes the number of samples while preserving the largest ground motion variability (Tarbali and Bradley, 2015).

The surface ground motion can also be considered an “outcropping bedrock” response, which is classically used in 1D site-effect and soil-structure interaction analyses and which may require deconvolution.

# 6 Limitations and perspectives

Since HEMEWS-3D is the first database that provides 3D ground motion, it is constrained by some hypotheses to control its size and allow machine learning applications.

First, the minimum S-wave velocity of 1071 m s⁻¹ is rather high when compared to S-wave velocities in soft sediments (typically of a few hundred metres per second) but coherent for hard sediments (Molinari and Morelli, 2011). V_s values in the HEMEWS-3D database are also in line with national velocity models with a low spatial resolution that display surface values of around 2000 m s⁻¹. The choice of minimum V_s implies in particular that the HEMEWS-3D dataset is not targeted towards site-effect analyses in sedimentary basins. One should also note that the vertical resolution of the geological models is 300 m, while very low V_s values are more commonly encountered in the first tens of metres. These low values would be averaged with higher deeper values in our models. In particular, this means that the upper velocity in the HEMEWS-3D database should be understood as the average velocity between 0 and -300 m. It is not comparable with the common notion of V_{s,30} (average V_s in the first 30 m). Reducing the minimum velocity poses no theoretical limitation but would increase the computational cost of the subsequent numerical simulations since it increases the number of mesh elements.

Second, the maximum S-wave velocity of 4500 m s⁻¹ corresponds to existing V_s values at the bottom of the Earth's crust that are often adopted in velocity models (Duverger et al., 2021; Molinari and Morelli, 2011). In addition, the bottom layer has a fixed thickness and value that originates from earlier works. Therefore, variability is considered only above this constant layer.

Third, we do not constrain the ordering of layerwise V_s values to provide a large database variability that is essential for machine learning perspectives. Similar choices were made in the widely used OpenFWI database (Deng et al., 2022), in which the FlatVel-B family is based on a random arrangement of layers. This means that some layer arrangements may be unphysical – for instance, if the mean values are linearly decreasing with depth. However, it is important to notice that the physics of wave propagation is still satisfied in those situations, which is the main concern of this work. One choice opposite to ours would be to impose layerwise values increasing with depth (as done in the OpenFWI FlatVel-A family for instance). However, this would remove all models with velocity inversion, which can be found in complex geological contexts. From the metadata provided, users can filter geological models with custom criteria to exclude those outliers from their study. For instance, one custom criterion could be the impedance contrast, defined as the ratio between impedance (Z = V_s × ρ) in one geological layer, Z_ℓ, and the layer above, Z_{ℓ-1}. Figure A2 illustrates the distribution of impedance contrasts and shows that they are
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The original geological model (a) contains two homogeneous layers and a circular basin inserted inside the top layer. The PCA reconstruction was obtained with 1000 PCA components (b).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** For two spatial points, the velocity field predicted by the MIFNO (dashed red line) is compared with the reference from the HEMEWS$^S$-3D database (black line).</figcaption></figure>

realistic. A total of 9028 geological models have a minimum impedance contrast smaller than 0.7.

Additionally, more diverse configurations could be designed by relaxing the assumption that all geological parameters depend on a single variable. This would imply, for instance, varying the $V_p/V_s$ ratio. Random anisotropic heterogeneities can also be generated for more diversity (Ta et al., 2010).

The domain size was limited to 9.6 km to prove that SciML was possible with a manageable dataset size. This size allows for reasonable local studies and is already larger than existing 2D databases (Table 1). Extending the spatial size is certainly of interest for some seismological applications and requires additional computational costs. In summary, with larger computational budgets and lesser memory constraints, it would become possible to

- consider larger and deeper geological models,

- design models with a higher spatial resolution,

- include lower minimum $V_s$ values,

- increase the frequency limit of the wave propagation simulations,

- increase the spatial sampling of virtual sensors,

- increase the temporal duration of signals to match the longer epicentral distances coming from larger models.

It should also be noted that numerical simulations are only valid for a frequency of up to 5 Hz due to the mesh design, with numerical pollution for frequencies larger than 5 Hz. We observed that it is crucial to apply a low-pass filter (with a cutoff frequency of 5 Hz) to the velocity fields before using machine learning models; otherwise, the model may try to fit numerical noise.

# 7 Code and data availability

The database is referred to as Lehmann (2023) and can be downloaded at https://doi.org/10.57745/LAI6YU. The wave propagation code, SEM3D, is available at https://github.com/sem3d/SEM (CEA et al., 2017). The code used to gener-
ate the HEMEWS³D database is given at https://github.com/lehmannfa/HEMEW3D (Lehmann, 2024).

## 8 Conclusions

We presented the HEMEWS³D (HEterogeneous Materials and Elastic Waves with Source variability in 3D) database that contains 30 000 geological models, source parameters, and time- and space-dependent surface wavefields generated by the propagation of seismic waves through each geological model. This database was conceived as a solution to the forward problem of wave propagation.

Geological models are built from horizontal layers that are randomly arranged, and they correspond to the velocity of shear waves ($V_S$). They represent a domain of size $9.6 \times 9.6 \times 9.6$ km discretized in 300 m wide elements. $V_S$ values are comprised between 1071 and 4500 m s⁻¹. Random fields are then added independently in each geological layer to create 3D heterogeneities. Their parameters (coefficients of variation and correlation lengths) vary widely to cover diverse geological configurations and are given as metadata. Geological models are provided as cubes with $32 \times 32 \times 32$ voxels.

Seismic waves propagate numerically from the earthquake source to the surface. Pointwise sources have a random position and orientation. Ground motion wavefields are synthesized at the surface of the propagation domain by a grid of $32 \times 32$ sensors recording for 8 s. Simulations are conducted with the high-performance computing code SEM3D and amount to a total computational time of $9 \times 10^5$ core hours. The dataset description shows that the 8 s time window covers most significant ground motion at the surface.

Ground motion characteristics differ strongly between samples. They were analysed in terms of relative significant duration (RSD), P-wave arrival time, peak ground velocity (PGV), and pseudo-spectral acceleration (PSA). In addition to quantifying the distributions of essential intensity measures in seismology, these analyses confirm expected relationships between physical parameters and ground motion characteristics. In particular, hypocentral distance, $V_S$ at the source location, and mean velocity were investigated. Comparisons with ground motion models (GMMs) show that PSA is comparable with estimates from recorded earthquakes. This indicates the usefulness of the HEMEWS³D database in complementing databases of recorded earthquakes.

Due to the size of individual samples in the HEMEWS³D database, one may wonder whether data could be represented with fewer parameters to reduce memory requirements. To this end, we explored different methods to estimate the data-intrinsic dimension, and we exemplified the well-known fact that they can lead to very different values. Taking the MLE as a lower bound, one can argue that the intrinsic dimension of the geological database is at least 30. In addition, the low values of the SSIM indicate that geologies are sparse and distant from each other in the HEMEWS³D database.

Concerning the velocity wavefields, the PCA and the MLE confirm the intuition that the intrinsic dimension is larger than the geological dimension since the source adds variability to the time arrival of wavefields as well as their location at the surface. In this situation, it is reasonable to consider that the intrinsic dimension of ground motion is at least on the order of 100. However, if data are decomposed with the PCA, then the number of principal components is a few thousand. The correlation dimension yields questionable estimates of the intrinsic dimension that contradict our intuition and the PCA and MLE outcomes.

By providing a large number of physics-based simulations, the HEMEWS³D database offers new perspectives from which the relationship between geological properties and surface ground motion can be studied. It led to the first neural operator that predicts 3D ground motion, but many applications in statistics, scientific machine learning, and deep learning are envisioned. We designed the database to be as generic as possible, and we believe that several scientific communities can benefit from it.
Appendix A: Dataset description

<figure><img src="image_14.png" /><figcaption>Figure A1. The seismic moment function in the HEMEWS-3D database.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure A2.** Distribution of impedance contrasts in the HEMEWS-3D database. The impedance contrast is computed as the ratio between impedance in one layer ($Z_ℓ = Vs_ℓ × ρ_ℓ$) and impedance in the layer above it ($Z_{ℓ-1}$).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure A3.** Frequency spectra corresponding to Fig. 2. Each row corresponds to a different sensor, located on a line parallel to the y axis at x = 8.85 km. Each column corresponds to one of the three components. The grey line denotes frequencies larger than 5 Hz where numerical simulations are not accurate.</figcaption></figure>
# Appendix B: Dimensionality of data

## B1 Principal component analysis

The intrinsic dimension based on the principal component analysis (PCA) components has been evaluated with the *scikit-dimension* package. Figure B1 illustrates the number of PCA components required to retain 95% of the data variance depending on the number of samples. It can be noted that, for computational reasons, the velocity fields are represented by a single component (the east-west component, parallel to the *x* axis) for all three methods.

For comparison purposes, the wavefield intrinsic dimension is also computed for a previous version of the database, where the source has a fixed position and orientation (HEMEW-3D database, https://entrepot.recherche.data.gouv.fr/dataset.xhtml?persistentId=doi:10.57745/LAI6YU&version=1.0 (last access: 12 July 2024)). With this database, the wavefield intrinsic dimension was around 3200. It is a reasonable fact that adding degrees of freedom with a random source increases the variability in data.

## B2 Correlation dimension

The correlation dimension is determined as the slope of the linear part in the log-log representation of $C_N$ (Fig. B3). This definition is subject to some interpretation since one should determine which portion constitutes the linear part. Nevertheless, we found that small variations in the linear part limits had very little influence on the slope estimate (less than one unit).

## B3 MLE-based intrinsic dimension

The intrinsic dimension based on the maximum likelihood estimator (MLE) has been computed with the *scikit-dimension* package. Figure B4 shows the evolution of the intrinsic dimension as a function of the number of samples for geological models and velocity fields.

<figure><img src="image_17.png" /><figcaption>**Figure B1.** Number of principal components (y axis) required to represent 95% of the variance in data as a function of the dataset size (x axis) for geological models (a) and ground motion wavefields (b). For ground motion, the HEMEW-3D database is used for the fixed source (black line) and the HEMEW<sup>S</sup>-3D database corresponds to the blue line.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

**Figure B2.** Correlation dimension (y axis) as a function of the dataset size (x axis) for geological models (a) and ground motion wavefields (b). For ground motion, the HEMEW-3D database is used for the fixed source (black line) and the HEMEW$^S$-3D database corresponds to the blue line.

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

**Figure B3.** The correlation dimension, $C_N(r)$, is computed from the number of samples that are at (Euclidean) distance lower than $r$ for different values of $r$ (Eq. 8). The correlation dimension is then obtained as the slope of the linear part in the log-log representation. (a) Correlation dimension for 30 000 geological models. (b) Correlation dimension for 30 000 ground motion wavefields with a fixed source (HEMEW-3D database). (c) Correlation dimension for 30 000 ground motion wavefields with a random source (HEMEW$^S$-3D database).
## B4 Structural similarity index

Figure B5 exemplifies two pairs of geological models with high similarity (SSIM of 0.6) but different properties. The first pair (Fig. B5a and b) has similar mean values but different heterogeneities, while in the second pair, geological models are almost homogeneous but exhibit different mean values (Fig. B5c and d).

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

**Figure B4.** Intrinsic dimension estimated by the MLE (y axis) as a function of the dataset size (x axis) for geological models (a) and ground motion wavefields (b). For ground motion, the HEMEW-3D database is used for the fixed source (black line) and the HEMEW$^S$-3D database corresponds to the blue line.

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

**Figure B5.** Two pairs of geological models with a high SSIM of 0.6. (a, b) Geological models with a SSIM of 0.6 and normalized distance of 0.03 (c, d) Geological models with a SSIM of 0.6 and normalized distance of 0.15.
## Appendix C: Multiple-Input Fourier Neural Operator (MIFNO)

The Multiple-Input Fourier Neural Operator (MIFNO) architecture is shown in Fig. C1. Details about the model are given in Lehmann et al. (2024).

<figure><img src="image_22.png" /><figcaption>**Figure C1.** The MIFNO is made of a geology branch that encodes the geology with factorized Fourier (F-Fourier) layers and a source branch that transforms the vector of source parameters ($x_s, \theta_s$) into a 4D variable $v_S$ that matches the dimensions of the geology branch output, $v_K$. Outputs of each branch are concatenated after elementary mathematical operations and the remaining factorized Fourier layers are applied. Uplift, *P*, and projection $Q_E, Q_N$, and $Q_Z$ blocks are made of two dense layers with 128 neurons each.</figcaption></figure>

**Author contributions.** FL, FG, and DC designed the study. FL conducted the analyses. MB and DC conceived the original idea. FL wrote the paper with input from all authors.

**Competing interests.** The contact author has declared that none of the authors has any competing interests.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

**Acknowledgements.** We used the code from https://github.com/bakerjw/GMMs/tree/master (last access: 9 July 2024, Baker, 2022) to compute the GMMs. The authors are grateful for the resources and human support of the Très Grand Centre de Calcul (TGCC, CCRT, France). They thank the anonymous reviewers for their numerous suggestions that helped to improve the paper.

**Review statement.** This paper was edited by Andrea Rovida and reviewed by two anonymous referees.

## References

Aki, K. and Richards, P. G.: Quantitative Seismology, Freeman, San Francisco, 1980.

Altindal, A. and Askan, A.: SIGMOID-TR: A Simulated Ground Motion Dataset for Turkey, Zenodo [data set], https://doi.org/10.5281/ZENODO.7007917, 2022.

Annon: Simulate CO₂ Flow with Open Porous Media, Github [code], https://github.com/microsoft/AzureClusterlessHPC.jl/tree/main/examples/opm, 2022.

Arias, A.: A Measure of Earthquake Intensity, in: Seismic design for nuclear plants, edited by: Hansen, R. J., The MIT Press, 438–483, ISBN 978-0-262-08041-5, 1970.

Arroucau, P.: A Preliminary Three-Dimensional Seismological Model of the Crust and Uppermost Mantle for Metropolitan Franc, Tech. Rep. SIGMA2-2018-D2014, https://www sigma-2.net/medias/files/sigma2-2018-d2-014-3d-velocity-model-france-approved-public.pdf (last access: 12 April 2022), 2020.

Atkinson, G. M.: Ground-Motion Prediction Equation for Small-to-Moderate Events at Short Hypocentral Distances, with Application to Induced-Seismicity Hazards, B. Seismol. Soc. Am., 105, 981–992, https://doi.org/10.1785/0120140142, 2015.

Atkinson, G. M. and Boore, D. M.: Earthquake Ground-Motion Prediction Equations for Eastern North America, B. Seismol. Soc. Am., 96, 2181–2205, https://doi.org/10.1785/0120050245, 2006.

Bahrampouri, M., Rodriguez-Marek, A., Shahi, S., and Dwood, H.: An Updated Database for Ground Motion Parameters for KiK-net Records, Earthquake Spectra, 37, 505–522, https://doi.org/10.1177/8755293020952447, 2021.

Baker, J. W.: GMM, Github [code], https://github.com/bakerjw/GMMs/tree/master (last access: 9 July 2024), 2022.
Bonev, B., Kurth, T., Hundt, C., Pathak, J., Baust, M., and Kashinath, K.: Modelling Atmospheric Dynamics with Spherical Fourier Neural Operators, in: ICLR 2023 Workshop on Tackling Climate Change with Machine Learning, https://www.climatechange.ai/papers/iclr2023/47 (last access: 24 August 2023), 2023.

Castro-Cruz, D., Gatti, F., and Lopez-Caballero, F.: High-Fidelity Broadband Prediction of Regional Seismic Response: A Hybrid Coupling of Physics-Based Synthetic Simulation and Empirical Green Functions, Natural Hazards, 108, 1997–2031, https://doi.org/10.1007/s11069-021-04766-x, 2021.

CEA, CentraleSupélec, IPGP and CNRS: SEM3D, Github [code], https://github.com/sem3d/SEM (last access: 17 April 2024), 2017.

Chaljub, E., Celorio, M., Cornou, C., Martin, F. D., Haber, E. E., Margerin, L., Marti, J., and Zentner, I.: Numerical Simulation of Wave Propagation in Heterogeneous and Random Media for Site Effects Assessment in the Grenoble Valley, in: He 6th IASPEI/I-AEE International Symposium: Effects of Surface Geology on Seismic Motion, 2021.

Chernov, L. A.: Wave Propagation in a Random Medium, Translated by Richard A. Silverman, Mineola, New York, dover publications Edn., 1960.

Chiou, B. S.-J. and Youngs, R. R.: Update of the Chiou and Youngs NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra, Earthquake Spectra, 30, 1117–1153, https://doi.org/10.1193/072813EQS219M, 2014.

Colvez, M.: Influence of the Earth's Crust Heterogeneities and Complex Fault Structures on the Frequency Content of Seismic Waves, Ph.D. thesis, Université Paris Saclay, Paris-Saclay, https://theses.hal.science/tel-03551874 (last access: 11 March 2023), 2021.

Convertito, V., De Matteis, R., Amoroso, O., and Capuano, P.: Ground Motion Prediction Equations as a Proxy for Medium Properties Variation Due to Geothermal Resources Exploitation, Sci. Rep., 12, 12632, https://doi.org/10.1038/s41598-022-16815-x, 2022.

de Carvalho Paludo, L., Bouvier, V., and Cottereau, R.: Scalable Parallel Scheme for Sampling of Gaussian Random Fields over Very Large Domains: Parallel Scheme for Sampling of Random Fields over Very Large Domains, Int. J. Numer. Meth. Eng., 117, 845–859, https://doi.org/10.1002/nme.5981, 2019.

De Martin, F., Chaljub, E., Thierry, P., Sochala, P., Dupros, F., Maufroy, E., Hadri, B., Benaichouche, A., and Holender, F.: Influential Parameters on 3-D Synthetic Ground Motions in a Sedimentary Basin Derived from Global Sensitivity Analysis, Geophys. J. Int., 227, 1795–1817, https://doi.org/10.1093/gji/ggab304, 2021.

Deng, C., Feng, S., Wang, H., Zhang, X., Jin, P., Feng, Y., Zeng, Q., Chen, Y., and Lin, Y.: OpenFWI: Large-scale Multi-Structural Benchmark Datasets for Full Waveform Inversion, in: Advances in Neural Information Processing Systems, edited by: Koyejo, S., Mohamed, S., Agarwal, A., Belgrave, D., Cho, K., and Oh, A., Curran Associates, Inc., 35, 6007–6020, https://proceedings.neurips.cc/paper_files/paper/2022/file/27d3ef263c7cb8d542c4f9815a49b69b-Paper-Datasets_and_Benchmarks.pdf (last access: 24 August 2023), 2022.

Ding, Y., Chen, S., Li, X., Wang, S., Luan, S., and Sun, H.: Self-Adaptive Physics-Driven Deep Learning for Seismic Wave Mod-

eling in Complex Topography, Eng. Appl. Artif. Intel., 123, 106425, https://doi.org/10.1016/j.engappai.2023.106425, 2023.

Duverger, C., Mazet-Roux, G., Bollinger, L., Trilla, A. G., Vallage, A., Hernandez, B., and Cansi, Y.: A Decade of Seismicity in Metropolitan France (2010–2019): The CEA/LDG Methodologies and Observations, Bulletin de la Société Géologique de France, 192, p. 25, https://doi.org/10.1051/bsgf/2021014, 2021.

El Haber, E., Cornou, C., Jongmans, D., Lopez-Caballero, F., Youssef Abdelmassih, D., and Al-Bittar, T.: Impact of Spatial Variability of Shear Wave Velocity on the Lagged Coherency of Synthetic Surface Ground Motions, Soil Dyn. Earthq. Eng., 145, 106689, https://doi.org/10.1016/j.soildyn.2021.106689, 2021.

Equinor: Sleipner 2019 Benchmark Model, CO2DataShare [data set], https://doi.org/10.11582/2020.00004, 2020.

Faccioli, E., Maggio, F., Paolucci, R., and Quarteroni, A.: 2d and 3D Elastic Wave Propagation by a Pseudo-Spectral Domain Decomposition Method, J. Seismol., 1, 237–251, https://doi.org/10.1023/A:1009758820546, 1997.

Feng, S., Wang, H., Deng, C., Feng, Y., Liu, Y., Zhu, M., Jin, P., Chen, Y., and Lin, Y.: EFWI: Multiparameter Benchmark Datasets for Elastic Full Waveform Inversion of Geophysical Properties, in: Thirty-Seventh Conference on Neural Information Processing Systems Datasets and Benchmarks Track, https://openreview.net/forum?id=3BQaMV9jxK (last access: 10 June 2024), 2023.

Fu, H., He, C., Chen, B., Yin, Z., Zhang, Z., Zhang, W., Zhang, T., Xue, W., Liu, W., Yin, W., Yang, G., and Chen, X.: 18.9-Pflops 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, SC '17, Association for Computing Machinery, New York, NY, USA, ISBN 978-1-4503-5114-0, https://doi.org/10.1145/3126908.3126910, 2017.

Gadylshin, K., Lisitsa, V., Gadylshina, K., Vishnevsky, D., and Novikov, M.: Machine Learning-Based Numerical Dispersion Mitigation in Seismic Modelling, in: Computational Science and Its Applications – ICCSA 2021, edited by: Gervasi, O., Murgante, B., Misra, S., Garau, C., Blečić, I., Taniar, D., Apduhan, B. O., Rocha, A. M. A. C., Tarantino, E., and Torre, C. M., Springer International Publishing, Cham, 34–47, ISBN 978-3-030-86653-2, 2021.

Gatti, F. and Clouteau, D.: Towards Blending Physics-Based Numerical Simulations and Seismic Databases Using Generative Adversarial Network, Comput. Method. Appl. M., 372, 113421, https://doi.org/10.1016/j.cma.2020.113421, 2020.

Grady, T. J., Khan, R., Louboutin, M., Yin, Z., Witte, P. A., Chandra, R., Hewett, R. J., and Herrmann, F. J.: Model-Parallel Fourier Neural Operators as Learned Surrogates for Large-Scale Parametric PDEs, Comput. Geosci., 178, 105402, https://doi.org/10.1016/j.cageo.2023.105402, 2023.

Grassberger, P. and Procaccia, I.: Measuring the Strangeness of Strange Attractors, Physica D, 9, 189–208, 1983.

Hartzell, S., Harmsen, S., and Frankel, A.: Effects of 3D Random Correlated Velocity Perturbations on Predicted Ground Motions, B. Seismol. Soc. Am., 100, 1415–1426, https://doi.org/10.1785/0120090060, 2010.

Heinecke, A., Breuer, A., Rettenberger, S., Bader, M., Gabriel, A. A., Pelties, C., Bode, A., Barth, W., Liao, X. K., Vaidyanathan, K., Smelyanskiy, M., and Dubey, P.: Petascale High Order Dy-
namic Rupture Earthquake Simulations on Heterogeneous Su-
percomputers, in: SC '14: Proceedings of the International Con-
ference for High Performance Computing, Networking, Storage
and Analysis, New Orleans, LA, USA, ISBN 2167-4337, 3–14,
https://doi.org/10.1109/SC.2014.6, 2014.

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Fleming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.-N.: The ERA5 Global Reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049, https://doi.org/10.1002/qj.3803, 2020.

Imperatori, W. and Mai, P. M.: Broad-Band near-Field Ground Motion Simulations in 3-Dimensional Scattering Media, Geophys. J. Int., 192, 725–744, https://doi.org/10.1093/gji/ggs041, 2013.

Jessell, M., Guo, J., Li, Y., Lindsay, M., Scalzo, R., Giraud, J., Pirot, G., Cripps, E., and Ogarko, V.: Into the Noddyverse: a massive data store of 3D geological models for machine learning and inversion applications, Earth Syst. Sci. Data, 14, 381–392, https://doi.org/10.5194/essd-14-381-2022, 2022.

Karimpouli, S. and Tahmasebi, P.: Physics Informed Machine Learning: Seismic Wave Equation, Geosci. Front., 11, 1993–2001, https://doi.org/10.1016/j.gsf.2020.07.007, 2020.

Khazaie, S., Cottereau, R., and Clouteau, D.: Influence of the Spatial Correlation Structure of an Elastic Random Medium on Its Scattering Properties, J. Sound Vib., 370, 132–148, https://doi.org/10.1016/j.jsv.2016.01.012, 2016.

Komatitsch, D. and Tromp, J.: Introduction to the Spectral Element Method for Three-Dimensional Seismic Wave Propagation, Geophys. J. Int., 139, 806–822, https://doi.org/10.1046/j.1365-246x.1999.00967.x, 1999.

Lehmann, F.: Physics-based Simulations of 3D Wave Propagation with Source Variability: HEMEW\S-3D, Recherche Data Gouv [data set], https://doi.org/10.57745/LAI6YU, 2023.

Lehmann, F.: Lehmannfa/HEMEW3D: Initial Version, Zenodo [code], https://doi.org/10.5281/ZENODO.13625608, 2024.

Lehmann, F., Gatti, F., Bertin, M., and Clouteau, D.: Machine Learning Opportunities to Conduct High-Fidelity Earthquake Simulations in Multi-Scale Heterogeneous Geology, Front. Earth Sci., 10, https://doi.org/10.3389/feart.2022.1029160, 2022.

Lehmann, F., Gatti, F., and Clouteau, D.: Multiple-Input Fourier Neural Operator (MIFNO) for Source-Dependent 3D Elastodynamics, arXiv [preprint], https://doi.org/10.48550/arXiv.2404.10115, 2024.

Levina, E. and Bickel, P.: Maximum Likelihood Estimation of In-
trinsic Dimension, in: Advances in Neural Information Process-
ing Systems, edited by Saul, L., Weiss, Y., and Bottou, L., MIT
Press, vol. 17, https://proceedings.neurips.cc/paper_files/paper/
2004/file/74934548253bcab8490ebd74afed7031-Paper.pdf (last
access: 17 April 2023), 2004.

Liu, B., Yang, S., Ren, Y., Xu, X., Jiang, P., and Chen, Y.: Deep-Learning Seismic Full-Waveform Inversion for Realistic Structural Models, Geophysics, 86, R31–R44, https://doi.org/10.1190/geo2019-0435.1, 2021.

Maechling, P. J., Silva, F., Callaghan, S., and Jordan, T. H.: SCEC Broadband Platform: System Architecture and Software Implementation, Seismol. Res. Lett., 86, 27–38, https://doi.org/10.1785/0220140125, 2015.

Mansoor, K., Buscheck, T., Yang, X., Carroll, S., and Chen, X.: LLNL Kimberlina 1.2 NUFT Simulations June 2018 (V2), https://doi.org/10.18141/1603336, 2020.

Michelini, A., Cianetti, S., Gaviano, S., Giunchi, C., Jozinović, D., and Lauciani, V.: INSTANCE – the Italian seismic dataset for machine learning, Earth Syst. Sci. Data, 13, 5509–5544, https://doi.org/10.5194/essd-13-5509-2021, 2021.

Moczo, P., Kristek, J., Bard, P.-Y., Stripajová, S., Hollender, F., Chovanová, Z., Kristeková, M., and Sicilia, D.: Key Structural Parameters Affecting Earthquake Ground Motion in 2D and 3D Sedimentary Structures, B. Earthq. Eng., 16, 2421–2450, https://doi.org/10.1007/s10518-018-0345-5, 2018.

Molinari, I. and Morelli, A.: EPcrust: A Reference Crustal Model for the European Plate: EPcrust, Geophys. J. Int., 185, 352–364, https://doi.org/10.1111/j.1365-246X.2011.04940.x, 2011.

Moseley, B., Markham, A., and Nissen-Meyer, T.: Solving the Wave Equation with Physics-Informed Deep Learning, ArXiv [preprint], https://doi.org/10.48550/arxiv.2006.11894, 2020.

Mousavi, S. M. and Beroza, G. C.: Machine Learning in Earthquake Seismology, Annu. Rev. Earth Pl. Sc., 51, 105–129, https://doi.org/10.1146/annurev-earth-071822-100323, 2023.

Mousavi, S. M., Sheng, Y., Zhu, W., and Beroza, G. C.: Stanford Earthquake Dataset (STEAD): A Global Data Set of Seismic Signals for AI, IEEE Access, 7, 179464–179476, https://doi.org/10.1109/ACCESS.2019.2947848, 2019.

Ovadia, O., Kahana, A., Stinis, P., Turkel, E., and Karniadakis, G. E.: ViTO: Vision Transformer-Operator, ArXiv [preprint], https://doi.org/10.48550/arXiv.2303.08891, 2023.

Paolucci, R., Smerzini, C., and Vanini, M.: BB-SPEEDset: A Validated Dataset of Broadband Near-Source Earthquake Ground Motions from 3D Physics-Based Numerical Simulations, B. Seismol. Soc. Am., 111, 2527–2545, https://doi.org/10.1785/0120210089, 2021.

Pathak, J., Subramanian, S., Harrington, P., Raja, S., Chattopadhyay, A., Mardani, M., Kurth, T., Hall, D., Li, Z., Azizzadenesheli, K., Hassanzadeh, P., Kashinath, K., and Anandkumar, A.: FourCastNet: A Global Data-driven High-resolution Weather Model Using Adaptive Fourier Neural Operators, ArXiv [preprint], https://doi.org/10.48550/arXiv.2202.11214, 2022.

Petrone, F., Abrahamson, N., McCallen, D., Pitarka, A., and Rodgers, A.: Engineering Evaluation of the EQSIM Simulated Ground-motion Database: The San Francisco Bay Area Region, Earthquake Eng. Struc., 50, 3939–3961, https://doi.org/10.1002/eqe.3540, 2021.

Poursartip, B., Fathi, A., and Tassoulas, J. L.: Large-Scale Simulation of Seismic Wave Motion: A Review, Soil Dyn. Earthq. Eng., 129, 105909, https://doi.org/10.1016/j.soildyn.2019.105909, 2020.

Qiu, H., Yang, Y., and Pan, H.: Underestimation Modification for Intrinsic Dimension Estimation, Pattern Recogn., 140, 109580, https://doi.org/10.1016/j.patcog.2023.109580, 2023.

Raissi, M., Perdikaris, P., and Karniadakis, G.: Physics-Informed Neural Networks: A Deep Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Par-
tial Differential Equations, J. Comput. Phys., 378, 686-707,
https://doi.org/10.1016/j.jcp.2018.10.045, 2019.

Rasht-Behesht, M., Huber, C., Shukla, K., and Karniadakis, G. E.: Physics-Informed Neural Networks (PINNs) for Wave Propagation and Full Waveform Inversions, J. Geophys. Res.-Sol. Ea., 127, e2021JB023120, https://doi.org/10.1029/2021JB023120, 2022.

Rekoske, J. M., Gabriel, A.-A., and May, D. A.: Instantaneous Physics-Based Ground Motion Maps Using Reduced-Order Modeling, J. Geophys. Res.-Sol. Ea., 128, e2023JB026975, https://doi.org/10.1029/2023JB026975, 2023.

Ren, P., Rao, C., Chen, S., Wang, J.-X., Sun, H., and Liu, Y.: Seis-
micNet: Physics-informed Neural Networks for Seismic Wave
Modeling in Semi-Infinite Domain, Comput. Phys. Commun.,
295, https://doi.org/10.1016/j.cpc.2023.109010, 2024.

Ross, D. A., Lim, J., Lin, R.-S., and Yang, M.-H.: Incremental
Learning for Robust Visual Tracking, Int. J. Comput. Vis., 77,
125-141, https://doi.org/10.1007/s11263-007-0075-7, 2008.

Rosti, A., Smerzini, C., Paolucci, R., Penna, A., and Rota, M.: Val-
idation of Physics-Based Ground Shaking Scenarios for Empiri-
cal Fragility Studies: The Case of the 2009 L'Aquila Earthquake,
B. Earthquake Eng., 21, 95-123, https://doi.org/10.1007/s10518-
022-01554-1, 2023.

Scalise, M., Pitarka, A., Louie, J. N., and Smith, K. D.: Ef-
fect of Random 3D Correlated Velocity Perturbations on
Numerical Modeling of Ground Motion from the Source
Physics Experiment, B. Seismol. Soc. Am., 111, 139-156,
https://doi.org/10.1785/0120200160, 2021.

Shahjouei, A. and Pezeshk, S.: Alternative Hybrid Empirical
Ground-Motion Model for Central and Eastern North America
Using Hybrid Simulations and NGA-West2 Models, B. Seismol.
Soc. Am., 106, 734-754, https://doi.org/10.1785/0120140367,
2016.

Shinozuka, M. and Deodatis, G.: Simulation of Stochastic Pro-
cesses by Spectral Representation, Appl. Mech. Rev., 44, 191-
204, https://doi.org/10.1115/1.3119501, 1991.

Smerzini, C., Paolucci, R., and Stupazzini, M.: Comparison of 3D,
2D and 1D Numerical Approaches to Predict Long Period Earth-
quake Ground Motion in the Gubbio Plain, Central Italy, B.
Earthquake Eng., 9, 2007-2029, https://doi.org/10.1007/s10518-
011-9289-8, 2011.

Song, C., Liu, Y., Zhao, P., Zhao, T., Zou, J., and Liu, C.: Simulating Multi-Component Elastic Seismic Wavefield Using Deep Learning, IEEE Geosci. Remote Sens. Lett., 20, 1-5, https://doi.org/10.1109/LGRS.2023.3250522, 2023.

Ta, Q.-A., Clouteau, D., and Cottereau, R.: Modeling of Random Anisotropic Elastic Media and Impact on Wave Propagation, Eur. J. Comput. Mech., 19, 241–253, https://doi.org/10.3166/ejcm.19.241-253, 2010.

Tarbali, K. and Bradley, B.: Ground Motion Selection for Sce-
nario Ruptures Using the Generalised Conditional Intensity Mea-
sures (GCIM) Method, Earthq. Eng. Struct. D., 44, 1601–1621,
https://doi.org/10.1002/eqe.2546, 2015.

Thompson, E. M., Baise, L. G., and Kayen, R. E.: Spatial Correlation of Shear-Wave Velocity in the San Francisco Bay Area Sediments, Soil Dyn. Earthq. Eng., 27, 144–152, https://doi.org/10.1016/j.soildyn.2006.05.004, 2007.

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

Wang, Z., Bovik, A., Sheikh, H., and Simoncelli, E.: Image Quality Assessment: From Error Visibility to Structural Similarity, IEEE T. Image Process., 13, 600–612, https://doi.org/10.1109/TIP.2003.819861, 2004.

Wen, G., Li, Z., Long, Q., Azizzadenesheli, K., Anandkumar, A., and Benson, S. M.: Real-Time High-Resolution CO₂ Geological Storage Prediction Using Nested Fourier Neural Operators, Energ. Environ. Sci., 16, 1732–1741, https://doi.org/10.1039/d2ee04204e, 2023.

Witte, P. A., Konuk, T., Skjetne, E., and Chandra, R.: Fast CO2 Saturation Simulations on Large-Scale Geomodels with Artificial Intelligence-Based Wavelet Neural Operators, Int. J. Greenhouse Gas Control, 126, 103880, https://doi.org/10.1016/j.ijggc.2023.103880, 2023.

Wu, Y., Aghamiry, H. S., Operto, S., and Ma, J.: Helmholtz Equation Solution in Non-Smooth Media by Physics-Informed Neural Network with Incorporating Quadratic Terms and a Perfectly Matching Layer Condition, Geophysics, 88, 1–66, https://doi.org/10.1190/geo2022-0479.1, 2023.

Zhang, T., Trad, D., and Innanen, K.: Learning to Solve the Elastic Wave Equation with Fourier Neural Operators, Geophysics, 88, 1–63, https://doi.org/10.1190/geo2022-0268.1, 2023.

Zhu, C., Riga, E., Pitilakis, K., Zhang, J., and Thambiratnam, D.: Seismic Aggravation in Shallow Basins in Addition to One-dimensional Site Amplification, J. Earthquake Eng., 24, 1477–1499, https://doi.org/10.1080/13632469.2018.1472679, 2020.

Zhu, W., Hou, A. B., Yang, R., Datta, A., Mousavi, S. M., Ellsworth, W. L., and Beroza, G. C.: QuakeFlow: A Scalable Machine-Learning-Based Earthquake Monitoring Workflow with Cloud Computing, Geophys. J. Int., 232, 684–693, https://doi.org/10.1093/gji/ggac355, 2022.