Aegean-style extensional deformation in the contractional southern Dinarides: incipient normal fault scarps in Montenegro

We describe two 5–7 km long normal fault scarps (NFSs) occurring atop fault-related anticlines in the coastal ranges of the Dinarides fold-and-thrust belt in southern Montenegro, a region under predominant contraction. Both NFSs show well-exposed, 6–9 m high, striated, and locally polished fault surfaces, cutting uniformly northeastward-dipping limestone beds at high angles and documenting active faulting. Sharply delimited ribbons on free rock faces show different colors, varying karstification, and lichen growth and suggest stepwise footwall exhumation, which is typical of repeated normal faulting during earthquake events. Displacements, surface rupture lengths, and geometries of the outcropping fault planes imply paleoearthquakes with M_w≈6 ± 0.5 and slip rates of ∼ 0.5–1.5 mm yr⁻¹ since the Last Glacial Maximum. This is well in line with (more reliable, higher-resolution) slip rates based on cosmogenic ³⁶Cl data from the scarps for which modeling suggests 1.5 ± 0.1 mm yr⁻¹ and 6–15 cm slip every 35–100 years during the last ∼ 6 kyr. The total throw on both NFSs – although poorly constrained – is estimated to ∼ 200 m and offsets the basal thrust of a regionally important tectonic unit. The NFSs are incipient extensional structures cutting (and postdating emplacement of) the uppermost Dinaric thrust stacks down to an unknown depth. To explain their existence in a region apparently under pure contraction, we consider two possibilities: (i) syn-convergent NFS development or – less likely – (ii) a hitherto undocumented propagation of extensional tectonics from the hinterland. Interestingly, the position of the extensional features documented here agrees with geodetic data, suggesting that our study area is located broadly at the transition from NE–SW-directed shortening in the northwest to NE–SW-directed extension to the southeast. While the contraction reflects ongoing Adria–Europe convergence taken up along the frontal portions of the Dinarides, the incipient extensional structures might be induced by rollback of the Hellenic slab in the southeast, whose effects on the upper plate appear to be migrating along-strike of the Hellenides towards the northwest. In that sense, the newly found NFSs possibly provide evidence for a kinematic change of a thrust belt segment over time. However, with a significantly higher probability, they can be regarded as second-order features accommodating geometrical changes in the underlying first-order thrust faults to which they are tied genetically.

1 Introduction

Active normal faults in the Mediterranean frequently develop bedrock normal fault scarps (NFSs). Their suitability as tools for paleoseismic analyses has been proven by many authors (see the following passages for differentiated references). We describe two 5–7 km long NFSs along the southwestern slopes of the Rumija mountains in the coastal ranges of the Dinarides fold-and-thrust belt in southern Montenegro on the western Balkan peninsula (Figs. 1–4; S1 and S2; Table S3). According to their positions between the eponymous towns, we refer to them as Bar (BFS) and Katërkolle (KFS) fault scarps. The primary aim of this paper is a first-time description and interpretation of the mentioned structures, including (i) verification that BFS and KFS are active NFSs at all, (ii) basic studies on the timing of NFS exhumation and slip rates, and (iii) a discussion on how these active NFSs may be embedded in the contractional geodynamic setting of the southern external Dinarides. Our work is based on a thorough mapping campaign, involving both ground truthing and remote sensing, followed by low-spatial-resolution ³⁶Cl cosmogenic dating of NFS free-face samples. Both paleoseismic interpretations based on overall NFS geometries and visual indicators (e.g., Armijo et al., 1992; Giraudi and Frezzotti, 1995; Roberts and Michetti, 2004; Papanikolaou et al., 2005, 2013; Papanikolaou and Roberts, 2007; Faure Walker et al., 2009; McCalpin, 2009; Grützner et al., 2013, 2016; Bubeck et al., 2015; Mason et al., 2016, 2017; Mechernich et al., 2018, 2022) as well as ³⁶Cl cosmogenic dating (e.g., Benedetti et al., 2002, 2013; Palumbo et al., 2004; Schlagenhauf et al., 2010, 2011; Roberts et al., 2014; Tesson et al., 2016; Cowie et al., 2017; Mechernich et al., 2018, 2022; Goodall et al., 2021; Iezzi et al., 2021) have previously been successfully applied for NFSs in the Mediterranean region, delivering valuable instructions and benchmarks for our study. BFS and KFS closely resemble the abundant and well-studied NFSs in the central Apennines (Italy; e.g., those ruptured during the 2016–2017 earthquake series) as well as the Greek part of the Hellenides in terms of geomorphology, structural setup, and length, showing evidence of segmentation and repeated co-seismic footwall exhumation during earthquakes (e.g., Papanikolaou et al., 2005, 2013; Grützner et al., 2016; Mason et al., 2016, 2017; Civico et al., 2018; Mechernich et al., 2018, 2022). While the formation of large-scale NFSs in the extensional geodynamic settings of Italy and Greece is barely surprising and the NFSs correlate with strong extensional earthquakes, it is puzzling that all known major instrumentally recorded earthquakes between southern Croatia and Albania (including Montenegro 1979, M_w≈ 7.1, and Durrës–Albania 2019, M_w≈ 6.4; e.g., Benetatos and Kiratzi, 2006; Papadopoulos et al., 2020) were exclusively contractional (Fig. 1). Two different explanations are discussed: (i) the NFSs are formed syn-convergently, a phenomenon that has been frequently observed in other locations (e.g., Philip and Meghraoui, 1983; Nábělek, 1985; Avouac et al., 1992; Bennett et al., 2012, Hicks and Rietbrock, 2015; Corrado et al., 2019; Riesner et al., 2021), or (ii) the NFSs indicate a possible temporal transition in the kinematic behavior of a mountain range – a phenomenon rarely documented in detail. The hitherto inchoate understanding of the newly discovered Montenegrin NFSs and the fact that (i) similar NFSs in Italy and Greece are associated with major earthquakes that caused many casualties and severe economic losses from destroyed medieval villages and modern infrastructure (e.g., Chiaraluce et al., 2017, Table S4) along with the fact that (ii) comparably small fault systems are generally underestimated in terms of their seismic hazard (Grützner et al., 2013) underline the urgent need to further analyze such structures. In the following, we present the methods and results of NFS mapping (structural and surface exposure indicators), cosmogenic ³⁶Cl sample collection, preparation, and AMS (accelerator mass spectrometry) dating, modeling of the ³⁶Cl data, and the resulting tectonic interpretations.

Figure 1 (a) Earthquakes in the study area: where available, fault plane solutions (FPSs) are scaled and color-coded according to magnitude and stress regime. All other earthquakes are color-coded by grey scales according to magnitude. SPFZ: Shkodra–Peja fault zone. (b) FPS projected onto profile A-A^′ within the range of the white box (Fig. 1a). FPS from EMSC, EMRCMT and Harvard catalogues, Pondrelli et al. (2006), and Grünthal et al. (2013).

2 Geological setting

The BFS and KFS are located in the extreme south of the Dinarides, slightly north of the Dinarides–Hellenides transition. The latter is marked by the ca. NE–SW-striking normal–transverse Shkodra–Peja fault zone (SPFZ, e.g., Handy et al., 2019, and references therein). The driver for the seismicity along the coast of Montenegro and Albania is the northward motion of the Adriatic microplate that bidirectionally subducts below the Balkan and Apennine peninsulas, creating almost mirror-imaged tectonic settings on both sides of the Adriatic Sea (e.g., Nocquet and Calais, 2004; Battaglia et al., 2004; Faccenna et al., 2014; Le Breton et al., 2017; Király et al., 2018). Both the Dinarides–Hellenides and Apennine fold-and-thrust belts are characterized by orogen-parallel, NW–SE-striking tectonic units. NE–SW-directed contraction along the deformation fronts is replaced by extensional domains in the hinterland (Fig. 1; e.g., D'Agostino et al., 2008; Nocquet, 2012) that are attributed to subduction rollback, gradually migrating towards Adria in both cases (e.g., Cavinato and de Celles, 1999; Dumurdzanov et al., 2005; Carminati and Doglioni, 2012; Handy et al., 2019). The restriction has to be made that the latter statements are not fully valid for the southern Apennines, where continental subduction recently stopped (Patacca and Scandone, 2007) to replace frontal thrust belt motion by strike-slip faults involving the lower and upper plate (Di Bucci et al., 2006). Despite all similarities between the Dinarides–Hellenides and Apennines, seismotectonic characteristics reveal major differences. In Italy, the most destructive earthquakes are commonly normal faulting events, creating distinct large-scale NFSs in a pronounced, continuous hinterland extensional domain (e.g., Chiaraluce et al., 2017; Galadini and Galli, 2000). On the Balkan peninsula, the greatest risk emanates from contractional earthquakes with epicenters close to the Montenegrin–Albanian coast (e.g., Pondrelli et al., 2006; Copley et al., 2009; Papadopoulos et al., 2020; Vittori et al., 2020; Fig. 1). Known extensional tectonics and related earthquakes are restricted to the internal Hellenides south of the SPFZ, resulting from a clockwise rotation of the Hellenides segment with respect to the Dinarides (e.g., Jouanne et al., 2012; Faccenna et al., 2014; Handy et al., 2015). NFSs are poorly developed in these regions of Albania (Handy et al., 2019), while they were hitherto entirely unknown in Montenegro. The existence of BFS and KFS (Figs. 2–4, S1 and S2; Table S3) is therefore striking in two respects: (i) they are a rare example of well-developed NFSs in the Dinarides–Hellenides north of Greece. (ii) They are – surprisingly – not located in the hinterland, where extensional focal mechanisms are well documented, but as close as 4 km from the coast in a fold-and-thrust belt segment solely characterized by horizontal contraction, as evidenced in both the geological structures and the predominance of reverse faulting focal mechanisms (Fig. 1). Structurally, the BFS and KFS cut the uniformly northeastward-dipping limestone beds of the Budva–Cukali, High Karst, and (to a very small extent) Kruja–Dalmatian units at high angles. The Budva–Cukali Unit – hosting the major proportion – is a regionally important tectonic nappe in the Dinarides–Hellenides consisting mainly of Mesozoic pelagic sediments topped by Paleogene synorogenic deposits (Fig. 4). Its emplacement onto the underlying Kruja–Dalmatian Unit is dated to between late Eocene and middle Miocene (e.g., Schmitz et al., 2020).

Figure 2 Panoramic view from aboard the Meteor research vessel cruise no. 86, leg 3 SE towards the Montenegro coast, including BFS_N. Basal thrusts separating the large-scale tectonic units and a choice of other landmarks are sketched.

Figure 3 Photographs of the northern section of the Bar fault scarp (BFS_N): (a) panoramic view, with the fault trace shown in the inset. (b) Repeated earthquake ribbons; hammer for scale. Arrows: blue – young; yellow – older. (c) Fault plane outcrop with characteristic corrugations and oblique-slip lineations; person for scale. (d) Cataclastic fault breccia in the footwall of the fault plane; hammer for scale. (e) Sampling site BFS_N, with a 2 m ruler for scale; for localization see Fig. 4.

In the study area, along the Rumija mountain front, only the uppermost part of its stratigraphic section appears as a ≤ 5 km wide corridor between the structurally underlying Kruja–Dalmatian and the overlying High Karst tectonic units. Furthermore, the Budva–Cukali Unit appears in remnants at the base of the High Karst Unit in two isolated nappe outliers (Fig. 4). For more detailed information on the regional geology of the area, the reader is referred to, e.g., Biermanns et al. (2019), Handy et al. (2019), Schmid et al. (2020), and Schmitz et al. (2020). Current horizontal shortening rates for the region lie in the range of 3–5 mm yr⁻¹ (Kotzev et al., 2008; Jouanne et al., 2012; Devoti et al., 2017), while vertical uplift rates were estimated to around 1 mm yr⁻¹ (Biermanns et al., 2019, and references therein).

Figure 4 (a) Localization of BFS and KFS on the geological map (compiled from the 1:25 k geological map sheet Vladimir 160-4-3 and 1:50 k geological map sheet Podgorica 3; Geological Survey of Montenegro, 2009). Stars indicate the sites of NFS profiles (Fig. 5). Stereoplots show orientations of the fault planes and striations for BFS_N (grey), BFS_S (orange), and KFS (green); semi-transparent coloring – range of main fault plane orientations; dots – mean orientations of the striation. Single great circles represent (i) smaller transfer faults that connect segments of the main fault plane and (ii) the fault branch at the western end of KFS. HK: High Karst, B-C: Budva–Cukali. (b) Cross section showing the hypothetical minimum fault slip derived from the reconstructed visible fault scarp height at sampling site KFS (Fig. 5, bottom). (c) Cross section showing the proposed maximum fault slip derived from the maximum offset of the Budva–Cukali marker between the outcrops at the Vladimir nappe outlier and the Rumija range.

3 Methods

3.1 Fieldwork

3.1.1 Structural and morphological mapping, site selection, and fault scarp profiling

The full extent of BFS and KFS was geologically and structurally mapped to gain details of fault morphology and orientation, fault striae, surface roughness, and lichen growth (Figs. 4; S5–S8). Particular attention was directed to the identification and mapping of striation-perpendicular, ca. slope-parallel horizons (“ribbons”) on the fault planes, interpreted to display co-seismic displacements per earthquake event (Figs. 3b and S6; more detail is described in Sect. 4.1). Based on the mapped fault lengths (Table S9) and an evaluation of the mentioned horizons, we estimated earthquake magnitudes after Wells and Coppersmith (1994). Furthermore, four representative sites (BFS_N, BFS_S1, BFS_S2, and KFS; Figs. 4 and S1; Table S3) were selected (i) to collect samples for ³⁶Cl dating (site BFS_N only; see following Sect. 3.1.2) and (ii) to construct topographic profiles across the fault scarps as a basis to estimate low-resolution long-term (post-Last Glacial Maximum, LGM) slip rates for comparison (Fig. 5, following examples of, e.g., Papanikolaou et al., 2005; Mason et al., 2016; see also Sect. 5.3). For this technique, it is assumed that the preservation of NFSs initiated ∼ 15 ± 3 kyr ago, following the demise of the LGM in the Mediterranean region (e.g., Benedetti et al., 2002; Papanikolaou et al., 2005; Giraudi and Frezzotti, 1995). Until then, periglacial conditions allowed slope-degrading processes to exceed fault throw rates. Post-LGM warming, waning freeze–thaw cycles, and slope stabilization by vegetation allowed fault throw to outpace slope degradation, thus forming pronounced NFSs (e.g., Papanikolaou et al., 2005). The four sites fulfill all essential requirements like minimum erosion and deposition, flat scarp base, intact scarp surface, and representative striations (e.g., Bubeck et al., 2015; Cowie et al., 2017; Mechernich et al., 2018). The fault planes were cleared from vegetation, followed by a thorough structural survey including the immediate surroundings. Profiles were measured by a broomstick and clinometer parallel to striations in 1 m steps, ∼ 50 m upslope and downslope of the NFS. The entire NFS height consists of two sections: (i) height of the distinct free-face and (ii) degraded NFS height, interpolated from hanging wall and footwall slope (Fig. 5). Both portions are incorporated in the calculation of a generalized post-LGM (15 ± 3 kyr) slip rate (Table S10; for discussion of LGM timing and interpretation of the degraded fault scarp, see Sect. 5.3). Since the calculation principle integrates the NFS height over a full post-LGM time period, the method only yields hypothetical constant slip rates but does not resolve actual phases of fault activity (e.g., with slip clusters) or quiescence. Despite these weaknesses and the availability of ³⁶Cl dating for site BFS_N (see Sect. 3.1.2), we still consider the comparison of topography-based slip rate estimations to be a reasoned benchmark in the framework of our study: (i) those sites not dated by the ³⁶Cl method are considered. (ii) The degraded scarp at site BFS_N (which is not included in ³⁶Cl dating) can be incorporated in considerations and related to the ³⁶Cl-dated free face. (iii) The four sites are distributed across different sections of BFS and KFS (Fig. 4), each of which has a distinct setup (structurally and in terms of exposition), allowing for a comparison and rating of these sites. (iv) The qualitative calculation normalized to a full post-LGM time frame sets a lower limit in terms of slip rates and earthquake recurrence intervals.

Figure 5 Profiles across the fault scarps at four selected sites (see Fig. 4 for locations). Slip rates are derived from the free-face and degraded scarp heights (compare Table S10, two right columns). Stereoplots show fault plane orientations (great circles) and striations (triangles) within ± 5 m of the study site. Sites are indicated in Fig. 4; sites BFS_N and BFS_S are shown in Figs. 6a–c and S2c.

3.1.2 ³⁶Cl dating: sampling procedure

Samples for ³⁶Cl dating were taken at site BFS_N (introduced in Sect. 3.1.1; Figs. 3e, 6, and S2; Table S3). For extraction of datable rock samples, a trace across the highest part of the outcrop wall was defined parallel to the visible striations (Fig. 6c and f; Mechernich et al., 2018). To achieve an adequate resolution for the reconstruction of long-term slip rates and seismic events, we chose a sample spacing of ∼ 50 cm parallel to slip direction whilst avoiding confounding factors (e.g., joints, Fig. 6f). This sample spacing is rather large compared to previous studies of cosmogenic fault scarp dating, and hence it is not possible to identify the stepwise ³⁶Cl concentration pattern (so-called cusps, e.g., Schlagenhauf et al., 2010). In our study, this is not a disadvantage, since the capable offsets of the small faults are lower than the ∼ 2 m co-seismic offset required for the generation of a stepwise ³⁶Cl concentration pattern. We take the low sample density into account in the modeling approach. As published, such resulting modeled slip rates are comparable to slip rates derived from dense (continuous) sampling (Beck et al., 2018; Iezzi et al., 2021). Generally, the multiparametric approach of our study compensates for the low sample amount, as conclusions and interpretations are not solely contingent upon ³⁶Cl-dating-based input data. The lowermost samples were collected ∼ 1 m below the scarp base by manual excavation of a trench (Fig. 6c). These buried samples are essential to characterize pre-exposure conditions (e.g., Cowie et al., 2017). The previously marked 15×5 cm sample blocks were carefully extracted with the help of an angle grinder, hammer, and chisel (Fig. 6d–f). Subsequently they were marked according to their distance from the scarp base and packed for shipping. To quantify the risk of insolation weathering at the sampling sites, exposure angles were determined in 10^∘ steps with the help of a clinometer.

Figure 6 (a) Location BFS_S1 as an example for which the criterion of a perfectly flat scarp base is met. (b) Difference between overgrown (a, b) and cleaned (e, f) fault plane at sampling location BFN_N. (c) Full view of the cleaned sampling location BFS_N including a trench below the scarp base. (d–f) Work in progress: sample blocks are marked and extracted with the help of an angle grinder, hammer, and chisel at sampling location BFS_N. Image (f) shows the trace of extracted sample blocks parallel to striations but avoiding disturbing factors such as joints.

3.2 ³⁶Cl dating: sample preparation and data modeling

3.2.1 Sample preparation

Based on a spacing of 100–200 cm (striation-parallel distance on the fault plane), six sample blocks from sampling site BFS_N were prepared at the Institute of Geology and Mineralogy of the University of Cologne. Weathered parts and pore surroundings were carefully removed with a rotary tool before crushing and sieving. The following chemical treatment and the measurement at the Cologne AMS facility were performed as described in Mechernich et al. (2018). Resulting ³⁶Cl $/$ ³⁵Cl, ³⁶Cl $/$ ³⁷Cl, and ³⁵Cl $/$ ³⁷Cl ratios were used to calculate the concentrations of ³⁶Cl and natural chlorine (Cl_nat). Their reliability is confirmed by the simultaneous preparation of CoCal–N ³⁶Cl standard material (Mechernich et al., 2019) and one blank in the respective batch. The blank subtractions were 0.8 %–1.7 % (Table S11). The calculated ³⁶Cl concentrations of the six analyzed samples range from ∼ 7 × 10⁴ atoms per gram of rock at 0.55 m below the scarp base to ∼ 2 × 10⁵ atoms per gram of rock at a height of 5.8 m above the scarp base. In general, the concentrations continuously increase with fault scarp height (Fig. 7). The Cl_nat concentrations are very low, ranging from 6 to 17 µg g⁻¹ (Table S11). One replicate sample was prepared and measured in Cologne (Table S11). An aliquot of each dissolved sample was analyzed by in-house ICP-OES (inductively coupled plasma–optical emission spectrometry) at the University of Cologne to determine the concentrations of the principal ³⁶Cl target elements, Ca, K, Ti, and Fe. The ICP-OES Ca concentrations of the BFS_N samples range from 38.9 % to 40.0 %, indicating local variabilities (Table S13) with a minor impact on the ³⁶Cl production rate. We used one non-treated free-face sample from the BFS_N site as a reference for the assumed thermal and epithermal neutron flux and thus constrain production of ³⁶Cl on ³⁵Cl. Equally, trace element analyses on this sample were used for the ³⁶Cl production estimate (Table S13). Both analyses were performed by Actlabs (Canada). For the hanging wall composition, we used the soil composition of the colluvium.

Figure 7 ³⁶Cl concentrations (1σ deviations) as a function of the height up the scarp (distance measured on the free face). (a) Modeled ³⁶Cl concentrations for constant stick-slip rates (1.5 ± 0.1 mm yr⁻¹) using the code of Schlagenhauf et al. (2010), with the highest likelihood. (b) The resulting correlation of age and scarp height of (a). (c) Modeled ³⁶Cl concentration for the “most likely” landslide–rockfall scenario. (d) The resulting correlation of age and scarp height of (b).

3.2.2 ³⁶Cl scarp modeling method and parameters

To determine earthquake ages and slip rates from the ³⁶Cl concentrations, we used the MATLAB® code of Schlagenhauf et al. (2010) that models synthetic ³⁶Cl concentrations while accounting for all influencing factors, i.e., the time-dependent variability of the fault scarp geometry, the chemical composition, and the respective amount and timing of progressive exhumation. The code was adapted to the large sampling spacing, and furthermore, the mapped ribbon heights were added to the input parameters. All input parameters are described in the following passages and in Tables S11–S13. Several parameters have an influence on the production of cosmogenic ³⁶Cl, which typically extends from several meters to tens of meters below the surface. In addition to the chemical composition and density of the bedrock scarp and the colluvial wedge (Table S12), the ³⁶Cl production rate depends strongly upon the rate at which the scarp is exhumed. Continuous accumulation of ³⁶Cl in the footwall rock occurs in both the shallow subsurface (inherited or pre-exposure component) and pre-eminently after sub-aerial exposure and growth of the fault scarp (e.g., Schlagenhauf et al., 2010; Mechernich et al., 2018). This typically leads to increasing ³⁶Cl concentrations with fault scarp height, although this is somewhat complicated by erosion of the scarp free face, whereby ³⁶Cl in the rock is reduced by weathering. Furthermore, production rates have to be scaled appropriately to the local and distant shielding of the site from cosmic rays. The average density of the limestone samples was determined using the sample weight and their volume by suppression in water, yielding 2.55 g cm⁻³. The density of the colluvium was estimated at ∼ 1.5 g cm⁻³; more specific measurements were not undertaken due to local variabilities in the clast occurrence and humidity impact. We used an annual ³⁶Cl production rate of 48.8 ± 3.4 atoms per gram of rock from Ca spallation (Stone et al., 1996) as it is derived from a similar latitude (39^∘ N) and a rather comparable altitude (1445 m a.s.l.), and it integrates over a time span of 17.3 kyr, which is appropriate for our postglacial focus. All further production rates used are given in Table S13. Scaling with respect to latitude and elevation was performed using the Stone (2000) scaling scheme assuming a constant geomagnetic field intensity. The geometry of the fault as derived from the topographic profile (Fig. 5) is used to calculate shielding factors for the time-dependent self-shielding during the progressive exhumation of the fault plane. Thereby, 33^∘ was used for the dip of the hanging wall, 56^∘ as dip of the fault plane, 35^∘ as dip of the footwall, and 22.2 m as the total displacement of the hillslope (Fig. 5). Additional topographic shielding does not occur since the mountains in sight occur just insignificantly above the horizon. There is significant local variation in the amount of weathering of the exposed fault plane, ranging from zero at locations with preserved slickensides (0–3 m of the free-face height) to ∼ 3 cm at the solution flutes. While the fault plane at the scarp base is smooth at a millimeter scale, rock surface relief at 8.0–8.8 m height is 2–8 mm around the sampling line. Assuming this as the minimum amount of erosion and using a preliminary age estimation of ∼ 15 kyr at 8.8 m height, this would correspond to a minimum erosion rate of 0.13 to 0.53 mm kyr⁻¹ at our sampling locations. If the whole fault scarp were of postglacial age (∼ 15 kyr at 22.2 m height, i.e., ∼ 6 kyr at 8.8 m height), the erosion rates would be 0.33 to 1.3 mm kyr⁻¹. Such a low rate was also observed on other fault planes in carbonates (Goldberg et al., 2016; Mechernich et al., 2018, 2022). In our further analysis, we focused on the likely value of 1 mm kyr⁻¹ of erosion. For the numerical modeling of the ³⁶Cl concentration, we applied the MATLAB® code modelscarp.m presented in Schlagenhauf et al. (2010). All parameters used are displayed in Table S14. The code was used by iteratively modeling constant slip rates of the BFS, which fit the measured ³⁶Cl concentrations best. The modeling of “constant” arbitrary slip rates is done by using simple scenarios of a stick-slip behavior with 6–15 cm co-seismic offsets (as suggested by the mapped ribbons; the approach is described in detail in the Supplement of Mechernich et al., 2018). The uncertainty of the constant slip rate is mainly based on the measured ³⁶Cl concentrations, the amount of co-seismic offset, and additionally on external effects which are not incorporated in the given slip rate uncertainty (e.g., deviations of production rate, shielding). We applied the criterion that all scenarios covering the 68 % uncertainties of at least three of the five ³⁶Cl samples are considered so that a 1σ in ³⁶Cl internal uncertainty of this constant slip rate is derived. To approach the lowest possible slip rates on the BFS, we additionally applied a “sliding scenario” for the fault scarp part above the sampled part of the free face, even though we did not find any indicators of such a “sliding event”. Here we also used the code modelscarp.m and iteratively tested what amount of sliding at what time results in the lowest slip rate for the sampled part of the free face. A hypothetical rockslide was modeled using a large offset within a small amount of time (i.e., 1 year). To calculate the ages of the last earthquake events, we used the resulting slip rate on the Bar fault together with the amount of co-seismic offset based on the ribbon locations. An uncertainty propagation was included by considering the uncertainties of the slip rate and the amount of co-seismic offset.

4 Results

4.1 Structure and morphology of the fault scarps

Both the BFS and KFS crosscut the ∼ 30–45^∘ northeastward-dipping beds of the Budva–Cukali, High Karst, and Kruja–Dalmatian tectonic zones at high angles (Figs. 4 and S5). Bedrock limestone in the footwall is juxtaposed against carbonatic, partly cemented colluvium and slope scree in the hanging wall, consisting of centimeter- to meter-sized clasts (Figs. 3b and c). The fault zone is marked by a cataclastic fault breccia of ≥ 1 m width (Fig. 3d). Fault dip directions, pervasive slickenside striations, slickenfibers, Riedel shears, well-developed triangular facets, and wine-glass-shaped valleys (Fig. 3a and c; e.g., Dramis and Blumetti, 2005) prove active normal faulting. Free faces are 1–9 m high in domains with negligible erosion (Figs. 3e, 5 and S2c) and up to 40 m in domains of strong hanging wall erosion, e.g., where gullies occur (Fig. 3a and c). Fault planes frequently exhibit ∼ 5–50 cm wide, horizontal, and sharply bound ribbons of distinct color and roughness, increasingly better preserved toward the scarp base (Figs. 3b and S6). Boundaries between the ribbons are ca. perpendicular to local fault plane striations and therefore mostly slope-parallel. Their distance to the present-day scarp base is nearly constant. Across these boundaries, widths of karstic solution flutes, surface roughness, lichen growth, and micro-karstification decrease stepwise toward the scarp base (Fig. S6). While the lower ribbons are partly correlatable over longer distances, the higher-up ones are often hardly distinguishable and defaced local occurrences. Individual ribbons can be correlated across several locations in terms of their width and habitus (Fig. S7). However, they are rather isolated occurrences that do not enable a gapless tracing along the fault planes. The highest density of perceptible single horizons was encountered on BFS_S. The ribbons are interpreted as single exhumation events (e.g., Mechernich et al., 2018, 2022, and references therein; for further discussions see Sect. 5.1) and therefore used as an input parameter for ³⁶Cl dating (see Sects. 3.2.2 and 4.2.2). The BFS dissects the slopes of Mt. Lisinj (Fig. 2). Approximately midway along its ∼ 5 km extent, it changes its mean fault plane orientation from moderately steep (∼ 55^∘) ca. NNW–SSE-striking in the north to steep (∼ 70^∘) ca. east–west-striking in the south (Figs. 4, S5, and S7a and b). The transitional area between the northern (BFS_N) and southern (BFS_S) section of BFS presents as rather complex, with an apparently diffuse array of multiple fault planes (Figs. 4, S1 and S8; Table S3). Despite fewer conclusive outcrops, the pervasiveness of faults in bare bedrock and relays at the transition between major branches of BFS_N and BFS_S (mapped by means of remote sensing) suggests a likely interconnection of the two BFS sections. The formation of distinct NFS outcrops (like in most other locations) is likely prevented by lower offset as a result of slip distribution among multiple fault branches. Ribbon abundance and widths are comparable along the full length of BFS (Fig. S7a and b). KFS follows the southern slopes of the Rumija mountains for >7 km and crosscuts thrusts at the base and top of the Budva–Cukali zone (Fig. 4). A connection between BFS_S and KFS is conceivable, as suggested by (i) a similar mean fault plane orientation and (ii) an interjacent penetrative step in terrain steepness (Fig. S8). However, a lack of intermittent outcrops for ∼ 3 km along-strike and less-abundant ribbons (Fig. S7c) render such correlation less certain. Along all NFS sections, fault planes reveal systematic undulations and corrugations with wavelengths up to several meters (Fig. 3c). The trends of striations follow the mean fault plane orientation, indicating dominant dip-slip kinematics. A tendency to increasing strike-slip components away from the section centers creates patterns of radially outward-diverging striations (Fig. S5). At several locations along BFS_N and KFS, transfer faults were mapped (Fig. 4a). These faults branch away obliquely from the main fault plane, sometimes at high angles. The free faces formed by these faults are commonly less high (∼ 1–4 m) compared with the main NFS, but otherwise show the same characteristics (slickenside striations, undulations, etc.). The western end of KFS is marked by a large-scale (tens of meters high; cleared by an adjacent stream) fault plane that deviates northward from the main KFS trend by almost 90^∘ (Fig. 4a) and features slickensides indicative of strike-slip movement. Geometric considerations of the thickness of the Budva–Cukali Unit and the degree of folding of its tectonic contacts yield a maximum throw of ∼ 200 m across KFS: Fig. 4b illustrates a hypothetical offset of 16 m (exclusively based on the observed throw at location KFS, Fig. 5), assuming that the outcropping width of the Budva–Cukali Unit at the front of the Rumija mountain range goes back to an unfaulted (∼ 600 m thick) state of the unit. Figure 4c shows a lower thickness of the Budva–Cukali Unit (∼ 380 m), which can only achieve the required outcropping width along the mountain front by artificially thickening the unit by faulting with an offset of ∼ 200 m; this represents the maximum geometrically conceivable throw. With a larger offset and constant thickness, the contact between the Kruja–Dalmatian and Budva–Cukali units would no longer be congruent with the geological map. The opposite (southwestward) dipping nappe contacts between the Kruja–Dalmatian, Budva–Cukali, and High Karst units at the Mt. Vladimir nappe outlier prove the existence of an eroded anticlinal hinge, which supports the scenario shown in Fig. 4c. Thus, an offset close to 200 m is likely. Obtained higher slip rates (Table S9) and topographic cross sections across BFS_N, where clearly perceptible knickpoints (∼ 850 m a.s.l.) mark the NFS ∼ 400 m below the overlying highest parts of the Rumija ridge (∼ 1250 m a.s.l), suggest that the offset across this NFS section may be slightly higher. However, a lack of convenient markers (such as those presented for KFS) prevents sufficiently reliable estimations.

4.2 Proxies for fault scarp exhumation

4.2.1 Slip rate and magnitude estimates derived from fault scarp profiling, surface rupture lengths, and ribbons

For all four NFS sampling sites (i.e., also those that were not dated by means of the ³⁶Cl method), we calculated post-LGM slip rates according to the simplistic procedures described in Sect. 3.1.1. The obtained rates range between 0.51 ± 0.12 (site BFS_S2) and 1.46 ± 0.37 (site BFS_N) mm yr⁻¹. Measurements of earthquake-related ribbons (see also Sect. 4.1) at a total of 48 sites (Fig. S7) revealed up to five horizons per location, with 15 cm average and 5–50 cm individual ribbon width. An average displacement of 15 cm per event at the studied representative NFS sites yields recurrence intervals of ∼ 100–400 years. The lower three horizons at site BFS_N (from bottom to top: 15 ± 1 cm; 11 ± 1 cm; 5.5 ± 1 cm) were used as input parameters for ³⁶Cl-dating-derived earthquake ages with significantly different results (see following Sect. 4.2.2). Rough magnitude calculations after Wells and Coppersmith (1994) were based on the input parameters (ribbon widths and fault lengths) presented in Table S9. As the connection between different NFS sections and NFS genesis are not trivial, we use different presumptions and calculation methods. Derived magnitudes range from M_w≈ 5.3 to 6.5.

4.2.2 Earthquake ages obtained from ribbons and ³⁶Cl dating

The modeling of the ³⁶Cl concentrations of the BFS_N free face highlights the fact that the measured ³⁶Cl pattern can be generated by a constant slip rate of 1.5 ± 0.1 mm yr⁻¹ (Fig. 7a). Since all five samples are aligned pretty well to fit this synthetical slip rate of ∼ 1.5 mm yr⁻¹, we feel that it is robust. Tests using different amounts of co-seismic offsets did not reveal significant deviations. The given uncertainty does not include external factors like the ∼ 10 % production rate uncertainty, which affects all ³⁶Cl samples in the same way. However, it covers the 68 % confidence interval for at least three of the five ³⁶Cl samples (Fig. 7a) and can hence be considered 1σ internal reliability. Based on this uncertainty we did the further calculations of free-face and fault scarp ages as well as earthquake ages. The retrieved slip rate suggests that the 8.8 m high free face was most likely exhumed within the last 5.9 ± 0.4 kyr (Fig. 7b), and the according fault scarp age is presumably 14.8 ± 1.0 kyr (Fig. 7b). The earthquake ages were estimated by injecting the mapped co-seismic slips of the ribbons around site BFS_N (15 ± 1 cm for EQ1, 11 ± 1 cm for EQ2 and 5.5 ± 1 cm for EQ3) in the ³⁶Cl model (this method is described in detail in Mechernich et al., 2018). The resulting earthquake ages are 100 ± 14 years (EQ1), 173 ± 24 years (EQ2), and 210 ± 29 years (EQ3; Fig. 8). These ages also consider the uncertainties of the ³⁶Cl production rates. Hence, the earthquake recurrence interval appears to be around 35–100 years. Since the slip rate of 1.5 ± 0.1 mm yr⁻¹ is higher than expected, we tested how to receive the lowest possible slip rate explaining the ³⁶Cl data. An erosion of the fault plane itself can be ruled out, since slickensides are visible at least within the lower ∼ 3 m of the free face, and also previous studies report that such faults indicate very low erosion rates (e.g., Goldberg et al., 2016). To be as open-minded as possible, we used any hypothetical scenario. The lowest slip rate on the sampled free face is achieved when minimizing its inherited component of ³⁶Cl, which is generated in the subsurface. The minimum amount in inheritance is caused by a very fast exhumation of the upper scarp, e.g., by a sliding event, which exhumes 14.5 m very fast (1 year in the modeling code; Fig. 7c and d). Before this hypothetical sliding event, the hillslope was exposed for an arbitrary 200 years (from 6.7 ± 0.6 to 6.5 ± 0.6 kyr; Fig. 7d). After the hypothetical sliding event, the stick-slip modeling as described above revealed the best fitting for a slip rate of 1.15 ± 0.10 mm yr⁻¹ during the last 6.5 ± 0.6 kyr (Fig. 7c). Hence, we highlight the fact that the slip rate of the BFS during the last ∼ 6 kyr was surely higher than ∼ 1.15 mm yr⁻¹ and presumably around 1.5 ± 0.1 mm yr⁻¹. The reliability and uncertainties are evaluated in the discussion in Sect. 5.2.

Figure 8 The exhumation history of the free face at site BFS_N based on the modeling results (see Fig. 7a and b). The slip rate of 1.7 ± 0.1 mm yr⁻¹ together with the co-seismic amount of offset based on the mapped earthquake horizons results in earthquake ages of 100 ± 14 years (EQ1), 173 ± 24 years (EQ2), and 210 ± 29 years (EQ3). All values are given within the 68 % (1σ) confidence interval.

5 Discussion

5.1 Interpretation of ribbons, surface rupture lengths, and earthquake magnitude estimation

The described ribbons on fault planes (see Sect. 4.1) are often correlatable across several locations based on their habitus and constant widths, while commonly showing sharp boundaries (Figs. 3b, S6 and S7a–c). We interpret these characteristics to exclude any kind of gradual or localized exhumation by erosion or gravitational processes (see also Sect. 5.4). For example, erosion from human or animal activity (despite being very unlikely in such remote, steep, and overgrown terrain) would not yield such constant ribbon widths over long distances, nor would it account for repeated centimeter-scale exhumation in relatively large time steps (as suggested by the different states of surficial weathering, see Sect. 4.1). Weather- and/or climate-related (i.e., precipitation-related) exhumation would show more gradual transitions between individual ribbons and/or create inconstant ribbon widths due to (topography-related) non-uniform surface runoff. Snow – as the only conceivable weather-related factor – is ruled out, since it is an extremely short-lived phenomenon in this coastal climatic setting at low elevation. Aside from local snow's non-existent capability of leaving distinctly visible marks, it would possibly never create such uniform ribbon widths, as snowdrifts would certainly yield variable thicknesses of the snow blanket. Gravitational sliding of material at the scarp base may occasionally occur. However, we regard this as a phenomenon restricted to particularly susceptible (“exposed”, steep) locations along the fault scarp. Therefore, such processes would equally not create uniform ribbon widths over long distances, but rather deface existing ones. Based on this argumentation, the attribution of the observed ribbons to individual seismic events is highly ascertained. This underlines the interpretation of the ribbons as earthquake horizons (e.g., Mechernich et al., 2018) and technically qualifies them as input parameters for magnitude calculations after Wells and Coppersmith (1994). Minor error sources are misinterpretations of displacement per event, as ribbons may be defaced and overseen. When using fault lengths as input data for the Wells and Coppersmith (1994) method, incorrect recognition of the actual fault lengths constitutes a similar minor error source (see also Sect. 4.2.1; Table S9). The most severe error source, however, is the application of the empirical approach itself. For our setting with short fault lengths and relatively low magnitudes (i.e., M_w<6), Wells and Coppersmith (1994) only present limited data, so the adequacy of this statistical method is questionable. Depending on the interpretation of our NFSs in a per se contractional setting (see Sect. 5.4), they may portray magnitudes ranging ∼ 1 M below the actual magnitude evoked by the rupture of first-order thrust faults – a phenomenon that can, e.g., be observed on Crete, where M_w≈ 8 (or higher) uplifted the western part, whereas onshore normal faults are much shorter but also seismogenic with M_w≈ 6 ± 0.5 (e.g., Grützner et al., 2016; Mason et al., 2016; Schneiderwind et al., 2017; Mechernich et al., 2018, 2022). The obtained values would therefore advance to magnitudes in the range of the Montenegro 1979 earthquake, which are likelier to produce crustal ruptures of such scale (e.g., McCalpin, 2009).

5.2 Reliability of slip rates and ages from ³⁶Cl dating (site BFS_N)

The applied forward modeling method accounts only for the analytical ³⁶Cl uncertainties and not for the uncertainties of the parameters introduced in Sect. 3.2.2 and Table S13. Changes in these input parameters would shift the modeled earthquake ages to older or younger values, without changing the relative recurrence interval (e.g., Mechernich et al., 2018). The largest effect of such a parameter change is related to the ³⁶Cl production rates from Ca spallation or muon capture. A change in these two rates in the framework of published uncertainties would systematically shift all ages and slip rates within ∼ 10 %. This shift is included in the age calculations but not in the slip rate calculations. Furthermore, the estimated parameters for the density of the colluvium, the erosion rate, and the apparent pre-exposure duration can cause similar shifts of the calculated ages. Changes in the erosion rate, e.g., using the minimum erosion rate of 0 mm kyr⁻¹, would result in 3 % younger ages at the top of the free face compared to the 1 mm kyr⁻¹ used, which was chosen based on the 2–8 mm of relief at the top of the free face. Due to the large degraded part of the fault scarp, the choice of the apparent pre-exposure duration has no impact on the restored slip history of the free face. According to the mapping results of the ribbons, we suppose a stick-slip behavior of the NFS as the most likely scenario for the interpretation of the ³⁶Cl data. The few data points would indeed leave a margin for other scenarios, e.g., a landslide or rockfall that exhumed the degraded part of the scarp ∼ 6.5 kyr ago, followed by free-face exhumation (Fig. 7c and d). However, this is ruled out by the fact that no indicators of landsliding were found in the hanging wall at all. Several more complicated scenarios are possible but they would require a larger number of samples and could be the topic of future studies. This study aimed to estimate the slip rate of the free face, which is likely to be 1.5 ± 0.1 mm yr⁻¹ owing to the clearly increasing ³⁶Cl concentrations with scarp height and the mapping of the ribbons and surrounding geology. Also, the very young past earthquake ages of 100 ± 14 years (EQ1), 173 ± 24 years (EQ2), and 210 ± 29 years (EQ3) are quite robust according to the combination of ³⁶Cl data with the mapping of exposure duration. A comparison of the obtained ages with records of (historical) strong earthquakes yields well-fitting and reasonable matches for the two younger proposed earthquakes: EQ1 (100 ± 14 years) may either correspond to the Shkodra 1905 (M_w≈ 6.6; e.g., Koçiaj and Sulstarova, 1980; Aliaj et al., 2010) or the Durrës 1926 (M_w≈ 6.2; e.g., Sulstarova et al., 1980; Aliaj et al., 2010) event. EQ2 (173 ± 24 years) is likely related to the 1855 earthquakes near Shkodra (up to 8^∘ on the MSK scale; Aliaj et al., 2010) or alternatively the 1869 (M_w≈ 6.2) $/$ 1870 (M_w≈ 6.7) events near Durrës (Papazachos and Papazachou, 2003). Concerning the upper part of the fault scarp, we have a lack of data. An extrapolation of the ∼ 1.5 mm yr⁻¹ of slip on the free face would result in a fault scarp age of 14.8 ± 1.0 kyr. This is a reasonable age also found for several fault scarps in the high altitudes of the Apennines (e.g., Cowie et al., 2017; Beck et al., 2018) and corroborated by our LGM-normalized calculations (see following Sect. 5.3). The estimated ∼ 200 m offset presented in Sect. 4.1 indicates incipient activity of the introduced normal faults. Extrapolation of the calculated slip rate of ∼ 1.5 mm yr⁻¹ and using a maximum offset of 200 m results in a potential initiation time of ∼ 135 kyr ago. It is stressed that this estimation does not take possible quiescence intervals into account in the long term. As an absolute maximum age for fault initiation, we define the early Oligocene (∼ 30 Myr). By that time, the Budva–Cukali Unit was conclusively emplaced on top of the Kruja–Dalmatian Unit – a milestone that certainly pre-dates NFS initiation, as BFS and KFS demonstrably overprint the Budva–Cukali Unit and cross-cut its basal thrust into the Kruja–Dalmatian Unit (see also Sect. 2).

5.3 Reliability of slip rates obtained from fault scarp profiling

As described in previous sections, ³⁶Cl dating was conducted for one site (BFS_N) only. To enable a comparison of the different (structurally and exposure-related distinct) sections of the fault scarps nonetheless – and to provide at least one benchmark for the obtained ³⁶Cl dating results – we invoke the rather simplistic technique of fault scarp profiling (see also Sects. 3.1.1 and 4.2.1) for slip rate derivation. The applied slip rate estimation relating free-face and NFS heights to a proposed post-LGM period of ∼ 15 ± 3 kyr holds three main error sources: (i) the exact timing of initiating fault scarp preservation – in the case that it was effectively impeded during glaciation, (ii) the local impact of LGM climate on erosion at all, and (iii) the interpretation of NFS geometry. Being aware that the exact timing of the (local) LGM is a matter of debate, we synthesized 15 ± 3 kyr as an adequate time frame for our location (based on available data and literature from surrounding areas in the Balkans, Greece, and Italy; e.g., Giraudi, 1995; Giraudi and Frezzotti, 1995, 1997; Allen et al., 1999; Kuhlemann et al., 2009; Papanikolaou et al., 2005, 2013). The large variety of ages presented in the literature (for comparison, e.g., Pope et al., 2017; Pavlopoulos et al., 2018) still causes us to admit a given degree of uncertainty concerning this value. The estimation of errors connected to (ii) and (iii) is equally non-trivial and related to both tectonic and climatic–erosional impacts. Even if the above-specified LGM timing is accurate for a broader region, small-scale local variations may exist, e.g., as a result of a location's elevation, exposition, and microclimate. When interpreting the NFS geometries at BFS and KFS, the transient zone between the characteristic free face and regional hillslope – the so-called “degraded scarp” – is striking. It forms part of the full NFS height and is constructed by an interpolation of free face and regional hillslope (Fig. 5). The full NFS height is typically used as an input value to calculate generalized post-LGM slip rates (e.g., Papanikolaou et al., 2005). Our studies based on ³⁶Cl dating suggest that the highest part of the BFS_N free face was exhumed ∼ 6 kyr ago, yielding slip rates of ∼ 1.5 ± 0.1 mm yr⁻¹. For the same site, our LGM-normalized and profile-based calculation yields an annual slip of 1.46 ± 0.4 mm along the full NFS height (Table S10) of 22.2 m. Since these values are strikingly well-fitting, we are confident that the hypothesis of a post-LGM formation of both the free face and degraded scarp as well as the assumed timing of initiated fault scarp formation (i.e., ∼ 15 kyr) are appropriate. Based on our results from the only ³⁶Cl-dated site BFS_N, we are largely convinced that our profiling-derived slip rates for BFS_S and KFS (Sect. 4.2.1; Table S10) are equally reliable. Since the general setup of these nearby NFSs is similar (in terms of ribbon widths, general NFS heights, and morphologies), all necessary preconditions are given. When eventually looking at the obtained results from KFS and particularly BFS_S, it is remarkable that slip rates are significantly lower and degraded scarps are generally less developed. Both observations are likely connected to the BFS_S and KFS sections' orientation, which significantly differs from that of BFS_N: (i) the orientation of BFS_S and KFS with respect to underlying first-order structures (“Dinaric strike”) could lead to a different (slower) slip behavior. While BFS_N best represents the overall Dinaric strike, BFS_S and KFS reproduce the strike of the majority of neighboring geological and geomorphic features (e.g., the southerly adjacent anticlines). (ii) The KFS and BFS_S are south-facing (i.e., not immediately sea-facing) and better protected by vegetation, which likely leads to a better preservation of intact free faces (and less-developed degraded scarps). This hypothesis is substantiated by the fact that, especially for BFS_S, the dissection by erosional gullies is minimal. Free faces within this section are steep and moderately high, with clearly visible and abundant earthquake ribbons (Fig. S7b). By contrast, BFS_N is more exposed to weathering (i.e., sea-facing, surrounded by less vegetation) and dissected by numerous gullies (Figs. 3 and S1). Here, more degraded scarps, high, shallowly dipping free faces, and fewer earthquake ribbons are observed (Fig. S7a). In a similar approach, the transition from well-defined free faces to degraded scarps could possibly be explained by changing climatic conditions favoring free-face preservation staring ∼ 6 kyr ago (e.g., Lambeck et al., 2011; Boulton and Stewart, 2015). Provided that the introduced quantitative NFS profiling method may be rather error-prone (see discussion above) and low-resolution, we again stress that it was used as an auxiliary tool complementing the ³⁶Cl dating method (see Sects. 3.1.2, 3.2.1, 3.2.2, 4.2.2, and 5.2).

Figure 9 (a) Apulia referenced GPS velocity field and (b) interpolated strain rate with a Gaussian–Voronoi cell weighting of a net reweighting threshold of W_t = 6 in the southwestern Balkans. (c) Swath topographic section with GPS velocity information through the working area. Panels (a)–(c) are modified after and reprinted with permission from D'Agostino et al. (2020); see their work for details. (d) Normal-faulting-related horizontal velocity component and its accelerating role in cross-regional convergence.

Figure 10 3D block diagram of the tectono-morphological features of the Montenegrin–Albanian coastal border region. Forward-modeled structural cross section through the external Dalmatian nappe stack, modified after Schmitz et al. (2020). The spatial proximity of the extensional NFS and contractional dry valleys as well as a supposedly tectonically deviated stream (Biermanns et al., 2019; Schmitz et al., 2020) support the conclusion of highly interactive tectonic regimes.

5.4 Formation mechanisms of the BFS and KFS normal fault scarps

The position of the BFS and KFS in the hinge of thrust-related anticlines within the nappe stacks of the Budva–Cukali and High Karst units suggests their possible origin along pre-existing planes of weakness (fold-related longitudinal fractures, e.g., Ramsay and Huber, 1987; Tavani et al., 2015, and references therein). Several indicators strongly suggest that the NFSs were formed within the framework of regional contractional tectonics and repeatedly ruptured as second-order structures during earthquakes released along first-order active basal thrust faults, where strain is partitioned in the upper plate or hanging wall (e.g., Hicks and Rietbrock, 2015). In this context, a possible explanation for initial NFS formation could be the underlying setup of ramp and flat geometries, e.g., leading to fold crest collapse as a result of material slipping beyond the upper end of a ramp segment of the detachment (comparable to the El Asnam 1980 earthquake; e.g., Avouac et al., 1992). A significantly less probable explanation is NFS formation as a direct effect of northwestward-migrated hinterland extensional tectonics. The boundary between foreland contraction and hinterland extension has evidently propagated westward since the late Eocene (e.g., Dumurdzanov et al., 2005; Reicherter et al., 2011; Handy et al., 2019). Recent geodetic studies show that the working area lies in the frontal-most part of the deformation zone, right at the tip of a northwestward-shifting line separating hinterland extension from foreland contraction (Figs. 1 and 9b; D'Agostino et al., 2021). In the hypothesis of a migrated extensional hinterland domain, the existence of the described NFSs would have implications for observed GPS-derived convergence rates: the fault slip rate of 1.5 mm yr⁻¹ along the normal fault plane with a dip of 60^∘ results in a horizontal extension of 0.75 mm yr⁻¹. In order to achieve the geodetically observed convergence of 3–4 mm yr⁻¹, the actual convergence must be 3.75–4.75 mm yr⁻¹ (Fig. 9d). Despite the proven proximity of a hinterland extensional domain, the number of arguments in favor of a still predominantly contractional regime is significant: (i) the local seismicity is characterized by exclusive occurrence of instrumentally recorded contractional earthquakes and an entire lack of extensional ones. In combination with geodetic data, this rather excludes the advance of hinterland extension into the debated domain. Such a scenario would require ceased activity of the basal orogenic thrust owing to dissection by the obliquely oriented normal faults, which is rather unlikely. Although the downdip depths of BFS and KFS cannot be precisely quantified, it is rather unlikely that they cross-cut the active basal thrust at depth. (ii) The accommodation of total convergence in the coastal area (Fig. 9c) is further substantiated by the existence of other recent geomorphological features in the study area, such as dry valleys and deflected river channels (Fig. 10; Biermanns et al., 2019; Schmitz et al., 2020). (iii) A full migration of hinterland extension into the study area would likelier – or at least additionally – create inland-facing fault scarps. Although graben-like characteristics have been recognized for nearby areas (such as the Shkodra Lake basin, Šasko Lake valley, Zadrima or Balldreni plain; e.g., Sulstarova et al., 1980; Aliaj et al., 2001; Aliaj, 2006; Uncu, 2012; Biermanns et al., 2019, and references therein), their geomorphic expression does not involve prominent fault scarps or other unequivocal evidence of large-scale normal faulting. The formation of BFS and KFS as a consequence of gravitational collapse, i.e., landsliding, can be ruled out. (i) To our knowledge, it has never been reported or proven that landslides create geomorphic landscapes with large-scale surface ruptures and/or fault planes in a similarly complex array (i.e., largely consistent overall characteristics of the rupture, despite changing orientation, bedrock, and several apparent outcrop gaps). Instead, exactly these characteristics apply to other verified examples of seismogenic NFSs (e.g., Mason et al., 2016, 2017, a study that is probably best comparable in terms of bedrock, fault scarp, and slope morphology). (ii) The same is valid for the slickenside striations on ∼ 10 m high(!) free faces that – although slightly diverging outward – are almost consistent with respect to the undulating but otherwise planar fault planes. (iii) No convincing (geomorphic) features typical of landsliding have been identified in the surroundings of the NFSs (for comparison, e.g., Highland and Bobrowski, 2008). Large portions of the hanging wall, particularly along the eastern ends of BFS_S and KFS, are clearly lacking any potential landslide deposits (i.e., Holocene colluvial deposits; Fig. 4a). (iv) Neither the described earthquake-related ribbons nor the mapped transfer faults or the western termination of KFS (indicative of strike-slip motion; see end of Sect. 4.1) are compatible with a landslide-related origin. (v) The described transition of BFS_N into BFS_S features fault traces in bare rock, forming relays. (vi) The presented setting is assessed to be scarcely landslide-prone with its barely water-saturated massive limestone bedrock and unhampered runoff (as particularly underpinned by the pronounced gullies across the NFS).

6 Conclusions

We report two previously unknown, active normal faults with well-preserved bedrock NFSs along the contractional front of the southern Dinarides fold-and-thrust belt. We propose a tectonic, co-seismic origin of these structures. Relations between fault orientation, striations, earthquake ribbons, and surrounding structures suggest that the normal faults are either the result of rollback-induced westward-migrating extensional tectonics or more likely second-order features linked to subjacent, higher-order thrusts capable of triggering earthquakes up to M_w≈ 7 ± 0.5. ³⁶Cl cosmogenic dating for one sampling site on the NFS suggests that the corresponding ∼ 9 m high free face was exhumed within the last ∼ 6 kyr. An annual slip rate of 1.5 ± 0.1 mm is released in 6–15 cm steps during earthquakes with recurrence intervals estimated to 35–100 years. To incorporate the exhumation history of the degraded part of the same fault scarp in our model and to allow statements about three additional NFS sites that were not dated by means of the ³⁶Cl method, we calculated sets of qualitative slip rates based on fault scarp profiles. Normalized to a full post-LGM period of 15 ± 3 kyr, this method yields slip rates of 0.5 ± 0.1 to 1.5 ± 0.4 mm yr⁻¹ and recurrence intervals for major earthquakes in the range of 100–400 years – values that compare well with the ³⁶Cl-derived ones. Altogether, all of the presented values appear realistic against the backdrop of available GPS rates and common earthquake magnitudes in the region. The normal faults are exactly located above the subsurface thrust fault and epicenter that was responsible for the M_w 7.1 Montenegro 1979 earthquake and hence suggest a relation. In any case, we regard the NFSs as a manifestation of repeated earthquake activity in the study area.