Mature karst systems consist in well developed channels of various morphologies within a relatively low permeable rock matrix, typically consisting of soluble rocks such as limestones, dolomites or evaporites. Different hydraulic parameters between channels and a matrix lead to the so-called duality of karst (Király, 1994, 1998). On the surface, the epikarst, consisting of a highly altered surface bedrock, is thought to enhance the concentrated infiltration within the karst conduits (Király, 1998). Bonacci (1987) carried out numerous analyses on karst and found that the effective porosity of karst host rock is between 0.1% and 1% on average, being locally much higher in the presence of large caves or conduits. This configuration, i.e., high permeable conduits and a low porosity matrix, explains satisfactorily karstic spring hydrographs following recharge events. Figure 1 illustrates the characteristic hydrographs of springs related to 1) a fractured rock mass Figure 1A, 2) a bare karst Figure 1B, and 3) a fractured rock mass with karst conduits. The shape of the hydrographs depends on the extent of the fractures within the host rock. The hydrograph resulting from the combination of a fractured host rock with karst conduits (Figure 1C) represents the global response of both the high permeable conduits and low permeable matrix in response to recharge events. In this scenario, the recharge event involves the inversion of the flow regime. During flood events, the evacuation capacity of the karst system is exceeded and results in a rapid and large water level increase within the channels. This leads to a flow from the channels to the matrix, recharging the matrix. During low-water periods, the water stored in the matrix flows back to the conduits, maintaining base flow levels. Such hydrographs are typical of karstic regions and are clear evidence of rapid hydraulic head increases. For example, Figure 2 shows the discharge of the Schlichenden spring draining the Hoelloch cave system (Muotathal, Switzerland). During an important recharge event, a cave located 540 m above the permanent water table was flooded Miller (2008).

Another important aspect to consider is the relationship that can exist between karst systems and fault zones. A fault zone’s hydraulic properties depend on various parameters and are subject to changes over time and space. Considered either as a barrier or conduits for fluid circulation, their role depends greatly on their anatomy (fault core and damaged zone) and composition (breccia, cataclasite, or ultracatalasite) (Caine et al., 1996; Faulkner et al., 2010), which will influence the spatial distribution of permeabilities (Caine et al., 1996; Evans et al., 1997). Additionally, the hydraulic properties of a fault zone can vary in time as has been observed in laboratory tests (Seront et al., 1998) and modelling experiments (Preisig et al., 2016), in which a shear slip on a fault increases the permeability in response to dilation. As karst development is driven by water flow and faults can favor flow paths, it is likely that karst systems develop in relation to faults and that the two systems are hydraulically connected. However, the exact role of the fault, i.e., whether fluids can pass through the fault zone or flow along, is determined by its internal structure and composition (Matonti et al., 2012). This potential link between karst and faults makes the latter particularly sensitive to rainfall and karst flooding events.

External changes in the stress regime might impact the stability of the fault depending on its orientation, initial stress regime, and pore pressure. We identify four potential mechanisms to explain the influence of rainfall on fault stability:

1. crustal loading,

The impact of crustal loading 1), i.e., the weight of precipitation as an additional load on the crust, depends on the stress regime. From Anderson’s stress regime classification (Anderson, 1951), an additional load in a normal faulting stress regime (NF) reduces the stability of the fault, increases stability in a thrust-faulting stress regime (TF), and leaves the stability unchanged in a strike-slip stress regime (SS). Additionally, the additional load on the crust modifies normal effective stresses due to increasing pore-fluid pressure through poroelastic deformation 2), reducing the shear strength and hence the stability of the fault. A rising fluid pressure at the surface might also initiate a pressure diffusion front 3), increasing pore pressure and reducing the effective normal stress gradually while the pressure front propagates. A rising water level and a direct hydrogeological connection through the karst between the surface and the fault plane 4) result in an instantaneous increase in pore pressure, impacting the stability of the fault.

We propose a conceptual model (Figure 3) that considers specific karst features and rainfall-induced mechanisms that potentially lead to fault instability. An important difference between these mechanisms depends on the relationship between the karstic network and the faulted system and whether a direct hydraulic connection exists (iv) or not (ii and iii). In cases in which the karstic network is connected to a highly permeable fault zone (iv), the groundwater pressure increase at a given depth is quasi-instantaneous. In cases in which the karstic network lies above a fault, and there is no direct hydraulic connection, the mechanism leading to stress changes is similar to the one described for reservoir impoundment, i.e., (ii) a rapid increase in pore pressure resulting from the elastic deformation due to the vertical load and undrained conditions and/or (iii) a delayed response due to fluid migration towards depth. Although an increasing groundwater pressure occurs in all scenarios, they can be differentiated by the timing, which is quasi-instantaneous for (ii) and (iv) and delayed for (iii).

According to the classical Mohr–Coulomb failure criterion, failure on a fault initiates when the following stress conditions are met: τ = c + μ σ n ′ (1) where τ is the shear stress, μ the coefficient of internal friction, the effective normal stress σ n ′ = σ n − p , with σ _n the normal stress, c the cohesion, and p the pore pressure in the fault zone. The coefficient of internal friction μ typically lies between 0.6 and 0.85 (Byerlee, 1978) and in the context of reactivation, faults are often considered to be cohesionless (c = 0) (Sibson, 1985). For constant fault properties, the external parameters affecting fault stability are the shear stress τ and the effective normal stress σ n ′ , which are both dependent on the fault’s orientation (angle θ between σ 1 ′ and the normal to the fault plane) and stress magnitude as expressed in two-dimensional plane strain conditions or when the intermediate stress is coplanar with the fault by the following relations (Jaeger et al., 2009): σ n ′ = σ 1 ′ + σ 3 ′ 2 + σ 1 ′ − σ 3 ′ 2 cos 2 θ (2) τ = σ 1 ′ − σ 3 ′ 2 sin 2 θ (3) For a fault favorably oriented (θ = 60°) and according to Eqs 1–3, rupture can be initiated by an increasing shear stress τ through an increase of σ 1 ′ and/or a decrease in σ 3 ′ or through a decrease of the effective normal stress σ n ′ .

To evaluate the proximity of a fault to failure, the Coulomb failure stress is a good indicator and is defined as follows (see, e.g., Beeler et al., 2000): C F S = τ − μ σ n − p − c (4) A negative value resulting from Eq. 4 indicates that failure has not been reached, as the shear strength ( c + μ σ n ′ ) is greater than the shear stress (τ). The change in Coulomb failure stress shows the impact of stress changes on the stability of the fault. Coulomb failure stress change is defined as follows (Beeler et al., 2000): Δ C F S = Δ τ − μ Δ σ n − Δ p (5) where Δσ _n is the change in normal stress and is positive for compression. Positive values of ΔCFS indicate that the fault is closer to failure, whereas negative values indicate that the fault is moving away from the rupture point.

If the fault zone is in direct hydraulic connection with the karstic network, and considering a simplified environment with no pressure losses, an increasing hydraulic head in the karst conduits due to heavy precipitation will generate a quasi-instantaneous and non-attenuated pressure signal. The increase in pore pressure (and the decrease of the effective principal stresses σ 1 ′ and σ 3 ′ ) will directly relate to the increase in the hydraulic head ΔH by the relation ΔP = ρgΔH, leading to a decreasing σ n ′ .

The variation in the hydraulic head ΔH [m] depends on the storage coefficient S [-], a considered area A [m²], and the volume of infiltrated water V _w [m³], which is controlled by the (effective) rainfall intensity i [m/d] falling over the area A during a given time dt [d] (Fetter, 2001): Δ H = V w S A (6) where V w = i d t A (7) Combining Eqs. 6, 7 gives: Δ H = i d t S (8) In a porous rock mass in which the water table can fluctuate freely, S is equal to the effective porosity, which can range from less than 0.01 to 0.30 in sedimentary karstic rocks. In a porous rock mass in which the water table cannot fluctuate freely, i.e., the piezometric level is above the top of the permeable unit (confined aquifer), S is in the order of 0.005 or less (Fetter, 2001).

Considering a karst network without a direct hydraulic connection with the fault zone, fault stability might be influenced by a supplementary load, i.e., the infiltrated water. In such a case, an underlying saturated fault system will behave as a poroelastic medium with increasing pore pressure (and decreasing effective stress) in response to the vertical loading.

An increasing hydraulic head will lead to fluid pressure diffusion from zones of high to zones of low hydraulic head. Considering a step increase in pressure at time t = 0, the simplified one-dimensional diffusion problem has the solution (Carslaw and Jaeger, 1959): p z , t = P a erfc z 2 D t (12) where erfc is the complementary error function, p is the pressure at depth z and time t resulting from the applied pressure P _a on a medium with diffusivity D, where D = T/S, with T the transmissivity and S the storage coefficient.

Now considering a step increase in pressure at time t = 0 and a linear decrease over a time τ from P _a to 0, the solution reads as follows (Perrochet P., personnal communication, presented in Supplementary Material S1): p z , t = P a erfc z 2 D t − m t + z 2 2 D erfc z 2 D t − z t π D exp − z 2 4 D t + H t − τ m t − τ + z 2 2 D erfc z 2 D t − τ − z t − τ π D exp − z 2 4 D t − τ (13) where H(u) is the Heaviside step function (H(u) = 1 when u > 0 and H(u) = 0 when u < 0) and m = P _a /τ. Note that Eqs. 12 and 13 do not consider the coupled poroelastic response.

In 2001, the University of Munich, together with the Bavarian State, installed a new seismic monitoring network in Bavaria, as well as a subnet of six stations around Mount Hochstaufen (Kraft et al., 2006). These were installed following preliminary observations made by Schwarzmann (1996), who analyzed 2 months of seismic records from a temporary network installed after a M _L 2.7 of 26 June 1995 in the region of Bad Reichenhall. The author suggested a possible relationship between seismicity and rainfall.

Mount Hochstaufen experienced two seismic swarms in 2002, both occurring after an extended period of rainfall in March and August. In August, ∼300 mm fell over the region of Bad Reichenhall (German Weather Service, DWD) in less than 1 week, with a maximum intensity of ∼120 mm/day and ∼110 mm/day on the sixth and the 11th of August, respectively. In all, 1,171 earthquakes were recorded with a maximum event magnitude of M _L = 2.4. The located events have epicenters distributed over the Staufen Massif with a focal depth mostly lying between 0 and 3 km depth. Figure 4 summarises the geologic context around the Staufen Massif, which is part of the Northern Calcareous Alps (NCAs), forming an east-west striking fold-and-thrust belt. The NCAs are part of the Upper Austroalpine unit and are commonly divided into three main nappes units: (1) Bajuvaric, (2) Tirolic, and (3) Juvavic. Mount Hochstaufen is part of the Tirolic nappes overthrusting the Bajuvaric Nappes, which are overthrusting the Helvetic and Penninic units (Figures 4A, B) (Tollmann, 1976; Henrich and Zankl, 1981). Mount Hochstaufen is located south of this major thrust fault oriented E-W and north of an important SW-NE sinistral strike-slip fault called the ISAM fault (Innsbruck-Salzburg-Amstetten). Contemporaneous stress filed measurements provided by the Stress Map Germany (Reiter et al., 2016) suggest an N-S maximum horizontal stress orientation.

The stratigraphic sequence of the Staufen Massif consists in units from the lower to middle Triassic with lithologies mainly dominated by limestones and dolomites. Figures 4C, D show the lithologies, along with their hydraulic properties (Bayerisches Landesamt für Umwelt, last accessed in January 2021) The higher part of the Staufen Massif consists of Wetterstein limestone (Tollmann, 1976). This formation, along with the Dachstein formation outcropping in the southern and eastern NCAs, is the most important karst host rock in the NCAs.

For the analytical modelling, we consider an average altitude of 1,500 m.a.s.l and an average focal depth of 1,000 m.b.s.l. The altitude of a perennial karstic spring located 4.5 km west of Hochstaufen is taken as the base water level. The horizontal distribution of epicenters simplified as a circle with a radius of 5 km is also considered.

In August 2005, the regions of Muotathal and Riemenstalden in central Switzerland experienced a succession of events over a 12-h period after 4 days of heavy rain with a total precipitation of ∼200 mm, with a maximum intensity of ∼75 and 90 mm/day on the 21 and 22 of August. The return period of such intensive rainfall periods is approximately 300 years (Husen et al., 2007). A total of 47 earthquakes were recorded and distributed in two clusters (regions of Muotathal and Riemenstalden). The focal depths are between 0 and 5 km. Muotathal and Riemenstalden are both located in the Helvetic unit in central Switzerland and are part of the Axen nappe (Pfiffner, 2011). The Axen nappe consists of a formation from the Triassic to the Paleogene. In the area of interest, the lithologies are dominated by Cretaceous limestones. Due to extensive thrusting, the units are stacked on top of each other, leading to an important thickness of Cretaceous units Figures 5A and B (Schmid et al., 2004).

A major thrust delimits the Axen nappe with the Drusberg nappe close to the villages of Muotathal and Riemenstalden. Numerous fault lineaments are visible on the geological map and show an N-S general trend as well as minor lineaments showing a WNW-ESE strike. Some of these faults show distinct horizontal offsets, indicating a strike-slip environment. Data provided by the World Stress Map (Reiter et al., 2016) suggest a ∼SE-NW orientation of the maximum horizontal stress for this region.

With a total length of more than 250 km of galleries (Wildberger et al., 2010), the region is a well-known karstic area. The Hölloch cave and Silberen system are part of this karstic network and are located east of the village of Muotathal. Karst developed mostly in the upper Cretaceous, in the Urgonian and Seewen limestones, and as a result of thrusting, the thickness of these karstified rocks can be doubled or tripled (Figure 5C) (Wildberger et al., 2010). The lithologies and their hydraulic properties are presented in Figure 5D.

For the analytical modelling, we consider an average altitude of 1,400 m.a.s.l. and an average focal depth of 2,000 m.b.s.l. The altitude of the perennial Schlichenden karstic spring, the outlet of the Hölloch cave system, is taken as the base water level. The horizontal distribution of epicenters simplified as a circle with a radius of 2.5 km is also considered.

In south-eastern France, in the Gard department, a catastrophic storm occurred in September 2002. In 2 days, accumulated precipitation reached ∼400 mm in Pujaut Basin, near Avignon, and even reached 650–700 mm further to the north, near the city of Alès. Seismicity increased during the 2 weeks following the storm. Eleven earthquakes were recorded within a temporary network centered around the Pujault Basin, and two events were located outside the network but still relatively close in distance. The focal depth mostly lies between 0 and 10 km (Rigo et al., 2008).

The Gard department located in the Provence region between the Alps in the northeast and the Pyrenees in the southwest. The regional geology consists of a sedimentary cover of Mesozoic and Cenozoic ages, varying in thickness from 2 to 7 km (Figures 6A and B) (Rigo et al., 2008). Three major faults exist in the area of interest: the Nîmes Fault, the Pujaut fault, and the Roquemaure fault. The Nîmes fault has an NE-SW orientation and a senestral strike-slip movement (Schlupp et al., 2001). The Pujaut fault, parallel to the Nîmes fault, corresponds to a normal fault steeply dipping towards the north and is the result of the Oligocene extension regime (Jolivet and Faccenna, 2000). The northern compartment of the Pujaut fault results in a depression called the Pujaut basin or graben. The Roquemaure fault is considered inactive or even non-existent (Rigo et al., 2008). Further east, the Salon-Cavaillon fault system strikes SE-NW and has a normal component possibly combined with a dextral sense of movement (Terrier et al., 2007). The present-day maximum horizontal stress orientation is thought to be N-S to NW-SE (Jolivet and Faccenna, 2000).

The area is well known for its deep karstic conduits that developed during the Messinian crisis (Julian and Nicod, 1984; Audra et al., 2004). Exploration by a remotely operated vehicle (ROV) in the Fontaine-de-Vaucluse brought evidence of former epiphreatic karst conduits down to 100 m below sea level (Bayle and Graillot, 1987). From seismic data, drilling, and field observations, evidence shows that canyons are lying approximately 1,000 m below sea level near the city of Avignon (Clauzon, 1982; Gargani, 2004). At the end of the Messinian crisis, the canyons were flooded and filled with terrigeneous sediments Clauzon (1982). Figure 6C shows the relationship between the deep karstic system and the ancient canyons and Figure 6D illustrates the lithologies in the area of interest.

For the analytical modelling, we consider an average altitude of 200 m.a.s.l. and an average focal depth of 3,000 m.b.s.l. The altitude of the Rhone River near Avignon is taken as the base water level. The horizontal distribution of epicenters simplified as a circle with a radius of 19 km is also considered.

We quantify the effect of each mechanism on fault stability for the different case studies by considering a semi-infinite environment with initial stress conditions of σ ₁, σ ₃ and the initial pore pressure p _init at depth z, with cohesion c and internal friction coefficient μ. The environment is saturated at depth and unsaturated towards the surface with an unsaturated layer of thickness z _unsat (see Figure 7). Further considered is a circular area on the surface of the saturated layer of radius r being loaded with pressure P, which is inferred from the rain intensity i.

Values for z, z _unsat , r, and P are based on the case studies and initial stress conditions are estimated using average lithostatic and hydrostatic pressure gradients of 0.0225 and 0.0098 MPa/m, respectively. The depth z is inferred from the focal depth, the unsaturated thickness z _unsat is calculated by subtracting the altitude of local springs from the average altitude of the area. The radius r corresponds approximately to the horizontal spatial distribution of seismicity. The pressure P is inferred from the rain intensity. In the case of a direct hydrogeologic connection, the pressure is quasi-equivalent to ΔH (Eq. 8). In the case that there is no hydrogeologic connection, the pressure load corresponds to the additional weight of the infiltrated water. For the stress conditions, we consider a strike-slip environment where σ _H = σ ₁, which is at least equal to σ _v , and σ _h = σ ₃, which is set to the minimum possible value to have a CFS < 0, meaning that rupture has not yet occurred. Cohesion c is thought to be negligible and a value of 0.6 for the internal friction coefficient μ is taken.

Input values for the analytical modelling are shown in Table 1. Precipitations are inferred from the case studies as shown in Figure 9. Note that we slightly simplify the rainfalls for the analytical modelling by considering only two rainfall events for each case study (i ₁ and i ₂). The storage coefficient S is set according to the observed water level fluctuations in the karstic network reported in the literature, which range between a few tens of meters (Lesparre et al., 2017) to several hundreds of meters (Jeannin, 2001). For Muotathal, a water level increase of 540 m for 200 mm of rain corresponds to a storage coefficient of 0.0004 (Eq. 8). As Muotathal is the only case study in which this value can be calculated, we decided to use a more conservative value of S = 0.001, which corresponds to a level fluctuation of 100 m for 100 mm of infiltrated water.

Using Eqs. 1–12 and data from Table 1, changes in pore pressure are calculated, as well as changes in the Coulomb Failure Stress (CFS). Results are separated between direct and delayed effects. Direct effects include changes in pore pressure resulting from the hydraulic head increase, as well as poroelastic deformation resulting from an increasing vertical stress, while the delayed effect is a pore pressure increase resulting from the diffusion process.

For poroelastic computations, we assume undrained conditions. We consider high rainfall intensity events, which implies that the drainage capacity of the karstic network is overcome during the rainfall. This results in a rapid increase in the hydraulic head within the karstic channels. For the delayed effect, simplified drained conditions are considered and represented by a linear decrease occurring over 12 days from the maximum pore pressure increase induced by ΔH at time t = 0 to 0 at time t = 12 days. The diffusion front is considered to migrate from the bottom of the karstic network towards focal depths. The distance d varies according to the extent of the karstic network below the water table. Two different values are tested. The first one, d ₁, considers that the pressure front starts migrating from the top of the saturated thickness, meaning that the distance d ₁ corresponds to the total depth z minus the unsaturated thickness z _unsat . The second, d ₂, considers a saturated karstic network with a thickness z _Ksat of 1,000 m, meaning that the distance d ₂ corresponds to the total depth z minus the unsaturated thickness z _unsat minus z _Ksat (Figure 7).

Note that the analytical modelling aims to compare the mechanisms with each other, not to conduct a sensitivity analysis of the different parameters.