Cliff recession geodynamics variability and constraints within poorly consolidated landslide-prone coasts in the southern Baltic Sea, Poland

Abstract

This study identifies the reasons for geodynamics variability of the coastal system within two cliff-shore sections of the southern Baltic Sea (SBS). The comparative analysis included distinct moraines and their foregrounds near the open sea (S1) and within the Gulf of Gdańsk (S2). Short-term trends indicate a direct link between landslide occurrence and increased cliff retreat. Long-term (total) values were obtained by developing the 4F MODEL for large-scale applications, based on the analysis of remote sensing and hydroacoustic data (to determine the extent of shore platforms), the modelling of higher-order polynomial functions describing their extent, followed by the integral calculus of the indicated functions within the open-source Desmos environment. The retreat dynamics for individual landslides (S1) was an order of magnitude higher (m/yr) than the average for the whole cliff section (0.17 ± 0.008 m/yr), which correlates well with medium- and long-term development tendencies and recession dynamics, revealed by the numerical modelling method, since approximately 8 ka b2k, years before 2000 CE (at S1 = 0.17 ± 0.020 m/yr, at S2 = 0.11 ± 0.005 m/yr). While the approach described in this paper can reveal, project, and simulate the dynamics of past and future trends within other cliffed coasts shaped in tideless conditions, it also proves stable moraine erosional responses to sea-level rise since the Mid-Holocene.

Introduction

The monitoring of coasts consisting mainly of poorly consolidated (soft rock) Pleistocene sediments is common in the Baltic Sea countries along the coasts of Germany (Kolander et al, 2013; Kuhn and Prüfer, 2014), Poland (Kostrzewski et al., 2015; Łabuz, 2017; Frydel et al., 2017; Winowski et al., 2022), Russia (Spiridonov et al., 2011), Lithuania (Bitinas et al., 2005; Bagdanavičiūtė et al., 2012), Latvia (Eberhards and Saltupe, 2006), and Estonia (Orviku et al., 2013). Identification of geohazards and their impact on coastal development is undertaken worldwide in the North Sea basin countries, that is, in the United Kingdom (Walkden and Hall, 2005; Poulton et al., 2006; Hurst et al., 2016), the Netherlands (van Rijn, 2011); in France (Medjkane et al., 2018); around the Mediterranean Sea (Montoya-Montes et al., 2012); in Portugal (Terefenko et al., 2018); in New Zealand (de Lange and Moon, 2005; Dickson and Perry, 2016); in China (Liu et al., 2011); and in the United States (Collins and Sitar, 2008; Young, 2018).

Local studies of coastal sea cliffs (Rosser et al., 2005; Hapke and Plant, 2010; Hackney et al., 2013; Earlie et al., 2015; Frydel et al., 2017) show that the pace of cliff retreat can be derived from light detection and ranging (LIDAR) observation data, which enable the mapping of coastal topography and designating areas susceptible to landslide (Poulton et al., 2006; Schulz, 2007; Baldo et al., 2009; Ventura et al., 2011; Dickson and Perry, 2016). However, a wide array of processes need analyses that would enable an assessment of the evolutionary tendencies of a coastal system and an understanding of its behaviour (Ashton et al., 2011). On a long-term scale, the model aimed at determining the response of soft rock coasts (Carpenter et al., 2014) to sea-level fluctuations and waves indicates that the sea-level rise will accelerate the recession of cliffs to the greatest extent, while the increase in storm frequency (by moderately large rather than extremely large waves) may account for further acceleration in retreat rate (Trenhaile, 2010). Yet in the southern Baltic Sea (SBS), which is considered virtually nontidal (Ostrowski and Stella, 2016), short-term morphogenetic factors (sea level and sea state, wind speed and direction, rainfall intensity, groundwater level) and the degree of human influence shape the coast most noticeably, mainly through the impact of extreme meteorological phenomena (Zawadzka-Kahlau, 1999; Bitinas et al., 2005; Kostrzewski et al., 2015; Zhang et al., 2011).

Research objectives

Research objectives included recognition of the coastal systems of two seaside moraines subject to mass wasting processes and investigation of the influence of lithology and hydrometeorological factors on landslide development. Multitemporal models were required to illustrate cliff deformation and erosion on a short timescale (decadal) and to compare hydrodynamic regimes and their influence on cliff erosion dynamics in the open sea (Swarzewo Moraine) and the Gulf of Gdańsk (Redłowo Moraine). A further aim was the identification of areas most susceptible to mass wasting and the determination of their development tendencies as an element of spatial planning available for application by maritime authorities to enhance human and infrastructure safety. An additional task comprised resolving uncertainties regarding discontinuities of chronostratigraphic units within the Redłowo Moraine (Kaulbarsz, 2005; Woźniak et al., 2018). Further palaeogeographic analysis was targeted at determining the extent of the palaeomoraines subject to erosion from the Mid-Holocene onwards. Yet the main objective facilitated introducing a method capable of revealing total cliff erosion dynamics within moraines containing boulder-rich glacial tills of the Pleistocene age in their profile, using numerical modelling incorporating the 4F MODEL (Frydel, 2022), a step forward to the quantification of geologic processes (coastal erosion and accumulation, marine transgression and regression, ice-sheet advance and decline) in an open-source spatiotemporal environment.

Regional setting

Due to the fractal-like geometry of the coastline (Mandelbrot, 1983) and its ephemeral course, the length of the Polish Baltic coast equals 498 km at a map scale of 1:200,000 (Uścinowicz et al., 2004). Accordingly, lithology and topographic features allow distinguishing barriers, wetlands, and cliffs that extend across approximately 86 km of the Polish coast and include approximately 62 km of active cliffs, 20 km of inactive (dead) cliffs no longer subject to wave action with limited mass wasting within the slopes, and approximately 4 km of cliffs encompassed by coastal engineering structures.

Origin of the Polish Baltic Sea coast

Within the Polish SBS coast, a pattern of valleys, plains, and shoals built of Holocene sands, silts, and biogenic sediments of marine, fluvial-lake, and glaciofluvial origin intertwine with common coastal moraines (morainic uplands) built from typical Pleistocene sediments, especially boulder-containing tills, sands, and gravels of glacial and glaciofluvial origin, as well as locally occurring Pleistocene and Neogene clays and silts of ice-marginal origin (Tomczak, 1995; Uścinowicz et al., 2004). Quaternary sediments were accumulated over approximately 1 Ma during cyclic periods of glaciation (stadials), when the Scandinavian Ice Sheet (SIS) moved southwards, causing the deepening of the Baltic basin and glacial erosion of bedrock and sediments from older glaciations, as well as depositing sediments carried by the ice sheet (Uścinowicz, 2011). During climate warming (interstadials), the ice sheet underwent melting and recession towards the north, when the accumulation of glaciofluvial sediments and the deposition of fine-grained sediments under reservoir conditions proceeded (Uścinowicz, 2011). The decline of the last SIS during the Gardno Phase (16.5 ka ago) resulted in the gradual filling of the reservoir with water from a melting glacier front (Tylmann and Uścinowicz, 2022). As a result of the interaction of isostatic factors and eustatic fluctuations of the ocean level, the level, extent, and salinity of the reservoir were variable in time and space, with periods of transgression and regression (Uścinowicz, 2003, 2006; Harff et al., 2017) characterised by varying dynamics of shoreline migration (Frydel, 2022; Tylmann and Uścinowicz, 2022).

Locations of the study sites

The study examines two seaside moraines (morainic uplands) and their marginal zones (cliffs) located in NE Poland (Fig. 1) near the open SBS (Site 1 [S1]) and near the semi-enclosed SBS, in the Gulf of Gdańsk (GOG) (Site 2 [S2]), subject to distinct hydrodynamics and diverse slope exposures, but representing fairly common lithology. S1 encompasses the Jastrzębia Góra city area covering the northwestern part of the Swarzewo Moraine with Jastrzębia Góra Cliff (JGC) (Supplementary Fig. 9); and the eastern part of the Karwia Lowland. Whereas S2 is located near Gdynia city, which includes the moraine marginal zone of the Gdynia Orłowo Cliff (GOC) in the eastern part of the Redłowo Moraine (Supplementary Fig. 10), with its foreground subject to presumably calmer wind–wave environment when compared with S1. GOC slope exposure ranging from northwestern to southeastern, in comparison to the latitudinal course of the JGC, is noteworthy. Anthropogenic impacts vary between the sites, with S1 subject to hydroengineering protection since the 1990s and S2 remaining virtually unprotected (Supplementary Text 1).

Figure 1. (A) Main types of Polish coasts with location of first (S1) and second (S2) research zones, modified from (Uścinowicz et al., 2004), and regional setting of the study areas: (B) S1, Jastrzębia Góra region; and (C) S2, Gdynia Orłowo region

Geology, geomorphology, and geodynamics

Within the Polish coast, soft Quaternary rock cliffs predominate, occasionally underlain by Neogene or Paleogene strata (Tomczak, 1995; Fig. 2). In littoral zones of the inspected regions, the thickness of the Quaternary-age deposits falls below 10 m and is amongst the lowest in the Polish coastal zone (Uścinowicz, 1995a), where more resistant carbonaceous clays (Neogene) occasionally appear (Kramarska et al., 2014). Considerable shallow-water areas lack sand coverage (up to 16 m below sea level [m bsl]), indicating storm wave influence resulting in erosion of the sea bottom, followed by the sediment transport predominantly to the east (Uścinowicz, 1995b).

Figure 2. Geologic cross section along the coast from Jastrzębia Góra (W) to Gdynia (SE); H, Holocene: sands and biogenic deposits; P, Pleistocene: tills, sands, silts, and clays; N, Neogene: sands, brown coal, silts, clays, and siltstones; Pg, Paleogene: siltstones, claystones, silts, sands. Thick black line stands for coastal engineering structures; notice the relationship between lithology and presence of coastal protection systems in places where the Holocene strata occur and within the edges of Pleistocene moraines. Modified from Tomczak (1995); Maritime Office in Gdynia (https://www.umgdy.gov.pl/mapy); Google Earth. ALS, airborne laser scanning; TLS, terrestrial laser scanning; SBES, single-beam echosounder; MBES, multibeam echosounder; SSS, side scan sonar.

The Swarzewo Moraine (S1) consists of up to three layers of Pleistocene tills from subsequent glaciations, fluvioglacial sands, and gravels overlaying the Miocene strata comprising silts and clays (Jurys et al., 2006) deposited in the glaciolacustrine environment (Olszak et al., 2008), which corresponds with the Grudziądz/Sassnitz interstadial (Marine Isotope Stage 3) of the Vistulian/Würm glaciation (Börner, 2010). By comparison, laminated silts and clays interbedded with fine and silty sands are floodplain-zone deposits (Sokołowski et al., 2019) favouring the development of mass wasting processes due to complex hydrogeologic conditions and glaciotectonic deformations (Lidzbarski and Tarnawska, 2015). However, the presence of glaciofluvial sand layers between clayey strata influences the groundwater drainage towards the sea. Yet landslides in the Polish coastal zone are not limited to land-based slumping but also include landslides of a subaqueous origin (Rudowski et al., 2016).

In the case of the Redłowo Moraine (S2) margin, three cliff sections were distinguished, that is, the northern, southern, and central parts (Orłowo Headland), which consist of tills (Kaulbarsz, 2005). As a result, rockfalls occur within the promontory, while glaciotectonically deformed silts and clays allow landslides and slumps to develop in the northern and southern parts of the GOC. Such clay clasts are characteristic of the fluvial system of meandering river deposits in northern Poland (Sokołowski et al., 2019). The Pleistocene tills, which contain a substantial number of boulders, are separated by a layer of sand and gravel. Near 81.3 km chainage, a sudden lack of continuity within two chronostratigraphic units O-1 and O-2 was mapped (Woźniak et al., 2018), also recognized as an overthrust (Kaulbarsz, 2005; Jurys et al., 2006), while another view on its origin is discussed here.

Numerous studies on cliff retreat in the Polish coastal zone have been undertaken based on transverse profiles (e.g., Subotowicz, 1982; Kwoczek, 2007; Florek et al., 2008, Łęczyński and Kubowicz-Grajewska, 2013; Kostrzewski et al., 2015). According to geodetic and photogrammetric measurements of the GOC, the recession pace equalled 0.45 m/yr in the 1970s, later increased to 0.9 m/yr in 1977–1990, to reach 1.2 m/yr in 1992–1997 and 0.78 m/yr in 1997–2007.

Environmental factors

In the Baltic Sea area, extreme short-term weather events such as storm winds and storm surges associated with cyclonic circulation (Zeidler, 1995) due to the North Atlantic Oscillation (NAO) may pose a high threat to coastal zones (Kostrzewski et al., 2015), while the impact of tides with an amplitude of several tens of millimetres (Müller-Navarra, 2003) is negligible. Conditions likely to generate storm winds over the Baltic Sea are observed mostly from October to January (Cyberski, 2011), though multiannual storm frequency and severity are irregular, but the SBS lacks a long-term pattern in terms of multidecadal variability and strength (Bärring and Fortuniak, 2009). However, intense storm winds frequently blow from the west or southwest (Kwiecień, 1990; Cyberski, 2011), demonstrating a link with positive NAO phases (Formela and Marsz, 2011) in northern Europe, while in negative NAO phases, westerly winds shift southwards into the western Mediterranean (Hurrell, 1995) leading to periods of increased storminess detectable at a millennial timescale (Sabatier et al., 2012). The foundation for environmental factor analysis was provided in Supplementary Text 2.

Materials and Methods

A range of techniques and datasets were employed to inspect impacts on the S1 and S2 coastal systems and analyse the behaviour of cliffs. To some extent, these issues have already been described for S1 by Kramarska et al. (2011), Uścinowicz et al. (2014), Karwacki (2021), and Kamiński et al. (2012, 2023); and for S2 by Subotowicz (1982), Kwoczek (2007), Kubowicz-Grajewska (2016), and Małka et al. (2017). New elements incorporated into this study encompass geohazard investigations with slope deformation analyses, consideration of environmental forcing factors, definition of multiscale erosion dynamics, and investigation of palaeogeographic aspects of moraine setup and erosion dynamics since the Mid-Holocene.

The landslide boundaries were confirmed during the fieldwork, following the national Landslide Counteracting System (SOPO) (Grabowski et al., 2008). Variability in cliff erosion dynamics was assessed across seasonal, annual, decadal, centennial, and millennial (since the Mid-Holocene) timescales. The methodology included qualitative and quantitative analysis of cartographic, remote sensing, hydroacoustic, and topographic LIDAR data, along with hydrometeorological models. Landslides and cliffs were profiled; deformations of slopes and hydro-technical structures were detected using GIS tools as well as 3D and 4D modelling (multitemporal digital terrain models [DTMs]).

Short-term erosion—seasonal, annual, and decadal—based on terrestrial laser scanning (TLS), and airborne laser scanning (ALS)

Overall, 11 TLS campaigns were undertaken between April 2010 and November 2017, alongside ALS campaigns from September 2008 to September 2017, to account for slope deformation and seasonal and annual rates of cliff erosion. TLS campaigns employed a Riegl VZ-400 laser scanner, a Trimble R8 GNSS, a controller (TSC2), a calibrated camera (Nikon D700), and real-time network point cloud positioning using the ASG-EUPOS reference station network. These tools were used to acquire consecutive scans (beam width: 0.04°–0.08°). As a result, each point cloud centre was specified with a precision of 10–30 mm and subjected to filtering before obtaining DTMs. Data capture, registration, filtration, and resampling followed the standard procedures (Poulton et al., 2006; Schulz, 2007; Baldo et al., 2009; Kramarska et al., 2011; Frydel et al., 2017) along with an examination of GNSS data integrity followed by error minimisation using the Multi Station Adjustment algorithm (introduced by Riegl in the RiScan PRO software). Afterwards, only single- and last-echo pulses were kept, which allowed the reduction of point cloud by 15–20%, and an octree filter (3D analogue of quadtree; Finkel and Bentley, 1974) was applied to decimate point cloud as much as 95%. Finally, DTMs were produced. A visualisation of the simplified TLS LIDAR data-handling scheme resulting in a 3D description of the coastal moraine/cliff can be reached via https://youtu.be/Mt1rtdDKfbQ (last accessed on September 12, 2023). Slope deformations at S1 were displayed using ALS data and at S2 using TLS data acquired by the Polish Geological Institute (PGI-NRI). A description of the procedures, including digital terrain modelling, volumetric calculations, and postprocessing, is provided in Supplementary Text 3. The laser scanning survey dates (TLS and ALS) for study areas S1 and S2 are described in Supplementary Table 1.

Medium-term erosion—decadal and centennial—based on maps

Morphodynamic analyses utilised rectified German historical maps (scale: 1:25,000; root mean square [RMS] = 15 m, Messtischblatter dated 1875 and 1908 for S1 and S2 areas, respectively) to examine centennial shore dynamics, while topographic maps (scale: 1:10,000; RMS = 8 m, dated 2000) were used to inspect cliff retreat rates during the first decade of the twenty-first century, following standard procedures. Large-scale maps (1:10,000) were presented in the PL-2000 coordinate system (EPSG2177). For smaller scales (1:25,000), the PL-1992 coordinate system (EPSG2180) was used. To bring remote sensing and hydrometeorological data into a common Cartesian coordinate system, altitudes determined based on GNSS measurements (PL-KRON86-NH) referenced to the height of the Pomeranian quasi-geoid (at 5 × 5) above the reference GRS80 ellipsoid.

Long-term erosion—millennial—based on numerical modelling

Palaeogeographic analyses utilised remote sensing (LIDAR) and hydroacoustic data (sound navigation and ranging [SONAR]; multibeam [MBES] and single-beam echosounder [SBES]). The total recession coefficients (C _{mid−H→b2k}) were revealed via numerical modelling of polynomial and (occasionally) linear functions in the Desmos environment (https://www.desmos.com/calculator?lang=en, last accessed on September 10, 2023) and validated in the GIS environment (https://www.esri.pl/produkty/arcmap, last accessed on September 10, 2023), in accordance with Frydel (2022). For coordinate system conversions and production of *.kmz maps, the Global Mapper suite was used (https://www.bluemarblegeo.com/global-mapper, last accessed on September 10, 2023). Modelling utilised hybrid DTMs based on hydroacoustics, TLS, and ALS acquired from the Maritime Office in Gdynia (MOiG) (https://www.umgdy.gov.pl/mapy, last accessed on September 12, 2023) and the Central Geodetic Documentation Centre (CODGiK), originating from the nationwide project Informatyczny System Osłony Kraju (ISOK). The resolution of respective data equalled 9 pts/m² (MOiG) and 4 pts/m² (CODGiK). The area between cliff edges and the shoreline was rendered using TLS data (100 pts/m², local maximum up to 400 pts/m²). The littoral zone at S1 was described using data acquired from the Maritime Institute in Gdańsk for the 4D Cartography Pilot Project (Kramarska et al., 2014) and the 4D Cartography Project (Uścinowicz et al., 2018), while at S2 it was described based on data provided by Pakszys (2014). The range of data input and an overview of the method are provided in Supplementary Text 4.

Onset time of erosion of SBS coastal moraines

According to Uścinowicz (2003, 2006), the relative sea-level (RSL) curve indicates that SBS shoreline erosion began in the initial phase of the Atlantic period after the tectonic uplift of terrain ceased and the moraines fully emerged, owing to the glacio-isostatic adjustment (Harff et al., 2017), a specifically glacio-isostatic rebound that followed the deglaciation (Uścinowicz, 2003). Accordingly, the sea-level curve of approximately 8.0–7.5 ka before the year 2000 (b2k = −50 yr BP) settled around −15.0 m bsl. This time frame corresponds to the beginning of cliff retreat in Germany (Hoffman and Lampe, 2007) and New Zealand (de Lange and Moon, 2005) also adopted by Meyer and Harff (2005) for modelling of coastline changes by comparing DTMs and RSL curves. Differences between the sea-level rise curves are evident between the sites due to distinct location, origin, and lithology of the coasts and glacio-isostatic and neotectonic factors, although they clearly show a steep increase in water level during the early Littorina stage of the Baltic Sea (Müller, 2001), which links with increased coastal erosion. Notably, data from the Little Belt region (Bennike and Jensen, 2011) and the Wismar Bay region (Lampe et al., 2011) indicate that the sea level in the central Holocene (approximately 8 ka b2k) was similar to that of northern Poland (Uścinowicz, 2003, 2006; Jegliński, 2013), that is, 12 or so metres lower than today.

The 4F MODEL adaptation to obtain total retreat dynamics coefficients $( \C _{{\rm mid}-{\rm H}\to {\rm b}2{\rm k}})$

Moraines found in northern Poland commonly consist of glacial deposits that contain boulders, which are visible on the cliff face and below the foot of the cliff. Additionally, the rocky residuum visible in the littoral zone is a result of gradual marine erosion of the same Pleistocene tills (Subotowicz, 1982). Thus, it is expected that the range of boulders in the coastal zone determines the area of formerly existing moraines ($A_{t_{mid-H}\to t_{b2k}}$), while the time when erosion began (t _mid−H = approximately 8 ka b2k) is estimated, by relating bathymetric data from hybrid DTMs to the RSL curves for the SBS at the maximum seaward extent of boulder residuum for each shore platform.

Such recognition, following Occam's razor, allows determining the cliff's total retreat dynamics, that is, the coefficient of moraine recession (C _{mid−H→b2k}) since its initial phase in the Atlantic period, which can be achieved by adjustment of the formula by Frydel (2022), to the form introduced herein:

$$C_{mid-H\to b2k} = A_{mid-H\to b2k}/tB_t = \displaystyle{{\mathop \smallint \nolimits_{x_l}^{x_u} ( \mathop \sum \nolimits_{i = 0}^n a_ix^i) dx-\mathop \smallint \nolimits_{x_l}^{x_u} ( \mathop \sum \nolimits_{i = 0}^n b_ix^i) dx} \over {\vert {x_u-x_l{\rm \;}} \vert {\rm \ast }\vert {t_{mid-H}-t_{b2k}} \vert }}{\rm \;}( {{\rm m/yr}} ) \eqno ( Eq. 1) $$

C _{mid−H→b2k}

coefficient of moraine recession dynamics since the Mid-Holocene (cliffs total average retreat rate), expressed in metres per year (m/yr);

A _{mid−H→ b2k}

an area of an abrasion platform since the Mid-Holocene up to b2k (years before 2000, −50 yr BP), expressed in square metres (m²);

t _mid−H

beginning of moraine erosion, approximately 8 ka b2k, expressed in years (yr);

t _b2k

“present” time, expressed in years (yr);

t

interval between t _m−1 and t _m, expressed in years (yr);

B _t = |AB|

length, expressed in metres (m);

x _u, x _l

upper and lower integral limits, expressed in metres (m);

The projected average range of cliff extent subject to erosion $( {F_{t_0\to t_1}} )$ in a definite number of years (p) equals (adjusted from Frydel, 2022):

$$F_{t_{m-1}\to t_m} = pC_{t_{m-1}\to t_m}, \;{\rm for}\;p\in {\rm N\;\;\;}( {\rm m} ) \eqno ( Eq. 2) $$

The key 4F numerical modelling principles for a typical SBS cliff coast area are illustrated in Figure 3, where the division of A _{mid−H→b2k}, which stands for the difference between two Riemann integrals$\;[ \mathop \smallint \limits_A^B f_1( x ) dx]{, \;\;\mathop \smallint \limits_A^B f_2( x ) dx} ]$, by time span t, multiplied by the integral limit B _t = |AB|, designates the average recession dynamics C _{mid−H→b2k} for each moraine (cliff) since the Littorina transgression per shore metre. Within the functional programming framework in the Desmos environment, polynomial regression of successive functions based on nodal points, followed by integral calculus and application of the formula (Eq. 1) revealed the moraines’ total erosion dynamics.

Figure 3. Southern Baltic Sea (SBS) coastal zone diagram illustrating the adopted dynamics modelling principles, while the determination of the extent of the abrasion platform allows the application of a numerical modelling method based on functional programming and thus the determination of erosion dynamics since the Mid-Holocene (developed from Frydel, 2022).

SBS moraines have been subjected to abrasion since the Atlantic period, throughout the Subboreal to the present day (b2k, −50 yr BP) (Supplementary Fig. 4), when approximately 8 ka b2k the sea level began to rise from about 12 metres below the current mean sea level (Uścinowicz, 2003, 2006; Bennike and Jensen, 2011; Lampe et al., 2011; Jegliński 2013). During the initial phase of the Littorina transgression, the connection to the ocean through the Danish Straits resulted in a gradual change to a brackish environment (Rößler et al., 2011). The extent of the sea at that time increased in a southward direction, affecting the erosion intensification of the uplands commonly found on the southern coasts (remnants of the Vistulian/Würm glaciation). As a result, glacial outcrops—residua containing a significant number of boulders—formed on top of the developing shore platform while the sea level gradually increased. These boulder fields are visible on the seabed to a depth of about 12 m bsl, thus marking the maximum extent of eroded palaeomoraines. Based on the uniformitarian principle, it can be conjectured that the process of shoreline destruction, similar to what occurs at present, was accompanied by the development of mass movements, including landslides, causing a local intensification of the erosion process across the marginal slopes.

Environmental forcing factors in the research areas

Heterogeneous temporal attribution of some factors is related to the limited availability of hydrometeorological and geomorphological data. Therefore, an approach taking into account direct causal relationships between coast conditions and erosion determinants, based on the work of Young et al. (2021), could not be applied. Instead, qualitative analyses were performed based on available quantitative data. Erosion and landslide triggers were examined based on time-series data including rainfall (Prec) and snowfall (Sn _d) characteristics acquired from the Institute of Meteorology and Water Management. Velocity (V _c) and direction of wind-driven barotropic surface currents and sea-level (storm surge) data resulted from the M3D_UG model (Kowalewski, 1997) readouts (http://model.ocean.univ.gda.pl, last accessed on September 13, 2023). Warning (>5.50 m) and alarm sea levels (>5.70 m) were defined with tide-gauge data based on the Kronstadt reference system (PL-KRON86-NH). Events exceeding 6.0 m, that is, above 1.0 m above the mean sea level (msl = 5.0 m at the tide gauge), formed a separate category. Values of H _s derived from the WAM model (http://www.satbaltyk.pl, last accessed on September 13, 2023) were calculated using the energy spectrum for location in the closest possible vicinity to each site: 7 km north of Jastrzębia Góra (S1) and 4 km east of the Orłowo Headland (S2) (Fig. 1).

On-site environmental forcing factors based on time-series data, including significant wave heights (H _s) (Supplementary Fig. 1), time interval December 01, 2015 – January 03, 2019; sea levels, storm surge, and surface current velocities (V _c), time interval February 9, 2011–December 31, 2018 (Supplementary Fig. 2); precipitation (daily, Prec _d; monthly, Prec _m; annual, Prec _a), time interval January 1, 1951–October 31, 2018 (Supplementary Fig. 3A–E); and snow cover (Sn _d), time interval January 1, 1951–October 31, 2018 (Supplementary Fig. 3F–G), are described in Supplementary Text 2.

Results

The current study allowed for multiscale erosion dynamics hindcasts and projections for cliffed coasts in the southern Baltic Sea, Poland. Cliff margins of inspected moraines include landslide-prone areas (Fig. 4; Supplementary Fig. 11) with a total of 31 landslides and compound landslides (~4/km) covering 385,700 m² (Uścinowicz et al., 2018) identified within the Swarzewo Moraine and five (~5/km) in the southern part of the Redłowo Moraine (26,000 m²), where environmental factors affect the degree and dynamics of coastal transformation, landslide development and increased cliff retreat which mostly occur due to autumn rainfalls and winter storm surges (Fig. 5A, Supplementary Figs. 1–3). Both cliffs are characterized by a comparable average cliff crown elevation: JGC = 30 m above sea level (m asl), σ = 2.4 m; GOC = 32 m asl, σ = 10.8 m; however, the highest section of the GOC (S2) is almost twice as high as that of the JGC (S1). The trends of short-term erosion dynamics for inspected margins are fairly similar and amount to C _2008→2017(S1) = 0.16 ± 0.0088 m/yr and C _2014→2017(S2) = 0.23 ± 0.011m/yr. On the scale of the last century, the cliff development pace equalled C _1875→2015(S1) = 0.19  ± 0.02 m/yr and C _1908→2015(S2) = 0.25 ± 0.02 m/yr, which also corresponds to total recession dynamics for the subsequent moraines. Errors were determined as a percentage contribution; for the short-term scale (TLS/ALS) based on the resulting GNSS measurement accuracies and point cloud generalisation errors during DTM modelling; for the medium-term scale (maps) taking into account the RMS errors due to rectification; and in the long-term based on R² and the deviation from the middle value of the time interval of the designated onset time of upland erosion (7–9 ka b2k).

Figure 4. Northwestern margin of the Swarzewo Moraine, the Jastrzębia Góra Cliff subject to mass wasting. Hybrid digital terrain model (DTM) originating from terrestrial laser scanning (TLS)/airborne laser scanning (ALS).

Figure 5. (A) Geomorphometric profile (I-I’) through a slumped promontory near 134 km, profiles: April 2010, November 2010, March 2011 (Kramarska et al., 2011) and November 2012 (Uścinowicz et al., 2014); lithology based on Jurys et al. (2006) and Uścinowicz et al. (2014). Multitemporal digital terrain models (DTMs) showing biannual (B–E), annual (F), and cumulative (G) deformations of the landslide-prone coast to the east of massive cliff stabilisation system (MCSS), based on airborne laser scanning (ALS) data.

Site 1—Jastrzębia Góra region (S1)

During the last decade (2008–2017), the most severe mass wasting impacted the western part of the cliff (Fig. 4), especially the promontory at 134.0 km, which was destroyed due to landslide development, resulting in local crest shifting by nearly 30 m (Fig. 5A) and material loss of C _2010→2017 = 3.7 ± 0.18 m³/m/yr. The degree of slope deformation (qualitative) and the amount of erosion (quantitative) in two landslides, in the vicinity of the massive cliff stabilisation system (MCSS) (Fig. 5) and its hinterland (Fig. 6) for the years 2008–2010 and 2010–2012, are justified by the results of the qualitative hydrometeorological data analyses (Supplementary Figs. 2B, D, and E, 3B and C). Overall, the retreat coefficient for the JGC in the last decade, C _2008→2017(S1) = 0.17 ± 0.008 m/yr, is comparable with the centennial rate of C _1875→2015(S1) = 0.19 ± 0.020 m/yr, yet the erosion pattern along the coast shows high spatiotemporal variability (Fig. 7, Table 1).

Figure 6. (A) Longitudinal geomorphometric profile through the massive cliff stabilisation system (MCSS) at 134.27 km, based on terrestrial laser scanning (TLS); lithology modified from Witkowski and Wolski (2015). Multitemporal digital terrain models (DTMs) showing biannual (B–E), annual (F), and cumulative (G) MCSS deformations due to landsliding based on airborne laser scanning (ALS); colour intensity indicates the size of subsidence (red) or uplift/seaward push (green)

Figure 7. (A) Land use in the zone particularly exposed to landslides, along with projections (F) of the cliff retreat based on the extrapolation of retreat coefficients (C): in 50 yr (average), in 10 yr (maximum). (B) Morphodynamics of the Jastrzębia Góra Cliff (JGC; S1) in 2008–2017, groups of peaks correspond to the most active landslides, left-axis: retreat rate (m/yr), top axis: chainage (km), bottom axis: x-coordinate of the PL-1992 coordinate system (EPSG 2180).

Table 1. Cliff retreat dynamics (m/yr) for the Jastrzębia Góra Cliff (JGC; S1) between 132.25 and 134.24 km chainage.

^a Cliff outline inaccuracies within inspected maps.

The broadest deformations (Fig. 5B–G) impacted the slope between April 2010 and November 2012, linked with storm winds at the end of July 2010 (up to 10° B), a storm surge in January 2012 along with extreme rainfalls in 2010 and 2012, followed by temporary slope equilibrium and the alignment of the cliff main scarp. Consequently, since 2012, the cliff retreat and the deformation pace of the MCSS (Fig. 6) were lower, which reflects beach nourishment in 2013, low annual rainfall totals, and a calm winter in 2014/2015. Vertical displacement (subsidence) and the horizontal component (seaward push) (Supplementary Figs. 5 and 6) of the MCSS edges exhibit significant distortions. The maximum subsidence in 2010–2015, equalled 3 m (3rd step, i.e., the third elevation level of the MCSS structure) in comparison to about 4.5 m of total subsidence. Seaward motion reached 3.5 m, while the maximum total horizontal displacement equalled 6 m from the initial state. The subsidence rate and push slowed down in 2012–2016 (Fig. 6), whereas increased deformations in 2016–2017 were associated with high rainfall totals (daily, monthly, and annual) and storm surges accompanied by high-velocity surface currents responsible for debris removal in summer 2016 and an increased H _s in December 2016 and October 2017. Deformations induced by the reactivation of landsliding processes in the cliff protected by the MCSS (Fig. 4) are almost unrelated to the direct influence of the sea, other than the regular vibrations produced by waves and transmitted en bloc through the ground and synchronized with the wave impacts upon the coast, as observed in numerous water seepages within the cliff and below the landslide foot.

The process of cliff retreat locally occurs rapidly due to the influence of landslides developing where quick clays are present (Figs. 5 and 6). Providing favourable hydrometeorological conditions occur, the gabion band at the cliff base in the vicinity of the most active landslide (134 km), equalling about 2 m asl (while farther to the east it is 3.5–4.0 m asl), slows down the erosion rate by limiting the impact of the storm waves to the landslide foot, but eliminates neither slope deformations associated with landslide development nor the crown retreat. This complex landslide associated with glaciotectonic deformations poses a significant threat to the community. Its side embayments near the promontories are connected with the natural process of maintaining the cliff equilibrium that eventually creates an even larger embayment, along with the mechanism described by Poulton et al. (2006). This is why the nearby promontory at 133.9 km and the viewpoint at the end of the seafront path (Światowid Promenade) are greatly threatened (Fig. 4). Therefore, the stability of the currently intact promenade ought to be carefully examined, whereas its state requires real-time geohazard monitoring and warning systems. The threatened buildings include the primary school, governmental infrastructure (“Baltic”), guesthouses, and residential buildings (Fig. 7) where the distance to the cliff edge ranges from 12 to ~25 m. Vulnerable parts of the Swarzewo Moraine include Rozewie Cliff (130.75–132.08 km) and Chłapowo Cliff (127.0–130.75 km), especially the landslides below the military base in Chłapowo and the lighthouse in Rozewie, despite protection with a concrete breakwater.

Short- and medium-term erosional trends indicate that the cliff retreat dynamics are average rather than extremely high, yet it needs to be considered that the slope margin lately (since the 1990s) is protected by hydro-technical facilities. Also, the lithology of the cliff hinterland suggests that erosion of the western part of the cliff is likely to exert a further negative impact on the local infrastructure in the upcoming decades and centuries. Therefore, widespread discussion concerning sustainable management of this area is necessary.

Swarzewo palaeomoraine—the Mid-Holocene extent and dynamics reconstruction

The hydroacoustic data from the littoral zone, including bathymetric maps (MBES), sea bottom typology based on side scan sonar (SSS) data, and seismic profiling, analysed within the study (Supplementary Fig. 7) allowed a palaeogeographic reconstruction to be achieved. Accordingly, the Mid-Holocene extent of the northwestern part of the Swarzewo Moraine reached up to 2 km from the current coastline (Fig. 8, compare with Supplementary *.kmz image). The western margin of the palaeomoraine is indicated by the seabed topography and lithology (Fig. 8, Supplementary Figs. 7 and 8). Adjacently to the Karwia Lowland, a structural boundary exists between the Pleistocene and Miocene deposits (Supplementary Fig. 7). Recognized relics of moraine relief of the Swarzewo Moraine and its foreshore include topolineaments associated with Miocene outcrops overlapping with residua (boulder fields) and consequently define the extent of the former palaeoshoreline in the Mid-Holocene (Fig. 8). Moreover, the azimuth of the western boundary of the Swarzewo Moraine forms a natural sub-bottom continuation that also precisely complies with the western boundary of the rocky residuum (Figs. 1B and 8). In such cases, the palaeomoraine extents are marked with a solid red line. Otherwise, the dashed red line indicates the location of expected (approximate) boundaries. Thus, moraine total recession dynamics coefficient C _{mid−H→b2k}(S1) = 0.17 ± 0.020 m/yr, revealed in the Desmos and validated in the GIS environment, with a calculation compliance of 88% between the given instances. Therefore, the total recession dynamics of the Swarzewo Moraine C _{mid−H→b2k}(S1) closely match short- and medium-term erosional trends for the JGC.

Figure 8. The extent of the northeastern part of the Swarzewo palaeomoraine in the Mid-Holocene during the Atlantic period, approximately 8 ka b2k, f ₂₀(x), [ R² = 1], f ₂₁(x), [R² =   0.87], f ₂₂(x), [R² = 0.89], dynamics coefficient C _{mid−H→b2k} = 0.17 ± 0.020 m/yr and R² statistics were calculated in the Desmos environment. The Desmos button contains a hyperlink to the map in Desmos environment, with unlocked polynomial nodes allowing for manual modelling.

Site2—Gdynia Orłowo region (S2)

The elevation of the GOC ranged between 12 and 60 m asl. Due to considerable lithological variability, the GOC is exposed to mass movements, including falls, topples, landslides, and earth flows. In the southern part of the cliff, a large landslide of bipartite topography (Figs. 9 and 10) imprints its signature of past wasting processes. Currently, most of the colluvium in its older, southwestern part (1^st) is absent due to debris removal through coastal and subaerial processes. The more recent northeastern part located below the 9- to 20-m-high secondary main scarp (2^nd) has developed into a rotational slump whose foot is subject to intense erosion (at 81.33–81.41 km), whereas the sharp lithological boundary that reaches below sea level near 81.33 km emerged due to landsliding. Remarkably, deposits within the displaced colluvial packet are not preserved in situ, as the actual position of strata differs by minus 9–20 m.

Figure 9. Multitemporal digital terrain models (DTMs) of the Gdynia Orłowo Cliff illustrating total volume changes in 2010–2015: (A) bird's-eye view; (B–D) perspectives.

Figure 10. Erosion intensity of the Orłowo Headland: (A) modified from Subotowicz (1982) by Zaleszkiewicz and Pikies (2007); (B) based on terrestrial laser scanning (TLS) data collected in 2010–2015.

From 2010 to 2015, cliff deformations (Fig. 9) affected colluvia near the 81.33 km and 81.55 km chainage and the secondary scarp retreated by up to 7.5 m due to after-slides within the primarily displaced landslide packet at 81.33–81.41 km. The presence of the latter landslide is associated with the spring discharge through an erosional channel (Supplementary Fig. 11). The Orłowo Headlands (OH) cliff face also experienced almost equally intense erosion, leading to occasional falls and topples due to the systematic development of a wave-cut notch at the cliff's foot.

Farther to the north, an intensive erosion of the OH was related to the low elevation of the shore platform at the promontory base, which drops below 0.5 m asl (H _min = 0.27 m asl). On April 20, 2017, the mean elevation of the GOC base (H _mean) equalled 0.93 m (σ = 0.21 m), while 67% of the cliff base length is at an elevation of 0.5–1.0 m asl, 32% at 1.0–1.5 m asl, and 1% below 0.5 m asl. Despite the headlands’ higher resistance (due to lithologic composition), impacts of even lesser waves below the warning state (<0.50 m) favour the wave-cut notch deepening. In 2010–2015, retreat dynamics of the OH equalled 9.6 ± 0.5 m³/m/yr, while the landslide to the south receded at a rate of 6.0 ± 0.3 m³/m/yr, whereas the maximum retreat of the headland equalled 4.3 m (Fig. 10). However, multitemporal DTMs (Fig. 9) show an uneven erosion pattern across the promontory. Therefore, the dynamics of cliff retreat differ significantly between cross sections. Volumetric calculations for the central and southern GOC (81.3–81.6 km) indicate that sediment delivery to the littoral zone in 2010–2015 equalled 6200 m³, which is about 4.1 ± 0.2 m³/m/yr of deposits. The average retreat dynamics of the OH (81.45–81.55 km) equalled C _2010→2015 = 0.23 m/yr, while an interannual retreat pace along the whole cliff equalled C _2014→2017 = 0.26 ± 0.013 m/yr and C _2000→2015 = 0.68 ± 0.1 m/yr (see Tables 1 and 2 for cliff outline inaccuracies within inspected maps). However, the cliff outline derived from topographic maps in comparison to LIDAR data, erroneously suggests local seaward advance (!) of the cliff top by 20 m. Yet an average recession during the last century equalled C _1908→2015 = 0.25 ± 0.02 m/yr; even though the historical German maps also include inaccuracies, these values correlate well with dynamics recognised within the most recent time frame (Table 2).

Table 2. Average retreat dynamics for the Gdynia Orłowo Cliff (GOC; S2) between 81.15 and 82.08 km chainage.

^a Cliff outline inaccuracies within inspected maps.

The approach utilised in this study illustrates the relationship between landsliding processes and increased dynamics of coastal erosion. Due to the geomorphologic features of the more recent part of the slump, its further retreat is expected until the landslide colluvium below the secondary scarp is completely removed. Consequently, the course of the educational path that leads across the main body of the slump and the viewpoint at 81.35 km require altering unless real-time geohazard monitoring and warning systems become discussed and established.

Redłowo palaeomoraine—Mid-Holocene extent and dynamics reconstruction

Furthermore, based on collected materials, a palaeogeographic reconstruction was accomplished (Fig. 11, compare with the Supplementary *.kmz image). Accordingly, the northern and southern limits of the research area, utilised for calculations, correspond with the current extent of the GOC. In its northern section, the seaward boundary of the residuum follows the course of the isobath 8 m bsl, while the southern part matches the isobaths representing depths of 8–10 m bsl, shallower than in the Swarzewo Moraine foreshore (Fig. 8). Therefore, the seabed adjacent to the Redłowo Moraine must have been subject to less substantial transformations and/or its abrasion started during the later phase of the Atlantic period. Consequently, moraine recession dynamics since the Mid-Holocene C _{mid−H→b2k}(S2)equalled 0.11 ± 0.005 m/yr. This parameter was calculated using split polynomial functions and raster base projection (in accordance with Frydel, 2022; Supplementary Material) in the Desmos environment and validated in the GIS environment, with a calculation compliance of 95% between the given instances.

Figure 11. Redłowo palaeomoraine range in the Mid-Holocene during the Atlantic period, approximately 8 ka b2k, dynamics coefficient C _{mid−H→b2k} = 0.11 ± 0.005, <R ² ∈  0.94, 0.97>, calculated in the Desmos environment upon the base map projected by 90° counterclockwise. The Desmos button contains a hyperlink to the map in Desmos environment, with unlocked polynomial nodes allowing for manual modelling.

Discussion

Electrical resistivity tomography models within S1 predict the extent of quick clays (fine glaciolacustrine deposits) behind the MCSS to reach at least 135 m to the south (Kamiński et al., 2012), with the top of the clay layer tilted towards the north, which additionally favours landslide development (Kamiński et al., 2023). Yet drill-hole data suggest an even wider extent of the zone susceptible to a landslide that stretches as far as 750 m along the coast (133.7–134.45 km) and 370 m to the south. Palaeo-reservoir deposits were not recognized in only two boreholes (P–1, P–2) (Fig. 7). Thus, the northern boundary of the zone corresponds with the maximum seaward extent of the landslide-prone area about 60–70 m from the current cliff edge, matching a recession rate of about 1 m/yr (Subotowicz, 2000). Within S2, fieldwork and modelling explained a discontinuity in Pleistocene deposits (acknowledged by Kaulbarsz [2005] and Woźniak et al. [2018]) and displayed geohazard zones within the cliff, with the main landslide location having a tourist viewpoint on its colluvium, highlighting the threat the landslide poses. Revealed headland retreat dynamics match values for previous decades described by Subotowicz (1982) and Kwoczek (2007), while the erosion amounts correspond with Bełdowska et al. (2016) for the period 2009–2013, when about 5.6 ± 0.3 m³/m/yr became displaced. The retreat dynamics of the JGC and the GOC shown in this study are comparable to the average values determined for the western part of the Polish coast in Trzęsacz, Wolin Island, and Chłapowo (Table 3); however, the pace is several times up to an order of magnitude lower than in Orzechowo near Ustka Bay (Frydel et al., 2017).

Table 3. Retreat dynamics of the coastal cliffs in northern Poland.

^a Age: P, Pleistocene; M, Miocene; N, Neogene.

^b Compiled from Frydel et al. (2017).

^c Compiled from Łabuz (2017).

^d Compiled from Kostrzewski et al. (2015).

^e Compiled from the current study.

Taking a broader perspective, the Estonian sandstone cliffs in Pakri also retreat at a pace of 0.25 m/yr, with the erosion beginning dating back to approximately 7 ka b2k (Orviku et al., 2013), while 8 ka b2k was used in the current study. The accuracy of the time range adopted here can be verified through the dating of the cosmogenic content of ¹⁰Be of an active shore platform that allowed determination of the rate of chalk cliff erosion in East Sussex (southern England) at 0.02–0.06 m/yr during most of the Holocene, followed by an increase to 0.22–0.32 m/yr over the last 150 yr (Hurst et al., 2016). In the case of the SBS littoral zone, further research is required to provide an answer to the question concerning the rate of lowering of an active shore platform as a result of storm wave erosion. Such characteristics could be verified by specifying equations of critical shear stress velocities (u* _cr-Shields) and critical current velocities (u _{cr-Hjulström}) (Ziervogel and Bohling, 2003) for the SBS environment.

Recognition of development trends of cliff shores can be accomplished using transects based on Bayesian modelling that takes into account the probability of a given scenario (Hapke and Plant, 2010). Importantly, a different approach, established on alongshore (not transect-based) explorations using functions, was adopted to identify coastal development trends in the Anthropocene (forecasting) and to determine and simulate the dynamics of coastal transgression and regression during the evolution of the Baltic Sea in the Holocene as well as the advance and retreat of the ice sheet in the Pleistocene (hindcasting), executed as part of 4F MODEL (Frydel, 2022). Reconstructions of the recession dynamics are also achievable following Tylmann et al. (2022). To some extent, hindcasting and forecasting of the dynamics of the aforementioned processes recognized within the 4F MODEL framework can be considered a competitive approach to Bayesian modelling; however, in its current form, it does not take into account the probability of occurrence of a given scenario.

Coastal zone development dynamics linked to global climate change resulting in melting glaciers and rising global ocean levels (HELCOM, 2013; IPCC, 2018; Supplementary Fig. 12) also depends on a spectrum of local and regional factors increasing coastal erosion. Evaluation of the response of the coast to the hydrometeorological conditions (Dudzińska-Nowak, 2017), including that of soft-rock cliffs to recent climatic changes (Carpenter et al., 2014), is among the most important factors influencing local economies, especially land-use planning of human settlements (Łabuz, 2012). In the case of the SBS, the application of the M3D_UG model highlighted disparities between the hydrodynamic conditions of the SBS and the GOG that reflect the relationship between the storm surge count and intensity, surface current velocities, and degree of slope deformation, especially in the areas impacted by landslides, corresponding to erosion dynamics in individual years. Models used to describe the properties of storms on the Polish coast (Kowalewski and Kowalewska-Kalkowska, 2011) and to predict future Baltic Sea levels (Kowalewski and Kowalewska-Kalkowska, 2017) portray a comparably less-severe wave regime and inferior significant wave heights (H _s) within the GOG. Calmer conditions occur due to moraine topography, longitudinal GOC course, and land use, resulting in a wind-shielding effect by the Redłowo Moraine protecting the coast from prevailing westerly winds. Results show consistency with the highest waves of 7.4 m (H _s) measured for the SBS (Tuomi et al., 2011), while the maximum modelled value may reach more than 9 m (Alari, 2013), which indicates that the highest wind waves characteristic for the relatively shallow SBS may equal half the height of the ones encountered within the world ocean (Hanafin et al., 2012).

Given that in the early phase of upland development in the Mid-Holocene (8 ka b2k), the SBS level rise was faster than in the last centuries and millennia, it is reasonable to suspect that, initially, the recession dynamics surpassed the estimated averages. Yet the results of short- and medium-term analyses indicate that the recognized trends match the retreat dynamics for the scale of millennia (C _{mid−H→b2k}). Moreover, the current results may erroneously suggest that the revealed retreat dynamics are average. However, given that in the long term, the rate of sea-level rise plays a predominant role in the rate of coastal erosion (Trenhaile, 2010) and the estimated coefficients (C _{mid−H→b2k}) determine the total average, recession dynamics in the last decade can be considered to be high. Lately (1885–2005), the sea-level rise rate equals 1.6 mm/yr (normalised readings from the Gdańsk tide gauge), which is twice as low as the 3–4.25 mm/yr in the initial phase (8–6 ka b2k) of the Littorina transgression (Uścinowicz, 2003, 2006; Bennike and Jensen, 2011; Lampe et al., 2011), and significantly higher than the 0.33–0.66 mm/yr for the subsequent 6 ka. Therefore, in light of the climate changes resulting in increased wind speed (Groll et al., 2017), sea-level rise, and higher frequency of extreme storms and heavy rainfalls (HELCOM, 2013) directly affecting moraine erosion, it can be assumed that the magnitude and speed of cliff recession will also gain momentum in the coming decades. Forecasts prepared by extrapolating the latest trends for the next decade indicate that the cliffs in Gdynia Orłowo and Jastrzębia Góra, but also in Chłapowo and for Wolin Island, are expected to recede by approximately 2 m inland, while in landslide-prone areas, the extent may be an order of magnitude higher.

Conclusions

The influence of variable environmental conditions resulted in relatively constant coastal erosion dynamics within the inspected uplands (Swarzewo Moraine, S1; Redłowo Moraine, S2) at all analysed timescales throughout the last 8 ka (Supplementary Fig. 12A). These results confirm that the gradual rate of sea-level rise since the mid-Atlantic period (Uścinowicz, 2006) plays a predominant role in the long-term rate of cohesive coast erosion (Trenhaile, 2010), also within the SBS. Short-term retreat dynamics of inspected margins, C _2008→2017 (S1) = 0.17 ± 0.020 m/yr and C _2014→2017 (S2) = 0.23 ± 0.01 m/yr, and medium-term pace, C _1875→2015 (S1) = 0.19 ± 0.02 m/yr, and C _1908→2015 (S2) = 0.25 ± 0.02 m/yr, matches long-term recession trends at the millennial scale (since 8 ka) for each site, equalling C _{mid−H→b2k}(S1) = 0.17 ± 0.020 m/yr and C _{mid−H→b2k} (S2) = 0.11 ± 0.005 m/yr (Supplementary Fig. 12), as revealed a through large-scale application of the 4F MODEL (Frydel, 2022). Given outcomes may enrich the GlobR2C2 database by Prémaillon et al. (2018).

For the S2, short- and medium-term trends suggest higher retreat dynamics as compared with the long-term trends. The former extent of the Swarzewo palaeomoraine reaches 1300–2000 m into the sea, which corresponds with the Usedom Island palaeomoraines (on the German coast) that stretch 1800–2400 m into the sea (Hoffman and Lampe, 2007). The maximum extent of the Redłowo palaeomoraine of 500–1200 m away from the present coastline and its total recession coefficient values are lower, which results from a less severe long-term hydrometeorological regime at S2 than at S1, although the coastal retreat pace at S2 surpassed S1 short- and medium-term trends, due to lack of anthropogenic impacts.

Yet the development of the inspected coastal system is influenced by factors on local, regional, and global scales (Supplementary Fig. 12B). In places where landslides and falls occur, the retreat of moraines progresses approximately an order of magnitude faster than the average value for the whole cliff coast at a given site. Therefore, landslides ought to be treated as one of the most vital components of the local cliff retreat mechanism. Development at a slower pace may proceed despite slope protection measures, whereas quick clays present in the subsurface and intense precipitation (daily totals >80 mm) contribute to increased infiltration and pore pressure that trigger landslide reactivation, resulting in the downslope transport of liquefied colluvia subject to subsequent removal by the wind–wave regime.

The existence of quick clays at the study sites coincides with similar floodplain, fine-grained sediments of meandering rivers (approximate longitudinal flow direction from the east to the west) as described by Sokołowski et al. (2019) for the southern part of the Swarzewo Moraine. Within its northeastern part, the river course with flow direction from the northeast to the southwest explains irregular occurrence of landslides along the coast and rapid cliff retreat over the space of several dozen metres to more than 100 m (equal to the width of river streams), consistent with Bennett and Glasser (2009). However, confirmation of this thesis requires further research.

Given the increasing threat to local communities and infrastructure caused by climate change and sea-level rise, especially in landslide-prone areas, systematic monitoring and warning systems for erosional tendencies in the most threatened zones need to be included as a key element of sustainable coastal management. This need especially applies to landslide-prone zones: at the S1 study site within the JGC, at the end of the Światowid Promenade, and in the central and southern part of the GOC. The application of the approach described in this article is particularly important to the ongoing monitoring of the Polish coastal zone—the 4D Cartography Project (Kramarska et al., 2014; Uścinowicz et al., 2018). However, the designed scope of investigation, which extends 2 km towards the sea, may be insufficient to determine the seaward limits of some shore platforms, thus preventing the recognition of the total recession dynamics of the moraines.