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

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^S-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}=\mathrm{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 frameworks, 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 Earthquake 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; Poursartip 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. Gadylshin 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 HEMEW^S-3D (HEterogeneous Materials and Elastic Waves with Source variability in 3D) database contains 30 000 high-fidelity simulations in 3D domains of $\mathrm{9.6}\times \mathrm{9.6}\times \mathrm{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.

(Jessell et al., 2022)(Deng et al., 2022)(Deng et al., 2022)(Deng et al., 2022)(Feng et al., 2023)

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.

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 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 𝔼^FWI 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^S-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∈ℝ³, that propagates in a 3D Euclidean space. We consider a truncated propagation domain, $\mathrm{\Omega }=[\mathrm{0},L{]}^{\mathrm{3}}$, with absorbing boundary conditions all around, except the traction-free top surface, and a solution, $\mathit{u}:\mathrm{\Omega }\times [\mathrm{0},T]\to {\mathbb{R}}^{\mathrm{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

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathit{\rho }{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{u}}{\partial t^{\mathrm{2}}}}& =\mathrm{\nabla }\mathit{\lambda }\left(\mathrm{\nabla }\cdot \mathit{u}\right)+\mathrm{\nabla }\mathit{\mu }\left[\mathrm{\nabla }\mathit{u}+{\left(\mathrm{\nabla }\mathit{u}\right)}^T\right]\\ & +\left(\mathit{\lambda }+\mathrm{2}\mathit{\mu }\right)\mathrm{\nabla }\left(\mathrm{\nabla }\cdot \mathit{u}\right)-\mathit{\mu }\mathrm{\nabla }\times \mathrm{\nabla }\times \mathit{u}+\mathit{f},\end{array}\end{array}$$

where $\mathit{\rho }:\mathrm{\Omega }\to \mathbb{R}$ is the material unit mass density; $\mathit{\lambda }:\mathrm{\Omega }\to \mathbb{R}$ and $\mathit{\mu }:\mathrm{\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 ρ, λ, and μ 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:

$$\begin{array}{}\text{(2)}& \mathcal{L}(V_{\mathrm{S}},\mathit{u})=\mathit{f}.\end{array}$$

3.2 Earthquake source

In our database, the forcing term $\mathit{f}(\mathit{x},t)=\mathrm{div}\mathit{m}\left(\mathit{x}\right)\cdot s\left(t\right)$ 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_{\mathrm{s}},y_{\mathrm{s}},z_{\mathrm{s}})\in \mathrm{\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{array}{}\text{(3)}& \begin{array}{rl}{x}_{\mathrm{s}}& \in [\mathrm{1.2},\mathrm{8.4}]\mathrm{km},\\ y_{\mathrm{s}}& \in [\mathrm{1.2},\mathrm{8.4}]\mathrm{km},\\ z_{\mathrm{s}}& \in [-\mathrm{9.0},-\mathrm{0.6}]\mathrm{km}.\end{array}\end{array}$$

In addition to the source position, the source is parametrized by the symmetric 3×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_{\mathrm{0}}=\mathrm{2.47}\times {\mathrm{10}}^{\mathrm{16}}$ N m, and the source time evolution is a spice bench given by $s\left(t\right)=\mathrm{1}-\left(\mathrm{1}+\frac{t}{\mathit{\tau }}\right)e^{-\frac{t}{\mathit{\tau }}}$, with τ=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 HEMEW^S-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₀. For the reference source time function, s(t), the solution, u(x,t), provided in HEMEW^S-3D can be written as the convolution of the Green function with the source time function, i.e. $\mathit{u}(\mathit{x},t)=\mathbf{G}(\mathit{x},t)*s\left(t\right)$.

Computing the solution, u₁, for a new source time function, s₁ (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, $\widehat{s}:=\mathcal{F}\left(s\right)$, and the solution, $\widehat{\mathit{u}}:=\mathcal{F}\left(u\right)$;

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

3. compute the Fourier transform of the new source time function, ${\widehat{s}}_{\mathrm{1}}$;

4. compute the new solution in the frequency domain, ${\widehat{\mathit{u}}}_{\mathrm{1}}=\widehat{\mathbf{G}}*{\widehat{s}}_{\mathrm{1}}$;

5. deduce the new solution in the temporal domain, ${\mathit{u}}_{\mathrm{1}}={\mathcal{F}}^{-\mathrm{1}}\left({\widehat{\mathit{u}}}_{\mathrm{1}}.\right)$

From these remarks, one should remember that ground motion wavefields in the HEMEW^S-3D database originate from pointwise sources with a different location, x_s∈ℝ³, and orientation, θ_s∈ℝ³, but the same source time function.

3.3 Heterogeneous geological models

The HEMEW^S-3D database contains samples $\mathit{\{}{V}_{\mathrm{S}}^{\left(i\right)},{\mathit{x}}_{\mathrm{s}}^{\left(i\right)},{\mathit{\theta }}_{\mathrm{s}}^{\left(i\right)},{\dot{\mathit{u}}}^{\left(i\right)}{\mathit{\}}}_i$ that satisfy Eq. (2) ($\dot{\mathit{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.

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.

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_{\mathrm{S},\max }=\mathrm{4500}$ m s⁻¹. The minimum S-wave velocity is $V_{\mathrm{S},\min }=\mathrm{1071}$ m s⁻¹. Above the bottom layer, the number of horizontal layers and their thickness are randomly chosen for each sample $V_{\mathrm{S}}^{\left(i\right)}$, 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, 𝒰([μ₁,μ₂]). All layerwise values are chosen independently. Values of ${\mathit{\mu }}_{\mathrm{1}}=V_{\mathrm{S},\min }/\mathrm{0.6}$ = 1785 m s⁻¹ and ${\mathit{\mu }}_{\mathrm{2}}=V_{\mathrm{S},\max }/\mathrm{1.4}=\mathrm{3214}$ m s⁻¹ were determined to ensure that most values remain bounded within $[V_{\mathrm{S},\min },V_{\mathrm{S},\max }]$ after the addition of random fields in each layer (see Table 2 for a summary of the parameters).

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.

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

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathit{\rho }& =\mathrm{1.6612}{V}_{\mathrm{P}}-\mathrm{0.4721}{V}_{\mathrm{P}}^{\mathrm{2}}+\mathrm{0.0671}{V}_{\mathrm{P}}^{\mathrm{3}}\\ & -\mathrm{0.0043}{V}_{\mathrm{P}}^{\mathrm{4}}+\mathrm{0.000106}{V}_{\mathrm{P}}^{\mathrm{5}}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(5)}& Q_{\mathrm{P}}=max\left({\displaystyle \frac{V_{\mathrm{P}}}{\mathrm{20}}},{\displaystyle \frac{V_{\mathrm{S}}}{\mathrm{5}}}\right),Q_{\mathrm{S}}={\displaystyle \frac{V_{\mathrm{S}}}{\mathrm{10}}}.\end{array}$$

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, $\left|\mathcal{N}\right(\mathrm{0.2},\mathrm{0.1}\left)\right|$, 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).

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, ℛ, can be decomposed as a sum of independent identically distributed random variables $(V_{\mathrm{S},n}{)}_{-N\le n\le N}$, with a uniform distribution over [0,2π].

$$V_{\mathrm{S}}\left(x\right)=\sum _{n=-N}^N\sqrt{\mathrm{2}\widehat{\mathcal{R}}\left(n\mathrm{\Delta }k\right)}\cos (n\mathrm{\Delta }k\cdot \times +V_{\mathrm{S},n}),$$

where $\widehat{\mathcal{R}}$ is the Fourier transform of the autocovariance function ℛ and Δk is the unit volume in ℝ³.

Finally, V_S values are clipped between $V_{\mathrm{S},\min }=\mathrm{1071}$ m s⁻¹ and $V_{\mathrm{S},\max }=\mathrm{4500}$ m s⁻¹. 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_{\mathrm{S}}^{\left(i\right)}$, is discretized over a grid of $\mathrm{32}\times \mathrm{32}\times \mathrm{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, ${\mathit{\lambda }}_{min}=V_{\mathrm{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{\mathrm{7}}{\mathrm{5}}\cdot \frac{V_{\mathrm{S},\min }}{f_{max}}=\mathrm{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. $\mathrm{2}\times f_{max}=\mathrm{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.

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.

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 $\mathrm{32}\times \mathrm{32}\times \mathrm{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.

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.

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.

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.

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

$$\begin{array}{}\text{(6)}& I_{\mathrm{A}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}g}}\underset{\mathrm{0}}{\overset{T}{\int }}{a}^{\mathrm{2}}\left(t\right)\mathrm{d}t,\end{array}$$

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 HEMEW^S-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-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 HEMEW^S-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).

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.

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.

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.

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

$$\begin{array}{}\text{(7)}& E_{\mathrm{s}}={\displaystyle \frac{M_{\mathrm{0}}\mathrm{\Delta }\mathit{\sigma }}{\mathrm{2}\mathit{\mu }}},\end{array}$$

where M₀ is the seismic moment, Δσ is the stress drop, and μ is the shear modulus at the fault location. Knowing that the shear wave velocity is written as $V_{\mathrm{S}}=\sqrt{\mathit{\mu }/\mathit{\rho }}$, Eq. (7) indicates that the seismic energy is inversely proportional to $V_{\mathrm{S}}^{\mathrm{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 HEMEW^S-3D database is between $\mathrm{2.3}\times {\mathrm{10}}^{-\mathrm{3}}$ and 81.2 g. It decreases to the interval between $\mathrm{5.8}\times {\mathrm{10}}^{-\mathrm{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.

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.

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 HEMEW^S-3D database is compared with four GMMs (all taken with a moment magnitude of 4.9 that corresponds to the seismic moment of the HEMEW^S-3D database):

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

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_{\mathrm{S},\mathrm{30}}=\mathrm{2000}$ m s⁻¹, 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 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.

Figure 8 shows that the mean PSA computed from the HEMEW^S-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 depth, especially for the smallest hypocentral distances. The HEMEW^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 $\mathrm{32}\times \mathrm{32}\times \mathrm{32}$ = 32 768 points, and ground motions contain $\mathrm{16}\times \mathrm{16}\times \mathrm{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 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.

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.

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 $\mathit{\{}{V}_{\mathrm{S}}^{\left(i\right)}{\mathit{\}}}_{\mathrm{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_{\mathrm{S}}^{\left(i\right)},V_{\mathrm{S}}^{\left(j\right)}{)}_{i\ne j}$ that are at a distance lower than r:

$$\begin{array}{}\text{(8)}& C_N\left(r\right)={\displaystyle \frac{\mathrm{2}}{N(N-\mathrm{1})}}\sum _{i=\mathrm{1}}^N\sum _{j=i+\mathrm{1}}^N\mathbb{1}\left(\Vert V_{\mathrm{S}}^{\left(i\right)}-V_{\mathrm{S}}^{\left(j\right)}\Vert \le r\right).\end{array}$$

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.

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

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

$$\begin{array}{}\text{(9)}& \mathrm{SSIM}(A,B)={\displaystyle \frac{(\mathrm{2}{\mathit{\mu }}_A{\mathit{\mu }}_B+C_{\mathrm{1}})(\mathrm{2}{\mathit{\sigma }}_{AB}+C_{\mathrm{2}})}{({\mathit{\mu }}_A^{\mathrm{2}}+{\mathit{\mu }}_B^{\mathrm{2}}+C_{\mathrm{1}})({\mathit{\sigma }}_A^{\mathrm{2}}+{\mathit{\sigma }}_B^{\mathrm{2}}+C_{\mathrm{2}})}},\end{array}$$

where μ_A and μ_B are the means of A and B, σ_A and σ_B are the unbiased estimators of the variance of A and B, σ_AB is the unbiased estimator of the covariance of A and B, and C₁ and C₂ 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 distance (Fig. B5a), while the second geologies have different mean values, leading to a higher Euclidean distance (Fig. B5b).

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.

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 HEMEW^S-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 different from the random fields contained in the HEMEW^S-3D database – for instance, embedding a basin shape.

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

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 HEMEW^S-3D database from a geometrical point of view. To match the design of the HEMEW^S-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 HEMEW^S-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).

Figure 11 For two spatial points, the velocity field predicted by the MIFNO (dashed red line) is compared with the reference from the HEMEW^S-3D database (black line).

For each geological model and source in the HEMEW^S-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 HEMEW^S-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 HEMEW^S-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 HEMEW^S-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 HEMEW^S-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 HEMEW^S-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_{\mathrm{S}}\times \mathit{\rho }$) in one geological layer, Z_ℓ, and the layer above, Z_ℓ−1. Figure A2 illustrates the distribution of impedance contrasts and shows that they are 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_{\mathrm{P}}/V_{\mathrm{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 generate the HEMEW^S-3D database is given at https://github.com/lehmannfa/HEMEW3D (Lehmann, 2024).

8 Conclusions

We presented the HEMEW^S-3D (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 $\mathrm{9.6}\times \mathrm{9.6}\times \mathrm{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 $\mathrm{32}\times \mathrm{32}\times \mathrm{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×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×10⁵ 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 HEMEW^S-3D database in complementing databases of recorded earthquakes. Due to the size of individual samples in the HEMEW^S-3D 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 HEMEW^S-3D 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 HEMEW^S-3D 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 A1 The seismic moment function in the HEMEW^S-3D database.

Figure A2 Distribution of impedance contrasts in the HEMEW^S-3D database. The impedance contrast is computed as the ratio between impedance in one layer $(Z_{\mathrm{\ell }}=V{s}_{\mathrm{\ell }}\times {\mathit{\rho }}_{\mathrm{\ell }})$ and impedance in the layer above it (Z_ℓ−1).

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.

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.

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^S-3D database corresponds to the blue line.

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 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 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 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 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 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,θ_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.