# Uncertainty quantification for acoustic wave propagation in a shallow water environment

Shahram Khazaie, Xun Wang, Dimitri Komatitsch, Pierre Sagaut

► To cite this version:

Shahram Khazaie, Xun Wang, Dimitri Komatitsch, Pierre Sagaut. Uncertainty quantification for acoustic wave propagation in a shallow water environment. Wave Motion, 2019, 91, pp.102390. 10.1016/j.wavemoti.2019.102390. hal-02467993

HAL Id: hal-02467993

https://hal.science/hal-02467993

Submitted on 13 Feb 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Uncertainty quantification for acoustic wave propagation in a shallow water environment

Shahram Khazaie<fnref n="1" />ª, Xun Wangᵇ, Dimitri Komatitsch<fnref n="2" />ᶜ, Pierre Sagautᵈ

ªInstitut de Recherche en Génie-Civil et Mécanique, UMR CNRS 6183, Université de Nantes, Saint Nazaire, France

ᵇDepartment of Civil and Environmental Engineering, Hong Kong University of Science and Technology, Hong Kong

ᶜAix-Marseille University, CNRS, Centrale Marseille, LMA, Marseille 13453, France

ᵈAix-Marseille University, CNRS, Centrale Marseille, M2P2, Marseille 13451, France

## Abstract

Sound wave propagation in a shallow water environment is complex due to e.g. the uncertainties of sound speed profile being inhomogeneous and imprecisely measured, the bottom reflections, etc. The propagation and influence of several uncertainty parameters are quantified in this paper. A four-layer model, which can approximately represent a wide range of shallow water environments, is considered; six parameters representing sound speed profile and water depth are considered as random variables. We investigate how the wave field (pressure) in this model is influenced by these uncertainties. For this purpose, the sound field is computed for different realizations of the random variables, when the medium is excited with sources whose frequencies are appropriate, for example, for marine seismic exploration applications. Since classical Monte Carlo methods require a huge sample size to converge, we use three surrogate modeling techniques (Kriging, Polynomial Chaos, and Polynomial Chaos-based Kriging). The proposed methods require much smaller sample sizes, which makes the uncertainty quantification (UQ) possible. Wavelength-to-depth ratio ($\lambda/d$) is introduced as the key parameter that defines the degree of interaction (reflection and transmission) of the sound waves with the boundaries of the shallow water waveguide. The results show that for small and large values of $\lambda/d$, the wave field is more sensitive to the variations of the water depth and the velocity of the bottom layer, respectively. The robustness (precision) of the surrogate models is shown to decrease for lower values of $\lambda/d$. The proposed UQ methodology can be used for more complicated underwater environments; it is even more advantageous because it can efficiently deal with a large number of model uncertainty parameters and identify the most influential ones.

**Keywords:** Uncertainty propagation, Wave propagation, Surrogate models, Shallow water

<footnote n="1">¹Electronic address: shahram.khazaie@univ-nantes.fr</footnote>

<footnote n="2">²Deceased on 21 January 2019.</footnote>
Introduction

Knowledge of the Green's function (impulse response) of the medium is essential in the majority of wave
propagation applications, for instance medical imaging, underwater communication and structural health
monitoring, among others. In order to model the acoustic wave field (pressure or velocity) [4, 6, 22, 30, 32, 33,
38], one should precisely define the values of the environmental parameters (typically velocity and density at
each point) and the boundary and/or initial conditions. The physical properties of the propagation medium
are inherently inhomogeneous. For instance, in ocean acoustics, depth, pressure, salinity and temperature
are among the influencing factors on the propagation velocity [21]. Drastic variations on the sound speed
profiles are observed particularly near the sea surface. Hence, the environmental parameters are spatially
varying and only their statistical information might be provided [24, 52]. Taking into account this spatial
variability plays an important role in modeling of both direct and inverse (source identification) problems.

The values of the input parameters of the above-described inhomogeneous medium are usually under-
sampled and poorly known. As a result, the uncertainties due to our lack of knowledge should be accounted
for in the modeling process. This adds a probabilistic dimension to the classical wave propagation equation.
The mentioned uncertainties exist also in boundary and initial conditions. In this regard, Khazaie et al.
[24] studied acoustic wave propagation, accounting for the uncertainties related to the source. However, the
uncertainties related to the boundary and initial conditions are out of the scope of the present paper and
the main focus is on inserting uncertainty on the physical properties of the propagation medium.

In this work we consider an axisymmetric model which is also known as 2.5-D and therefore it can
describe a 3-D inhomogeneous shallow water environment with an underlying semi-infinite ocean layer. The
spatial variations on the density are considered as negligible compared to that of the sound speed so that it
is considered as a constant for each layer. The inhomogeneous shallow water sound speed is considered to be
only depth-dependent (horizontally invariant) and is characterized with two constants along with a (quasi-
)linear variation in between. The velocity and density of the ocean layer are both constants. This model
was elaborated in Wang et al. [55] and was shown to efficiently approximate shallow water environments in
particular during summer [27].

The uncertainty on the aforementioned sound speed profile of the whole system is modeled via six random
variables with given lower and upper bounds. These random variables correspond to: (i) two constant sound
speeds of the shallow water; (ii) the sound speed of the ocean bottom; and (iii)three layer depths in the
shallow water. One of the objectives of this work is to propagate the input uncertainties through the system
and to classify the input parameters in terms of their contribution to the variance (uncertainty) of the wave
field via a sensitivity analysis. For instance [36] used the same approach for wave propagation through
poroelastic environments. It should be noted that hereinafter the wave field means the absolute value of
the Fourier transform of the temporal pressure field filtered around the dominant frequency of the source.
The relatively high number of random variables involved in this problem forces us to solve deterministic
realizations of the problem at least several thousands times to obtain a wave field with converged statistics.
On the other hand, the computational cost of the wave field modeling for each particular realization of such
a medium is a *priroi* high, particularly in three-dimensional problems. To cope with this problem, we use
the fast algorithm proposed in [55] to solve the wave equation. This method is based on the wavenumber
integration [6, 33, 54] and provides a semi-analytical Green's function for the system that enables drastic
reductions of the computational costs. As a result, we can use classical sampling methods such as Quasi
Monte Carlo (QMC) to calculate the statistics of the random wave field.

Several studies have been devoted to assessing the uncertainty effects on the pressure wave field. For
instance, Lermusiaux et al. [31] studied the sensitivity of the transmission loss (TL) to several factors such
as uncertainties in geoacoustic parameters, seabed thickness and bathymetry. Spectral stochastic method
has been used to predict the sound pressure in an uncertain waveguide [10–12, 25]. In these works, the UQ
analysis is done using a Polynomial Chaos Expansion (PCE) method without any further discussion about
its robustness nor about the output sensitivity with respect to the input randomness. However, Gerdes and
Finette [13] studied the relative importance of tidal forcing, source depth and seabed parameters on the
spatial structure and variance of intensity using a PCE method. The PCE method has been used in Pekeris
waveguides to identify the sound field uncertainties arising from the environmental ones [19]. The former
has also been used as a substitute for the perturbation technique which is limited to weak uncertainties [37].
The Kriging surrogate model is also used to quantify the sound pressure uncertainty [17]. Wang et al. [53]
performed a two-parameter UQ using the Kriging method in the context of the sound source localization in
shallow water environments. Our paper aims at introducing the randomness on the parameters of the four
layer model introduced in Wang et al. [55] and then doing a full UQ analysis using and comparing three
approaches: the PCE, Kriging and PCE-based Kriging. It should be noted that the methodology could be
used for more complicated environments but here for comparison with the QMC results, we consider the
four layer model for which the semi-analytical Green's function is derived in [55]. Kriging yields a local
interpolation of the random pressure and provides a local error estimate. While PCE approximates the
global behavior of the computational model. The PCE-based Kriging has both of these features and can
thus be expected to generate better results. The robustness of the results obtained by these methods is
found to be determined by the wavelength-to-depth ratio. The latter also shows to be influential on the sensitivity indices.

The outline of the paper is the following. In Section 2, we recall the model of the propagation medium and the calculation of its Green's function. We then specify the parameters that will be modeled as random variables based on some given statistical parameters. In Section 3, the classical uncertainty quantification methods used in ocean acoustics are briefly recalled. In Section 4 we then present the numerical case studies and the results obtained based on the UQ methods. Conclusions are drawn in Section 5.

## 2. Sound propagation in a four-layer shallow water environment

The sound speed in shallow water region {$\mathbf{r} = (x, y, z) : z \in [0, d]$}, where $d \in \mathbb{R}^+$ is the water depth, is horizontally stratified and vertically follows a three-layer profile consisting of two separate constant layers, with a layer varying linearly with respect to $k^2(z)$ sandwiched in between, as shown in Figure 1. Here, $k(z) = 2\pi f/v(z)$ is the wavenumber, $v(z)$ is the sound speed, and $f$ is the frequency. The bottom is assumed to be an infinite fluid halfspace {$\mathbf{r} = (x, y, z) : z \in [d, \infty)$} where the sound speed is constant. The sound speed function with respect to the depth $z$ is concluded as follows:

$$v(z) = \begin{cases} v_1, & z \in [0, d_1) \\ \sqrt{\frac{1}{az+b}}, & a > 0, z \in [d_1, d_2) \\ v_2, & z \in [d_2, d) \\ v_3, & z \in [d, \infty) \end{cases}, \qquad (1)$$

where the coefficients $a$ and $b$ are

$$a = \frac{v_1^2 - v_2^2}{v_1^2 v_2^2 (d_2 - d_1)}, \quad b = \frac{v_2^2 d_2 - v_1^2 d_1}{v_1^2 v_2^2 (d_2 - d_1)}, \qquad (2)$$

such that the sound speed is continuous at $z = d_1$ and $z = d_2$. In both the water column and the ocean bottom regions, the density is constant:

$$\rho(z) = \begin{cases} \rho_1, & z \in [0, d] \\ \rho_2, & z \in (d, \infty) \end{cases}. \qquad (3)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Shallow water environment with a four-layer sound speed profile $v(z)$.</figcaption></figure>

The Helmholtz equation for the displacement potential, being a function of the spatial coordinate $\mathbf{r} = (x, y, z)$ and the frequency $f$, reads

$$ (\nabla^2 + k^2(z)) \psi(\mathbf{r}, f) = \delta(\mathbf{r} - \mathbf{r}^0)S_f, \quad (4) $$

where $S_f$ is a deterministic function and $\mathbf{r}^0 = (x^0, y^0, z^0)$, $z^0 \in (0, d)$ is the source position. The sound pressure is obtained from the displacement potential as $p(\mathbf{r}, f) = \rho\omega^2\psi(\mathbf{r}, f)$ where $\omega = 2\pi f$ is the angular frequency. The Helmholtz equation (4) can be solved via a wavenumber integration approach [55], whose computation is fast but still needs a one-fold numerical integration. This method is used in this paper to compute the sound propagation in the four-layer model.

### 3. Surrogate modeling and sensitivity analysis for uncertainty quantification

In this section, the uncertainty quantification methods used in this paper are introduced. Consider the random pressure field (at any point $\mathbf{x}$) or any other scalar stochastic output quantity of interest (QoI) $Y$ as a function of a vector of $M$ input random variables $\mathbf{X} = \{X_1, X_2, ..., X_M\}^\top \in \mathbb{R}^M$, i.e. $Y = \mathcal{P}(\mathbf{r}, \mathbf{X})$. Hereinafter the dependence on $\mathbf{r}$ is omitted for clarity purposes. Classical sampling techniques such as Monte Carlo (MC) simulation, which takes random samples from $\mathbf{X}$, could be used to estimate the probabilistic content of the QoI. However, for large stochastic dimensions $M$, a relatively high number of samples is needed to guarantee the convergence. This implies high computational costs which makes MC methods intractable, even if the solution of the wave equation is semi-analytical (i.e. only a one-fold numerical integration is required). Hence, meta-models could alternatively be used to build fast-to-evaluate surrogate models from a limited number of calls to the numerical solver. In this work, Kriging [40, 41, 43, 49] and Polynomial Chaos Expansion (PCE) [1, 3, 15, 50] as two popular surrogate models will be used. The PCE-based Kriging
method [45, 46] will also be employed as a variant of the Kriging method. These surrogate models will be briefly discussed in Sections 3.1–3.3. Based on the three surrogate modeling methods, we further introduce the Sobol index to quantify the global sensitivities in Section 3.4.

### 3.1. Kriging (Gaussian process regression)

Kriging is based on the decomposition of the QoI into a deterministic mean which is referred to as trend or regression part, and a residual stochastic term that is assumed to be a Gaussian process [43]:

$$ \mathcal{P}^K(\mathbf{x}) = \mathbf{a}^\top \mathbf{f}(\mathbf{x}) + \sigma Z(\mathbf{x}; \ell) \quad (5) $$

wherein $\mathbf{x}$ is a realization of the random vector $\mathbf{X}$, $\mathbb{E}[\mathcal{P}^K] = \mathbf{a}^\top \mathbf{f}(\mathbf{x})$ is the mean of the QoI, represented by a set of predefined shape functions $f_i, i \in \{0, ..., P-1\}$ and coefficients vector $\mathbf{a}$, and $Z(\mathbf{x}; \ell)$ is a zero-mean, unit-variance stationary Gaussian random field. The standard deviation of the Kriging’s residual is denoted by $\sigma$. Note that the superscript $\top$ represents the transpose operation. The random field $Z(\mathbf{x}; \ell)$ is fully described by its auto-correlation function (ACF) $R(\mathbf{x}, \mathbf{x}'; \ell) = \mathbb{E}[Z(\mathbf{x}; \ell)Z(\mathbf{x}'; \ell)]$. The ACF is often considered as a function of the lag distance $d = |\mathbf{x} - \mathbf{x}'|$ and a vector $\ell = \{\ell_1, ..., \ell_M\}^\top$ also called the correlation distance, i.e. $R(\mathbf{x}, \mathbf{x}'; \ell) = R(d; \ell)$. The vector $\ell$ is oftentimes unknown and should be estimated from the pressure measurements. Some frequently used ACF models are introduced in [23, 26].

The objective is to evaluate the unknown parameters using a sample of size $N$ (aka experimental design), $\mathcal{X} = \{\mathbf{x}_1, ..., \mathbf{x}_N\}$ and the corresponding pressure measurements $\mathcal{Y} = \{y_1, ..., y_N\}^\top$ in which $y_j = \mathcal{P}(\mathbf{x}_j)$, $j \in \{1, ..., N\}$ is the recorded pressure for $j^{th}$ realization of input parameters. The Kriging predictor at any point (realization of input parameters) $\mathbf{x}$ is a Gaussian random variable $\hat{Y}(\mathbf{x}) \sim \mathcal{N}(\mu_{\hat{Y}}(\mathbf{x}), \sigma_{\hat{Y}}(\mathbf{x}))$. The least-squares estimate of the $P \times 1$ coefficients vector $\mathbf{a}$ is $\hat{\mathbf{a}} = (\mathbf{F}^\top \mathbf{R}^{-1} \mathbf{F})^{-1} \mathbf{F}^\top \mathbf{R}^{-1} \mathcal{Y}$ in which the $N \times P$ matrix $\mathbf{F}$ is defined as $F_{ij} = f_j(\mathbf{x}_i)$, $i \in \{1, ..., N\}$, $j \in \{1, ..., P\}$ and $\mathbf{R}$ is the correlation matrix whose components are $R_{ij} = R(|\mathbf{x}_i - \mathbf{x}_j|; \ell)$, $i, j \in \{1, ..., N\}$. Note that the dependence of the correlation matrix $\mathbf{R}$ on $\ell$ is omitted for the sake of clarity. The mean and variance of the Kriging estimator could be obtained as a function of $\ell$ (see [43] for more details):

$$ \mu_{\hat{Y}}(\mathbf{x}) = \mathbf{f}(\mathbf{x})^\top \hat{\mathbf{a}} + \mathbf{r}(\mathbf{x})^\top \mathbf{R}^{-1}(\mathcal{Y} - \mathbf{F}\hat{\mathbf{a}}), \quad (6) $$

$$ \sigma_{\hat{Y}}^2(\mathbf{x}) = \sigma^2 (1 - \mathbf{r}(\mathbf{x})^\top \mathbf{R}^{-1} \mathbf{r}(\mathbf{x}) + \mathbf{u}(\mathbf{x})^\top (\mathbf{F}^\top \mathbf{R}^{-1} \mathbf{F})^{-1} \mathbf{u}(\mathbf{x})) \quad (7) $$

in which $\mathbf{r}(\mathbf{x})$ specifies the correlation between each realization point $\mathbf{x}_i$ and the prediction point $\mathbf{x}$, i.e.
$\mathbf{r}_i(\mathbf{x}) = R(|\mathbf{x} - \mathbf{x}_i|; \ell), i \in \{1, ..., N\}$ and the $P \times 1$ vector $\mathbf{u}(\mathbf{x})$ is defined as $\mathbf{u}(\mathbf{x}) = \mathbf{F}^\top \mathbf{R}^{-1} \mathbf{r}(\mathbf{x}) - \mathbf{f}(\mathbf{x})$. Note also that the mean Kriging estimator interpolates the data, i.e. $\mu_{\hat{Y}}(\mathbf{x}_i) = \mathcal{P}(\mathbf{x}_i)$ and $\sigma_{\hat{Y}}^2(\mathbf{x}_i) = 0, i \in \{1, ..., N\}$. The main steps of the Kriging method are summarized below.

The first step consists in choosing the trend function. Constant, linear or higher order polynomials are often used. We employ an unknown constant so that $\mathbf{a}^\top \mathbf{f}(\mathbf{x}) = a_0$ is to be determined. This variant of the Kriging method is referred to as ordinary Kriging. The second step consists in choosing a model for the ACF, $R$. The von Kármán ACF will be used in this work [26], i.e. $R(\eta = h/\ell, H) = 2(\sqrt{H\eta})^H K_H(2\sqrt{H\eta})/\Gamma(H)$ wherein $H$ is the shape parameter, $K_H$ is the modified Bessel function of the second kind and of order $H$, $\Gamma$ is the Gamma function and $\ell = \ell_i, i \in \{1, ..., M\}$. We use a shape parameter of $H = 1.5$. The next step consists in setting an optimization problem and solve it for the vector $\ell$ using the experimental design $\mathcal{X}$ and the values of the QoI $\mathcal{Y}$. For this purpose, maximum likelihood and cross-validation are frequently used [7, 35]. In this work, the maximum likelihood method will be used that leads to the following optimization problem for the vector $\ell$ [28]:

$$ \ell = \arg \min_{\mathcal{D}_{\theta}=\mathbb{R}^{+M}} \left( \frac{1}{2}\log(\det(\mathbf{R})) + \frac{N}{2}\log(2\pi\sigma^2) + \frac{N}{2} \right). \quad (8) $$

In this study, a hybrid genetic algorithm [9] with upper and lower bounds of [0.001, 10] for each $\ell_i$ ($1 \le i \le M$) is used to solve the optimization problem (8) via Matlab's built-in optimization functions. Finally, one should use the optimal values of $\ell$ to estimate the coefficients of the trend $\hat{\mathbf{a}}$, the mean $\mu_{\hat{Y}}(\mathbf{x})$ and the variance $\sigma_{\hat{Y}}^2(\mathbf{x})$ of the Kriging predictor.

A leave-one-out (LOO) cross-validation error estimator $\epsilon_{LOO}$ is frequently used to assess the quality of the meta-models and/or to choose the best meta-model among several candidates. It is based on the creation of $N$ surrogate models each of which lacks one point of the experimental design, i.e. $\mathcal{X}^{(-j)} = \{\mathbf{x}_i, i = 1, ..., N, i \neq j\}$. One should then compare the corresponding predictions with the real value of the QoI at the excluded point [2]. Dubrule [8] derived an analytical solution for $\epsilon_{LOO}^K$ without computing the $N$ meta-models.

## 3.2. Polynomial Chaos Expansion (PCE)

PCE belongs to the spectral stochastic methods [14, 29] and is based on a functional representation of the random QoI. It consists in expanding any second-order random variable into an infinite linear combination of orthogonal polynomials [44]. The orthogonality is with respect to the PDFs of the input parameters. A truncated form of the PCE of the random pressure field in which polynomials of total degree not larger than
p are used, reads:

$$ \mathcal{P}_p^{\text{PCE}}(\boldsymbol{x}) = \sum_{\boldsymbol{\alpha} \in \mathcal{A}_{M,p}} a_{\boldsymbol{\alpha}} \Psi_{\boldsymbol{\alpha}}(\boldsymbol{\zeta}) = \sum_{i=0}^{P-1} a_i \Psi_i(\boldsymbol{\zeta}) = \boldsymbol{a}^{\top} \Psi(\boldsymbol{\zeta}), \quad \text{where} \quad \mathcal{A}_{M,p} = \left\{ \boldsymbol{\alpha} \in \mathbb{N}^M : 0 \le \sum_{j=1}^{M} \alpha_j \le p \right\}, \qquad (9) $$

which includes $P = \binom{p+M}{p}$ polynomials and coefficients that are put on the column vectors $\Psi$ and $\boldsymbol{a}$, respectively. Hence, $P$ grows rapidly with $p$ and $M$ implying the necessity of large number of model resolutions. This issue makes this method less interesting for large stochastic dimensions (say, $M > 10$ which is not the case in this study) [39]. The standard random vector $\boldsymbol{\zeta} = (\zeta_1, ..., \zeta_M)$ is related to the random input vector via an isoprobabilistic transformation $\mathcal{T}$, i.e. $\mathbf{X} = \mathcal{T}(\boldsymbol{\zeta})$. In (9), the PCE is expressed in terms of the $M$-tuples $\boldsymbol{\alpha} = (\alpha_1, ..., \alpha_M)$ along with the simple indices ranging from 0 to $P-1$. Moreover, $\Psi_{\boldsymbol{\alpha}}$ ($\boldsymbol{\alpha} \in \mathcal{A}_{M,p}$) are multidimensional orthonormal polynomials under the inner product defined with the joint PDF of the random vector $\mathbf{X}$ as the weight function, i.e. $\langle \Psi_{\boldsymbol{\alpha}}, \Psi_{\boldsymbol{\beta}} \rangle_{f_{\mathbf{X}}(\boldsymbol{x})} = \delta_{\boldsymbol{\alpha}\boldsymbol{\beta}}$ being 1 when $\boldsymbol{\alpha} = \boldsymbol{\beta}$ and vanishes otherwise. The family of the orthonormal polynomials will be chosen based on the marginal PDFs of the input random variables. For instance, Hermite, Legendre, Laguerre and Jacobi polynomials are used for Gaussian, uniform, Gamma and beta random variables. For independent input random variables, the basis functions are constructed by tensorizing one dimensional bases, i.e. $\Psi_{\boldsymbol{\alpha}}(\boldsymbol{\zeta}) = \Pi_{i=1}^{M} \psi_{\alpha_i}(\zeta_i)$. Different methods to calculate the expansion coefficients are summarized in [29] from which we adopt a non-intrusive least squares approach. The latter gives the estimation $\hat{\boldsymbol{a}} = (\mathbf{F}^{\top}\mathbf{F})^{-1}\mathbf{F}^{\top}\mathbf{y}$ in which the $N \times P$ matrix $\mathbf{F}$ is defined as $F_{ij} = \Psi_j(\zeta_i)$. A sample size larger than $(M-1)P$ is empirically recommended in order to prevent the matrix $\mathbf{F}^{\top}\mathbf{F}$ to be singular. Once the coefficients are calculated, the mean and variance of the random QoI could be estimated directly via $\mathbb{E}[\mathcal{P}_p^{\text{PCE}}] = a_0$ and $\text{Var}[\mathcal{P}_p^{\text{PCE}}] = \sum_{i=1}^{P-1} a_i^2$, respectively.

Since the least-squares minimization method is employed to calculate the coefficients, one can show that similar to the Kriging method, it is unnecessary to construct $N$ meta-models to calculate the LOO error.

To choose the appropriate PCE degree $p$, we use a degree-adaptive PCE approach with degrees varying from 1 to 5. Since the least-square minimization approach allows to get an estimation of the LOO error, this approach adaptively returns the best value for the polynomial degree (in the chosen interval) that minimizes the error. Apart from the aforementioned parameters of the QoI, the sensitivity indices could also be calculated directly from the coefficients of the PCE. This will be discussed in Section 3.4.

### 3.3. PCE-based Kriging

The Kriging method interpolates the local variations of the QoI as a function of the neighboring experimental design points. On the other hand, the PCE captures its global behavior using a set of orthogonal
polynomials. The PCE-Kriging method combines these two meta-modeling techniques. This allows to approximate the global behavior of the QoI using orthonormal polynomials in the trend part of the Kriging method and the local variabilities using the Gaussian process [45, 46]. The main assumption of this method is that one can directly use the optimal set of polynomials obtained in the context of pure PCE as the optimal trend for the PCE-based Kriging model. This meta-model could thus be written as

$$ \mathcal{P}^{PCEK}(\boldsymbol{x}) = \boldsymbol{a}^{\top}\boldsymbol{\Psi}(\boldsymbol{\zeta}) + \sigma Z(\boldsymbol{x}; \boldsymbol{\ell}). \qquad (10) $$

There are two main steps to construct a PCE-based Kriging meta-model. First, determination of the optimal set of polynomials in the trend part and second, the estimation of the vector $\boldsymbol{\ell}$ as well as other Kriging parameters $\{\boldsymbol{a}, \boldsymbol{\sigma}\}$. The readers are referred to Schöbi et al. [45] for further details about different numerical approaches for these estimations. In this work, we use a sequential PCE-Kriging approach that allows to sequentially calculate the set of polynomials and the Kriging meta-model. The main hypothesis of this approach is that the optimal set of polynomials obtained in the context of pure PCE can be directly used as the optimal trend for the PCE-based Kriging model. The calibration of the meta-model as a classical Kriging model is subsequently carried out.

### 3.4. Sensitivity analysis

A summary of some efficient methods to calculate the sensitivity indices are introduced in [42]. In this study we compute the Sobol sensitivity indices [48] through quasi Monte Carlo estimators [42] for Kriging and PCE-Kriging methods. However, the PCE model allows to obtain analytical formulas for sensitivity indices explicitly in terms of the PCE coefficients.

#### 3.4.1. QMC-based sensitivity indices for Kriging and PCE-Kriging

Let $\Xi_1$ and $\Xi_2$ be two independent quasi Monte Carlo samples of random parameters whose dimension is $N \times M$, $N$ and $M$ are the number of samples and the number of uncertain parameters, respectively. If we define $\Xi_1^i$ the matrix $\Xi_1$ whose $i$-th column is replaced by the $i$-th column of $\Xi_2$, an estimation for the Sobol index of the random variable $X_i$ ($1 \le i \le M$) writes:

$$ S_i = \frac{\sigma_i^2}{\sigma_P^2} \approx \frac{1}{N\sigma_P^2} \sum_{j=1}^{N} \hat{P}_j(\Xi_2) (\hat{P}_j(\Xi_1^i) - \hat{P}_j(\Xi_1)), \quad 1 \le i \le M \qquad (11) $$
wherein $\sigma_P^2 = \text{Var}[\hat{P}]$ is the estimate of the total variance of the random wave field and $\sigma_i^2 = \text{Var}[\mathbb{E}[\hat{P}|X_i]]$ is the partial variance with respect to the random variable $X_i$ or the share of variance of $\hat{P}$ that is due to $X_i$. The second order index $S_{ij}$ describes the sensitivity of the wave field (QoI) to the interactions between the random variables $X_i$ and $X_j$ ($1 \le i \ne j \le M$) without their marginal effects, i.e.

$$S_{ij} = \frac{\sigma_{ij}^2}{\sigma_P^2} = \frac{\text{Var}[\mathbb{E}[\hat{P}|X_i X_j]]}{\text{Var}[\hat{P}]} - S_i - S_j, \quad 1 \le i \ne j \le M \quad (12)$$

The higher order interactions could be defined by similar coefficients which all lie between [0, 1] and the higher values imply higher contribution to the wave field uncertainty (variance). There are $2^M - 1$ Sobol indices that can be defined in which $M$ is the stochastic dimension of the problem. In the numerical examples considered in this paper in Section 4, $M = 6$ in this work which yields 63 Sobol indices. Since large stochastic dimensions imply high number of Sobol indices, Homma and Saltelli [18] introduced the total Sobol indices to decrease the computational costs and make the interpretation easier. These total indices describe the total contribution of any random variable $X_i$ (marginal and joint) to the variance of the random output (wave field) and are defined as:

$$S_{T_i} = S_i + \sum_{i \ne j} S_{ij} + \sum_{(j<k) \ne i} S_{ijk} + \dots, \quad 1 \le i \le M \quad (13)$$

An estimation for the total Sobol index reads:

$$S_{T_i} \approx \frac{1}{2N\sigma_P^2} \sum_{j=1}^{N} (\hat{P}_j(\Xi_1^i) - \hat{P}_j(\Xi_1))^2, \quad 1 \le i \le M \quad (14)$$

### 3.4.2. Analytical sensitivity indices for PCE

For PCE, one can analytically obtain the Sobol sensitivity indices of any order, instead of taking Quasi Monte Carlo samples of the underlying meta-model. The analytical formulas are derived in Sudret [50] in terms of particular sums of squared PCE coefficients. The first-order and total Sobol indices are given as [16, 50]:

$$S_i = \sum_{\alpha \in A_i} a_{\alpha}^2 / \sum_{\alpha \in A, \alpha \ne 0} a_{\alpha}^2 \quad \text{where} \quad A_i = \{\alpha \in A : \alpha_i > 0, \alpha_{j \ne i} = 0\}, \quad 1 \le i \le M, \quad (15)$$
and

$$S_{T_i} = \sum_{\alpha \in \mathcal{A}_i^{\text{tot}}} \frac{a_{\alpha}^2}{\sum_{\substack{\alpha \in \mathcal{A}, \\ \alpha \neq 0}} a_{\alpha}^2} \quad \text{where} \quad \mathcal{A}_i^{\text{tot}} = \{\alpha \in \mathcal{A} : \alpha_i > 0\}, \quad 1 \le i \le M. \qquad (16)$$

# 4. Numerical results

## 4.1. Setup

In this section, the characteristics of the sound velocity profile are modeled by six independent random variables. For the sake of simplicity, the uncertainty is put only on the sound speed profile and the water depth. However, the extension to more complicated cases, for instance when the density of the propagation medium or the properties of the bottom layer are uncertain, is straightforward. In this study, the densities of the water column and the bottom layer are considered as homogeneous and equal to $\rho_1 = 1000 \text{ kg.m}^{-3}$ and $\rho_2 = 1800 \text{ kg.m}^{-3}$, respectively. Table 1 summarizes the range of variation for each uncertain parameter that are taken from [51].

<table><caption>Table 1: Lower and upper bounds for the random parameters considered in this paper.</caption><thead><tr><td>Characteristic of the velocity profile</td><td>min value</td><td>max value</td></tr></thead><tbody><tr><td>Medium depth d [m]</td><td>97.5</td><td>102.5</td></tr><tr><td>First layer depth d<sub>1</sub> [m]</td><td>17.5</td><td>22.5</td></tr><tr><td>Third layer depth d<sub>2</sub> [m]</td><td>57.5</td><td>62.5</td></tr><tr><td>First layer sound speed v<sub>1</sub> [m.s<sup>-1</sup>]</td><td>1535</td><td>1540</td></tr><tr><td>Third layer sound speed v<sub>2</sub> [m.s<sup>-1</sup>]</td><td>1507.5</td><td>1512.5</td></tr><tr><td>Bottom layer sound speed v<sub>3</sub> [m.s<sup>-1</sup>]</td><td>1750</td><td>1850</td></tr></tbody></table>

Since the only available information about these parameters is their lower and upper bounds, based on the maximum entropy principle [20], the best probability density function (PDF) to be used is the uniform distribution. As such, these parameters are assumed to be independent uniform random variables. For instance, the depth of the propagation medium follows the uniform distribution $d \sim U(97.5, 102.5)\text{m}$.

In this paper, two different scenarios are considered and the influence of the random parameters are discussed for each case. These two cases correspond to large and small values of the wavelength to depth ratio $\lambda/d = v/df_0$ which is the key factor in this work. Indeed, $\lambda/d \to \infty$ (or $\lambda \gg d$) implies no reflection at the top and bottom boundaries which means a smoothly varying wave field following $x$ direction. On the contrary, the asymptotic limit $\lambda/d \to 0$ (or $\lambda \ll d$) means high degrees of reflection at the boundaries which implies in particular a highly fluctuating wave field along the path of the wave propagation. To distinguish these cases, the sound waves are generated via a source located at $(x^0, y^0, z^0) = (0, 0, -50)\text{m}$ with two different frequencies of $f_0 = 10 \text{ Hz}$ and $f_0 = 50 \text{ Hz}$. One of the major applications of this range of frequencies is
<figure><img src="image_3.png" /><figcaption>Figure 2: The configuration of the source and the microphones.</figcaption></figure>

in marine seismic exploration where an image of the underlying geology (elastic bottom layer of our model which could also be considered as multiple layers) is constructed by emitting the acoustic energy through the water and then recording the energy that is reflected and refracted by the bottom rock layer(s) back to the receivers. In this case, the conventional acquisition systems work with frequencies under 100 Hz [34]. 401 measurement stations are assumed at (25(*m* − 1), 0, −49) m, *m* = 1, 2, ..., 401, to record the wave field. The length of the medium in *x* direction is thus 10 km. This configuration is depicted in Figure 2. The simulation of sound wave propagation is realized by the wavenumber integration method implemented by a numerical code developed by the authors, its accuracy is justified in [55].

4.2. Results of the quasi Monte Carlo method

Following the distributions given in Table 1, 50000 realizations of the random velocity profile are drawn using a quasi Monte Carlo (QMC) sampling method. The random wave field recorded at the point **r** = (*x*, *y*, *z*) is denoted by *P*(**r**) which is a positive, second-order, mean-squared continuous stochastic field indexed on **R**⁺ and defined on a probability space (Ω, **F**, P). The ensemble average of *P* at each point **r** is thus defined as E[*P*(**r**)] = ∫Ω *P*(**r**)dP. Figure 3 shows the convergence of the normalized mean E[*P*(**r**)]/E[*P*(**0**)] (plots a and d), normalized standard deviation √Var[*P*(**r**)]/√Var[*P*(**0**)] (plots b and e) and the coefficient of variation δ(**r**) = √Var[*P*(**r**)]/E[*P*(**r**)] (plots c and f) of the wave field for the frequencies f₀ = 10 Hz (plots a,b and c) and f₀ = 50 Hz (plots d, e and f). The convergence of these statistics are represented in terms of the number of QMC samples for the wave field recorded at 4 different stations located at
(x ∈ {1, 2500, 5000, 10000}, y = 0, z = -49) m. These figures reveal that one needs large sample sizes (larger
than 10⁴) to guarantee the convergence of the first and second order statistics of the recorded wave field.
This justifies the importance of using surrogate models in which we aim at approximating the probabilistic
content of the random wave field using significantly smaller sample sizes.

<figure><img src="image_4.png" /><figcaption>Figure 3: Convergence of the normalized average (a,d), normalized standard deviation (b,e) and coefficient of variation (c,f) of the wave field recorded at x = 1 m (solid lines), x = 2500 m (dashed lines), x = 5000 km (dotted lines) and x = 10000 m (dash-dot lines) for the frequencies of 10 Hz (a,b,c) and 50 Hz (d,e,f).</figcaption></figure>

4.3. Results of the surrogate models

In the sequel, we use the methods introduced in Section 3 in order to construct meta-models of the
random wave field at each recording point. Kriging, PCE and PCE-Kriging methods are thus used to
estimate the statistics of the wave field and the results are compared to those obtained from the QMC
method. We used in-house Matlab-based scripts to construct our meta-models along with the subsequent
global sensitivity calculations. The readers can use software like UQLab (available at uqlab.com) and
OpenTURNS (openturns.org) in order to reproduce the results.

A sample size of 500 realizations based on Sobol sequences (which are an example of quasi-random
low-discrepancy sequences as explained in [47]) is taken from the six input random variables to build each
surrogate model. Figure 4 shows for both frequencies, the mean, standard deviation and coefficient of
variation of the wave field in terms of the station position. The black, blue, red and green curves correspond
to the QMC (with a sample size of N = 50000), Kriging, PCE and PCE-Kriging methods (each with N =
500), respectively. As it can be observed from the plots (a,b,c), all surrogate models appropriately estimate
the first and second order statistics of the wave field when the source's central frequency is f₀ = 10 Hz. In
this case, the ratio $\lambda/d \approx 1.53$ is large enough so that the influence of the reflections and transmissions at the boundaries on the wave field's fluctuations become very small (see the smooth variation of the mean field in plot (a)). The black curves (QMC results) in plots (c) and (f) show a reduction in the fluctuation amplitudes after some travel distance. This reduction could be justified by the fact that at short travel distances the wave field is reflected by and transmitted into the interface. By contrast, at long travel distances, the energy of the initial waves is attenuated due to the transmissions at the interface and absorption at the level of the PMLs and the wave front is like a wave packet whose shape remains almost unchanged. Figure 4(d, e and f) shows the same results for $f_0 = 50$ Hz in which the ratio $\lambda/d \approx 0.31$ is much lower compared to the previous case so that the interactions of the waves with the boundaries of the medium become substantial and thus the wave field fluctuates more. In this case, deviations from the QMC results become apparent for the second order statistics. These deviations start from the source-station distances above 4500 m and 2500 m for Kriging and PCE methods, respectively. On the contrary, the PCE-Kriging method appears to yield the results closer to the QMC (green curves in the plots (e,f)). Figure 5 depicts, for the both frequencies, the snapshots of the wave field (pressure) near and far away from the source which is located at $\mathbf{r}^0 = (0, 0, -50)m$. Note that these simulations are carried out using the open-source code SPECFEM2D [5] which is based on the spectral element method and is available at geodynamics.org. This software simulates the sound wave propagation in the time domain.

To assess quantitatively the meta-models, Figure 6 shows the LOO error estimates (in percent) in terms of the recording point position for both cases. In the case $f_0 = 10$ Hz, the PCE-Kriging method gives errors less than 1% almost everywhere while the other two methods result in relatively larger errors. For $f_0 = 50$ Hz, Figure 6(b) reveals acceptable errors only within the interval (0, 2000)m. Above 2000 m, the LOO error of the PCE quickly becomes very large and the PCE-Kriging provides overall smaller errors compared to the Kriging method. This drastic increase on the surrogate model errors compared to the case $f_0 = 10$ Hz results from the fact that smaller values of the ratio $\lambda/d$ imply more interactions between the waves and the boundaries of the medium. Hence, in the recorded signals at long propagation distances, the direct waves are highly mixed with the reflected waves. This will result in noticeable fluctuations of the wave field as it could be observed from Figure 4(d) (compared to Figure 4(a)). This irregularity of the wave field at long distances is the main reason that the surrogate models could not yield reasonable errors. However, as far as the first and second order statistics are considered, the Kriging and PCE-Kriging methods approximate well the QMC results.

In order to compare the results more precisely, the PDFs of the wave field normalized by the corresponding
<figure><img src="image_5.png" /><figcaption>(a) $f_0 = 10$ Hz, mean</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) $f_0 = 10$ Hz, standard deviation</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) $f_0 = 10$ Hz, coefficient of variation</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) $f_0 = 50$ Hz, mean</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(e) $f_0 = 50$ Hz, standard deviation</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(f) $f_0 = 50$ Hz, coefficient of variation</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(a) $f_0 = 10$ Hz, near field</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) $f_0 = 10$ Hz, far field</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) $f_0 = 50$ Hz, near field</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) $f_0 = 50$ Hz, far field</figcaption></figure>

<caption>Figure 5: Snapshots of the wave field near (a,c) and far (b,d) from the source for $f_0 = 10$ Hz (a,b) and $f_0 = 50$ Hz (c,d).</caption>

mean value, i.e. $P(r)/E[P(r)]$, are computed at two different distances from the source, $|r-r^0| \in \{1, 5000\}$m and the results are displayed in Figure 7 for both source central frequencies. The black, blue, red and green curves correspond to QMC, Kriging, PCE and PCE-Kriging methods, respectively. This figure reveals that the near-field is almost deterministic (see for instance the Dirac delta-like function in plot b) which becomes more randomized in the far-field. This figure also shows that larger source-station distances imply lower accuracy of the surrogate models, especially at higher source frequencies. For $f_0 = 50$ Hz, the PCE method fails to satisfyingly estimate the wave field. One solution is to employ the experimental designs with significantly larger sizes which is against the main interest of the meta-modeling.
<figure><img src="image_15.png" /><figcaption>Figure 6: LOO errors (in percent) in terms of different surrogate models in terms of the source-station distance for $f_0 = 10$ Hz(a) and $f_0 = 50$ Hz(b).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 7: PDFs of the normalized wave field recorded at two different distances from the source whose central frequency is 10 Hz(a) and 50 Hz(b). Black, blue, red and green curves represent the results obtained by QMC, Kriging, PCE and PCE-Kriging methods, respectively. Solid and dashed dotted lines correspond to $|\mathbf{r} - \mathbf{r}^0| = 1$ m and $|\mathbf{r} - \mathbf{r}^0| = 5000$ m.</figcaption></figure>

### 4.4. Sensitivity analysis

The objective of this section is to identify how the uncertain parameters describing the sound source profile along with the depths, i.e. $d$, $d_1$, $d_2$, $v_1$, $v_2$ and $v_3$ contribute to the variability of the wave field recorded at different distances from the source. Some of the objectives of such a study are: (i) identification of less influential parameters in order to reduce the stochastic dimension of the problem and (ii) reduction of the variance of the recorded wave field (in the forward (direct) problem) by minimizing the variance of the most influential parameters. This variance reduction is crucial for getting better estimations when solving the inverse source localization problem is regarded.

Figures 8 and 9 display for $f_0 = 10$ Hz the variation of the Sobol and total Sobol sensitivity indices in terms of the source-station distance obtained by the QMC (black curves), ordinary Kriging (blue), PCE (red) and PCE-Kriging (green) methods. It turns out that the most influential factor is the velocity of the bottom layer $v_3$. Then, the depth of the shallow water $d$ plays the most important role on the variance of the recorded
pressures. The other 4 random variables ($d_1, d_2, v_1, v_2$) could therefore be considered as uninfluential. Thus, one should try to decrease the uncertainty on $v_3$ and $d$ in order to reduce the uncertainty on the sound pressure. Moreover, the difference between the single and total Sobol indices shows that the interaction terms do not contribute to the output's uncertainty. An important result that could be drawn from these figures is the good accuracy of the PCE method to estimate the Sobol indices. As stated in Section 3.4.2, the Sobol indices could be analytically derived from the PCE coefficients which is not the case for the Kriging and PCE-Kriging methods.

Figures 10 and 11 show the same results for the source's central frequency of $f_0 = 50$ Hz. In this case, the PCE method fails to appropriately estimate the Sobol indices. The reason is that the wave field is such irregular that it could not be approximated by a polynomial function. It can be observed that the PCE-Kriging method returns the best results. Similar to the case where $f_0 = 10$ Hz, the most influential parameters are $v_3$ and $d$. However, the total sensitivity indices indicate that the interaction terms now play an important role on the uncertainties of the parameters $v_1$ and $v_2$.

<figure><img src="image_17.png" /><figcaption>Figure 8: Variation of the Sobol sensitivity indices of all random variables in terms of the source-station distance for $f = 10$ Hz. Black, blue, red and green curves represent the results obtained by QMC, Kriging, PCE and PCE-Kriging methods, respectively.</figcaption></figure>

# 5. Conclusions and perspectives

The robustness of the Kriging, PCE and PCE-based Kriging methods for uncertainty quantification (UQ) in acoustic wave propagation through shallow water environments has been investigated. It is shown
<figure><img src="image_18.png" /><figcaption>Figure 9: Variation of the total Sobol sensitivity indices of all random variables in terms of the source-station distance for $f = 10$ Hz. Black, blue, red and green curves represent the results obtained by QMC, Kriging, PCE and PCE-Kriging methods, respectively.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 10: Variation of the Sobol sensitivity indices of all random variables in terms of the source-station distance for $f = 50$ Hz. Black, blue, red and green curves represent the results obtained by QMC, Kriging, PCE and PCE-Kriging methods, respectively.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 11: Variation of the total Sobol sensitivity indices of all random variables in terms of the source-station distance for $f = 50$ Hz. Black, blue, red and green curves represent the results obtained by QMC, Kriging, PCE and PCE-Kriging methods, respectively.</figcaption></figure>

that the accuracy of these methods in the case of a shallow water with a layered sound velocity profile highly depends on the wavelength-to-depth ratio. Indeed, the latter defines the degree of interaction between the waves and the boundaries. Larger (resp. smaller) values of this ratio yield better (resp. worse) accuracies of the UQ methods in terms of the LOO error. Overall, the PCE-based Kriging has been shown to give better results since it has both the local and global interpolation properties of the Kriging and PCE methods. The variance-based sensitivity analysis has also shown the importance of the depth of the shallow water and the velocity of the bottom layer on the uncertainties of the recorded sound pressure.

In this paper, an underwater wave propagation model with six independent uncertain parameters is considered. The proposed UQ methods can be tested for more complicated cases, e.g. 3D wave propagation in range-dependent environment with irregular sound speed profile, non-flat water bottom and more layer numbers. In these cases, more uncertain parameters with interdependence must be considered, which is more challenging for UQ tools.
Acknowledgments

This work has been carried out thanks to the support of the A*MIDEX grant (reference number: ANR-
11-IDEX-0001-02) funded by the French Government "Investissements d'Avenir" program and Guangxi
Natural Science Foundation (No. 2018GXNSFBA281013).

References

[1] M. Berveiller, B. Sudret, and M. Lemaire. Stochastic finite element: a non intrusive approach by regression. *European Journal of Computational Mechanics/Revue Européenne de Mécanique Numérique*, 15(1-3):81-92, 2006.

[2] G. Blatman and B. Sudret. An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis. *Probabilist. Eng. Mech.*, 25(2):183–197, 2010.

[3] S.-K. Choi, R. V. Grandhi, R. A. Canfield, and C. L. Pettit. Polynomial chaos expansion with Latin hypercube sampling for estimating response variability. *AIAA J.*, 42(6):1191–1198, 2004.

[4] C. A. Clark and K. B. Smith. An efficient normal mode solution to wave propagation prediction. *IEEE J. Ocean. Eng.*, 33(4):462–476, 2008.

[5] P. Cristini and D. Komatitsch. Some illustrative examples of the use of a spectral-element method in ocean acoustics. *J. Acoust. Soc. Am.*, 131(3):EL229-EL235, 2012.

[6] F. R. DiNapoli and R. L. Deavenport. Theoretical and numerical Green's function field solution in a plane multilayered medium. *J. Acoust. Soc. Am.*, 67(3):92-105, 1980.

[7] V. Dubourg, B. Sudret, and J. M. Bourinet. Reliability-based design optimization using kriging surrogates and subset simulation. *Struct. Multidiscip. Optim.*, 44(5):673–690, 2011.

[8] O. Dubrule. Cross validation of kriging in a unique neighborhood. *J. Int. Ass. Math. Geol.*, 15(6):687-699, 1983.

[9] T. A. El-Mihoub, A. A. Hopgood, L. Nolle, and A. Battersby. Hybrid genetic algorithms: A review. *Engineering Letters*, 13(2):124-137, 2006.

[10] S. Finette. Embedding uncertainty into ocean acoustic propagation model. *J. Acoust. Soc. Am.*, 117(3):997-100, 2005.

[11] S. Finette. A stochastic representation of environmental uncertainty and its coupling to acoustic wave propagation. *J. Acoust. Soc. Am.*, 120(5):2567-2579, 2006.

[12] S. Finette. A stochastic response surface formulation of acoustic propagation through an uncertain ocean waveguide environment. *J. Acoust. Soc. Am.*, 126(5):2242-2247, 2009.

[13] F. Gerdes and S. Finette. A stochastic response surface formulation for the description of acoustic propagation through an uncertain internal wave field. *J. Acoust. Soc. Am.*, 132(4):2251-2264, 2012.

[14] R. G. Ghanem and P. D. Spanos. Stochastic finite element method: Response statistics. In *Stochastic finite elements: a spectral approach*, pages 101–119. Springer, 1991.

[15] R. G. Ghanem and P. D. Spanos. *Stochastic finite elements: a spectral approach*. Dover publications, 2003.

[16] L. L. Gratiet, S. Marelli, and B. Sudret. Metamodel-based sensitivity analysis: Polynomial chaos expansions and Gaussian processes. *Handbook of Uncertainty Quantification*, pages 1–37, 2016.

[17] XP. Guo, L. Zhang, D. Wang, Z. Lu, Y. Wang, and J. Zhang. Efficient computation and uncertainty analysis of underwater
acoustic propagation based on kriging surrogate model. *DEStech Transactions on Computer Science and Engineering*, (csae), 2017.

[18] T. Homma and A. Saltelli. Importance measures in global sensitivity analysis of nonlinear models. *Reliab. Eng. Syst. Safety*, 52(1):1–17, 1996.

[19] K. R. James and D. R. Dowling. Pekeris waveguide comparisons of methods for predicting acoustic field amplitude uncertainty caused by a spatially uniform environmental uncertainty (L). *J. Acoust. Soc. Am.*, 129(2):589–592, 2011.

[20] E. T. Jaynes. Information theory and statistical mechanics. *Phys. Rev.*, 106(4):620, 1957.

[21] F. B. Jensen, W. A. Kuperman, M. B. Porter, and H. Schmidt. *Computational Ocean Acoustics, Second Edition*. Springer, 2011.

[22] J. B. Keller. Rays, waves and asymptotics. *Bull. Am. Math. Soc.*, 84(5):727–750, 09 1978. URL http://projecteuclid.org/euclid.bams/1183541137.

[23] S. Khazaie, R. Cottereau, and D. Clouteau. Influence of the spatial correlation structure of an elastic random medium on its scattering properties. *J. Sound Vib.*, 370:132–148, 2016.

[24] S. Khazaie, X. Wang, and P. Sagaut. Localization of random acoustic sources in an inhomogeneous medium. *J. Sound Vib.*, 384:75–93, 2016.

[25] Y. Y. Khine, D. B. Creamer, and S. Finette. Acoustic propagation in an uncertain waveguide environment using stochastic basis expansions. *J. Comput. Acoust.*, 18(04):397–441, 2010.

[26] L. Klimeš. Correlation functions of random media. *Pure Appl. Geophys.*, 159(7-8):1811–1831, 2002.

[27] W. A. Kuperman and J. F. Lynch. Shallow-water acoustics. *Phys Today*, October:55–61, 2004.

[28] C. Lataniotis, S. Marelli, and B. Sudret. Gaussian process modelling using UQLab. *arXiv preprint arXiv:1709.09382*, 2017.

[29] O. Le Maître and O. M. Knio. *Spectral methods for uncertainty quantification: with applications to computational fluid dynamics*. Springer Science & Business Media, 2010.

[30] D. Lee and A. D. Pierce. Parabolic equation development in recent decade. *J. Comput. Acoust.*, 3(02):95–173, 1995.

[31] P. F. J. Lermusiaux, J. Xu, C.-F. Chen, S. Jan, L. Y. Chiu, and Y.-J. Yang. Coupled ocean-acoustic prediction of transmission loss in a continental shelfbreak region: Predictive skill, uncertainty quantification, and dynamical sensitivities. *IEEE J. Ocean. Eng.*, 35(4):895–916, 2010.

[32] S. J. Levinson, E. K. Westwood, R. A. Koch, S. K. Mitchell, and C. V. Sheppard. An efficient and robust method for underwater acoustic normal-mode computations. *J. Acoust. Soc. Am.*, 97(3):1576–1585, 1995.

[33] W. Y. Luo, X. L. Yu, X. F. Yang, and R. H. Zhang. Analytical solution based on the wavenumber integration method for the acoustic field in a Pekeris waveguide. *Chin. Phys. B*, 25(4):044302, 2016.

[34] A. Mànuel, X. Roset, J. Del Rio, D. M. Toma, N. Carreras, S. Shariat Panahi, A. García-Benadí, T. Owen, and J. Cadena. Ocean bottom seismometer: design and test of a measurement system for marine seismology. *Sensors*, 12(3):3693–3719, 2012.

[35] A. Marrel, B. Iooss, F. Van Dorpe, and E. Volkova. An efficient methodology for modeling complex computer codes with Gaussian processes. *Comput. Stat. Data Anal.*, 52(10):4731–4744, 2008.

[36] A. Mesgouez, S. Buis, G. Lefeuve-Mesgouez, and G. Micolau. Use of global sensitivity analysis to assess the soil poroelastic parameter influence. *Wave Motion*, 72:377–394, 2017.

[37] Y. Mi, H. Zheng, and H. P. Lee. A domain decomposition method for stochastic analysis of acoustic fields with hybrid
and localized uncertainties. *Wave Motion*, 83:121-133, 2018.

[38] C. L. Pekeris. Theory of propagation of explosive sound in shallow water. *Geol Soc Am Mem*, 27:1-116, 1948.

[39] S. Rahman. Wiener-Hermite polynomial expansion for multivariate Gaussian probability measures. *J. Math. Anal. Appl.*, 454(1):303-334, 2017.

[40] C. E. Rasmussen. Gaussian processes in machine learning. In *Advanced lectures on machine learning*, pages 63-71. Springer, 2004.

[41] J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn. Design and analysis of computer experiments. *Statistical science*, pages 409-423, 1989.

[42] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and S. Tarantola. Variance-based sensitivity analysis of model output. Design and estimator for the total sensitivity index. *Comput. Phys. Commun.*, 181(2):259-270, 2010.

[43] T. J. Santner, B. J. Williams, and W. I. Notz. *The design and analysis of computer experiments*. Springer Science & Business Media, 2013.

[44] E. Savin and B. Faverjon. Computation of higher-order moments of generalized polynomial chaos expansions. *Int. J. Numer. Meth. Eng.*, 111(12):1192-1200, 2017.

[45] R. Schöbi, B. Sudret, and J. Wiart. Polynomial-chaos-based kriging. *Int. J. Uncertain. Quantif.*, 5(2), 2015.

[46] R. Schöbi, B. Sudret, and S. Marelli. Rare event estimation using polynomial-chaos kriging. *ASCE ASME J. Risk Uncertain. Eng. Syst. A Civ. Eng.*, 3(2):D4016002, 2016.

[47] I. M. Sobol'. On the distribution of points in a cube and the approximate evaluation of integrals. *Zhurnal Vychislitel'noi Matematiki i Matematicheskoi Fiziki*, 7(4):784-802, 1967.

[48] I. M. Sobol'. Sensitivity estimates for non-linear mathematical models. *Mathematical Modeling and Computational Experiment*, 1(4):407-414, 1993.

[49] M. L. Stein. *Interpolation of spatial data: some theory for kriging*. Springer Science & Business Media, 2012.

[50] B. Sudret. Global sensitivity analysis using polynomial chaos expansions. *Reliab. Eng. Syst. Safety*, 93(7):964-979, 2008.

[51] J. Tabrikian and J. L. Krolik. Robust maximum-likelihood source localization in an uncertain shallow-water waveguide. *J. Acoust. Soc. Am.*, 101(1):241-249, 1997.

[52] X. Wang, S. Khazaie, and P. Sagaut. Sound source localization in a randomly inhomogeneous medium using matched statistical moment method. *J. Acoust. Soc. Am.*, 138(6):3896-3906, 2015.

[53] X. Wang, S. Khazaie, L. Margheri, and P. Sagaut. Shallow water sound source localization using the iterative beamforming method in an image framework. *J. Sound Vib.*, 395:354-370, 2017.

[54] X. Wang, S. Khazaie, and X. Chen. Linear approximation of underwater sound speed profile: Precision analysis in direct and inverse problems. *Appl. Acoust.*, 140:63-73, 2018.

[55] X. Wang, S. Khazaie, D. Komatitsch, and P. Sagaut. Sound-source localization in range-dependent shallow-water environments using a four-layer model. *IEEE J. Ocean. Eng.*, 44 (1):220-228, 2019.