Earthquake Swarms Triggered by Groundwater Extraction Near the Dead Sea Fault

In 2013 and 2018, earthquake swarms with a maximum moment magnitude of 4.5 occurred ~5 km from the northern section of the Dead Sea Transform Fault. Here we show that aquifer pressure data, interferometric synthetic aperture radar surface deformation time series, and seismic monitoring suggest that groundwater withdrawal triggered these earthquakes. Continuous groundwater extraction from several wells located ~10 km west of the swarms has accelerated since 2010 and resulted in a total decrease of ~50 m of the groundwater level at the time of the 2018 earthquake swarm. The withdrawal also corresponds to surface subsidence of ~10 mm/year based on repeat interferometric synthetic aperture radar measurements. The temporal correlation, extensive subsidence, anomalous swarm characteristics, and normal faulting orientation suggest a connection between the groundwater withdrawal and recent earthquakes. Poroelastic modeling demonstrates that pumping‐induced pore pressure decrease west of the earthquake could have caused significant dilatational stresses that led to normal faulting events outside the aquifer.


Introduction
The Lake Kinneret (Sea of Galilee) basin is located on the northern section of the Dead Sea Transform (DST) separating the Arabian plate and the Sinai subplate ( Figure 1) and has accommodated~100 km of sinistral slip motion since the Miocene . The unique archeological (Ellenblum et al., 1998(Ellenblum et al., , 2015 and geological record in Northern Israel (Marco et al., 2005;Wechsler et al., 2014) indicates that this segment of the Dead Sea Tranform has ruptured several times since the Iron Age (980 BCE). The archeological records allow a relatively accurate determination of timing and displacements of earthquake rupture, including a left lateral strike-slip of 1.6 m of slip in 1202 and a 0.5 m of slip in 1759. From the beginning of the instrumental record in 1985 until 2013, seismicity in the lake was limited to sporadic events with only one M W ≥ 3 earthquake every few years (Figures 1 and 2) and most of these earthquakes could be associated with the main transform fault ( Figure S1 in the supporting information).
The situation changed suddenly in 2013 when two seismic swarms affected Lake Kinneret. The first began in September 2013 and included five M W > 3 earthquakes and the second occurred in July and August 2018 with 15 M W > 3 earthquakes (Figures 1 and 2). The 2013 and 2018 seismic swarms are approximately collocated. In both cases, the majority of the seismicity was located in the northwest section of the lake and earthquakes of magnitude M W > 3 were localized along NNW-SSE strands with focal depths of only 0-10 km in contrast to the deeper seismicity north of the lake ( Figure S1). A closer inspection shows that the 2013 swarm is slightly to the west of the 2018 swarm (see Text S1 in the supporting information for location method details).
In this paper we explore a possible relationship between the decline of groundwater level within the aquifers and the seismic activity. Multi-well aquifer pressure monitoring, surface deformation time series based on satellite imagery, and detailed seismic monitoring capabilities are all used to constrain pumping-induced deformation in this region and its relationship to the earthquakes. We use the data to build on the temporal correspondence of the earthquakes noted above with three additional observational lines of evidence suggesting that the swarms are human induced. We also combine the information into a poroelastic model that demonstrates a plausible relationship between groundwater withdrawal and earthquake swarm generation next to the Dead Sea Transform fault.

10.1029/2019GL083491
Key Points: • Two anomalously shallow earthquake swarms occurred below Lake Kinneret following periods of rapid groundwater extraction • Water level decrease and extensive ground subsidence captured by InSAR are correlated with the earthquake swarms • Poroelastic modeling supports earthquake triggering by groundwater extraction Supporting Information: • Supporting Information S1   (Sneh & Weinberger, 2014), earthquake catalog between 1985 and September 2018. Black arrows represent the relative plate motion along the main Dead Sea Tranform. The moment tensor solutions for the two M W ≥ 4 events are plotted (see more information in the supporting information). Shaded relief by Hall (1994). (c) Focal mechanism solutions based on full moment tensor inversion, Time Domain Moment Tensor, are presented for selected events of the 2013 (orange) and 2018 (blue) swarms. Focal mechanism of previous studies (Hofstetter et al., 2007;Van-Eck & Hofstetter, 1990) are presented in gray. (d) The temporal occurrence of earthquakes along Lake Kinneret is shown for events of M L ≥ 3 by the colored dots. Lower magnitude earthquakes are plotted in gray. DST = Dead Sea Transform.

Hydrological Information of the Aquifer
The timing of the earthquake swarms is particularly interesting in light of the hydrological history of the region. Lake Kinneret used to be the source for one third of the domestic water of Israel. The increase in demand for water for domestic uses and the continuous decline in rain have changed the water resources management. The sharp decline in the lake level and the concern that the salinity of the lake might increase to very high levels (Markel & Shamir, 2002) has resulted in a sharp drop in the pumping rates from the lake from a level of 500 MCM/year (million cubic meters per year, http://www.water.gov.il) in the 1990s to 60 MCM/year in recent years ( Figure S3). An alternative resource for the water was groundwater around the lake and especially to the northwest of the lake. Increasing pumping rates in the regional wells (Figures 2a and S4) caused the groundwater level to drop by more than 50 m (0.5 MPa) since the 1990s. During that period water level is monitored in the wells 3-12 times per year. The measurements in the wells are performed hours after pumping was paused to allow the recovery of the dynamic levels during pumping to static levels. The actual time required to reach the static level changes between the wells depends on the hydraulic properties and pumping rates. The measured static water levels are used to construct water level maps that are used to manage the water budget at the aquifers.
The two periods of fastest groundwater decline rate were 2007-2013 and 2016-2018. These two periods both correspond to earthquake swarms at the lake, which were located 2-10 km from the pumping wells.

Abnormal Earthquake Sequences
The first line of evidence for the unusual nature of the earthquakes comes from the sense of motion. Focal mechanism solutions of the 2013 and 2018 clusters (listed in Table S1) show that most of the seismic events are dominated by normal faulting striking NW to NNW ( Figure 1c). We solve the focal mechanism solutions for all earthquakes of M w ≥ 3.4 using the full-waveform Time Domain Moment Tensor technique by Dreger and Helmberger (1993) implemented at Geological Survey of Israel. Events with magnitude below M W 3.4 resulted in poor convergence of the inversion, due to low signal-to-noise at the stations. The algorithm inverts three component waveforms to estimate moment tensor in a point-source approximation and decomposed into the scalar seismic moment, a double-couple moment tensor and a compensated linear vector dipole moment tensor. The isotropic component in this study is set to be zero. Synthetic seismograms are represented as linear combination of 10 fundamental Green's functions that form the basis functions for any arbitrarily oriented double-couple moment tensor mechanism. In this study, the Green's functions are computed by using the frequency-wavenumber integration code (FKPROG) of Saikia (1994) for six-layered, regionally calibrated velocity structures of Gitterman et al. (2005). We calculate the Green's functions at 5-km horizontal intervals (from 5 to 400 km) and at 1-km vertical intervals between 1 to 40 km. Focal depth is constrained from the inversion run over various depths searching for the maximum in variance reduction annual pumping rate at the regional wells (black dots), and the yearly rate of earthquakes above magnitude 3 (pink histograms). The pumping rate of the individual wells is shown in Figure S4. Two arrows plotted on the x axis are pointing to the increase of seismicity during the two swarms. (b) The b value calculated for seismicity before and after the first swarm of 2013 (see supporting information section for more information, Figure S13).
(fit between the waveforms and the synthetics). We use velocity waveforms from broadband and shortperiod three component stations provided by the Israeli Seismic Network ( Figure S2). Data are extracted in 100-s windows starting 60 s before the event origin time, corrected for instrument response and the horizontal components are rotated to great circle path. Finally, data and Green's functions are filtered by applying a band-pass filter in the frequency band 0.05-0.1 Hz for earthquakes with M W < 4 ( Figure S5), and 0.02-0.1 Hz for earthquake with magnitude 4 ≤ M W ≤ 5 ( Figure S6), in order to capture the low frequency seismic energy content (Scognamiglio et al., 2009). The best result is achieved by performing a grid search on the depth, and choosing the moment tensor solution and centroid depth for which the variance reduction is at maximum.
The local stress field that corresponds to these normal faulting mechanisms (Figure 1b) is distinct from earlier focal mechanisms and the dominant sense of motion on the Dead Sea Fault (Palano et al., 2013). The discrepancy might suggest a local masking of the regional stress field by a local and more dominant stress field. There is a lack of agreement about the fault architecture of the basin, and we suggest that faulting occurred on either one of the subparallel east dipping (60°± 15°) faults (e.g., Rosenthal et al., 2019), or unmapped low-angle west dipping fault planes (30°± 12°).
By itself, the fault orientation would not qualify as a major line of evidence for human involvement in the earthquake triggering, however, the normal faulting is suggestive in context of the other two lines of evidence discussed below. It is also noteworthy that the sense of motion is consistent with that expected from withdrawal as will be shown below by the modeling results.
Additional evidence for the unusual nature of the seismicity comes from its swarm-like character. Swarms are earthquakes sequences that lack the usual mainshock-aftershock behavior and instead have a disproportionate number of small earthquakes. Quantitatively, swarms can be distinguished by anomalous Gutenberg-Richter b values, which are used to measure the relative number of small to large earthquakes. The swarms decreased the Gutenberg-Richter b from 1.0 to 0.7 (Figure 2b), compatible with swarm-like induced seismicity (Goebel et al., 2016;Shelly et al., 2016).

Surface Deformation Response to Groundwater Pumping
An additional line of evidence is subsidence that correlates temporally and spatially to the groundwater withdrawal, thus providing direct evidence of a significant deformation of the solid Earth. We produced a time series of deformation from 11 October 2014 to 17 February 2019 using Sentinel-1 SAR data acquired on the descending (Track 21), orbits processed using the Sentinel stack processor of Jet Propulsion Laboratory's ISCE software (Fattahi et al., 2017). We take 5 and 15 looks in azimuth and range direction, respectively. For each acquisition we form four sequential interferograms and use the routine workflow of the MintPy InSAR time-series software package to invert for the displacement time series. Reference point is positioned at the town Kafar Kana (32.7458°N, 35.3436°E), located on the Mount Scopus Group (Senonian Paleocene). The horizontal north-south tectonic displacements associated with the Dead Sea fault system larger scale than the area considered and are unlikely to result in the local, relative line-of-sight changes. The expected changes are so small in part because north-south trending faults are difficult if not impossible to study with polar-orbiting InSAR satellites because of the viewing geometry. From the InSAR time series of each cell, we define the behavior of surface vertical motion by fitting displacement time series to a parabolic function capturing the accelerated subsidence rate of the surface. The displacement time series are fitted using the MATLAB@ nonlinear curve-fitting (data-fitting) solver in least squares sense. We solve the regression for the parabola and sinusoidal oscillation with native seasonal frequency of one cycle per year searching for four free parameters: acceleration (a 1 ), U 0 (a 2 ), gain (a 3 ), and phase (a 4 ; Figure S7), where t is the time in years. To keep the vertex of the parabola at the beginning of the time series, we do not use the linear term. In general, the subsidence field presents an area that coincides with the location of the wells (Figure 3) with the seasonal fluctuations corresponding with the mean regional monthly precipitation reported by the Israel Meteorological Service (http://www.ims.gov.il/IMSENG/All_Tahazit/homepage. htm). The subsided area is west of the swarm with the strong subsidence in a region elongated east-west. Minor uplift is visible at regions without wells and any known withdrawal. Several regions indicate accelerated subsidence with~3 mm/year 2 with mean root-mean-square of 25 mm. The accelerated subsidence can be translated to mean of −12 mm/year or a total of~50 mm over the four and a half years of our time frame.

Conceptual Model of Induced Earthquakes
In most of the cases of fluid-induced seismicity, the triggering mechanism involves hydraulic connection between the well and the fault with an increase in pore pressure reducing the frictional strength and promotes failure at the regional ambient stress (Raleigh et al., 1976). In the present case, groundwater pumping results in a decrease in pore pressure and so an alternative model is required.
To obtain more insight into potential triggering mechanisms, we calculated production-induced stress changes in a 2D-layered poroelastic model. The model is informed by the extensive hydrogeological information available for the region. The three main aquifers west of the lake include the Cretaceous Judea group, the Lower Cretaceous Kurnub group, and the Jurassic Arad group with a total depth of 4 km ( Figure S8a). All of these aquifers are hydraulically connected and were shown to affect the salinity of the springs around the lake (Gvirtzman et al., 1997). All of the groundwater is extracted from the shallowest Cretaceous Judea aquifer. The basin fill sediments (fluvial and lacustrine sediments) beneath the lake are 5-8 km thick and have very low permeability (Ben-Avraham et al., 1996Hurwitz et al., 2000;Marcus & Slager, 1985).
The model ( Figure S8b) consists of five geological units including upper and lower Cretaceous formations, as well as Eocene, Jurassic, and Triassic sediments (Table S2). The different geological units are hydraulically connected so that ground water is expected to migrate horizontally and vertically during fluid pumping. The x") are shown at the six panels below (blue dots), in which the gray lines are the fitted seasonal curves, and the dashed lines are the 95% residuals, and the black solid line is the parabolic acceleration without the sinusoidal parameters. The mean monthly precipitation for the region is shown (Panel 1), and water level at selected well are shown for the nearby reference point (Panels 2 and 3). base of the model represents the top basement that is composed by impermeable granitic rocks. We built the model using the Comsol Multiphysics finite-element package (COMSOL, 2017). We used an adaptive mesh to accurately resolve poroelastic stress changes close to the ground water wells and across the aquifer boundaries. Element sizes vary between 50 m close to the well and across the high-low permeability transition at the edge of the aquifer and 500 m at the center of the aquifer. We benchmark our model ( Figures S9 and S10) using analytical solutions from Helm, 1994 for shallow pumping in a vertically confined aquifer (Helm, 1994) and through comparison with analytical solutions for poroelastic compaction in a layered aquifer model (Leake & Hsieh, 1997).
Changes in pore fluid pressure and solid elastic stresses are fully coupled via Comsol's poroelasticity module and boundary conditions and model geometry are shown in Figure S8b. We prescribe the measured drop in hydraulic head at a pumping well of about 50 m on the left boundary within upper and lower Cretaceous at a depth of 1,000 m from 2007 to 2017 ( Figure 4). For this simplified model, we assume a constant rate of head decrease during this 10-year period that encompasses the time of most rapid groundwater extraction. The fluid production zones are bounded to the east by the impermeable basin fill.

Discussion and Conclusions
The simplified poroelastic model predicts subsidence rates within the range of 7 to 100 mm/year at different distances from the well ( Figure S11a), which is roughly comparable to the observed subsidence from InSAR measurement (Figures S11b and S11c). Groundwater pumping in our model results in large subsidence due to poroelastic coupling, where the largest subsidence is predicted close to the pumping wells as a result of production-induced pore pressure drop (Figures 4and S11a). At the edge of the aquifer, the model results suggest large pressure and displacement gradients across the bounding low permeability sediment-fill underneath the lake. These gradients lead to dilatational horizontal stresses of up to 0.1 MPa, which might promote failure in the form of normal faulting events on preexisting faults. We perform several tests to evaluate the robustness of these observations. These tests include varying the hydraulic conductivity values outside the aquifer ( Figure S12), changing far-field boundary conditions and fault dip angle between 60°and 90°. We find that displacement gradients and horizontal strains are generally present across permeability boundaries at the edge of the aquifer.
Previous studies suggested that anthropogenic stress changes of~0.1 MPa are sufficient to trigger observable earthquakes with typical regional seismic coverage (e.g., Keranen et al., 2014;Sumy et al., 2014). Similar thresholds are commonly suggested for naturally triggered earthquakes (King et al., 1994;Saar & Manga, 2003). Other studies indicate that there is no physical lower limit to triggering stresses as faults are distributed at different points in their loading cycles (Van Der Elst & Brodsky, 2010). An observed triggering threshold is best thought of as corresponding to sufficient perturbation to produce a large enough earthquake rate change to be observable on a given network (Brodsky & van der Elst, 2014).
Although there are many examples of subsidence due to excess withdrawal of groundwater (e.g., Bawden et al., 2001;Cabral-Cano et al., 2008), there are few cases where groundwater extraction has triggered earthquakes (Amos et al., 2014;González et al., 2012). Previous studies proposed that earthquakes triggered by fluid extraction can be caused by crustal unloading and stress transfer processes (González et al., 2012). Crustal unloading from groundwater extraction would produce thrust faulting in the geometry studied here, which is at odds with the observed focal mechanisms.
Here we use poroelastic models similar to those previously utilized for evaluating energy extraction activities (Segall, 1989;Segall & Fitzgerald, 1998) and evidence for a more localized effect. Modeling results suggest that differential subsidence due to a strong lateral contrast between the aquifer and the impermeable basin sediments can significantly amplify extraction-induced stress changes. Such Horizontal stresses at depth of 2.5 km. Relatively high dilation is generated at the impermeable basin sediments under the lake. stress changes extend beyond the produced aquifer and can induce seismic events at larger distances with magnitudes that are not constrained by the size of the aquifer. Note that other mechanisms for propagating fluid pressures to the far field such as utilizing the fault damage zones as hydraulic conduits (e.g., Ortiz et al., 2019) are not appropriate for the pore pressure decreases driving the seismicity here.
The aquifer, deformation and seismic data analyzed here imply that groundwater withdrawal is producing earthquake swarms along the northern Dead Sea Tranform fault. Continuous groundwater extraction accelerated since 2010 and resulted in a total water level decrease of~50 m at the time of the 2018 earthquake swarm. Poroelastic stresses appear to extend the impact of the human activities well beyond the length of the aquifer. Individual earthquakes have been as large as magnitude 4.5. The geological history of the fault combined with the close proximity to populated areas suggest that future pumping in the region should be closely monitored.