The size distributions of fractures and earthquakes: implications for orogen-internal seismogenic deformation

Pre-existing geological discontinuities such as faults and fractures represent structural and mechanical discontinuities in rocks which influence earthquake processes. As earthquakes occur in the subsurface, seismogenic reactivation of pre-existing fracture networks is difficult to investigate in natural settings. However, it is well-known that there exists a physical link between both fractures and earthquakes since an earthquake's magnitude is related to the ruptured fracture area and therefore fracture length. Furthermore, fractures and earthquakes exhibit similar statistical properties, as their size distributions follow power laws.

In this study, we exploit the relation between the size (or length) distributions of pre-existing fractures and earthquakes to decipher the seismic deformation processes within the exhumation-related orogen-internal setting of the southwestern Swiss Alps, which due to its well-monitored seismic activity and the excellent outcrop conditions provides an ideal study site. Characterizing exhumed fracture networks from different tectonic units based on multi-scale drone-based mapping, we find that power law exponents of 3D fracture size distributions generally range between 3 and 3.6. Comparing these values with the depth-dependent exponents of estimated earthquake rupture lengths, we observe significantly larger values of 5 to 8 for earthquake ruptures at shallow depths (< 3 km below sea level (b.s.l.)). At intermediate crustal depths (∼ 3 to 9 km b.s.l.), the power law exponents of fractures and earthquakes appear to be similar. These findings imply depth-dependent differences in the seismogenic reactivation of pre-existing fractures in the study region: while partial rupturing is the prevailing deformation mechanism at shallow depths, fractures are more likely to rupture along their entire length at intermediate crustal depths. Therefore, the present-day near-surface differential stresses are likely insufficient to rupture entire pre-existing fractures seismogenically. Our findings have direct implications for seismic hazard considerations, as earthquakes that rupture along entire fractures appear to become less likely with decreasing depth.

1 Introduction

Within the Earth's crust, the occurrence of earthquakes is commonly attributed to the reactivation of pre-existing geological discontinuities, commonly referred to as fractures. Nevertheless, in regions characterized by extensive fracturing and scattered seismic activity, like in collisional orogens such as the Alps, the task of assigning earthquakes to specific fracture segments becomes notably complex. This challenge limits our ability to decipher the current deformation processes and discern the significance of such intricate pre-existing fracture networks in influencing the underlying seismic events.

Despite their complexity, fractures in nature often appear qualitatively self-similar at different scales. This qualitative self-similarity can be quantified by power law distributions, which inherently implies scale invariance and thus the absence of a characteristic length scale. This means that, strictly speaking, fracture segments do not have a unique length but that a measure of the length of a fracture is always dependent on the scale of observation. Previous studies have shown that a power law is best suited for describing the size (or length) distribution of natural fracture traces (Odling, 1997; Bonnet et al., 2001; Torabi and Berg, 2011; Yielding et al., 1996; Davy et al., 2010; Scholz, 2019), even if other types such as exponential or log-normal distributions have been proposed (Cowie et al., 1995, 1993; Ackermann et al., 2001; Nicol et al., 1996). A power law distribution defines the number n of features of a given length l as

$$\begin{array}{}\text{(1)}& n\left(l\right)=c{l}^{-\mathit{\alpha }},\end{array}$$

where α represents the power law exponent (i.e., the slope of the linear trend in a log–log plot), and c denotes the density constant. The power law exponent α is of main interest herein, since it measures the relative proportion of short and long features (Davy, 1993; Pickering et al., 1995). As fracture data are often acquired at a single mapping scale, a major problem for the identification of power law distributions is the limited dynamic length range of the fracture traces and the uncertainties related to sampling bias, which can lead to size distributions that appear to deviate from power laws (Odling, 1997; Scholz, 2007, 2019; Torabi and Berg, 2011).

Similar to fractures, it is well-known that earthquake magnitudes follow power law distributions, often referred to as the Gutenberg–Richter law (Gutenberg and Richter, 1944) or magnitude–frequency distribution, implying that also earthquakes are scale invariant. As proposed by Hatton et al. (1993), “the observed fractal nature of both fracture length distributions and earthquake magnitude–frequency distributions suggest that there may be a relationship between the structure of active fracture systems and the resulting seismicity.” This is supported by the fact that the statistical properties of micro-fracturing at the laboratory scale are similar to those of earthquakes (e.g., Scholz, 1968; Mogi, 1962). Furthermore, fractures represent the geometrical anisotropies that provide planes of weakness within the rock masses and are thus most likely to be reactivated seismically if suitably oriented with respect to the local stress field. Fractures and earthquakes can thus be regarded as the long- and short-timescale phenomena of brittle tectonics (Scholz, 2019), and there should exist an inherent relation between the size distributions of fractures and earthquakes (Scholz, 1998; Turcotte, 1997). However, the relationship between these size distributions in natural datasets has so far received little attention, mainly due to the challenging acquisition of reliable quantitative information of fracture networks due to limited outcrop conditions and small dynamic range. The southwestern Swiss Alps, with (i) a long-lasting exhumation history of about 20×10⁶ years (e.g., Egli et al., 2017; Herwegh et al., 2020; Boutoux et al., 2016; Cardello et al., 2024), (ii) well-preserved extensive Neogene faulting in high-Alpine regions with formation depths similar to the present-day earthquakes (Cardello and Mancktelow, 2015; Huggenberger and Aebli, 1989; Cardello et al., 2024), (iii) an increased and exceptionally well-recorded seismic activity (e.g., Diehl et al., 2021a; Lee et al., 2023), and (iv) ongoing vertical tectonics (e.g., Ustaszewski and Pfiffner, 2008; Brockmann et al., 2012; Herwegh et al., 2024; Piña-Valdés et al., 2022) implying a similar stress field to the time of fracture formation, provide a unique natural laboratory to investigate the link between the size distributions of fractures and earthquakes. Owing to the abundant occurrence of thermal springs, the region is also a potential target for the exploitation of geothermal energy, where knowledge of faulting at depth and constraints on the potential risks of induced seismicity are critical.

Herein, we explore the relation between the size distributions of fractures (fracture length) and earthquakes (rupture length) and its implications for orogen-internal seismogenic deformation in the exhumation-related (vertical) tectonic setting of the Rawil depression region (southwestern Swiss Alps) (Fig. 1). We use field observations and drone-based fracture trace mapping on multiple scales to quantitatively characterize the size and orientation distributions of fracture networks in different tectonic units. We employ a multi-scale power law fitting approach on the fracture network data, which allows us to expand the limited dynamic length range of single-scale observations. We then exploit the relationship of the fracture orientation and size distributions with similar properties of the recent seismic activity and discuss the implications for seismogenic fracture reactivation processes in the subsurface.

Figure 1 Regional setting of the investigated seismotectonic domain z21 of Wiemer et al. (2016) (black polygon), with the four fracture mapping locations (A–D) shown as circles. (a) Tectonic map (1:500 000; © swisstopo). RSF: Rhône-Simplon Fault Zone. (b) Earthquakes from 1984 to 2020 (Diehl et al., 2021b) and available focal mechanism data (see main text for references). Blue arrows show the recent crustal stress field orientations after Kastrup et al. (2004), with contrasting stress regimes north and south of the RSF. The shaded colors indicate the different subdomains A–D. RFZ: conceptual Rawil Fault Zone. (c) Recent crustal displacement rates based on GNSS surface velocities in the vicinity of the study area with reference to stable Europe (data: EUREF Permanent GNSS Network (Bruyninx et al., 2019); swisstopo). Digital elevation model: EU-DEM v1.1 (© European Union, Copernicus Land Monitoring Service 2023, European Environment Agency (EEA)).

2 Regional setting

In the framework of a nation-wide seismic hazard assessment, Switzerland was zoned into different seismotectonic domains (Wiemer et al., 2016). We herein focus on seismotectonic domain z21 in the southwestern Swiss Alps, mainly due to its well-exposed rock outcrops and the exceptionally well-monitored natural seismic activity (Fig. 1). In the study area, several external crystalline massifs (ECMs) are exposed: the Aiguilles Rouges massif in the west and the Aar and Gastern massifs in the east (Fig. 1a). Between their exhumed parts, these basement units form a saddle-like structure, commonly referred to as the Rawil depression, which is overthrusted by the Helvetic limestone nappes (Burkhard, 1988; Dietrich, 1989; Dietrich and Casey, 1989; Ramsay, 1989, 1981). The Helvetic nappe system is composed of different, heavily folded nappes which have been stacked on top of each other during an early thin-skinned deformation phase of the Alpine orogen (thrusts in Fig. 1a; e.g., Musso Piantelli et al., 2022; Escher et al., 1993; Pfiffner 1993; Burkhard, 1988; Ramsay, 1981, 1989). In a late-Alpine stage, dominated by thick-skinned vertical tectonics, both the crystalline basement units and the Helvetic limestones were heavily affected by late-Alpine deformation, resulting in pervasive and complex faulting that affects all tectonic units and crosscut the older thrusts (faults in Fig. 1a; Cardello and Mancktelow, 2015; Gasser and Mancktelow, 2010; Huggenberger and Aebli, 1989; Herwegh et al., 2020, 2024; Cardello et al., 2024). Some of the predominantly subvertical late-Alpine faults have been linked to the regional seismic activity (Cardello and Mancktelow, 2015; Pavoni, 1980a; Ustaszewski et al., 2007; Ustaszewski and Pfiffner, 2008; Pavoni and Mayer-Rosa, 1978; Cardello et al., 2024). With the exception of the Rhône-Simplon Fault Zone (RSF), a steeply south-dipping dextral strike-slip fault system along the southern margin of the study area (e.g., Campani et al., 2014, 2010; Egli and Mancktelow, 2013; Mancktelow, 1985), regional-scale faults are largely absent at the surface of the study area. The nappes of the Helvetic units are separated by shallow-dipping thrust faults which are arguably not seismically active nowadays (Lee et al., 2023), so we herein do not discuss these thrusts further.

The presently dense seismic monitoring network in the region resolves the seismic activity exceptionally well, detecting earthquakes with a magnitude of completeness M_C of about M_L≈ 1.5 until 2016 and about M_L≈ 1.0 thereafter (Diehl et al., 2021a; Lee et al., 2023). The earthquakes predominantly occur in a NE–SW-striking, elongated corridor in the center of the Rawil depression (Fig. 1b), which is referred to as the conceptual Rawil Fault Zone (RFZ) (Lee et al., 2023). The seismicity at depth arguably reactivates a complex system of subvertically oriented faults (Diehl et al., 2021a, 2018; Lee et al., 2023; Truttmann et al., 2023). This is supported by the variety of previously published, predominantly strike-slip-type focal mechanisms shown in Fig. 1b (Baer et al., 2005, 2003, 1997; Deichmann et al., 2012, 2006, 2002, 2000; Delacou et al., 2005; Diehl et al., 2021a, 2018, 2013; Jimenez and Pavoni, 1983; Maurer, 1993; Maurer and Deichmann, 1995; Pavoni et al., 1997; Truttmann et al., 2023), suggesting an overall transtensional regime (Kastrup et al., 2004; Houlié et al., 2018; Maurer et al., 1997).

The southwestern Swiss Alps are further characterized by their long-lasting vertical tectonics. During the latest stage of the Alpine orogeny, thick-skinned tectonics prevailed and led to the exhumation of the ECMs associated with intense faulting (faults in Fig. 1a; e.g., Herwegh et al., 2020, 2024; Boutoux et al., 2016). Recent uplift rates in the range of 1 mm per year (Brockmann et al., 2012; Piña-Valdés et al., 2022) and the pronounced seismic activity document the ongoing exhumation (Fig. 1b and c). As a consequence of these long-lasting vertical tectonics, a similar stress field likely prevailed during both the formation of the exhumed fracture networks and today's earthquakes.

In summary, the investigated seismotectonic domain reveals complex and pervasive fracture networks, well-recorded seismic activity, and long-lasting vertical tectonics. In our view, the Rawil depression region in the southwestern Swiss Alps thus provides a unique natural laboratory to elaborate on the statistical link between fractures and earthquakes in an orogen-internal setting.

3 Data and methods

3.1 Fracture networks

3.1.1 Fracture data

We selected four sites along the Rawil depression to conduct a detailed analysis of the fracture networks in the ECM units of (A) the Aiguilles Rouges massif and (D) the Gastern massif, as well as in the Helvetic limestones of (B) the Diablerets and (C) the Wildhorn nappes (Fig. 1). As pronounced glacial and fluvial erosion, vegetation, and Quaternary cover greatly limit the outcrop conditions and thus affect reliable sampling of fracture networks, we selected sites with little weathering, vegetation, and Quaternary cover, as well as flat terrain, to minimize topographic effects of fracture lengths and directions (Baumberger et al., 2022), which were encountered only at high-altitude sites.

At each site, we performed a detailed characterization of the fracture networks (Fig. 2a). We carried out drone-based surveys (DJI Mavic 2 Pro; Hasselblad L1D-20c camera) at two different observation heights to produce orthorectified images (orthophotos) and digital elevation models (DEMs) with Agisoft Metashape (v1.6.6) at millimeter and centimeter ground resolutions (see Table 1). Additionally, we used orthophoto and DEM products at DM ground resolution from the Federal Office of Topography swisstopo. We then manually digitized the visible fracture traces at fixed mapping scales of 1:10, 1:100, and 1:1000 within circular mapping areas, preferring circular extents to minimize potential orientation bias (e.g., Bonnet et al., 2001) (Table 1, Fig. 2a). Fracture traces are herein defined as curvilinear features that represent the intersections between a fracture plane and the Earth's surface (Baumberger et al., 2022; O'Leary et al., 1976), which serve as a proxy for the pre-existing geometrical anisotropies. Because such pre-existing anisotropies are likely to be reactivated also in today's stress field if suitably oriented and because the power law exponents are insensitive to fracture type (Bonnet et al., 2001; Marrett et al., 1999), we chose not to distinguish between fracture types in our analysis. Unambiguous differentiation of fractures from other discontinuities such as bedding structures based on remote sensing alone is often ambiguous, which is why we systematically validated the fracture maps in the field to ensure their representativeness of the real fracture populations. We also measured both the 3D orientations of the fractures and where possible the kinematic indicators (slickenfibres, slickolites) in the field.

Table 1 Overview of the orthoimages and digital elevation models (DEM) used that serve as a basis for the fracture trace maps. Orthoimages and DEMs of the 1:1000 mapping scale (DM resolution) are publicly available from swisstopo.

Figure 2 (a) Workflow of the fracture network analysis conducted separately for each location, combining observations from three different mapping scales of 1:1000, 1:100, and 1:10 (for details see “Data and methods”). (b) Examples of size distributions of the fracture networks, plotted with the non-normalized single-scale approach, treating the data from each mapping scale individually (upper panel), and the multi-scale approach, combining all information across scales by normalization with the fractal area (lower panel).

Despite thorough data collection, the mapping of fracture traces includes several limitations. The detectability of fracture traces by remote sensing depends on both the fracture length, the image resolution, the mapping scale, and the degree of erosion of the related tectonites (i.e., the resulting morphological incision). Furthermore, image artifacts caused by illumination or shadow effects in low-contrast domains or by partial coverage of the bedrock by Quaternary debris or vegetation can prohibit the appropriate identification of fractures. These factors lead to deviations from the true fracture trace length which in turn influences the quantitative analysis (e.g., Cao and Lei, 2018). To minimize such potential artifacts, we carefully selected sites with little Quaternary cover and topographic effects and acquired the images in similar light conditions. Areas covered by debris or vegetation were excluded as “no data” areas. The influence of the remaining uncertainties is further discussed in Sect. 5.1.

3.1.2 Statistical analysis of fracture networks

We used the data from the fracture trace maps to quantify both the orientation and size distributions of the fracture networks. To describe their 2D orientation distribution, we derived length-weighted rose diagrams from the fracture traces for each individual dataset presented in Table 1 (Fig. 2a). Additionally, we used field measurements to characterize the 3D orientation distribution of the fracture networks, as well as their kinematics. The 2D and 3D distributions capture different aspects of the fracture orientations: while the 3D distribution yields information about the dip range of the fracture sets, the 2D orientation distribution from the fracture trace maps represents a better-suited quantitative measure of the dominant strike orientations.

In a next step, we evaluated the size distribution of the fracture networks from the four sites. We assume an underlying power law distribution, which is most commonly used to model fracture size distributions (e.g., Bonnet et al., 2001; Bour et al., 2002; Odling, 1997; Scholz et al., 1993; Torabi and Berg, 2011; Yielding et al., 1996; Davy et al., 2010; Scholz, 2019). In nature, all power laws must have upper and lower limits (Bonnet et al., 2001; Torabi and Berg, 2011; Turcotte, 1997, 1989). For fractures, the upper limit is likely related to the thickness of the crust or the stratigraphic layering, while the lower limit is constrained by a physical length scale (e.g., grain size). Defining these limits is nontrivial, and deviations from a power law in the tails of the size distribution of fracture networks commonly occur presumably independent of the physical limits of the system. Instead, sampling bias such as truncation and censoring effects can result in size distributions that appear to be exponential or log-normal (Bonnet et al., 2001; Odling, 1997; Odling et al., 1999; Torabi and Berg, 2011). Due to the underestimation of the frequency of small fractures, these structures are systematically underrepresented due to the limited resolution of the orthophotos and DEMs, which is commonly referred to as truncation (Bonnet et al., 2001; Bour et al., 2002; Torabi and Berg, 2011) (Fig. 2b). The truncation length l_trunc defines the length above which fractures can be reliably detected. Similarly, large fractures are more likely to lie partially outside the boundaries of the sampling area, so their lengths and frequencies are often systematically underestimated (Bonnet et al., 2001; Pickering et al., 1995). This leads to a steepening of the curve at the upper tail of the distribution, commonly called the censoring effect (Fig. 2b). Note that the fracture lengths at which truncation and censoring effects occur do not necessarily correspond to the true physical limits of the power law distribution (Bonnet et al., 2001). Since both truncation and censoring effects in natural fracture datasets are often rather pronounced, it is a major challenge to identify the scaling range where the distribution follows a power law. It is therefore crucial to ensure that the power law fitting is based on a large range of fracture lengths, ideally spanning values over 2 to 3 orders of magnitude (Bonnet et al., 2001). Practical issues and mapping bias (see previous section), however, prevent the acquisition of high-resolution orthophotos over large areas that would be required for this. On a single observation scale, the length range is thus limited, rendering the identification of power law distributions challenging (Bonnet et al., 2001; Bour et al., 2002; Davy et al., 2006; Pickering et al., 1995; Scholz, 2007).

To overcome these limitations, we herein used a multi-scale approach, combining the observations from different mapping scales (Fig. 2b). The basic idea is to normalize the fracture density with the fractal area, which allows us to greatly extend the length range (e.g., Bonnet et al., 2001; Heffer and Bevan, 1990; Odling, 1997). The normalized density size distribution is given as (Bour et al., 2002; Davy et al., 1990)

$$\begin{array}{}\text{(2)}& n\left(l,L\right)=c{l}^{-{\mathit{\alpha }}_{\mathrm{F}}}{L}^D,\end{array}$$

with l as the fracture trace length, c representing the fracture density, α_F as the power law exponent of the fracture networks, and L^D denoting the area term. The latter is defined by the diameter of the study area L and the fractal dimension D of the fracture trace barycenters (i.e., the centroids of the fracture traces) and fixes the number of structures of a given length per unit fractal area independently from the observation scale and is thus a scale-independent measure (Davy et al., 1990; Bour et al., 2002). This allows for comparison of the data obtained at different scales.

The fractal dimension D was derived with the two-point correlation function C₂ (Hentschel and Procaccia, 1983; Vicsek, 1992), which describes the spatial correlation of fracture trace barycenters and has been suggested as an appropriate measure for D of fracture networks (Bonnet et al., 2001; Bour et al., 2002; Bour and Davy, 1999; Davy et al., 1990). It is defined as

$$\begin{array}{}\text{(3)}& C_{\mathrm{2}}\left(r\right)={\displaystyle \frac{\mathrm{1}}{N^{\mathrm{2}}}}{N}_d\left(r\right),\end{array}$$

where N is the total number of fracture trace barycenters and N_d is the number of pairs whose barycenter interdistance is less than r (Hentschel and Procaccia, 1983). For fractal populations, C₂ scales with r as r^D, with D as the correlation dimension. Following Bonnet et al. (2001), we estimated the fractal dimension by identifying a plateau of the local slopes of D.

To then fit the fracture trace data with Eq. (2), we only extracted the unbiased data from each fracture trace map. To account for truncation effects, we determined the truncation length l_trunc with the statistical approach of Alstott et al. (2014) and Clauset et al. (2009), using the Kolmogorov–Smirnov (KS) distance as a measure for the goodness of fit between the fitted power law and the data. Hereby, the minimum KS distance defines l_trunc (Fig. S1 in the Supplement). We then removed all fracture traces with lengths smaller than l_trunc. We also accounted for censoring effects by removing fractures that intersect the sampling window (Lei et al., 2015). The cleaned fracture densities of each mapping scale n(l) were then normalized with the area term L^D. We finally fitted the normalized dataset with a power law using maximum likelihood estimation (MLE). This is commonly preferred over linear regression (Clauset et al., 2009; Goldstein et al., 2004). We finally obtained power law size distributions which, owing to the multi-scale approach, span several magnitudes of fracture lengths (Fig. 2b). As the number of fracture traces is critical for the reliable extraction of power law exponents (Zeeb et al., 2013) and the number of fractures is rather small for non-dominant orientations, we chose to not incorporate orientation-dependent α_F values but rather derived a single bulk value of α_F for the respective site.

3.2 Earthquakes

We analyzed the recent seismicity within four subdomains of the seismotectonic domain z21 of Wiemer et al. (2016), each associated with one of the fracture mapping sites A to D (Fig. 1b). In analogy to the fracture data, we characterized both the orientation and size distributions of the earthquakes within these subdomains.

The orientation of the rupture plane activated during an earthquake often remains unknown, and the determination of the orientation distribution of seismic ruptures is thus challenging. However, focal mechanisms (e.g., determined by the polarities of the first arriving P waves) give a first-order constraint on the orientation of the ruptured plane for earthquakes with significant energy release. Such first-motion focal mechanisms yield an ambiguous solution with two possible rupture plane orientations. Additional constraints such as the spatial distribution of hypocenters can in some cases be used to decipher the activated ruptured plane (e.g., Truttmann et al., 2023), but this information is often lacking. For this study, we compiled previously published focal mechanism data as shown in Fig. 1b (Baer et al., 2005, 2003, 1997; Deichmann et al., 2012, 2006, 2002, 2000; Delacou et al., 2005; Diehl et al., 2021a, 2018, 2013; Jimenez and Pavoni, 1983; Maurer, 1993; Maurer and Deichmann, 1995; Pavoni et al., 1997; Truttmann et al., 2023), which partially incorporates information about the active rupture plane deciphered from the spatial distribution of fore- and aftershocks. Based on these focal mechanisms, we evaluated the 3D orientation distribution of earthquake ruptures within each subdomain by plotting the nodal planes in stereographic projections and incorporate, where available, the information about the effectively ruptured plane.

We consecutively analyzed the size distribution of earthquakes, using the data from the earthquake catalog published by Diehl et al. (2021b) only considering natural earthquakes (Fig. 1b). Indications exists that the M_L magnitudes in the SED bulletin include inconsistencies between earlier and modern epochs due to changes in instrumentation and determination procedures of magnitudes (Staudenmaier et al., 2018). We therefore only considered events after 2005, for which the M_L magnitudes are derived from digital data in a relatively consistent manner. Since we are interested in a measure that is directly comparable to the size distribution of fractures, here we derive the size distribution of earthquake ruptures, where size refers to length l_R, defined as the diameter of a circular earthquake rupture, rather than magnitude or seismic moment, which serves as a proxy for the recently activated part of the fracture network at depth. This approach is justified by the fact that the earthquake magnitude is directly related to the ruptured fracture length (e.g., Wells and Coppersmith, 1994). Therefore, we translated the local magnitudes M_L to the earthquake rupture lengths l_R. After converting the recorded local earthquake magnitude M_L to the moment magnitude M_w after Goertz-Allmann et al. (2011) using an empirical M_L−M_W scaling relationship, we derived the earthquake rupture area A_R, given as

$$\begin{array}{}\text{(4)}& A_{\mathrm{R}}={\mathrm{10}}^{\frac{(M_{\mathrm{w}}-a)}{b}},\end{array}$$

where a and b represent constants derived from empirical scaling relations (e.g., Wells and Coppersmith, 1994). We herein used three different scaling relations determined for strike-slip earthquakes in stable continental regions with a= 4.18 and b= 1 (Leonard, 2014), a= 3.49 and b= 0.94 (Thingbaijam et al., 2017), and a= 3.98 and b= 1.02 (Wells and Coppersmith, 1994). From A_R, we then calculated the earthquake rupture length l_R under the assumption of a circular rupture plane. Based on this new dataset containing a rupture length for each earthquake, we then employed the MLE power law fitting procedure of Alstott et al. (2014) and Clauset et al. (2009) using Eq. (1) to derive the power law exponent α_R of the earthquake ruptures. α_R represents a modified b value that is the commonly used power law exponent in seismology. Similar to the fitting procedure of the fracture data, we used the minimum KS distance as a measure for the optimal truncation length (Fig. S2), which is comparable to M_c often used by seismologists.

4 Results

4.1 Fracture networks

Based on remote sensing, we mapped between 1608 and 10 894 fracture traces at each of the three different mapping scales on the four study sites (Table 1, Fig. 3). In the following, we use these data to characterize the fracture networks in terms of their orientation and size distributions.

Figure 3 Fracture trace maps at different mapping scales (rows) for all four study sites (columns). The maps show both fracture traces (black lines) and no data areas (grey areas). The length-weighted rose diagrams to the top right of each map show the 2D orientation distribution of the fracture traces (strike directions). N indicates the number of mapped fracture traces.

4.1.1 Orientation distribution

Site A is characterized by three main sets of subvertical fractures identified in the field (Fig. 4a): a NE–SW-striking set dipping towards the SE (A_I) and a similarly striking but NW-dipping set (A_II), as well as an almost perpendicular NW–SE-striking set (A_III). Set A_I consistently indicates normal movements, while A_II fractures reveal a reverse shear sense (Fig. 4b). For set A_III, we only observed a single kinematic indicator, revealing oblique sinistral movements. The length-weighted rose diagrams for site A suggest that the NW–SE-striking A_III fractures are the dominant set (Fig. 3).

Figure 4 Field-based characterization of the fracture populations. The stereoplots in the upper row show the 3D orientation distributions (black dots: poles to planes on a lower-hemisphere projection) with Kamb density contours in red (2σ interval) for all study sites. The lower row shows the respective kinematic indicators from fracture observations (lower-hemisphere projection). The dots represent the slip vectors, colored after fracture type and shear sense.

At site B, the measured fracture planes reveal a similar 3D orientation distribution with a slight clockwise rotation compared to site A (Fig. 4c). The NW–SE-striking set B_I, dipping SE, is oriented parallel to fracture set B_II and perpendicular to the NW–SE-striking set B_III. Both B_I and B_II exhibit dextral strike-slip movements (Fig. 4d), while fractures of set B_III mostly indicate normal movements. The rose diagrams in Fig. 3 show that the NE–SW-striking B_I and B_II fractures are the most prominent sets at all scales.

Also at site C, NE–SW striking and steeply SE fractures (C_I) and NW-dipping fractures (C_II) are the dominant sets (Fig. 4e). Set C_I exhibits large variations in strike directions, ranging from NNE–SSW to E–W. Both sets predominantly imply dextral strike-slip movements (Fig. 4f). The rose diagrams indicate that E–W to ENE–WSW are the dominant fracture strike orientations, corresponding to both fracture sets C_I and C_II (Fig. 3).

The fracture network at site D is characterized by two main sets, with a diffusely oriented, approximately NE–SW-striking set (D_I) and a dominant subvertical, NW–SE-striking fracture set (D_II) (Fig. 4g). For neither set D_I nor D_II do the kinematic indicators reveal a clear pattern (Fig. 4h). The rose diagrams in Fig. 3 reveal a bimodal distribution, with both sets D_I and D_II equally abundant.

In summary, we document a change in dominant fracture orientations, from NW–SE striking in the Aiguilles Rouges massif (site A), to E–W and NE–SW dominated in the limestone nappes in the center of the Rawil depression (sites B and C), and to a NE–SW and NW–SE bimodal fracture distribution in the Gastern massif (site D).

4.1.2 Size distribution

To combine the fracture trace data from the different scales, we first derive the fractal dimension D for each fracture trace map. Using the two-point correlation function C₂, we obtain values of D between 1.64 and 1.83 (Fig. 5), suggesting a relatively homogeneous spatial fracture distribution, as D= 2 would imply a perfectly homogeneous spatial distribution. This is in good agreement with the qualitative observation of widespread and pervasive fracturing in the area. Since the datasets from different scales reveal consistent values of D, we argue that the fractal dimension is scale-independent.

Figure 5 Estimation of the fractal dimension D for each individual fracture network based on the two-point correlation function C₂ (colored line). The black line shows the local slopes of C₂ (derivative of C₂; right-hand axis). The plateau in the local slopes within the grey shaded area approximates the fractal dimension D, which is represented by the dashed black lines.

Using the obtained values of D for the area normalization (see Sect. 3.1.2), we find that the normalized fracture size data follow a linear trend in a log–log plot for all study sites, supporting the assumption of an underlying power law distribution (Fig. 6). The power law scaling seems to hold over a length range of 3 orders of magnitude for fracture lengths between ca. 10⁰ to 10³ m. However, as can be seen for the data of site C, for example, individual datasets slightly deviate from the power law behavior. Such deviations can be caused by sampling bias such as the choice of the sampling area of the consecutive mapping scale, since subareas with anomalously high or low fracture densities can lead to an offset in the multi-scale size distributions. We will further discuss potential uncertainties in Sect. 5.1.

The obtained fracture power law exponents α_F, which are of main interest for comparison with the earthquake data, range between 2.43 and 2.73 (Fig. 6), which suggests that little regional variations exist. As the fracture networks are sampled in both crystalline basement rocks and the Helvetic limestone units, α_F seems to be relatively insensitive to differences in lithologies. For the dimensionless fracture density term c, we observe somewhat larger variations. With a value of $c=\mathrm{1.57}\times {\mathrm{10}}^{-\mathrm{4}}$, site B reveals the lowest fracture density term, while sites A, C, and D show larger c values of up to $\mathrm{5.76}\times {\mathrm{10}}^{-\mathrm{4}}$. Even though all c values range in the same order of magnitude, the fracture density term varies by a factor of 3 for the different sites, however, without any clear regional trends.

Figure 6 Normalized size distributions of the fracture networks for all four study sites, derived from the multi-scale approach. Circles denote the unbiased data, while crosses show the truncated and censored parts of the data. The different color shades represent data from the three mapping scales. The respective fracture density term c and the power law exponent α_F are given at the top right.

4.2 Earthquake ruptures

Next, we assess the orientation and size distributions of seismic rupturing at depth, derived from focal mechanism data and earthquake rupture lengths for the different subdomains denoted in Fig. 1b. We therefore use a compilation of 69 previously published focal mechanisms and 4612 earthquakes of the earthquake catalog of Diehl et al. (2021b) that lie within the seismotectonic domain z21 shown in Fig. 1 (see Sect. 3.2).

4.2.1 Orientation distribution

The data from the focal mechanism catalog reveal consistent rupture orientations within the different subdomains, with known earthquake ruptures predominantly occurring along subvertical, E–W- to NE–SW-striking planes (Fig. 7). While in subdomain A in the Aiguilles Rouges massif, the known active planes dip steeply towards the south (Fig. 7a), steeply north-dipping earthquake ruptures are observed additionally in the central and eastern subdomains C and D (Fig. 7c and d). For the data of subdomain B, currently no active rupture planes are provided in the published focal mechanism data, but ENE–WSW-striking nodal planes render similar active rupture orientations possible (Fig. 7b). Together with the observation of the NE–SW-directed clustering of the seismicity (Fig. 1b), this implies that the majority of the seismicity in the Rawil region occurs along subvertically oriented E–W- to NE–SW-striking fractures.

Figure 7 Three-dimensional orientation distributions of earthquake ruptures based on the nodal planes from published focal mechanism data (see text for references). Solid colored lines: known active fracture planes with slip vectors plotted as circles. Focal mechanisms without a known active fracture plane are represented by dashed grey lines (planes) and triangles (slip vectors). Secondary nodal planes for focal mechanisms with a known active plane are not shown in the figure.

4.2.2 Size distribution

As shown in Fig. 8, the size distribution of earthquake ruptures for the different subdomains all follow power laws at length scales between about 10² and 10³ m. For the data of subdomains A and D, we observe a deviation from power law behavior in the upper tail of the distribution, likely related to individual larger earthquakes. In subdomain A, this can be attributed to the M_L= 4.9 Vallorcine earthquake in 2005 (Fréchet et al., 2011), while in subdomain D, it relates to the M_L= 4.1 Salgesch earthquake in 2016 (Diehl et al., 2018). In subdomains B and C, we observe well-defined power law distributions for the earthquake rupture lengths. The three empirical M_w−A_R scaling laws used all yield similar values of α_R, which indicates that a significant dependency from the empirical scaling relations used on the obtained power law exponents appears unlikely (Fig. 8). In general, the derived earthquake rupture power law exponents α_R all range between values of 3.88 and 4.81 regardless of the earthquake rupture scaling law used (L14, T17, and WC94), implying little along-strike variations in the α_R values. As the size distribution of earthquakes often shows variations with depth, we further investigate a possible depth dependency in Sect. 5.4.

Figure 8 Three-dimensional size distribution of earthquake ruptures for the four subdomains. The earthquake magnitudes M_w were converted into rupture lengths L_R with the scaling laws of Leonard (2014) as L14, Thingbaijam et al. (2017) as T17, and Wells and Coppersmith (1994) as WC94 represented by the different grey scales, yielding three slightly different size distributions and α_R values. The values in the bottom left indicate the respective values of α_R.

Overall, the α_R and α_F exponents exhibit significant differences that are unlikely to be solely attributed to uncertainties. These differences will be further explored in the discussion section.

5 Discussion

5.1 Uncertainties in fracture size distributions

In the following, we discuss two of the major uncertainties related to data sampling that can affect the fracture size distributions: (i) the finite mapping resolution of fracture trace lengths and (ii) the choice of the subsampling mapping area as well as the potential bias due to manual mapping.

Mapping fracture traces on a fixed mapping scale on imagery with a given finite resolution necessarily leads to a bias in the derived fracture lengths, as the terminations of the fracture traces can only be mapped with a certain accuracy that is dependent on the image resolution. This potentially shifts the data along the horizontal axis of a log–log plot, leading to variations in the derived power law exponent α_F. We here use the rather conservative assumption that the fracture trace lengths are mapped with errors ϵ of 0.05 m (1:10), 0.5 m (1:100), and 5 m (1:1000) at the different mapping scales. We then randomly perturb all measured fracture traces with length values between [−ϵ, ϵ] with a uniform distribution and calculate the α_F value for 100 Monte Carlo (MC) simulations. As shown in Fig. 9a, the normalized density distributions show relatively little variation. The derived mean power law exponents α_F between 2.14 and 2.80 are similar to the values reported in Fig. 6, and the standard deviations of up to 0.15 suggests that the influence of uncertainties related to fracture length estimation on α_F is small (Fig. 9b).

Figure 9 (a) Assessment of the fracture trace length uncertainties. The plots show the normalized density size distributions with fracture trace lengths perturbed with values 0.05, 0.5, and 5 m. (b) Histogram of the derived α_F values. The values report the mean α_F value and the respective standard deviation. (c) Assessment of the subsampling area uncertainties, showing the normalized density size distributions with the number of fractures in each dataset randomly altered between −25 % and +25 %. (d) Histogram of the derived α_F values of the sampling area uncertainties. The values report the mean α_F value and the respective standard deviation.

As the multi-scale approach used in this study relies on fracture datasets from different scales, it is critical to assess the influence of the choice of different subsampling areas. The derived fractal dimension values below 2 suggest that the fracture density is not a perfectly homogeneous property in the investigated fracture networks (Figs. 3 and 5). Therefore, the spatial choice of the subsampling area for the next detailed scale leads to uncertainties in the size distribution, as differences in fracture densities shift the data along the vertical axis, affecting the calculation of α_F. Additionally, fracture traces may be overlooked or linear features not representing fractures may be interpreted as fractures during manual mapping of fracture traces, leading to a similar effect. When selecting subsampling areas, we tried to avoid specific geometric features (i.e., fracture linkages, fracture terminations) and focused on minimizing topographic effects and Quaternary cover, since these factors can significantly bias the power law exponent estimates, leading to an overestimation of α_F (Cao and Lei, 2018). To assess the influence of changes in fracture densities due to the arbitrary choice of the subsampling area, we randomly altered the number of mapped fractures from −25 % to +25 % and derive α_F for 100 MC simulations. Compared to the fracture trace length uncertainties, the effect of the subsampling area uncertainties is more pronounced (Fig. 9c and d). However, the derived α_F values still vary in a small range between 2.44 and 2.86, with standard deviations of up to 0.29.

We conclude from this analysis that both uncertainties related to the fracture trace length estimation and the choice of the subsampling area and manual mapping have a limited influence on the derived power law exponents α_F. We thus expect the true power law exponents of the fracture networks α_F to range between 2.5 and 3, with little regional and lithological variation.

5.2 Fracture network characteristics

The fracture networks from the four different study sites across the Rawil depression illustrate the pervasive but systematic fracturing in the study area (Fig. 3). Interestingly, the 2D orientation distributions for an individual site appear to be consistent across scales (Fig. 3), which is in agreement with the findings of Odling (1997). This implies that the fracture orientations are self-similar, and thus scale-independent, suggesting that a single-scale mapping should be sufficient to capture accurate information on the orientation distribution within a fracture network. Comparing the fracture networks of the different sites, however, we observe variations in the dominant strike orientations (Fig. 3). While in the ECMs, at both site A in the Aiguilles Rouges massif and site D in the Gastern massif, NW–SE- and NE–SW-striking fractures prevail, the limestone nappes overlying the central part of the Rawil depression are characterized by predominantly E–W- to NE–SW-striking fractures. NE–SW-striking fractures are, however, also present in both sites A and D, implying that these fractures are regional features, roughly following the direction of the RSF and the seismogenic RFZ (Fig. 1a).

The α_F values for the four different study sites are all in the same order and range between 2.5 and 3 (Fig. 6). As shown in Table 2, these results are in good agreement with previously published values that generally lie between 2 and 3, which suggests that this range can be seen as a first-order universal law of fracture size distributions, independent of the host rock lithology. Also the derived fractal dimension D values of around 1.7 are consistent with values reported in the literature (Bonnet et al., 2001, and references therein).

Table 2 Compilation of published power law exponents α_F from fracture size distributions sorted by the analysis approach used (single- vs. multi-scale power law fitting), with complementary information deduced where available. As both cumulative (CDF) and probability density functions (PDF) are used for fitting fracture length data, we have to make the CDF-based values comparable to our PDF-based α_F values by using α_CDF+1 (Bonnet et al., 2001). N: total number of data; α_F (published): published power law exponent and the distribution used; α_F (PDF): power law exponents converted to PDF values (Bonnet et al., 2001); lower and upper limit: approximate bounds of the power law scaling range. Bold font emphasizes PDF-converted values (directly comparable to our values).

The multi-scale approach used herein has the benefit to greatly expand the length range of the individual fracture trace datasets compared to single-scale analysis (Figs. 2b and 6). The power law distributions derived in this study cover fracture length scales of 3 orders of magnitude between about 10⁰ and 10³ m, and it is likely that the power law holds also for smaller and larger fractures. In nature, however, all power laws must have upper and lower physical limits (Bonnet et al., 2001; Torabi and Berg, 2011; Turcotte, 1997, 1989). The detected lower scaling limit of 10⁰ m likely does not reflect the actual physical length scale (e.g., caused by grain size effects), and we thus speculate that our analysis does not incorporate the lower limit of the power law scaling. Identifying the true limit of the power law would require a similar analysis of high-resolution datasets at smaller scales. The derived upper limit of around 10³ m is way below the assumed thickness of the seismogenic part of the crust of around 15 km, as well as the thickness of the Helvetic nappe system of a few kilometers (Burkhard, 1988; Levato et al., 1994; Pfiffner et al., 1997; Steck et al., 1997). However, stratigraphic layering effects within the Helvetic nappes (e.g., Pfiffner, 1993) could influence the upper limit of the power law size distributions (Odling, 1997), at least for sites B and C. Since the upper scaling limits do not significantly differ between the different lithologies, we consider this effect to be of minor importance within the observed length range. We thus assume that our data only cover a limited length range of the actual fracture size distribution and that the power law scaling holds for smaller and larger fractures than studied herein. This interpretation is supported by the fact that the previously reported α_F values for lengths between 10⁻¹ and 10⁵ m all lie in a similar range (Table 2).

5.3 Potential links to the seismicity

The strong vertical components due to the exhumation-related tectonics during the latest stage of the Alpine orogeny, accompanied by enhanced glacial surface erosion in the past 2 Myr (Häuselmann et al., 2007; Glotzbach et al., 2010; Valla et al., 2012; Fox et al., 2015; Sternai et al., 2019), led to the exposure of the fracture networks studied herein (Egli et al., 2017; Herwegh et al., 2020; Cardello et al., 2024). Based on observations of both exhumed fracture patterns, earthquake hypocenter locations, and focal mechanisms, previous authors suggested that the exhumed fracture networks in the region of the Rawil depression were formed under a similar stress field to that which prevails today (Cardello and Mancktelow, 2015; Maurer et al., 1997; Pavoni, 1980a, b). As earthquakes and fractures are inherently related, they are expected to exhibit similar orientation and size distributions (Bonnet et al., 2001; Scholz, 2007, 1998, 1997). In the following, we attempt to compare these statistical properties of the exhumed fracture networks with the same properties of the seismicity in the study region.

Comparing the orientation distributions of the fracture networks (Figs. 3 and 4) with the earthquake ruptures (Fig. 7) in a first step, we observe a correlation only in the center of the Rawil depression for study sites and subdomains B and C. The E–W- to NE–SW-striking fracture set C_I (Fig. 4e and f) correlates with the active fracture planes from the respective subdomain C (Fig. 7c). The dextral strike-slip kinematics of the NE–SW-striking fractures reveals striking similarities to the kinematics of the earthquake ruptures (Figs. 4f and 7c). In subdomain B, where no active fracture planes are deciphered (Fig. 7b), the average orientation of the fracture planes generally renders earthquake rupturing along E–W- to NE–SW-striking planes possible, which would correspond to the reactivation of fracture sets B_I and B_II (Figs. 4c and d). However, recent high-precision hypocenter datasets suggest that NW–SE-striking fracture planes are reactivated at least in the northern part of the Rawil depression (Diehl et al., 2024). In the ECMs on both sides of the Rawil depression, the correlation between the orientation distributions is less clear. In subdomain A, earthquake ruptures strike ENE–WSW (Fig. 7a), likely corresponding to the similarly oriented A_I fractures (Fig. 4a). It thus seems that the recent earthquakes mainly reactivate prevailing A_I fractures as strike-slip faults (Fig. 7a). The E–W-striking active rupture planes in subdomain D (Fig. 7d) are not reflected by the bimodal fracture network of site D at all (Fig. 4g), which is therefore likely not a good proxy for the E–W-directed seismicity in subdomain D. This might be due to the fact that the chosen study site lies slightly outside the main seismic corridor. This difference could be explained by mechanically detached Aar and Gastern massifs, which also exhibit significantly lower seismic activity (Fig. 1b). The mapped NW–SE-striking fractures could thus represent the eastern lateral boundary of the main E–W- to NE–SW-striking seismogenic RFZ (e.g., Lee et al., 2023). Overall, earthquake ruptures in the Rawil depression and the RFZ seem to occur mainly along subvertical E–W- to NE–SW-striking planes, and the exhumed fracture networks show similar directions at least in the central part of the Rawil depression above the main earthquake corridor (Fig. 1b). This implies that the active fractures in the basement are consistent with exhumed fracture networks in the Helvetic nappes, indicating a link across various structural levels and different lithologies. Combined with the observation that α_F values reveal little sensitivity to lithological variations, we argue that the fracture networks of sites B and C serve as a good first-order analog of the seismically reactivated fracture networks at depth.

5.4 Implications for orogen-internal settings

As earthquakes usually occur several kilometers below the surface where direct information about the prevailing fracture networks is sparse, the underlying processes of seismogenic deformation often remain uncovered. However, (i) the good agreement of fractures and earthquakes in terms of their orientation distributions, (ii) the insensitivity of α_F to lithological variations across different structural levels, and (iii) the observation that both features arguably formed in a similar stress field (Cardello and Mancktelow, 2015; Maurer et al., 1997; Pavoni, 1980a, b) imply that the exhumed fracture networks in our study area might serve as a good first-order analog for the fracture networks at depth. We therefore attempt to compare the size distributions of the fracture networks and earthquake ruptures in a next step. As the power law exponent of fracture networks is dependent on the dimensionality of the data (e.g., Bonnet et al., 2001), we have to consider that the power law exponents were obtained from 2D (fractures) and 3D (earthquake ruptures) observations. Stereological considerations have demonstrated that the differing dimensionalities can be corrected for by employing the equation ${\mathit{\alpha }}_{\mathrm{3}D}=\mathrm{1.28}\cdot {\mathit{\alpha }}_{-D}-\mathrm{0.23}$ (Bonnet et al., 2001; Borgos et al., 2000; Hatton et al., 1993). Consequently, the 3D power law exponents of the fracture networks α_F(3D) range somewhere between 3 and 3.6.

The power law exponent of the earthquake size distributions, commonly referred to as b value by seismologists, has been proposed to decrease with increasing differential stress and thus depth (Mori and Abercrombie, 1997; Spada et al., 2013; Scholz, 2019, 2015). We therefore assess the depth-dependent variations in the frequency and size distribution of earthquake ruptures. As no significant lateral variations in α_R values were observed within the investigated seismotectonic domain (Fig. 8), we calculate α_R values for depth slices that incorporate the entire seismotectonic domain. To minimize potential bias due to the data selection, we randomly iterate through different depth ranges and only fit a power law for datasets with more than 200 entries. Since the majority of vertical uncertainties in the modern era of the earthquake catalog are in the subkilometer range (Diehl et al., 2021b; Lee et al., 2023), we do not consider hypocenter uncertainties.

Figure 10 Depth dependency of the seismicity in terms of (a) the number of earthquakes, (b) the rupture length l_R of the individual earthquakes, and (c) the derived power law exponents α_(3D) for both the earthquake ruptures (greyscale) and the fracture networks (colored bars and red box). Error bars denote 1σ uncertainties of the power law exponents (vertical bars) and the depth interval used (horizontal bars).

As shown in Fig. 10a, earthquakes in the study area occur most frequently at depths between 3 and 8 km. The maximum rupture lengths lie, with few exceptions, on the order of 10² m for depths shallower than 3 km, with significantly larger ruptures up to 10³ m occurring at depths between 3 and 9 km (Fig. 10b). At greater depths, the seismicity becomes sparse and maximum rupture lengths generally decrease, which could be due to the increasing component of temperature-dependent viscous deformation representing the frictional-viscous transition at the lower depth end (e.g., Wehrens et al., 2016). As expected α_R values generally decrease with depth (Fig. 10c), in line with previous studies (e.g., Mori and Abercrombie, 1997; Spada et al., 2013). However, the values for shallow (< 3 km) and intermediate (3–9 km depth) crustal earthquakes differ significantly: while shallow earthquakes reveal α_R values of around 5 to 8, values for intermediate earthquakes constantly decrease with depth to values of around 3.5. Comparing the α_R values with α_F(3D) of the exhumed fracture networks, we notice that for shallow earthquake ruptures a significant discrepancy in the α_F(3D) values exists (Fig. 10c). For intermediate earthquakes, α_R and α_F(3D) values converge towards similar values at around 8 km depth. This observation could potentially be explained by the following two hypotheses.

First, it could be stated that similarly to α_R, also α_F(3D) values decrease with depth. This cannot be strictly proven, as fracture networks at depth are not directly accessible. The enhanced earthquake frequency at depths > 3 km, however, suggests that the majority of brittle fractures are formed at intermediate depth levels (Fig. 10a). Even though fracturing certainly occurs at shallow depths < 3 km during exhumation as well, the significantly lower earthquake frequency implies that the majority of fractures in the exhumed fracture networks were formed at intermediate depth levels. This interpretation is supported by the proposed formation depths of exhumed fractures in the surroundings of the Rawil depression derived from geological field observations (Cardello and Mancktelow, 2015; Gasser and Mancktelow, 2010; Ustaszewski et al., 2007; Cardello et al., 2024). Together with the absence of significant lithological variations in α_F and the consistency with previously reported values from a large variety of tectonic settings (Table 2), we consider this first hypothesis as unlikely.

As a second hypothesis, we thus propose that the α_F(3D) values are relatively constant with depth, and α_F(3D) and α_R converge towards similar values at intermediate depths. According to this argumentation, the decrease in α_R values therefore records variations in the reactivation potential of pre-existing fractures at different depth levels due to increasing differential stresses. As the differential stresses are smaller at shallow depths, rupture initiations are more likely to stop before growing into large earthquakes, which leads to larger α_R values (Mori and Abercrombie, 1997). In contrast, with increasing differential stress, the likelihood of a rupture initiation (i) to grow into a large earthquake that reactivates the full available fracture length and (ii) to jump across multiple fracture segments is increasing (Manighetti et al., 2007), leading to a convergence of α_R and α_F(3D) values at intermediate depths. We propose that the similarities in the power law exponents at depths of around 6 to 8 km reveal that earthquakes are more likely to rupture along the entire fracture length at these depths.

Apart from the influence of differential stress, previous studies have shown that also other factors such as the fracture type (Gulia and Wiemer, 2010; Petruccelli et al., 2019a, b; Schorlemmer et al., 2005) or material heterogeneities (Mori and Abercrombie, 1997; Mogi, 1962; Goebel et al., 2017) can lead to changes in the earthquake power law exponents. As the dominant orientations of both the fracture networks and earthquake ruptures are in good agreement (Figs. 4 and 7) and the focal mechanisms mostly imply transtensional movements (Kastrup et al., 2004; Houlié et al., 2018), we argue that the fracturing style has a minor influence on the observed variations in our case. Material heterogeneities, such as differences in lithologies, could potentially explain the increased α_R values at shallow depths above 3 km, as this depth correlates with the proposed basement-cover contact below the Rawil depression (Burkhard, 1988; Levato et al., 1994; Pfiffner et al., 1997; Steck et al., 1997; Lee et al., 2023). Stratigraphic layering effects due to the finite thickness of the Helvetic limestone units may limit the size of earthquake ruptures, increasing the α_R values at this lithological boundary (Ouillon et al., 1996). However, since the derived α_F(3D) values of the limestone and crystalline basement units are rather similar, and the discrepancy between the size distributions of exhumed fractures and earthquake ruptures occurs at shallow depths, we argue that the elevated α_R values at shallow depths cannot be explained by differences in material heterogeneities alone, in line with the findings of Scholz (1968). As an additional factor, the existence of fluids impacting pore pressure and friction can weaken the strength of the fractures in the uppermost crust, preventing accumulation of larger amounts of stresses.

In summary, the depth-dependent comparison of the size distributions of fractures and earthquakes suggests that the likelihood of an earthquake rupture to reactivate the entire fracture length is higher at intermediate crustal depths, while for shallower earthquakes differential stresses are rather small, rendering earthquakes that rupture the full available fracture length less likely, which is similar to the interpretations of Tormann et al. (2014) and Hetényi et al. (2018). As illustrated in Fig. 11, partial seismic rupturing thus seems to be the common seismogenic deformation mode at shallow depth in the Rawil depression region. However, it should be noted that other factors, such as lithological changes and the presence of fluids, can influence fracture strength. Consequently, the observed depth-dependent fracture reactivation potential may differ in other regions. To assess the universality of this observation, further investigations in different settings are required.

Figure 11 Schematic block diagram illustrating the depth-dependent seismogenic fracture reactivation patterns observed herein by relating the power law exponents of the fracture networks (α_F) to the exponents of recent earthquakes (α_R). Full seismogenic fracture reactivation appears to occur predominantly at depths between 3 and 9 km below sea level (b.s.l.).

6 Conclusion

The comparison of the statistical properties of fractures and earthquakes points towards differences in the seismogenic deformation processes in the investigated orogen-internal setting, characterized by pervasive fracturing and distributed seismicity. Our findings show that the potential of pre-existing fractures to reactivate seismically seems to vary with depth, likely associated with depth-dependent changes in differential stresses (Fig. 11): while partial seismogenic rupturing of fractures is the dominant deformation mode at shallow depths (< 3 km b.s.l.), earthquakes are more likely to reactivate the entire fracture length at intermediate depths (ca. 3 to 9 km b.s.l.) in the Rawil depression region. Initiating earthquake ruptures thus reactivates the full available fracture length predominantly at intermediate crustal depths.

The earthquake activity at shallow and intermediate depths furthermore documents ongoing brittle deformation in the Rawil depression region. As a consequence, recent fracturing processes likely enhance the permeability of pre-existing fractures, which could be of great interest for the exploration of the geothermal potential within such orogen-internal settings. In light of potential induced seismicity, our findings indicate that shallow fracture systems are less likely to rupture over the entire length than deeper ones. Furthermore, quantitative information about orientation and size distributions of fracture networks such as that obtained herein is crucial for accurate modeling of recent fluid flow within intricate fracture networks. Similar investigations incorporating the statistical properties of both fracture networks and earthquake datasets could therefore potentially be used for innovative geothermal exploration in orogen-internal settings.