Exploring the Potential for Multivariate Fragility Representations to Alter Flood Risk Estimates

In flood risk analysis, limitations in the multivariate statistical models adopted to model the hydraulic load have restricted the probability of a defense suffering structural failure to be expressed conditionally on a single hydraulic loading variable. This is an issue at the coastal level where multiple loadings act on defenses with the exact combination of loadings dictating their failure probabilities. Recently, a methodology containing a multivariate statistical model with the flexibility to robustly capture the dependence structure between the individual loadings was used to derive extreme nearshore loading conditions. Its adoption will permit the incorporation of more precise representations of a structure's vulnerability in future analyses. In this article, a fragility representation of a shingle beach, where the failure probability is expressed over a three‐dimensional loading parameter space—water level, wave height, and period—is derived at two localities. Within the approach, a Gaussian copula is used to capture any dependencies between the simplified geometric parameters of a beach's shape. Beach profiles are simulated from the copula and the failure probability, given the hydraulic load, determined by the reformulated Bradbury barrier inertia parameter model. At one site, substantial differences in the annual failure probability distribution are observed between the new and existing approaches. At the other, the beach only becomes vulnerable after a significant reduction of the crest height with its mean annual failure probability close to that presently predicted. It is concluded that further application of multivariate approaches is likely to yield more effective flood risk management.


INTRODUCTION
Over the past two decades, focus, particularly in Europe, has shifted from flood prevention to flood risk management. Flood risk analysis (FRA) is the only quantitative tool available for guiding invest- or continent (Alfieri, Feyen, Dottori, & Bianchi, 2015) have been developed, although each are fundamentally based on a Monte Carlo framework first set out in a flood risk context by USACE (1996) where the economic consequences associated with a large number of plausible events are examined in order to estimate the flood risk, an idea that was conceptualized by the source-pathway-receptor model (Dawson & Hall, 2002;Defra/Environment Agency, 2003) (Fig. 1). The structural performance of the flood defense assets protecting a given area can have a large bearing on these estimates (Apel, Thieken, Merz, & Blöschl, 2004;Eleuterio, Payraudeau, Mosé, & Rozan, 2014;Gouldby, Sayers, Panzeri, & Lanyon, 2010). Most of the key factors that influence a structure's reliability, such as its structural and material properties, exhibit inherent variability in both space and time (aleatory uncertainty) (van Gelder, 2000); however, the majority of the uncertainty in a defense's ability to resist structural failure primarily stems from a lack of data collected on such structural characteristics (epistemic uncertainty). When combined with further (epistemic) uncertainty induced by a lack of knowledge of the exact mechanics behind the failure mechanisms themselves, the rational for adopting a probabilistic description of a flood defense asset's structural vulnerability (that is able to account for all of this uncertainty) in order to achieve accurate flood risk estimates becomes evident. The probability that an asset fails when subjected to a specific loading is referred to as its fragility.
Within an FRA where a large number of plausible loading events are considered, the probability of the defense failing will be required for a range of plausible hydraulic loading scenarios. Until recently, restrictions in the modeling of the dependence structure between the hydraulic loading variables have meant that the loading was considered in terms of a single variable. In order to incorporate the fragility of coastal structures into the analysis, fragility representations have taken the form of a curve, which are referred to as fragility curves. There are a number of methods available for deriving fragility curves. Due to the shortage of data of failures necessary to apply a statistical model and the potential bias introduced by relying solely on expert judgment, processmodel-based curves have become the standard in FRA (Schultz, Gouldby, Simm, & Wibowo, 2010). They are based on limit state equation(s) (LSE) that are established off the back of a large body of engineering knowledge, including field observations, laboratory experiments, and full-scale testing. LSEs describe the relationship between the loadings to which a structure is subjected, collectively denoted by an S, and the resistance, R, to these loadings provided by the structure (Melchers, 1999). They therefore take the following general form: where x 1 , . . ., x m are the basic resistance variables, e.g., material strength and/or crest elevation and x m+1 , .. . ., x n are the basic loading variables, e.g., water level. Once the loadings are greater than the resistance, i.e., once Z < 0, the structure is deemed to have failed. For a given hydraulic loading event, the state of all of the flood defense assets-failed or nonfailed-are obtained by referring to the relevant fragility representations. This determines whether during the event water is able to enter the flooding area through a particular defense (in the case of failure) or simply overtop the defense (when the defense is in a nonfailed state). Once the volume of water entering the flooding area is known, its distribution over the flooding area is estimated by running an inundation model (Lhomme et al., 2008). The assets in the flood plain are grouped into homogeneous classes relating, for instance, to their land use and size. Damage functions that relate the characteristics of the flood, primarily the inundation depth, to the economic damage for a given set(s) of asset classes (Merz, Kreibich, Schwarze, & Thieken, 2010) are then applied to find the economic damage incurred as a result of the hydraulic loading event. By repeating the procedure for a large number of plausible events, an estimate of the flood risk expressed in terms of expected annual damage (EAD) can be obtained.
For coastal structures such as embankments and seawalls, the prevalence of each failure mechanism is dictated by the type and location of breaking. The type of wave breaking can be identified by the Iribarren number (Battjes, 1974) (or surf similarity parameter) ξ , which is expressed as: where tan is the tangent function, α is the bed slope, H s is the significant wave height, and L m is the deep water wavelength. Plunging breakers are present when 0.5 < ξ < 3.3, they arise after the base of the wave begins to feel the sea bed causing it to slow down while the crest continues at its original speed. This action results in the wave steepening and eventually overturning, causing a mass of water to "plunge" down into the water ahead of the wave. By releasing all of their energy during a single localized impact, they are the breaker type most commonly associated with the structural failure of coastal structures such as seawalls and breakwaters (Müller et al., 2008). Although plunging breakers play an important role in shaping the morphology of gravel beaches, it is surging breakers that are most closely associated with the erosion and overwashing of such beaches. They appear once ξ > 3.3, typically on steep slopes where the base of the wave is able to catch up with the wave crest, leading to the wave traveling up the beach face as a wall of water. For ξ < 0.5, spilling breakers emerge when a wave becomes unstable after its base comes into contact with the bed. This contact causes the wave to steepen and a turbulent mixture of air and water to "spill" down the front of its crest. In addition to the breaking conditions, the exact destructive potential of a wave will depend on the location of breaking. In the nearshore, breaking is primarily depth induced. Given the beach slope, the type and location of breaking is therefore almost exclusively controlled by water level and the wave climate (a combination of wave height and wave period). If the fragility is to be accurately represented it is thus essential that the probability of failure is expressed explicitly on each of these variables. Gouldby, Méndez, Guanche, Rueda, and Mínguez (2014) proposed a methodology for deriving extreme nearshore hydraulic loading conditions. The first step involves applying the Heffernan and Tawn model (Heffernan & Tawn, 2004) (herein referred to as HT04) to several of the offshore hydraulic loadings that characterize an extreme sea state. The HT04 model is composed of a series of nonlinear regression models. Each captures the shape of the dependence between the variables when a particular variable is deemed to be extreme after the data are standardized (by invoking the probability integral transform) to possess common marginal distributions. In contrast with the existing copulabased models, the HT04 model does not prescribe a specific dependence structure to the joint tail regions, thus it removes the assumption of asymptotic dependence or independence between the variables required by earlier models. A Monte Carlo simulation procedure outlined in more detail in Lamb et al. (2010) and Wyncoll and Gouldby (2013), among others, is used to generate samples from the fitted model. It commences by conditioning a specified variable to be extreme. A joint residual from those obtained during the fitted procedure when the specified variable was conditioned to be extreme is then sampled. The values of the remaining variables at the transformed scale can then be obtained using this residual and the relevant regression model. The sample is subsequently transformed back to the original scale by applying the inverse probability integral transform to each realization. The sample is, however, rejected if the conditioning variable is not the most extreme on the transformed scale. The procedure is repeated, conditioning on each variable in turn, until the proportion of samples where each variable is most extreme, on the transformed scale, is consistent with the empirical joint distribution.
In the method put forward in Gouldby et al. (2014) the HT04 model is applied to model H s , wind speed (U), and surge (S). The wave steepness st which is dependent on and physically restricted by H s , is also a function of the wave period (T m ). A regression model of steepness given H s was therefore developed to generate the T m , which in deep water is related to the wavelength through the dispersion relation: The wind direction (θ U ) and wave direction (θ H s ) are sampled independently from their empirical distributions, conditional on H s and U. The total water level h is composed of the astronomical tidal level and S. Given the magnitude of the simulated S, to find the astronomical tidal level, the month in which the surge is assumed to have occurred is first determined by sampling from the conditional distribution of the month given the surge. This conditional distribution is estimated from the empirical surge data. A year from the 18.6-year tidal cycle is also independently sampled. The astronomical tidal level is then obtained by randomly sampling an astronomical tidal level from those observed in the sampled month and year. Numerical wave models, such as SWAN (Booij, Ris, & Holthuijsen, 1999), are commonly adopted to propagate wave conditions from the offshore to the nearshore where waves are subjected to several transformational processes such as refraction and shoaling. FRA generally requires the simulation of a large number of plausible sets of loadings. The computational intensity of numerical wave models can make their adoption in such applications impractical. In Gouldby et al. (2014), a maximum dissimilarity algorithm is therefore applied to the set of simulated offshore conditions in order to obtain a subset of conditions with sufficient coverage of the parameter space. The subset of conditions is propagated to the nearshore using SWAN, while a meta-model is applied to the input and output of the SWAN modeling and used to propagate the remaining conditions to the nearshore, thus ensuring the methodology is computational tractable for implementation within an FRA. The Gouldby et al. (2014) methodology has been used to generate extreme nearshore loading conditions at a 1-km resolution around the U.K. coast (Gouldby et al., 2017).
It is thus now possible to generate a long synthetic record, which captures the dependence structure between the individual loadings, of the offshore conditions that characterize an extreme sea state in the source component of an FRA. Consequently, multivariate fragility representations that express the probability of the defense failing simultaneously conditional on each of the loadings can now form the pathway component of a system-wide FRA. For a more thorough discussion of the role of fragility in FRA and the rationale for multivariate fragility representations, see Jane et al. (2016). The adoption of multivariate fragility representations has already been shown to result in substantial differences in the estimated structural vulnerability of seawalls compared with when the representations presently adopted in nationwide U.K. FRA are assumed (Jane et al., 2015(Jane et al., , 2018. Gravel beaches are composed of sediment with a medium-size diameter, D 50 , between 2 mm and 256 mm. They are commonplace along the world's high-latitude coastlines (e.g., the United States, Japan, and northern Europe). Constituting approximately a third of the coastline of the United Kingdom (Fuller & Randall, 1988), they are highly effective at dissipating wave energy either alone or as part of a broader flood defense scheme, which often includes structures such as seawalls or offshore breakwaters. Consequently, they form a vital part of the U.K. flood defense infrastructure. Laboratory studies have shown that, similarly to embankments and seawalls, their structural vulnerability depends on the combination of wave height, period, and water level (Obhrai, Powell, & Bradbury, 2008); thus a multivariate representation is required if their fragility is to be accurately represented. As they are so widespread, any changes in the assessment of a structure's vulnerability by adopting new representations has the potential to not just alter regional flood risk estimates but also impact on national flood estimates.
The aim of the article is to examine the changes in the estimated structural vulnerability of a gravel beach when the potentially more robust fragility representations, derived according to the methodology proposed in this article are adopted, compared to when the representations used in present-day national U.K. FRA are assumed. First, the measures available for estimating the response of the gravel beach to an extreme event are evaluated for the purpose of application in FRA (Section 2). The case study site containing the localities where the fragility representations will be applied is subsequently introduced (Section 3). The existing representation of fragility for a gravel beach used in U.K. FRA are then shown (Section 4) before the method for deriving a new representation that explicitly considers each of the variables, site-specific information as well as a more robust method of predicting the structural failure of a beach, is explained (Section 5). Once these new representations are derived for a case study site in the United Kingdom, the changes in the vulnerability of the gravel beach are analyzed (Section 6) before appropriate conclusions are drawn (Section 7).

REVIEW OF GRAVEL BEACH FAILURE PREDICTORS
A comprehensive overview of the classification, modeling, and management of barrier beaches is given by Stripling, Bradbury, Cope, and Brampton (2008). However, at present, information on profiles during storm events is sparse. Poate, McCall, Masselink, Russell, and Davidson (2012) undertook an analysis of pre-and post-storm profiles of several beaches that lie along the south coast of the United Kingdom. They found that their responses during the storm period were seemingly uncorrelated with the nearshore conditions, thus indicating that if there is a relationship, then it is a complex interaction between the initial profile of the beach and the nearshore conditions. The lack of field data has led to empirical models based predominantly on laboratory experiments, occasionally incorporating limited field data and, more recently, process-based numerical models to predict the short-term response of shingle beaches to storm events.
Given a sufficient supply of sediment, over time the ridge will adapt to the stationary water level and wave conditions to approach an equilibrium profile. On the basis of tests at the 1:17 scale, where coal was used to represent shingle of four material sizes and gradings, Powell (1990) derived a model to describe the profile. In each test, waves from a JONSWAP wave spectrum (Hasselmann et al., 1973) were grouped in terms of steepness. The wave conditions in each grouping were applied to the beach in order of increasing severity for a burst of 3,000 waves, which was found to be a sufficient number of waves to achieve a near stable state. The final profile of the beach was defined by a set of three curves primarily determined by the sediment size and loading conditions to which the beach was subjected, through a set of functional relationships. The effect of a permeable core was also investigated. This was found to reduce the volume of sediment able to adapt to the loading conditions, which in turn increases the mobility of the active sediment. It also reduces a beach's ability to dissipate wave energy through infiltration. Both of these factors increase a beach's susceptibility to failure. When assessing the effect of an impermeable core under identical conditions to those adopted in the other tests (but only applying bursts of 1,000 waves on each occasion), it was discovered that the compatibility of the core has a significant influence on the horizontal displacement of the profile. On the basis of these results, Powell (1990) provided a correction factor to account for this effect when predicting the beach profile.
The Powell (1990) model has previously been adopted in reliability analysis (Buijs, Simm, Wallis, & Sayers, 2007) where failure is considered to have occurred once the distance that the crest retreats by is greater than the original crest width. However, there are concerns regarding its application in reliability analysis due to the assumption that a constant supply of sediment is available, which is likely to be unrealistic, particularly in storm conditions. Furthermore, the short-lived nature of such events means that there will not be sufficient time for the profile to reach equilibrium. In addition, as the maximum sediment size used in the laboratory experiments corresponds to D 50 of 30 mm at full scale, the model can only be considered for shingle beaches composed of relatively fine sediment.
Overwashing is the primary process behind breaching of a beach. It is therefore often employed as a conservative surrogate for predicting breach. An empirical model for predicting whether overwashing and thus a breach will ensue for a given set of hydraulic loadings from laboratory tests and fieldwork based on Hurst Spit (Hampshire, UK) is given by Bradbury (2000). During the modeling, it was deduced that as well as the height of the barrier above the still water level influencing the chance of overwash occurring, for a given freeboard, a barrier with a large cross-sectional area is less likely to fail than one with a small cross-sectional area. This led to the introduction of the barrier inertia parameter B i , a function of barrier freeboard R c and area above the still water level denoted B A (m 2 ) nondimensionalized by the significant wave height, H s , By plotting the results of the tests by a regression of the inertia parameter on the wave steepness H s L m , where L m is the deep water wavelength, it was discovered that overwash can be expected for  (2000) Obhrai, Powell, and Bradbury (2008).
Although the results can be used to assess whether a barrier will breach during a given shortterm loading event, the results from the tests only covered a limited range of wave steepnesses, 0.015 < H s L m < 0.032. In addition, Bradbury, Cope, and Prouty (2005) found that the model performed poorly when applied to sites other than Hurst Spit.
A series of tests were undertaken at HR Wallingford to examine the threshold for overwashing for both swell and storm conditions in order to extend the valid prediction range of the Bradbury model (Obhrai, Powell, & Bradbury, 2008). The sediment sizes ranged from a D 50 of 16-57 mm. Once an equilibrium profile was established, the wave height was increased incrementally for sets of 1,000 waves until the barrier failed. It was concluded that the existing formula overpredicts the threshold for breaching under swell wave conditions and underpredicts the threshold for steep waves. The original Bradbury model was therefore reformulated incorporating the new data to provide a more accurate model for the overwashing threshold of a shingle barrier as follows: valid for 0.01 < H s L m < 0.06. Uncertainty remains in the exact location of the threshold. To account for this, factors are applied, upon an inspection of the Bradbury (2000) and Obhrai, Powell, and Bradbury (2008) data, to both the strength and loading terms in the model. The resulting likelihood of ex-  ceeding the overwashing threshold given the barrier inertia parameter and wave steepeness are shown in Fig. 2 alongside Equation (6) and the data sets used in its derivation. The investigators also found that the level of the hinterland influenced the likelihood of overwash since it determines whether sediment remains part of the structure once it is over the crest or whether it is deposited further inshore contributing to a reduction in both height and width of the structure, thus increasing its vulnerability. It was concluded that any future model should incorporate these effects as derived on the basis of laboratory tests and as with all physical model tests may be subject to scale effects. Originally developed to predict the response of sandy beaches to hurricanes, X-beach (Roelvink et al., 2009) is a process-based model for computing nearshore hydrodynamics and morphological responses of the nearshore, beaches, and dune systems during storm events. It has been adapted for gravel beaches where it is referred to as X-Beach-G   . It includes an extension to the standard model to account for the influence of short waves in intermediate and shallow water depths on the beach profile. These effects are much more pronounced for steep reflective gravel beaches than at their sandy counterparts. It also includes the addition of a groundwater model to capture the infiltration and exfiltration that occur on such beaches. Although overtopping and overwash cannot be successfully simulated without morphodynamic updating of the bed level, which is not presently possible in the modified X-beach model, the type of response can be inferred from the simulated hydrodynamics on the initial barrier profile. Mccall et al. (2013) used the model to predict whether overwash occurs for 22 historical storm impacts on five gravel beaches in the United Kingdom and three BARDEX physical model experiments. In these cases, the model was shown to offer more accurate predictions than the empirical models albeit at a greater computational cost and with greater uncertainty in overwash predictions near the empirical threshold for overwash. In addition, the depth of the gravel beach toe and the gravel beach slope were found to have a large influence on the threshold criteria for overwash; these are unaccounted for in the empirical models. Although these early results are promising, further validation of the model will be required to establish whether the improvement in the accuracy of the predictions, compared with those given by the empirical model, warrant the extra computational effort of its adoption. The Obhrai (2008) reformulation of the Bradbury (2000) model will therefore be adopted to analyze the fragility of gravel beaches in this article.

CASE STUDY SITE
Chesil Beach is a coarse (D 50 = 50 mm at Chiswell (Carr, 1969)) pure gravel beach situated on the south coast of the United Kingdom (Fig. 3). The point of inflection of the rear slope on the average profile corresponds to a maximum on the gradient plot. The elevation of the average profile at this chainage (6.561 mODN) is the elevation at which the base of the rear and front slopes will be defined when fitting "idealized" profiles to those recorded at the site.
Although debate surrounds its origin, it is commonly believed to have formed as a result of the landward migration of an offshore bar, composed of excavated sediment and fluvial-glacial deposits, as sea levels rose during the holocene marine transgression (Bray, 1990;Carr & Blackley, 1973). Subsequently, the volume of the coarse sediment that is believed to lie on top of a more impermeable core (May, 2003) increased due to the longshore transport of sediment eroded from cliffs to the west of the structure. However, due to human intervention reducing its sediment supply both directly, as well as through contributing to the emergence of natural headlands such as Golden Cap, the beach can now be considered a closed system (Bray, 1997). This has culminated in reduction in its total volume and as a consequence left the beach more vulnerable during storm events. A beach of scientific significance as well as recreational and esthetic value, in common with many shingle beaches situated along the coastlines of northern Europe, it plays an essential role in defending the adjacent land from flooding. Just as with man-made structures, such natural defenses when exposed to prolonged extreme swell conditions or storm events can fail. For shingle beaches, failure is typically considered to have occurred once the structure is breached. Breaching is generally defined as the short-term lowering of crest height as a result of wave-induced overwashing of the beach (Bradbury, 2000). Flooding, to a lesser extent, can also occur by water overtopping and/or percolating through the structure.
The wave climate of the south coast of the United Kingdom is dominated by Atlantic swell and locally generated wind waves. Bimodal sea states, where long-period swell waves coincide with wind waves, are not uncommon and have been known to exacerbate the impacts of an event (Mason, Bradbury, Poate, & Newman, 2009). On the basis of a series of laboratory tests, Polidoro, Pullen, Eade, Blanco, and Mason (2017) recently derived an empirical model for predicting the morphological response of a gravel beach to bimodal sea states. However, the authors stressed that the model is not designed as a means of predicting breach and it was therefore not considered in the earlier review of gravel beach failure predictors.
Hurricane strength winds during the "Great Storm" on November 22-23, 1824, heralded the arrival of one or two tsunami-like waves and an accompanying storm surge, which led to the destruction of many houses and the deaths of 50-60 people (Haslett & Bryant, 2009). Although by far the most severe event on record at Chesil, the eastern end of Chesil Beach, where this study focuses, has a long history of flooding (West, 2017a). After flooding events in 1942 and 1956, construction commenced in 1958 of a vertical concrete seawall to protect the most exposed section of the village's frontage. The wall was eventually completed in 1965. Further flood defense measures were deemed necessary after the village experienced a succession of flooding in the late 1970s (Williams & Hardwick, 1996). In 1981, a 150-m gabion crest protection mattress was installed immediately to the north of the wall to strengthen the beach, with improvements to the wall including the addition of a wave return wall also undertaken. The Flood Alleviation Scheme implemented in 1986 saw a 550-m long drainage channel constructed behind the ridge to collect any water that seeped through it, as well as the drilling of an 18-m deep steel sheet pile into the back of the ridge in order to prevent water flowing through the sand and shingle at the foot of the beach (West, 2017a). On February 5, 2014, substantial damage was incurred by the beach, including the destruction of a number of crest protection mattresses around 100 m to the east of the sites studied in this article. This was the latest of a series of storm events during the winter of 2013-2014 that ultimately resulted in around 150, 000 m 3 of sediment being removed from the beach (Environment Agency, 2015). By February 3, 2014, an Atlantic depression had reached the southwest tip of Ireland; the resultant storm surge coincided with south-southeasterly wind waves generated by a storm in the English Channel ahead of the depression (Sibley, Cox, & Titley, 2015). This culminated in the failure of the seawall at Dawlish on the opposite side of Lyme Bay the following day. The long-period swell generated by the Atlantic depression arrived on the morning of February 5. The coupling of the high potential for offshore transportation of sediment at Chesil owing to its southwesterly orientation resulting in storm wave generally approaching parallel to the beach and the arrival of the swell coinciding with high tide, gave rise to Fig. 10. Scatter plots of the (normalized) variables that define the geometry of the beach (lower panels) are given along with their accompanying Kendall's τ coefficient and its associated p-value (upper panels). a highly destructive storm event (Masselink et al., 2016). It resulted in structural damage to the gravel beach as well as to the toe of the seawall protecting the village of Chesilton. This once again shone a spotlight on the vulnerability of the beach and the important role it plays in defending the adjacent land from flooding.

EXISTING FRAGILITY CURVES
When implementing FRA at the national scale in the United Kingdom, due to the expense of obtaining detailed site and structure-specific information for every defense, each flood defense asset is represented by one of 61 generic structures. The generic structures can broadly be classified into four fluvial defenses (Vertical Walls/Slope or Embankment/ High Ground/Culvert) and three coastal defenses (Vertical Seawall/Sloping Seawall or Dike/Beach) with each class further subdivided according to structural characteristics such as width or material composition (Environment Agency, 2004). Each defense is assigned a condition grade that ranges from 1 -Ex-cellent to 5 -Very poor on the basis of a visual inspection of the structure and linguistic descriptions of the condition grades (Environment Agency, 2012). A set of five generic fragility curves therefore exist for each asset, one for each condition grade.
To derive a particular curve, any structural characteristics that influence a defense's propensity to failure are either modeled by the relevant distribution or assigned an appropriate fixed value. It is the parameter of these distributions and constants that is adjusted to delineate defenses belonging to the different condition grades. A large number of samples comprising a single realization from each distribution are then simulated. For each sample, a set of loadings that produce a particular overtopping discharge are also generated and the LSE subsequently evaluated. The conditional probability of failure is thus simply the proportion of samples for which the LSE is less than zero. The process is repeated for sets of conditions corresponding to different overtopping discharges to produce the curve.
In U.K. FRA, a gravel beach is assumed to have failed once the distance of crest retreat is greater than the initial crest width of the beach, i.e., once Z < 0 in the following LSE: where w is the width of the beach and p c is the distance of crest retreat according to Powell (1990). The values assumed for the structural variables that were selected to characterize a shingle beach in each condition grade are given in Table I with the corresponding fragility curves shown in Fig. 4. Although none of the structural variables are modeled probabilistically, due to differing extents of crest retreat manifesting from combinations of the water level, wave height, and wave period that produce identical overtopping discharges, the fragility curves appear probabilistic.

DERIVING MORE ROBUST FRAGILITY REPRESENTATIONS
Profiles of the beach for three locations labeled as 6a00131, 6a00140, and 6a00144, shown in Fig. 5, are available at approximately six-month intervals between 2007 and 2013 (Fig. 6). It is evident that at each location along the beach the rear face of the shingle ridge remains almost wholly unaltered throughout this period. The beach experiences both erosion and accretion, causing the front face to migrate seaward and retreat, predominantly altering the position and width of the crest as well as the front slope and, ultimately, the vulnerability of the ridge in the process. In the Obhrai, Powell, and Bradbury (2008) model, the shape of the profile is characterized by a small number of structural parameters relating to an "idealized" profile of a shingle ridge (Fig. 7). They are the crest width c w , front slope N f , and rear slope N r , as well as the crest height c h through the consideration of the relative freeboard R c . To incorporate the observed profiles into the fragility calculations, a procedure for deriving their corresponding "idealized" profiles, which are completely defined by three straight lines representing the rear slope, front slope, and the crest, respectively, was developed. To initiate the process, an estimate of the average elevation of the base of the ridge is required. The base is defined to lie at the point of inflection of the rear slope to ensure that a conservative estimate of the cross-sectional area of the backface of the ridge is obtained. The chainage of this point corresponds to a maximum when the gradient between the beach's average profile and point with its maximum elevation is plotted against the chainage (Fig. 8). Once the elevation has been obtained, when fitting an "idealized" cross-section to a profile recorded at the site the chainages located landward and seaward with the same elevation are taken as the location of the base of the rear and front slopes, respectively. The intersection of the lines that define these two slopes and the horizontal line defining the crest that minimizes the root-mean-squared error between the lines and the observed profile is defined as the profile's corresponding "idealized" profile. Once fitted, the structural parameters required for evaluating the LSE follow through trivial calculation. Examples of the "idealized" cross-sections fitted to the observed profiles of a gravel beach are shown in Fig. 9. Subsequently, to capture the uncertainty in the ridge profile at a point in time, the parameters are modeled by relevant statistical distributions as is standard in FRA (Simm et al., 2008).
Scatter plots of the structural variables that define the geometry of the shingle ridge extracted from the observed profiles taken at the locations along the beach are shown in Fig. 10. As may have been anticipated, the crest width was found to be negatively correlated with both the elevation of the crest and the gradient of the front slope, which consequently were themselves positively correlated. At each of the sites, the rear slope remained relatively fixed and therefore, with the exception of the crest height, exhibited no discernible associations with the other structural variables. In order to quantitatively assess the degree of (linear) correlation between the pairs of variables, Kendall's τ correlation coefficient was computed. The τ coefficient is a parametric measure whose value is equivalent to the difference between the probability of concordance and probability of discordance for a pair of variables. Consequently, it lies between [−1, 1] where 1 corresponds to perfect positive linear correlation between the variables and −1 corresponds to a perfect negative linear correlation. The Kendall's τ coefficient for each pair of the variables is given in the upper panels of Fig. 10 along with their associated p-value (≤0.01) when testing the null hypothesis τ = 0, i.e., that there is no correlation against the alternative hypotheses τ = 0-that some correlation is present. All of the p-values, with the exception of those associated with the correlation between the rear slope and both the front slope and crest width, are negligible and so the null hypothesis of independence between the pairs of variables is rejected. Thus the possibility of the structural variables being dependent cannot be ruled out. So, when modeling the structural input parameters this dependence between these pairs of variables will have to be accounted for if an accurate description of the fragility is to be achieved for the shingle ridge at Chesil Beach.

Copula
The d-dimensional copula is a multivariate distribution function on I ∈ [0, 1] d with uniform [0,1] marginals. Sklar's theorem (Sklar, 1959) provides the foundation for their application, it states that if F is a  Evolution of the fragility surface at site 6a00140 as the crest height is progressively lowered. In each case, the mean water level is set at 1.345 mODN, the mean water level of the extreme events (events with a return period >1 in 12.5 years). decomposed into a copula and a set of marginal distributions. Consequently, it enables the modeling of the dependence structure to be carried out independently from that of the marginal distributions. This offers a more flexible approach than traditional multivariate modeling, where the choice of joint distribution is restricted by the choice of possible marginal distributions. For a more in-depth introduction to copulas, consult the classical texts by Joe (1996) and Nelsen (2006).
Since the copula only allows uniform margins the observations are transformed to [0,1] applying the nonparametric transformation, given below, to obtain the (normalized) ranks (Genest & Favre, 2007), where X i, j is the ith observation out of n of the jth variable, respectively. Now that only information relating to the dependence structure between the variables remains, the copula can be fitted. The limited number of available profiles has the potential to introduce bias into the copula parameter estimates and, ultimately, the overall assessment of the beach's vulnerability at these sites. A bootstrapping procedure will therefore undertake to assess the bias and uncertainty that is a product of the limited number of observed profiles. The Archimedean and elliptical are the most well-known copula families. Archimedean copulas include the Gumbel, Frank, and Clayton copulas. They use a single parameter to define the dependence structure between a set of variables. This can inhibit their ability to accurately capture the dependence structure between a set of variables, particularly when considering a large number of variables where the dependence structure often becomes complex. More flexibility can be achieved by arranging multiple bivarite copulas in a tree structure when constructing higher dimensional copulas, e.g., hierarchical/nested Archimedean copulas (Joe, 1997;Mc-Neil, 2008) or vine copulas (Joe, 1996). However, due to the desire to keep the proposed methodology as accessible as possible they are not considered here. The elliptical class of copula comprise the Gaussian and Student's t copula; both are radially symmetric. The Gaussian copula is the most commonly applied in practice. In the d-dimensional setting, it is given by: where R denotes the joint distribution function of the d-variate standard normal distribution function with linear correlation matrix R, and −1 denotes the inverse of the distribution function of the univariate standard normal distribution. Just as the Student's t distribution is a generalization of the standard Gaussian distribution, the Student's t copula is a generalization of the Gaussian copula. In an analogous way to the univariate case, the Student's t copula assigns more probability density to the tails of the copula than the Gaussian copula, in particular to the upper and lower tails, thus increasing the probability of occurrences of joint extreme events. The Student's t copula contains an additional parameter compared to the Gaussian copula, ν, the degrees of freedom, which controls the magnitude of the density shifted to the tails, which increases as ν decreases, with the Student's t copula approaching the Gaussian copula as ν → ∞. Both copulas were fitted to the transformed crest heights, crest widths, and rear and front slopes observed at the sites. The Akaike information criterion (AIC) (Burnham & Anderson, 2002), a goodness-offit criterion that provides an assessment of the quality of model fit while penalizing model complexity was adopted to choose between the competing copulas. The Gaussian copula was found to achieve the better AIC of the two copulas and so was selected to model the dependence between the variables.

Fragility Surfaces
A Monte Carlo approach similar to that put forward in Simm et al. (2008) is adopted to derive  the site-specific multivariate fragility representations. The method proposed in this article commences in an identical manner with the identification of the failure mechanism(s) that a structure is susceptible to, and the selection of the most appropriate means of modeling them, as determined for a gravel beach in Section 2. To incorporate information relating to the specific structure and to express the failure conditionally on each of the hydraulic loadings, subsequent steps of the approach in Simm et al. (2008) require modification. According to the LSE nominated to model the structural vulnerability of a beach, the hydraulic loading conditions required to estimate its likelihood of breaching are H s , T m , and h, the latter measured relative to the same datum as beach's crest elevation. A brief outline of the sequence of steps required to obtain the probability of breaching given a set of these hydraulic loadings once the copula is fitted to capture the dependence between the structural variables at a particular site is given below. Fig. 18. Scatter plots of the (normalized) variables that define the geometry of the ridge at the East Beach site (lower panels) are given along with their accompanying Kendall's τ coefficient and its associated p-value (upper panels).

STEP 1: Simulate a large number of profiles.
Simulate uniform random variables with the observed dependence structure by sampling from the copula. A smooth estimate of the probability density function of each variable at a given site was provided by a Gaussian kernel density estimate (Parzen, 1962;Rosenblatt, 1956) with a bandwidth selected according to the rule of thumb given in Silverman (1986). By combining the relevant marginal distribution of each of the variables with the inverse probability integral transform, these uniform random variates can then be transformed to give the structural variables (i.e., c h , c w , N f , and N r ) that fully describe a particular cross-section of the gravel beach at Chesil (see Figs. 11-13). STEP 2: Account for uncertainty in the LSE. Generate a random variate from the distributions associated with the model uncertainty in the strength (γ R ) and loading (γ S ) elements of the LSE, for each simulated profile.
STEP 3: Evaluate LSE for each profile. The LSE, Z, associated with the breach of a gravel beach according to the Obhrai, Powell, and Bradbury (2008) reformulation of the Bradbury (2000) model, while also accounting for model uncertainty is: where B A = R c (c w + 1 2 R c [N f + N r ]) and R c = c h − h. STEP 4: Estimate failure probability. The estimated failure probability is then given by the proportion of simulated profiles for which Z < 0, i.e., where and Z i is Z corresponding to profile i ∈ {1, . . ., N}, where N is the total number of simulated profiles.
To generate a multidimensional fragility representation, these steps are repeated for every plausible set of (H s , T m , h). A summary of the structural parameters that describe the profile of the gravel beach at each of the sites is provided in Table II. In comparison with the corresponding values assigned when deriving the generic curves, it is quickly apparent that the gravel beach at the selected sites is steeper and wider than was previously assumed. Given the latter result in particular and the observations in the laboratory by Powell (1990) and Obhrai, Powell, and Bradbury (2008), it is reasonable to expect that the gravel beach may be more resilient to extreme storm events than would have previously been assumed. Multidimensional fragility representations for the site at 6a00140 are given in Fig. 14.

COMPARISON OF FRAGILITY REPRESENTATIONS
The annual failure probability (AFP) of a defense is the probability that it suffers one or more structural failures within a single year. In order to quantitatively compare the change in the description of the fragility of the gravel beach between the newly proposed and existing fragility representations, their associated AFPs are then computed. The AFPs are estimated by integrating the fragility representations across the distribution of the hydraulic loadings. This is implemented in practice by considering the structural response of a defense to hydraulic loading events with a range of return periods within a Monte Carlo framework. In this article, a long-term synthetic record of the nearshore loading conditions, spanning a full range of return periods, is obtained according to the approach given in Gouldby et al. (2014). The probability of failure for a set of simulated nearshore loadings, i.e., p f |x where x = [H s , T m , h], can be obtained directly from the new representations where the probability of failure is expressed conditionally on each of these basic loadings. The univariate fragility representations, however, require overtopping rates to determine the probability of failure for a corresponding set of nearshore conditions. The nearshore conditions were therefore propagated to the toe of the structure and the overtopping rates calculated, i.e., x = q according to the formulas given in Pullen et al. (2007). For a given simulated event x i , the defense state, D, can be considered either failed D i = 1 or nonfailed D i = 0, according to D i ∼ Ber noulli ( p f |x i ). After generating a defense state for each of the simulated events, N, under the strong assumption of independence between the n events arising in a given year, the AFP can subsequently be computed as follows: By repeating the process of determining the defense state for each simulated event and subsequently computing the AFP, the distribution of the AFP of a defense associated with a particular fragility representation can be built up. All of the fragility representations concurred that the present-day AFP for each of the sites was zero. This may appear surprising since there are many flood events within living memory; however, these, along with previously documented events, have almost exclusively been attributed to overtopping and/or percolation of water through the beach rather than due to breaching. There is evidence of washover fans that arise as sediment is transported and subsequently deposited further inland, often in a shape resembling that of a fan, which were possibly created here during the "Great Storm" of 1824 (West, 2017a). Although there is no conclusive evidence of a breach in the recent geological past, the likely occurrence of overwashing during this period, when the beach is likely to have been in a form similar to that which it currently takes, means that further investigation into the likelihood of breach as its "strength" is reduced is warranted.
The mean AFP was recalculated as the crest height of each defense was progressively lowered in order to explore how the vulnerability of the defense changes as its ability to resist structural failure is reduced. By repeating the procedure for each fragility representation, a greater understanding of the magnitude of the variation in the estimated vulnerability of the defense between the different representations can simultaneously be achieved. The results are shown in Fig. 15 with the crest heights for which significant mean AFPs are first exceeded provided in Table III. Both illustrate that according to the newly derived fragility representations, site 6a00144 is typically the least vulnerable, as expected since it is the widest section of the beach. Its fragility is approximately equivalent to that of the generic beach when in an "Excellent" condition. The beach at site 6a00140 is straddled between the estimated vulnerability of the generic beach in various condition grades. For lower crest elevations (<7.2 mODN), site 6a00140 is the least vulnerable of the three sites, whereas for higher crest elevations, its vulnerability is very similar to that of the beach at site 6a00131. The structural vulnerability of the beach at site 6a00131 is estimated to lie between those of the other two sites at lower crest heights, resembling that of the generic beach in condition grades 1 and 2.
According to the new fragility representations, the crest height at site 6a00144 can be allowed to drop lower than the other sites before the AFP become equivalent. This suggests that site 6a00140 and, in particular, site 6a00131 should be prioritized over site 6a00144 if any intervention work were to be carried out on the beach. This highlights a practical application of these fragility representations as part of an approach to managing the maintenance of coastal structures. By relating the AFP to the condition of the structure, in this case solely described by the crest height, it provides a quantitative basis upon which decision making relating to the adoption of relevant intervention options can be made. Although the crest height is a key variable in determining the fragility of a beach for a robust analysis, each of the structural characteristics that determine the profile of the beach would have to be considered. For instance, the mean AFP could be displayed conditionally on each of the structural characteristics in a form similar to that of the surfaces in Fig. 14. In addition, once the structure deteriorates to such an extent that a significant AFP threshold is exceeded, such as an AFP of 1%, then this can act as a trigger for any intervention to be carried out. At present, the timing of interventions is often linked to the condition of the structure (Wallis, Whitehouse, & Lyness, 2009); however, relating them directly to the structural vulnerability of a defense can be expected to further increase their effectiveness (Ogunyoye, 2011).
The proposed methodology was subsequently applied to East Beach, a gravel beach located immediately eastward of the harbor at West Bay, at the eastern end of Chesil Beach (Fig. 16). A site of special scientific interest (SSSI), the beach affords protection to approximately 254 residential and 35 commercial properties. In response to numerous flood events over the preceding years, including a particularly severe event in 1974, a shingle bank was built up to a height of 7.5 mODN in 1986 (Dorset County Council, 2016). In order to maintain its capabilities as a flood defense, the beach is periodically reprofiled and if any significant reductions in the volume of sediment on the beach are observed, beach recycling is also carried out (Dornbusch & Cargo, 2011). Despite this, the beach is narrower with a lower crest elevation as well as shallower rear and front slopes than the sites at the opposite end of Chesil Beach included in this study (see Fig. 17 and Table IV). The values of the structural variables synthesized from beach profiles indicated that the crest height was not significantly correlated with the crest width or either of the rear or front slopes (Fig. 18). A significant correlation was found between the crest width and both slopes; therefore, a trivariate copula was fitted to model the dependence between these variables. For a set of structural variables simulated from copula and transformed back to their original scale, realizations of the crest height are generated (independently) from the kernel density estimate of its probability density function to complete the description of the shape of the beach. Sets of these simulated variables that define the shape of the profile along with the observed values are plotted in Fig. 19.
The mean AFP of the site at East Beach was found to decrease from 5.301 × 10 −3 to 1.478 × 10 −3 when the new representations were adopted (see Fig. 20). Since the original figure lay outside the 95% confidence interval of estimate of the AFP given by the newly proposed approach, it can be considered a significant decrease. This is in contrast to the results observed at the sites at the other end of Chesil Beach where the initial results suggested that the estimates of a beach's structural vulnerability given by the existing representations may closely emulate those produced by their new multivariate site-specific counterparts. The magnitude of the disparities in the AFPs will be dependent on the condition grade assigned to the beach when implementing the existing approach and its geometry, in combination with the nearshore conditions to which it is exposed. Thus, there are unlikely to be any universally applicable rules for predicting the size of the disparity at a particular beach.
To ascertain the proportion of the uncertainty in the estimate of the beach's AFP that can be at-tributed to the small number of profiles recorded along transect 6a00688, a large number of bootstrap resamples of these original profiles were drawn. The process of fitting a copula to model the dependence between the sampled profile's structural variables through to estimating the AFP was repeated for each sample. The estimates of the fitted Gaussian copula's parameters for each resample are shown in Fig. 21. Although the bias introduced by the small number of profiles was negligible, it was responsible for around 45% of the overall variation in the AFP. Hence, increasing the frequency at which profiles are recorded along a transect is certainly recommended as a means of providing potentially substantial reductions in the uncertainty associated with the AFP.
When these new fragility representations are subsumed into the existing FRA framework, potentially more accurate flood estimates will automatically be produced for applications where simple identification of the risk is required such as when conveying flood risk to the public or for the setting of premiums by insurance companies. These estimates are also used by flood risk managers to guide investment planning and decision making, where the aim is often to achieve some predetermined reduction in risk. When implemented at the national level, the potentially more accurate quantification of the risk can be expected to instigate more effective allocation of resources. At the more regional scale, it may bring about an improvement in the selection of the mitigation measures for the assets that defend a given area. At the same scale, it can also precipitate more suitable defense schemes in future. Given the wide array of geometries of shingle beaches around the U.K. coast at other locations, for instance, at a beach of a similar height to that at Chesil but with a narrower width, the existing fragility representation may underestimate its propensity to breach, thus necessitating (further) strengthening of the beach or the implementation of other flood defense measures. On the contrary, at sites similar to East Beach the existing fragility representation may have provided a conservative estimate of its vulnerability. In other words, the beaches may potentially provide a more significant role in the defense of a stretch of coast than previously assumed, making other structures obsolete (which may have been considered when designing a new coastal defense system) or at least necessary at a reduced scale, and thus offering potential savings for the agency responsible for the management of the coastline.

CONCLUSION
Assessing the fragility of coastal structures has become an established procedure when conducting FRA. The results of the assessment have been shown in the past to heavily influence flood risk estimates. It is therefore essential that defense fragility is accurately assessed. To date, in the majority of FRAs, the probability of a defense failing is considered conditionally on a single loading variable. At coastal locations in particular where multiple loadings influence a defense's performance, simple intuition suggests that a better description of a defense's fragility will be obtained if each variable is considered explicitly rather than if they are amalgamated into a single composite loading as is currently the adopted approach. To gauge the impact on the estimated structural vulnerability of a gravel beach by considering each loading explicitly, incorporating site-specific information and an improved method of predicting breach, new fragility representations were derived for several sites at the eastern end of Chesil Beach. Both the existing and new fragility representations concurred that in its current form the AFP of the beach at the three analyzed sections was zero. Furthermore, as the crest height was progressively lowered and the AFP recalculated, the mean AFPs obtained from the new fragility representations were approximately bounded by those given by the existing representations. The consistency in the estimated vulnerability of the beach at these sites between the fragility representations provided confidence in the robustness of those adopted in current practice. When the methodology was applied to a site at the western end of Chesil Beach that possessed a substantially different geometry from the sites further east, considerable changes were observed in the estimated fragility compared with those obtained when the existing representations were assumed. Considering the intuitive superiority of the fragility representations derived in the article, it is reasonable to expect that the description of the beach's fragility provided by the new representations is likely to be closer to the true fragility of the beach than those given by their predecessors.
The relationship between the "hazard" and "consequence" components of risk are complex. Nevertheless, the substantial differences in a structure's estimated vulnerability between the new and existing fragility representations that have been shown in this article mean that it is likely that significant changes to the existing flood risk estimates would be observed at certain locations if the new representations were subsumed into an FRA. The observed differences in a structure's vulnerability according to the different representations thus provide justification for the adoption of more complex and computationally intense multivariate representation of the hydraulic load and fragility of coastal defenses in FRA. It is envisaged that more accurate FRA will have the capacity to deliver more efficient allocation of resources as well as improvements in the selection of mitigation measures for existing ones and in the designing of new coastal defense schemes.