Towards a method for forecasting earthquakes in Iceland using changes in groundwater chemistry

The study investigates whether changes in groundwater chemistry can forecast M ≥ 5 earthquakes in northern Iceland several months in advance. Building on a 2014 study that observed precursory chemical changes before two consecutive 2012–2013 earthquakes, the authors hypothesize that similar signals in the following decade (2014–2023) could have enabled forecasts within 4–6 months. This work places groundwater geochemistry in the broader context of decades of reported pre-seismic hydrological and geochemical anomalies and aims to quantify forecasting skill using sensitivity (TPR), false positive rate (FPR), and positive predictive value (PPV). The purpose is to evaluate the feasibility and potential utility of a site-specific, medium-term forecasting approach that complements longer-term probabilistic hazard models and short-term aftershock forecasts.

Prior research over more than 50 years has documented groundwater chemical and hydrologic anomalies before earthquakes, including changes in isotopes (e.g., δ2H) and elements (e.g., Na) and radon anomalies (e.g., Kobe 1995). Foundational studies discussed hydrologic responses to seismic waves and empirical estimates of preparation zone sizes (e.g., Dobrovolsky et al. 1979; Rikitake 1988). The authors’ earlier Icelandic work (Skelton et al., 2014) found statistically significant associations between groundwater isotopic/chemical changes and the 2012–2013 M ≥ 5 earthquakes near Húsavík (HA01 well). Subsequent studies linked Na peaks to non-stoichiometric dissolution of analcime along microfractures and documented water–rock interaction processes (e.g., analcime formation, zeolites) supported by petrographic evidence. Broader literature addresses microfracturing mechanics and hydrological responses (Gudmundsson et al., Anders et al.) and reanalyses of oscillatory pre-seismic chemical signals (e.g., Kobe), supporting the plausibility of geochemically detectable preparatory processes.

Study area and sampling: Groundwater from well HA01 at Hafralækur, penetrating basaltic rocks and sediments, was sampled weekly since October 2008. The well is cased to 35 m with inlets at 65, 82, and 96 m, discharging ~7.7 l/s. Groundwater is hot (73–76 °C) and alkaline (pH ~10.2 at 25 °C). Samples for stable isotopes and dissolved elements (Na, Ca, Si, K) were filtered (0.2 μm) into acid-washed HDPE bottles; HNO3 was added to elemental samples. Samples were cooled and shipped for analysis. Analytical methods: δ18O and δ2H were measured at Stockholm University by cavity ring-down spectroscopy (Los Gatos Research LWIA), normalized to VSMOW-SLAP differences (−55.5‰ for δ18O; −428‰ for δ2H). Element concentrations were determined by ICP-OES (Stockholm University: Varian Vista AX until 2011; Thermo ICAP 6500 in 2012–2013; University of Iceland: Spectro Ciros Vision thereafter). Systematic shifts (Si, Na, Ca, K) due to silicate precipitation on glass nebulizers were mitigated at the University of Iceland using a geothermal-water protocol (short runs, Teflon nebulizer, frequent HF cleaning). Archived samples (four per year) were reanalyzed to correct earlier data. Parameters and time series: In addition to δ18O and δ2H, two composite parameters were computed to separate effects of water–rock interaction and meteoric mixing: - d = δ2H − 8×δ18O (deuterium excess; proxy for water–rock interaction; more negative implies greater interaction or certain climatic effects). - d′ = δ2H + (δ18O/8) (proxy for mixing between Holocene and pre-Holocene meteoric sources; more negative implies larger pre-Holocene component). Cumulative means and ±2σ envelopes were constructed from 1 Oct 2008 start to each measurement date. Anomalies were defined as deviations ≥2σ from the cumulative mean. Earthquake data and derived fields: Earthquake locations, magnitudes, and focal mechanisms were obtained from USGS NEIC and Global CMT; nodal planes matched onshore fault geometries. Tapered faults were modeled in Coulomb 3.3 to estimate dilational/compressional strains and seismic energy densities at HA01. Local Icelandic Meteorological Office locations were compared; maximum discrepancy (12 km) was too small to affect event selection. Magnitudes reported are surface-wave M. Statistical testing (association): A binomial test compared the occurrence of earthquakes within 4–6 months of groundwater chemistry changes against a null of random timing, using p = tf/t where t is study duration (63 months in 2014 study; 172 months in present study) and tf the 4–6 month forecast window. p-values were computed for the original (2008–2013) and extended (2008–2023) series. Forecast analysis (skill metrics): For 2014–2023, forecasts were defined as occurrences where d′ deviated ≥2σ for 1, 2, or 3 consecutive measurements. Forecast windows were 4, 5, or 6 months. Skill metrics were computed: - TPR = TP/(TP+FN) - FPR = FP/(TN+FP) - PPV = TP/(TP+FP) TP (anomalies followed by an earthquake within the window), FP (anomalies not followed by an earthquake within the window), TN (non-anomalous points not followed by an earthquake), FN (non-anomalous points followed by an earthquake). Counts and rates are reported in Table 2. Spectral analysis: FFT (Excel) on pre-event windows (441 days before 2020; 889 days before 2012, detrended) assessed periodicities in d and d′. Element concentrations: Trends in Na, Ca, Si, K were assessed relative to d and d′ behavior, with statistical tests for anomaly significance at earthquake times.

- Forecastability: Of the M ≥ 5 earthquakes affecting northern Iceland (2014–2023), only the strongest (M 6.0, June 2020 sequence) could have been forecast based on statistically significant groundwater chemistry changes 4–6 months prior; the M 5.0 (2018) and M 5.2 (2022) events were not forecastable by this criterion. - Skill metrics (2014–2023, using d′ anomalies ≥2σ): • TPR (sensitivity): 0.28–0.32 (4 months), 0.24–0.30 (5 months), 0.20–0.27 (6 months). • FPR: 0.01–0.03 (4 months), 0.01–0.02 (5 months), 0.01–0.02 (6 months). • PPV: 0.62–0.80 (4 months), 0.72–0.85 (5 months), 0.76–0.85 (6 months). Representative Table 2 counts include (examples): for 4 months window, 1, 2, 3 consecutive anomalies yield TP/FN/FP/TN of (18/39/11/417), (17/40/6/422), (16/41/4/424), respectively. - Statistical association: Binomial tests reject the null of random coincidence of earthquakes with groundwater changes: p = 0.02–0.04 for 2/2 events (2012, 2013) in 2008–2013; p = 0.003–0.01 for ≥3/5 events (2012, 2013, 2020) in 2008–2023. - Timing of anomaly: For d′ and δ18O, oscillatory maxima/minima exceeded 2σ about 5 months before the 2020 event. - Oscillatory behavior: Spectral analysis indicates an ~75-day periodicity in d and d′ before the 2020 earthquake (and similar, less well-defined periodicity before 2012–2013), interpreted as cyclic expansion–contraction of the groundwater source region. - Element concentrations: While Na showed statistically significant maxima coincident with the 2012 and 2013 events (attributed to non-stoichiometric analcime dissolution along microfractures), apparent anomalies around 2020 were not statistically separable from background; thus major elements were not useful for forecasting the 2020 event. - Mechanical context: Co-seismic dilational/compressional strains at HA01 were small (10−10 to 10−7), but seismic energy densities (10−1 to 101 J m−3) were within ranges associated with co-seismic groundwater responses elsewhere. All five major events considered occurred within an empirically inferred region where precursory groundwater changes might be detectable at HA01.

The results support the hypothesis that groundwater chemistry changes at HA01 can, in some cases, precede significant regional earthquakes by several months, enabling a site-specific forecasting method with moderate PPV and low FPR but limited sensitivity. The approach complements other probabilistic models: it targets a medium-term window (4–6 months), potentially bridging long-term hazard forecasts and short-term aftershock probabilities. Underlying mechanisms likely involve microfracturing associated with crustal dilation during stress build-up. Microfractures can strongly influence groundwater chemistry due to high fracture densities at small scales, promoting mixing between distinct meteoric sources (captured by d′) and affecting water–rock interaction (d). Oscillatory anomalies with ~75-day periodicity suggest cyclic expansion and contraction of the groundwater source region, potentially coupled with fracture mineralization processes. Na maxima coincident with the 2012–2013 events and petrographic/geochemical evidence support microfracture-related processes (e.g., analcime dissolution). The less pronounced Na signal in 2020 may reflect prior depletion or preferential re-opening of existing microfractures. Collectively, these findings indicate that while not all M ≥ 5 events are forecastable using this method, when anomalies occur they are informative (high PPV) and rarely produce false alarms (low FPR). Implementation at other sites would require prolonged baseline monitoring and suitable hydrogeologic and lithologic conditions to resolve mixing and interaction signals.

The study advances a site-specific method for forecasting M ≥ 5 earthquakes in northern Iceland using groundwater isotopic chemistry, particularly the d′ parameter, with a 4–6 month forecast window. Retrospective analysis indicates: sensitivity of 20–32%, low false positive rates (1–3%), and high positive predictive values (62–85%). The method could have forecast the 2020 M 6.0 event but not the weaker 2018 and 2022 events. Observed oscillations and periodicities support a mechanism of microfracturing, mixing, and mineralization within the groundwater source region during stress accumulation. Future work should involve prospective testing to validate performance, expansion to additional sites with long-term monitoring, refinement of anomaly definitions and thresholds, integration with seismic and geodetic observations, and further mechanistic studies linking hydrochemistry to crustal stress evolution.

- A posteriori retrospective design; requires prospective validation to confirm forecasting skill and avoid overfitting. - Site-specific results (HA01, northern Iceland); generalizability depends on local hydrogeology, isolation from surface influences, and distinct mixing end-members. - Limited event sample size (five M ≥ 5 sequences over 14 years), which constrains statistical power and sensitivity estimates. - Sensitivity is modest (20–32%), implying a high probability of missed events (68–80%). - Major element concentration anomalies (e.g., Na, Si, Ca) were not statistically significant for the 2020 event; reliance on isotopic parameters may limit robustness where isotopic signals are noisy. - Potential climatic influences on deuterium excess (d) require caution in interpretation despite use of d′ to mitigate confounding. - Dependence on continuous, long-term, high-quality sampling and analysis; instrumental biases required post hoc corrections in earlier elemental data. - Earthquake location uncertainties are mostly small but could affect strain/energy estimates; magnitudes are reported as surface wave M but conclusions are intended to be magnitude-scale agnostic.