Remote Triggering of Damage Followed by Healing Recorded in Groundwater Pressure

Abstract

Water levels in three adjacent water wells in the Yarmouk Gorge area have all responded to the 2020 Elazığ Mw 6.8 teleseismic earthquake. Water levels in two aquifers exhibited reciprocal behavior: during the first eight days after the earthquake, water level decreased by 40 cm in the deeper highly confined aquifer, and increased by 90 cm in the shallower less confined aquifer. The recovery of the water levels in both aquifers continued for at least three months. We interpret these observations as reflecting the increase in damage along the fault at the Yarmouk Gorge. Ground shaking increased the damage and permeability of this fault, temporarily connecting the two aquifers, allowing flow from the deep aquifer to the shallow one. Model results showing decreased permeability suggest that the fault healed by one order of magnitude within three days. This is the first documentation of decrease in permeability in a fault zone within such short time scales.

1. Introduction

Hydraulic and mechanical properties of rocks change during and after remote earthquakes due to damage and healing/sealing processes [1]. When the fault is active during seismic-induced shaking, the permeability of the fault increases due to fracture damage enhancement while the mechanical properties decrease [2]. Upward migration of fluids from an aquifer with higher pressure to another one through activated faults is expected at this stage [3]. Over time, faults recover their mechanical properties through fracture healing/sealing and the permeability decreases as a result of a combination of chemical and mechanical processes [4,5,6].

Aquifer permeability has been observed to be enhanced by transient stresses induced by distant earthquakes [7,8,9]. Other studies show the co-seismic effects of earthquakes on fault permeability [10,11] or long time permeability variations [12]. Wang et al. [13] showed that a large earthquake creates vertical permeability in groundwater systems by breaching the interlayered aquitards. They hypothesized that breaching of aquitards by large earthquakes is achieved by the formation of many hydraulically conductive cracks across the aquitards.

Groundwater near active faults is affected by stress accumulation associated with tectonic loading, which implies strong feedback between fault strength and permeability [14]. The processes responsible for changes in the strength and permeability of faults include healing and sealing [15] whereby healing involves closure of microscale cracks driven by surface energy [16], and sealing of cracks is carried out by precipitation of minerals transported from distant or local sources [17]. Subsurface fluid flow is expected along faults, whenever earthquakes create pathways (by increasing the permeability). In reactive systems, narrow channels are quickly blocked by mineral precipitation and the permeability is reduced [18]. Mineralized fluids with elevated temperatures accelerate both healing and sealing processes reducing the permeability and terminating the flow.

Dieterich [19,20] performed frictional experiments with different values of normal stress and hold times up to 105 s, and obtained a time-dependent logarithmic increase of the static coefficient of friction, during several hours. The results were interpreted as reflecting enlargement of the real contact area with time, as a result of indentation creep around geometrical asperities [21,22]. Recent laboratory studies have shown that micro-cracks can last for hours to months under a range of hydrothermal conditions [16,23,24,25,26]. Geophysical evidence for various durations of post-seismic healing processes, has been observed in Landers [27], Parkfield [4], Tokai region, Japan [28], and Loma Prieta [29] earthquakes. The recovery of P and S wave velocities (healing) following the 1992 Landers earthquake was recorded to last for years [27]. Cochran et al. [30] documented reduced elastic moduli along the Calico fault using seismic travel times, trapped waves, and InSAR observations that represent the cumulative mechanical damage from past earthquake ruptures. They inferred that the long-living fault damage persisted for hundreds or thousands of years with no significant healing.

The various healing mechanisms cannot be fully isolated from each other in recovery experiments on natural rock samples, and the recovery of damaged rocks in nature is a complex interplay between various healing mechanisms [15]. Therefore, we will not attempt to distinguish between the mechanisms and refer to the healing/sealing as a whole. While fracture, damage and permeability can heal relatively quickly in the experiments, in real fault zones, fast healing/sealing has yet to be documented [31]. This might be due to the fact that the tools used to monitor healing/sealing process in the field do not have the resolution for time scales smaller than weeks. In this paper, we show evidence for the remote earthquake-induced (24 January 2020, Elazığ Mw 6.8) sudden permeability increase of a fault zone, connecting two sub-aquifers followed by the gradual healing/sealing within a few days.

2. Hydrological Setup

The Yarmouk River and the groundwater resources in its watershed are a strategic transboundary freshwater resource of Syria, Jordan, and Israel. The scattered springs of Hammat Gader (Figure 1) discharge hot pressurized groundwater (45.5 MCM/year up to 50 °C) through extensive deep faulting. Groundwater emerges with different hydraulic heads, chemical compositions, and temperatures [32]. Recent studies attempted to explain the mechanisms responsible for the outflow of hot springs and high-pressure groundwater. Siebert et al. [32] distinguished the source of the springs and wells in the area. Groundwater originates at three different sources: groundwater flows from Golan Heights (north) and Ajloun-Plateau (south) towards the Yarmouk Gorge, where both flows are hydraulically separated by a major fault. Some of the springs and wells in the lowermost Yarmouk Gorge recharge at the Syrian Hauran Plateau (east).

Groundwater pressure data from three wells (

Table 1. Basic data of the three wells.

Well Name	Depth (m)	Casing Length (m)	Diameter of Well (m)	Diameter of Bare Hole (m)	Aquifer Lithology
Meizar 1	1250	954	0.34–0.71	0.31	dolomite
Meizar 2	807	448	0.34–0.51	0.31	dolomite
Meizar 3	336	80	0.38–0.51	0.35	chalk, marl

) is analyzed in this work (see Figure 1 for location). The Meizar 1 and 2 deep wells represent water from the upper Cretaceous Judea group and the shallow Meizar 3 well represent water from the Senonian Ghareb-Mishash formations. Meizar 2 and 3 are located 20 m apart and Meizar 1 is located three km to the North of Meizar 2 and 3 (Figure 2). Meizar 2 and 3 are artesian wells with an over-pressure of ~1.7 and ~0.7 bar, respectively. Meizar 1 is not an artesian well due to its high elevation (160 m higher than Meizar 2 and 3).

3. Pumping Tests

Meizar 2 and 3 wells are artesian and pumping tests were performed by opening the valves and measuring the discharge and pressure. The hydraulic connection between the Judea and Mishash-Ghareb aquifers was very limited as observed during pumping tests (Figure 3). However, the aquifers were elastically connected and pumping from Meizar 2 well caused a temporary small increase in the water pressure in Meizar 3. This phenomenon is called the reverse water-level effect [33,34] and it emphasizes the lack of hydrological connection between two aquifers. Water in the wellbore was hot after the pumping (42 °C) and cools down with time to 20 °C, decreasing the pressure after the valves are closed.

It has been observed in many field tests that when groundwater is pumped from an aquifer, the water level may momentarily increase in adjacent aquifers and aquitards [34,35,36,37,38]. This reverse water-level response was first observed in Noordbergum, The Netherlands, where the phenomenon earned its name for a special case where the water-level rise occurs in an aquifer separated by a confining unit from the pumped aquifer [35].

The loss of water due to pumping involves a decrease of volume, producing radial solid displacements towards the well. The displacement field (loading) compensates for the loss of water. For deeper layers, this displacement field causes compression, and hence an increase of pore pressure [35]. This loading information is transferred immediately throughout the elastic medium because of stress compatibility, while water flow and re-distribution of the pore pressure due to pumping, requires time.

4. Hydrological and Seismological Data

Water pressure is monitored using PA 33X in the artesian wells and PA 36XW in the non-artesian wells pressure transmitters (Keller-Druck, Winterthur, Switzerland), configured to a range of 0–0.1 MPa and 0–0.3 MPa, respectively. The accuracy and precision of these transmitters are 50–150 Pa and 2–6 Pa, respectively, depending on the pressure transmitter measuring range. Temperature dependencies and nonlinearities of the pressure transmitters are compensated by the transmitters’ microprocessor. Water pressure is recorded at 40 samples per second. We do not record temperature or electrical conductivity measurements because lowering sensors to the uncased section of the wells is very difficult due to the artesian conditions. Seismic data are obtained from station GLHS (Figure 1), of the recently upgraded Israel Seismic Network [39]. GLHS is a station with collocated sensors, recording ground accelerations, velocities, and displacements. The velocity data used in this study was measured by a Trillium 120 PH broadband seismometer (Nanometrics), located within a posthole, at a depth of 5.5 m below the surface, and recorded by a 6-channel Centaur (Nanometrics), at 200 samples per second, and 40 V peak to peak sensitivity.

The seismic waves generated by the remote 2020 Elazığ Mw 6.8 earthquake [40,41] are clearly seen in the seismic and hydrological data (Figure 4). During the passage of the waves, there was an average water pressure increase in Meizar 1 well. This sustained increase was also observed in previous cases [42,43]. The highest amplitudes occur after the arrival of the S waves; in Meizar 2 and 3 wells, the amplitudes are 3 m and 1.2 m, respectively; the amplitudes in Meizar 1 are relatively small (5.5 cm), due to its non-artesian nature. The attenuation in the high-frequency range in non-artesian wells is an effect of wellbore storage, where water is actually displaced during oscillations; in closed artesian wells, the attenuation is much smaller [44]. Similar trends are observed in the seismic data and in Meizar 2 and 3 wells, emphasized for specific bandwidths. Such data obtained by applying a Butterworth low-pass filter (with four poles at a cutoff frequency of 0.02 Hz) are provided in Figure 5.

Immediately after the passage of the seismic waves, pressure started to change in all wells. The pressure in Judea aquifer started to decrease and the pressure in the Ghareb-Mishash aquifer started to increase (Figure 6). Unfortunately, industrial pumping in Meizar 2 began soon after the event. Therefore, we discuss here only the water pressure changes in Meizar 1 and 3 wells. During the first 10 days, the water level in Meizar 1 well decreased by 40 cm and then recovered towards its pre-event value. In Meizar 3 well, the water level increased by 90 cm and then started to decay back to its pre-event value (Figure 6). We interpret these changes as fault activation and damage increase in a relatively narrow zone connecting two aquifers, followed by gradual healing/sealing of this zone. Several floods in the nearby Yarmouk River, that recharged the shallower Mishash-Ghareb aquifer without affecting the deeper Judea aquifer, illustrate the separation between the aquifers.

5. Groundwater Modeling

5.1. Model Formulation

Fault activation and damage heeling is incorporated into the groundwater modeling. The results of Dieterich’s [19,20] healing experiments demonstrate the logarithmic increase of the static coefficient of friction f_s:

$$f_s=f_0+Alo{g}_{10}\left(1+Bt\right)$$

(1)

where t is the duration of the hold time in seconds, f₀ ≈ 0.6–0.8, A ≈ (1–3)10⁻² and B ≈ 1–2 s⁻¹. The non-dimensional A-coefficient defines the magnitude of the static friction change with time and the attribute of B provides a time scale for the evolution of the static friction coefficient. The relation of the rate and state dependent friction model [22] connects the steady-state value of the friction coefficient, f_ss, with the sliding velocity at slip distances larger than a certain critical value:

$$f_{ss}=f_0+\left(a-b\right)ln\left(\frac{V}{V_0}\right).$$

(2)

The a and b non-dimensional values control the transient stage, while a defines the direct effect to the instantaneous slip velocity change and b defines the recovery to a new steady-state friction value corresponding to the slip velocity V normalized to the reference velocity V₀. The difference (a − b) defines the overall change in the steady-state friction value. Numerous experimental studies demonstrated that both the a and b parameters are of the order of 10⁻¹ to 10⁻² and the sign of the (a − b) defines the stability of the sliding contact [22]. Elevated values of a, b ~0.15 have been reported for wet frictional experiments in hydrothermal conditions [45,46]. By comparing laboratory results and microphysical models Van den Ende et al. [47] suggested that the b parameter increases with depth (temperature and pressure) approaching b~0.2. Larger values of b up to 0.35 were reported for gouge-bearing faults [48], and in general agreement with recent micromechanical models [23]. These elevated values were also used in dynamic rupture modeling [49,50].

Motivated by the observed logarithmic increase of the static coefficient of friction, Lyakhovsky et al. [51] used a damage-dependent function for the kinetics of healing and suggested the rate of healing under constant compaction strain, ε_cmp, to be of the form:

$$\frac{d\alpha }{dt}=-C_1{\epsilon }_{cmp}^{2}exp\left(\frac{\alpha }{C_2}\right).$$

(3)

where α is the damage parameter, and $C_1$, $C_2$, are two positive coefficients. Under static load conditions, the damage healing Equation (3) predicts logarithmic-in-time healing:

$$\alpha \left(t\right)={\alpha }_0-C_{2}ln\left[1+\frac{C_1}{C_2}{\epsilon }_{cmp}^{2}exp\left(\frac{{\alpha }_0}{C_2}\right)t\right],$$

(4)

with the same structure as of Equation (1). Lyakhovsky et al. [52] also demonstrated that the model with exponential healing (Equation (3)) reproduces the main features of the rate and state friction and connected the parameters $C_1$ and $C_2$ to the static and steady-state friction coefficients. According to their derivations, $C_2 ~b/f_{ss}$ is between 10⁻¹ and 10⁻². The value of the coefficient $C_1~B{C}_{2}exp\left(\frac{{\alpha }_0}{C_2}\right)/{\epsilon }_{cmp}^2$ strongly depends on the value of the coefficient $C_2$ and varies within several orders of magnitude. For two different values of $C_2$ = 0.03, 0.05, and B~1 s⁻¹, ε_cmp ~10⁻², the coefficient $C_1$is ~10⁻¹² s⁻¹, ~10⁻⁶ s⁻¹, respectively. Fault healing at seismogenic depth of 5–10 km and pressure of about 200 MPa is very fast with these parameters during the initial few hours after an earthquake and the rate of healing significantly decreases with time [53]. With healing, the value of damage (α) decreases, leading to permeability reduction. Experimental studies [54,55,56] confirm the correlation between the log-scale permeability variation and micro-crack density. In this model, the damage is a measure for the micro-crack density and thus the enhanced fluid flow and poro-elastic coupling is incorporated using the exponential permeability–damage relation [55]:

$$\kappa \left(\alpha \right)={\kappa }_0 exp\left(D\alpha \right),$$

(5)

To account for the several orders of magnitude change in the permeability $\kappa \left(\alpha \right)$ between intact rock ($\kappa \left(0\right)={\kappa }_0$) and heavily damaged fault zone ($\alpha \to 1$), the exponential factor D should vary between 10 and 20 [57].

5.2. Numerical Results

Modeling of the response of the hydrological system to the remote earthquake was done by using a two-dimensional, vertical cross-section, implemented into FeFlow [58]—a finite-element simulator that solves the coupled variable density groundwater flow and solute transport equations. Layers’ properties were obtained from the pumping tests analysis, and the marl and chalk layer separating the two aquifers were assigned with lower hydraulic conductivity. The length of the numerical cross-section was set to 30 km to avoid the influence of the right boundary, and a mesh of 228,090 elements and 114,752 nodes was used to capture the geometrical boundaries of the model. Flow boundary conditions are hydraulic heads at the right boundary and no flow at the top, bottom, and left boundaries (Figure 7). The fault is set to 300 m away from Meizar 3 well. The increased damage of the fault by the passage of seismic waves is simulated by an abrupt increase of the hydraulic conductivity of the fault from 0.1 to 1 m/d. Healing and sealing of the fault reduces the hydraulic conductivity back towards its original value. The amount of time required for healing/sealing processes was evaluated according to this setup. We used different $C_1$, $C_2$ values (Equation (4)) and D = 15 (Equation (5)) as a representative value and run the groundwater modeling. The model was run for 100 days to represent the monitoring data between 24 January and 1 May (Figure 6).

Figure 8 shows the expected damage decrease under ε_cmp = 3 × 10⁻⁴ corresponding to 15 MPa confining pressure (~600 m depth) for 50 GPa bulk modulus, $C_2$= 0.1, 0.25, 0.5, and $C_1$= 0.05, 0.3, 0.4 1/s. The time scale of healing with these parameters is days. Equation (5) predicts the permeability reduction with time (Figure 3b), with these damage-rate parameters. The permeability reduction of one order of magnitude occurs within 0.75, 5, or 13 days depending on the healing parameters ($C_1$ and $C_2$).

The groundwater flow simulations, using these time-dependent fault zone permeability profiles, predict water level variations (Figure 9). Results of the numerical simulations show that the hydraulic head in the deeper and pressurized aquifer decreased due to the drainage of the groundwater through the fault while the hydraulic head at the shallower and less pressurized shallower aquifer increased. Simulated water levels are sensitive to the time required for healing/sealing of the fault. The fault healing parameters that best fit the observed water levels are C₁ = 0.3 s⁻¹ and C₂ = 0.25 (Figure 9). It is possible to estimate rate and state friction parameter—b with these damage-healing parameters. Using the obtained C₂ ~0.25 and $f_{ss}$ = 0.6–0.8 results with $b~C_2\times f_{ss}$ ~0.15–0.2. This value of the rate and state friction model parameter falls slightly above most of the reported values for hard rocks, but still meets the laboratory-based values in gouge-bearing faults. This is compatible with the hydro-geological settings, suggesting hot water flow from the deep aquifer through the activated fault zone.

6. Discussion

Fault zone permeability can be used as a proxy for fracturing and healing. Groundwater observations and modeling presented in this study demonstrate significant permeability increase (above two orders of magnitude) induced by remote Elazığ earthquake, and relatively fast recovery, with a duration of several days (Figure 9). There are many reports of permeability changes of regional aquifer following earthquakes that cause water level changes and increased spring discharge [59,60,61]. To the best of our knowledge, there is no field documentation of fault zones, showing rates of post-seismic permeability change at time-scale of hours to days. Contrary to the present study, previous observations reported relatively small permeability variations, representing the time window well after the earthquake (Figure 10). For example, Xue et al. [62] used the tidal response of water level in a deep borehole drilled into a ruptured fault zone to track permeability. They observed healing rates that reduced the permeability by a factor of 1.5 during a few months. Similar results were also reported by Shi and Wang, [63,64]. Injection experiments showed a decrease of 50% of the permeability in three years following the Kobe earthquake [65]. In Nojima fault, the permeability stabilized eight years after the occurrence of the earthquake [66]. Tidal analysis requires that the length of the recorded observations should be at least a few days for meaningful property monitoring. Pumping tests cannot be carried out continuously because, depending on the hydraulic properties, relatively long periods are required for recovery. Neither of the methods (tidal analysis and pumping tests) can monitor permeability at a frequency greater than days.

The gradual material recovery in laboratory experiments and the reported post-seismic velocity changes after strong earthquakes are often modeled with a logarithmic fit [19,20,67,68,69,70,71]. The logarithmic function predicts relatively fast healing rate immediately after the event with a gradual decrease in its rate. As noted by Qiu et al. [72], the logarithmic function suffers from an overshooting issue in curve fitting over long periods. Strong healing is not consistent with long-living wide fault zone damage persisted for hundreds or thousands of years with no significant healing. In order to resolve this contradiction, Hobiger et al. [73] suggested using an exponential function to fit the velocity variation curves. Qiu et al. [72] noted that the post-seismic recovery process could be described by either an exponential or logarithmic function.

Here we suggest that relatively short-term healing might be described by a Dieterich type logarithmic relation (Equation (1)) or damage decrease according to the damage-dependent kinetics (Equation (3)). Neglecting any time-dependent stress variations, the damage healing equation (Equation (3)) predicts logarithmic-in-time healing (Equation (4)) similar to the Dieterich equation (Equation (1)). Such constant stress conditions are typical for laboratory experiments and could be utilized for relatively short periods of time in field observations. The complete damage kinetic equation considers three-dimensional structure of the stress or strain tensor and includes the condition for the transition between damage growth and recovery [51]. Even under constant confining pressure, shear stress significantly decreases the rate of material recovery and may even lead to its repeated accumulation when the onset of damage condition is satisfied. Such a decrease in the healing rate leads to the damage stabilization schematically shown in Figure 10. In this case, we start with almost constant load during the initial stage leading to logarithmic healing and reduce the strain function leading to damage stabilization at the level between 0.2 and 0.3, which are values reported for the study area [43].

7. Conclusions

The water level in water wells in the Lower Yarmouk Gorge responded to the remote 2020 Elazığ Mw 6.8 earthquake. The response included a sustain change during the passage of the seismic waves and a post-seismic diffusive response to an activation of a nearby fault. Immediately after the passage of the seismic waves, pressure started to change in all wells. The high pressure in the deep Judea aquifer started to decrease and the lower pressure in the shallower Ghareb-Mishash aquifer started to increase. During the first 10 days, the water level in Meizar 1 well decreased by 40 cm and then recovered towards its pre-event value. In Meizar 3 well, the water level increased by 90 cm and then started to decay back to its pre-event value. The increased damage of the fault allowed for water flow from the deep-pressurized aquifer to the shallow aquifer. As a result, water level decreased in the deep aquifer and increased in the shallow aquifer. Numerical modeling suggest that the hydraulic conductivity of the fault decreased by an order of magnitude within three days. To the best of our knowledge, such fast healing/sealing was never observed by hydrological processes in the field.