Global Sensitivity Analysis of Groundwater Related Dike Stability under Extreme Loading Conditions

Abstract

With up to 15% of the world’s population being protected by dikes from flooding, climate-change-induced river levels may dramatically increase the flood risk of these societies. Reliable assessments of dike stability will become increasingly important, but groundwater flow through dikes is often oversimplified due to limited understanding of the important process parameters. To improve the understanding of these parameters, we performed a global sensitivity analysis on a comprehensive hydro-stability model. The sensitivity analysis encompassed fifteen parameters related to geometry, drainage conditions and material properties. The following three sensitivity settings were selected to characterize model behavior: parameter prioritization, trend identification and interaction qualification. The first two showed that dike stability is mostly dependent on the dike slope, followed by the type of subsurface material. Interaction quantification indicated a very prominent interaction between the dike and subsurface material, as it influences both groundwater conditions and dike stability directly. Despite our relatively simple model setup, a database containing the results of the extensive Monte Carlo analysis succeeded in finding most of the unsafe sections identified by the official inspection results. This supports the applicability of our results and demonstrates that both geometry and subsurface parameters affect the groundwater conditions and dike stability.

1. Introduction

Over 45 major flood events occurred in Europe between 1950 and 2005 that each resulted in more than 70 fatalities or a collected economic damage of EUR 7.6 × 10¹⁰ [1]. As a result, many flood prone areas have an extensive network of artificially elevated levees or dikes, which, along Europe’s major rivers, add up to a length of approximately 60,000 km [2]. To ensure the safety of people living behind dikes, continuous maintenance and reinforcements are needed to warrant the stability of dikes and their proper functioning during high water events. Climate change, e.g., earlier snow melt or an increase in extreme precipitation events in the upstream drainage area [3], poses a new threat that may increase the risk of a society to flooding [4]. To maintain safety levels under changing climatic conditions, major investments are needed for dike maintenance and reinforcement, of which the costs for the latter are in the order of EUR 1–20 million per kilometer [1]. Improved knowledge of the processes during and following a high-water event that can result in dike failure is crucial for more cost-effective dike reinforcements, which may reduce the total expenditures on dike reinforcements substantially and can support more societally acceptable flood defense measures [5].

Many dike failure mechanisms are related to local groundwater conditions and pore pressures in the dike body. In response to elevated river stages, changing groundwater conditions may increase the pore pressure and, thus, reduce the effective normal strength, while, at the same time, the lateral load of river water pushing against the dike is increased. Therefore, parts of its inner or outer slope may slip, or the dike may slide along its base (soil slip sensu lato), threatening the structural integrity of the dike. Accordingly, when analyzing a dike failure hazard in relation to high groundwater levels and river stages, multiple failure mechanisms must be considered.

Although we acknowledge that critical groundwater heads in dikes are primarily driven by the occurrence and nature of high-water events, their variation in space and time also depends on surface geometry and subsurface properties. Previous research on this topic can roughly be divided into the following three categories: research focusing purely on hydrology, research focusing on single cases or research assessing variability in either surface geometry or subsurface properties. An extensive analysis on only the hydrology near a river dike was provided [6,7]. Research on a single case [8] was often also focused on the effect of artificial reinforcements [9,10]. Attempts including a local sensitivity analysis investigated the influence of material properties [11] or geometry [12] on the stability of embankments. In sum, none of these previous studies conducted a full analysis that considers both variations in hydrological parameters as factors influencing the stability of a dike.

Such a full analysis has already been widely applied in landslide probability modelling [13,14,15]. Nonetheless, whereas slope hydrology is mostly dependent on rainfall infiltration, flow as a result of elevated river water levels mostly occurs horizontally and is often affected by intersecting aquitards, resulting in very different patterns of groundwater flow and pore pressure buildup. These differences inhibit the direct application of results from landslide modelling to river dike failure scenarios. Nonetheless, to assess variations in hydrological parameters as factors influencing the stability of a dike, a comprehensive hydro-stability system needs to be modelled.

To quantify the model and parameter uncertainty, local sensitivity approaches estimating the partial derivatives of the model at a specific point in the parameter space are no longer sufficient [16]. Alternatively, a global sensitivity analysis can handle nonlinearity and local variations expected in more complex models. A global sensitivity analysis considers the entire variation of the input factors [17]. Global sensitivity analysis recently gained interest in environmental modelling [18], and the different goals and methods related to global sensitivity in environmental models have been extensively reviewed [19,20]. Hydrological models have seen a similar rise in interest for global sensitivity analysis, both from a methodological point of view [21,22,23] and for analyzing geo-hydrological systems [24] and slope stability uncertainty [25,26,27]. Though some attempts have been made to analyze the global sensitivity of dike stability based on the uncertainty in its internal characteristics [11,25], no complete sensitivity analysis covering both geometrical as subsurface characteristics has yet been made.

Thus, to assess both geometrical and subsurface characteristics of dikes, we created a coupled high-resolution groundwater model and a limit equilibrium stability analysis. To constrain our results and to highlight first-order relationships, we evaluated the stability under the most critical loading conditions and the maximum pore water pressures. Three failure mechanisms that affect the macro-stability of a dike were considered, being inner slope stability, outer slope stability and basal sliding, as their occurrence is directly linked to the geometry of a dike and its composition.

To add to previous research on both river dike hydrology and global sensitivity analysis, the goal of our analysis was to identify the overall stability of a dike in terms of its factor of safety F under different hydrological loading conditions, subsurface geometries and material properties, and including pore pressure calculations. In this work, our research goal translates into the following three sensitivity settings to characterize model behavior: parameter prioritization, trend identification and interaction qualification. We aim to provide insights in each of these settings, while maintaining a reasonable computation time. For parameter prioritization we used the Elementary Effect [28] and the delta-importance measure [29]. As we are mostly interested in those factors that could lead to unstable dikes, we used a regional sensitivity analysis (RSA) to perform trend identification and identify regions in the parameter space with a safety factor below one. For interaction qualification, we focused on the subsurface, and used response surfaces to analyze the interaction activate between hydrology, material characteristics and dike stability. Moreover, the outcome of this global sensitivity analysis can be used to inform semi-qualitative assessments of dike stability as often applied in regional inventories, and the conducted set of model runs is used for a direct comparison to a case study site in the Netherlands.

2. Materials and Methods

2.1. Case-Study Schematization

We applied the global sensitivity analysis on a cross-section from the river to the hinterland behind the dike. Fifteen parameters (

Table 1. Name, symbol and range of the model parameters. A visualization of each of the parameters is shown in Figure 1. Layer type descriptions are found in Table 2.

Parameter	Symbol	Range	Unit
Dike Height	$D_h$	3–10	m
Dike crest width	$D_w$	2–5	m
Dike slope	$D_s$	0.2–1	m/m
Dike type	$D_{typ}$	C-CL-L-SL-S	-
Upper layer thickness	$U_{thck}$	0.3–1.9	m
Upper layer type	$U_{typ}$	C-CL-L-SL-S	-
Lower layer thickness	$L_{thck}$	5–10	m
Lower layer type	$L_{typ}$	C-CL-L-SL-S	-
Foreland width	$F_w$	0–100	m
Drainage depth	$D{r}_d$	0.1–2	m
Drainage spacing	$D{r}_s$	1–20	m
Riverbed slope	$R_d$	0.33	m/m
River depth	$R_d$	$0.9*(U_{thck}$$+ L_{thck}$)	m
Flood height	$H$	$Dh$	m
Drawdown time	$T_d$	1	days

) describe the cross-section, subdivided in the following three groups: topographical parameters, subsurface parameters and human management parameters. The topography is described by the following six parameters: the dike height ($D_h$), dike crest width ($D_w$), dike slope ($D_s$) and floodplain width ($F_w$), riverbed slope ($R_s$) and river depth ($R_d$) (Figure 1). The subsurface is described by five parameters, which divide the subsurface in three units with uniform characteristics. The dike is schematized by its material type ($D_{typ}$), in addition to the previously mentioned geometry parameters $D_{typ}$, $D_w$ and $D_s$. The upper subsurface layer is schematized by a thickness ($U_{thck}$) and material type ($U_{typ}$). The same applies to the lower subsurface layer ($L_{thck}$, $L_{typ}$). Two parameters describe human management by specifying the drainage conditions behind the dike (Figure 1), which are drain spacing ($D{r}_s$) and drainage depth ($D{r}_d$).

2.2. Coupled Hydrology-Stability Model

2.2.1. Hydrological Model Setup

We added a stability module to a groundwater model to examine the stability of the schematized cross-section. The MODFLOW 6 software [30,31], a Modular Three-Dimensional Finite-Difference Groundwater Flow Model, simulates the groundwater conditions, which are included in the structural stability using the Generalized Limit Equilibrium Method (GLEM) [32,33]. The MODFLOW 6 hydrological model is constrained by the river water level and drainage depth. On the river side, the imposed river stage at the top of the dike constitutes a head-controlled boundary condition that is regulated by the MODFLOW river package. A head-controlled boundary on the inner side of the dike is regulated by the MODFLOW drain package, which creates outflow-only seepage points [30]. Seepage is possible if the hydraulic head at the surface is higher than the surface elevation. In addition, at a distance of $D{r}_s$ behind the dike, a ditch is located with a depth of $D{r}_s$ (Table 1), enabling faster drainage of deeper layers. The cell size is 0.5 m in all directions, which enables the assessment of small-scale spatial variation while retaining the computational efficiency needed. The model first performed a steady-state simulation, in which the river stage ($H$) was set at the dike crest elevation. It was assumed that under these conditions, the pore pressures reach their most critical values. After the steady-state simulation, a rapid drawdown of $H$ from the dike crest to the dike toe in a time period ($T_d$) of one day was transiently simulated with a 3-hour timestep. As pore pressures in the dike do not immediately follow river water level changes and the stabilizing effect of the high river water levels is absent, these conditions might provoke outer slope failure. An exploratory sensitivity analysis showed a time step of 3 h does not significantly impact results when compared to smaller time steps.

2.2.2. Dike Slope Stability Model Setup

Dike stability is expressed by the factor of safety ($F$), which is calculated separately for inner slope failure, basal sliding and outer slope failure using the Generalized Limit Equilibrium Method [32,33], resulting in the following three safety factors: $F_{inner}$, $F_{lat}$, $F_{outer}$. This method solves both moment and force equilibrium on a slip surface for different ratios between the vertical and horizontal inter-slice shear forces. The relationship between the magnitude of the inter-slice shear and normal forces is assumed to be constant [34]. The factors of safety presented in this paper always represent the factor of safety of the most critical circular slip surface, derived by an effective critical slip surface minimization technique adapted from [35]. To constrain the slip surface for inner and outer sliding the slip surface is forced to enter on the dike crest or on the corresponding dike side. For basal sliding an infinite slump radius is assumed, which results in a horizontal slip surface and enables the calculation of $F$ using only force equilibrium. To ignore very small slumps not causing dike breaches, a minimum cross-sectional slip surface area of 5 m² is imposed.

2.3. Workflow and Parameters for Global Sensitivity Analysis

As conducting the entire global sensitivity analysis with all parameters was not a feasible option due to increasing computation times, we first screened our inputs using the Elementary Effect test (EE) [28,36]. Hereby, we identified which input parameters have a small contribution to the output variation in high-dimensional models and can, therefore, be set to a fixed value. Subsequently, the global sensitivity analysis was performed using a Monte-Carlo (MC) approach, as it captures the entire range of possible combinations, while facilitating parameter interaction understanding (Figure 2). The Delta Moment-Independent measure (DMI) [29] was used to quantify the sensitivity of the factors contributing significantly to the variation in the model output.

All variable parameters defining the cross-section are screened using the Elementary Effect test (EE). As no information is available a priori, the parameters are sampled from a uniform distribution within their possible range (Table 1). The material types ($D_{typ}$, $U_{typ}$, $L_{typ}$) each represent a single lithological class (Figure 3), having deterministic attributes used in the coupled hydrology-stability model. We selected the values for these attributes (

Table 2. Subsurface types used in the model, related to the $D_{typ}$, $U_{typ}$ and $L_{typ}$ parameters. The subsurface type is linked to the hydraulic conductivity ($K_{sat}$ ), drained cohesion ($c\prime$ ), effective friction angle ($\varphi$ ), the bulk unit weight ($\gamma$ ) and the saturated bulk unit weight (${\gamma }_{sat}$ ).

Subsurface Type	Abbreviation	${\mathit{K}}_{\mathit{s}\mathit{a}\mathit{t}}(m{d}^{-1})$	$\mathit{c}\prime \left(kPa\right)$	$\mathit{\varphi }{(}^o)$	$\mathit{\gamma }\left(kN{m}^{-3}\right)$	${\mathit{\gamma }}_{\mathit{s}\mathit{a}\mathit{t}}\left(kN{m}^{-3}\right)$
Clay	C	0.13	5.0	17.5	17	17
Clay-Loam	CL	0.18	4.0	22.5	18	18
Loam	L	0.19	1.0	30.0	20	20
Sandy Loam	SL	1.54	0.5	31.25	19.5	19.5
Sand	S	8.94	0.0	32.5	18	20

) based on characteristic values in the literature. The hydraulic conductivity ($K_{sat}$) is derived from the geometric mean of multiple laboratory $K_{sat}$ measurements of soils with a relatively high density [37]. The cohesion ($c\prime$), the effective friction angle ($\varphi$), the bulk unit weight ($\gamma$) and the saturated bulk unit weight (${\gamma }_{sat}$) are in line with the European standardized characteristics for soil stability [38].

2.4. Parameter Prioritization

This sensitivity setting focusses on identifying the input factor (parameter) that has the largest effect on the model output, e.g., the input factor that, when fixed, decreases output variability the most. For this, we used the Elementary Effect test (EE) and the Delta Moment-Independent measure (DMI).

The Elementary Effect test [36] is an effective sensitivity analysis (SA) method that is widely used for screening practices, as it provides relatively good results at small sample sizes [28]. This method is basically a One-At-a-Time analysis, which is extended to the full input factor space. The original method [28] measures sensitivity in terms of $\mu$, indicating the first order influence or elementary effect, and $\sigma$, indicating second order influences, being nonlinearity or interaction effects. For each input factor, random baseline points are selected from which the other are varied (Figure 4A). Given function “$y$” (in our case the calculation of $F$), step size $\Delta$ and a random baseline sample $X,$ the elementary effect of input factor $X_i$ is given by

$$E{E}_i\left(X\right)=\frac{y\left(X_1,\dots ,X_{i-1},X_i+\Delta ,X_{i+1},\dots ,X_k\right)-y\left(X_1,\dots ,X_{i-1},X_i,X_{i+1},\dots ,X_k\right)}{\Delta }$$

(1)

The final $\mu$ for any input factor is the mean of the $E{E}_i$ at all baseline points $X$. Non-monotonic models result in both positive and negative $EE$’s for a given input factor, which average out when taking the mean. Therefore, [36] introduced ${\mu }^{*}$, which is the mean of all absolute elementary effects and is found to be suitable for input factor ranking.

$${\mu }_i^{*}=\frac{1}{r}{\displaystyle \sum }_{j=1}^r|E{E}_{i,j}\left(X\right)|$$

(2)

where $\left|E{E}_i\left(j\right)\right|$ is the elementary effect for input factor $X_i$ using the $j$-th step with step size $\Delta$, with r being the number of steps in the parameter space. A ${\mu }_i^{*}$ value near zero indicates that the parameter has a small general effect on the output. This measure is used in the Insensitive Factor fixing procedure.

The Delta Moment-Independent measure (DMI) [29] is based on shifts in the probability density function $f_Y\left(y\right)$, contrary to most SA techniques, which are variance-based. Variance-based sensitivity, according to classical utility theory, is not suitable to describe uncertainty in case of a non-normal probability distribution and in case of a non-quadratic utility function. Moreover, the probability density function provides a more complete overview of sensitivity than variance-based techniques. DMI returns the measure (${\delta }_i$), which represents the non-overlapping area between an unconditional input vector $X_i$, including all parameter values, and a conditional input vector $X_i$, consisting of a subset of parameter values (Figure 4B). Mathematically, it is expressed as:

$${\delta }_i=\frac{1}{2}{E}_{X_i}\left[s\left(X_i\right)\right]$$

(3)

with

$$E_{X_i}\left[s\left(X_i\right)\right]=\int f_{X_i}\left(x_i\right)[\int \left|f_Y\left(y\right)-f_{Y|X_i}\left(y\right)\right|dy]d{x}_i$$

(4)

which shows that the input factor specific ${\delta }_i$ depends on the shift in probability density function for multiple conditional inputs $X_i$ and on the underlying probability of that shift to occur. As in our method, a uniform probability function is used, and the mean of the separate shifts represents the final ${\delta }_i$.

2.5. Factor Fixing Procedure

Factor fixing is often used as a SA setting in itself, where the goal is to simplify model and prevent overparameterization [23]. In this research, it was just a means by which we aimed to keep MC simulation runs to an acceptable level, i.e., by fixing the least influential parameters to some nominal value. To provide evidence for the identification of the least influential parameters for the model output, we used an iterative version of the SA repeatability test [40], previously successfully adapted for environmental models [41,42]. This approach focusses on testing the predictive capacity of parameters.

First, 1200 samples of all model parameters are created. The test then consists of the comparison of two conditional input samples, $X2$ and $X3,$ to the previously created unconditional sample, $X1$. Set $X2$ fixes the input factors deemed insensitive at a predetermined value, while $X3$ fixes the input factors deemed sensitive. Afterward, the unconditional result $F\left(X\right)$ is compared with the conditional results $F\left(X2\right)$ and $F\left(X3\right)$. If a correct classification of important and non-important parameters was used, the correlation coefficient ($r^2$) between $F\left(X2\right)$ and $F\left(X1\right)$ approaches 1, while the correlation coefficient of $F\left(X3\right)$ and $F\left(X1\right)$ approaches 0, as the parameters fixed in $X3$ should have a small influence on the results. We iteratively applied this approach, starting with only the most important factor classified as sensitive, and consecutively also classifying the next important parameter as sensitive, until the correlation coefficient $r_{F\left(X2\right)|F\left(X1\right)}^2$ exceeds a certain threshold. A threshold $r_{F\left(X2\right)|F\left(X1\right)}^2$ of 0.95 has been successfully applied [41] to limit the dimensionality of a problem, while retaining sufficient model variability.

To initially rank the parameters from sensitive to insensitive, we used the enhanced Elementary Effects method [36] on the initial sample X1. The iteration was performed for each failure mechanism separately, but we used an inclusive selection strategy, indicating that only those parameters that were found to be insensitive for all failure mechanisms were excluded from the MC-parameters. Although the inclusive approach increases the dimensionality of our problem, it also enables an easy comparison between the different failure mechanisms. As we used a threshold for $r_{F\left(X2\right)|F\left(X1\right)}^2$ of 0.95 as closing criterion for the iterative Factor Fixing procedure, this threshold being lower than 1 indicated that the factors to be fixed still influenced the model outcome, though their influence was limited. As this research was investigating worst case scenarios, any factor to be fixed should have been set at a value that resulted in relatively low safety factors. This value is selected from the input vector $X\left(allparameters\right)$ that results in the globally minimized $F\left(X\right)$ within the specified parameter ranges (

Table 3. Minimum $F\left(X\right)$ per failure mechanism and parameter values of input vector X resulting in that minimum $F\left(X\right)$. See Table 1 for the abbreviations and units of the parameters.

	${\mathit{F}}_{\mathit{m}\mathit{i}\mathit{n}}$	${\mathit{D}}_{\mathit{h}}$	${\mathit{D}}_{\mathit{w}}$	${\mathit{D}}_{\mathit{s}}$	${\mathit{D}}_{\mathit{t}\mathit{y}\mathit{p}}$	${\mathit{U}}_{\mathit{t}\mathit{h}\mathit{c}\mathit{k}}$	${\mathit{U}}_{\mathit{t}\mathit{y}\mathit{p}}$	${\mathit{L}}_{\mathit{t}\mathit{h}\mathit{c}\mathit{k}}$	${\mathit{L}}_{\mathit{t}\mathit{y}\mathit{p}}$	${\mathit{F}}_{\mathit{w}}$	$\mathit{D}{\mathit{r}}_{\mathit{d}}$	$\mathit{D}{\mathit{r}}_{\mathit{s}}$
Basal sliding	0.69	9.97	2.14	1:1	Sand	0.31	Clay	5.07	Clay	50	−1.05	20.0
Inner slope stability	0.00	9.74	3.50	1:1	Sand	1.10	Clay-Loam	7.50	Clay-Loam	50	−1.05	10.5
Outer slope stability	0.26	9.34	4.61	1:1	Sandy Loam	1.82	Loam	5.25	Clay-Loam	85	−0.62	20.0

) for each failure mechanism separately. Initiated at a random starting point, a modification of Powell’s method [43] performs the minimization operation. This method performs a bi-directional search in one dimension, meaning it searches for the local minimum $F\left(X\right)$ by changing only one input parameter. The input parameter is updated to the value resulting in the minimum $F\left(X\right)$, and the bi-directional search is applied to the next input factor. After minimizing all input factors, the intermediate model output is stored, and the first factor is again selected. When the difference between the previous and current intermediate output is lower than a given threshold, the globally minimized input vector $X$ is found. Though we acknowledge that fixing only a single input factor to value in the globally minimized input vector does not necessarily result in the local minimum $F\left(X\right)$ at any given point in the parameter space, we believe that it results in a safety factor near the real minimum at that point.

2.6. Factor Sampling

As neither the real parameter probability density functions of the selected model parameters were known and no information was available on their correlation, we used a uniform uncorrelated sampling strategy. The sensitive factors were sampled using a discrete uniform distribution in their possible range (Table 1), which suggested a known, finite number of outcomes that were equally plausible. We used five steps ($n$) at which to sample discrete values from the minimum value ($a$) to the maximum value ($b$) in the possible range given parameter $p$, i.e.,

$$x_p=a,a+\frac{b-a}{n-1},a+2\frac{b-a}{n-1},\dots ,b$$

(5)

Afterwards, each possible combination of the parameter values was selected for the MC analysis, resulting in a number of $P^n$ model runs, with $P$ being the number of selected parameters.

2.7. Trend Identification and Interaction Qualification

For trend identification and to increase our understanding of model sensitivity, we explored whether any change toward higher parameter values would also lead to a larger safety factor. To this end, we performed a regionalized sensitivity analysis (RSA), which aimed to identify regions in the inputs space that result in an output in a specified zone [23]. In our case, the selected zone was $F\le 1$, as it was intuitively the most interesting region of model outcomes; that of dike failure. To indicate this regional effect, we used $p\left(F\le 1\right)$ as a measure, which is the probability for a fixed parameter to result in $F\le 1$ given the variations of the other parameters. This measure can be easily calculated as a result of the uniform discrete sampling distribution. Interaction qualification uses response surfaces, which directly show the correlation between material properties, groundwater and dike stability. As we focused on a qualitative description and interpretation of these interactions, no quantitative statistical measures were used.

3. Results

3.1. Exemplary Results of the Hydro-Stability Model

The typical model results that highlight the interplay between the model input factors, groundwater hydrology and macro-stability are presented in Figure 5. Focusing on changes in groundwater conditions, the phreatic surface in the dike after the steady-state simulations ($t_0$) seems to be more influenced by dike geometry than material properties. Due to increased drainage, steeper and narrower dikes seem to result in a generally lower phreatic level, as the hydraulic heads at the center of the dike are more directly affected by gradient changes near the surface. Lower river water levels after draw-down also affect the phreatic level, but the amount of lowering depends on the possibility to drain the excess pore pressure. This is reflected in lower groundwater levels at steeper and thus smaller dikes, at dikes on a sandy substrate and in case of longer drawdown times.

For inner slope stability, these results have a more direct effect, as the minimum factors of safety $F$ are always found at $t_0$, when the pore pressures are the highest. For basal sliding, this is the case as well, although this is also influenced by the still high river water levels, which apply a lateral force. As a result, the dike’s stability related to basal sliding increases to infinity at $t_e$, as the river water level is at the dike toe and no driving force is exerted by it. For outer slope stability, the opposite argument applies, and due to a decrease in the lateral river water pressure, the lowest $F$ is most often found at $t_e$ (Figure 5).

These example results clearly show that not only the safety factor, but also the slip surface area and location are dependent on the hydrological conditions. Higher pore pressures are likely to result in larger slumps that result in a greater chance of dike breach, in addition to lower stability. The results also indicate an effect of drawdown time, which is related to the flood wave shape, but in this study the drawdown time ($T_d$) was kept constant.

3.2. Parameter Priorization

3.2.1. Factor Fixing on Globally Minimized Input Vector

The values used to fix the insensitive factors are derived from the input vector $X$ that results in the globally minimized $F\left(X\right)$. The globally minimized factors of safety ($F_{min}$) are 0.69, 0.00 and 0.26 for basal sliding, inner slope stability and outer slope stability, respectively, of which the corresponding input vectors are shown in Table 3. The iterative factor fixing is performed in the order of the mean absolute elementary effect (${\mu }^{*}$) of each input factor per failure mechanism (Figure 6A). If the dike slope, dike height, dike material type, upper subsurface type, lower subsurface type, upper subsurface thickness, drainage spacing and foreland width are classified as sensitive, the $r_{F\left(X2\right)|F\left(X1\right)}^2\ge 0.95$ for all the failure mechanisms. These eight factors (Figure 6) are, thus, selected for the MC-analysis. Using the final $F\left(X2\right)$ and $F\left(X1\right)$, a mean error in $F$ of −2.76·${10}^{-1}$, 7.86·${10}^{-3}$ and −1.81·${10}^{-2}$ and a $r_{F\left(X2\right)|F\left(X1\right)}^2$ of 0.965, 0.950 and 0.986 are observed for basal sliding, inner slope stability and outer slope stability, respectively, and a combined $r_{F\left(X2\right)|F\left(X1\right)}^2$ of 0.985 across all the failure mechanisms. Based on these small errors and the high correlation coefficient, the drainage depth, lower layer thickness and dike crest will be fixed. Furthermore, most points in Figure 6B are above the 1:1 line, indicating that fixing these input factors either results in equal or lower factors of safety, which suits our goal of performing a sensitivity analysis under the most critical conditions and justifies the use of the global minimum as the fixing point.

3.2.2. Delta Moment Independent Measure

The eight sensitive parameters are selected for the global sensitivity analysis, which results in 390,625 combinations for each failure mechanism. The DMI method (Figure 7), based on differences in output probability density functions, clearly indicates that for each mechanism, the dike slope ($D_s$) is most influential. The other dike parameters, namely its height ($D_h$) and material properties ($D_{typ}$), are also among the more important factors. Furthermore, some input factors have a different importance for the different failure mechanisms. For example, the foreland width $F_w$ mainly influences the outer slope stability, while the spacing of drainage ($D{r}_s$) mainly influences the inner slope. Other important differences include the relatively high influence of the $U_{typ}$ for basal sliding, as the sliding surface at the dike base is in continuous contact with this upper subsurface layer. Two remarkable differences are the smaller influence of $D_{typ}$ on outer slope stability, and the high standard deviation related to upper layer thickness regarding the inner slope stability (Figure 7). While these results indicate the importance of several input parameters, they do not provide further information about their relation to dike failure and local variability.

3.3. Trend Identification

To identify the trend, the probability of an unstable dike $p\left(F\le 1\right)$ given a certain input factor value is determined (Figure 8). For basal sliding, only a small fraction of the input factors’ combinations result in unstable dikes, making it the least important process of dike failure. For inner and outer slope stability, the effect is often similar, although critical safety factors are more often found for inner slope stability. Figure 8 reflects many of the DMI results (Figure 7), for example, a large influence of dike slope and dike height, where less steep sleep slopes and lower dikes are, in general, more likely to remain stable, and the small influence of $U_{thck}$. However, new insights become apparent too, as it is shown that an increase in foreland width ($F_w$) increases the outer slope stability only if the dike is less than 25 m away from the river. A similar effect is seen regarding drainage spacing ($D{r}_s$), which mostly decreases the inner slope stability at small distances from the dike (<$\sim$10 m). The most important local variability, however, is observed for the material types, where a unidirectional change in the input factor does not necessarily result in a unidirectional change in dike stability. Where any shift towards a sandier subsurface will, on average, still result in a smaller $p\left(F\le 1\right)$ for outer slope stability, this is not the case inner slope stability, which has its lowest $p\left(F\le 1\right)$ at $L_{typ}$ sandy loam. An even more striking trend is observed for $D_{typ}$, where the lowest failure probability $p\left(F\le 1\right)$ for outer slope stability and the highest failure probability for inner slope stability both coincide at a sand dike. Although basal sliding results in failure less often, a similar local variability in the response can be observed.

3.4. Subsurface Interaction Qualification

To explain these results in more detail, the response surfaces of the average dike stability $F$ at a given combination of material properties is compared with the corresponding pore pressure (Figure 9). The results are in line with the $p\left(F\le 1\right)$, which, for example, shows that sand as $D_{typ}$ results in a low $F$ for inner slope stability, but reaches a high $F$ for outer slope stability (Figure 9). The results also show that low stability is strongly linked to high pore pressures, in addition to material parameters related to strength, such as cohesion or friction angle. This is mostly the case for basal sliding, where the low factors of safety are reached at both a clay dike/sand cover layer and a sand dike/clay cover, while high pore pressures are only observed in the latter situation.

However, for inner slope stability, there is a clear coincidence of high pore pressures, sand dikes and low $F$ values. As shown in Figure 5, the dikes are least stable at $t_0$, which, due to the initial steady state, shows maximum groundwater heads for a dike type of sand, as it is more permeable. Sand dikes generally have low cohesion and rely largely on the frictional strength, but friction is lost as the high pore pressures reduce the effective normal stress. As the largest part of the slip surfaces intersects the dike and only slightly touches the cover layer, subsurface types $U_{typ}$ and $L_{typ}$ only have a minor influence on the dike’s stability on the inside.

This is not the case for outer slope stability, which often finds its most unstable condition at the end of river water drawdown. In this case, the dike’s stability is mostly dependent on the hydrological conditions, as lower stability coincides with higher pore pressures. Strikingly, low stability occurs most often for $D_{typ}$ sandy loam, for which conductivity is high enough to become fully saturated during prolonged high river water levels, but low enough to prevent rapid drainage when the river level falls. This effect becomes more prominent when the material of the dike is more permeable than that of the underlying layers; therefore, the excess pore pressure during river level drawdown cannot dissipate to the underlying layers.

There are no apparent differences between shallow and deeper subsurface layers, although the explanation differs. In case of a shallow subsurface layer ($U_{typ}$), the lower permeability directly inhibits drainage into this layer, while an impermeable deeper subsurface layer ($L_{typ}$) inhibits flow to the lateral drainage channels installed behind the dike.

4. Discussion and Practical Application

4.1. Global Sensitivity Indices for Groundwater Induced Dike Failure

Using the relations derived from the MC simulations, a database is constructed including the input factors and related factor of safety. For each failure mechanism, 390,625 parameter combinations are in this database, of which $p\left(F\le 1\right)$ was 0.02, 0.443 and 0.379 for basal sliding, inner slope stability and outer slope stability, respectively. It follows that our assumptions of an infinitely long period of high-water level at the dike crest followed by a rapid drawdown both favor inner and outer slope failure. Nonetheless, the inner slope is generally more unstable, as its minimum stability is reached during the infinitely long high water, whereas outer slope failure only occurs during the drawdown. This rapid but transient drawdown decreases pore pressure, which thus increases friction and results in more stable conditions than on the inner slope. The larger fraction of simulations in the database with inner slope failures matches with practical experience [44], which suggests the modelling assumptions in this paper result in groundwater conditions that are representative of real-world conditions. Combining all the three mechanisms results in a total $p\left(F\le 1\right)$ of 0.48 out of all the combinations, i.e., based on simulations that result in $F\le 1$ for one or more mechanisms.

All the parameters involved in the analysis have at least a small effect on the dike stability. In general, the least influential parameters are $F_w$, $L_{thck}$ and $D{r}_d$. However, the foreland width becomes influential at small values (Figure 8). Furthermore, the inner and outer slope stability are insensitive for changes in the dike width given the steady-state groundwater conditions, while the dike width does affects basal sliding. In addition, drainage spacing influences inner slope stability mainly if drainage occurs close to the dike. As such, many parameters have a considerable effect on a part of their parameter space and a given failure mechanism, but few are important over their entire parameter space and for all the failure mechanisms. This suggests that although global sensitivity analysis is a useful tool in determining the importance of parameters, the resulting sensitivity indices should be handled with care. In some cases, sensitivity indices such as the DMI [29] underestimate the general importance of a parameter, as they provide a mean parameter sensitivity that does not necessarily reflect the abundance of strong local effects.

4.2. Toward Flooding

The probability of flooding is not only related to the probability of dike instability, but also to the water levels at the time of failure and the volume of the instability. To analyze the occurrence of these properties, we select the set of model runs for dikes with a 1:3 slope, a height of 6.5 m and a cover layer thickness of 1.5 m. Of the 625 selected scenarios that remain in the database with these parameters fixed, the slip surface shape and phreatic water level are related to the upper substrate type ($U_{typ}$) for outer slope stability (Figure 10). In this example set, the area of failure decreases from 59.2 to 26.4 m², while the mean phreatic level in the upper subsurface layer decreases from 2.36 to 1.93 m when the substrate material changes from clay to sand. Accordingly, the average safety factor increases from 1.06 to 1.29. This indicates not only that dike failure as a result of macro-instability is less likely due to occur with a sand subsurface, but also that potential instabilities are less threatening for a dike breach, as their volume is much smaller.

It should be noted that the groundwater model assumes a stationary response to the flood wave at the dike crest, but a transient response on the drawdown of this flood wave. From a hydrological perspective, the simulations thus start at a situation of minimal internal strength and, hence, minimal dike stability. However, as the high river water levels act as a stabilizing external force on the outer slope, this aspect of failure generally occurs after river water level drawdown. The breaches that subsequently occur create vulnerable dikes but, owing to the falling flood water levels, this is unlikely to lead to major dike beaches and flooding. Obviously, if this instability is followed by a second flood wave, the situation might become critical. Thus, additional research is needed in the transient dike response under very common multi-peak flood waves [45].

4.3. Limitations Regarding Sampling, Hydrology and Subsurface Uncertainty

First, this research used a uniform sampling strategy, not taking into account any possible correlations. We suggest, however, that by sampling each parameter uniformly over a range of possible outcomes and not taking account of possible but unknown correlations between parameters, we are likely to overestimate the possible range of outcomes. This conservative estimate of possibilities is deemed a rational choice, if no information is available about the a priori joint probability distribution of parameters.

Second, we used the most adverse hydrological conditions for dike slope stability calculations, being steady-state conditions for inner slope stability. We acknowledge that this might not be the case when considering flood duration and transient groundwater conditions, as the infiltration curve might not reach the inner toe and the safety factor becomes more favorable. A similar argument applies to basal sliding.

Third, this study used many possible scenarios of subsurface material and surface geometry combinations. Still, the subsurface material in this study is assumed to be homogeneous and deterministic in each of the layers, and likewise the layer thickness and surface profile are constant per section. Due to the long history and continuous improvement of many dikes, their interior is presumably very heterogeneous [6,7]. In reality, the subsurface properties, induced by previous river systems, are also known to have a large spatial variability [46,47]. These heterogeneous topographic and subsurface properties have a large influence on both hydrological conditions [48] and stability [49], while still ignoring the 3D slope effects [50]. In addition, animal burrows and human measures inhibiting or enhancing groundwater flow may also affect the pore pressure evolution.

4.4. Suggestions for Further Research

Nonetheless, the simulations show that for all the mechanisms, the dike type and upper layer type have a large effect on the dike stability. Despite their importance for dike stability and the valid assumption that they are heterogeneous, these two subsurface parameters are often partly unknown. As such, extending this research with variable material properties [11] may provide an even more extensive analysis of groundwater-related dike stability. In this way, the uncertainty of the material properties can be assessed too, which is important in real-world cases.

In addition to assessing the uncertainty, decreasing it is of major importance when assessing dike stability [51]; this is achieved by having more detailed subsurface data available. Thus, incorporating large scale subsurface heterogeneity (in 3D) is an important step in actively incorporating groundwater calculations in dike stability calculations, although we are confident that the large contrast in dike and subsurface materials used in this study, to a great extent, covers the large variation across many dikes.

Furthermore, additional research is needed in the transient groundwater response given a certain flood wave instead of steady-state conditions. Finally, the database constructed and explored in this paper could be used for mapping those regions where factors of safety might reach critical values.

4.5. Case Study: Application of the Database for Fast High-Resolution Dike Safety Assessment

As a first attempt for identifying those regions, the created database was applied to a real-world case. This case concerns the area near the village of Ameide in The Netherlands (51.954594 N, 4.963298 E; Figure 11), for which an official preliminary assessment of the dike stability was made [52]. The failure probability of the dike was calculated using the Bishop’s modified method [53], and based on the characteristic values of material strength parameters. The phreatic level was simplified in the schematization and was roughly a straight line between the river water level and the ditch water level. In the assessment, larger segments with reasonably similar characteristics were tested against a failure probability of 1/360 years and assigned the final judgement. In a later stage, the precision of these judgements was drastically improved using local schematizations and locally derived characteristics for the important parameters; hence, the current assessment was seen as preliminary.

4.5.1. Case Study Methodology

The dike in this area was subdivided in segments with a length of 100 m, of which several failed to meet the expected failure probability. Here, we compared the values of the factor of safety for comparable situations in terms of dike geometry and composition from the Monte Carlo simulations with the assessment in order to test whether unsafe conditions were revealed by our approach. We hypothesized that using a pre-constructed database with factors of safety in combination with actual dike characteristics could provide a quick a priori analysis of the dike stability. To compare the official assessment with the database, those dike parameters should be selected from a database that corresponds to the actual dike. The parameters concerning geometry and composition of the dike and the subsurface were determined at an interval of 10 m along the entire dike crest. The dike height, crest width and slope were automatically derived from the high-resolution lidar-based AHN3 surface elevation model. The properties of the subsurface, being layer thickness and lithology, were derived from the GeoTOP subsurface model [54]. An approximation of the dike material was made from publicly available cone penetration tests (BRO) using a simple but effective method proposed by [55]. When these parameters were assembled, the corresponding safety factor was selected from the database.

4.5.2. Case Study Results

First, as the acquisition of surface geometry and subsurface composition can have a high spatial resolution of 10 m, it provides more detailed information than the official assessment, which is conducted on a 100-meter resolution. On visual inspection, the calculated safety factors already clearly coincide with the official dike assessment (Figure 11B), although the variation of the calculated values is much higher, as the official safety assessments are carried out only per 100 m section, whereas the factor of safety is calculated here every 10 m. There is also a clear difference visible between the inner and outer slope stability. Those dike sections that are assessed as sufficiently safe have an average safety factor of 1.53 ± 0.31 and 1.44 ± 0.13 for their inner and outer slope stability, respectively. The insufficiently safe sections according to the official assessment have a calculated average safety factor of 1.41 ± 0.38 and 1.34 ± 0.20, respectively. Thus, according to our method, most of the sections found to be insufficiently safe have an F $>$ 1 (Figure 11A). In the official assessment, 69.3% of the dike was found insufficiently safe against 10.8% in our analysis. Of this 10.8%, 9.5% is found on the insufficiently safe sections and 1.3% is found on those sections that were found to be safe enough.

4.5.3. Case Study Discussion and Conclusions

These false positive assessments, all occurring on the inner dike slope, are most dangerous. Their lower safety factors are caused by a combination of steeper dike slopes and the presence of sandy dike material. Thus, in addition to showing high spatial variability in the expected factor of safety, the analysis also clearly shows those sections that, according to our calculations, are the most critical. These differences seem to be largely related to variations in dike material but can also be related to some of the parameters (drainage, dynamic river level) that are not included in our analysis. Moreover, it is likely the cause of different definitions of failure and the use of different stability calculation methods, such as drained or undrained loading [39]. In conclusion, the high-resolution database comparison could help to focus further research and data assembly, by indicating areas to improve local schematizations and derive local characteristics for the important parameters in a more detailed stage of dike reinforcement design.

5. Conclusions

In this study, an extensive global sensitivity analysis was carried out for dike stability using calculated groundwater heads and resulting pore pressures that represent a worst-case scenario. The following three sensitivity settings formed the basis of this analysis: parameter prioritization, trend identification and interaction qualification. The results show that each of the three studied failure mechanisms, namely, basal sliding, inner slope stability and outer slope stability, can possibly result in dike failure, where failure on the inner slope of the dike is the most likely.

In the parameter prioritization settings, eight parameters were determined to be influential for any of the failure mechanisms, being dike slope, dike height, dike material type, upper layer material type, lower layer material type, upper layer thickness, drainage spacing and foreland width. In contrast, the dike crest width, drainage depth and aquifer thickness had a negligible effect on the stability. According to the Delta Moment-Independent measure (Equation (4)), the dike slope was found to have the largest effect on dike stability, as it has both a direct effect on slope stability (Figure 5) and indirectly affects the slope stability by changing the pore pressure conditions (Figure 7). We conclude that the delta sensitivity index provides a clear indication of parameter importance. The Delta Moment-Independent measure does not provide information on parameter interaction and the underlying mechanisms.

A regional sensitivity analysis by means of the probability of failure $p\left(F\le 1\right)$ was used to show local and process dependent variability. This measure shows that dike slope and dike material are most influential for resulting in unstable dikes. Furthermore, it indicates local variability, such as that drainage spacing only affects the inner slope stability if the drainage location is close to the dike. Combining the Delta Moment-Independent measure indices and the $p\left(F\le 1\right)$ showed that while basal sliding is most sensitive to changes in most parameters, it is least likely to result in dike failure.

In addition, we qualified the interaction between different material types, hydrology and stability. Most strikingly, a combination of a relatively permeable dike and impermeable subsurface inhibits the dissipation of pore pressures to the lower layers. Therefore, pore pressures remain high and dike instabilities are more likely to occur. This effect is more prominent for outer slope stability, as dissipating high pore pressures is mostly important during river level drawdown. The area of failure is often also larger in case of an impermeable subsurface, increasing the chance of severe flooding after a slope instability.

Applying the database containing geometry parameters, subsurface properties and safety factors to a real case resulted in high-resolution estimates of dike stability. These estimates show that the unsafe segments derived from the database are mostly on those segments also found to be unsafe by the official assessment, although there are some false positives, e.g., segments that are estimated to be unsafe but classified as safe in the assessment. Overall, the database estimates provide a more differentiated picture and allows for more targeted analyses and measures. Although three-dimensional subsurface buildup and variable flood waves can improve the simulation results, the comparison gives confidence that our results provide useful insights in the process of groundwater-related dike stability and that the underlying database can be used to focus additional local research. The analysis of groundwater-related dike failure with global sensitivity methods clearly shows the importance of high-resolution groundwater modelling for estimating dike slope stability.