Evaluating residual dyke resistance using the Random Material Point Method

Abstract Due to a lack of large deformation dyke assessment models, primary failure mechanisms, such as inner slope failure, are often used as a proxy to assess the probability of failure of a dyke. However, a dyke continues to fulfil its main function unless, or until, flooding occurs. The Random Material Point Method (RMPM) is used here to investigate residual dyke resistance, which is the resistance against flooding after initial failure. RMPM combines random fields with MPM in a Monte Carlo simulation and has been extended here to include the effects of an external hydrostatic pressure on a dyke’s outer slope. The residual resistance of an idealised dyke (computed using RMPM) is shown to reduce the probability of flooding by 25% with respect to the initial failure. A lower degree of anisotropy of the spatial variability increases the residual dyke resistance. RMPM simulates, as expected, a lower residual dyke resistance for larger initial failures and/or a higher water level. A ‘safe’ remaining geometry has not been found, since even small initial failures can result in an unacceptable probability of flooding, highlighting the importance of modelling the entire failure process.


Introduction
"What is failure of a structure?" is a simple question to ask, but can be difficult to answer in practice. Failure of a structure is defined as the structure becoming unfit for service, which either means a loss of safety under extreme circumstances, i.e. the Ultimate Limit State (ULS), or a loss of function under normal circumstances, i.e. the Serviceability Limit State (SLS) (ASCE, 2005;CEN, 2002). To define ULSs and SLSs for a flood protection structure, the consequences of flooding must be considered. Indeed, the potential damage to the flood protection itself is usually insignificant in comparison. So, while the Eurocode defines a loss of equilibrium of the structure or significant damage to the structure as ULS failure, for flood protections in particular, a loss of equilibrium or significant deformation can be accepted so long as the structure retains its primary function, i.e. preventing flooding of the hinterland (CEN, 2002(CEN, , 2004. In such circumstances, a loss of equilibrium and/or significant deformation of a flood protection can be considered as SLS failure. In the Netherlands, the Water Law (Waterwet, 2009) has defined ULS dyke failure as flooding. Therefore, for reliability assessments, the probability of ULS failure of a dyke is given by the probability of flooding. For most failure processes, the initiation of the mechanism is generally considered by guidelines to always lead to flooding (MIM, 2017a). For example, to determine the probability of flooding due to inner slope failure, i.e. failure of the landward slope of the dyke, the loss of equilibrium is usually considered as the limit state. However, several theoretical case studies and dyke failure tests have shown that after macro-instability flooding can still be prevented, i.e. residual dyke resistance, or 'reststerkte' in Dutch, was available (ENW, 2009; van Montfoort, 2018;'t Hart et al., 2016;Klein Breteler et al., 2009). For other mechanisms, such as (wave) overtopping or piping, residual dyke resistance is also often available (Calle, 2002). Hence, including this additional resistance can lead to more efficient designs, and to dyke reinforcement taking place only where it is most necessary.
Some guidelines acknowledge the existence of residual dyke resistance and propose methods to include this additional resistance, especially for wide dykes, i.e. dykes where the critical failure surface only affects either the inner berm or slope, and dykes with a large foreland (Blinde et al., 2018;ENW, 2009;MIM, 2016). Almost all of these methods reduce the probability of flooding (or increase the factor of safety) compared to the probability of inner slope instability based upon the remaining dyke geometry. However, since existing methods for the assessment of inner slope instability, for example the Limit Equilibrium Method (LEM) or Finite Element Method (FEM), generally stop at the start of failure, both the remaining geometry and potential secondary mechanisms are evaluated based on crude approximations. This paper further develops the understanding of, and modelling capabilities for, the failure processes using the Material Point Method (MPM), a tool that is able to model the entire failure process. MPM is a large deformation variant of FEM, that makes use of both a material point discretisation for storing the information of the material and a mesh discretisation on which the governing equations are solved. MPM has been used to model large deformations for a range of geotechnical applications (Andersen et al., 2010;Fern et al., 2019;Martinelli et al., 2017;Phuong et al., 2016;Wang, et al., 2016b), and is able to model the entire dyke failure process, from initial macro-instability to flooding (Fern et al., 2017;Leclercq, 2020;Zabala and Alonso, 2011;Zuada Coelho et al., 2019).
A further consideration in dyke stability assessments is the variability of soil properties, which can have a large impact on the likelihood of initial failure as well as on the failure process. The Random Finite Element Method (RFEM) combines random fields, for modelling the spatial variability of soil properties, with FEM, and has been used successfully to quantify the impact of variability of soil properties on the probability and geometry of inner slope instability (Griffiths & Fenton, 2004;Hicks & Li, 2018;Hicks & Samy, 2002;Hicks & Spencer, 2010;van den Eijnden & Hicks, 2018). Similarly, the Random Material Point Method (RMPM) combines random fields with MPM to investigate the impact of spatial variability on the entire failure process. This paper extends the implementation of RMPM, as presented in Wang et al. (2016a) and Remmerswaal et al. (2019b), to make it applicable to dyke safety assessments incorporating residual strength. Section 2 describes an idealised dyke cross-section that has been used throughout the paper to demonstrate the approach. In Section 3, a probabilistic framework is presented for assessing residual dyke resistance, and this framework is implemented using RMPM in Section 4. Sections 5 and 6 present a parametric investigation into the probability of initial failure and subsequent flooding.

Example dyke
An idealised dyke cross-section is used throughout, to demonstrate the capabilities of the method. The dyke, shown in Fig. 1, has an initial height H i = 5 m, an inner (right-hand-side) slope of 1 in 3 and an outer (left-hand-side) slope of 1 in 1. The outer slope angle is steeper than that of a usual dyke, in order to increase the necessity of an external water load on the outer slope for maintaining stability and thereby test the veracity of the external loading implementation. In addition, the width of the crest W c = 10 m is relatively large in order to more fully explore the effect of residual dyke resistance. The water level h is varied and the dyke is underlain by additional soil layers. Fig. 2 shows potential failure mechanisms related to macro-instability. A deep rotational slide is shown in Fig. 2a, where flooding occurs immediately. Deep slides usually result in limited residual dyke resistance, due to the size of the failure. Conversely, residual dyke resistance is usually seen when the dyke is founded on a stiff layer (as shown in Fig. 2b-e). In this case, a shallow rotational slide (as shown in Fig. 2b) causes a reduction of the width of the dyke W (defined at the water level), while the maximum height of the dyke H remains unchanged. A rotational slide can be followed by a retrogressive failure (Fig. 2c), which may or may not result in flooding. Multiple successive retrogressive failures can also lead to flooding. The external water loading can also push the entire dyke section sideways to cause a horizontal sliding mechanism (Fig. 2d). The height of the dyke tends to slowly decrease during such a slide, while the width stays more or less constant. In the longitudinal direction of the dyke, a horizontal slide may only occur along a small stretch, so that water may flow through a gap created by the slide. However, since the simulations performed here are 2 dimensional, this 3-dimensional effect is neglected. Finally, especially for lower water levels, failure of the outer slope may occur, see Fig. 2e, after which retrogressive failure in the opposite direction is possible. Combinations of the previously mentioned failure mechanisms may also occur; for example, a horizontal sliding mechanism after inner or outer slope instability.

Risk framework
To design or assess dykes in a risk based framework, both the probability of failure P(F) of a dyke and the consequence of failure C(F) are required, since the risk of failure is defined as  Calle (2002) and 't Hart et al. (2016).
Since flooding is usually considered as the dominant limit state for dykes, the probability of flooding and the consequence of flooding should be considered. This paper focuses on computing the probability of flooding, since the consequence of flooding is often included separately in safety standards. Therefore, a maximum allowable probability of flooding P s , based upon the minimisation of the combined cost of dyke reinforcement and the consequence of flooding, is defined: The major failure categories of a dyke are considered to be overtopping, piping, macro-instability of the inner slope, damage to the cover of the outer slope and failure of structures which are part of the dyke (MIM, 2017b). Assuming, conservatively, that all failure categories are mutually exclusive, the probability of flooding can be computed by summing up the probability of occurrence of the failure categories. Herein, to demonstrate MPM, only macro-instability is considered, although the risk framework concepts can be extended to other failure categories.
A limit state function Z can be defined, by considering the resistance against failure R and the loads triggering failure L (Schneider, 2006). Hence, for failure to occur, so that the probability of flooding becomes R and L are modelled and/or measured, and both may contain model uncertainty (M), material (or 'quantity') uncertainty (Q) and geometric uncertainty (D), i.e.
The material uncertainty can be further split into measurement uncertainty and transformation uncertainty where indirect measurements are taken. In addition, the resistance and loads may be functions or models that are dependent on further variables, each of which may have further uncertainties. Combining Eq.
(3) and Eq. (5), Z becomes a function of all the uncertainties: Geometric uncertainties D R and D L tend to be limited, due to the relatively high measurement accuracy, and are, from here on, considered to be negligible (see Varkey et al. (2021) for an investigation of the relative importance of geometric uncertainties in dyke assessments). Moreover, model uncertainties are ignored, since they have not yet been established for MPM. Hence, this paper focusses on uncertainties relating the spatial variability of material properties and uncertainties relating to the external loading (due to the external water level), i.e. Q R and Q L , respectively. The limit state function Z for macro-instability, and most dyke failure categories, is usually defined according to a loss of equilibrium, i.e. the start of slope instability, for which the resistance R and load L can be determined using for example LEM or FEM. However, to include residual dyke resistance, the process from slope instability to flooding should be evaluated, and Z must be defined accordingly.
A series of events leading to a flood always ends with a breach event, i.e. erosion of the dyke by flowing water (Calle, 2002). This occurs once the maximum height of the dyke drops below the external water level. Therefore, the limit state function Z can be defined as where the dyke height H can be seen as the resistance R against flooding and the water level h as the load L. The probability of flooding is therefore given by After a primary macro-instability failure, flooding may occur directly ( Fig. 2a) or by way of secondary mechanisms; e.g. retrogressive failure (Fig. 2c), internal erosion (i.e. micro-instability), or wave-overtopping. For simplicity, in this paper the occurrence of wave-overtopping or internal erosion after macro-instability are not considered, and only retrogressive failure, the secondary mechanism most closely related to the primary mechanism, is considered.
Considering an uncertain water level h, Eq. (8) may be expanded using fragility curves, i.e. probabilities of flooding given a specific water level P(F|h) (Simm et al., 2009), such that where f(h) is the probability density function of h given the hydraulic conditions imposed on the dyke, i.e. all conditions influencing the external water level, such as wind, upstream conditions, tidal currents, and so on. The integral can be numerically approximated from the calculation of discrete values and interpolation, i.e. the probability of flooding can be computed separately for specific water levels, weighted using the likelihood of the water levels, and then summed across the water levels.
By evaluating all possible macro instability events, the probability of flooding given a specific water level can be evaluated as where S is one of the possible failure processes, which may or may not lead to flooding. Each failure process S may be split into a sequence of separate failure mechanisms, starting with an intact dyke S 0 affected by a primary macro-instability failure S 1 , followed by a potentially infinite number of successive failures S 2 … S ∞ . In recent research (van der Krogt et al., 2019), Eq. (10) has been evaluated assuming that the probability of flooding can be computed by following the most critical failure process S c , i.e. only evaluating the most likely failure mechanisms S c 1 ⋯S c k using LEM. This assumes that the sequence of critical failure mechanisms is a good representation of all possible failure mechanisms and thereby the probability of flooding. However, sub-critical failure mechanisms might have a larger probability of flooding, as the failure mechanisms may be larger, reducing residual dyke resistance (van der Krogt et al., 2019). In some cases, a sub-critical primary failure mechanism might even lead to flooding directly. To avoid assuming critical failure mechanisms, Eq. (10) is here modelled completely, i.e. all possible failure mechanisms are considered.

Modelling method (RMPM)
The entire failure process is here modelled using the Material Point Method (MPM). Spatial variability of the material properties is considered by way of random fields, and an ensemble approach is utilised such that uncertainties are included via use of the Monte Carlo Method, resulting in the so-called Random MPM (RMPM).

Material Point Method
MPM was proposed by Sulsky et al. (1994) as an extension of the Particle in Cell method, and can also be considered an extension of the Finite Element Method (FEM). In this method, the integration (sampling) points are replaced by material points (MPs) which can move freely through the background grid. MPM can model large deformations and strains, without encountering the mesh distortion problems associated with FEM. Wang et al. (2016b) showed the method to be capable of modelling retrogressive failure, similar to the failure process observed in dykes, in (shear) strain softening clay materials. In this study, a single phase, total stress, plane strain MPM is used. The undrained shear strength softening model described by Wang et al. (2016c) is used to represent strength reduction due to pore pressure build up during failure.
Implicit MPM is used, so that larger time steps can be used to reduce the computation time (Wang et al., 2016c). In the standard MPM formulation, stress oscillations occur as MPs move through the mesh and especially when they cross element boundaries. In order to limit such oscillations, the DM-G technique has been used in this paper (González Acosta et al., 2020Acosta et al., , 2021. This method uses extended shape functions (and shape function gradients) to reduce the cell crossing error, as well as FEM Gauss point stiffness integration to reduce the error due to movement within the background grid cells. The Gauss point properties are acquired by the mapping of material point properties via the background grid nodes.

Random Fields
Random fields are used to model the spatial variability of soil properties, and have herein been generated using Local Average Subdivision (Fenton & Vanmarcke, 1990). The random field values are mapped onto the MPs in their initial positions at the start of the analysis, and the material properties of each MP remain constant throughout a single realisation. Spatial variability has only been considered for the shear strength properties of the material, with all other material properties being assumed deterministic. Specifically, the initial and residual undrained shear strengths, c i and c r , respectively, are considered to be spatially variable. They are assumed to be fully correlated, so that the spatial variability of both properties can be constructed from a single random field. Partially correlated random fields for different properties are also possible (e.g. Vardon et al., 2016), but are not considered here for reasons of simplicity.
The spatial variability is defined by the point statistics, i.e. the mean μ and standard deviation σ, and by the spatial statistics, i.e. the horizontal and vertical scales of fluctuation, θ h and θ v , respectively. The point statistics of the initial and residual undrained shear strengths, μ ci , σ ci and μ cr , σ cr respectively, are different, whereas the same scales of fluctuation have been assumed for both properties.

External water load
Applying boundary conditions is a relatively simple task when the nodes of the mesh coincide with the material boundary. However, the application of external loads in MPM is more difficult, since the location of the boundary is (in general) not known a priori (Bing et al., 2019;Cortis et al., 2018;Remmerswaal, 2017). As the boundary represents the surface (or edge) of the material and the outer layer of MPs are, by definition, locally at the centre of the material they represent, the boundary cannot simply pass through the outer layer of MPs. Hence, the boundary of the dyke must be located outside the MPs, and herein it is located using boundary detection (Remmerswaal, 2017).

Boundary detection
The boundary is detected using the Proximity Field Method (PFM) (Remmerswaal, 2017;Remmerswaal et al., 2019a). PFM, which is based on the level set method, assigns kernel functions to each MP to represent its influence domain (Sethian, 1996). The standard level set method describes the boundary location as the zero level set of an auxiliary field (Sethian, 1996). Here, the proximity, i.e. distance, to nearby MPs is used as the auxiliary field. This proximity field (PF) is created by summing the kernel functions assigned to the MPs, where each kernel function represents the proximity to its MP, see Fig. 3a. The PF is evaluated only at the background grid nodes, by summing the kernel functions at each node as shown in Fig. 3a. The standard level set method would define the boundary as the zero level of the PF. However, for a smoother boundary a higher threshold is beneficial. Therefore, the boundary is located by determining the locations on the grid edges where the PF is equal to a threshold value, i.e. points belonging to the boundary are determined (Fig. 3b).
The size of the kernel functions is generally larger than the material volume they represent, so as to facilitate a smooth boundary. In this paper, Epanechnikov (parabolic) kernel functions have been chosen (Fig. 3a), which have a circular basis in the initial position in 2 dimensions. The strains of the MPs can be used to transform the circular kernel functions into ellipsoids to improve the boundary detection (Remmerswaal, 2017), but this has not been used here. In 1 dimension, the boundary points are the boundary; i.e. they denote the 2 ends of each 1D domain as illustrated for the 2 1D domains in Fig. 3b. In 2 dimensions, the complete boundary is created by connecting the points using Composite Bézier Curves; i.e. a set of connected polynomial boundary segments C contained within background grid cells (Fig. 3c). Multiple sets (C 1 until C m ) might be needed if the dyke splits into multiple pieces or holes occur within the dyke after deformation. Here linear boundary segments are used, which usually give a good representation if the mesh discretisation is small compared to the size of the dyke (Fig. 3c).

Application of external load
Hydrostatic pressure is applied on the boundary segments from left to right, i.e. from the outer side towards the inner side of the dyke, until one segment is above (or partially above) the water level (Fig. 3c). Gauss integration is used along the segments to integrate the hydrostatic pressure (Bing, 2017;Bing et al., 2019;Remmerswaal, 2017), as illustrated in Fig. 3c for one segment (which extends from y = -2 m to y = -1.2 m). The hydrostatic pressure and the direction of the applied pressure at the Gauss point are computed using the depth below the water level and the direction normal to the boundary, respectively. Gauss integration along the boundary allows the pressure to be converted to external nodal loads.

Detection of S i and F from residual dyke geometry
In each time step, the external geometry is determined using the boundary detection method. Two values are extracted; namely, the width W and height H, as shown in Fig. 2. H is the maximum height of the dyke, whereas W is the width of the dyke at the water level h. H(t) and W(t) are used to detect the moments when the initial failure mechanism S 1 and flooding F occur. After an initial slide, W will be smaller than the original width of the dyke, i.e. W < W i as shown in Fig. 2a and 2b. Moreover, W decreases to 0 and H < h when flooding occurs (Fig. 2c). The residual width W r and residual height H r are the width and height at the end of the simulation. The probability of S 1 with a time of failure (t S1 ) before a time t can be computed from the realisations in a Monte Carlo simulation, i.e. P(S 1 |h, t S1 < t) ≈ N S1 (t) where N S1 (t) = N W(t)<Wi is the number of realisations where an initial failure mechanism has occurred before time t and N is the total number of realisations. Similarly, the probability of flooding with a time of flooding (t F ) before t is computed as where N F (t) = N H(t)<h is the number of realisations in which flooding occurred before t.

Analyses
The dyke presented in Section 2 has been simulated. A weak strain softening clay, with a mean initial undrained shear strength μ ci of 13 kPa and a mean residual undrained shear strength μ cr of 5 kPa, is used as the dyke material. These material properties were chosen to represent a dyke with a relatively high probability of failure, in order to enable reasonably quantitative results within a manageable number of realisations for investigating the deformation processes during dyke failure. For the efficient computation of smaller, more realistic, probabilities of initial failure and flooding, RMPM can be combined with subsetsimulation as has already been done successfully with RFEM ( Table 1 and Table 2, respectively. A mesh of 4 noded square elements and an element size Δx = Δy of 0.25 m is used, together with 4 equally spaced MPs per element (giving 6400 MPs in total). Along the sloping faces of the dyke, the number of MPs per element is adjusted and the MPs redistributed, in order to match the geometry. The nodes along the bottom boundary are completely fixed, so as to model a dyke resting on a firm foundation (as illustrated in Fig. 2b to 2e). The random field is generated with a cell size corresponding to the material point domain, i.e. ½Δx by ½Δy, such that the properties of a single cell correspond to one material point. Hence, as the random field cell size is 4 times smaller than the smallest scale of fluctuation (i.e. θ v ), the spatial variability is adequately captured (see, for example, Huang & Griffiths (2015)).
The time step Δt is 0.01 s, which is significantly larger than would have been required using explicit MPM (González Acosta et al., 2020Acosta et al., , 2021. The maximum simulation time t max is 40 s. The in-situ stresses are generated quasi-statically using gravity loading, assuming an elastic material, whereas plasticity and dynamics are considered during the simulation itself. At the start of the simulation, the plasticity generates unbalanced forces at the nodes, which may lead to an initial failure. To reduce the vibrations caused by suddenly switching on plasticity at the start of the simulation, damping is used. The unbalanced forces are artificially reduced (damped) by 100% at the start and then gradually increased until they have been fully applied after 1 s; thereafter, no damping is used. The realisation is terminated if (a) flooding occurs, (b) the dyke does not fail at all, or (c) the dyke, having experienced slope  Each simulation comprises 10 000 realisations, which have been performed in parallel on a new grid computing system called Spider, at SURFsara, a national computing centre in the Netherlands. This number is higher than for typical RFEM simulations, due to the greater range of possible failure mechanisms arising from the evolving failure process being modelled in the simulations.

Analysis 1: Base Case
The final dyke geometries for 6 realisations of Analysis 1 are presented in Fig. 4. Realisations without any failure are by far the most common type, as illustrated by the example in Fig. 4a. In contrast, Fig. 4b shows an initial failure mechanism, initiating through a weak layer approximately 1.5 m above the base. Due to the weak layer, a slip circle is formed, along which the soil reduces in strength due to strain softening as the failure progresses, as indicated by the dark slip surface. Moreover, due to the large deformations, the sliding mass breaks up into 2 distinct blocks separated by a secondary failure surface inclined at approximately 45 • with the slip circle. Since the depth through which the failure surface forms is limited, the initial failure does not develop further and reaches the stable configuration shown in the figure. Fig. 4c presents the most likely initial failure mechanism, which occurs along a weaker layer close to the base of the dyke and resembles failure in homogenous material. Again, the sliding mass breaks into 2 blocks separated by a secondary failure surface at approximately 45 • . The figure shows some softening in the weak layers directly behind the slip circle, but enough resistance remains to prevent a second (retrogressive) failure and the possibility of flooding.
Retrogressive failure does occur in Fig. 4d and 4e. After an initial slide, a second slide fully develops in both figures; in addition, a third slide partially develops in Fig. 4d and fully develops in Fig. 4e. Once again, breaking up of the sliding masses can be observed. While the height of the dyke is almost unaffected in Fig. 4b and 4c, the height of the dyke reduces in both Fig. 4d and 4e. The final height of the dyke in Fig. 4d is slightly higher than the water level and, in reality, the dyke may flood due to wave-overtopping. However, under the current modelling assumptions, the dyke remains stable and flooding was prevented. In contrast, a flood is triggered in Fig. 4e. In reality, a breach would then occur due to erosion, causing further (catastrophic) damage to the remainder of the dyke; but, since the model immediately stops after flooding, the geometry in Fig. 4e is the final configuration.
A different type of mechanism is shown in Fig. 4f, where an initial rotational slide triggers a horizontal secondary slide along a very weak layer near the base of the dyke, eventually resulting in flooding. This highlights that RMPM is capable of approximating both rotational and horizontal sliding mechanisms. It is seen that the leading tip of the primary slip surface in Fig. 4b-f is suspended above the remaining slope or foundation, most likely due to the element discretisation adopted in MPM. This error may be removed by including a contact algorithm between bodies (González Acosta et al., 2021). However, in these analyses this effect only occurs towards the end of the simulation and therefore does not affect the results.
The probability of failure, i.e. the probability of initial failure P(S 1 ) and flooding P(F), of Analysis 1 over time and the probability of failure before a given time are shown in Fig. 5. The probability of flooding is lower than the probability of initial failure (i.e. ~ 9% compared to ~ 12%), indicating that residual dyke resistance is present in this analysis. Failure occurs faster than would normally be expected due to the low residual undrained shear strength, the absence of damping and the instantaneous triggering of failure due to the loading conditions.   A deterministic MPM simulation based on the mean strength properties (μ ci = 13 kPa and μ cr = 5 kPa) did not lead to initial failure, and gave a Factor of Safety (FOS) ≈ 1.3. In a deterministic homogeneous MPM simulation with a reduced strength, such that FOS is just below 1 (c i ≈ 10 kPa and c r ≈ 3.85 kPa), failure occurs after 4 s, and flooding occurs 3.5 s later due to retrogressive failure. By integrating the probability distribution of c i until 10 kPa, the probability of FOS < 1.0 has been determined as 0.13. Therefore, the estimated probability of initial failure based only on the point statistics (13%) is similar to the probability of initial failure found with RMPM for this problem (12%). However, the difference becomes much larger with smaller values of θ h , as will be demonstrated by Analysis Set 2. Indeed, with realistic soil properties, when the probability of failure of the dyke is much lower than the one computed here, RFEM (or RMPM) generally computes a much lower probability of initial failure compared to the estimation based only on the point statistics . In this case, a reduction in the probability of flooding compared to the probability of initial failure is not observed in the homogeneous MPM simulation, since flooding occurred at the same strength as initial failure. A further reduction in the probability of flooding is observed using RMPM. Therefore, RFEM can be used to reduce the conservatism regarding initial dyke failure by accounting for spatial variability, while RMPM can further reduce this conservatism by also accounting for residual dyke resistance.
6.2. Analysis Set 2: Effect of horizontal scale of fluctuation 6.2.1. Probability of failure Fig. 6 presents the probability of an initial failure mechanism and the probability of flooding for various horizontal scales of fluctuation, i.e. the results of Analysis Set 2. Fig. 6a shows that, if it occurs, the initial failure mechanism triggers 2 to 5 s after the start of the simulation. As for the Base Case, flooding occurs later if the residual resistance is overcome, see Fig. 6b. The times over which flooding may occur are more widely spread than the times of initial failure. Fig. 6a and 6b clearly show that, as expected, the probability of flooding is lower than the probability of initial failure for all horizontal scales of fluctuation. Moreover, according to Fig. 7a, the largest absolute decrease in the probability of flooding compared to the probability of initial failure occurs at an intermediate value of θ h . However, more importantly, Fig. 7b shows that the residual resistance is highest at small horizontal scales of fluctuation, since the probability of flooding given the occurrence of an initial failure is the smallest. Subset simulation could be used to further investigate residual dyke resistance at lower probabilities of initial failure, since the geometry of the failure mechanism may then be different (van den Eijnden & Hicks, 2017). For large horizontal scales of fluctuation, the probability of an initial failure increases and the residual dyke resistance decreases, see Fig. 7. Failures are more likely to occur in weak layers in the presence of larger horizontal scales of fluctuations and retrogressive failure has a greater tendency to occur through these same weak layers, thereby reducing the residual resistance. Fig. 8 shows histograms of the residual geometry, H r − h and W r , for Analysis Set 2. The histograms in Fig. 8a and 8b exclude completely stable dykes, which all have W r ≈ 10.5 m and H r − h = 0.25 m, and flooded dykes, which all have W r = H r − h = 0, see Fig. 8c. The histograms in Fig. 8a and 8b have been normalised with respect to the total number of included dykes. Horizontal scales of fluctuation below θ h = 4 m are not presented, because the number of included dykes becomes very small.

Residual dyke geometry
Due to the limited height above the water level of the dyke, most failures with residual dyke resistance, i.e. do not flood, do not affect the height of the dyke (H r − h ≈ 0.25 m). Therefore, a peak is visible at H r − h ≈ 0.25 m (Fig. 8a). The average height above the water table H r − h is slightly higher for smaller horizontal scales of fluctuation, i.e. 0.22 m for θ h = 4 m against 0.20 m for θ h = 48 m.
The residual width tends to be smaller with larger horizontal scales of fluctuation (Fig. 8b). The approximately bi-model nature of the histograms of residual width are the result of the expected size of the initial and retrogressive failures. Since rotational failures are most likely and significantly influenced by depth (due to the depth-independent statistics) the initial failure mechanism is most likely to result in W r between 4 and 8 m. Secondary (i.e. retrogressive) failures further reduce W r to roughly 2-3 m. Secondary failures are more likely for larger horizontal scales of fluctuation, i.e. the peaks for larger horizontal scales of fluctuation are located slightly more to the left in Fig. 8b.
Guidelines (e.g. Blinde et al., 2018;ENW, 2009;MIM, 2016) sometimes positively correlate the dyke width after initial failure (W s1 ) with Fig. 8. Histograms of the residual geometry for several horizontal scales of fluctuation expressed as a) height above the water level (H r -h), b) dyke width at the water level W r and c) indication of excluded and included realisations. Fig. 9. Histogram of the probability of flooding given initial failure, i.e. residual dyke resistance, as a function of the width of the initial failure. residual dyke resistance, i.e. the probability of flooding is reduced or nullified for large W s1 . The results of all simulations with sufficient results (θ h ≥4 m) have been combined in Fig. 9 (again, stable dykes have been excluded from this figure), which shows the distribution of probability of flooding for a given W s1 , i.e. P(F|S 1 , W s1 ). A negative correlation between the remaining width after initial failure and the probability of flooding is present, suggesting that the reduction in the probability of flooding based on W s1 according to the guidelines makes sense. However, even though the positive correlation between W s1 and residual dyke resistance is clear, flooding can still be likely after initial failure even at high W s1 . Therefore, guidelines which specify a 'safe' remaining geometry, i.e. P(F|S 1 , W s1 ) = 0 for W s1 larger than a specified minimum, seem risky, especially when the mechanisms after initial failure are not yet fully understood. However, note that the likelihood of flooding after initial failure is heavily influenced by the adopted strength and softening properties, and by larger horizontal scales of fluctuation (which are more prominent in Fig. 9), as well as by the geometry of the dyke, and is probably much lower in practice. 'Safe' remaining dyke geometries may therefore be more reasonable in practice, although further research is needed to investigate this.

Analysis Set 3: Effect of water level
The effect of different external water levels, i.e. Analysis Set 3, is presented in Figs. 10-12. At first, the probability of initial failure reduces as the pressure on the dyke decreases (0.25 ≤ H i − h ≤ 1.0 m, Fig. 11a). However, once the water level reduces further, outward slope failure, as shown in Fig. 10, becomes more likely and the probability of initial failure increases. Of course, outward slope failure can be prevented by using a more gentle slope. The initiation time of the outward slope   failure in this study seems to be faster in most cases than for inward slope failure, probably due to a smaller failed volume.
The initial width of the dyke W i is larger for lower water levels and therefore, for a larger width, more or larger (retrogressive) failures are necessary to trigger flooding. Moreover, the crest of the dyke needs to settle further to trigger flooding for lower water levels. Consequently, the probability of flooding is comparatively lower for lower water levels ( Fig. 11b and 12). Fragility curves would therefore be useful for the macro-instability failure category, as this would inherently take into account the lower probability of flooding at lower water levels, thereby reducing the over conservatism.
Analysis Set 3 assumes a constant water level for each set of realisations, thereby ignoring that an outward slope failure due to a period of low water, followed by a period of high water, can still be dangerous if the required repair works have not yet been performed. Fig. 13 shows the distributions of the residual geometry for the various water levels; once again, the histograms exclude completely stable dykes and flooded dykes, and have been normalised). The horizontal axes show the residual geometry normalised with respect to the initial geometry. In general, an H i -h of 1.5 m or more offers less initial resistance than an H i -h of 1.0 m or lower due to a switch in the dominant failure mechanism from inner to outer slope failure. However, Fig. 13 shows that a decrease in water level increases the residual resistance, due to the fact that the crest can settle more without triggering flooding (Fig. 13a), and the remaining width at the water level is larger (Fig. 13b). Moreover, Fig. 13b shows that the width of the dyke at the water level may actually increase due to outward slope failure, causing a further increase of residual resistance.

Conclusions
A probabilistic framework for Ultimate Limit State (ULS) dyke failure, i.e. flooding, has been developed and an example analysis presented. The Random Material Point Method (RMPM), a numerical method capable of simulating the full failure mechanism and the effect of spatial variability of the material properties, has been used to investigate the residual resistance of an idealised dyke cross-section. This method accounts for the influence of weaker zones on dyke stability, and it therefore leads to an increased probability of initial failure, and consequently to an increased probability of flooding, compared to standard MPM based on mean strengths.
The residual dyke resistance is relatively small in the example, due in part to the chosen softening parameters. This is especially so for materials with layered spatial variability, where a weak layer which triggers a primary mechanism is also likely to trigger secondary mechanisms, i.e. retrogressive failure. RMPM agrees with existing guidelines which increase residual dyke resistance for smaller initial failures. However, the RMPM analyses do not present a 'safe' remaining geometry, i.e. flooding can always occur due to secondary mechanisms, even if the remaining geometry is large, in part due to the spatial correlation of material properties. The analyses (Analysis Set 3) demonstrate that the water level plays a central role in the dyke assessment, especially for the probability of flooding. Fig. 13. Histograms of the residual geometry for various outer water levels in terms of a) height above the water level (H r -h) and b) dyke width at the water level (W r ).