Characterising thermal water circulation in fractured bedrock using a multidisciplinary approach: a case study of St. Gorman’s Well, Ireland

A hydrogeological conceptual model of the source, circulation pathways and temporal variation of a low-enthalpy thermal spring in a fractured limestone setting is derived from a multidisciplinary approach. St. Gorman’s Well is a thermal spring in east-central Ireland with a complex and variable temperature profile (maximum of 21.8 °C). Geophysical data from a three-dimensional(3D)audio-magnetotelluric(AMT) survey are combined with time-lapse hydrogeological data and information from a previously published hydrochemical analysis to investigate the operation of this intriguing hydrothermal system. Hydrochemical analysis and time-lapse measurements suggest that the thermal waters flow within the fractured limestones of the Carboniferous Dublin Basin at all times but display variability in discharge and temperature. The 3D electrical resistivity model of the subsurface revealed two prominent structures: (1) a NW-aligned faulted contact between two limestone lithologies; and (2) a dissolutionally enhanced, N-aligned, fault of probable Cenozoic age. The intersection of these two structures, which has allowed for karstification of the limestone bedrock, has created conduits facilitating the operation of relatively deep hydrothermal circulation (likely estimated depths between 240 and 1,000 m) within the limestone succession of the Dublin Basin. The results of this study support a hypothesis that the maximum temperature and simultaneous increased discharge observed at St. Gorman’s Well each winter is the result of rapid infiltration, heating and recirculation of meteoric waters within a structurally controlled hydrothermal circulation system. Supplementary Information The online version contains supplementary material available at 10.1007/s10040-021-02393-1.

positioned in a nearby barn. Water level measurements are presented in metres above an arbitrary datum; in this case the datum is the position of the logger downhole.

S1.1.1 Hydrochemical measurements and analysis
In the course of the IRETHERM project, samples were collected and analysed from the six thermal springs shown in Fig. 1c. These hydrochemical analyses are the subject of a detailed study in Blake et al. (2016). Geographical coordinates, geological setting, maximum temperatures and a brief description of each spring is provided here in Table S1. Data were recovered for analysis over five seasons to assess the temporal variation in the spring chemistry and to provide some seasonal overlap for a more robust analysis. The springs were sampled in July/August and October 2013, and in January, May and August 2014 (see Fig. 6 of main article). Temperature, electrical conductivity, and pH were recorded at each spring prior to sampling using a calibrated YSI 556 portable multi-probe and a Hanna HI 98130 Combo meter (use of both instruments facilitated cross-checking of results).
The samples were analysed for major and minor ions by ELS Ltd., Cork, Ireland, and also for a suite of 70 trace elements by Acme Ltd. (now trading as Bureau Veritas Ltd.), Vancouver, Canada.
Samples were collected using a low-flow peristaltic pump and an in-line 0.45 μm filter to obtain samples that were most representative of the formation water (see Henry, 2014). The method used was based upon the technical standards ASTM D6452-99 and ISO 5667-11, and the methods of Barcelona et al. (1994) and Puls and Barcelona (1996). Samples for trace element analysis were acidified to a pH of < 2 with trace metal grade nitric acid (John, 2000;Huang et al., 2013) prior to transport. All spring samples were duplicated for each round and a system of blanks using ultrapure water (from a Milli-Q® water purification system) was analysed by each laboratory for each round (as recommended by USGS, 2006). The blanks showed no contamination of the samples. The samples were stored in a cool box during transport and storage prior to shipment to the laboratories for analysis.
At ELS Ltd. the analyses for total concentrations of major ions were measured using inductively coupled plasma mass spectrometry (EM130 ICP-MS), total alkalinity was measured using a Titralab EW153, and sulfate and chloride concentrations were determined using EW154M-1 AQ2-UP2 EW015/016 Autoanalyser Spectrophotometry. At Acme Ltd. trace element concentrations were determined using ICP-MS. Details of the analytes, including limits of quantification (LOQ), are presented in Table S2.
As a check on the accuracy of the analytical results, ionic balance errors were calculated using PhreeqC (version 2.18) (Parkhurst and Appelo, 1999) with the minteq.dat database. The majority of samples had calculated errors below the recommended standard of ±5% (Freeze and Cherry, 1979), with 14 % of the samples having elevated errors of between ±5% and ±10%. All samples were retained for further analysis as the ionic balance error of 10% was deemed acceptable (e.g., Cloutier et al., 2008;King et al., 2014   Assuming all of the groundwater comes from a karst limestone aquifer in the Waulsortian Limestone Fm., and an effective porosity of 0.01 (based upon lower end of range of values for regional karstic aquifers), the thermal spring has a contributory area of 145 km 2 , which could feasibly be contained within the Blackwater (Longwood) sub-catchment (area of 181 km 2 ).

S2.1 Dimensionality analysis of AMT data
The dimensionality of the data was analysed by investigating the Z and T responses independently of each other. For the Z responses, the dimensionality analysis was performed by examining the phase tensors (Caldwell et al., 2004), which have the advantage of being unaffected by galvanic distortion of the electric fields. Figure S1 shows the calculated phase tensor for each frequency for each station, depicted as an ellipse. For a 1-D scenario the phase tensor will be represented by a circle, and for a 2-D case the phase tensor will be represented by a symmetrical ellipse, with the orientation of the major axis aligned either parallel or perpendicular to the regional geoelectrical strike direction. For 3-D cases the phase tensor will be non-symmetrical, necessitating the use of an additional angle, β, to characterise the tensor. In Figure S1, the ellipses representing 3-D conditions are coloured depending upon the magnitude of β normalized by the corresponding error, following the approach of Campanyà et al. (2016). All stations in Figure S1 show coloured ellipses for some frequencies, indicating 3-D conditions for the survey area. For the T responses, induction arrows (Schmucker, 1970) following the Parkinson criteria (i.e., the real arrows tend to point towards current concentrations in conductive anomalies (Jones, 1986)) were used. Figure S2 shows the induction arrows for each station and each frequency (station 33 has no induction arrows because the T data quality was poor for this station).
For a 1-D scenario the length of the induction arrows will be less than the threshold length of the assumed errors as there is no induced vertical magnetic field. For a 2-D scenario the induction arrows will point in the same or exactly opposite directions for all periods and stations. In a 3-D scenario, real and imaginary induction arrows will point in different (oblique) directions at any one frequency for any station (as can be seen in Figure S2). The results from Figures S1 and S2 indicate the existence of a 3-D scenario beneath the survey area.

S2.2 3-D inversion of AMT data
AMT data from 28 frequencies (excluding frequencies in the dead-band, particularly between 800 Hz and 2,000 Hz) were prepared for the inversion; these data were subsequently re-edited on a station-bystation basis to remove particularly noisy frequencies. The data were inverted using the ModEM 3-D inversion code (Egbert and Kelbert, 2012;Kelbert et al., 2014). The vertical magnetic transfer functions (T) were inverted alongside the four components of the impedance tensors (Z) to improve the resolution of the subsurface resistivity values (e.g., Siripunvaraporn and Egbert, 2009) No correction or compensation was applied to the data to account for galvanic distortion, which is a tractable problem in 2-D cases, but far less practicable in 3-D (see Jones, 2011). An examination of the apparent resistivity curves revealed no particular "problem areas" for galvanic distortion. As a 3-D modelling approach was used, with a fine parameterization in the uppermost part of the model, it was expected that the model would not be greatly affected by near-surface galvanic distortion effects at our target depths (e.g., Sasaki and Meju, 2006;Farquharson and Craven, 2009;Meqbel et al., 2014).
Galvanic distortion may affect the very shallowest layers of the model, but at depth, particularly beneath 100 m where every part of the model is sampled by numerous stations, the conductivity shows a smooth and consistent distribution and any unresolved features near to the surface appear to have been assimilated by the model. Hence, at depths greater than 100 m the effects of galvanic distortion in the model should be negligible. Also, the inversion of T alongside Z should decrease the susceptibility of the model to the effects of galvanic distortion. As T does not involve the electric field directly, it is not subject to the same galvanic distortion as Z. However, it can become distorted if the deflection of electric currents by in-phase electric fields alters the vertical magnetic fields (Booker, 2014). The resulting models do not show obvious artefacts (i.e., site-correlated model structures), which commonly indicate the presence of static shifts.

Figure S1
: Phase tensor dimensionality analysis using Z responses. White-grey colours indicate frequencies affected by the presence of 1-D or 2-D structures. Other colours represent frequencies affected by 3-D structures. Stations are arranged from W to E in three panels to correspond with the boxes outlined in the map of the survey area (see Fig. 2 in main article).

Figure S2
: Induction arrow dimensionality analysis using T responses, following the Parkinson criteria. Stations are arranged from W to E in three panels to correspond with the boxes outlined in the map of the survey area (see Fig. 2 in main article).

Figure S3
: (cont. overleaf) Representation of the data fit to the model responses from the final 3-D inversion of AMT data. Data from both components of T and all four components of Z are represented (real and imaginary parts). The stations are grouped to reflect the three boxes in the inset: the stations are arranged in order of their appearance from W to E. The coloured scale represents the difference between the data and the model response divided by the error for each frequency at each station.