Evidence of slow millennial cliff retreat rates using cosmogenic nuclides in coastal colluvium

The erosion of rocky coasts contributes to global cycles of elements over geological times and also constitutes a major hazard that may potentially increase in the future. Yet, it remains a challenge to quantify rocky coast retreat rates over millennia – a time span that encompasses the stochasticity of the processes involved. Specifically, there are no available methods that can be used to quantify slow coastal erosion (< 1 cm yr⁻¹) averaged over millennia. Here, we use the ¹⁰Be concentration in colluvium, corresponding to the by-product of aerial rocky coast erosion, to quantify the local coastal retreat rate averaged over millennia. We test this approach along the Mediterranean coast of the eastern Pyrenees (n=8) and the desert coast in southern Peru (n=3). We observe a consistent relationship between the inferred erosion rates and the geomorphic contexts. The retreat rates are similar, 0.3–0.6 mm yr⁻¹ for five samples taken on the Mediterranean coast, whereas two samples from vegetated colluvium have a lower rate of ∼ 0.1 mm yr⁻¹. The coastal retreat rate of the Peruvian site currently subject to wave action is similar to the Mediterranean coast (0.5 mm yr⁻¹), despite Peru's more arid climate. The other two Peruvian sites, which have not been subjected to wave action for tens of thousands of years, are eroding 20 times more slowly. The integration periods of the two slowest Mediterranean coast erosion rates may encompass pre-Holocene times, during which the sea level and thus the retreat rate were much lower. We explore here this bias and conclude that the associated bias on the inferred retreat rate is less than 80 %. These data show that rocky coasts are eroding 1 to 20 times faster than catchments in the same regions on average over the last few thousand years. We anticipate that this new method of quantifying slow rocky coastal erosion will fill a major gap in the coastal erosion database and improve our understanding of both coastal erosion factors and hazards.

1 Introduction

Rocky coasts, which represent approximately 50 % of the world's coastline (Young and Carilli, 2019; Regard et al., 2022), are eroding by wave action (Sunamura, 1992; Trenhaile, 2002) through processes not fully understood yet (Prémaillon et al., 2018). Coastal erosion likely contributes to global cycles of elements (Regard et al., 2022) and also constitutes a hazard with implications for infrastructures and economy that is possibly amplified by global climate change (e.g. Ashton et al., 2011). Both aspects require documentation of coastal cliff erosion over different time spans, including millennial timescales; this is something that remains a challenge today. Coastal erosion occurs irregularly and varies over several timescales, namely daily, depending on wave energy (tide); annually, depending on the frequency and intensity of storms; over decadal to centennial periods because of rising sea levels linked to climate change; and over multi-millennial periods, depending on relative sea level changes. Landslides can also affect coastal erosion and participate in the stochastic nature of the retreat over decadal to centennial time spans (Prémaillon et al., 2018). Thus, there is currently a growing effort to quantify the erosion rates of rocky coasts (Prémaillon et al., 2018). The two most commonly used methods are the comparison of aerial photos (e.g. Dornbusch et al., 2008; Hapke et al., 2009; Letortu et al., 2014) and, for the past 15 years, the comparison of 3D point clouds collected by photogrammetry or lidar (see Rosser et al., 2013; Dewez et al., 2013; Prémaillon et al., 2021; Swirad and Young, 2022). Lidar offers very good resolution, but its maximum time span does not integrate all the stochastic natures of erosion, which in some cases is achieved by collapses with a long (> 100 years) return time (Dewez et al., 2013; Dornbusch et al., 2008). Aerial photography can increase the temporal range to more than 50 years (Letortu et al., 2014; Dornbusch et al., 2008) but is then limited by the resolution, and therefore, this technique is ineffective for capturing low retreat rates.

The lack of a reliable method to document millennial rates has led to the development of retreat rate measurements based on the concentration of cosmogenic isotopes in rocks from the shore platform (Choi et al., 2012; Regard et al., 2012; Rogers et al., 2012). These measurements have proven effective in quantifying millennial velocities between 10 and 300 mm a⁻¹ but require accurate measurements of the platform shape (Regard et al., 2012; Hurst et al., 2016; Swirad et al., 2020; Duguet et al., 2021; Shadrick et al., 2021; Clow et al., 2023). In addition, such a measurement has never been attempted on rocks more resistant to wave action for which shore platforms are narrower with complex geometry, limiting the possibility of providing evidence of slow coastal retreat rates. Given that the available methods are adapted for rapid retreat rates, it is possible that the current database of rocky coast erosion rates is biased towards these rapid values (Prémaillon et al., 2018). As a result, complementary methods are much needed.

In this work, we introduce a new method for quantifying millennial cliff retreat by analysing cosmogenic isotopes in colluvium. Cosmogenic isotopes have already been used to quantify the erosion of escarpments, either from river sediments to derive an average erosion rate for the catchment draining the escarpment (e.g. Wang and Willett, 2021; Stokes et al., 2023) or from local sampling of outcropping bedrock (e.g. Cockburn et al., 2000; Heimsath et al., 2006). The main difference with these previous methods is that the method developed here applies to a coastal scarp and that this scarp is not dissected by rivers. This method is derived from a method that is already being implemented to measure the retreat rate of alluvial valley sides (Zavala et al., 2021). With this method, we obtain slow cliff retreat rates between 0.05 and 0.5 mm a⁻¹. The integration time of this method is of the order of the inverse of the retreat rate, namely between 2 and 20 kyr. These periods are longer than the rarest storms and collapses and therefore average out the stochastic phenomena. We therefore believe that this approach can provide a new tool for quantifying the erosion of rocky coastlines over these time periods.

2 Method and sampling

The proposed method is based on the assumption that 1 kg of colluvium sampled at the foot of a rocky coast includes grains detached all along the coastal scarp above it (Fig. 1). The sampled colluvial wedge is active, indicating the ongoing erosion of the cliff. The mean ¹⁰Be concentration of this colluvium sample can be converted into the mean erosion rate of the rocky coast, assuming a secular steady state in this concentration. This principle is nothing more than one of the assumptions that underlines the widely used method using ¹⁰Be in river sands to quantify the average erosion rate of catchments (von Blanckenburg, 2005).

Figure 1 (a) Shore platform/retreating cliff sketch showing the sampling strategy on colluviums at the foot of a coastal cliff. The sketch illustrates a multi-point sampling of the colluvial deposits described in the text. The sampled colluviums are active, suggesting that the erosion of the cliff is ongoing. (b) View from the southwest of the coastal cliff at the BRAV3 site. The cliff is 80–100 m high. (c) Sampling of one of the colluvium samples at the BRAV3 site. (d) Sampling close-up view; the finest part of the grains is sampled until reaching roughly 1 kg of grains from several colluvium samples.

Each study site is composed of a shore erosion platform backed by a retreating cliff (Fig. 2; see the Supplement for complete site descriptions). The cliff or escarpment face has a roughly constant slope. The erosion of the cliff produces colluvium that lies at its base. In detail, the cliff is divided into segments along the coast, each with a specific source area that feeds the colluvium, which is defined between the cliff top and the cliff foot. Three criteria were used to select the sampling sites. (i) The lithology must guarantee the presence of quartz grains in the colluvium. (ii) The geomorphological context must limit, as much as possible, the contribution of sediment coming from the areas above the cliff, which allows us to constrain the source of the colluvium. To do this, we selected portions of the coast where the cliff top also corresponds to a drainage divide or which are located at the front of ridges between two rivers. In the latter case, the probability of a grain coming from above the linear face of the escarpment is minimal. (iii) Last, it must be possible to access the foot of the cliffs in order to carry out the sampling.

Figure 2 Sampling sites of this study. Top: general location of the series of sampled cliffs (Mediterranean eastern Pyrenees on the left and southern Peru on the right) using © Google Earth views. (a) Location of the VERM series of samples on the 5 m horizontal-resolution DEM (French Institut Géographique National). (b) Location of the BRAV series of samples on the 5 m horizontal-resolution DEM (Institut Cartogràfic i Geològic de Catalunya). (c) Location of the COSTA series of samples on the 30 m horizontal-resolution DEM (SRTM1). Note the presence of the uplifted platform at the foot of the cliff for COSTA2 and COSTA3. Contour lines are traced every 20 m in all topographic maps.

The studied escarpments are a few dozen metres high on the Pyrenean coast and a few hundred metres high in Peru. The cliff surface is covered by a regolith in Peru, and in the Pyrenean coast, a thin regolith alternates with outcropping bedrock. Sampling was conducted using the method described in Zavala et al. (2021). For each sampling site, about five colluvium samples were collected at the surface of debris wedges along a 50 m stretch of the cliff and then mixed together to obtain approximately 1 kg of material (Fig. 1). This protocol was implemented so that our sampling was as representative as possible of sediment sources of coastal escarpment at each site. We took care to avoid sampling sediment that might have been deep in colluvium before being excavated very recently. We give here two examples of this concern in sampling. (i) In Peru, we were careful to sample well above the road cut and its associated disturbance (Fig. S23). (ii) At the foot of coastal escarpments, we sampled colluvium high enough above the sea to avoid contamination by wave-borne material (Figs. S2–S24). To do this, we systematically sampled debris wedges that covered any slope break at the toe of the escarpment (usually located a few metres above the sea), so that the sediment sampled necessarily came from higher up (i.e. no contamination by sand brought in by waves). In most cases, the debris wedges are located at the outlets of shallow debris channels eroding the escarpment, which increases the likelihood that the sediment collected statistically come from the entire escarpment (Figs. S2–S24). We also collected sediment a few metres above sea level to avoid any contribution from pelagic sand.

The 0.5–1 mm sand fraction was then chemically prepared, following the protocol described in Zavala et al. (2021). The concentration in ¹⁰Be C was then measured at AMS ASTER (CEREGE, Aix-en-Provence, France).

To calculate the average ¹⁰Be production rate, we first determined a polygon bounded downstream by the sampling line and upstream by the cliff crest line, which also constitutes a limit for sediment sources. The average production rate (Table 1) was then calculated by averaging its value over all the pixels of the digital elevation models (DEMs; Fig. 2) contained within this polygon. The French airborne lidar-derived RGE ALTI 5 m DEM with 5 m resolution was used for France (https://geoservices.ign.fr/, last access: 21 January 2023) and the Shuttle Radar Topography Mission (SRTM) DEM (about 30 m of resolution) was used for Peru (https://www.earthdata.nasa.gov/, last access: 21 January 2023). For each pixel, the ¹⁰Be production rate is determined by

$$\begin{array}{}\text{(1)}& {\displaystyle }P=P_{\mathrm{sp}}+P_{\mathrm{sm}}+P_{\mathrm{fm}},\text{(2)}& {\displaystyle }{P}_{\mathrm{sp}}=P_{\mathrm{SLHL}}{f}_{\mathrm{sp}}{S}_{\mathrm{sp}},\text{(3)}& {\displaystyle }{P}_{\mathrm{sm}}=P_{\mathrm{SLHL}}{f}_{\mathrm{sm}}{S}_{\mathrm{sm}},\text{(4)}& {\displaystyle }{P}_{\mathrm{fm}}=P_{\mathrm{SLHL}}{f}_{\mathrm{fm}}{S}_{\mathrm{fm}},\end{array}$$

where P_sp, P_sm, and P_fm are the ¹⁰Be production rates of a given cosmogenic nuclide (CN) at the Earth's surface by spallation (subscript “sp”), slow muon capture (subscript “sm”), and fast muon interactions (subscript “fm”), respectively (Braucher et al., 2003). P_SLHL is the total sea level/high-latitude production rate of the considered nuclide (P_SLHL= 4 atoms g⁻¹ yr⁻¹). f_sp, f_sm, and f_fm are the fractions of this production rate due to spallation, slow muon capture and fast muon interactions occur (f_sp= 0.9886, f_sm= 0.0027, and f_fm= 0.0087). S_sp, S_sm, and S_fm are the respective scaling factors, depending on latitude and elevation based on Stone (2000).

P can be decreased by topographic shielding and increased by the shorter cosmic ray distance between the scarp surface and any given point inside the rock on steep slopes (DiBiase, 2018). In order to evaluate the factor that should multiply P, we used the MATLAB code provided by DiBiase (2018), and we computed this factor for different scarp slopes between 20–85°, corresponding to the range of our sites. We found that the factor increases slightly from 1.0006 for 20° to 1.107 for 60° and then sharply from 1.17 for 65° to 3.1 for 85°. The factor is always greater than 1 because the effect of the shorter cosmic ray paths dominates over the topographic shielding, as explained by DiBiase (2018). We corrected all the ¹⁰Be production rates calculated in the following by this factor, which remain on the order of 1.1, as the slopes are smaller than 60°.

The mean erosion rate ϵ was calculated assuming a steady state and neglecting the radioactive decay as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathit{ϵ}={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }C}}(P_{\mathrm{sp}}{\mathrm{\Lambda }}_{\mathrm{sp}}+P_{\mathrm{sm}}{\mathrm{\Lambda }}_{\mathrm{sm}}+P_{\mathrm{fm}}{\mathrm{\Lambda }}_{\mathrm{fm}}),\end{array}$$

where C is the sample ¹⁰Be concentration; the rock density is ρ=2.6 g cm⁻³; and where Λ_sp=160 g cm⁻², Λ_sm=1500 g cm⁻², and Λ_fm= 4320 g cm⁻² (Braucher et al., 2003). Neglecting radioactive decay is justified here because the estimated integration times are much shorter than the ¹⁰Be half-life.

The erosion rate uncertainty corresponds to the propagated analytical uncertainty and the 15 % uncertainty in P. Using the neutron attenuation length $L={\mathrm{\Lambda }}_{\mathrm{sp}}/\mathit{\rho }$ of 0.6 m, we calculated the integration time $\mathit{\tau }=\frac{L}{\mathit{ϵ}}$ for each sampled cliff (Table 1).

For catchments with a strong 3D curvature including an escarpment, Wang and Willett (2021) proposed a method for quantifying the horizontal retreat of the escarpment by considering that the flux of eroded material occurred through a vertical surface. In our case, the escarpments have an almost constant slope, α. Thus, we simply turn ϵ into a horizontal cliff retreat rate, i.e. $r=\frac{\mathit{ϵ}}{\tan \left(\mathit{\alpha }\right)}$. In the following examples, r is similar to ϵ because the slopes of the cliffs are all close to α= 45° (Table 1).

Table 1 Features and results for each sample. The mean slope, α (°), of the cliff was measured in the field. The mean production rate, P (atoms g⁻¹ yr⁻¹), of cosmogenic ¹⁰Be was calculated for the (narrow) catchment upstream each sampling site using the DEMs of each location; it includes a correction for the slope, following DiBiase (2018) (see the text). The vertical erosion rate, ϵ (mm yr⁻¹), was calculated from the concentration, C, in ¹⁰Be with a neutron attenuation length, L, of 0.6 m and assuming 15 % of uncertainty in P. The horizontal retreat rate, r, (mm yr⁻¹) was calculated from α, ϵ, and propagating the uncertainty in ϵ. The integration time, τ (years), was calculated from L and ϵ (see text and Supplement for more details).

We tested the consistency of our results by comparing our results with the geomorphic context. Three study areas were sampled, namely the Côte Vermeille in southern France, close to the Spanish border, to the north of Banyuls-sur-Mer (VERM samples; Fig. 2a); the Costa Brava coast around Sant Feliu de Guixols and Tossa de Mar in Spain (BRAV samples; Fig. 2b); and around Atico in the Peru coastal desert (COSTA samples; Fig. 2c).

The cliffs of the VERM series are 15 to 40 m high and consist of pelites and sandstone pelites (Argelès-sur-Mer geologic map, 1:50 000, French Bureau de Recherches Géologiques et Minières, BRGM; Llançà geologic map, 1:50 000, Institut Cartogràfic i Geològic de Catalunya). The colluvium of the VERM1 and VERM2 sites is relatively fresh, while that of the VERM3 site appears older and is partly stabilized by herbaceous plants (it can be estimated that this vegetation needs a decade to become established).

The cliffs of the BRAV series are 50 to 120 m high and consist of leucogranites (BRAV1, BRAV2, BRAV3, and BRAV4), granodiorites, and alkaline granites (BRAV5; Baix Empordà and Selva geologic maps, 1:50 000, Institut Cartogràfic i Geològic de Catalunya). The colluvial wedges of the BRAV1, BRAV2, BRAV3, and BRAV4 sites are relatively fresh, while those of the BRAV5 site are largely stabilized by abundant vegetation consisting of trees. BRAV2 is the only sample for which part of the sampled sand could come from above the coastal escarpment. We discuss this further below. In this area, the sea level rose rapidly just after the Last Glacial Maximum, and then at a decreasing rate, to reach a more or less constant rise by around 0.4 mm yr⁻¹ over the last 6 millennia (Vacchi et al., 2021).

The sample sites for the COSTA series were selected to test the effect of a drier climate. Furthermore, some parts of the coast have been protected from the action of waves by the uplift of a shore platform related to the subduction of the Nazca plate beneath South America at a rate of about 0.45 mm yr⁻¹ over the Pleistocene (see Regard et al., 2021, 2010; Melnick, 2016; Malatesta et al., 2022; Saillard et al., 2017). The emergence of such a platform at the base of the cliff stops direct wave action at the bottom of the cliff (emergence since 200 ka after cliff foot elevation and a constant uplift rate of 0.45 mm yr⁻¹). Erosion rates measured on these cliffs should therefore be lower than those measured on active cliffs. COSTA1 was sampled at the base of an active cliff, whereas COSTA2 and COSTA3 were sampled over an uplifted shore platform spanning approximately 300 m and 1 km wide, respectively (Fig. 2c). The cliffs are between 200 and 300 m high and are made of intrusive rocks (coastal batholith). The arid climate prevents the development of vegetation on these colluvial wedges.

3 Results

Table 1 and Fig. 3 show the obtained erosion rates, similar to the retreat rates, and their respective integration times. For the VERM and COSTA series, five samples yield similar low-erosion rates between 0.3 and 0.6 mm yr⁻¹. Two other samples (VERM3 and BRAV5), gathered from vegetated colluvium, have possibly lower rates between 0.05 and 0.1 mm yr⁻¹ that are within 1σ uncertainties in the faster-eroding cliffs. An erosion rate similar to VERM and COSTA active sites is obtained for the COSTA1 sample in arid Peru, subject to wave erosion. The COSTA2 and COSTA3 samples, gathered on a palaeo-cliff protected from wave action, are eroding ∼ 20 times more slowly (Table 1; Fig. 3). The integration times range from 1.1 to 33.6 kyr (Table 1l Fig. 3), i.e. from late Holocene to pre-LGM (Last Glacial Maximum) periods.

Figure 3 Results from the cosmogenic ¹⁰Be abundance measurement in colluviums. Red: VERM series. Green: BRAV series. Blue: COSTA series. Circles: vertical erosion rates (ϵ) calculated from the C concentration of ¹⁰Be. Squares: horizontal retreat rates (r) deduced from the mean cliff slope. Columns: integration time (τ) of the measurement. The colluviums at the VERM3 and BRAV5 sites are vegetated, and the BRAV2 site is a cape, which is consistent with the lower-erosion rates. Due to uplifted platforms at the base of the cliffs, coastal erosion no longer occurs on COSTA2 and COSTA3.

4 Discussion

The obtained erosion rates are remarkably consistent with their geomorphic setting. For the BRAV and VERM samples, six samples that could be considered part of a repeatability test show similar erosion rates but with large uncertainties. The possibly lower-erosion rates of VERM3 and BRAV5 are consistent with the vegetated nature of colluvial wedges either because the vegetation protects the coast from erosion or because lower erosion has allowed vegetation to develop.

The COSTA1 sample in Peru corresponds to an erosion rate similar to the values for the Mediterranean coasts, despite a more arid climate. It is beyond the scope of this paper to discuss the factors that control the coastal erosion in further depth, especially since it would require an analysis of the wave amplitude and frequency distribution. Nevertheless, the similar coastal erosion obtained in Peru and for the Mediterranean coast suggests that different combinations of climates, wave frequency distributions, and vertical uplift rates may result in similar erosion rates.

Furthermore, the Peruvian site COSTA1, where waves attack the base of the cliff, yields, as expected, a 20 times higher erosion rate than the other sites in front of uplifted marine terraces which are protected from wave action. One obvious limitation of our study is the lack of an alternative method for comparing the erosion values we obtain averaged over several millennia, which is often the case with any new method. We note, however, that the erosion rates of the two sites in Peru preserved from wave action (0.018–0.025 mm yr⁻¹) have values similar to the average erosion rates of the catchments draining the Andes near these sites (0.015–0.02 mm yr⁻¹; Starke et al., 2020). In the Pyrenees, the minimum erosion rates of the coastal escarpment (VERM3 and BRAV5) are similar to millennial erosion rates in the Pyrenees nearby (0.06–0.14 mm yr⁻¹; Molliex et al., 2016), as well as other small French catchments draining towards the Mediterranean in the Southern Alps (Mariotti et al., 2021) and Corsica (Molliex et al., 2017) (0.01–0.24 mm yr⁻¹). However, these comparisons should be treated with caution as the erosion processes are different. Overall, the rates of cliff retreat are around 1 to 20 times higher than the average erosion rates in the nearby catchment areas.

As we have obtained coastal erosion values that are slower than those documented worldwide (Prémaillon et al., 2018), in the following we discuss several biases that could lead to an underestimation of retreat rates in our study.

One bias could be the contribution of grains eroded from the land surface above the rocky coast. Although we paid attention to selected sites between rivers that ensure a negligible contribution of the land surface above the cliff, half of the surface of the catchment draining toward BRAV2 is a gentler slope upstream of the cliff itself. If the upstream zone erodes much more slowly than the coastal escarpment, it is possible that this contribution partly explains the slightly lower-erosion rate of BRAV2. Nevertheless, even an extremely low-erosion rate of the land surface above the coast could not artificially decrease the retreat rate by a factor larger than 1.5 because the contribution of this area to the flux of quartz grains would be negligible in that case (see the Supplement).

It has been shown that outcropping bedrock can erode at a smaller rate than surrounding loose material in a catchment (e.g. Bierman and Caffee, 2001; Heimsath et al., 2006; Lodes et al., 2023). We may wonder if this difference can affect the estimation of the mean scarp erosion rate from the ¹⁰Be concentration in colluvium. This is unlikely because the coastal scarp retreats horizontally, so that outcropping bedrock contributes as much as the other sources that provide the sampled colluvium over millennia. However, if a slower-eroding part of the bedrock was not sampled at all by our approach, then our erosion rate would be overestimated and not underestimated. For the Peruvian samples, the coastal scarp is almost entirely covered by a regolith so that the distinction between bedrock and loose sediment does not hold.

Another source of bias could be the shallow landslides that may feed the colluvium with grains that have lower-¹⁰Be concentrations. However, in that case, the bias would increase the apparent erosion rate. Furthermore, the similar erosion rates obtained for five samples in the VERM and BRAV series in the same geomorphic context indicate that these stochastic processes have a negligible effect.

Another bias could arise from the delayed adjustment of the ¹⁰Be concentration in response to an increase in the erosion rate. The coastal erosion rate probably increased once the current sea level had been established 6 kyr ago on average (e.g. Lambeck, 1997; Bintanja and van de Wal, 2008; García-Artola et al., 2018); this timescale is shorter than the integration time of the samples VERM3 and BRAV5. The cosmogenic signal adapts to a changing erosion rate with a delay. We thus wonder if the ¹⁰Be signal could be inherited from a former low-erosion rate period, leading to estimations of erosion rates smaller than 1 or 2 orders of magnitude. In order to quantify this bias, we carried out end-member simulations for which the erosion rate is constant (0.05 or 0.5 mm yr⁻¹) during 100 kyr and then multiplied by 10, 100, or 200 in the last 6 kyr. In the worst scenario, the erosion rate averaged over the last 6 kyr is underestimated by 80 % (see Fig. S1 in the Supplement). Therefore, this bias can reduce the real value by half, but it cannot change the order of magnitude of these two erosion rates. In other words, if the retreat rate had been 1 or 2 orders of magnitude faster in the last 6 kyr, the ¹⁰Be concentration for VERM3 and BRAV5 should have recorded it. For the other samples subject to erosion by sea waves, the integration time is shorter than 6 kyr, although the erosion rates are still low (< 1 mm yr⁻¹), which gives confidence in these low values.

5 Conclusions

Our approach provides a new method to quantify coastal erosion rates less than 1 cm yr⁻¹ over millennia. These rates are typically averaged over integration periods of millennia, with some of them (the highest integration times) encompassing the current highstand and the period beforehand when the sea level was lower and the waves did not reach the cliff foot. The bias on the coastal erosion rate associated with this variable erosion should not exceed −80 %, thus giving a valuable order of magnitude for coastal retreat rates. As there is no limitation to reproduce this approach where colluvium is present, we anticipate that it will fill a significant gap in the rocky coast retreat rate database (retreat rates ranging 0.05–5 mm yr⁻¹), improve our understanding of controlling factors, and provide a temporal benchmark to evaluate current and future rocky coast erosion hazards.