Resilient Coastal Protection Infrastructures: Probabilistic Sensitivity Analysis of Wave Overtopping Using Gaussian Process Surrogate Models

Abstract

This paper presents a novel mathematical framework for assessing and predicting the resilience of critical coastal infrastructures against wave overtopping hazards and extreme climatic events. A probabilistic sensitivity analysis model is developed to evaluate the relative influence of hydrodynamic, geomorphological, and structural factors contributing to wave overtopping dynamics. Additionally, a stochastic Gaussian process (GP) model is introduced to predict the mean overtopping discharge from coastal defences. Both the sensitivity analysis and the predictive models are validated using a large homogeneous dataset comprising 163 laboratory and field-scale tests. Statistical evaluations demonstrate the superior performance of the GPs in identifying key parameters driving wave overtopping and predicting mean discharge rates, outperforming existing regression-based formulae. The proposed model offers a robust predictive tool for assessing the performance of critical coastal protection infrastructures under various climate scenarios.

1. Introduction

Coastal regions are densely populated, with over 2.4 billion people (40% of the world’s population) living within 100 km of the coast [1]. More than 600 million people (10% of the world’s population) live in coastal areas less than 10 meters above sea level, making them highly vulnerable to flooding and inundation during extreme climatic events. Over the past century, global sea levels have risen at an accelerating rate, with projections suggesting a rise of over a meter by 2100. Such increases contribute to a host of adverse outcomes, including intensified erosion, storm surge flooding, inundation, contamination of freshwater, and the loss of coastal lowlands and wetlands [2,3]. The rise in sea level and associated risks affect two-thirds of the world’s 23 most densely populated cities [4]. In recent years, coastal flooding triggered by extreme events has had devastating socioeconomic impacts. The scientific prediction for a high-emission climate change scenario (+2 °C temperature rise) estimates an annual $14 trillion cost of flooding worldwide [5], with significant damage to critical infrastructures.

Coastal defences play a vital role in protecting coasts from flooding and erosion. The effective and optimal design of coastal protection structures are important for safeguarding coastal communities from storm surges and extreme climatic events [6]. In the past decades, several types of flood protection schemes, including hard-engineered, soft nature-based, and hybrid solutions, have been developed and tested [7]. Hard defences, in particular, vary widely in terms of shape, structural design, and materials. Extensive research has been conducted to evaluate the performance of these defences under changing hydrodynamic, geomorphological, and structural conditions. Wave overtopping is a crucial design parameter for coastal defences as it directly affects both the structural integrity of the defence and the level of protection it provides. Given the direct link between wave overtopping and the vulnerability of coastal areas to flooding, it is vital to develop robust methods for evaluating and predicting wave overtopping under different hydrodynamic and structural conditions. However, accurately predicting the performance of critical coastal infrastructure remains a significant challenge. This difficulty stems from the complexity of nearshore processes, which operate at varying spatial and temporal scales. It is further compounded by the non-linear interactions between hydrodynamics and geomorphology, the complex geometry of nearshore areas and coastal defences, and the uncertainties associated with the incident wave climate. Furthermore, a lack of comprehensive data and field observations exacerbates these challenges [8,9].

The interaction between nearshore wave dynamics and defence structures leads to complex hydrodynamic responses, which are influenced by factors such as water depth h, incident wave height H, and period T. Previous studies undertook extensive laboratory physical modelling and field measurements to investigate the effects of wave structure interactions on the structural response to wave overtopping. The majority of these investigations led to empirical formulae for overtopping; however, these formulae are typically valid only within specific hydrodynamic ranges or for particular structural configurations. In recent decades, efforts have been made to develop integrated tide-surge and wave flood models, often using RANS-VOF models and shallow water equations [10,11]. Despite their utility, these models are limited by their depth- and time-averaged nature of flow solvers and the inherent simplifications in turbulence closure models. More recently, Lagrangian particle-based models have been employed to simulate wave interactions with coastal infrastructures and overtopping processes ([12,13]). These models are advantageous because they account for depth-varying hydrodynamics in wave–structure interactions. However, they are often case-specific, computationally expensive, and time-consuming. To address the challenges associated with understanding wave overtopping processes from laboratory and numerical models, machine learning approaches have gained traction as a reliable method for predicting wave run-up and overtopping at coastal defences [6,14]. Machine learning techniques have shown superior predictive capabilities compared with traditional overtopping formulae based on a regression analysis of physical model data. Nonetheless, existing predictive tools for wave overtopping are generally limited by their focus on specific geometries and hydrodynamic conditions. There remains a critical need to develop a more comprehensive predictive framework capable of providing robust overtopping estimates across a wide range of defence structures and hydro-meteorological conditions.

We now explore the complex processes driving overtopping. A thorough understanding of these processes will provide insights for selecting input parameters to develop an efficient model and improving predictive capabilities.

2. Wave Overtopping Processes

Overtopping is a complex, multi-faceted phenomenon influenced by the interaction between nearshore wave processes, seabed topography, and the configuration of coastal defence structures. The complicated interactions between the incident wave and coastal defences can generate overtopping flows that are sudden and do not follow the typical wave up-rush (run-up) processes. To fully understand wave–structure interactions, it is necessary to study the nearshore wave processes. Waves approaching coastal protection structures can be classified into ‘breaking’ or ‘non-breaking’ waves. Extensive research has been conducted to explore the distinct interactions between these wave types and coastal defences (e.g., see [15,16]). Although these two terms may not be entirely precise, they are widely used in the literature to characterise wave overtopping behaviours. This section provides an overview of the key concepts related to wave overtopping.

Wave breaking on beaches and gently sloping structures is commonly characterised by the Iribarren number, also known as the surf similarity parameter (shown in Equation (1)). The Iribarren number defines four distinct breaking regimes, including spilling (${\xi }_{op}<0.4$), plunging ($0.4<{\xi }_{op}<2.3$), collapsing ($2.3<{\xi }_{op}<3.2$), and surging (${\xi }_{op}>3.2$). These regimes are frequently used in design manuals [17] to assess armour stability for sloped structures.

$${\xi }_{op}={\displaystyle \frac{tan\left(\alpha \right)}{{\left(S_{op}\right)}^{0.5}}}$$

(1)

where $\alpha$ is the beach slope and $S_{op}$ is the wave steepness.

For vertical walls and steep structures, the impulsive and non-impulsive breaking is defined by the slope and the incident wave’s wavelength. Besley et al. [18] introduced the wave breaking parameter, $h_{*}$, which is based on the depth at the toe of the wall ($h_s$) and nearshore incident wave conditions (Equation (2)).

$$h_{*}={\displaystyle \frac{h_s}{H_{si}}}\left({\displaystyle \frac{2\pi h_s}{g{T}_m^2}}\right)$$

(2)

where $H_{si}$ is the inshore significant wave height and $T_m$ is the averaged wave period determined from spectral moments or a zero-crossing analysis. Besley et al. [18] suggested that impulsive wave conditions occur at the wall when $h_{*}\le 0.3$ and that pulsating conditions occur when $h_{*}>0.3$.

Wave overtopping rates must remain below tolerable thresholds under design and operational conditions to ensure the safety of people and property on or behind coastal defence structures [19,20,21]. The volume of overtopping is significantly influenced by the nature of wave interactions with the defence structure. When waves break onto or over the structure, overtopping tends to produce relatively continuous volumes of water, known as ‘green water’. In contrast, when waves break seaward of the structure, overtopping occurs in the form of fine droplets—a phenomenon referred to as ‘splash overtopping’. This type of overtopping is typically carried by the wave’s momentum or enhanced by onshore winds [22]. Onshore winds can play a critical role in intensifying overtopping, especially when reflected waves (e.g., from steep walls) interact with incoming waves, creating local clapotis effects. Although onshore wind significantly influences spray overtopping, particularly from vertical walls, limited data are quantifying its impact. Research indicates that while onshore wind has minimal effect on ‘green overtopping’, it can amplify spray overtopping from vertical walls by a factor of up to three [23].

3. Database

The overtopping data utilised in this study are derived from both two-dimensional and three-dimensional physical modelling tests, as well as full-scale field-based (prototype) measurements reported in the CLASH dataset [24]. This comprehensive dataset comprises 163 test series with a total of 10,532 data points and 31 parameters. The dataset includes overtopping records for a range of hydrodynamics and structural configurations, covering both small-scale and large-scale laboratory studies, as well as field measurements.

Table 1. Summary of dataset parameters used in this study and the experimental range.

Parameter	Range	Parameter	Range	Description
Structural parameters		Hydrodynamic parameters
$h_{deep}\left[m\right]$	(0, 100)	$H_{m0deep}\left[m\right]$	(0.003, 5.920)	$=4\sqrt{m_0}$
$m[-]$	(6, 1000)	$T_{pdeep}\left[s\right]$	[0.545, 15]
$h\left[m\right]$	(0.029, 9.32)	$T_{mdeep}\left[s\right]$	(0.454, 12.5)	= $m_2/m_0$
$h_t\left[m\right]$	(0.025, 7.78)	$T_{m-1,0}deep\left[s\right]$	(0.495, 13.636)	= $m_{-1}/m_0$
$B_t\left[m\right]$	(0, 10)	$\beta {[}^o]$	(0, 80)
${\gamma }_f[-]$	(0.35, 1)	$H_{m0toe}\left[m\right]$	(0.003, 3.8)	=4$\sqrt{m_0}$
$cot{\alpha }_d[-]$	(0, 7)	$T_{ptoe}\left[s\right]$	(0.545, 16.4)
$cot{\alpha }_d[-]$	(−5, 9.706)	$T_{mtoe}\left[s\right]$	(0.454, 11.881)	$=m_2/m_0$
$cot{\alpha }_{excl}[-]$	(−1.533, 8.144)	$T_{m-1,0toe}\left[s\right]$	(0.495, 10.64)	$=m_{-1}/m_0$
$cot{\alpha }_{incl}[-]$	(−1.533, 12.821)	$q[m^3/s.m]$	(0, $1.65\times {10}^{-1}$)
$R_c\left[m\right]$	(0, 8.345)	$P_{ow}[-]$	(0, 81)
$B\left[m\right]$	(0, 8)
$h_b\left[m\right]$	(−0.208, 1.175)	General parameters
$tan{\alpha }_B[-]$	(0, 0.125)	$RF[-]$	(1, 4)
$B_h\left[m\right]$	(0, 8)	$CF[-]$	(1, 4)
$A_c\left[m\right]$	(0, 7.87)
$G_c\left[m\right]$	(0, 5.6)

summarises the parameters recorded in the dataset. The homogeneous database includes measurements from various types of hard-engineered coastal protection structures, including vertical structures, rubble mound breakwaters (with rock or concrete armour), dikes, berm breakwaters, and composite structures. Notably, the dataset also features tests with zero overtopping, making it particularly valuable for developing predictive models. De Rouck et al. [25] compared the dataset to empirical overtopping prediction formulae for vertical structures [22], sloping structures (${\gamma }_f=0.5$, as suggested in [26]), and dikes [26]. This comparison aimed to identify outliers within the dataset. A general form of the overtopping formulae proposed in the EurOtop manual is presented in Equation (3) [27].

$${\displaystyle \frac{q}{{\left(g{H}_{m0}^3\right)}^{0.5}}}=Aexp\left(-{\displaystyle \frac{B{R}_c}{\gamma H_{m0}}}\right)$$

(3)

where $H_{m0}$ denotes the significant wave height based on spectral analysis, $R_c$ is the structure crest freeboard relative to still water line, $\gamma$ is a correction factor for the roughness and angle of wave approach or structure geometry, and A and B are constants (fitting coefficients) that are determined for different types of structures based on a regression analysis.

All data used in this study are classified with a reliability factor (RF) and complexity factor (CF), based on the screening study conducted in [28]. Table 1 presents the range of RF and CF for the CLASH dataset. Tests with an RF/CF of four are deemed unreliable and are excluded from the modelling process in this study. Figure 1 shows a schematic of the hydrodynamic and structural parameters considered in this study for the probabilistic sensitivity analysis and for developing a predictive model. To present a holistic overview of the data used for this study and the range they cover, Figure 2 plots dimensionless wave overtopping rates measured against a relative crest freeboard. The dataset effectively covers the range of ${10}^{-6}\le q\le {10}^{-1}$ and $0.3\le R_c/H_{m0}\le 3.5$. Notably, outliers are present, manifested as large, unexpected overtopping values. These outliers are primarily attributed to experimental conditions, such as very small wave steepness and shallow foreshores. Additionally, the unrealistically low wave overtopping discharges observed for cases with small relative freeboards are associated with tests involving high and wide-crested rubble mound armour.

The range of wave steepness as a function of the wave height was determined for all the data and is presented in Figure 3. Wave steepnesses ($S_{op}$) exceeding 0.07 are physically unrealistic, while values lower than 0.005 are challenging to generate reliably in controlled conditions. Therefore, tests with $S_{op}>0.07$ or $S_{op}<0.005$ are considered unreliable and are excluded from further analysis in this study.

Figure 4 illustrates the combination of upper and down slopes for structures with and without berms. The majority of the tests involve uniform sloping structures. A negative upper slope corresponds to structures with a large return wall, while data with $cot\left({\alpha }_{up}\right)=0$ represent structures with a vertical upper section and a sloping lower section. Data points along the vertical axis correspond to structures with a vertical down slope and sloping upper section. Further analysis of the dataset reveals that most data points fall within the range of $0<\mathrm{width}/H_{m0,toe}<10$ and $-5<\mathrm{level}/H_{m0,toe}<5$.

4. Method

This section outlines the use of Gaussian process regression (GPR) to probabilistically model the complex relationships between wave overtopping (q) and the input parameters listed in Table 1. Before introducing the GPR approach, the variance-based and emulator-based methods are described as efficient probabilistic sensitivity analysis (SA) techniques. These methods are used to determine the relative importance and influence of each input parameter on wave overtopping. By building a model that focuses on the most influential input parameters, we aim to reduce model complexity, prevent overfitting, and improve predictive accuracy.

4.1. Probabilistic Sensitivity Analysis

In this paper, a global SA of the model output is performed, which evaluates the relative importance of input parameters when they are varied extensively, accounting for their uncertainties over a broad range. One approach to global SA is the analysis of variance of the model response originally proposed by [29]. This method identifies the contribution of individual inputs or groups of inputs to the overall variance in the model’s output. It also assesses the total effect of each input on output variance, including both its marginal influence and its interaction with other inputs. Several computational techniques can be applied to carry out this SA, as detailed in [29,30,31]. This study adopts the emulator-based method outlined by [32] to determine the sensitivity measures.

To conduct the SA, we consider how a function $f\left(\mathbf{x}\right)$ depends on its input variables. In this study, f typically represents the function that computes wave overtopping as a function of a vector of input parameters illustrated in Table 1. Some key notations are introduced below. We define a d-dimensional random vector as $\mathbf{X}=(X_1,\dots ,X_d)$, where $X_i$ denotes the ith element of $\mathbf{X}$. The sub-vector $(X_i,X_j)$ is represented as ${\mathbf{X}}_{i,j}$. More generally, if p represents a set of indices, then ${\mathbf{X}}_p$ indicates the sub-vector of $\mathbf{X}$ comprising elements with those indices. ${\mathbf{X}}_{-i}$ is defined as the sub-vector of $\mathbf{X}$ containing all elements except $x_i$. Similarly, $\mathbf{x}=(x_1,\dots ,x_d)$ represents the corresponding observed random vector $\mathbf{X}$. In this context, $\mathbf{X}$ serves as an input vector consisting of all input parameters outlined in Table 1, while q (wave overtopping) is regarded as the output variable and is denoted by Y.

4.1.1. Function Decomposition for Main Effects and Interactions

Sobol [29] demonstrates that any function $f(·)$ with quadratic integrability can be expressed through its main effects (MEs) and interactions, as follows:

$$y=f\left(\mathbf{x}\right)=z_0+{\Sigma }_{i=1}^d{z}_i\left(x_i\right)+{\Sigma }_{i<j}{z}_{i,j}\left({\mathbf{x}}_{i,j}\right)+\dots +z_{1,2,\dots ,d}\left(\mathbf{x}\right).$$

(4)

In the context of a relationship between y and $\mathbf{x}$, expressed by $y=f\left(\mathbf{x}\right)$ ( where $f(.)$ represents a function involving uncertain quantities $\mathbf{x}$), the expected value of which is represented by $z_0=E\left[f\left(\mathbf{X}\right)\right]$. The function $z_i\left(x_i\right)$, as seen in Equation (5), stands for the “main effect” of the ith variable, $x_i$, and is expressed as:

$$z_i\left(x_i\right)=E[f\left(\mathbf{X}\right)\mid x_i]-E\left[f\left(\mathbf{X}\right)\right].$$

(5)

The ME, $z_i\left(x_i\right)$, is the function solely dependent on $x_i$ that provides the optimal approximation of $f(.)$ by minimising the variance when averaged over all other variables [32,33].

The function $z_{i,j}\left({\mathbf{x}}_{i,j}\right)$, as presented in Equation (6), characterises the first-order interaction between variables $x_i$ and $x_j$. Likewise, $z_{i,j,k}\left({\mathbf{x}}_{i,j,k}\right)$ represents the second-order interaction among $x_i$, $x_j$, and $x_k$, with this pattern extending to higher-order interactions.

$$z_{i,j}\left({\mathbf{x}}_{i,j}\right)=E[f\left(\mathbf{X}\right)\mid {\mathbf{x}}_{i,j}]-z_i\left(x_i\right)-z_j\left(x_j\right)-E\left[f\left(\mathbf{X}\right)\right].$$

(6)

In the context of this study, the main effect (ME) represents the expected change in wave overtopping when input parameter i takes a specific value $x_i$, while the uncertainty in the remaining parameters is still accounted for.

It is common for the functions in Sobol’s decomposition to be pairwise orthogonal. The definitions of the main effects and interaction terms in Sobol’s decomposition, as given in Equation (4), depend on the distribution of the input parameters, $\mathbf{X}$ denoted by G (further details are discussed in Section 4.2.2).

The sensitivity metrics evaluated in this study help identify which input parameters in $\mathbf{X}$ most significantly contribute to the uncertainty in $f(.)$. Examining the MEs and, where applicable, first-order interaction terms provides valuable insights. Their corresponding plots offer a useful visual tool for understanding how individual inputs affect the model’s output and how these inputs interact to influence the overall behaviour of the model.

4.1.2. Variance-Based Methods

Variance-based methods assess the sensitivity of the output, in this case, the overtopping volume, $Y=f\left(\mathbf{X}\right)$, by analysing how changes in model input parameters affect the variance of Y. A comprehensive review of this approach can be found in [30]. Two key sensitivity measures for the model output Y with respect to an individual input $x_i$ are introduced. The variance of the ME indicates the potential reduction in the overall variance of $f(·)$ if $x_i$’s were known. The first of these measures is defined in Equation (5).

$$V_i=var\left\{E(Y\mid X_i)\right\}.$$

(7)

The second variance-based SA measure, proposed by [34], can be written as:

$$V_{T_i}=var\left(Y\right)-var\left\{E(Y\mid {\mathbf{X}}_{-i})\right\}$$

(8)

which is the remaining uncertainty in Y that is unexplained after everything has been learnt except $x_i$.

These two measures (Equations (6) and (7)) can be converted into scale-invariant measures by dividing by $var\left(Y\right)$, as follows:

$$S_i={\displaystyle \frac{V_i}{var\left(Y\right)}},S_{T_i}={\displaystyle \frac{V_{T_i}}{var\left(Y\right)}}=1-S_{-i}$$

(9)

where $S_i$ is the ME index of $x_i$ and $S_{T_i}$ is the total effect index of $x_i$.

The variance measures are linked to the Sobol decomposition when the parameters are independent. The total variance of f can be represented as the sum of the variances for each term as given in Equation (7) (see [32,33] for further details on the Sobol decomposition).

4.2. Emulator-Based Sensitivity Analysis

Theoretically, if the function $f\left(\mathbf{x}\right)$ was not extremely complex, the sensitivity measures discussed earlier could be computed analytically. However, given the complexity of the models in this study, these measures cannot be evaluated analytically. Instead, a computationally efficient and robust model is required to compute the sensitivity measures discussed in Section 4.1.

When $f\left(\mathbf{x}\right)$ is computationally inexpensive and can be quickly computed and easily evaluated for many different inputs, standard Monte Carlo (MC) methods are sufficient for estimating $var\left(Y\right)$ and other sensitivity measures introduced in Section 4.1. However, the computation techniques proposed by [29,30] require many thousands of function evaluations, making them impractical for more expensive functions.To address this computational complexity, the promising methodology outlined in [32], based on the Bayesian paradigm, enables the estimation of all necessary quantities for sensitivity analysis in modeling and predicting wave overtopping for coastal defenses.

Because the functional relationship $f(.)$ between wave overtopping and the input parameters given in Table 1 is unknown for any specific input configuration $\mathbf{x}$ until the model is actually run for those inputs, we need to make some assumptions. Within a Bayesian framework, it is appropriate to define a prior distribution for the values of $f\left(\mathbf{x}\right)$ at various $\mathbf{x}$ points. This prior distribution is then updated using Bayesian methods, with data $\mathcal{D}=\{({\mathbf{x}}_i,y_i):y_i=f\left({\mathbf{x}}_i\right),i=1,\dots ,n\}$ being generated from a sequence of model simulations. The outcome is a posterior distribution for $f(·)$, enabling formal Bayesian inferences regarding the sensitivity analysis measures described earlier.

While uncertainty persists regarding the function $f(·)$ at input or parameter values where it has not been evaluated, considering the correlation between function values at different points can help reduce this uncertainty. Typically, the expected value of the posterior distribution serves as a point estimate for $f(·)$. In SA, two distinct distributions are utilised: first, the distribution G, representing uncertainty in model inputs/parameters $\mathbf{x}$, which is propagated to output values via $f(·)$; and second, the posterior distribution on $f(·)$, serving a purely computational role. This latter can be refined by increasing the number of training points $\mathbf{x}$, but it has no operational interpretation.

4.2.1. Gaussian Process Emulators

Consider the function or complex model under investigation as a deterministic code that yields an output $y=f\left(\mathbf{x}\right)$ for a given input vector $\mathbf{x}$. A GP emulator offers a prior representation of our uncertain understanding of this function’s values before the deterministic code is executed. Evaluating the code at various input configurations generates the data required to construct a posterior distribution. The process for defining the prior distribution for $f(·)$ involves two stages: initially, a GP is chosen as the prior distribution based on specific hyperparameters. Subsequently, these hyperparameters are estimated to define the mean vector and covariance matrix of the GP model. The key requirement for employing the GP is that $f(·)$ should be a smooth function. Thus, knowing the value of $f\left(\mathbf{x}\right)$ should provide insight into the value of $f\left({\mathbf{x}}^{\prime }\right)$ for $\mathbf{x}$ near ${\mathbf{x}}^{\prime }$. This smoothness assumption gives GPs a significant computational advantage over Monte Carlo (MC) methods, which often ignore the expected similarity of function values at proximate points. In this study, we select high-reliability points (see Section 5.1) to minimise noise in our model and select a subset of data that only includes overtopping measurements in the presence of a berm, thereby eliminating a distinct class of events that could violate the continuity assumption inherent in GPs. One of the key strengths of GPs is their ability to provide uncertainty quantification in predictions, even when the smoothness assumption is partially violated. While incorporating non-stationary kernels could allow us to handle potential discontinuities more explicitly, such an approach lies beyond the scope of this study. Instead, we rely on careful data selection to maintain model validity and enhance prediction robustness.

Using a GP prior for $f(·)$ means that the uncertainty about $f\left({\mathbf{x}}_1\right),\dots ,f\left({\mathbf{x}}_n\right)$, given any set of points ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n$, can be expressed as a multivariate normal distribution. Consequently, we need to make feasible prior assumptions about the mean and covariance. The mean of $f\left(\mathbf{x}\right)$, given the hyperparameters $\beta$, is modelled as

$$E\left[f\left(\mathbf{x}\right)\right|\mathit{\beta }]=\mathbf{h}{\left(\mathbf{x}\right)}^T\mathit{\beta }$$

(10)

where $\mathbf{h}(·)$ is a vector of q known functions of $\mathbf{x}$ and $\beta$ is a vector of coefficients. The selection of $\mathbf{h}(·)$ is flexible, but it should be made to incorporate any prior beliefs we have about the form of $f(·)$. The covariance between $f\left(\mathbf{x}\right)$ and $f\left({\mathbf{x}}^{\prime }\right)$ is given by:

$$cov(f\left(\mathbf{x}\right),f\left({\mathbf{x}}^{\prime }\right)|{\sigma }^2)={\sigma }^{2}c(\mathbf{x},{\mathbf{x}}^{\prime })$$

(11)

where $c(·,·)$ is a monotone correlation function on ${\mathbb{R}}^{+}$ with $c(\mathbf{x},\mathbf{x})=1$, and it decreases as $|\mathbf{x}-{\mathbf{x}}^{\prime }|$ increases. Furthermore, the function $c(·,·)$ must ensure that the covariance matrix of any set of outputs $\{y_1=f\left({\mathbf{x}}_1\right),\dots ,y_n=f\left({\mathbf{x}}_n\right)\}$ is positive semi-definite. Throughout this paper, we use the following correlation function, which satisfies all the conditions mentioned above and is widely used for its computational convenience:

$$c(\mathbf{x},{\mathbf{x}}^{\prime })=exp\{-{(\mathbf{x}-{\mathbf{x}}^{\prime })}^T\mathbf{B}(\mathbf{x}-{\mathbf{x}}^{\prime })\},$$

(12)

where $\mathbf{B}$ is a diagonal matrix composed of positive smoothness parameters ${\left\{{\left(\sqrt{2}{b}_i\right)}^{-2}\right\}}_{i=1}^d$, where d represents the dimension of $\mathbf{x}$. It should be noted that $\mathbf{B}$ functions to re-scale the distance between $\mathbf{x}$ and ${\mathbf{x}}^{\prime }$, thereby determining the proximity required between two inputs $\mathbf{x}$ and ${\mathbf{x}}^{\prime }$ for the correlation between $f\left(\mathbf{x}\right)$ and $f\left({\mathbf{x}}^{\prime }\right)$ to achieve a specific value.

The normal inverse gamma distribution for $(\mathit{\beta },{\sigma }^2)$ was proposed in [32] for fixed hyperparameters $\mathbf{z},V,a$ and d, as represented by

$$p(\mathit{\beta },{\sigma }^2)\propto {\left({\sigma }^2\right)}^{-{\displaystyle \frac{1}{2}}(d+q+2)}exp\{-\{{(\mathit{\beta }-\mathbf{z})}^T{V}^{-1}(\mathit{\beta }-\mathbf{z})+a\}/\left(2{\sigma }^2\right)\}$$

The function $f(·)$ yields outputs at n predetermined design points ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n$, generating the dataset as $\mathbf{y}=\{f\left({\mathbf{x}}_1\right),\dots ,f\left({\mathbf{x}}_n\right)\}$. Unlike MC methods, these points are intentionally chosen to provide informative insights about $f(·)$. Typically, these design points are strategically distributed across the input space $\mathcal{X}$ of $\mathbf{X}$, the unknown variables, guided by the probability distribution $G\left(\mathbf{X}\right)$. Consequently, the selection of design points is influenced by $G(·)$, as detailed in [35]. The standardised posterior distribution of $f(·)$, conditioned on $\mathbf{y}=\{f\left({\mathbf{x}}_1\right),\dots ,f\left({\mathbf{x}}_n\right)\}$, is then determined as:

$${\displaystyle \frac{f\left(\mathbf{x}\right)-m^{*}\left(\mathbf{x}\right)}{\widehat{\sigma }\sqrt{c^{*}(\mathbf{x},{\mathbf{x}}^{\prime })}}}\mid \mathbf{y}\sim t_{d+n},$$

(13)

where $t_{d+n}$ stands for a Student’s t-distribution with $n+d$ degrees of freedom.

The mean of the resulting posterior is then given by

$$m^{*}\left(\mathbf{x}\right)=\mathbf{h}{\left(\mathbf{x}\right)}^T\widehat{\mathit{\beta }}+\mathbf{t}{\left(\mathbf{x}\right)}^T{A}^{-1}(\mathbf{y}-H\widehat{\mathit{\beta }}),$$

(14)

the updated correlation function described in Equation (12) given the observed data can be written as:

$$\begin{array}{cc}\hfill c^{*}(\mathbf{x},{\mathbf{x}}^{\prime })& =c(\mathbf{x},{\mathbf{x}}^{\prime })-\mathbf{t}{\left(\mathbf{x}\right)}^T{A}^{-1}\mathbf{t}\left({\mathbf{x}}^{\prime }\right)+(\mathbf{h}{\left(\mathbf{x}\right)}^T\hfill \\ \hfill & -\mathbf{t}{\left(\mathbf{x}\right)}^T{A}^{-1}H){\left(H^T{A}^{-1}H\right)}^{-1}{(h{\left(x^{\prime }\right)}^T-\mathbf{t}{\left({\mathbf{x}}^{\prime }\right)}^T{A}^{-1}H)}^T\hfill \end{array}$$

(15)

and

$$\begin{array}{c}\hfill \mathbf{t}{\left(\mathbf{x}\right)}^T=(c(\mathbf{x},{\mathbf{x}}_1),\dots ,c(\mathbf{x},{\mathbf{x}}_n)),\\ \hfill H^T=({\mathbf{h}}^T{\left({\mathbf{x}}_1\right)}^T,\dots ,{\mathbf{h}}^T{\left({\mathbf{x}}_n\right)}^T),\end{array}$$

(16)

$$A=\left(\begin{array}{cccc}1& c({\mathbf{x}}_1,{\mathbf{x}}_2)& \dots & c({\mathbf{x}}_1,{\mathbf{x}}_n)\\ c({\mathbf{x}}_2,{\mathbf{x}}_1)& 1& & ⋮\\ ⋮& & \ddots & \\ c({\mathbf{x}}_n,{\mathbf{x}}_1)& \dots & & 1\end{array}\right)$$

(17)

$$\begin{array}{cc}\hfill \mathit{\beta }& =V^{*}(V^{-1}\mathbf{z}+H^T{A}^{-1}\mathbf{y}),\hfill \\ \hfill {\widehat{\sigma }}^2& ={\displaystyle \frac{\{a+{\mathbf{z}}^T{V}^{-1}\mathbf{z}+{\mathbf{y}}^T{A}^{-1}\mathbf{y}-{\widehat{\mathit{\beta }}}^T{\left(V^{*}\right)}^{-1}\widehat{\mathit{\beta }}\}}{(n+d-2)}}\hfill \\ \hfill V^{*}& ={(V^{-1}+H^T{A}^{-1}H)}^{-1}.\hfill \end{array}$$

The outputs for any set of inputs will follow a multivariate t-distribution, with the covariance between any two outputs being defined by Equation (13). The t-distribution arises as the marginal distribution for $f(·)$ after integrating the hyperparameters $\mathit{\beta }$ and ${\sigma }^2$. In practice, additional hyperparameters, known as smoothness parameters $\mathbf{B}$, are involved in modelling the correlation function, $c(·,·)$. It is often impractical to give $\mathbf{B}$ a fully analytical Bayesian treatment as integrating the posterior distribution analytically with respect to these parameters is generally impossible. One straightforward approach is to keep $\mathbf{B}$ fixed. Alternatively, numerical methods, such as Markov chain Monte Carlo (MCMC) sampling, can be used to integrate the posterior distribution, though this is computationally intensive. A practical and robust approach is to estimate the hyperparameters of $c(·,·)$ from the posterior distribution and substitute these estimates into $c(·,·)$ in the relevant formulae [36]. These estimates can be obtained using the posterior mode combined with a cross-validation approach [37]. The GEM-SA tool can estimate the smoothness parameters using either method.

4.2.2. Analysis of Main Effects and Interactions

This section explains how sensitivity metrics discussed earlier can be estimated using the GP posterior distribution obtained in Section 4.2.1. An important insight from [32] is that inferences about $f(.)$ can be used to derive information about the main and interaction effects of $f(.)$. This is because these effects are linear functions of $f(.)$ and $t_{d+n}$ after standardisation, as illustrated in (13). As a result, the derived posterior for the main and interaction effects will also be $t_{d+n}$. In particular, if the posterior mean of $f(.)$ is expressed as shown in Equation (14), subsequently for

$$E(Y\mid {\mathbf{x}}_p)=\int {\chi }_{-p}f\left(\mathbf{x}\right)d{G}_{-p\mid p}({\mathbf{x}}_{-p}\mid {\mathbf{x}}_p)$$

(18)

(where ${\chi }_{-p}$ refers to the input space corresponding to ${\mathbf{x}}_{-p}$, while $G_{-p\mid p}({\mathbf{x}}_{-p}\mid {\mathbf{x}}_p)$ represents the conditional distribution of ${\mathbf{x}}_{-p}$ given ${\mathbf{x}}_p$ under G), the posterior mean of this quantity can be written as:

$$E_{\mathrm{post}}E(Y\mid \mathbf{x}p)=Rp\left(\mathbf{x}p\right)\widehat{\mathit{\beta }}+Tp\left({\mathbf{x}}_p\right)\mathbf{e}$$

(19)

where

$$R_p\left({\mathbf{x}}_p\right)={\int }_{{\chi }_{-p}}h{\left(\mathbf{x}\right)}^{T}d{G}_{-p\mid p}({\mathbf{x}}_{-p}\mid {\mathbf{x}}_p),$$

(20)

$$T_p\left({\mathbf{x}}_p\right)={\int }_{{\chi }_{-p}}t{\left(\mathbf{x}\right)}^{T}d{G}_{-p\mid p}({\mathbf{x}}_{-p}\mid {\mathbf{x}}_p)$$

(21)

and $\mathbf{e}=A^{-1}(\mathbf{y}-H\widehat{\mathit{\beta }})$.

Similarly, the mean of posterior for ME or interaction can be derived in the following way:

$$E_{post}\left\{z_i\left(x_i\right)\right\}=\{R_i\left(x_i\right)-R\}\widehat{\mathit{\beta }}+\{T_i\left(x_i\right)-T\}\mathbf{e}.$$

(22)

In a similar manner, the standard deviations of the MEs and interactions can be derived; see [32] for more details on the computational aspects.

The posterior mean of the ME $E_{pos}\left(z_i\left(x_i\right)\right)$ can be plotted against $x_i$, with bounds representing, for example, plus and minus two posterior standard deviations. By standardising the input variables, we can visualise $E_{pos}\left(z_i\left(x_i\right)\right)$ for $i=1,\dots ,d$ on a single plot. This provides a concise graphical summary of the influence that each input variable exerts on the model’s output. Section 5 will present this plot using the dataset introduced in Section 3, offering insights into the impact of different parameters on wave overtopping predictions.

Direct posterior inference for the variance-based measures, specifically $V_{T_i}$ discussed in Section 4.1.2, presents greater challenges due to their nature as quadratic functionals of the underlying function $f(.)$. These measures are more complex than simple linear functionals like the main effects. To handle this computational complexity and derive these sensitivity measures within a Bayesian framework, advanced techniques are required. An in-depth exploration of these methods, including how Gaussian process (GP) emulators can be utilised to compute such measures, is detailed in [32].

5. Results

5.1. Data Preparation and Initial Examination of the CLASH Dataset

This section describes the pre-processing steps taken prior to implementing a novel probabilistic sensitivity analysis approach and predictive modelling for the wave overtopping dataset (Section 3). Initial examination and visualisation of the CLASH dataset in Section 3 revealed a highly non-linear and complex relationship between variables. As a result, standard linear regression models are not suitable for modelling wave overtopping data with such a complex dependency structure. Previous studies have attempted to model these intricate relationships between variables in wave run-up and overtopping using artificial neural networks (ANNs) [38,39].

The emulator-based sensitivity analysis, outlined in Section 4.2, will be used to perform an SA of the parameters influencing wave overtopping. This probabilistic method is developed based on the GP regression as a computationally efficient non-parametric Bayesian machine learning technique. To train the GP model required for computing the SA measures described in Section 4.1.2 and Section 4.2.2, and to address the other computational objectives discussed in Section 4, the following steps were implemented:

• Database cleaning and the selection of a highly reliable subset.
• Perform an exploratory analysis of the subset.
• Fit a Gaussian process regression model for the selected subset of the dataset.
• Compute the SA measures, including the variance-based indexes and the main effects.
• Illustrate the corresponding SA plots.
• Interpret the SA results, perform an uncertainty analysis, and draw conclusions about the most influencing input parameters affecting the wave overtopping.

The CLASH wave overtopping database contains observations and measurements collected from several sources with varying levels of reliability. Numerical modelling has been implemented to fill the gaps in the dataset, and in order to discern between observed and modelled data, each row was given a reliability rating ranging between 1 and 4, with 1 being for “highly reliable” data points and 4 for “highly unreliable” data.

Selecting only the highest reliability entries in the CLASH dataset offers two significant benefits: it reduces the training time for the Gaussian process regression (GPR) model and increases accuracy due to the improved quality of the training points. This initial screening of the dataset, focusing on high-reliability data, results in a dataset of size $3385\times 29$. After removing any entries with missing values, the final dataset is further reduced to 3208 observations, consisting of 27 feature variables and 1 output ($q[m^3/s.m]$, which represents the wave overtopping discharge per unit width. The choice of using GPR for the model development is well justified because GPR is known for its efficiency in constructing reliable models with far fewer training points compared with methods like artificial neural networks (ANNs). ANNs, as employed in previous studies, typically require large training datasets to generalise effectively. In contrast, GPR’s ability to incorporate uncertainty and work with smaller, higher-quality datasets makes it a more appropriate choice for this study. This further supports the adoption of the GPR approach for modelling wave overtopping, especially given the non-linear and complex dependency structure present in the data.

The dataset screening process was further refined by eliminating non-informative feature variables. This was accomplished through an initial data analysis involving the regression of feature variables against the output variable. During this analysis, the scatter plot (Figure 5 left) and histogram (Figure 5 right) of ‘Width of Berm (B)’ revealed that the recorded values for this variable were predominantly zeros. This presented a potential issue for the modelling process. Not only does a variable with mostly zero values contribute little information to the model, but it can also cause numerical challenges during the training phase, particularly for data-driven models like the GPR approach employed here. To mitigate these issues, this study opted to focus on analysing a smaller subset of the dataset, consisting of 330 high-reliability data points. By selecting this high-quality subset, the aim was to probabilistically identify the most suitable model linking wave overtopping to the remaining 27 input parameters. This step ensures that the model is developed using informative and reliable data, thus improving the overall accuracy and robustness of the predictive modelling process.

5.2. GP-Based Sensitivity Analysis for the Wave Overtopping Dataset

The SA of the overtopping parameter, q, with respect to the changes in the input parameters described in Section 5.1, has been conducted. The results, expressed in terms of variance-based measures and total effects, are presented in

Table 2. The emulator-based SA of the overtopping parameter with respect to the changes in other input parameters. The seven parameters contributing the greatest variance are in bold.

Parameters	Variance (%)	Total Effect
Signif wave height (${\mathbf{H}}_{\mathbf{m},\mathbf{deep}}$)	10.93	11.22
Peak period in the deep (${\mathbf{T}}_{\mathbf{p},\mathbf{deep}}$)	5.53	5.83
Mean period m2/m0, deep ($T_{m,deep}$)	0.79	0.89
Mean period, deep (${\mathbf{T}}_{\mathbf{m}-\mathbf{1},\mathbf{deep}}$)	7.63	7.77
Off-shore Water depth, ($H_{deep}$)	1.93	2.04
Slope of foreshore (m)	1.28	1.38
Angle of wave attack ($\beta$)	0.85	0.96
Water depth at toe (h)	1.59	1.69
Signif wave height at toe (${\mathbf{H}}_{\mathbf{m},\mathbf{toe}}$)	8.11	8.22
Peak period, toe (${\mathbf{T}}_{\mathbf{p},\mathbf{toe}}$)	15.75	15.86
Mean wave period, toe (${\mathbf{T}}_{\mathbf{m},\mathbf{toe}}$)	8.48	8.59
Spectral wave period at toe ($T_{m-1,toe}$)	4.06	4.17
Water depth on toe ($h_t$)	1.27	1.40
Toe width ($B_t$)	0.95	1.06
Roughness/perm factor (${\gamma }_f$)	0.06	0.17
Cot downward slope, berm ($cot{\alpha }_d$)	4.01	4.30
Cot upward slope, berm ($cot{\alpha }_u$)	0.65	0.68
Cot slope, excl berm ($cot{\alpha }_{\mathbf{excl}}$)	16.91	17.02
Cot slope, incl berm ($cot{\alpha }_{incl}$)	1.07	1.37
Crest freeboard ($R_c$)	0.74	1.04
Width of berm (B)	2.15	2.22
Water depth on berm ($h_b$)	2.51	2.65
tan of slope of berm ($tanab$)	0.32	0.43
Width of horizontally schematised berm ($B_h$)	1.50	1.62
Width of crest ($G_c$)	0.39	0.50
Armour crest freeboard ($A_c$)	0.20	0.31
Total variance (%)	99.64
Estimated mean output	0.00779018
Estimated variance output	$1.35091\times {10}^{-6}$

.

To identify the variables that are most influential on the overtopping parameter, the correlation coefficients between the input variables and the wave overtopping discharge parameter, q, were examined. Figure 6 illustrates the correlation matrix of these variables. Several strong positive correlations are apparent. For example, significant correlations are observed between $(H_{(m,deep)},q)$, $(T_{(p,deep)},q)$, and $(T_{(m,deep)},q)$. In contrast, the correlations between $(m,q)$ and $(b,q)$ are very weak. Despite Figure 6 being useful for assessing the strength and direction of correlations between the input variables and q, it does not provide insight into the underlying functional relationships between the input parameters and q.

Figure 7 illustrates the scatter plots between several input parameters (mainly structural and geomorphological parameters, as listed in Table 1) and the wave overtopping discharge parameter, q. Due to space constraints, not all input parameters could be included in this figure. From the scatter plots and the computed correlations between each input variable and q, it becomes clear that the relationships between these variables and the wave overtopping discharge are highly non-linear. Moreover, some input parameters, such as $cot{\alpha }_{excl}$ or $cot{\alpha }_{incl}$, exhibit distributions with multiple modes, suggesting that their relationships with q and even other input parameters defy interpretation via any known mathematical functional forms. Consequently, it would be misleading to rely solely on correlation coefficients for the SA of this complex, non-linear system, as highlighted by [40].

Using correlation coefficients in such a scenario could obscure the true effects of input variables on the prediction of the wave overtopping discharge parameter. Additionally, conventional SA methods, including Markov chain Monte Carlo (MCMC)-based approaches [41], may struggle to capture the complex, non-linear dependencies and interactions inherent in this model.

To examine the variation of model output with respect to the uncertainty in the input variables, we further develop the probabilistic SA approaches described in Section 4.1. The GP emulator helps to overcome the computational challenges involved in calculating the SA indices for the complex, non-linear system considered in this study. These SA indices allow us to evaluate the contribution of each input variable to the overall variability of the system. Specifically, by leveraging the GP emulator, we efficiently compute these indices and quantify the influence of each variable on the wave overtopping model’s output. The probabilistic SA enables us to assess how uncertainty in the input variables propagates through the model and affects the overtopping discharge parameter. This approach is crucial for understanding which variables are the most influential in predicting wave overtopping, allowing for a more focused analysis and potential simplifications in the model complexity without compromising accuracy. The SA indices thus provide a powerful tool for identifying key drivers of uncertainty and guiding decision making in model refinement and resource allocation.

Table 2 provides the results of the emulator-based SA for the wave overtopping parameter with respect to changes in the input parameters discussed in Table 1. The results highlight that a significant portion of the variability in wave overtopping discharge, approximately 20%, can be attributed to ‘significant wave height’ measurements. Furthermore, the results indicate that $H_{m,deep}$ accounts for 10.93% and $H_{m,toe}$ accounts for 8.11% of the total variance. Further investigation into the relationship between these hydrodynamic parameters and wave overtopping discharge reveals strong correlations. Figure 8 illustrates that significant wave height parameters ($H_{m,deep},H_{m,toe}$) are also highly correlated. This correlation indicates that significant wave height is a key factor influencing wave overtopping and suggests it could be pivotal in simplifying the final model or reducing model dimensionality. Given these findings, a significant wave height should be included as one of the primary features in any predictive model for wave overtopping. This inclusion will enhance model accuracy and focus, leveraging the substantial variance contribution from these hydrodynamic parameters.

The most influential input parameter identified by the models is the ‘mean cotangent of structure slope without contribution to the berm’, denoted by $cot{\alpha }_{excl}$, which contributes approximately 17% to the total variance of the wave overtopping output. This finding of the model aligns with fundamental underlying theories, which emphasise the importance of the slope in determining the wave breaker index (Iribarren number). The Iribarren number provides an important role in the wave breaking process, affecting turbulent kinetic energy and momentum transfer, both of which are crucial factors in wave overtopping at coastal defence structures.

In addition to slope, several wave period-related parameters also significantly influence wave overtopping. The ‘Peak wave period at toe’, $T_{p,toe}$, contributes 15.75%, while the ‘off-shore peak wave period in the deep water’, $T_{p,deep}$, and ‘off-shore spectral wave period’, $T_{m-1,deep}$, contribute 5.53% and 7.63% of total variance contribution toward the output’s variance, respectively. These findings are consistent with existing theories as wave periods directly impact the frequency and intensity of wave impacts on coastal structures.

Furthermore, parameters such as ‘Roughness/permeability factor for the structure (${\gamma }_f$)’, ’Width of the structure crest ($G_s$)’, and ‘Armour crest freeboard ($G_s)$)’ contribute relatively little to the overall variance in the wave overtopping model (see Table 2). These parameters, with their relatively insignificant contributions to variance, suggest that they are not critical for predicting the wave overtopping discharge. Notably, 99% of the variability in the output is attributed to MEs, which further underscores the suitability of using GPR in this analysis.

Figure 9 illustrates the estimated MEs, $E\left(q\right|X_i)$, where $X_i$, defined in the general case in Equation (18), is an input parameter (detailed in Table 1). These MEs are approximated using the Gaussian process (GP) model as described in Section 4.2.2. This approach allows us to evaluate the sensitivity of the wave overtopping parameter with respect to individual input variables while maintaining computational efficiency and flexibility. To compute the SA measures reported in Table 2 and the MEs for each input parameter, the GP model was trained using only 330 high-reliability data points after cleaning the dataset, as described in Section 5.1.

The SA measures proposed in this study are computed by developing a custom code in R. To facilitate computations, the input variables and the output (wave overtopping discharge) were standardised prior to fitting the GP emulator. Standardising the data not only helps to avoid numerical issues but also makes the computations more consistent across different models. Once the GP model was trained, the results were then transformed back to their original scale to make them interpretable and applicable to real-world scenarios. This process, as supported by [42], showcases the advantages of using standardised data, particularly when working with complex models like the GP emulator. This ensures numerical stability while maintaining flexibility across a wide range of computational implementations. The final model results can then be transformed back to the original scale.

Figure 9 shows the estimated MEs of the overtopping parameter in response to changes in various input parameters. The figure highlights that the overtopping parameter is most sensitive to uncertainties in terms of $cot{\alpha }_{excl},T_{p,toe},H_{m,deep}$, $T_{m,toe},H_{m,toe},T_{m-1,deep}$, and $T_{p,deep}$, respectively. The width of the uncertainty bands depicted in Figure 9 represents the uncertainty of the GP emulator linked to each input. This uncertainty quantification provides insight into the reliability of predictions for each parameter. For instance, when $H_{m,deep}$ is fixed at a value of 0.5, the corresponding point on the graph represents the expected value of the overtopping discharge parameter, q, derived by averaging across the remaining input parameters.

The emulator-based SA method allows for systematic exploration by varying small groups of input variables, while others are held fixed at their default values. This approach enables the examination of the sensitivity of the output (overtopping parameter, q) to particular input variables while controlling for other factors. When comparing the thickness of the ME plots, it becomes clear that there is less uncertainty associated with certain key variables, such as $cot{\alpha }_{excl},T_{p,toe},H_{m,deep},T_{m,toe},H_{m,toe},T_{m-1,deep}$, and $T_{p,deep}$, which are illustrated in blue, than the rest of the input parameters. The thickness of the ME plots is the result of simulating multiple (in this case, 200) realisations from the posterior distribution of the output (overtopping parameter q). These realisations are computed at a regularly spaced grid of input points, as described in Section 4.2.1. The thicker ME bands indicate more uncertainty, while thinner bands suggest a higher degree of confidence in the model’s predictions for those variables. The variability in the spread of ME lines illustrates the combined uncertainty present in both the model and input parameters. Thus, these sensitivity measures themselves are subject to uncertainty, and generating multiple realisations helps to reflect that uncertainty in a probabilistic framework.

To enhance the clarity of the visualisations of the ME plots illustrated in Figure 9, we now focus on the ME plots of the most sensitive input parameters, namely $cot{\alpha }_{excl},T_{p,toe},$$H_{m,deep},T_{m,toe}$, $H_{m,toe},T_{m-1,deep}$, and $T_{p,deep}$. The improved visualisations are presented in Figure 10, which provides a more detailed focus on the most influential parameters. In each ME plot, the blue solid line represents the posterior mean of the ME, calculated as described in Section 4.2.2. This line indicates the central tendency of the model’s prediction for the effect of each input parameter on the output given the available data. The red dotted lines show the 95% confidence interval around the ME (see Section 4.2.2 for details).

This confidence interval is derived from the posterior distribution of the emulator and illustrates the uncertainty in the estimated ME. It is worth noting that this confidence interval offers an alternative method for assessing uncertainty compared with the approach of using multiple realisations simulated from the posterior distribution of the overtopping parameter, q. Both methods provide insight into the uncertainty surrounding the estimated ME, though they approach the analysis from slightly different perspectives.

One final point to note about the proposed method in this study is that the machine learning Gaussian process (GP) sensitivity analysis results for the wave overtopping discharge parameter are consistent with the findings derived from the variance-based sensitivity analysis (SA) measures, as reported in Table 2. This agreement reinforces the robustness of the GP-based approach, indicating that it effectively captures the influential input parameters and their contributions to the variability of the wave overtopping discharge parameter. The concordance between the two approaches—one probabilistic (GP-based) and the other variance-based—validates the applicability of the GP method for complex, non-linear systems like the one studied here. It suggests that the GP emulator can be reliably used as a computationally efficient alternative to traditional variance-based methods for sensitivity analysis, especially in cases where the underlying functional relationships between inputs and outputs are highly non-linear or difficult to model with simpler techniques.

6. Conclusions

This study outlines a novel probabilistic machine learning approach for sensitivity analysis of complex, non-linear, multi-variable datasets on wave overtopping from coastal defences. The overtopping data used in this study are based on 163 two- and three-dimensional physical modelling tests as well as full-scale field-based (prototype) measurements reported in the CLASH dataset. The homogeneous database of this study includes measurements for all common types of hard-engineered coastal protection structures and captures data on the hydrodynamics, geomorphological and structural parameters influencing wave overtopping discharge.

The analysis revealed that the overtopping process is highly complex, driven by non-linear interactions between wave kinematics, geomorphological features, and structural parameters. For the first time, this paper proposes a mathematically robust and computationally efficient framework based on the GP emulator to perform a sensitivity analysis of the multi-faceted wave overtopping problem. This framework analyses how variations in key input variables influence the wave overtopping from coastal defences.

The GP model developed in this study was successfully tested and validated using data from the CLASH database. The proposed method allows for an effective sensitivity analysis with significantly fewer model runs compared with conventional SA approaches, including MCMC-based methods [33].

The machine learning-based sensitivity analysis model presented in this paper is crucial for predicting wave overtopping discharge from critical coastal infrastructures, a key factor in forecasting coastal flooding. The results of the sensitivity analysis can be directly applied to improve predictive models by highlighting the parameters with the most significant impact on wave overtopping, thereby enhancing the accuracy of assessments.

Using the GP emulator, computationally expensive sensitivity measures, such as variance-based analyses and main effects (MEs) of the overtopping parameters, were efficiently computed. The study found that significant wave height features ($H_{m,deep},H_{m,toe}$), slope without berm contribution ($cot{\alpha }_{excl}$), peak wave period at the toe ($T_{p,toe}$), offshore peak wave period ($T_{p,deep}$), and offshore spectral wave period ($T_{m-1,deep}$) are the most influential factors determining the intensity of wave overtopping discharge.

The results from the GP-based SA methods were leveraged to simplify and develop a robust predictive model for wave overtopping and coastal flooding based on the CLASH database. Figure 11 illustrates the initial results of predictive modelling for wave overtopping discharges from vertical defence structures (e.g., seawalls). The GP-based model demonstrates an efficient and reliable prediction of wave overtopping volumes, showing potential to outperform existing methods, such as those based on artificial neural networks, as reported in [38,39]. Additionally, due to the capability of GPs to make accurate predictions with far fewer data points, compared with ANNs, the proposed method allows for SA and modelling on smaller data subsets, enabling rigorous analyses, even with limited data availability.

The mathematical framework presented in this paper, along with the machine learning Gaussian process (GP)-based model developed, can significantly contribute to creating a reliable decision support tool. This tool would empower coastal scientists and engineers to better evaluate and predict the performance of key coastal infrastructures under the growing threats posed by climate change.

This study provides a foundation for future work. While our method exhibits impressive predictive power on the CLASH dataset, many real-world coastlines are under-represented in these data. Future research should investigate the generalisability of this approach, which would require an expanded data collection. This paper primarily focused on topological and hydrological parameters, yet other sources of variability—such as wind direction and speed—are known to affect overtopping. The proposed methodology could be extended to explore these additional parameters, although collecting such data may be expensive. This enhances the value of our data-efficient machine learning approach.

With the availability of more data, future work should explore methods to enhance the computational efficiency of the GP model, particularly through approaches such as sparse GPs [43].

Although we introduced a novel method and demonstrated its promising performance relative to the current state of the art, future research should undertake a thorough comparison with other existing methods. Such work could help clarify the specific use cases where our approach excels, contributing to a deeper understanding of when and where to apply various predictive techniques.

One particularly impactful area for future exploration is spatio-temporal forecasting. Given sufficient data, GPs could be employed to predict the effects of climate change on coastal regions, forecasting the likelihood of wave overtopping in the coming decades. These predictive insights could help identify critical risks to vulnerable coastal areas well in advance, enabling timely and effective preventive measures.