The Modelling of Coastal Cliffs and Future Trends

Abstract

About 80% of the world’s oceanic shorelines include diverse types of cliffed and rocky coasts: plunging cliffs, bluffs backing beaches and rocky shore platforms. In combination, approximately 60% of the world’s population lives within 60 km of the coast. Rapidly retreating soft cliffs may be found worldwide and are particularly vulnerable to changes in the forcing factors. The study and analysis of the rate of change in shoreline position through time is important or even imperative for coastal management. The development of cliff erosion predictive models is mainly limited to geomorphological data because of the complex interactions between physical‐chemical processes acting simultaneously in time and space that result in large scale variations. Current historical extrapolation models use historical recession data, but different environments with the same historical values can produce identical annual retreat characteristics despite the potential responses to a changing environment being unequal. For that reason, process‐response models (PRMs) are necessary to provide quantitative predictions of the effects of natural and human‐induced changes that cannot be predicted using other models. Several models are explained and discussed, including a process‐response model, based on real data at Holderness Coast (UK).

1. Introduction

Approximately 60% of the world’s population lives within a band 60 km wide from the coastline towards coastal hinterland: the coastal zone. The “coastal zone” is defined in a World Bank publication as “the interface where the land meets the ocean, encompassing shoreline environments as well as adjacent coastal waters. Its components can include river deltas, coastal plains, wetlands, beaches and dunes, reefs, mangrove forests, lagoons and other coastal features”. The increasing world population, from 6 billion in 2000 to 9–10 billion by 2100, will lead to accelerating occupation of coastal zones even if these areas are seriously threatened by erosive processes [1–3], if so, an integrated management of resources and spaces is urgently required. Coastal recession and instabilities increases the risk on properties, people and their economic activities. Because of these and to improve coastal management and engineering intervention, it is critical to understand the coastal change pattern. The shoreline is dynamic in nature with different phenomena that acts at different temporal and spatial scales [4]. This area is also very sensible to changes in the conditions like sea‐level rise (SLR) or extreme weather events, which are mostly related to the acceleration of coastal retreat [5–7]. Some rapidly retreating areas like the Great Lakes of North America, England (north‐eastern coasts), and parts of the United States (north‐eastern coasts) are highly exposed to changes in these factors [8–12]. On rapidly retreating coasts, it is not just the geographic position of the cliff face that is important, but the entire interaction of the land with the constantly changing hydrodynamic regime creating a complex coastal system [4, 7, 13, 14]. Therefore, coastal managers and engineers must also consider the on‐ and offshore environment; weather; climate; the geotechnical properties of the slope forming materials; the influence of man‐made structures and also the effects of climate change and the associated sea level rise. This means that many factors control the complex behaviour and plan or profile shape of the coastline.

As well as assisting with the prediction of social implications of land loss, increasing the understanding of cliff erosion is of great importance in the wider science that underpins coastal engineering. However, in order to improve different coastal problems, the critical data are the shoreline position through time and/or its rate of change. These problems can be hazard zoning, developing appropriate management responses and coastal zone governance, assessing the possible effects of sea‐level rise and conceptual or predictive modelling of coastal morphodynamics [15–19]. Hence, appropriate and practical tools are needed to help forecast future shoreline evolution, especially where the shoreline responds in a non‐linear fashion due to non‐uniform and non‐homogeneous variations in the geology, the environment, the hydrodynamic regime and the changing climate. The tools, in terms of computational models, that are designed should be able to calculate an accurate approach of the coastal cliff state (position and shape) on the meso‐timescale (1–100 years), although longer term forecast may be useful or required for more long‐lasting structures. Although predictions have been based on historical rates of recession and on a variety of empirical and probabilistic methods [13, 18, 20], historical rates of change cannot be assumed to continue into the future, especially under changing conditions like geology, environment, hydrodynamic regime and climate [7]. Process‐response models (PRMs) [21–23] are needed to address these issues and provide quantitative predictions of the effects of natural and human‐induced changes, which cannot be predicted from statistical analysis of historic recession data or other models such as the Bruun‐rule model. For example, Bray and Hook [13] based on empirical methods like the “modified Bruun” model (also called historical trend analysis—[24]) related the erosion only with the sea‐level rise. The application of the latter model provides deterministic predictions of recession, without truly reflecting the uncertainty and process variability.

There are few published reliable process‐response models for cliff recession. These are based on functional relationships between the dominant physical processes covering the shore‐face, beach and cliff [21–23]. However, additional models are needed to simulate basal erosion and the resultant effects transmitted by gravitational movements up the cliff to the backscar, and this forms a major research gap [13], although progress is possible by adapting geotechnical stability analyses. Ashton et al. [25] developed a model called CEM (Coastline Evolution Model) focused on non‐cliffed coasts that includes a simplified geology, wave directions and sediment currents, and simply human activities. Walkden and Hall [22, 23] designed the first process‐response model (PRM) that includes individual modules to represent the evolution of soft rock‐shore profiles. These modules are, for example, rock erosion, beach protection, near shore wave transformation or sea level changes. In the model called SCAPE, at every time step, marine conditions are read from input files, and the new profile determines the system state. This whole view (holistic) is necessary to extract the information of the behaviour of complex systems. Trenhaile [9, 21] created new modules to explicitly incorporate waves and tides erosive mechanism like shear, abrasion or direct impact erosion. He also includes the beach morphodynamic on the erosive processes. Recently, the authors of this chapter [11] created a new PRM based on the idea of Kamphuis [26], and later Castedo et al. [27] developed another version of the marine treatment based on the work of Trenhaile [9, 21]. Castedo’s models can solve cliff recession problems, including geotechnical stability analysis, in real profiles with variable lithology.

Usually, quality of the model results is established against field measurements and historic catalogues. Since the recession processes are regularly slow on a human scale, the designated study areas, where the results will be validated, must have a sufficiently extensive data record to provide the necessary contrast and allow model error estimation. As long as the catalogue has more registered values, the quality of the estimation error is improved. Also as the recession rate is higher, the sample rate (or field campaign to get the data) should also be higher.

In some cases, for certain lithology, it might think that the speed of recession is a criterion for discriminating those coasts that are not cliffs. However, this is not so. Although clay dominated coastlines change more rapidly than most rocky coasts, they share the fundamental characteristics that distinguish them from beaches; those are that erosion and coastal retreat are irreversible. In this case, we do not suppose the traditional distinction that has been made between clay and rocky coasts because it is untenable [21].

2. Coastal recession processes

In coastal geomorphology, the cliffs (Figure 1a) are defined as a geographical feature in the form of denuded coastal rock or cohesive soil moulded by the action of two processes acting at the same time. The marine processes act as a transport and erosion force below the water level. The sub‐aerial processes determine the resistance force that produce the mass wasting phenomena [4, 28, 29], that is, landslides, topples, rockfalls, slumps, among others, of varying size distributions [20, 30]. Cliffs can be classified as follows: (a) steep to vertical slopes that expose hard or soft rock outcrops to the sub‐aerial and marine weathering actions; (b) the bluffs that are cliffs of smoothed and rounded profile in which the exposed materials have a lower resistance and cover with vegetation and soil the deep rocky basement. According to Sunamura [4], three major morphologies have been distinguished on rocky coasts: the type‐A shore platform or sloping type, the type‐B shore platform or horizontal or sub‐horizontal type, and plunging cliffs. The loss and erosion of the toe cliff by the wave’s action are responsible for the growth of the two types of shore platforms, though no appreciable recession occurs on plunging cliffs.

Figure 1.

(a) Cliff morphological features; (b) cliff behavioural unit (CBU) features and related processes and (c) spatial discretization and functions for process response numerical model.

The evolution of coastal cliffs is based on a balance between resisting and assailing forces. The lithology of the cliff material, principally described by its mechanical strength, determines the resisting force. The material resistance can be described by its compressive, tensile, cohesive or shear strength, being the first the most useful one [4]. The material resistance depends also on its formation conditions (bedding planes, cracks, joints, etc), the structure of the outcrops and latter of the human activity. The deterioration of its mechanical properties due to sub‐aerial weathering (physical, chemical and biological) or wave attacks can enhance the cliff erosion. Although the marine weathering processes, such that abrasion, strain, solution, bioerosion, freezing and thawing, occur on a small scale and are relatively slow, their results shape and modify the coastal topography. These actions are most important at the base of the cliff (Figure 1b), where the impact of the waves and the abrasive action of the water loaded with suspended solids and fragments of rock is frequent. In this way, a basal notch grows and deepens, the overhanging rock mass increases concentrating the shear stresses in the bottom of the notch. Once the strength of the rock is exceeded, cracks and fractures arise mobilizing, detaching and collapsing the upper rock mass. The repercussion of the slow erosion processes on a small scale is translated on a large scale and suddenly down‐slope movements in different ways (falls, topplings, slides, slumps, flows). Therefore, changes on coastal cliff morphology are not easily predicted because coastline recession is the cumulative result of successive numerous interacting gravitational phenomena, which are mostly coupled to marine hydrodynamics, intermittent and sporadic.

Observed complex behaviour arises naturally from several feedback patterns that control coastal erosion and recession, in the same way as has been usually presented in simple dynamical systems from physics theory. The system comprises two main components: the shoreline sub‐system and the cliff sub‐system. The influences on a coastal cliff system can be organized into a generalized hierarchy:

• Environmental activation mechanism: These are the elements that mainly control the recession problem like eustatic, vegetation, tectonic activity, isostasy, anthropic action, geology and climate

• Primary responses: They are the reactions (measurable) of the activation mechanisms that can produce changes in the coast like rainfall, wind regime, tidal‐wave conditions, sea‐level changes, morphology and geomechanical characteristics.

Coastal landsliding is a result of a combination of cliff toe erosion and geotechnical processes within the cliff, primarily connected with pore pressure distributions within the slope that influences stability [13]. Variations in groundwater level, which change the mobilized strength in soft materials, are significant contributors to the development of instability. The importance of groundwater level on the stability of coastal cliffs may be a significant uncertainty for predictions in view of a rapidly changing climate. In addition to long‐term changes in groundwater levels, potential changes in surges size and frequency are an important consideration in model development [9, 22, 31]. The more rapid removal of protective debris aprons and changing water pressure distributions contributes to more rapid instability and erosion. Therefore, any model must accommodate changing climatic conditions in order to provide adequately resilient assessments of hazard.

Cliff recession phenomena release debris deposits at the cliff toe, which, at least temporarily, provides a natural protection against the sea forces or further erosion. Refs. [32, 33] have highlighted the importance of this protective colluvium effect on the coasts of California (USA) and Holderness Coast (UK), respectively. These authors emphasized the importance of long shore currents, which may remove such material from the cliff base towards embayment or stagnation areas where could be deposited. The losses of this front unconsolidated detached material allow erosion to continue. Front detached material, like debris or densely fissured rocks to soils, has its critical limit strength reduced due to the fracturing caused by combination of the mass collapse and the weathering. Additionally, cracks and specific surface increase because size distribution of rock and earth fragments is moved to smaller sizes, so it increases their erodability. Therefore, less front material is consolidated and size reduced, less is the time to be eroded and spent at the cliff toe.

The backshore, the foreshore and the nearshore are all affected by the processes of coastal cliff recession; this can be grouped into a single element called “Cliff Behaviour Unit” (CBU). Each CBU unit consists of a 3D block of cliff coast that can be conceptually simplified and represented as a vertical profile (Figure 1c) with uniform homogeneous geologic and oceanographic conditions. This concept provides an important framework for cliff management in order to tackle the study and modelization of the complex recession process. Each vertical section is a reflection of the interrelationships between the morphodynamic processes and resultant changes in form over time in the coastline, the backshore, the foreshore and the nearshore. Subsequently, the recession process in the shoreline can be computationally solved into different CBU units. Then, for each CBU is reasonable to calculate the balance of forces and sediments. Later, the relation of both balances can be used to make a quantitative analysis of the complex erosional‐sedimentary dynamic regimes between CBUs. These units (CBUs) can be coupled to adjacent CBU’s within the framework to form littoral cells, which will allow the study of several areas together with more accurate results. The relations between CBU’s are important regarding the mass and energy flow exchange effects, that is, coastal landslides affecting the susceptibility of adjacent cliffs to further recession, sediments transport and others. The activation mechanism and their primary responses that determine recession in a CBU can be schematically represented through cause‐effect connections diagram (Figure 2) as a draft idea of the conceptual model acting on the whole CBU. This flow chart of interactions comprises soft and hard rocky coasts.

Figure 2.

Schedule of the cyclic evolution of a cliff profile due to the persistent incidence of the assailing forces of waves Fw, the erosion in the wave cut notch cause Fw is larger than the resisting rock force of cliff forming rocks FR, and the sudden instability and subsequent mass wasting which produces the cliff face recession.

In summary, coastline recession on unprotected cliffs or bluffs is characterized by a cyclic process of marine weathering and abrasion of rock or soil particles through multiple physical chemical mechanisms. As consequence, the removal of material from the cliff toe, and the increase of depth at the wave cut notch. The deeper the notch, the more important the stress concentration around it [32, 34]. Then, if the material strength is exceeded, slope becomes unstable, and a mass failure occurs that retreats the coastline at the cliff top. This can reduce the coastal slope and deliver debris to the front beach or shore platform at the cliff toe. The resultant deposit may only be partly removed by marine action allowing an accumulation of colluvium to form which provides a passive support to the slope, and some, albeit small erosion protection. As cliff recedes, erosion penetrates though zones of weakness in the cliff face and top such as cracks or joints, cutting clefts and crevices that may develop into caves, blowholes, which probably will evolve into bridges, arches, promontories, pillars, as common natural rocky coastal geomorphologies. The complete removal of the debris deposit or protective structures at the foot or front of the slope results in the cycle beginning again.

3. Previous coastal erosion‐recession models

In countless branches of science, computer models have helped in improving our understanding of dynamical behaviour in the earth physical processes. In the near shore, the design of computer algorithms tries to include process‐related information to understand flow dynamics, particulate and chemical distributions, and the overall oceanography of coastal environments. Whether cliff morphodynamics is considered, most of the main characteristics and activation mechanisms are illustrated in Figure 3. Most of them have been described individually or jointly in the literature, that is, coastal cliff recession system, groundwater load, geotechnical characteristics, waves and tides characteristics [21, 35–37]. Some of the very complex characteristics and their primary responses have been solved mathematical and numerically [21–22, 26, 38]; individually or coupled, side by side (Figure 2). However, as presented before, a coastal recession phenomenon is a consequence of the interaction between the ground and sea forces at the point where they meet. Scientific results that adequately reproduce such a complex land‐sea interaction have been scant and simplistic. Modelling results often presuppose that the recession pattern is the result of a single, unique distribution of waves, weather and environmental conditions. However, as different conditions, due to the range of behaviours found in different cliff morphologies (different CBUs), can cause diverse recession scenarios, the application of a simple prediction method cannot be universal [21]. Therefore, to develop a process‐response model, it is crucial to take into account a wide range of information and existing feedback relationships and multi‐scale coupling between the activation mechanisms and their primary responses that control the morphodynamics of the coastal recession events.

Figure 3.

Flowchart with details of the activation mechanisms and their primary responses.

3.1. The earliest models

First attempts at modelling this complex interaction of processes were made by Sunamura [4] and Kamphuis [26]. Both authors considered the erosion at a point on the cliff as an average rate that quantifies the amount of material loss by coastal recession over time (usually one year). Yielded values, as the rate (δ(x,t)) between marine driving and material resisting forces, could be considered as a safety factor from the geomechanical point of view [36, 37, 39] of the coastal erosion processes (Figure 5.2—chapter 5—[4]). Through laboratory experiments, some authors [40–42] validated closely first numerical approaches and realize that cliff change depends on spatial processes and time at different scales. Based on these data, and isolating the large number of elements involved in the recession event, an earliest model was developed. In this model, the change in position per unit time at a point x in the cliff face (erosion) is a function of: (i) time (t)—that control the number of tidal cycles that act a certain event; (ii) the assailing forces (F_W)—hydrodynamic and hydrostatic; and (iii) the resisting forces (F_R)—that depends on the lithology and structural conditions. The model is then based on the force ratio during a certain time (t) and quantifies the displacement of each point of the coastal section as, Eq. (1):

$\delta (x,t)=f(F_W(x),F_R(x),t)\leftrightarrow \frac{dx}{dt}\approx \frac{F_W(x,t)}{F_R(x,t)}$E1

Various attempts have been made to associate F_W and F_R with measurable, sea and rock, respectively, characteristics. Among others, the work of Craig [26] provided an expression in which the material strength is represented by a factor to be calibrated experimentally. Subsequently, the model by Sunamura [4], derived from experiments on a concrete wall in an artificial wave tank, suggested that F_R is directly proportional to the uniaxial compressive strength (UCS) of exposed rock. Also Sunamura [4] suggests that the erosion rate at a certain time δ(x,t) is proportional to the logarithm of the ratio between both forces, F_W/F_R, applied at the position x. Resistance forces in this context means the strength and hardness of rocks attacked by the physical forces of marine erosion, or their durability in the face of physical, chemical and biological weathering process. Rock resistance depends on several factors such that lithology, joint or fracture density, and the existence of bedding or cleavage planes which decrease rock hardness or soil cohesion. Identifying the most appropriate parameter to represent the strength of the material to erosion has been investigated in numerous laboratory studies suggest the compressive strength σ_c (Japan island: andesite lava, sandstone, clays, gravels, among others—[4]; Cilento area: sand, calcilutites and siltstones—[38]) the most appropriate parameter. Although less frequently, the tensile strength (for idealized cliffs in chalk—[34]), cohesion (Calvert Cliffs, Maryland: alternation of diatomaceous earth, sand, clay and marl—[43]), shear strength (Lake Erie: glacial silt/clay overlain by lake sand—[26]; Grates Lakes: cohesive clay coast—[21]) or Young’s modulus (Fukushima coast: gravel, sandstone, mudstone—[44]) have also been proposed. Nevertheless, the model that most completely describes the resisting forces is the Budetta et al. [38]. The authors introduce into Sunamura’s [4] model the Rock Mass Index (RMi), taking into account not only the UCS value of the intact rock but also taking into account the block volume and the fracturing state (joint condition factor). Their results show this parameter capture the relationship between the resisting forces and their geotechnical controls. Incorporating a suitable representative form of rock mass (i.e., the rock material plus the fractures) was an important step forward in the prediction of recession.

3.2. Historical recession based models

Historical trend analysis [24, 45] uses the expression based on the ‘Bruun Rule’ [46] and modifications to the Bruun rule [13], on two‐dimensional profile models, to predict future shoreline erosion depending on projected future sea‐level rise scenarios and measured shoreline recession rates. The main idea behind the Bruun Rule is the statement that the upper shore face (the morphologically active—years—part of the continental shelf) has a profile that is in equilibrium with the hydrodynamic forcing (i.e., wave driven processes mainly, but also density driven currents close to river or estuarine entrances) and that for any perturbation in the forcing, the reaction of the profile will be moderately fast. When sea level rise occurs the shoreline retreats and a new equilibrium profile will form at the new shoreline position by moving sediments to deeper water, or in other words, the profile is shifted to an upward and landward position. Bruun [46] developed a model that evolved as a result of sea level rise (S) until it reached the shoreline retreat equilibrium, R = S×L/(h+B). Where L is the cross‐shore width of the active profile, h is the closure depth, and B is the elevation of the beach or dune crest. Even though the Bruun Rule is very often used, it has a number of limitations, and care must be taken when using this tool to predict the coastline response to sea level rise. It basically tends to overestimate the future recession rate.

The following equation [24, 45], based on the Bruun Rule, predicts future shoreline erosion by using projected future SLR scenarios and measured shoreline erosion rates R₂ = S₂×(R₁/S₁), where R₂ is future shoreline recession rate (m/year); S₂ is projected future sea‐level rise rate (m/year); R₁ is measured shoreline erosion rate (m/year), and S₁ is sea‐level rise rate (m/year) during the period of measured shoreline erosion. The R₁ could be obtained from the Digital Shoreline Analysis System (DSAS) [47] fitting a least‐squares regression line to all shoreline points for a particular transect. The historic rate‐of‐change of position of cliff‐top is the slope of the linear regression (m/year), also known as historic rate of change. Recent approaches to historical shoreline analysis have been made using the Digital Shoreline Analysis System (DSAS) version 4.3 [47]. Additionally, statistical models [14, 24] should be used to predict future shoreline locations using historical information extracted from the DSAS. Within DSAS, the linear regression method is usually applied to obtain rates of coastline change statistics, computationally done in a GIS environment, for a time series of shoreline vector data that describe temporal evolution of the coast from multiple shoreline positions [47]. This method is selected because it includes the following features: first, all the data are used; secondly, it is also the most commonly applied statistical technique for expressing rates of change [48]; and thirdly, it includes additional statistical information essential to be compared with the yield results from other analytical methods (i.e. [14, 24]).

3.3. Probabilistic forecast methods

A key feature of these methods is the recognition of the episodic nature of cliff recession at many coastal sites. In other words, cliff recession proceeds primarily via occasional landslide episodes, non‐uniformly distributed in size, followed by periods of relative inactivity, which may last for more than 100 years on some coastlines [49]. The process is complex and far from purely random, so maybe probabilistically distributed in some cases.

This methods use to present significant limitations, but on the contrary, the projection of historical rates is probably the most evident and straight forward approach. However, the historical record uses to consist of a less than 10 measurements over the last 100 years. This is insufficient to explain the pattern of recession events, and so the cumulative process of erosion, rock tension, and finally failure. In this case, the historical record can only reveal a limited picture of the past events.

The scarcity of historical cliff position data can limit the usefulness of many conventional statistical methods, such as linear regression. One approach to addressing this problem is the development of probabilistic models [50] to simulate the recession process. These models consider the cliff recession as an episodic random process and use a Monte Carlo sampling tool. This model can be used to simulate synthetic time series of recession data, which conform statistically to the cliff recession measurements. If the cliff recession process is plotted against time, its episodic form is reflected in a stepped function. If the simulation is repeated several times, a probability distribution of cliff recession can be built. The simplest model can be defined by two distributions: a description of the timing of recession events, and an event size distribution describes the magnitude of recession events in terms of the mean size and their variability.

For different cliff settings, a different probabilistic method can be used. The selection of the method should be based on: the knowledge on cliff recession of the study area, the nature of available data, if the primary responses are expected to remain constants and the amount of investigation and analysis which can be justified.

The type of information available is influenced by the nature of the cliff recession process. Also the results of any model will be influenced by the information available. For these reasons, more than one method must be used to indicate, in some way, the robustness of the predictions. Even in some cases, as the quality and quantity of data increases, more sophisticated analysis must be re‐done.

3.4. Process‐response models

The factors and phenomena contributory to development of many different geological processes can be systematized in term of process‐response (PR) or cause‐effect models. In the case of sea cliff systems, they are based on experimental‐theoretical knowledge and mathematical description of the interactions between the hydrodynamic characteristics, the coastal processes, geotechnical properties and mechanical behaviour of the rock mass. Such models are still at an early stage of development, and most were originally purely deterministic. However, in order to incorporate uncertainty in the characteristics and variability of physical phenomena, a parametric probabilistic framework can be included addressing the complexity of the system using simple stochastic models of each considered process. In these models, as in previous ones, historical records were used for calibration. Also, in this type of model, most of the physical interactions are studied at laboratory scale so upscaling them remains problematic.

Walkden and Hall [22] developed a numerical model, SCAPE, to simulate the processes on a cohesive shore profiles. SCAPE is mainly based on the model of [26]. However, they added the effects of waves by a cumulative effect of an erosion function per wave developed through a tidal period, and the effects of coastal development platform, including changes in shore slope. The mathematical solution adopted means that the shore slope is limited to 10°in order to avoid numerical instabilities that produce non‐natural roughness. Their erosive function was based on erosion experiments by Refs. [41, 42]. The weakness in this model is that it does not consider the characteristics of the rock, so that the resistive capacity of the ground is adjusted through the calibration constant based on the availability of historic recession data. Similarly, the failure criterion used is a deterministic approach, after 10 recession events, the model removes any overhanging material from the most eroded point in the notch, which produces a vertical cliff. This severely limits the use of the model, since the time when a break occurs or not, greatly affects the final appearance of the erosion profiles.

Trenhaile [21] presented a model to examine the relationships between variables that affect the erosion and development of cohesive clay coasts. Trenhaile introduces a relationship to consider the effect of the abrasion due to the suspended material. The erosion rate within the model is introduced as the amount by which the applied bottom shear stress produced by the waves, exceeds the critical shear stress for erosion. This model produces shore‐normal submarine profiles retreat but it not reproduces the upper part of the section that is the cliff face.

4. Coastal cliff model development

The process‐response models (PRM) are the most modern and updated software to incorporate critical activation mechanisms and subsequently trigger the primary responses in a CBU. In order to reduce the design complexity in the model, the process in the cliff‐shoreline system is represented by single modules as in Figure 3. Each one is described mathematically by a single numerical approach to improve computational efficiency, high accuracy of results and assist understanding of emergent properties of the systems behaviour.

The recession model described here is validated in coasts with rocks with a compressive strength σ_c that varies from very soft (geotechnically named very weak, i.e. <1.25 MPa) to moderately hard (i.e., weak to moderately strong, 1.25–50 MPa). These kinds of rocks are also known as soft erodible materials. After the eroding process explained previously, the crest of the cliff moves inland a quantity R(t), which usually varies at each failure event. This is the point of interest in determining the risk to cliff‐top assets. So, at this point here is located control, and monitoring observations are regularly collected to determine the instantaneous position of the cliff, which is the point. For more details of model applicability and development see Refs. [11, 12].

The model explained here is the first one that explicitly includes a slope failure mechanism using a limiting equilibrium procedure in a dynamic manner, Fellenius’ method of slices [39]. Previous models like [21, 22] uses an expert judgement criteria where the authors assume that each failure event occurs at every ten recession events, or a probability distribution to describe the possible cliff top location following a failure event [51].

4.1. Erosive model details

As outlined previously [11, 12] in this numerical model, the rock profile, y (z, t + T) represented as a column of horizontally aligned layers of height Δz = 0.05 m, moves a quantity δy(z,T) in the points z affected by the erosive processes after one tidal period T. This quantity can be defined as the erosion of each element of the cliff front and depends on H_b, which is the breaking wave height; T_b, which is the wave period; K, which is a calibration term representing hydrodynamic constants (100 m^13/4s^7/2/kg for the application site; see Refs. [11, 12]); σ_c(z), which is the uniaxial compressive strength (USC) of the rock mass; ∂_zy(z,t)⁻¹, which is the slope of each rock element z (which changes through simulation time in response to the calculated erosion and consequently to profiles evolution); p_w(z,t), which is the shape function (which include from experimental data—[41, 42]—the erosion of different type waves) and w_t(t), which is the tidal expression introduced as a sinusoid that oscillates about mean sea‐level with an amplitude A_m and period T (both can be obtained through governmental marine agencies). The rock profile evolution can be described as seen in Eq. (2):

$y(z,t+T)=y(z,t)\pm \frac{H_b^{13/4}\left(t\right)T_b^{3/2}\left(t\right)}{K{\sigma }_c\left(z\right)}{\left(\frac{\partial y(z,T)}{\partial z}\right)}^{-1}{\int }_0^T{p}_w\left(w_t\right(t)-z)dt$E2

The in situ cohesive material can be represented by its lithological characteristics if a field survey has been carried out. They are incorporated into the model using the value of the uniaxial compressive strength σ_c‐mass(z) of rock or soil mass calculated at each elevation z. In order to estimate the value of UCS, σ_c‐mass(z) of rock mass, the Hoek‐Brown failure criterion was used based on rock mass rating (RMR) and estimated values of the constants m and s [52]. Alternatively if the coefficient of internal friction angle φ’, and the cohesive strength c’ values for the Mohr‐Coulomb constitutive model are known it is possible to estimate the value of UCS σ_c‐mass(z) of the rock mass (Eq. (2) from [52]). In addition, the original material of the cliff composes the colluvium, but their original soil structure can be modified controlling the mechanical properties of this new structure. So the considered value of UCS for the colluvium σ_{c‐weathered}(z) can be obtained through the UCS for remoulded materials σ_cR, the internal friction angle ${{\varphi }^{\prime }}_w$, and the cohesive strength ${c^{\prime }}_w$, following the same methodology as before. In both cases, when the necessary data to obtain the values of UCS σ_c‐mass(z) or σ_{c‐weathered}(z) following the method presented by [52] are unknown, the considered values are σ_c and σ_cR, respectively.

4.2. Extra modules (primary responses)

Erosion, cliff stability and sediment transport are calculated one per tidal cycle T which is considered as one time step in the computational procedure of the global simulation time. As the foot erosion progresses, the notch becomes deeper and deeper increasing the volume of the overhanging material. Each tidal cycle geotechnical equilibrium conditions are evaluated to determine the stability of the exposed rock mass. This stability will be higher or lower as a function of the fracture grade of the rock mass, the material strength and the groundwater level. The first two parameters are introduced into the model by the RMR. When the retreating forces overcome the resisting ones, the rock strength is exceeded, resulting in the material critical failure (that can be deposited or removed) causing the loss of some coastline and consequently coastal recession.

Usually, for practical purposes, the groundwater conditions are based on the position of the groundwater table, ranging from fully drained to saturated. It is defined by the ratio where the phreatic surface coincides with the ground surface at a distance behind the toe of the slope among the slope height [37]. In this model, the water table elevation is described as the Dupuit parabola based on the upstream head, measured from the impermeable base, and the downstream head data, see Table 1 for details [35]. The inland pressure is considered here as a difference of the depth from the surface to the seepage point b, minus the depth from the surface to the groundwater level a. The downstream pressure is considered as the reference zero‐level at the spring point, see Figure 4.

Figure 4.

(a) Schematic cliff section showing the adapted Dupuit equation solution for flow through a porous media; (b) schematic visualization of the factor of safety calculation following the Fellenius method; (c) search criteria for circular arc slip surfaces from (b); (d) schematic visualization of the maximum bending stress situation.

Factor of safety to slumps (FSslump)	$F{S}_{slump}=\frac{{\sum }_{i=1}^n({c^{\prime }}_i\Delta l_i+(W_i\cos {\theta }_i-u_i\Delta l_i\left)\tan {{\phi }^{\prime }}_i\right)}{{\sum }_{i=1}^n{W}_i\sin {\theta }_i}$
Variation in groundwater elevation (h)	$h\left(y_D\right)=+\sqrt{(y_D/L){(b-a)}^2}$
Maximum bending stress (σt max)	${\sigma }_{t\max }=M/Z$
Bending moment (M)	$M=1/2\gamma d{h}_{failure}{L}_n^2$
Modulus of the section (Z)	$Z=1/6d{h}^2{}_{failure}$

Table 1.

Equations to calculate the stability at every calculation time step T [11, 12].

where n is the number of slides, i is the i slide, c_i’ and φ_i’ are the effective cohesion and friction angle, respectively, Δl_i is the length of the slide, W_i is the weight, θ_i is the inclination angle of the bottom of the i slide, u_i is the pore water pressure at the centre of the base of the i slide, y_D is the horizontal distance to the seepage point to any point in the parabola, L the horizontal distance between the seepage point and the point where the phreatic surface was measured, b the seepage point, a the phreatic level measurement,γ is the material unit weight, d is the distance along the bluff face of the overhanging material, h_failure is the height, and L_n is the notch depth.

Clay cliffs are susceptible to deep‐seated rotational landslides, which occur where basal erosion is rapid enough to remove the debris and steepen the underlying coastal slopes [28]. This model utilizes limit equilibrium techniques; for simplicity, accuracy and computational time, the ordinary method of slices (OMS) technique is used [39]. The OMS method was developed by Fellenius in 1936 and is sometimes referred to as “Fellenius’ Method.” Details of the equation used and the schematic representation can be seen in Figure 4 and Table 1. In the OMS method implemented in this model, the failure circle is defined by the position of a point through which the circle must pass at the toe of the slope and the tangent angle (α) to the circle plane at this point. Also is necessary to define the point at which the circle intersects the cliff top (point 2 in Figure 4c). The variation of the second point (2 to 2′) in the x‐axis allows the code to look for unstable circles. The maximum depth, from the cliff edge, to look for an unstable circle is h_talus based on physical observations [36, 39]. Also, different circles are created by separating 2 to 2′ a discrete quantity Δc selected by the user. If the value of factor of safety (FS) for each slip surface is less than 1.0, slope is unstable. The factor of safety obtained with the “Fellenius’ Method” is the lowest when comparing with other techniques, so the results obtained with this method are on the safe side [39].

In addition to the OMS method, this recession model includes a notch collapse stability analysis based on the ideas proposed by Refs. [53, 54]. The cantilever beam model derives the distribution of stress inside the cliff assuming that tensile stress acts on the upper section and compressive stress on the lower section, see Figure 4d. The maximum bending stress (σ_t‐max) is calculated as shown in Table 1. The tensile strength (σ_t) is specified and limited, which in many analyses is taken to be 10% of the rock mass cohesion [37], for isotropic materials. Then, a collapse occurs when the maximum bending stress (σ_t‐max) exceeds the tensile strength of the rock (σ_t).

Following a failure event disintegrated collapse material is deposited at the foot of the cliff, which acts as a natural protecting material. This talus material now located at the foot of the slope is more altered and fractured than the original one, and, therefore its resistance to erosion is reduced, this effect becomes more noticeable at smaller values of σ_{c‐weathered} (z). It is usually quickly removed by the cumulative action of waves and currents. In the model, debris produced from the cliff face are deposited on the colluvium slope and covers the foot of the cliff, in the proportion defined by the user. Three different solutions are developed: the colluvium shape is determined by the friction angle for weathered materials ${{\phi }^{\prime }}_w$ [36, 37, 39]; the angle of reach for landslides, defined from the detaching point with respect to the horizontal, as the ratio of vertical height differences between the source point and the furthest point reached in distance to the horizontal [55, 56]; and based on expert field observations where the talus is created from the half cliff height (h_talus/2) with a small slope [7, 28].

5. Model application: Holderness Coast (UK)

This model has been applied within the cliffs of the Holderness Coast (UK), which stretches from Flamborough Head in the north to Spurn Head in the south, presenting a smooth “S‐shape.” The area is sparsely populated, and land use is predominately agricultural and tourism; however, critical infrastructures (transport and communication networks) are located near the coast [33] (Figure 5).

Figure 5.

Location and geological sketch of the Holderness Coast. Note that, the numbers are to situate the lithological columns through the area.

This area is one of the youngest natural coastlines of England, a 61 km long stretch of low glacial drift cliffs ranging from 3 to 35 m in height. The Holderness Coastline is made up of soft boulder clays (tills) left after the retreat of the Devensian ice sheets about 12,000 years ago. These soft, recent deposits sit on a platform of chalk which slopes away gently to the east originating from between the Santonian (85.5 Ma) and the Maastrichtian (65 Ma) stages from the Upper Cretaceous period. The chalk is located at a depth of at least 30 m below the commonly used UK datum (Ordenance Datum Newlyn—OD), across the vast majority of the coastline [57]. The most recently and accepted division of the glacial materials was summarized by [58], who divided the area into: Basement, Skipsea, and WithernseaTills. The oldest unit (Basement Till) is the lowest of the stratigraphic series. The following of the series is the SkipseaTill, and finally, the Withernsea Till is found. In addition, the majority of the coastline is also capped by a weathered till deposit. The cliff stability alongside the Holderness Coast has been the subject of study and controversy over the last 30 years. Field and laboratory testing have shown that the tills display essentially similar geotechnical characteristics used as an input parameter (see Refs. [11, 59–61]).

At Holderness, the environment is mainly dominated by wind wave development. The dominant wind direction is north‐easterly, determining also the wave direction, and creating a north‐south orientated net long shore current. During storm events, highest waves can reach up to 4 m. In addition to this, tidal range at the Holderness coast is very high and can reach up to 7 m.

Since the Roman times, many villages and parishes have been lost [8], with an estimated area of more than 200 km². Along the entire coast [33] an average recession rate of 1.55 m/year has been recorded, making this coast one of the fastest eroding coastlines of Europe. Due to this erosion rates, the East Riding of Yorkshire Council designates important economic resources to collect cliff top retreat and coastal cliff profiles measurements within the entire coastline since 1951. A common conclusion for the rapid erosional behaviour of the Holderness Coast is the low cohesive values of the glacial till and the groundwater.

The model requires input values for oceanographic data, sea level rise estimations, retained material that form a talus after a failure event, the geotechnical data and the groundwater level. Wave height in deep water (H₀ = 0.98 m) and mean breaker period (T_b = 7.66 s) are entered as a mean values from 3 years data registered by the Hornsea Directional Waverider Buoy, deployed on the 05/06/2008. The tidal amplitude (A_m = 3 m) and tidal period (T = 12.46 h) are obtained from The National Oceanography Centre. The hydrodynamic constant (K) had a value of 10² for the entire Holderness coastal area. The groundwater level is also an important factor in the entire system with special influence on the failure mechanism. For modelling purposes of the Holderness area, the groundwater level used is 2.5 m under the surface at about 250 m landward. The seepage point can vary from one profile to another and usually require specific site investigation [33].

A changing climate and particularly an accelerated sea‐level rise are believed to impact cliff retreat rates [3] with an expectation that shoreline retreat rate will generally accelerate in the future [11–13, 16, 17]. According to United Kingdom Climate Impacts Programme (UKCIP), relative net sea level will rise between 0.002 m/year and 0.01 m/year for Yorkshire by the 2080s for the low emissions and high emissions Scenarios, respectively [62]. However, Walkden and Rossington [63] recommend that for planning purposes, a relative rise of 0.006 m/year is taken into account.

5.1. Historical analysis based on DSAS

Table 2 shows the historical recession rates of change calculated using the real profiles provided by the East Riding of Yorkshire Council [64] and the artificial profiles created in ArcGIS to carry out this work for the 159 year period 1852–2011. The minimum R² value of these regressions is 0.89, see Ref. [12] for more details. The recession rates vary from one profile to another with values from 0.95 to 1.83 m/year, making the erosion trend so irregular. These variations can be due to numerous factors such as sub‐horizontal or vertical discontinuities, cliff height disparities or shore platform irregularities.

Profile	LRR (m/year)	LCI (m/year)
P21	1.83	0.63
22N	1.73	0.45
P22	1.76	0.24
23N	1.73	0.31
P23	1.65	0.44
24N	1.57	0.70
P24	1.46	0.54
25N	1.43	0.45
P25	1.51	0.60
26N	1.41	0.56
P26	1.46	0.69
27N	1.37	0.84
P27	1.20	0.62
28N	1.10	0.51
P28	1.13	0.68
29N	1.08	0.77
P29	1.05	0.57
30N	1.13	0.41
P30	1.05	0.39
31N	1.17	0.34
P31	1.14	0.44
32N	1.18	0.43
P32	1.10	0.24
33N	0.95	0.26
P33	1.03	0.26
34N	1.10	0.23
P34	1.13	0.30
35N	1.13	0.25
P35	1.10	0.23
36N	1.02	0.32
P36	0.95	0.28
37N	1.07	0.44
P37	1.03	0.32

Table 2.

Historical recession rates for several profiles measured in the field (noted as P followed by a number that follows the numbering provided by Ref. [64]) and artificial profiles obtained from DSAS analysis (that follows the same numeration as the real profiles with the addition of a letter N) for the 159 year period 1852–2011.

LRR, linear regression rate; LCI, standard error of LRR at 99.

5.2. Models output

Another way of checking the goodness of the PRM developed is comparing its output, with shorelines predicted by Lee and Clark [14] and Leatherman [24] models. Both models are used with the data extracted from the DSAS software. The PRM has been run with two different scenarios, one with a SLR of 0.002 mm/year to cover the lower prediction, and the other one with a SLR of 0.006 mm/year to cover the opposite (Figure 6). After that, future recession values are calculated in several epochs (2011, 2021, 2031, 2041 and 2051), and a shoreline was estimated in each one.

Figure 6.

Predicted cliff lines for PRM: (a) Skipsea Sands; (b) Far Grange region. The cliff edge marked line represents the 2011 cliff‐top position.

The Leatherman model [24], also called historical trend analysis, uses the equation presented in the point 3.2: R₂ = S₂×(R₁/S₁). Here, R₁ is the average recession rate (LRR) and S₁ the historical sea level rise equal to 0.002 m/year [1]. Finally, the S₂ is the projected future SLR considering two scenarios: 0.002 m/year (lower bound) and 0.006 m/year (upper bound), see Figure 7 for details.

Figure 7.

Predicted cliff lines for Leatherman model: (a) Skipsea Sands; (b) Far Grange region. The cliff edge marked line represents the 2011 cliff‐top position.

The Lee‐Clark model [14, 17] calculates the future recession as follows: recession by year A = (LRR ±LCI) × T. where LRR and LCI values are obtained from the DSAS analysis and for the study case can be found in Table 2. To cover the upper and lower bounds, the LCI can be subtracted or added to the LRR value, respectively. T is the time to extend the simulation. The output can be seen in Figure 8.

Figure 8.

Predicted cliff lines for Lee‐Clark model: (a) Skipsea Sands; (b) Far Grange region. The cliff edge marked line represents the 2011 cliff‐top position.

5.3. Models comparison

The study presented here has quantified the recession of part of the Holderness Coast shoreline, from profile P21 to P37, located between Bridlington and Hornsea. The recession values obtained with the three models presented here (Figure 9) are in accordance with data obtained by recent analysis based on historical data for the same location [10, 65]. The Leatherman model presents the highest recession rates and higher uncertainties as the shadowed region (between lower and upper bound) is the biggest. The Lee‐Clark model, based on historical analysis is the most optimistic in terms of low erosion; however, it is based on hundred years of data. In the middle, we found the PRM model which is less sensitive to the change of the sea level as higher values can decrease the cliff face height increasing its geotechnical stability. Also, this latter model respond better to other changes that may arise during time like changing in profile shape, marine conditions, etc. This empirical model can provide deterministic predictions of recession but in reality cannot reflect information of the process variability and uncertainty. Obviously, this is a problem as cliffs with the same recession rates can have a totally different behaviour. On the contrary, these models are of friendly use and “only” need historical data, while the PRM model needs geological, marine, topographical and also historical recession data.

Figure 9.

Predicted cliff lines for PRM, Leatherman and Lee‐Clark models at year 2051: (a) Skipsea Sands; (b) Far Grange region. The cliff edge marked line represents the 2011 cliff‐top position.

6. Final statements

The estimation of a future shoreline is a complex problem in engineering, due to the stochastic nature of cliff retreat, which is apparent from the irregular cliff lines. The most widely models used to forecast erosion is through use and analysis of historical records. The methods to predict the future cliff‐top evolution exposed here rely on the knowledge of historical recession rates, especially the empirical models, such as the Lee‐Clark model [14] and the Leatherman model [14]. The model presented here has been applied so far on shores with soft, erodible material, where erosion is dominated by wave action; however, with the introduction of the rock mass UCS, the model range applicability can be extended to hardest materials. The review and results provided here can thus be used to help the developers in the process of deciding how to manage the potential losses (social, ecological and economic), and when and where to intervene in the erosion processes at least for rapidly eroding coasts like Holderness, UK.