# Neural Network-Based Metamodel of synthetic seismograms: Application for uncertainty quantification

Mohamad Ali Noureddine, Florent de Martin, Rani El Meouche, Muhammad Ali Sammuneh, Fakhreddine Ababsa, Mickaël Beaufils

## ► To cite this version:

Mohamad Ali Noureddine, Florent de Martin, Rani El Meouche, Muhammad Ali Sammuneh, Fakhreddine Ababsa, et al.. Neural Network-Based Metamodel of synthetic seismograms: Application for uncertainty quantification. Engineering Applications of Artificial Intelligence, 2025, 151, pp.110613. 10.1016/j.engappai.2025.110613. hal-05037051

HAL Id: hal-05037051

https://brgm.hal.science/hal-05037051v1

Submitted on 16 Apr 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>
Research paper

# Neural Network-Based Metamodel of synthetic seismograms: Application for uncertainty quantification

Mohamad Ali Noureddine<sup>a,b,c,*</sup>, Florent De Martin<sup>b</sup>, Rani El Meouche<sup>a,D</sup>,
Muhammad Ali Sammuneh<sup>a,D</sup>, Fakhreddine Ababsa<sup>c,D</sup>, Mickael Beaufils<sup>b</sup>

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

## ARTICLE INFO

### Keywords:

- Neural networks

- Metamodel

- Uncertainty quantification

- Elastic wave propagation

- Earthquake

## ABSTRACT

We developed time and frequency domains neural network surrogate models of synthetic seismograms stemming from the resolution of the three-dimensional equation of motion. The surrogates predict the time or frequency-series knowing the input variables of the physical model. The surrogate models are used to quantify epistemic uncertainty in ground motion prediction through global sensitivity analysis. We first generate a dataset of time domain synthetic seismograms in a seven-dimensions uncertain space using the spectral-element method. To validate the surrogate model, we evaluate its ability to reproduce at least 80% of the bootstrap resamples. Additionally, the $R_2$ regression coefficient between the simulations generated by the spectral-element method code and those predicted by the neural network is 0.94 for the validation set, confirming the accuracy of the surrogate model. These surrogates allow fast predictions of velocity time-series or Fourier amplitude spectra where spectral-element simulations are not done (neural networks compute about 100,000 surrogates per second, while a single spectral-element simulation longs approximately 7 h on 48 cores). To quantify the uncertainty of the physical system under study, a global sensitivity analysis is undertaken to better understand how uncertain parameters affect the predicted state of the system. We present two sampling-based estimation methods, the so-called “pick-freeze” Sobol method and the Li and Mahadevan method, to quantify this uncertainty. The Sobol method requires approximately 500,000 simulations to achieve stability, whereas the Li and Mahadevan method requires only 30,000 simulations. Using the metamodel, both methods require only a few seconds to produce results, although the Li and Mahadevan method analyzes 10,000 simulations per second, compared to 3,000 for the Sobol method. The results indicate that the shear wave velocity in the physical system’s layer is the most influential parameter affecting the ground speed on the physical system. In contrast, at a reference station (without considering the geological properties of the physical system), the results show that the shear wave velocity in the first and third deep layers are the most influential parameters.

## 1. Introduction

Many physical systems are mathematically modeled by sets of partial differential equations (PDEs). Examples include diffusion processes (Meyer et al., 1974) or wave propagation (Chapman, 2004), among others. A major challenge for predicting reliable physical quantities arises from the presence of uncertainties in the systems due to systems' inherent variability and/or to our limited knowledge of them. The need to describe and quantify uncertainties in complex systems makes the field of uncertainty quantification (UQ) (Berthier, 2005) a fundamental component of predictive science, enabling robust predictability of PDE models for complex systems. UQ problems are

often formulated in a probabilistic framework, where uncertain parameters of the models are represented by probability distributions and random processes. Solving the UQ problem then amounts to solving a system of stochastic PDEs with random parameters (Du et al., 2023; Kolla et al., 2022). In the current work, we are concerned with direct UQ problems involving multiple realizations of the models using a computationally intensive solver based on finite-difference or finite-element methods. Specifically, given the distribution of the system's uncertain parameters, we compute the statistics of certain desired quantities of interest (QoIs). Various approaches exist to address direct UQ, among which global sensitivity analysis (Iooss, 2011) stands out.

<footnote><sup>a</sup> Institut de Recherche, ESTP, 28 Avenue du Président Wilson, Cachan, F-94230, France</footnote>

<footnote><sup>b</sup> Bureau de Recherches Géologiques et Minières, Orléans, 45100, France</footnote>

<footnote><sup>c</sup> PIMM, Arts et Métiers ParisTech, CNRS, CNAM, HESAM University, Paris, 75013, France</footnote>

<footnote>* Corresponding author at: Institut de Recherche, ESTP, 28 Avenue du Président Wilson, Cachan, F-94230, France.
*E-mail address: mnoureddine@estp.fr (M.A. Noureddine).</footnote>
It allows us to determine which model parameters contribute the most to a given quantity of interest predicted by the model, which have no influence, and which interact within the model. Hereafter, we will refer to the input variables as “inputs” and the numerical model responses as “outputs”. The sensitivity analysis of a model has long been based on evaluating the impact on output values of small perturbations in input values around a nominal point. This deterministic approach involves calculating or estimating indices based on the partial derivatives of the model at a specific point (Kleijnen, 1997). Thanks to adjoint approaches that allow for rapid handling of models with a large number of input variables, this type of approach is commonly used in solving large environmental systems (climatology, oceanography, hydrogeology) (Cacuci, 1981; Castaings, 2007). From the late 1980s, to overcome the limitations of local methods (assumptions of linearity and normality, local variations), a new generation of methods developed in a statistical framework emerged. Unlike local methods, these are termed “global” as they consider the entire domain of possible input variable variations (Saltelli et al., 2006). The objectives of these methods are to identify and prioritize the most influential inputs. In other words, if it is possible to obtain a sufficiently large sample of “input-output” simulations, it is easy to obtain sensitivity indices of the output concerning the inputs by applying linear regression techniques in the case of a linear model or functional variance decomposition in the case of nonlinear models (Dai et al., 2022). These indices quantify the influence of each input's uncertainty on the output. In this context, one of the main challenges in global sensitivity analysis is the high number of model evaluations needed to estimate parameter sensitivities, especially when models are expensive to evaluate (Sochala et al., 2020; De Martin et al., 2021; Schmitz et al., 2023). This becomes an obstacle in traditional approaches, such as the Monte Carlo method (Hammersley, 2013) or the variance-based method (Fajraoui, 2014), which often require a large number of repeated model calls for each parameter. In recent years, neural networks (Beale et al., 2010) have demonstrated remarkable success in solving a variety of artificial intelligence problems (Schmidhuber, 2015). For example, the authors in Cao et al. (2020) use a spectral temporal graph neural network to model temporal and spatial dependencies in multivariate time-series forecasting, while in the article (Song et al., 2020), the authors employ Long Short-Term Memory (LSTM) networks to capture long-term dependencies for predicting well performance based on historical time series. Our approach differs significantly from those of Cao et al. (2020) and Song et al. (2020), our method involves training a neural network on an existing dataset to predict time series for new configurations, capturing the complex relationship between the input parameters of the wave propagation model and the resulting time series, whereas the other methods focus more on modeling complex temporal and spatial dependencies.

In neural networks, uncertainty quantification primarily addresses two types of uncertainties: epistemic uncertainty (Hüllermeier and Waegeman, 2021), which is associated with a lack of knowledge or insufficient data, and aleatoric uncertainty, which arises from the intrinsic noise in the data (Kendall and Gal, 2017). Addressing these uncertainties not only enhances the interpretability of models but also ensures safer and more reliable usage in practical contexts. For example, various approaches have been proposed to manage uncertainties and delays in complex systems. The authors in Cao et al. (2024b) explored exponential state estimation for time-delay competitive neural networks subject to stochastic sampled controls and Markov jumps, while accounting for actuator failures. The authors also explored input-state stability in genetic networks with stochastic uncertainties and Markovian transitions (Cao et al., 2024a). However, in seismic risk analysis, UQ helps quantify the variability in ground motion predictions caused by geological uncertainties, which is crucial for designing resilient infrastructures. In this context, the neural networks can, for instance, be used to construct surrogate models for physical quantities

of interest, thus solving epistemic uncertainty quantification problems (Hu et al., 2022; Tan et al., 2023). The data necessary to construct (or train) a neural network as a surrogate model for a quantity of interest must be obtained by computing realizations of this quantity, where each realization involves solving a PDE problem via numerical discretization. In this framework, researchers have proposed different UQ methods and examined their performance in applications across diverse domains, such as computer vision (Kendall and Gal, 2017), image processing (Lambert et al., 2024). These applications highlight the broad importance of UQ in enhancing the robustness and reliability of machine learning models. This growing body of research on UQ in deep learning, as reviewed in Abdar et al. (2021), underscores its significance in ensuring that neural network models provide reliable and interpretable predictions, especially in high-stakes decision-making scenarios.

In this paper, we focus on predicting ground velocity field at the Dumanoir dam located in Guadeloupe, France. This field is described by the partial differential equation of motion (Chapman, 2004), and solved using EFISPEC3D code (De Martin, 2011) applying a spectral-element method (Komatitsch, 1997). To be solved, the equations of motion require velocity structure, seismic source information, as well as specific geological information about the region under study. Regretably, the uncertainty surrounding these parameters makes it difficult to accurately reproduce wave propagation in complex geological media. To quantify the impact of this uncertainty, a global sensitivity analysis (Cacuci, 1981) is undertaken to estimate the sensitivity indices of specific parameters, helping us better understand how uncertain parameters influence the predicted state of the system. To estimate these indices, we use two sampling-based methods: the first developed by Sobol (Soboł, 1993; Sobol, 2001; Jacques, 2011) and the second by Li and Mahadevan (2016). To obtain precise estimates of sensitivity indices, these methods are costly in terms of the number of model evaluations. According to Jiménez Rugama and Gilquin (2018), estimating a Sobol index would require 10,000 evaluations of the direct model to achieve a precision of 10%, for each input variable. To address this problem, a regression-based neural network is used to replace the computationally expensive direct model.

The structure of the paper is as follows. In Section 2, we present the site of study and its numerical replica. In Section 3, we explain the concept of neural networks and demonstrate their ability to replace the direct solver. In Section 4, we introduce global sensitivity analysis, detailing the two methods used, and then present the results of the sensitivity analysis using the bootstrap technique (Wehrens et al., 2000) to assess the stability of the estimates. Finally, conclusions and future perspectives are presented. A block diagram of the proposed research (step by step) is given by Fig. 1.

## 2. The Dumanoir dam: overview and numerical model

### 2.1. Site of study

We focus on simulating the displacement field generated by earthquakes of various magnitudes for the Dumanoir dam (see Fig. 2) located in Guadeloupe, in the French Antilles. This is a critical hydraulic infrastructure that regulates watercourses and contributes to energy production in the region. The area is situated amidst a combination of geological formations and is subject to various terrestrial phenomena, including the volcanic activity of La Grande Soufrière and tectonic movements resulting from the subduction of the Caribbean plate beneath the American plate. The dam is monitored by four stations distributed as shown in Fig. 3, a photo of a sensor is presented in Fig. 4. In the following, we will compare the signals recorded at the three site stations with the signals recorded at the reference station, which will allow us to deduce the effect of the dam on ground motion.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Block diagram of the proposed research.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Location of the Dumanoir dam (Guadeloupe). Maps data: ©Google. Photo: ©AnteaGroup.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Seismic sensors around the Dumanoir dam. Sensor 1: Site station, Sensor 2: Site station, Sensor 3: Site station, Sensor 4: Reference station.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Photo of the sensor deployed on the Dumanoir dam (Sensor 1).</figcaption></figure>

## 2.2. Theoretical model and numerical method

The ground displacement field is described by the equation of motion, which governs the propagation of seismic waves in a continuous medium. This mathematical equation is a partial differential equation (Chapman, 2004) formulated as follows:

$$ \rho \frac{\partial^2 u}{\partial t^2} = \nabla \cdot T + f \quad (1) $$

In this equation, $\rho$ represents the density of the medium, $u$ is the displacement field, $\nabla \cdot T$ is the divergence of the stress tensor, and $f$ is the excitation source. The elastic properties of the medium and its deformations are implicitly represented in the stress tensor $T$ due to a linear relationship between the stress tensor and the strain tensor $\epsilon$. This relationship is given by the generalized Hooke's law (Lucibella, 2014): $T = C : \epsilon$ where $\epsilon : \epsilon$ denotes the tensor product, $C$ is the elasticity tensor, and $\epsilon$ is the strain tensor. The elasticity tensor contains the elastic coefficients of the medium, which depend on geo-physical properties. According to the wave Eq. (1), the main parameters influencing the displacement field are the speeds of the compression wave (P-wave) and the shear wave (S-wave) involved in the stress-strain relationship, as well as the density. The parameters of the seismic source also play a significant role. The uncertainty associated with these parameters will be discussed in Section 3.

To solve the wave Eq. (1), we use the EFISPEC3D code,<fnref n="1" /> a three-dimensional parallel open-source computational program. The accuracy of the program has been verified with international benchmarks (Maufroy et al., 2015; Chaljub et al., 2015) and semi-analytical solutions (De Martin, 2011; Matsushima et al., 2014). EFISPEC3D employs a spectral finite element method to solve the weak formulation of Eq. (1) and utilizes an explicit time-stepping scheme constrained by a CFL (Courant-Friedrichs-Lewy) condition (De Moura and Kubrusly, 2013). Specifically, the three-dimensional study domain is modeled with a mesh size of 30 × 30 × 50 km, consisting of 700,000 hexahedral elements distributed over three layers to conform to the model proposed by Dorel et al. (1979) (see Fig. 5). Element sizes are adjusted across

the entire simulation domain to propagate a minimum wavelength of 260 m without numerical dispersion for Dorel layers and 62.5 m for the dam layer. The goal is to avoid excessively large element sizes relative to the desired wavelength, as this would result in underrepresentation of wave details. Conversely, overly small element sizes can significantly increase the mesh element count, thereby increasing computation time.

## 3. Surrogate model

The most common way to define a surrogate model is to refer to it as a model of the model. Surrogate models can also be referred to by other names such as emulator, approximator, simplified model, or metamodel (Misener and Biegler, 2023). In our study, the objective is to replace the solution given by the spectral method using the EFISPEC3D code for solving the elastodynamic wave equation. Additionally, we will use this surrogate model to calculate sensitivity indices, as the EFISPEC3D code is highly time-consuming, a single spectral-element simulation longs approximately 7 h on 48 cores. The construction of the surrogate model involves executing a limited number of simulations through an experimental design (or experience plan) that represent combinations of different parameters.

### 3.1. Experimental design: Uncertain parameters

In the context of uncertainty quantification, the experimental design consists of $N$ models with $N_u$ uncertain parameters. Based on the wave Eq. (1), the primary parameters influencing the displacement “u” are the velocities of the P and S waves involved in the stress-strain relationship, as well as the density. The characteristics of the seismic source (location, focal mechanism, etc.) are also influential parameters.

In this study, we assumed for simplicity that the focal mechanism of the seismic source and the magnitude are known, with fixed values to avoid a large number of uncertain parameters. Similarly, due to the empirical correlation between P-wave and S-wave velocities and density (Brocher, 2005), we chose to neglect the uncertainty related to P-wave velocities and density. Instead, we focus solely on the uncertainty of S-wave velocities ($V_s$) in the three-layer model of Dorel et al. (1979), the S-wave velocity of the dam, and the position of the seismic source. Our hypothesis is based on the fact that shear wave velocities ($Vs$) play a predominant role in the assessment of seismic effects, as they are generally the most critical for the studied infrastructures. Concerning the source time function, any temporal structure can be changed retroactively (by a convolution with the desired spectrum). These parameters are modeled with uniform distributions within specific intervals of confidence, as indicated in Table 1. This choice is based on the reason that, uniform distributions serve as a good initial approximation when detailed information about the exact parameter distributions is unavailable, as in the case of Guadeloupe. So, we chose not to complicate the analysis with other distributions such as Gaussian distribution that require more precise data. In total, seven parameters are chosen uncertain.

To create an effective surrogate model, it is crucial that the $N$ models with uncertain parameters are sufficiently well-distributed across the specified ranges shown in Table 1. In order to optimize the multi-dimensional distribution of the experimental designs, we use the Latin hypercube method (Stein, 1987).

### 3.2. Response surface

To construct the metamodels of time or frequency series, a dataset consisting of $N$ observations containing $N_u = 7$ uncertain parameters is generated together with the corresponding output variables of the model to be modeled by the surrogates. The $N_u$ uncertain parameters will be represented as a matrix $X \in M_{N,N_u}(R)$, generated as explained in the previous steps. The output variable to be explained will be

<footnote n="1">¹ Available at http://efispec.free.fr, co-development at https://gitlab.com/brgm/seismic-modeling/efispec3d.</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 5. 3D Mesh Containing the Dumanoir Dam. The mesh size in the third layer of Dorel (purple) ranges between 612 m and 890 m. In the second layer (yellow), the mesh size ranges between 248 m and 263 m, and in the first layer (green), it ranges between 132 m and 246 m. Within the dam (blue), the mesh size ranges between 14 m and 30 m.</figcaption></figure>

<table><caption>**Table 1**

Intervals of confidence of the uncertain parameters. Top: mean value. Bottom: interval of confidence. U: uniform probability distribution.</caption><thead><tr><th>V<sub>s,dam</sub></th><th>V<sub>s1</sub></th><th>V<sub>s2</sub></th><th>V<sub>s3</sub></th><th>X<sub>s</sub></th><th>Y<sub>s</sub></th><th>Z<sub>s</sub></th></tr></thead><tbody><tr><td>375 m/s</td><td>1400 m/s</td><td>2400 m/s</td><td>4000 m/s</td><td>14 568 m</td><td>18 053 m</td><td>-28020 m</td></tr><tr><td>U[250,450]</td><td>U[1150,1650]</td><td>U[2150,2650]</td><td>U[3750,4250]</td><td>U[14 068,15 068]</td><td>U[17 553,18 553]</td><td>U[-28 520,27 520]</td></tr></tbody></table>

represented as a matrix $Y \in M_{N,N_i}(R)$ where $N_i = 301$ represents the number of time steps.

$$X = \begin{bmatrix} x_1^{(1)} & x_2^{(1)} & \cdots & x_{N_u}^{(1)} \\ x_1^{(2)} & x_2^{(2)} & \cdots & x_{N_u}^{(2)} \\ \vdots & \vdots & \ddots & \vdots \\ x_1^{(N)} & x_2^{(N)} & \cdots & x_{N_u}^{(N)} \end{bmatrix} \qquad Y = \begin{bmatrix} y_1^{(1)} & y_2^{(1)} & \cdots & y_{N_t}^{(1)} \\ y_1^{(2)} & y_2^{(2)} & \cdots & y_{N_t}^{(2)} \\ \vdots & \vdots & \ddots & \vdots \\ y_1^{(N)} & y_2^{(N)} & \cdots & y_{N_t}^{(N)} \end{bmatrix}$$

The matrix $Y$ is generated by running the experimental design through the model $M$ where $M$ is computed by the spectral-element method code EFISPEC3D:

$$Y = M(X) \qquad (2)$$

in order to construct a response surface that contains the model response to each set of uncertain parameters. The Fig. 6(a) shows a response surface containing $N = 300$ simulations generated by the model $M$ in three directions (East-West, vertical, North-South) at the on-site station. Similarly, Fig. 6(b) depicts a response surface containing $N = 300$ simulations generated by the model $M$ in three directions (East-West, vertical, North-South) at the reference station. We have chosen to divide our analysis into three steps: the first step for time-domain signals, the second step for the amplitudes of the Fourier transform of these signals (see Figs. 7(a) and 7(b)), and the final step for the phases of the Fourier transform (see Figs. 8(a) and 8(b)). However, due to the complexity of processing phases, we will create a metamodel for the real and imaginary parts of the Fourier transform, which will then allow us to reproduce the phases.

### 3.3. Neural network

In the previous sections, we presented our initial model given by Eq. (2), where $Y$ represents the response surface, $X$ the experimental design containing various distributions of input parameters, and $M$ the numerical model computed by the spectral-element method code EFISPEC3D for solving the motion equation. The objective of this part is to create a substitution model called NN for the solution of the wave equation as follows:

$$Y \approx \mathcal{N}\mathcal{N}(X) + \epsilon \qquad (3)$$

where $\epsilon$ is the modeling error (difference between synthetic values generated by $M$ and values given by the substitution model). The creation of substitution models involves minimizing the modeling error.

There are numerous substitution models in the literature, including polynomial chaos (Baudin and Martinez, 2010), neural networks, and Gaussian processes (Fernique, 1971). In this article, we focus on using neural networks, which are often constructed from reasonably sized experimental designs.

#### 3.3.1. Multi-layer perceptron (MLP)

The neural networks implemented in this study is called "multi-layer perceptron (MLP)" (Popescu et al., 2009), consist of 6 layers: one input layer, four hidden layers, and one output layer (see Fig. 9). The input layer aggregates all uncertain input variables from our original model, which are then transmitted to neurons in the hidden layers where they are weighted by multiplicative weights $W_{ij}$. Finally, this value is modulated by an activation function before being transmitted to neurons in the output layer. There are several activation functions in the literature (Farhadi, 2017), among which the function used in this document is the "ReLU" function, represented by the following expression:

$$\mathrm{relu}(x) = \max(0, x).$$

The output layer is either the time-series or the frequency-series which are represented by discrete signals rather than continuous ones, it will be represented by 301 output neurons (see Fig. 9). In our context, the network inputs include the shear wave velocities $V_s$ in each layer of the Dorel model, the velocity in the dam layer, and the position of the seismic source, requiring seven neurons in this layer. It is worth noting that we have constructed nine neural networks, each with the same characteristics, to approximate all components (three directions for displacement, three for velocity, and three for acceleration). This article focuses solely on the velocity field.

### 3.4. Training and optimization

In the context of surrogate modeling, our objective is to train a neural $f_{NN}$ using the architecture and activation function defined in the previous section, to approximate the model $M(X)$, which can be represented by an unknown function $f: \mathbb{R}^d \to \mathbb{R}^{n_d}$ where $d = 7$ and $n_d = 301$. Specifically, we aim to find a set of network parameters $\Theta_{NN} := \{(W', b')'_{L+1} \text{ where } L \text{ represents the number of hidden layers such that } f_{NN}(X; \Theta_{NN}) \text{ approximates well } f(X) \text{ within a small tolerance. We performed a grid search to optimize the hyperparameters of the neural network for surrogate modeling. The search space included}$
<figure><img src="image_8.png" /></figure>

(a) A response surface containing 300 time-domain signals generated by the model $M$ depicting the variation of the three components of velocity over time at the on-site station

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

(b) A response surface containing 300 time-domain signals generated by the model $M$ depicting the variation of the three components of velocity over time at the reference station

<figure><img src="image_10.png" /><figcaption>Fig. 6. Time-domain signals.</figcaption></figure>

(a) Variation of the amplitude of the Fourier transform of the time-domain signals at the reference station with respect to frequency.

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

(b) Variation of the amplitude of the Fourier transform of the time-domain signals with respect to frequency at the site station.

<caption>Fig. 7. Fourier transform of the time-domain signals at the reference and the site station.</caption>

different values for the number of layers (2, 3, and 4), the number of neurons (100, 200, and 300), and the learning rate (0.01, 0.001, and 0.0001) (Jepkooch et al., 2021). This process allowed us to identify the optimal configuration for approximating the Fourier spectrum amplitude generated by model $M$, leading to the choice of 4 hidden layers for the site station and 3 for the reference station, based on the performance of the loss function. Based on Tables 2 and 3, we can conclude that the minimum value of the loss function at the last epoch is reached with 4 hidden layers at the site station and 3 hidden layers at the reference station. Here, abusing notation, we write $f_{NN}(X; \Theta_{NN})$ to emphasize the network's dependency on the parameter set $\Theta_{NN}$. A network formed in this way will then serve as a surrogate model for the EFISPEC3D code represented by the function $f$. To this end, we first collect the dataset of $N$ input-output pairs $(X^i, f(X^i))_{i=1}^{300}$ represented by the “Experimental Design-Response Surface” set generated

in Sections 3.1 and 3.2, and then formulate an optimization problem

$$ \arg \min_{\Theta_{NN}} \frac{1}{N} \sum_{i=1}^{N} L(f_{NN}(X^{(i)}; \Theta_{NN}), f(X^{(i)})) $$

where $L(a, b)$ is a cost (or loss) function measuring the distance between vectors $a$ and $b$. In this article, the quadratic cost function $C(a, b) = \|a-b\|^2$ is used to quantify the difference between predictions made by the function $f_{NN}$ and the EFISPEC3D-generated data to be approximated by the function $f$. The optimization problem can be solved using a gradient-based method, such as stochastic gradient descent (Driancourt, 1994) or Adam (Kingma and Ba, 2014) (in this article, we used the Adam algorithm, which has shown more success in machine learning than other optimization algorithms Pan and Li, 2023). In these methods, the gradient of the cost function with respect to the network parameters is typically computed using the chain rule and a technique
<figure><img src="image_12.png" /></figure>

(a) Variation of the amplitude of the Fourier transform of the time-domain signals at the reference station with respect to frequency.

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

(b) Variation of the amplitude of the Fourier transform of the time-domain signals with respect to frequency at the site station.

<figure><img src="image_14.png" /><figcaption>Fig. 9. Neural network architecture composed of one input layer, 4 hidden layers, and one output layer. The output layer is either the time-series or the frequency-series. The input layer contains the parameters of the wave propagation model shown in Table 1.</figcaption></figure>

called backpropagation (Boutin, 1994). For further details, refer to the review article (Bottou et al., 2018).

The performance of a substitute neural network, which is the accuracy of the approximation $f(X) \approx f_{NN}(X; \Theta_{NN})$, depend on two factors:

* The network structure, which includes the number of layers and neurons, the activation functions used, as well as the optimization/regularization techniques and parameters used to train the network,

* The choice of the dataset $\{(X^{(i)}, f(X^{(i)}))\}_{i=1}^N$

Currently, there is no theory addressing the precise dependence of approximation accuracy on the aforementioned factors (Sharma et al., 2017; Gomes and Ludermir, 2013). Therefore, it is generally recommended to explore and test multiple network architectures, different learning parameter configurations, and various datasets to determine the setup that provides the most accurate approximation for the function $f$. In Fig. 10, we present three subfigures showing the variation of the loss function with epochs for different configurations of hidden neurons and learning rates. On the left, for a neural network with two hidden layers; in the middle, for a neural network with three hidden layers; and on the right, for a neural network with four hidden layers. For all three graphs, we can observe that the Loss function converges for the learning rates [0.0001, 0.001], which is not the case for a learning rate of 0.01. To find the minimum among these three graphs, we present in Table 2 the value of the loss function at the last epoch for the three graphs. We observe that the global minimum is achieved

<caption>Table 2</caption>

Value of the loss function at the last epoch for the minimum of the three graphs in Fig. 10.

<table><caption>Table 3</caption><thead><tr><td>Nb of hidden layer</td><td>Nb of hidden neurons</td><td>leaning rate</td><td>Loss value at last epoch</td></tr></thead><tbody><tr><td>2</td><td>200</td><td>0.001</td><td>9.23e-5</td></tr><tr><td>3</td><td>300</td><td>0.001</td><td>8.49e-5</td></tr><tr><td>4</td><td>100</td><td>0.001</td><td>8.44e-5</td></tr></tbody></table>

Value of the loss function at the last epoch for the minimum of the three graphs in Fig. 11.

<table><thead><tr><td>Nb of hidden layer</td><td>Nb of hidden neurons</td><td>leaning rate</td><td>Loss value at last epoch</td></tr></thead><tbody><tr><td>2</td><td>300</td><td>0.001</td><td>0.00125</td></tr><tr><td>3</td><td>300</td><td>0.001</td><td>0.00064</td></tr><tr><td>4</td><td>300</td><td>0.001</td><td>0.00081</td></tr></tbody></table>

with three hidden layers, 300 neurons, and a learning rate of 0.001. This is done to construct a neural network to reproduce the Fourier transform amplitudes for the reference station. We conducted the same study for the site station, see Fig. 11 and Table 3. The results show that the minimum is achieved with four hidden layers, 100 neurons, and a learning rate of 0.001.

After setting up our “input-output” database, we proceeded to divide this dataset into two distinct parts for training and validating the network. 90% of the data was reserved for training the networks, totaling 270 simulations. As mentioned earlier, these simulations were
<figure><img src="image_15.png" /><figcaption>Fig. 10. Variation of the loss function with respect to epochs for different configurations of hidden neurons and learning rates (Site station). on the left for a neural network with two hidden layers; in the middle, for a neural network with three hidden layers; and on the right, for a neural network with four hidden layers.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 11. Variation of the loss function with respect to epochs for different configurations of hidden neurons and learning rates (reference station). on the left for a neural network with two hidden layers; in the middle, for a neural network with three hidden layers; and on the right, for a neural network with four hidden layers.</figcaption></figure>

used to adjust the network weights using the Adam optimization algorithm to minimize a loss function L. The remaining 30 simulations, constituting 10% of the data, were used to validate our networks. Unlike the training stage where target values were provided to the network for learning to reproduce expected outputs, during validation, only input vectors were fed into the networks. Thus, the network was evaluated on scenarios it had not encountered during training, testing its ability to generalize and predict new geological configurations.

### 3.5. Validation

In this section, we evaluated the performance of neural networks as a substitute model for seismic wave modeling compared to the model M. The objective was to demonstrate that neural networks can reproduce seismic signals while significantly reducing computation time. To assess the convergence error, we used the mean squared error (MSE). Fig. 12 shows the variation of MSE across epochs on the validation dataset. The figure illustrates how the error between Time-domain signals generated by M (target outputs) and by the neural networks predictions decreases during training. To evaluate the robustness of the surrogate model, we implemented a 5-fold cross-validation procedure (Lyu et al., 2022). This approach ensures that the model is trained and tested on different subsets of the dataset,

providing a comprehensive assessment of its performance. The results demonstrate that the model achieves a high coefficient of determination ($R^2$) of 0.88, indicating its strong ability to accurately reproduce the time-frequency signals generated by the wave equation model across all folds. To visualize the discrepancy between predictions made by the neural networks and those from the model M on a test simulation, we plot figures of the three-dimensional velocity field (for time-domain signals, see Figs. 13, 14, and for amplitudes, see Figs. 15, 16). The discrepancy between results (for Fourier transform amplitudes at both stations) is negligible, demonstrating the reliability of the neural networks in replacing M. However, concerning time-domain signals, we can conclude that our metamodel is capable of accurately reproducing those generated at the reference station but performs poorly for those recorded at the site station. Regarding the phases, as previously mentioned, we created a metamodel to reproduce the real part first and the imaginary part for the reference station (Figs. 17(a), 17(b)). The results show that the neural network is capable of reproducing the real and imaginary parts for the reference station. Similarly, for the site station, the results in Figs. 18(a), 18(b) show good performance of the metamodel in all three directions at low frequency, with less accuracy at high frequency. These results allow us to construct a phase metamodel (show Figs. 19(a), 19(b)) using the following formula:

$$ \phi = \arctan \left( \frac{\Im}{\Re} \right) $$
<figure><img src="image_17.png" /><figcaption>Fig. 12. Variation of the loss function for time-domain signals generated at the reference station in each direction over epochs.</figcaption></figure>

where:

* $\phi$ is the phase

* $\mathbb{R}$ is the real part of the signal

* $\mathfrak{F}$ is the imaginary part of the signal.

Figs. 20(a) and 20(b) shows the overall error between predictions and target outputs across the entire test set (for amplitude of Fourier spectrum at the site and reference station).

In addition to the evaluation of the model's performance using standard metrics, we applied the bootstrap technique to assess the stability of the surrogate model. Bootstrap sampling involves generating multiple new datasets by randomly resampling the original dataset with replacement, which allows us to estimate the variability in the model's predictions. By evaluating the model's performance on these bootstrap samples, we gained insight into the robustness of the surrogate model, as evidenced by the low standard deviation of the mean, represented in Figs. 21(a), 21(b). To assess the robustness of our surrogate model in handling rare and extreme seismic events, we focused on identifying the test simulations most distant from the median signal using the mean squared error (MSE). These simulations were selected as they represent atypical, extreme seismic scenarios. We then compared the model's prediction to the corresponding SEM signal generated by the wave equation model. The results showed that the model achieves a coefficient of determination ($R_2$) of 0.7, demonstrating its ability to accurately reproduce the rare signals. All the techniques previously employed have consistently validated the ability of our surrogate model to accurately reproduce both the time-series and frequency-domain signals generated by the wave equation model. These results confirm the robustness and reliability of the metamodel in capturing the essential features of the physical phenomena, ensuring its effectiveness for the simulation of seismic events under uncertainty.

<table><caption>**Table 4**
Comparison of performance metrics for different machine learning metamodels applied to seismic signal prediction.</caption><thead><tr><td>Metamodel</td><td>Mean absolute error</td><td>Regression coefficient ($R_2$)</td><td>Mean square error</td></tr></thead><tbody><tr><td>Convolutional Neural Network (CNN)</td><td>1.81 × 10<sup>-9</sup></td><td>0.9</td><td>2.6 × 10<sup>-9</sup></td></tr><tr><td>Random Forest</td><td>5.4 × 10<sup>-9</sup></td><td>0.82</td><td>5.4 × 10<sup>-9</sup></td></tr><tr><td>MLP</td><td>0.94</td><td>0.99</td><td>1.36 × 10<sup>-9</sup></td></tr></tbody></table>

Exploration of other neural network architectures and machine learning methods

In this part, we explore alternative approaches to surrogate modeling, including Convolutional Neural Networks (CNNs) and Random Forests. These models are widely applied across various applications, including image and signal processing, classification, and regression tasks (Li et al., 2021; Wang et al., 2021). We compare their performance to our MLP approach in terms of their ability to reproduce the time signals generated by the wave equation model based on input parameters. The performance of each model is evaluated based on several metrics, such as mean absolute error (MAE), to determine which provides the best results for predicting seismic signals as shown in Table 4.

Based on the results presented in the table, we can conclude that the Multi-Layer Perceptron (MLP) is the most effective surrogate model. It demonstrates superior performance in reproducing the time signals generated by the wave equation model, achieving the lowest Mean Absolute Error and the highest Coefficient of Determination ($R_2$). This indicates its strong ability to capture the complex relationships between input parameters and the corresponding seismic responses.
<figure><img src="image_18.png" /><figcaption>**Fig. 13.** Time-domain signals generated by $\mathcal{M}$ and by the neural networks for a test simulation at the reference station. The input vector used to generate this time series is $X = (V_{x,dem} = 284.9 \text{ m/s}, V_{y,1} = 1157.9 \text{ m/s}, V_{y,2} = 2634 \text{ m/s}, V_{y,3} = 3829.3 \text{ m/s}, X_{\text{source}} = 14242.6 \text{ m}, Y_{\text{source}} = 18076.8 \text{ m}, Z_{\text{source}} = -28091.4 \text{ m})$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 14.** Time-domain signals generated by $\mathcal{M}$ and by the neural networks for a test simulation at the site station. The input vector used to generate this time series is $X = (V_{x,dem} = 284.9 \text{ m/s}, V_{y,1} = 1157.9 \text{ m/s}, V_{y,2} = 2634 \text{ m/s}, V_{y,3} = 3829.3 \text{ m/s}, X_{\text{source}} = 14242.6 \text{ m}, Y_{\text{source}} = 18076.8 \text{ m}, Z_{\text{source}} = -28091.4 \text{ m})$.</figcaption></figure>
<figure><img src="image_20.png" /></figure>

(a) Amplitude of Fourier spectrum at the reference station for $X=(V_{s_{dam}} = 284.9m/s, V_{s_1} = 1157.9m/s, V_{s_2} = 2634m/s, V_{s_3} = 3829.3m/s, X_{source} = 14242.6m, Y_{source} = 18076.8m, Z_{source} = -28091.4m$. Blue curve: model $M$, Orange curve: neural network surrogate.

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

(b) Amplitude of Fourier spectrum at the reference station for $X=(V_{s_{dam}} = 323m/s, V_{s_1} = 1586.5m/s, V_{s_2} = 2356m/s, V_{s_3} = 4123m/s, X_{source} = 14750.5m, Y_{source} = 18346.6m, Z_{source} = -27631.1m$. Blue curve: model $M$, Orange curve: neural network surrogate.

<caption>Fig. 15. Amplitude of Fourier spectrum generated by $M$ and those predicted by the neural networks for 2 test simulation at the reference station. The x-axis is limited to 4 Hz, corresponding to the maximum valid frequency range for the simulations.</caption>

## 4. Global sensitivity analysis

When using our substitution model, sensitivity analysis (SA) methods are valuable tools. Several input parameters of this model may be subject to uncertainties that can affect simulation results. Sensitivity analysis allows us to understand the influence of these uncertain parameters and quantify their impact on the model response.

### 4.1. Sensitivity analysis methods

There are several sensitivity analysis methods in the literature (Iooss, 2011), and their selection often depends on the model complexity and the objectives of the analysis. In this article, we distinguish between two types of methods: screening methods and importance measures.

#### 4.1.1. Screening method

Screening can be a useful approach to assess and manage uncertainties. Once the sources of uncertainty are identified, screening helps prioritize or classify them based on their importance to the final outcome. Some uncertainties may be deemed less critical, while others may be considered high-uncertainty factors requiring particular attention. Among the various screening techniques, group screening is notable (De Zanet and Viquerat, 2022). It involves creating several groups of inputs and identifying the most influential ones. By gradually repeating this process while retaining the influential groups, one eventually extracts the most influential inputs. This technique assumes no interaction between inputs and that the output variation is monotonic with respect to each input. The objective of screening is to identify inputs that have little influence on the output and fix them at their nominal values, thereby reducing the number of uncertain parameters by eliminating those with negligible influences.

#### 4.1.2. Importance measures

These methods assume the ability to obtain N Samples of input-output pairs from simulations $(X, Y)$ of sufficient size where

$$X = \begin{bmatrix} x_1^{(1)} & x_2^{(1)} & \dots & x_d^{(1)} \\ x_1^{(2)} & x_2^{(2)} & \dots & x_d^{(2)} \\ \vdots & \vdots & \ddots & \vdots \\ x_1^{(N)} & x_2^{(N)} & \dots & x_d^{(N)} \end{bmatrix}, \quad Y = \begin{bmatrix} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(N)} \end{bmatrix} = \begin{bmatrix} K(X^{(1)}) \\ K(X^{(2)}) \\ \vdots \\ K(X^{(N)}) \end{bmatrix} \quad (4)$$

$K$ is a function that represents the metamodel created using the neural network. The $y^{(i)}$ correspond to the responses of the metamodel to the inputs $X^{(i)}$. Once this plan is obtained, the importance measures allow for a precise ranking of the influence on the output $y^{(i)}$ of all inputs $x_j^{(i)}$, unlike screening techniques which aim to detect non-influential inputs.

In the general context of a nonlinear and non-monotonic model, such as in our case, we can estimate the importance of inputs on the model output using the functional variance decomposition technique (Touzani, 2011).

#### 4.1.3. Functional variance decomposition

Consider a mathematical model $f$ associated with $d$ independent random parameters $x = (x_1, ..., x_d)$ leading to a random response $f(x)$. (Sobol', 1990) demonstrated that if $f(x)$ belongs to the space of square integrable functions, $L_2$, then it can be expanded as follows:

$$f(x) = f(0) + \sum_{i=1}^{d} f_i(x_i) + \sum_{i<j}^{d} f_{ij}(x_i, x_j) + \dots + f_{1,...,d}(x_1, ..., x_d) \quad (5)$$

where $f(0)$ represents the mean value of $f(x)$ and $f_i, f_{ij}, ...$ are orthogonal functions. This representation is generally referred to as variance
<figure><img src="image_22.png" /></figure>

(a) Amplitude of Fourier spectrum at the site station for $X=(V_{s_{dam}} = 284.9m/s, V_{s_1} = 1157.9m/s, V_{s_2} = 2634m/s, V_{s_3} = 3829.3m/s, X_{source} = 14242.6m, Y_{source} = 18076.8m, Z_{source} = -28091.4m$. Blue curve: model $M$, Orange curve: neural network surrogate.

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

(b) Amplitude of Fourier spectrum at the site station for $X=(V_{s_{dam}} = 284.9m/s, V_{s_1} = 1157.9m/s, V_{s_2} = 2634m/s, V_{s_3} = 3829.3m/s, X_{source} = 14242.6m, Y_{source} = 18076.8m, Z_{source} = -28091.4m$. Blue curve: model $M$, Orange curve: neural network surrogate.

<caption>Fig. 16. Amplitude of Fourier spectrum generated by the physical model $M$ and those predicted by the neural networks for two validation test at the site station. The x-axis is limited to 4 Hz, corresponding to the maximum valid frequency range of the simulations.</caption>

decomposition (ANOVA). St et al. (1989). In fact, by squaring and taking the expectation, we obtain the variance decomposition of $f(x)$:

$$V = \sum_{i=1}^{d} V_i + \sum_{i<j}^{d} V_{ij} + \dots + V_{1,...,d} \quad (6)$$

where $V$ represents the total variance of $f(x)$, $V_i$ is the partial variance due to the individual input $x_i$ and $V_{(i_1,...,i_s)}$ is the partial variance due to interactions among inputs belonging to the subset of parameters $(x_{(i_1)},...,x_{(i_s)})$. Importance measures are then defined as follows:

$$S_{i_1,...,i_s} = \frac{V_{i_1,...,i_s}}{V} \quad (7)$$

These coefficients are called “variance-based importance measures” or simply “Sobol indices”. They quantify the contribution of an input variable to the output variance, either individually through first-order indices $S_i$, or through interactions with other input variables via higher-order indices. $S_{ijk}$, $S_{ijkk}$, ...

## 4.2. Method of calculation

We begin by calculating the first-order indices $S_i$ which measures the contribution of $x_i$ solely on the variance of $y$, without interaction with the other inputs. The calculation of $S_i$ is based on the following formula:

$$S_i = \frac{V_i}{V(y)} = \frac{V_{x_i}(E_{x_{-i}}(y|x_i))}{V(y)} \quad (8)$$

where $V(y)$ represents the variance of $y$ and $V_i$ represents the variance due to $x_i$ alone, $x_{-i}$ represents all parameters other than $x_i$ and $E(y|x_i)$ is the conditional expectation representing the average of $y$ when $x_i$ is fixed.

According to Eq. (8), calculating $S_i$ analytically is complex because, $E_{x_{-i}}(y|x_i)$ involves multidimensional integration. Directly computing $S_i$ via Monte Carlo simulation (SMC) is also computationally expensive (Saltelli et al., 2008). The numerator in Eq. (8) involves a double SMC loop: the inner loop $E_{x_{-i}}(y|x_i)$ computes the average of $y$ using $n_1$ random samples of $x_i$; and the outer loop computes $V_{x_i}(E_{x_{-i}}(y|x_i))$ by iterating the inner loop $n_2$ times over different values of $x_i$. Additionally, $n_3$ more SMC iterations are needed to calculate $V(y)$. The computational cost of this sample-based method, defined as the total number of model evaluations to compute all $S_i$ ($i = 1,...,7$), is $7n^2 + n$ if $n_1 = n_2 = n_3 = n$. This cost scales with $n$ and $k = 7$ (number of uncertain parameters), and becomes prohibitive if a single model evaluation is costly in terms of computation time, since $n$ is often much larger than 1000 in many practical applications. Various algorithms have been proposed to reduce the computational cost of Sobol' indices calculations. These algorithms can be categorized into analytical methods and sample-based methods. In analytical methods, the original model $y = f(x)$ is typically approximated by a special form surrogate model, allowing the multidimensional integral to be converted into multiple univariate integrals that can be easily computed analytically or numerically. Zhang et al. (2015) approximated the original model $y = f(x)$ by a product of univariate functions; subsequently, the univariate integral was computed using Gaussian quadrature. Sudret (2008) proposed that if the original model is approximated by a polynomial chaos expansion (PCE), the Sobol' index can be calculated by post-processing the PCE coefficients. In this study, we detail two sample-based methods: the first is based on quasi-Monte Carlo sampling to estimate all indices using the Sobol method, and another method proposed by Li and Mahadevan to estimate first-order Sobol' indices to confirm our estimations.
<figure><img src="image_24.png" /></figure>

(a) Real part of Fourier spectrum at the reference station for $X=(V_{s_{dam}} = 284.9m/s, V_{s_1} = 1157.9m/s, V_{s_2} = 2634m/s, V_{s_3} = 3829.3m/s, X_{source} = 14242.6m, Y_{source} = 18076.8m, Z_{source} = -28091.4m$. Blue curve: model $M$, Orange curve: neural network surrogate.

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

(b) Imaginary part of Fourier spectrum at the reference station for $X=(V_{s_{dam}} = 284.9m/s, V_{s_1} = 1157.9m/s, V_{s_2} = 2634m/s, V_{s_3} = 3829.3m/s, X_{source} = 14242.6m, Y_{source} = 18076.8m, Z_{source} = -28091.4m$. Blue curve: model $M$, Orange curve: neural network surrogate.

<caption>Fig. 17. Real and imaginary part of Fourier spectrum generated by the physical model $M$ and those predicted by the neural networks for one validation test at the reference station. The x-axis is limited to 4 Hz, corresponding to the maximum valid frequency range of the simulations.</caption>

### 4.3. Protocol for calculation and evaluation of sobol indices

First, we will detail the two methods for calculating Sobol indices: the Li and Mahadevan method, and the Sobol method. The Li and Mahadevan method focuses on a more analytical approach, whereas the Sobol method uses a functional decomposition to evaluate the variance. Next, we will present the results obtained with the first method and compare them with those obtained by the second method. We will also calculate the total Sobol indices using these two approaches to determine the overall influence of uncertain parameters on ground speed. Subsequently, we will discuss the number of simulations required to evaluate the stability of each method, specifying the exact number of simulations needed. Finally, we will introduce the Bootstrap technique to assess the stability of the Sobol indices estimation. This technique will allow us to calculate the mean and variance of the indices across several generated Bootstrap samples, thus providing a robust evaluation of the stability of the estimates.

#### 4.3.1. Li and Mahadevan method

In this section, we present a method proposed by Li and Mahadevan (2016) to estimate the first-order Sobol index using only input-output samples from a model, without requiring additional model evaluations. The proposed method includes two algorithms for calculating the indices. In this paper, we will use the first algorithm of this method to calculate the Sobol indices.

In Section 3.3, the EFISPEC3D code is approximated by the neural function $f_{NN}(x; \Theta_{NN})$ where $\Theta_{NN}$ are parameters calculated by optimization. For simplicity, let us denote $y = f_{NN}(x; \Theta_{NN}) = f(x)$ where $x = \{x_1, ..., x_7\}$ is the vector containing the uncertain parameters.

Let us take up Eq. (8) again. The first-order Sobol index is represented by the expression,

$$S_i = \frac{V_{x_i}(E_{x_{-i}}(y|x_i))}{V(y)} \qquad (9)$$

We divide $x$ into two sets: ($x = x_i, x_{-i}$) where $x_{-i}$ are the inputs other than $x_i$. $E_{x_{-i}}(y|x_i)$ is a function of $x_i$, and it can be shown that the mean value of $y$ over ($x = x_i, x_{-i}$) is equal to the mean of $E_{x_{-i}}(y|x_i)$ over $x_i$, that is $E_{x_i}(E_{x_{-i}}(y|x_i)) = E(y)$, (see Weiss et al., 2006). Assuming that $x_i \in \varphi$ where $\varphi$ is a bounded interval, based on the extreme value theorem and since $E_{x_{-i}}(y|x_i)$ is a continuous function of $x_i$, it must have a maximum and a minimum value over $\varphi$. Its mean value, that is $E_{\varphi}(E_{x_{-i}}(y|x_i))$ for $x_i \in \varphi$ must be between the minimum and the maximum, thus

$$\min_{x_i \in \varphi} E_{x_{-i}}(y|x_i) \le E_{\varphi}(E_{x_{-i}}(y|x_i)) \le \max_{x_i \in \varphi} E_{x_{-i}}(y|x_i)$$

and using the intermediate value theorem, we can assume that there exists a $x_i^* \in \varphi$ such that

$$E_{\varphi}(E_{x_{-i}}(y|x_i)) = E_{x_{-i}}(y|x_i^*) = E_{\varphi}(y). \qquad (10)$$

The interval $\varphi$ is generated using stratified sampling, which divides the sampling space of $x_i$ into $L$ equally probable intervals $\varphi^j$ ($j = 1, ..., L$). Each $x_i^j$ thus represents a random sample selected within one of these intervals $\varphi^j$. Reconsidering the inner loop, $E_{x_{-i}}(y|x_i)$ and for $\varphi = \varphi^j$, Eq. (10) shows that there exists an $x_i^j \in \varphi^j$ such that

$$E_{x_{-i}}(y|x_i^*) = E_{\varphi}(y). \qquad (11)$$

The outer loop $V_{x_i}(*)$ requires fixing $x_i$ at different positions, and these selected positions are samples from the distribution of $x_i$. Based on stratified sampling, the set of these unknown but existing points $x_i^* =$
<figure><img src="image_26.png" /></figure>

(a) Real part of Fourier spectrum at the site station for $X=(V_{sdam} = 284.9m/s, V_{s1} = 1157.9m/s, V_{s2} = 2634m/s, V_{s3} = 3829.3m/s, X_{source} = 14242.6m, Y_{source} = 18076.8m, Z_{source} = -28091.4m$. Blue curve: model $M$, Orange curve: neural network surrogate.

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

(b) Imaginary part of Fourier spectrum at the site station for $X=(V_{sdam} = 284.9m/s, V_{s1} = 1157.9m/s, V_{s2} = 2634m/s, V_{s3} = 3829.3m/s, X_{source} = 14242.6m, Y_{source} = 18076.8m, Z_{source} = -28091.4m$. Blue curve: model $M$, Orange curve: neural network surrogate.

<caption>Fig. 18. Real and imaginary part of Fourier spectrum generated by the physical model $M$ and those predicted by the neural networks for one validation test at the site station. The x-axis is limited to 4 Hz, corresponding to the maximum valid frequency range of the simulations.</caption>

$\{x_i^*, ..., x_i^{L*}\}$ within the equally probable intervals $\varphi = \{\varphi^{1*}, ... \varphi^{L*}\}$ can represent the distribution of $x_i$.

Finally, the first-order Sobol index is approximated by

$$S_i \approx \frac{\mathbf{V}_{\varphi}(\mathbf{E}_{\varphi}^j(y))}{\mathbf{V}(y)} \qquad (12)$$

### 4.3.2. Sobol method

To estimate sensitivity indices using the Sobol method (Jacques, 2011), consider a sample of size $N$, denoted by $X_{(N)} = (x_{k1}, ..., x_{k7})_{k=1,...,N}$ representing $N$ realizations of the input variables $(X_1, ..., X_7)$. The expectation of $Y$, $E[Y] = f_0$ and the variance of $Y$, $V[Y] = V$ are estimated by:

$$\hat{f}_0 = \frac{1}{N} \sum_{k=1}^{N} f(x_{k1}, ..., x_{k7}), \quad \text{et} \quad \hat{V} = \frac{1}{N} \sum_{k=1}^{N} f^2(x_{k1}, ..., x_{k7}) - \hat{f}_0^2. \qquad (13)$$

Estimating Sobol indices requires estimating the variance of the conditional expectation represented by the numerator in Eq. (8). This variance $V_i$ can be expressed as:

$$V_i = V(E[Y|X_i]) = E[E[Y|X_i]^2] - E[E[Y|X_i]]^2 = U_i - E[Y]^2, \qquad (14)$$

Sobol (1993) proposes to estimate the quantity $U_i$ given in Eq. (14), which is the expectation of the square of the conditional expectation of $Y$ given $X_i$, as a standard expectation. This involves conditioning on $X_i$ by varying all variables except $X_i$ between two calls to the function $f$. This requires two sets of input realizations, denoted $X_{(N)}^1 = (x_{k1}^1, ..., x_{k7}^1)_{k=1,...,N}$ and $X_{(N)}^2 = (x_{k1}^2, ..., x_{k7}^2)_{k=1,...,N}$:

$$\hat{U}_i = \frac{1}{N} \sum_{k=1}^{N} f(x_{k1}^1, ..., x_{k(i-1)}^1, x_{k1}^1, x_{k(i+1)}^1, ..., x_{k7}^1) \qquad (15)$$

$$\times f(x_{k1}^2, ..., x_{k(i-1)}^2, x_{k1}^1, x_{k(i+1)}^2, ..., x_{k7}^2). \qquad (16)$$

The first-order Sobol indices are estimated by:

$$\hat{S}_i = \frac{\hat{V}_i}{\hat{V}} = \frac{\hat{U}_i - \hat{f}_0^2}{\hat{V}}. \qquad (17)$$

For second-order indices, based on Eq. (7), they are represented by:

$$S_{ij} = \frac{V_{ij}}{V} \qquad (18)$$

In a similar manner to first-order indices, the variance $V_{ij}$ can be expressed as:

$$V_{ij} = V(E[Y|X_i, X_j]) - V_i - V_j = U_{ij} - E[Y]^2 - V_i - V_j \qquad (19)$$

and so, the estimation of second-order Sobol indices amounts to estimating the quantity $U_{ij} = E[E[Y|X_i, X_j]^2|]$. Sobol proposes to estimate the quantity $U_{ij}$, which is the expectation of the square of the conditional expectation of $Y$ given $X_i$ and $X_j$, as a standard expectation. This is done by conditioning on $X_i$ and $X_j$ while varying all variables except $X_i$ and $X_j$ in alternating calls to the function $f$. Therefore, the quantity $U_{ij}$ is estimated by:

$$\hat{U}_{ij} = \frac{1}{N} \sum_{k=1}^{N} f(x_{k1}^1, ..., x_{k(i-1)}^1, x_{k1}^1, x_{k(i+1)}^1, ..., x_{k(j-1)}^1, x_{k1}^1, x_{k(j+1)}^1, ..., x_{k7}^1) \times f(x_{k1}^2, ..., x_{k(i-1)}^2, x_{k1}^1, x_{k(i+1)}^2, ..., x_{k(j-1)}^2, x_{k1}^1, x_{k(j+1)}^2, ..., x_{k7}^2)$$

And so on, to compute all higher-order indices. The number of sensitivity indices constructed in this way, ranging from first to $p$th order (herep = 7, the number of uncertain parameters), is equal to $2^p - 1$. When the number $p$ of input variables is large, the number of sensitivity indices increases significantly. Estimating and interpreting all
<figure><img src="image_28.png" /></figure>

(a) Phase of Fourier spectrum generated by the physical model M and those predicted by the neural networks for one validation test at the site station for $X=(V_{s_{dam}} = 284.9m/s, V_{s_1} = 1157.9m/s, V_{s_2} = 2634m/s, V_{s_3} = 3829.3m/s, X_{source} = 14242.6m, Y_{source} = 18076.8m, Z_{source} = -28091.4m$. Blue curve: model M, Orange curve: neural network surrogate.

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

(b) Phase of Fourier spectrum generated by the physical model M and those predicted by the neural networks for one validation test at the reference station for $X=(V_{s_{dam}} = 284.9m/s, V_{s_1} = 1157.9m/s, V_{s_2} = 2634m/s, V_{s_3} = 3829.3m/s, X_{source} = 14242.6m, Y_{source} = 18076.8m, Z_{source} = -28091.4m$. Blue curve: model M, Orange curve: neural network surrogate.

<caption>**Fig. 19.** Phase of Fourier spectrum generated by the physical model M and those predicted by the neural networks for one validation test at the reference station. The x-axis is limited to 4 Hz, corresponding to the maximum valid frequency range of the simulations.</caption>

these indices quickly becomes impractical. Homma and Saltelli (1996) introduced total sensitivity indices to address this issue. These indices express the total sensitivity of the variance Y to a variable, capturing its sensitivity in all forms (both its direct effect and its interactions with other variables).

The total sensitivity index $S_{T_i}$ for the variable $X_i$ is defined as the sum of all sensitivity indices related to $X_i$:

$$S_{T_i} = \sum_{k \neq i} S_k$$

Here, $k \neq i$ denotes all subsets of indices that include the index i.

For instance, in a model with three input variables, $S_{T1} = S_1 + S_{12} + S_{13} + S_{123}$. Moreover, these total sensitivity indices can be estimated directly without calculating all higher-order indices. Indeed, it is straightforward to observe that the total sensitivity index can be expressed

$$S_{T_i} = 1 - \frac{V(E[Y|X_{\sim i}])}{V(Y)} = 1 - \frac{V_{\sim i}}{V}$$

where $V_{\sim i}$ represents the variance of the conditional expectation of Y given all variables except $X_i$, $V_{\sim i}$ is estimated similarly to $V_i$, but with a distinction: instead of varying all variables except $X_i$, we only vary $X_i$.

To estimate $V_{\sim i} = E[E[Y|X_{\sim i}]^2] - E[E[Y|X_{\sim i}]]^2 = U_{\sim i} - E[Y]^2$, Sobol estimates $U_{\sim i}$ as follows:

$$\hat{U}_{\sim i} = \frac{1}{N} \sum_{k=1}^{N} f(x_{(1)k1}, \dots, x_{(1)k(i-1)}, x_{(1)ki}, x_{(1)k(i+1)}, \dots, x_{(1)kp}) \times$$

$$f(x_{(1)k1}, \dots, x_{(1)k(i-1)}, x_{(2)ki}, x_{(1)k(i+1)}, \dots, x_{(1)kp})$$

And thus, we obtain $S_{\hat{f}_i} = 1 - \frac{\hat{U}_{\sim i} - \hat{f}_0^2}{V}$.

### 4.3.3. Results

In this part, first, we will compare the results obtained without and with the use of the metamodel to assess the effect of the number of simulations on the calculation of the indices. Next, we will present the results of the two methods (Li and Mahadevan, and Sobol) applied to the two stations and in the three directions (East-West, North-South, Vertical). Finally, we will discuss the results of the total indices and the application of the Bootstrap technique to evaluate the stability of these estimates.

#### Results using the Li and Mahadevan method

In Fig. 22(a), we present the results of the global sensitivity analysis using the Li and Mahadevan method without utilizing the neural network. Based on Eq. (6), we can see that the Sobol indices are values between 0 and 1; however, Fig. 22(a) shows that there are negative indices and indices greater than 1, which indicates that the number of simulations (300 simulations) is insufficient to calculate the Sobol indices. In Fig. 22(b), we present the results of the global sensitivity analysis using the Li and Mahadevan method and utilizing the neural network after generating 100,000 simulations, which shows a sum of indices less than 1.
<figure><img src="image_30.png" /><figcaption>(a) Error between neural network surrogate predictions and the model $M$ predictions at the reference station. Blue curve: model $M$, Orange curve: neural network surrogate. Red curve: Error</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(b) Error between neural network surrogate predictions and the model $M$ predictions at the site station (Sensor 1). Blue curve: model $M$, Orange curve: neural network surrogate. Red curve: Error</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(a) Bootstrap analysis of the surrogate model at the site station.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(b) Bootstrap analysis of the surrogate model at the site station.</figcaption></figure>

<caption>**Fig. 21.** Bootstrap analysis of the surrogate model at the site and reference station. The low standard deviation of the mean demonstrates the stability and robustness of the model across bootstrap samples.</caption>

### Results using the Li and Mahadevan method on the Fourier transform amplitudes

After generating 100,000 simulations using the neural network, we performed global sensitivity analysis in the three directions (East–West, North–South, vertical) on the Fourier transform amplitudes. This was

done using Algorithm 1 of the Li and Mahadevan method, dividing the distribution space of each uncertain parameter into $L = 10$ equally probable intervals. We also conducted tests with an increased number of sub-intervals, specifically 20, to assess the impact on the Sobol indices. However, we found that while this increased the computational
<figure><img src="image_34.png" /></figure>

(a) Influence of each uncertain parameter on the Fourier transformation of the East-West component of velocity at the site station (Sensor 1) without using neural network.

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

(b) Influence of each uncertain parameter on the Fourier transformation of the East-West component of velocity at the site station (Sensor 1) using neural network.

<figure><img src="image_36.png" /><figcaption>Fig. 23. Influence of each uncertain parameter on the three components of the Fourier transform amplitudes at the site station.</figcaption></figure>

time, it did not yield any substantial differences in the results when compared to those obtained using 10 sub-intervals. This suggests that 10 sub-intervals provide an adequate balance between computational efficiency and the precision of the sensitivity analysis. Fig. 23 shows the results for the site station, indicating that the shear wave velocity at the dam level is the most influential parameter on the amplitudes in the N–S direction. In the EW-direction, the result indicate that the shear wave velocity of the first and third layers of the Dorel model are the most influential parameter at low frequencies [0 Hz, 2 Hz], while the shear wave velocity at the dam level is most influential at high frequencies [2 Hz, 4 Hz]. Conversely, Fig. 24 shows the results for the reference station, revealing that the velocities of the first layer of the Dorel model is the most influential parameter in the N–S and V direction. In the E–W

direction, we can conclude that the shear wave velocity in the first and third layer of the Dorel model are the most influential parameters. We can also conclude that the position of the seismic source has a slight effect on the ground velocity recorded at the dam (site station) and also at the reference station in all three directions.

### Prediction interpretation using SHAP

To complement the global sensitivity analysis, we employed SHAP (SHapley Additive exPlanations), an explainable AI technique, to interpret the predictions of the surrogate model. SHAP provides insights into how input parameters contribute to the model's predictions by calculating the marginal contribution of each feature based on cooperative game theory. In this study, SHAP was applied specifically to the peak
<figure><img src="image_37.png" /><figcaption>Fig. 24. Influence of each uncertain parameter on the three components of the Fourier transform amplitudes at the reference station.</figcaption></figure>

ground velocity (PGV) in the East–West direction. Figs. 25(a), 25(b) shows the results for the site and reference station. Each point on the graph represents the impact of a specific observation on the model’s prediction. So, the farther a point is from the vertical axis, the greater the contribution of the feature to that instance’s prediction. Also, if a feature has a strong influence on the model’s output, the points will show a greater spread. The results (Figs. 25(a) and 25(b)) indicate that on the site station, the shear wave velocity (Vs) in the dam layer is the most influential parameter affecting the peak ground velocity of the dam. On the reference station, the shear wave velocities (Vs) in the first and third layers are the most impactful. Also, the global sensitivity analysis indicate that the Sobol index of the shear wave velocity in the third layer of Dorel (Vs3) is 0.25. This is the highest Sobol index among all parameters, indicating that Vs3 is the most influential parameter affecting the peak ground velocity recorded on the reference station. In the same way, the sobol index of the shear wave velocity in the dam layer is 0.4, the highest Sobol index among all parameters which shows that this parameter is the most influential parameter. These consistent findings validate the robustness of both techniques in identifying the key parameters influencing seismic behavior, further strengthening the reliability of the model’s predictions.

## Results using the Sobol method

In this section, we will compare the results obtained by both methods for each parameter, and then we will calculate the total sensitivity index for each parameter and compare it to its first-order index. Using a Monte Carlo sample of size *N*, the actual number of simulations required to estimate sensitivity indices for input variables is *N* × 2<sup>*p*</sup>, since this estimation involves two sets of simulations. However, estimating only first-order and total indices requires only *N* × (2*<i>p</i>* + 1) evaluations. Therefore, we will first estimate the first-order and total indices. If there are significant differences between these indices, it indicates that interaction effects are notable, and it may be useful to estimate intermediate-order indices. Conversely, if the differences are minimal, the input variables’ effects are primarily first-order, and intermediate-order indices may not be necessary. In this context, we generated two experimental designs of size 100,000 using quasi-Monte

<caption>Table 5</caption>

First order index and total order index for the shear wave velocity. Top: mean value of the first order index. Bottom: mean value of the total order index.

<table><thead><tr><th></th><th>V<sub>s,dam</sub></th><th>V<sub>s1</sub></th><th>V<sub>s2</sub></th><th>V<sub>s3</sub></th><th>X<sub>s</sub></th><th>Y<sub>s</sub></th><th>Z<sub>s</sub></th></tr></thead><tbody><tr><td>Mean (First order index)</td><td>0.2</td><td>0.59</td><td>0.3</td><td>0.46</td><td>0.22</td><td>0.3</td><td>0.2</td></tr><tr><td>Mean (Total order index)</td><td>0.01</td><td>0.1</td><td>0.05</td><td>0.2</td><td>0.01</td><td>0.05</td><td>0.03</td></tr></tbody></table>

Carlo sampling (based on Sobol sequences). Subsequently, we used our neural network surrogate model to construct our response surface. The comparison is specifically made for the EAST–WEST component of the Fourier amplitude transformation of velocity.

## Comparison between methods on first-order Sobol indices for the east–west component at the site station and reference station

Figs. 26 and 27 show the Sobol indices calculated by the two methods at the reference station and the site station, respectively. We can conclude that the discrepancy between the two methods is negligible for most parameters except for the source’s z-position. Indeed, the Sobol method strongly depends on the number of simulations. In Fig. 28, we present how the results of both methods vary with the number of simulations. We observe that the Li and Mahadevan method becomes stable after 30,000 simulations, whereas the Sobol method requires more than 250,000 simulations to stabilize.

## Total indices

Regarding the total indices, we present in Figs. 29 and 30 the variation of these indices with respect to frequencies, as well as the variation of first-order indices. The results show that there are significant differences between these two indices for all parameters. In Table 5, we show the maximum value of the first and the total order index for each uncertain parameter (in the reference station). This shows that there is a significant difference between the two indices, and therefore, it is interesting to calculate the intermediate-order indices to determine the influence of the interaction between the different uncertain parameters on the ground velocity recorded at the dam and at the reference station.
<figure><img src="image_38.png" /></figure>

(a) SHAP summary plot showing the influence of input parameters on the predicted peak ground velocity (PGV) in the East-West direction at the site station.

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

(b) SHAP summary plot showing the influence of input parameters on the predicted peak ground velocity (PGV) in the East-West direction at the reference station.

<figure><img src="image_40.png" /><figcaption>Fig. 26. Comparison between the results of the two methods at the reference station.</figcaption></figure>

## The bootstrap technique

In general, we aim to obtain both the estimation $\hat{S}_i$ and measures of its precision, such as bias and standard error. For many statistics, like a mean, these measures can be calculated analytically. However, for more complex statistics, analytical formulas are often unavailable, as in our case where both calculation methods used are based on sampling estimation methods. In such cases, bootstrap technique (Wehrens et al., 2000) can be employed. The Bootstrap method allows for the estimation of the uncertainty associated with Sobol indices by generating

multiple resampled datasets from the available data (experimental design and response surface). This produces an empirical distribution of the Sobol indices (mean, variance, and standard deviation of the mean) without needing to repeat *r* estimations to obtain these measures for each index estimation. In practice, bootstrap works by generating multiple samples from the original population through random sampling with replacement. For each bootstrap sample created, the estimator of interest ($S_i$) is recalculated, thereby observing the inherent variability in the estimation. This allows for assessing whether the Sobol indices are sensitive to variations in the data. This technique validates the
<figure><img src="image_41.png" /><figcaption>Fig. 27. Comparison between the results of the two methods at the site station.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>Fig. 28. Variation of results from both methods as a function of the number of simulations. From top left to bottom right panels: 1920, 3840, 7680, 15360, 30720, 61440, 122880, 245760, 491520 simulations, respectively.</figcaption></figure>

results of Sobol indices by providing a measure of their variability and verifying the stability of the conclusions obtained. The bootstrap technique is illustrated in Fig. 31, where 10 replicas were generated by random sampling with replacement. Using the Sobol method, we

estimated the Sobol index for the last uncertain parameter ($Z_{source}$) on each replica. After estimating the Sobol index for each replicate (10 results), we can evaluate the variability of these results (see Fig. 32) by calculating their mean, variance, and standard deviation of the mean.
<figure><img src="image_43.png" /><figcaption>Fig. 29. The influence of each uncertain parameter on the EST-WEST component of the Fourier transform without interaction with other parameters (First-order index) and the total influence with interaction with all other parameters (Total index) at the reference station.</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Fig. 30. Influence of each uncertain parameter on the EST-WEST component of the Fourier transform without interaction with other parameters (First-order index) and total influence with interaction with all other parameters (Total index) at the site station.</figcaption></figure>

The figure shows that the variance of the results is negligible ($10^{-6}$), and the standard deviation of the mean is also very low, indicating that the estimates are relatively close to the mean. In other words, there is little dispersion around the mean. This suggests a certain stability in the estimations.

## 5. Conclusion

In this article, we created a metamodel for the solution of the wave equation using neural networks and applied it to uncertainty analysis of seismic wave propagation in the Guadeloupe region, particularly at the
<figure><img src="image_45.png" /><figcaption>Fig. 31. The bootstrap technique. Original sample (top), Bootstrap samples (middle), and GSA in each sample for the Z-source position (bottom).</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>Fig. 32. Distribution of the Sobol index estimator: Mean (top), variance (middle), and standard deviation of the mean (bottom).</figcaption></figure>

location of the Dumanoir dam. We used the EFISPEC3D numerical code to solve the wave equation, which required constructing a 3D mesh of the domain. To simplify the analysis and account for the correlation among the geophysical parameters of the medium, we selected S-wave velocities as uncertain parameters, along with the source position, totaling 7 uncertain parameters. The model from Dorel et al. (1979) was used, highlighting the presence of three distinct layers, each with specific velocity characteristics.

We constructed a metamodel of the solution to the wave equation generated by the EFISPEC3D code based on neural networks, allowing us to generate a large set of simulations in reduced computational time. This metamodel offers significant computational efficiency. It can predict time or frequency-series nearly instantly compared to time-consuming spectral-element simulations, with up to 100,000 predictions per second. However, the study is focused on the Dumanoir dam in Guadeloupe, which limits the direct applicability of the results to other regions with different geological conditions. The surrogate
model may need significant re-tuning for other locations, particularly
in areas with different tectonic and geological features. In fact, general-
izing the surrogate model to different locations requires regenerating a
new dataset that reflects the specific geological and geophysical prop-
erties of the target region. This involves defining new input parameters
corresponding to the new site conditions and using the wave equation
model to generate the corresponding output signals. Once this dataset
is established, the neural network must be retrained to capture the
new relationships between inputs and outputs. This process ensures that
the surrogate model remains physically consistent with the site-specific
characteristics of the new region. The neural network constructed was
a multilayer perceptron regression model, consisting of an input layer
containing the uncertain parameters, with n_h hidden layers (we demon-
strated that n_h can vary depending on the objective; for example, the
results showed that for the site station, which records ground motions
considering the geological properties of the dam, a neural network
with four hidden layers is sufficient, whereas for the reference station,
which does not account for these properties, a neural network with
three hidden layers was used), and an output layer with 301 neurons
representing ground velocities at each time step and 149 neurons rep-
resenting the frequency-domain signals. We divided our approach into
three stages: the first to reproduce the time-domain signals, the second
to reproduce the amplitudes of the Fourier transform, and the third for
the phases of the Fourier transform. These three stages were performed
for both stations. For the reference station, our metamodels for the
time-domain signals and the amplitudes of the Fourier transform suc-
cessfully reproduced the signals generated by the EFISPEC3D code with
very good performance. For the phases, due to their complexity, we suc-
cessfully created metamodels for the real and imaginary parts, which
then allowed us to accurately reproduce the phases. However, for the
site station, our metamodel poorly reproduced the time-domain signals,
as the signals recorded at the dam are more complex and require more
training data. For the amplitudes of the Fourier transform, the neural
network successfully reproduced the frequency-domain series, while for
the real and imaginary parts, we observed good performance of the
metamodel at low frequencies, with less precision at high frequencies.
This issue can potentially be resolved by generating more training data
using EFISPEC3D code specifically from the dam area. By doing so,
the neural network can better learn the relationship between uncertain
parameters and the recorded signals at the dam. This approach aims
to improve the network's ability to generalize and accurately predict
seismic responses under different conditions, accounting for the specific
geological and structural characteristics that affect signal propagation
at the site station.

Next, we applied the Sobol indices method for a global sensitivity
analysis to determine the influence of each uncertain parameter on the
amplitude of the Fourier transformation. We employed two sampling-
based estimation methods: the first developed by Sobol to calculate all
indices (first and higher order), revealing that the interaction effects
among uncertain parameters were not negligible. The second method
developed by Li and Mahadevan to calculate only the index of first
order, confirmed our estimations by comparing results with those from
the Sobol method. The findings indicated that velocities in the first
and third layers of the Dorel model exerted the most influence on
seismic wave propagation recorded in the reference station. However,
the results showed that at the site station, the shear wave velocity
of the dam layer is the most influential parameter on the Fourier
transform amplitudes. The study simplifies the seismic source and P-
wave velocities by treating them as known parameters. This limits
the generalizability of the model, as real-world seismic scenarios of-
ten involve significant uncertainties in these factors. Expanding the
uncertain parameter space add computational complexity. Finally, we
employed the bootstrap technique to assess the stability of our esti-
mations. Results showed that the estimates were relatively close to
a mean value, suggesting stability in the estimations. Looking ahead,
our future research aims to develop a robust and efficient digital twin

capable of providing realistic and responsive seismic predictions for the
Dumanoir dam. Several perspectives will be also explored to further
enhance the model's accuracy and generalizability. First, we plan to
introduce fluctuations into the velocity structure, particularly in the
most influential layers, to improve the accuracy of predictions when
compared with real data. This will allow us to refine the model's sen-
sitivity to variations in geological properties. Additionally, we aim to
incorporate other uncertain parameters, such as P-wave velocities and
the earthquake source mechanism, into the model. This expansion will
help us better generalize the model, providing more robust predictions
across a broader range of scenarios.

CRediT authorship contribution statement

**Mohamad Ali Noureddine:** Writing – review & editing, Visual-ization, Validation, Supervision, Project administration, Funding acquisition. **Florent De Martin:** Writing – review & editing, Writing – original draft, Validation, Supervision, Funding acquisition. **Rani El Meouche:** Writing – review & editing, Writing – original draft, Validation, Supervision, Project administration, Funding acquisition. **Muhammad Ali Sammuneh:** Writing – review & editing, Writing – original draft, Supervision, Project administration, Funding acquisition. **Fakhreddine Ababsa:** Writing – review & editing, Writing – original draft, Validation, Supervision, Project administration, Funding acquisition. **Mickael Beaufils:** Writing – review & editing, Writing – original draft, Validation, Supervision, Project administration, Funding acquisition.

Declaration of generative AI and AI-assisted technologies in the
writing process

During the preparation of this work M.A.N used ChatGPT in order
to enhance clarity and improve the quality of the text. After using this
service, the authors reviewed and edited the content as needed and
takes full responsibility for the content of the publication.

Declaration of competing interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Acknowledgments

This research was conducted and funded within the framework of
the “Digital Twins of Construction and Infrastructure in their Envi-
ronment” research chair at ESTP. We sincerely thank Almohamad, D.,
Agastin L., Bascle A., Colombain A., and Antoinet, E. for deploying the
sensors and collecting the data (not presented in this article). Link to
the research chair: https://www.estp.fr/en/digital-twins-construction-
and-infrastructure-their-environment.

Data availability

Data will be made available on request.

References

Abdar, M., Pourpanah, F., Hussain, S., Rezazadegan, D., Liu, L., Ghavamzadeh, M., Fieguth, P., Cao, X., Khosravi, A., Acharya, U.R., et al., 2021. A review of uncertainty quantification in deep learning: Techniques, applications and challenges. Inf. Fusion 76, 243–297.

Baudin, M., Martinez, J.-M., 2010. Polynômes de chaos sous scilab via la librairie NISP. In: 42èmes Journées de Statistique.

Beale, M.H., Hagan, M.T., Demuth, H.B., 2010. Neural network toolbox. User's Guid. MathWorks 2, 77–81.
Berthier, C.-H., 2005. Quantification des incertitudes des débits calculés par un modèle pluie-débit empirique. Université Paris Sud XI, Paris.

Bottou, L., Curtis, F.E., Nocedal, J., 2018. Optimization methods for large-scale machine learning. SIAM Rev. 60 (2), 223-311.

Boutin, M., 1994. Sur Certains Algorithmes D'apprentissage De Type Retropropagation Dans Les Reseaux De Neurones. Université de Sherbrooke.

Brocher, T.M., 2005. Empirical relations between elastic wavespeeds and density in the Earth's crust. Bull. Seismol. Soc. Am. 95 (6), 2081-2092.

Cacuci, D.G., 1981. Sensitivity theory for nonlinear systems. I. Nonlinear functional analysis approach. J. Math. Phys. 22 (12), 2794-2802.

Cao, Y., Chandrasekar, A., Radhika, T., Vijayakumar, V., 2024a. Input-to-state stability of stochastic Markovian jump genetic regulatory networks. Math. Comput. Simulation 222, 174-187.

Cao, Y., Subhashri, A., Chandrasekar, A., Radhika, T., Przybyszewski, K., 2024b. Exponential state estimation for delayed competitive neural network via stochastic sampled-data control with Markov jump parameters under actuator failure. J. Artif. Intell. Soft Comput. Res. 14 (4), 373-385.

Cao, D., Wang, Y., Duan, J., Zhang, C., Zhu, X., Huang, C., Tong, Y., Xu, B., Bai, J., Tong, J., et al., 2020. Spectral temporal graph neural network for multivariate time-series forecasting. Adv. Neural Inf. Process. Syst. 33, 17766-17778.

Castaings, W., 2007. Analyse de sensibilité et estimation de paramètres pour la modélisation hydrologique: potentiel et limitations des méthodes variationnelles (Ph.D. thesis). Université Joseph-Fourier-Grenoble I.

Chaljub, E., Maufroy, E., Moczo, P., Kristek, J., Hollender, F., Bard, P.-Y., Priolo, E., Klin, P., De Martin, F., Zhang, Z., et al., 2015. 3-D numerical simulations of earthquake ground motion in sedimentary basins: testing accuracy through stringent models. Geophys. J. Int. 201 (1), 90-111.

Chapman, C., 2004. Fundamentals of Seismic Wave Propagation. Cambridge University Press.

Dai, H., Zhang, F., Ye, M., Guadagnini, A., Liu, Q., Hu, B., Yuan, S., 2022. A computationally efficient method for estimating multi-model process sensitivity index. Water Resour. Res. 58 (10), e2022WR033263.

De Martin, F., 2011. Verification of a spectral-element method code for the Southern California Earthquake Center LOH.3 viscoelastic case. Bull. Seismol. Soc. Am. 101 (6), 2855-2865.

De Martin, F., Chaljub, E., Thierry, P., Sochala, P., Dupros, F., Maufroy, E., Hadri, B., Benaichouche, A., Hollender, F., 2021. Influential parameters on 3-D synthetic ground motions in a sedimentary basin derived from global sensitivity analysis. Geophys. J. Int. 227 (3), 1795-1817, arXiv:https://academic.oup.com/gji/article-pdf/227/3/1795/40084089/ggab304.pdf. URL https://doi.org/10.1093/gji/ggab304.

De Moura, C.A., Kubrusly, C.S., 2013. The courant-friedrichs-lewy (cfl) condition. AMC 10 (12).

De Zanet, G., Viquerat, A., 2022. Screening methods for sensitivity analysis applied to thin composite laminated structures. Thin-Walled CStruct. 172, 108870.

Dorel, J., Eschenbrenner, S., Feuillard, M., 1979. Coupes sismiques des structures superficielles dans les petites antilles—I: Guadeloupe. Pure Appl. Geophys. 117, 1050-1069.

Driancourt, X., 1994. Optimisation Par descente de Gradient Stochastique de Systèmes Modulaires Combinant Reenseaux de Neurones Et Programmation Dynamique. Application a La Reconnaissance de La Parole (Ph.D. thesis). Paris 11.

Du, W., Ma, J., Zhou, C., Yan, Y., Wriggers, P., 2023. Uncertain dynamic characteristic analysis for structures with spatially dependent random system parameters. Materials 16 (3), 1188.

Frajraoui, N., 2014. Analyse de sensibilité globale et polynômes de chaos pour l'estimation des paramètres: Application aux transferts en milieu poreux (Ph.D. thesis). Strasbourg.

Farhadi, F., 2017. Learning Activation Functions in Deep Neural Networks. Ecole Polytechnique, Montreal (Canada).

Fernique, X., 1971. Régularité de processus gaussiens. Invent. Math. 12, 304-320.

Gomes, G.S.d.S., Ludermir, T.B., 2013. Optimization of the weights and asymmetric activation function family of neural network for time series forecasting. Expert Syst. Appl. 40 (16), 6438-6446.

Hammersley, J., 2013. Monte Carlo Methods. Springer Science & Business Media.

Homma, T., Saltelli, A., 1996. Importance measures in global sensitivity analysis of nonlinear models. Reliab. Eng. Syst. Saf. 52 (1), 1-17. http://dx.doi.org/10.1016/0951-8320(96)00002-6, URL https://www.sciencedirect.com/science/article/pii/0951832096000026.

Hu, R., Monebhurrun, V., Himeno, R., Yokota, H., Costen, F., 2022. An uncertainty analysis on finite difference time-domain computations with artificial neural networks: improving accuracy while maintaining low computational costs. IEEE Antennas Propag. Mag. 65 (1), 60-70.

Hüllermeier, E., Waegeman, W., 2021. Aleatoric and epistemic uncertainty in machine learning: An introduction to concepts and methods. Mach. Learn. 110 (3), 457-506.

Iooss, B., 2011. Revue sur l'analyse de sensibilité globale de modèles numériques. J. Soc. Française Stat. 152 (1), 3-25.

Jacques, J., 2011. Pratique de l'analyse de sensibilité: comment évaluer l'impact des entrées aléatoires sur la sortie d'un modèle mathématique. Lille: Sn.

Jepkoech, J., Mugo, D.M., Kenduiywo, B.K., Too, E.C., 2021. The effect of adaptive learning rate on the accuracy of neural networks. Int. J. Adv. Comput. Sci. Appl. 12 (8).

Jiménez Rugama, L.A., Gilquin, L., 2018. Reliable error estimation for Sobol' indices. Stat. Comput. 28 (4), 725-738.

Kendall, A., Gal, Y., 2017. What uncertainties do we need in bayesian deep learning for computer vision? Adv. Neural Inf. Process. Syst. 30.

Kingma, D.P., Ba, J., 2014. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.

Kleijnen, J.P., 1997. Sensitivity analysis and related analyses: a review of some statistical techniques. J. Stat. Comput. Simul. 57 (1-4), 111-142.

Kolla, H., De, S., Jones, R.E., Hansen, M.A., Plewes, J.A., 2022. Comprehensive Uncertainty Quantification (UQ) for Full Engineering Models by Solving Probability Density Function (PDF) Equation. Technical report, Sandia National Lab.(SNL-CA), Livermore, CA (United States); Sandia National ....

Komatitsch, D., 1997. Méthodes Spectrales et Éléments Spectraux Pour l'équation de l'élastodynamique 2D et 3D en milieu hétérogène (Ph.D. thesis). Institut de physique du globe de paris-IPGP.

Lambert, B., Forbes, F., Doyle, S., Dehaene, H., Dojat, M., 2024. Trustworthy clinical AI solutions: a unified review of uncertainty quantification in deep learning models for medical image analysis. Artif. Intell. Med. 102830.

Li, Z., Liu, F., Yang, W., Peng, S., Zhou, J., 2021. A survey of convolutional neural networks: analysis, applications, and prospects. IEEE Trans. Neural Netw. Learn. Syst. 33 (12), 6999-7019.

Li, C., Mahadevan, S., 2016. An efficient modularized sample-based method to estimate the first-order Sobol' index. Reliab. Eng. Syst. Saf. 153, 110-121.

Lucibella, M., 2014. Robert hooke (1635-1703). Reflets Phys. (40), 24-25.

Lyu, Z., Yu, Y., Samali, B., Rashidi, M., Mohammadi, M., Nguyen, T.N., Nguyen, A., 2022. Back-propagation neural network optimized by K-fold cross-validation for prediction of torsional strength of reinforced concrete beam. Materials 15 (4), 1477.

Matsushima, S., Hirokawa, T., De Martin, F., Kawase, H., Sánchez-Sesma, F.J., 2014. The effect of lateral heterogeneity on horizontal-to-vertical spectral ratio of microtremors inferred from observation and synthetics. Bull. Seismol. Soc. Am. 104 (1), 381-393.

Maufroy, E., Chaljub, E., Hollender, F., Kristek, J., Moczo, P., Klin, P., Priolo, E., Iwaki, A., Iwata, T., Etienne, V., De Martin, F., Theodulidis, N., Manakou, M., Guyonnet-Benaize, C., Pitilakis, K., Bard, P.-Y., 2015. Earthquake ground motion in the mygdonian Basin, Greece: The E2VP verification and validation of 3D numerical simulation up to 4 Hz. Bull. Seismol. Soc. Am..

Meyer, P., Priouret, P., Spitzer, F., Priouret, P., 1974. Processus de diffusion et équations différentielles stochastiques. In: Ecole d'Eté de Probabilités de Saint-Flour III-1973. Springer, pp. 37-113.

Miserer, R., Biegler, L., 2023. Formulating data-driven surrogate models for process optimization. Comput. Chem. Eng. 179, 108411.

Pan, Y., Li, Y., 2023. Toward understanding why adam converges faster than sgd for transformers. arXiv preprint arXiv:2306.00204.

Popescu, M.-C., Balas, V.E., Perescu-Popescu, L., Mastorakis, N., 2009. Multilayer perceptron and neural networks. WSEAS Trans. Circuits Syst. 8 (7), 579-588.

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Tarantola, S., 2008. Global Sensitivity Analysis: the Primer. John Wiley & Sons.

Saltelli, A., Ratto, M., Tarantola, S., Campolongo, F., et al., 2006. Sensitivity analysis practices: Strategies for model-based inference. Reliab. Eng. Syst. Saf. 91 (10-11), 1109-1125.

Schmidhuber, J., 2015. Deep learning in neural networks: An overview. Neural Netw. 61, 85-117.

Schmitz, V., Arnst, M., El kadi Abderrezzak, K., Pirotton, M., Erpicum, S., Archambau, P., Dewals, B., 2023. Global sensitivity analysis of a dam breaching model: To which extent is parameter sensitivity case-dependent? Water Resour. Res. 59 (6), e2022WR033894.

Sharma, S., Sharma, S., Athaiya, A., 2017. Activation functions in neural networks. Towar. Data Sci. 6 (12), 310-316.

Sobol', I.M., 1990. On sensitivity estimation for nonlinear mathematical models. Mat. Model. 2 (1), 112-118.

Sobol, I., 1993. Sensitivity estimates for nonlinear mathematical models. Math. Model. Comput. Exp. 1, 407.

Sobol, I.M., 2001. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math. Comput. Simulation 55 (1-3), 271-280.

Sochala, P., De Martin, F., Le Maitre, O., 2020. Model reduction for large-scale earthquake simulation in an uncertain 3D medium. Int. J. Uncertain. Quantif. 10 (2), 101-127.

Song, X., Liu, Y., Xue, L., Wang, J., Zhang, J., Wang, J., Jiang, L., Cheng, Z., 2020. Time-series well performance prediction based on long short-term memory (LSTM) neural network model. J. Pet. Sci. Eng. 186, 106682.

St, L., Wold, S., et al., 1989. Analysis of variance (ANOVA). Chemometr. Intell. Lab. Syst. 6 (4), 259-272.

Stein, M., 1987. Large sample properties of simulations using latin hypercube sampling. Technometrics 143-151.

Sudret, B., 2008. Global sensitivity analysis using polynomial chaos expansions. Reliab. Eng. Syst. Saf. 93 (7), 964-979.

Tan, A.R., Urata, S., Goldman, S., Dietschreit, J.C., Gómez-Bombarelli, R., 2023. Single-model uncertainty quantification in neural network potentials does not consistently outperform model ensembles. Npj Comput. Mater. 9 (1), 225.
Touzani, S., 2011. Méthodes de surface de réponse basées sur la décomposition de la variance fonctionnelle et application à l'analyse de sensibilité (Ph.D. thesis). Université de Grenoble.

Wang, F., Wang, Y., Zhang, K., Hu, M., Weng, Q., Zhang, H., 2021. Spatial heterogeneity modeling of water quality based on random forest regression and model interpretation. Environ. Res. 202, 111660.

Wehrens, R., Putter, H., Buydens, L.M., 2000. The bootstrap: a tutorial. Chemometr. Intell. Lab. Syst. 54 (1), 35–52.

Weiss, N.A., Holmes, P.T., Hardy, M., 2006. A Course in Probability. Pearson Addison Wesley Boston, MA, USA:

Zhang, X., Zhang, Y., Pandey, M.D., 2015. Global sensitivity analysis of a CNC machine tool: application of MDRM. Int. J. Adv. Manuf. Technol. 81, 159–169.