TY - JOUR AU - Pasquali,, Riccardo AB - SUMMARY We present modelling of the geophysical data from the Newcastle area, west of Dublin, Ireland within the framework of the IRETHERM project. IRETHERM's overarching objective was to facilitate a more thorough strategic understanding of Ireland's geothermal energy potential through integrated modelling of new and existing geophysical, geochemical and geological data. The Newcastle area, one of the target localities, is situated at the southern margin of the Dublin Basin, close to the largest conurbation on the island of Ireland in the City of Dublin and surrounds. As part of IRETHERM, magnetotelluric (MT) soundings were carried out in the highly urbanized Dublin suburb in 2011 and 2012, and a description of MT data acquisition, processing methods, multidimensional geoelectrical models and porosity modelling with other geophysical data are presented. The MT time-series were heavily noise-contaminated and distorted due to electromagnetic noise from nearby industry and Dublin City tram/railway systems. Time-series processing was performed using several modern robust codes to obtain reasonably reliable and interpretable MT impedance and geomagnetic transfer function ‘tipper’ estimates at most of the survey locations. The most ‘quiet’ 3-hr subsets of data during the night time, when the DC ‘LUAS’ tram system was not operating, were used in multisite and multivariate processing. The final 2-D models underwent examination using a stability technique, and the final two 2-D profiles, with reliability estimations expressed through conductance and resistivity, were derived. In the final stage of this study, 3-D modelling of all MT data in the Newcastle area was also undertaken. Comparison of the MT models and their interpretation with existing seismic profiles in the area reveals that the Blackrock–Newcastle Fault (BNF) zone is visible in the models as a conductive feature down to depths of 4 km. The investigated area below Newcastle can be divided into two domains of different depths, formed as depth zones. The first zone, from the surface down to 1–2 km, is dominated by NE–SW oriented conductors connected with shallow faults or folds probably filled with less saline waters. The conductors are also crossing the surface trace of the BNF. The second depth domain can be identified from depths of 2–4 km, where structures are oriented along the BNF and the observed conductivity is lower. The deeper conductive layers are interpreted as geothermal-fluid-bearing rocks. Porosity and permeability estimations from the lithological borehole logs indicate the geothermal potential of the bedrock, to deliver warm water to the surface. The fluid permeability estimation, based on Archie's law for porous structures and synthetic studies of fractured zones, suggests a permeability in the range 100 mD–100 D in the study area, which is prospective for geothermal energy exploitation. Electrical properties, Fracture and flow, Permeability and porosity, Hydrothermal systems, Magnetotellurics, Numerical modelling 1 INTRODUCTION Geothermal energy as a clean and renewable resource can be used for electricity production or direct heating purposes (Olasolo et al. 2016), and geophysics has played a key role in the description of geothermal systems and estimation of their potential for energy utilization (e.g. Domra Kana et al. 2015). Within the framework of Ireland's IRETHERM project (IREland's geoTHERMal potential), geophysical and geological studies were carried out to develop an understanding of Ireland's low enthalpy geothermal energy potential through joint-inversion/integrated modelling of new and existing data (Jones et al. 2015; Blake et al. 2016). IRETHERM comprises three broad geothermal target types, Assessment of the geothermal energy potential of Ireland's radiogenic granites. Assessment of the geothermal energy potential of Ireland's deep sedimentary basins (Hot Sedimentary Aquifer—HSA). Assessment of the geothermal energy potential of warm springs. Linking these three together was subsurface imaging using both controlled-source (CSEM) and natural-source (magnetotellurics, MT) electromagnetic methods. Electrical conductivity, being a transport property, is a proxy for permeability and appropriate porosity-permeability relations are being developed. One of the chosen target localities for HSA investigation was the Newcastle area, situated at the southern margin of the Dublin Basin in the Greater Dublin Area, due to its proximity to Dublin and its significantly elevated geothermal gradient (over 30 °C km-1) with a high potential for district space heating of nearby urban areas (Fig. 1), and possibly even electricity generation (Goodman et al. 2004). The potential for geothermal energy development lies in permeable faults deep within the Dublin Basin below Newcastle that bring warm fluids to shallower depths (GT Energy 2009). The regional structural conditions provide large scale fractures (O'Neill & Pasquali 2005) extending over a significant distance along the Blackrock–Newcastle Fault (BNF). Figure 1. Open in new tabDownload slide Location of the MT stations, at which measurements were performed in 2011 and 2012, and wells in Newcastle area. Green lines represent railways tracks with DC traction, which are one of the main sources of electromagnetic noise in the area. (Inset) Survey areas in Ireland to investigate geothermal targets overlain on modelled temperatures at 2.5 km depth (Goodman et al. 2004). Our study is marked by the black square in the east. Figure 1. Open in new tabDownload slide Location of the MT stations, at which measurements were performed in 2011 and 2012, and wells in Newcastle area. Green lines represent railways tracks with DC traction, which are one of the main sources of electromagnetic noise in the area. (Inset) Survey areas in Ireland to investigate geothermal targets overlain on modelled temperatures at 2.5 km depth (Goodman et al. 2004). Our study is marked by the black square in the east. The identification of structures with higher primary porosities and permeability, which is one of the requirements for reservoir productivity, is the main objective for geothermal energy potential estimation in any prospective geothermal basin area. Particularly for low-enthalpy geothermal hydrothermal systems, where a large amount of fluid is required, such as in the Dublin Basin, the focus is made on the porous sandstones and limestones, the basal formations of sedimentary successions (Jones et al. 2015). The areas along major and minor crustal faults and shear zones associated with the main BNF, where the secondary porosity is high, can allow sufficient fluid flow circulation and increase geothermal source productivity. Knowledge of the BNF geometry within the Carboniferous sedimentary strata plays an important role in understanding the circulation system. This information can help to answer the overarching question; whether the presented lithologies and their depositional environment in this area are suitable to deliver sustainable geothermal heat for industrial and/or residential utilization or not (Muffler & Cataldi 1977). The geophysical MT surveying method plays an important role in the understanding of complex geothermal reservoirs properties (Asaue et al. 2006; Newman et al. 2008; Muñoz 2014). The fluids and permeability of rocks are the key factors for the transport of heat within the Earth (i.e. important characteristics of geothermal systems), and horizontal and vertical electrical conductivity variations in sedimentary rocks represent information about ionic transport properties within the basin sequences. Electrical conductivity is sensitive to fluid content and interconnectivity and to the movement of ions. Thus, saline-filled lithological units with greater porosity and permeability exhibit higher conductivity (lower resistivity). Of all physical properties, electrical conductivity is the most suitable proxy for porosity/permeability determinations observable from the surface and for understanding processes in the geothermal system. Including additional information about vibrational mechanical medium properties, provided by seismic data, and bulk density properties, aids constraining features and allows mapping of the system for structural and lithological development, such as faults, units thickness variation and fractures. There are few previous surveys of the BNF and Dublin Basin, and they are represented primarily by seismic investigations and two boreholes in the area (Strogen & Somerville 1984; Jones et al. 1988). The most important information about rock properties and geophysical parameters come from two relatively-recently drilled boreholes, NGE1 and NGE2 completed by GT Energy Ltd (GT Energy 2009). The geophysical surveys include Vertical Seismic Profiling (VSP), 2-D seismic reflection profiles and a regional gravity survey. Active seismic data were acquired with Vibroseis sources along two lines targeting depth to a maximum of 4 km, and final migrated sections were produced. Gravity measurements were acquired to map the distribution of sedimentary, meta-sedimentary and igneous rocks within and adjacent to the Dublin Basin (Readman et al. 1997; GT Energy 2009), and DC resistivity profiling was performed, but only to shallow depth of 100 m. In the most recent investigation of the southern margin of the Dublin Basin and the BNF, teleseismic receiver functions (RFs) were used to investigate the seismic velocities, anisotropy and the open cracks filled by fluids (Licciardi & Agostinetti 2017). In this study, we present multidimensional MT modelling of the IRETHERM MT data collected in the Newcastle area during two campaigns. The investigated area is very difficult for MT data acquisition due to the contamination of the MT natural signal by artificial electromagnetic noise from industrial parks, the DC railway system and the nearby Dublin city itself. Remote reference analyses of MT data were applied as standard technique (Gamble et al. 1979), which was included in different processing codes based on robust statistics. There are several successful case studies published on long term problem of noisy data (Jones et al. 1989) and related reviews (Szarka 1987; Junge 1996). The final data are modelled on 2-D profiles crossing the BNF and by 3-D inversion codes, and are compared with seismic models and borehole information from GT Energy survey. The final geoelectrical images are used for evaluation of low-enthalpy geothermal resources, which are suitable for the doublet system (Mottaghy et al. 2011) in the Newcastle area. In particular, we focus on the estimation of depth, dimension and permeability of fractured zones along the southern margin of the Dublin Basin. 2 GEOLOGICAL BACKGROUND The fault-bounded Dublin Basin (Fig. 1) is a Carboniferous basin comprising rocks of the Calp Limestone Formation deposited in deeper water. The limestones were deposited as a thick sequence by turbidity currents and are interbedded with conglomerates containing clasts of granite and Lower Palaeozoic rocks (Holland & Sanders 2009). Our study is focused on the southern margin of the Dublin Basin and its bounding fault structure near the adjoining Newcastle area. This part of the structure is known as the BNF line (Fig. 1), which separates the Carboniferous domains of the Dublin Basin and Kildare Shelf to the north and west, and the Lower Palaeozoic Leinster Massif to the south (McConnell & Philcox 1994). The area to the south of the BNF is characterized by Silurian rocks. The group is subdivided into several subgroups and consists of greywackes and shales, deposited primarily as turbidites (Brück et al. 1974). Further to the east, towards the Leinster Granites of the Wicklow Mountains, the area is characterized by the Aghfarrell Formation from the Ribband Group, which consists of thinly bedded greywacke siltstones and slates (probably Lower Ordovician) (McConnell & Philcox 1994) surrounded by dark blue-grey slates and grey quartzites. Andesites and andesitic tuffs are present in the areas closest to the Leinster Granites, which were intruded during the Caledonian Orogeny (McConnell & Philcox 1994). The Blackrock to Newcastle shear zone (BNF) has an approximately east–west trend and cuts the northeast trending structures originated from the Killcullen and Ribband Groups. Geological studies of the Newcastle area were carried out with focus on fracture systems, trend of folds and mineralization (Jones et al. 1988; Phillips et al. 1988; Hitzman 1999). The main faults were formed during the Caledonian Orogeny and during emplacement of the Leinster Batholith (475–405 Ma, Chew & Strachan 2013; Fritschle et al. 2018). Later Variscan deformation caused tight folding, indicating northward rotation of the normal faults during Variscan deformation, and creating vertically dipping faults (Strogen & Somerville 1984; Hitzman 1999). In the western part of the study area is the Wheatfield Fault (WFF, Fig. 1), which is oriented east–west and splaying off from the BNF. The rotation of the fault zone resulted in possible formation of a dilation zones (Rothery & Phillips 1983; Phillips et al. 1988), important for mineralization processes within the structures as well as geothermal exploitation. The approximately east–west trending dextral shear zone of the BNF intersects the north–south extensional faults obliquely and is likely to create permeability at depth that could potentially host fluids (GT Energy 2009). A series of shallow and deep exploration boreholes were completed to map the BNF and basin's the Carboniferous succession (Strogen & Somerville 1984; GT Energy 2009). Recent geothermal exploration comprised four boreholes with multiple geophysical surveys to investigate structural conditions for large scale dilation zones along the BNF. The deepest borehole NGE1 (location marked on Fig. 1), which is approximately 1400 m deep, showed that the Carboniferous succession thickness exceeds this depth. The lowest lithological sequence contains quartz rich sandstones and interbedded shale units with high porosity and presence of fluids in the lower part of the Dublin Basin from two fractures (GT Energy 2009). The overall south to southeast dip of the structures is observed for most of the borehole by a televiewer technique, except for a northwest dip in deepest part of the well. The presence of folded structures in the upper part of the borehole is observed, as was expected, rather than major faults. The NGE2 borehole, which is situated less than 1 km to the north of the NGE1, shows similar lithologies but with thicker Calp Limestone Formation indicating that the succession is thicker towards the centre of the Dublin Basin. The upper limestones have low porosity with limited fault structures. These structures form thermal barrier resulting in an elevated geothermal gradient of 32.38 °C km-1 with temperature of 46.2 °C at the bottom of the NGE1 borehole, that is approximately depth 1400 m. 3 DATA ACQUISITION AND PROCESSING The MT method uses plane polarized, natural electromagnetic waves, which are vertically incident on the surface of an inhomogeneous Earth, to obtain information about subsurface distribution of electrical conductivity (and its inverse value, resistivity) within the Earth (Tikhonov 1950; Cagniard 1953; Chave & Jones 2012). Measured time-varying magnetic and induced electric (telluric) fields on the surface of the Earth are used to estimate the frequency domain MT impedance and vertical magnetic (tipper) response functions. The response transfer functions varying with period can be modelled/inverted to derive multidimensional models of subsurface resistivity structure. The resistivity models provide complementary information to other geophysical methods about geological units with conductivity contrast between structures from scales of 0.1–100 km. Determining the three-dimensional variation of electrical conductivity can play a significant role in developing an understanding dynamic, compositional and transport properties of the subsurface geology. Dry compact rocks, or areas that lack interconnected conductive paths, exhibit low conductivity. Conductive features, structures and anomalies indicate just the opposite: the structures of compact conductive rocks, domains of fluid accumulations, faults or increase in temperature (Chave & Jones 2012). At 40 sites geophysical soundings of both short period audio-magnetotelluric (AMT) and broadband period MT (BBMT) were acquired in the highly urbanized Dublin suburb in summers of 2011 and 2012 with MTU systems from Phoenix Geophysics. The spacing between stations was about 1 km with a remote reference station situated approximately 25 km from the known target area. The remote reference station was used in order to reduce the effect of local noise contamination in data through cross-correlation type methods through extracting coherent signal and removing incoherent noise. The MT data are highly disturbed by electromagnetic noise from surrounding infrastructure, Dublin's industry and nearby DC tram system (0.8 km from the closest site) and Irish Railway (0.6 km from northernmost site). The noise is caused by non-harmonic impulse-like artificial sources, such as engines during acceleration in the DC traction system or electric fences where high voltages impulses are applied every 2 min (Adam et al. 1986; Szarka 1987). The source plane wave and orthogonality assumption for the theoretical definition of the MT method are violated by close localization of sources (Boerner et al. 1993). The magnitude and lateral distribution of the noise varies with the resistivity of shallow structures (Lowes 2009). The subsurface structures in the investigation area are resistive rocks with a very thin (only up to tens of metres) conductive sedimentary layer on top (resistivity survey by GT Energy 2009). Best practice to improve the quality of MT data is to avoid areas with dominant noise by moving the site location (Oettinger et al. 2001). Unfortunately, our study area is restricted to the 10 × 15 km size with randomly distributed noise sources, and consequently the number of noise free spots is significantly reduced. The signal-to-noise ratio in our MT data varies with locality, and is generally very low at all measured sites. We overcame this noise problem with two methodical steps. In the first step, we derived spectrograms of the electromagnetic field to quantify the effects of unwanted man-made noise in each electromagnetic component. The spectrograms are heavily affected by distortion for periods greater than 1 s. The most quiet periods of each day, visible directly from spectrograms and also the pre-processed transfer functions, are night periods from 0:30 am to 3:30 am (see Fig. S1). During this time window, most industry and train systems in the surrounding region are not operational. For the BBMT sites, we were able to select the lowest disturbed time window from three nights (Fig. S2) and select the best information. The AMT short-period measurements were carried out for just one night with no possibility to select the best night window. Therefore, the transfer functions for periods below 1 s were estimated from fewer amounts of quality data and for long periods we have to use sufficient but more distorted data. The selection of a proper processing method and its application to heavily distorted data is crucial for our study. There are several studies comparing various estimations of statistical MT parameters (Weckmann et al. 2005) and different processing methods (Jones et al. 1989). In the second step we used three different processing approaches to compare the stability of MT impedance and tipper transfer function estimates. The multiple noise sources with random lateral distribution and character did not allow detailed data analyses and filtering, so we decided to determine the noise effects by considering the differences in applying different processing methods. The first tool was the robust processing package EMTF from Egbert (1997), which is based on multivariate statistical methods. The second tool was the code by Smirnov (2003) based on one-step reduced M-estimator with Huber weights. All processing tasks were performed on 6 different time-series in Phoenix instrument format (AMT with sampling frequencies 24 kHz, 2.4 kHz and 150 Hz; MT with sampling frequencies 2.4 kHz, 150 Hz and 15 Hz) and estimated transfer functions were merged to one final station file. The third tool was the processing code from Phoenix Geophysics provided with their MT systems, and that is based on Jones's Least Trimmed Squares (LTS) robust code (Jones & Jodicke 1984). Generally, signal-to-noise ratios are at a minimum for the dead-band zones frequency ranges (Fig. S3). We were not able to estimate decent quality transfer functions for AMT dead-band zone (between 800 and 3 kHz) at any station (Garcia & Jones 2002). For final modelling, the results from Phoenix Geophysics processing code were used, as they visually were the most superior (Fig. 2 with data example of a good and bad site). The maximum periods for most of the sites (15 sites) used in modelling reached approximately 0.1–1 s. Only for four sites were there greater maximum periods, of about 100 s (see misfit plots in Figs S5–S8). However, the penetration depth of these short periods at most of the sites is much greater than depths of interest of the project, due to very resistive structures in the study area (Fig. 3). The shallowest penetration depth is in the most northern site of profile NC12 (1 km) and deepest site (∼100 km) at southern area, where only five sites are below maximal depth of our models, 5 km. The skin depth for each site is also shown as a dashed line in 2-D models in the following sections of the manuscript. Figure 2. Open in new tabDownload slide Example of good sounding curves for site NCW12 and bad sounding curves for BUS13 obtained by Phoenix commercial processing tool. The top row shows apparent resistivity with period and the bottom row shows phase shifts between electric and magnetic field. Shorter periods indicate shallower information about resistivity and longer periods refer to deeper information (approximately from 0.05 to 70 km). Data points crossed by red X were not taken to the analysis. Figure 2. Open in new tabDownload slide Example of good sounding curves for site NCW12 and bad sounding curves for BUS13 obtained by Phoenix commercial processing tool. The top row shows apparent resistivity with period and the bottom row shows phase shifts between electric and magnetic field. Shorter periods indicate shallower information about resistivity and longer periods refer to deeper information (approximately from 0.05 to 70 km). Data points crossed by red X were not taken to the analysis. Figure 3. Open in new tabDownload slide Plot of individual mean geoelectrical strike angle estimates for each site and both profiles NC11 and NC12. The colour of the vectors indicates the maximum skin depth for the MT data. The colour of the squares indicates the ratio between the maximum Niblett–Bostic depth from the Berdichevsky average of the impedance tensor (Berdichevsky & Dmitriev 1976) and shown 5 km depth of profile, that is 20 per cent means depth 1 km. Figure 3. Open in new tabDownload slide Plot of individual mean geoelectrical strike angle estimates for each site and both profiles NC11 and NC12. The colour of the vectors indicates the maximum skin depth for the MT data. The colour of the squares indicates the ratio between the maximum Niblett–Bostic depth from the Berdichevsky average of the impedance tensor (Berdichevsky & Dmitriev 1976) and shown 5 km depth of profile, that is 20 per cent means depth 1 km. 4 METHOD TO TEST THE 2-D MODELLING RELIABILITY—EXAMPLE FOR SYNTHETIC DATA There are numerous studies focused on quantification of data noise, model mesh, inversion parameters and a priori assumptions of the influence on the final geoelectrical inversion model uncertainty (Muñoz & Rath 2006; Kalscheuer et al. 2010; Menke 2012; Schnaidt & Heinson 2015). Methods that use a probabilistic approach for data inversion have posterior probability distributions for the model parameters, and they directly indicate uncertainties in the model. In the code used in this study, the non-linear inverse problem is addressed by an iterative linearized least-square solution, where the calculated sensitivity matrix relates changes in the model parameters with changes in the response. The underestimated new variances of impedances from decomposition are replaced by larger fixed error floors. Therefore, our noisy data are hidden amongst the good quality data and we need to estimate model parameter bias following from the unwanted inclusion of noisy data that can act as leverage points in the inversion process. We expect that poor data subgroups are uncorrelated and produce distinct distorted models. Therefore, we implemented a simple stability approach based on a statistical jackknife method, where we resampled the inversion input data set and analysed the variability of the resulting models. Statistical testing is not focused on the problem of sensitivity of the model to input MT data, but on the problem of information inconsistency in data due to noise among stations and the distortion effects on our final geoelectrical models. Note that for the most accurate data inversions, for example the synthetic data inversions, the variance estimation by our method should be close to zero. To demonstrate the advantages of our reliability estimation method for inconsistent input data sets, we generated synthetic data with noise for a 2-D model composed of layers of conductive and resistive structures (three layers with regularly distributed bodies), as presented in Fig. 4. Some parts of synthetic data were replaced by data sets obtained from reduced models (RMs), where one or two of the layers were removed. RMs generated responses with missing information about some of the layers, and we studied their effect on model variance in combination with the full data set. Two versions of input inversion data sets v1 and v2 were tested, where every second (v1) or fourth (v2) site was replaced by simpler data set with missing information about the deepest layer (reduced model—RM1), the shallowest layer (RM2) and two deepest layers (RM3). With these three different synthetic models we demonstrate effects with structures in different depths. The inversion parameter settings, that is error floors and regularization parameters, were identical for all calculations. The results show the highest variance and overall increase of variance within the model for the most reduced data set v1 and RM3. The RM2 (with missing shallowest layer) exhibits highest variance in much deeper structures and not for shallow structure. The elevated variance is situated under a conductor in the lowest constrained area. This feature is due to the static shift effects correction in the inversion algorithm being disabled. The data set with a high amount of missing data led to rough structures with blobby variance. The higher regularization parameter smooths out the variance in conductivity within the model, that is the amplitude of the lateral changes within the model and also between particular separated inversions are restricted. Therefore the lower nominal variance is lower with higher smoothing. It is interesting that the type v1 (50 per cent of input data set is replaced with synthetic data from reduced models) and v2 (25 per cent replaced data) noise variance differ for shallow and deeper depths (Fig. 4). Shallower variances are slightly higher for v1 type and v2 type reaches higher values in bottom part of models. So with less inconsistence in the input data the final model is more stable for shallow structures, but the position of a bad data site can affect resolution of deeper structures. Figure 4. Open in new tabDownload slide Conceptual scheme of the reliability method testing (Jackknife method) for different synthetic noisy data sets v1 [each 2nd site are with responses from reduced model (RM)] and v2 (replaced each 4th site), different smoothing parameters in objective function for data set v1 (global smoothing parameter), and different synthetic input data sets, where certain parts of original synthetic model (top left) are missing (RM1, RM2 and RM3). More description and discussion in Section 4. Figure 4. Open in new tabDownload slide Conceptual scheme of the reliability method testing (Jackknife method) for different synthetic noisy data sets v1 [each 2nd site are with responses from reduced model (RM)] and v2 (replaced each 4th site), different smoothing parameters in objective function for data set v1 (global smoothing parameter), and different synthetic input data sets, where certain parts of original synthetic model (top left) are missing (RM1, RM2 and RM3). More description and discussion in Section 4. 5 DATA AND STRIKE ANALYSIS The expected dominant east–west orientation of geological structures in the Newcastle area, represented by the BNF, led to the MT survey being designed as a 2-D problem with sites distributed mostly along profiles oriented in a NNE–SSW direction perpendicular to the BNF (Fig. 1). The pre-modelling step for 2-D studies is dimensionality analysis of impedance transfer functions. The input impedance data were heavily edited with the help of transfer functions from different processing codes. We used inconsistency between transfer function estimates from the different codes as an indicator for poor data. Only 23 stations out of 40 had sufficient good quality data to be used in the modelling, that is all four components could be used in dimensionality analysis. To prepare the input data for 2-D inversion, assumptions must be made about the orientation of geoelectrical 2-D strike of the two MT modes for each station, the mode where the electric field is along strike (TE mode), and the mode where the electric field is across strike (TM mode). The distortion decomposition of impedance tensor to strike direction allows removing some level of static shift distortion and provides data that are appropriate for 2-D modelling. Several approaches are used in regional surveys to reduce dimensional complications and distortion in the data (see a comparison in Jones 2012b). Carefully selected and edited processing results were analysed using the McNeice & Jones (2001) multisite and multifrequency extension to the Groom–Bailey MT tensor decomposition technique (Groom & Bailey 1989). In this technique, a global minimum is sought to determine the most appropriate regional 2-D strike direction and telluric distortion parameters for the range of frequencies (or approximate Niblett–Bostick (Jones 1983) depth ranges) and a set of sites. The penetration depth varies significantly among the sites at the same period, so we decided to use more relevant a depth‐related method. The analyses were performed for a depth range from the surface to 10 km and, as the phases are unaffected by static shifts, the maximum phase instead of the minimum resistivity was use to define the most conductive direction. The depth range was calculated from arithmetic (Berdichevsky) average of MT responses (Berdichevsky & Dmitriev 1976). We performed analyses independently for each site (Fig. 3) and in the final step estimated one strike direction for each profile. The statistics performed in this section were undertaken using a circular von Mises distribution with a 90° repetition frequency, following Mardia (1972). In Fig. 3 we also characterize maximum depth penetration for each site (square colour) and robustness of depth range estimation (vector colour). The dimensionality and strike direction were also checked by phase tensor and skew angle beta parameter representation (Fig. S4) of the primary impedance data (Caldwell et al. 2004; Booker 2014). The phase tensor is graphically depicted as an ellipse, where a circle represents 1-D data, and for the 2-D case the major or minor ellipse axis is aligned with geoelectrical strike. The magnitude beta significantly greater than 3° indicates 3-D data. From Fig. S4 we can estimate predominantly 2-D and 3-D character of data with strikes approximately in an east–west strike and they smoothly vary with periods. The main conclusion from these analyses is the lower penetration to the north and 3-D character near offset in the BNF in the middle of NC12 (contact point with Wheatfield Fault). The resulting MT transverse electric (TE, along strike e-field) and transverse-magnetic (TM, across strike e-field) data mode are off-diagonal components of impedance tensor decomposed to a newly defined regional azimuth angle. As a final step of decomposition, quality checks of the decomposed TM and TE data modes were performed by testing the consistency of the apparent resistivities and phases with theoretical predictions from the best fitting Rho+ model (Parker & Booker 1996), as demonstrated by, for example Spratt et al. (2005), followed by last round of editing/deleting to eliminate poor data. Due to the results and different level of the distortion in measurements and the lower effect of possible 3-D character of data, we divided the sites into two profiles grouped by year of collection (year 2011–9 sites with names NEW*, profile named NC11; year 2012–14 sites with names NCW* and BUS*, profile named NC12). The sites collected in 2011 are better quality time-series, however they are closer to Dublin city. The main reason for better 2011 data is that the DC tram system (LUAS) started to operate at a nearby station just after our 2011 measurements. The traffic on this line is more intense (each couple of minutes) compare to Irish Railway in the north. Regional geoelectrical strike azimuths of N125°E for line NC11 and N110°E for profile NC12 were determined. To obtain general 2-D structures for whole area the global NC 2-D models with all quality sites was also derived and strike angle N116°E was adopted. The estimated geoelectrical strike directions are very close to our expected regional structures represented by the surface trace of the BNF. The progressive change of strike from the east to west does not follow the expected rotation of the fault to a more northwest–southeast orientation. It is probably caused by the presence of Wheatfield Fault below NC12 with an east–west orientation, which slightly affects the strike direction for NC12. 6 2-D MODELLING The derived electrical resistivity models of the regional 2-D structure of the BNF and surrounding area facilitated undertaking comparative studies with seismic reflection profiles in the same area. The distribution of stations in the area defines NNE–SSW oriented models crossing the BNF, but with insufficient lateral coverage of the area that is necessary for precise 3-D modelling. Compared to 3-D inversion, 2-D inversion can be run effectively on a single profile and is far faster (by orders of magnitude) and requires far less computing power (laptop compared to cluster) than 3-D MT inversion; it can be quickly rerun to test the reliability and resolution hypotheses. The principal advantage of MT 2-D inversion is that it allows for controlling dimensional distortion that creates artefacts in the MT data models, and it could produce less distorted models, closer to our expected structures. The strike direction decomposed data from each profile were inverted using an 2-D inversion code (Baba et al. 2006), based on the nonlinear conjugate gradient algorithm of Rodi & Mackie (2001). Isotropic conductivity modelling of the impedance and tipper data was performed. The different sets of regularization parameters for model smoothness were tested to find the optimal balance between model roughness and final misfit. The starting model was same for all models and it was 100 Ωm halfspace with a fine mesh (112 horizontal and 83 vertical cells for NC11, 95 × 66 for NC12, thickness of first layer is 13 m, horizontal cells size are tens of meters and depends on the site positions on the profile). We also allowed the code to invert for static shift effects (with default settings). After decomposition, the new error bars of the data (apparent resistivity and phases) were derived based on statistical realizations, and they are usually too small. Therefore, the error floors of inverted data-types play an important role in the final inversion process and models. The applied error floors for 2-D principal modes represent new data error estimation. The phase error floors are set to 5 per cent (i.e. 2.8°), and the apparent resistivity error floors set to 20 per cent. The final presented set of 2-D models is based just on the inversion of impedance data, that is the tipper data were not included. In the two most northern sites (NEW01 and 02) in NC_All modelling, the longer period data were removed to improve data fit. The geomagnetic transfer functions (GTF) are much more affected by artificial noise, so these data were excluded in the final inversion runs. The normalized RMS misfits (the reduced chi-square value) between the observed and modelled data are 2.32 for NC_All, 1.9 for NC11 and 2.6 for NC_12 (should be 1 for a fit by 68 per cent of the data to within errors, or 2 for 95 per cent to within double the errors). Detailed information about the fit is presented in the supplementary part of this paper (Figs S5–S7). Due to the distorted input data for 2-D inversion and higher resulting RMS fit, we decided to perform a stability technique examination of the models. We show the 2-D reliability plots based on simple formulae for model conductance (calculated for each cell) and resistivity variance. The averaged 2-D models of the same structures from several independent inversions for different input data were constructed. In Fig. 5 we present the example for the shortest simple profile NC_12 with nine sites, where we dropped data from one station at a time and inverted the remaining eight stations. The same process was repeated for each of the sites and we obtained nine different inversion models based on 8 sites. The presented median model calculated from subinversions cannot be used as the final model, because it does not represent the best fitting model. However, the resulting variability between the models from each subinversion mainly comes from incompatible information about the subsurface among the nine inverted sites, and it can be used as the estimation of incoherent subsurface information. The presented version is the logarithmic standard deviation of the 2-D model resistivity in each cell from all calculated submodels. A standard deviation is used as the resistivity error input for later estimations of permeability or integrated modelling with other geophysical data. Reliability models indicate which conductive bodies tend to influence or modify the total conductance. Figure 5. Open in new tabDownload slide The example of reliability test of NC11 line, where 2-D model is based on 9 independent selection of inverted points (left-hand panel) with estimation of uncertainty of averaged model (right-hand panel). Figure 5. Open in new tabDownload slide The example of reliability test of NC11 line, where 2-D model is based on 9 independent selection of inverted points (left-hand panel) with estimation of uncertainty of averaged model (right-hand panel). All 2-D models (Fig. 6) are reliable in the areas with significant conductive structures, which is important for identification of water-bearing layers. The easternmost model NC11 has the highest values of variance exactly in the northern area, where the parallel profile NC12 shows conductive structures (conductive structure C1). The variance values for both profiles are the highest up to a depth of 2 km, while the combined NC_All profile exhibits an elevated variance in deeper structures and in aforementioned shallow structures in the northern area. Low variance in the shallower middle part of the model is due to the very dense layout of stations projected onto the profile line. Some of the stations are more than 2 km away from profile line, which is why the variance is underestimated. Figure 6. Open in new tabDownload slide 2-D geoelectrical models (NC11, NC12 and NC_All) with reliability standard deviation estimation, where red means less reliable part of the model. The labels C1–C4 indicates positions of main conductive features discussed in the Section 6. Figure 6. Open in new tabDownload slide 2-D geoelectrical models (NC11, NC12 and NC_All) with reliability standard deviation estimation, where red means less reliable part of the model. The labels C1–C4 indicates positions of main conductive features discussed in the Section 6. From the models we see that the geometries of the units close and below the BNF are complex and complicated, probably because the area has a far greater number of faulted and fractured horizons than known and mapped at the surface. For structures further into the Dublin Basin (more to the north from the BNF), we observe a simplified layered character with a low angle of apparent dip towards the north. Structures south of the BNF are more compact and resistive with a lack of deeper conductive layers (C2). The models exhibit relatively resistive structures in the survey area, with more conductive features (about 10 Ωm) to the north of the BNF (C1). This could indicate that the structures in our models could have underestimated conductivity, what should count in the final interpretation presented below. Generally, the models from the NC11 and NC12 profiles display very different conductivity distributions within the structures despite their closeness. However, they exhibit similar conductive structures near the BNF (C3, C4), that is the conductor C3 north to the BNF up to depth 1 km and the conductor C4 at depths from 1 km (1.3 km for NC12) to 2.2 km. The difference in conductivity for the rest of the structures is caused by 3-D character of structures, mixing of perpendicular faults and structures, and different strike direction. The general NC_All model does not resolve the discrepancies between eastern and western 2-D model and exhibits new conductive structures below 2 km with higher reliability variance. Therefore, we decided to pursue 3-D modelling of the data even though the sites were distributed for 2-D modelling in two profiles. 7 3-D MODELLING 2-D modelling alone is not recovering the higher complexity of geological structural geometries in the area, and 3-D modelling is required to reveal correctly all structures and their properties. We see significant changes of structures in the east–west direction from the two parallel north–south 2-D models in the previous section. In 3-D MT inversion more parameters in input data sets and higher dimensionality of the model bring more degrees of freedom than 2-D inversions. The input data, where the diagonal components of impedance tensor are included, require more careful editing to avoid unrealistic results in the final inversion model. Unfortunately, even the latest forward and inverse 3-D codes are not able to reproduce the full real electromagnetic field components collected in the field with all distortion effects accounted for (Jones 2011), with the exception of the x3di code (Avdeev & Avdeeva 2009; Avdeeva et al. 2015) which was not used in this study. The original MT transfer functions did not undergo any distortion removal processing and were directly used as input into the 3-D inversion code. The only modifications performed were manual removal of spurious data and selection of data for four periods per decade to ascertain evenly spaced period numbers. The best quality 23 MT sites in the Newcastle area were inverted using two different 3-D inversion codes. The inversions were performed by one code based on a data space inversion algorithm (Siripunvaraporn et al. 2005) and the second modelling code, which uses a model space inversion algorithm (Egbert & Kelbert 2012; Kelbert et al. 2014). The data space 3-D inversion model that we present and discuss here was obtained by the newest version of the 3-D MT inversion program WSINV3DMT, based on the data-space variant of the OCCAM scheme (DeGroot-Hedlin & Constable 1990), which includes modifications to allow inversion of the geomagnetic transfer functions, or ‘tippers’. The inversion models from the model space ModEM code exhibited little sensitivity to deeper structures for our data and the final RMS misfits were higher. The ModEM code also allows inversion with HMTF or phase tensor data type, but their inclusion in the inverted data did not improve the final model. We selected an inversion model from WSINV3DMT code as our final model that was based on the inversion of full impedance tensor. The parallelized code version was used to facilitate its implementation on computer clusters and improve the time of computation. The mesh size depended on the selected sites included in the inversion and was 110 × 90 × 54 (110 in the N-direction, 90 in the E-direction and 54 vertically). The horizontal size was 250 × 250 m for central part of the model, with a thickness of 20 m for the first layer, and the total depth of the model with numerical padding zone is 605 km. Initialization parameters input into WSINV3DMT included the error floor for the impedances, defined as a percentage of square root of multiplied off-diagonal xy and yx impedance components; an impedance error floor 5 per cent was used for the final model. However, two other sets of error floors were used to test different data weighting in the inversion between diagonal and off-diagonal components (percentage of square root of xx and yy impedances for diagonal components) or all four components alone, but the resulting models had worse fit without any improvement in information about the conductive structures. The inversion of impedance is also sensitive to the global resistivity of the starting model (usually a homogeneous half space), and a poor choice of starting model can result in being trapped in local minima of the misfit penalty function. We used a halfspace start model with a resistivity 100 Ωm. To avoid an unrealistic, rough model, different smoothing parameters for the inversion algorithm were tested. Topography was not included into models; it is not a significant issue as we are concerned with structures on scales of kilometres and change in topography in the area is of order a few tens of metres. In Fig. 7 we present horizontal slices through our final 3-D inversion model, derived after several rounds of inversions, where the starting model is the best fitting model from the previous inversion run. The RMS misfit for the final model is 2.7. The 3-D RMS misfit is higher than 2-D modelling misfits, which are presented above. However, the error bars applied for 3-D RMS misfit calculation are smaller than in the 2-D case, and also the misfit contribution of less stable diagonal (xx and yy) components of impedances is included in this aggregate. Due to the high memory requirement of WS3DMT, a coarser mesh was necessary and therefore some detail in 3-D was sacrificed compared to the 2-D modelling. The higher or lower smoothing parameters were tested without significant simplification of the structures in the resulting final model or with significant improvement of the RMS value. The horizontal distribution of misfit, summed over all frequencies per site, is presented in Fig. 8 with the full impedance matrix and each component fit. The fit corresponding to derived impedance modules and phases of all four components are also presented in the supplementary information (Fig. S8). Figure 7. Open in new tabDownload slide Horizontal sections of final 3-D model for different depths with shown surfaces traces of major faults and dip symbols in the modelled area (GSI online maps). Figure 7. Open in new tabDownload slide Horizontal sections of final 3-D model for different depths with shown surfaces traces of major faults and dip symbols in the modelled area (GSI online maps). Figure 8. Open in new tabDownload slide Horizontal distribution of misfit between original and predicted data for presented 3-D model. Figure 8. Open in new tabDownload slide Horizontal distribution of misfit between original and predicted data for presented 3-D model. The final model presented in Fig. 7 exhibits much higher conductive structures in comparison to the 2-D models, with similarly resistive hosting rocks. The lateral boundaries of shallow conductive structures, down to a depth of 1 km, have a north–south elongation correlated with the surface traces of faults or structural boundaries according to the GSI maps of bedrock stratigraphic and structural lines, which are perpendicular to the BNF. Deeper structures become more oriented to regional geoelectric strike (approximately N110°E) similar to the 2-D regional strike estimated in the previous section, and less conductive, as observed in the 2-D modelling. The only exception is a strong conductor in the southern edge of the investigated area, but this conductor is only based on responses from one or two sites and is uncorrelated with investigated dilation zones around the shear zone along the BNF, so is suspect. 8 MT MODELS DISCUSSION The main rationale behind presenting both the 2-D and 3-D models is to combine the simplification of real geoelectrical structures produced by 2-D modelling with the identification of possible 3-D bodies or structures distributed along the designed profiles. The area exhibits high resistive background structures, and AMT and MT data from the highest frequencies down to 1 Hz are adequate to obtain sufficient information to depths of more than 5 km. The conductive bodies could be associated with local dilation zones with enhanced secondary porosity and permeability, which are some of the desired targets for geothermal exploitation. Based on geological mapping in the area, 2-D structures along the Newcastle Blackrock shear zone were expected, but both models indicate a much more complicated situation, particularly for shallower structures. The MT inversion study suffers from some unreliable primary data, and both types of modelling can help to extract and strengthen robust, consistent information about conductive structures. We have completely removed one quarter of all stations in the area and almost half of all periods (mostly longer periods) at the stations we retained. We did this in order to be sure that we are not influenced by spurious MT data in the inversion modelling. There are well known problems of galvanic and induction effects of 3-D structures that cannot be removed by classic decomposition techniques (Ledo 2005; Jones 2012a) and the distribution of them has an effect on the two modes of 2-D MT data. Fortunately, the geology of the area is dominated by high resistivity rocks and any conductive zone can be inferred, with high probability, to be faults or fractured zones. The comparison of all available geophysical information is presented in Fig. 9. The 3-D inversion model shows a reasonably consistent match with the 2-D model for depths down to 5 km. The dissimilarity comes from different fundamental geometries (2-D versus 3-D) of the models. There are several studies focused on comparison of 2-D and 3-D approachs in modelling of MT data (Berdichevsky et al. 1998; Park & Mackie 2000; Ledo 2005). In our case, the 2-D models average all local conductors against very resistive background along the strike direction, which leads to underestimated conductivities within the model. The small change in strike direction, which is determined by a statistical average over stations and MT data depth penetration on profile, can cause a rapid change in the along-strike average of conductive structures projected to the profile line. On the other hand, the 3-D geoelectrical model can easily localize non-regional anomalies and the model tends to fragment compact conductive structures, particularly for surveys with site coverage or noisy data insufficient for 3-D modelling. Figure 9. Open in new tabDownload slide Left-hand panel: vertical section from 3-D model versus general 2-D model NC12 with included NGE1 borehole information and with indicated zones of no penetration (maximal skin depth). Isolines from the seismic velocity models and thin dashed line as representation of the BNF are based on Licciardi & Agostinetti (2017). The conductive features were label as C* following Fig. 6. Right-hand panel: comparison of borehole logs with conductivity information from 2-D and 3-D geoelectric models. Figure 9. Open in new tabDownload slide Left-hand panel: vertical section from 3-D model versus general 2-D model NC12 with included NGE1 borehole information and with indicated zones of no penetration (maximal skin depth). Isolines from the seismic velocity models and thin dashed line as representation of the BNF are based on Licciardi & Agostinetti (2017). The conductive features were label as C* following Fig. 6. Right-hand panel: comparison of borehole logs with conductivity information from 2-D and 3-D geoelectric models. The typical example of behaviour described above can be identified in resistivity logs from the NGE1 borehole (Fig. 9, right-hand panel). There are no MT stations with good quality data within a 1 km radius of the borehole, therefore, the resolution of the MT models for comparison with the borehole vertical resistivity is poor, especially at shallow depths. The 2-D model, where lateral resolution is not decreasing in the along-strike direction and is determined only by distance of sites along profile, has better correlation with the depth of conductive layers in the borehole, but underestimates conductivity and exhibits higher resistivities of underlying structures. On the other hand, the resistivity distribution from the 3-D model bears a general resemblance in absolute resistivity values, but underestimates the depth of the bottom boundary of conductor. Based on the assumption that we omit well-known poor depth resolution of the bottom boundary of a conductive anomaly, there is a lateral resolution limitation from the locations of the sites of 3-D models. To visualize the model ‘no-resolution’ zones, we masked certain parts of models in horizontal sections (from Fig. 10). Figure 10. Open in new tabDownload slide Horizontal sections at depth 1 km and 3 km with indicated main features in the model and their possible interpretation. The purple lines in the 1 km section highlight strong lateral boundaries that are correlated with surface traces of faults. The 3 km section shows deeper conductive structures that are oriented along the BNF (purple semi-transparent areas). The high conductive structure in the south is interpreted as unreliable (black arrow). Figure 10. Open in new tabDownload slide Horizontal sections at depth 1 km and 3 km with indicated main features in the model and their possible interpretation. The purple lines in the 1 km section highlight strong lateral boundaries that are correlated with surface traces of faults. The 3 km section shows deeper conductive structures that are oriented along the BNF (purple semi-transparent areas). The high conductive structure in the south is interpreted as unreliable (black arrow). The conductivity of fluids and their interconnectivity, that is porosity and permeability of the host rocks, which determines conductivity anomalies are mapped by direct measurement in NGE1 borehole. The fluid conductivity of artesian flows in the borehole is increasing from shallow groundwater, through slightly saline water in middle depths, to conductive (∼1 S m–1) saline brines presented in the lowest parts (below 1000 m). The geology and borehole information suggest the presence of a fractured resistive matrix, that is secondary porosity, with conductive fluids at drilled depths. The less conductive structures from the MT models down to the borehole terminal depth can result from the vertical averaging of several conductive thin layers that cannot be resolved by MT inversion due to the low vertical resolution of the method at these depths. The investigated area below Newcastle can be divided into two different depth zones. The first zone down to 1–2 km is dominated by NNE–SSW oriented conductors connected with shallow faults or structural fold boundaries probably filled with saline waters. The conductors cross the surface trace of the BNF in an approximately perpendicular direction. The shallow conductivity layer (down to the hundred metres depth) north of the BNF is interpreted as the limestone of the Upper Calp Formation. The second zone can be identified from depths of 2 km down to maximum 4 km, where structures are oriented along the BNF and have lower conductivity (Fig. 10). The 3-D model exhibits a change in regional strike, shown by NC11 and NC12, from 126° to 110° from east to west. The part of the NC_All model where 2-D conductive structures below the BNF are smeared to greater depth are not taken into consideration, because these areas have also high standard deviation and are thus less reliable. We have relatively good correlation between the 2-D and 3-D models for this depth range (Fig. 9). Structures in the centre of our investigated area at about 4 km depth and deeper are more resistive, and we cannot identify any significant conductors associated with the BNF; either a water-bearing layer or other conductive faults. The strong conductor in the south of area probably originates from the distortion at one of the southernmost sites, and is considered unreliable. The final visualization of suggested interpretation is shown schematically in Fig. 11. Figure 11. Open in new tabDownload slide Final visualization of structures orientation in the Newcastle area for proposed two depth zones. The red dipping units in shallower zones represents conductors correlated with folding structures. And the green structures in deeper zone show dilation zones with higher conductivity oriented along the BNF. Figure 11. Open in new tabDownload slide Final visualization of structures orientation in the Newcastle area for proposed two depth zones. The red dipping units in shallower zones represents conductors correlated with folding structures. And the green structures in deeper zone show dilation zones with higher conductivity oriented along the BNF. 9 DILATION ZONES MODELLING One of the goals of the survey was to provide the geoelectrical information on the base of the Carboniferous sequence at a depth of 4 km. The focus of the IRETHERM collaborator, GT Energy, is the permeable dilation zones and associated fracture networks created by shear reactivation governing the deformation on existing older faults. Since the primary porosity observed during the drilling (GT Energy 2009) is small and the rocks are tightly cemented, permeability is established by secondary porosity of fractures, which we image through the volumetric electrical connectivity of fluid filled fractures that is mapped in this study. The most recent geophysical studies of dilation zones around the BNF were carried out by a teleseismic receiver function study (Licciardi & Agostinetti 2017). The observed anisotropy and low anomaly in S-wave velocity is tentatively interpreted by those authors as fluid filled aligned cracks at depth 2.3 km. The results are promising, due to the good correlation with the zone of high secondary porosity encountered at the bottom of the borehole NGE1. However, the resolution of the seismic data is lower than that of our MT data, and the robustness of the anisotropic properties is influenced by the low number of stations and their wide spacing. The comparison of velocity models and cross-sections from MT models is presented in Fig. 9, where the BNF is identified as the lateral seismic velocity and conductivity contrast in proximity to the basin margin. The comparison models shows that the 3-D MT model has good correlation with seismic model. However, the MT 2-D model is only acceptable for shallow depths and the possible dilation zone is shifted to the south. From all models we can estimate that the deepest, and largest, dilation zone is situated near the BNF in the south–southwest direction from NGE1 borehole and east from the junction area of the Wheatfield fault and the BNF. To obtain a better characterization of thermal transport properties of investigated area we used porosity and resistivity data from the borehole NGE1 to estimate the relationship between porosity/permeability and electrical conductivity. The formulae are based on a generalized Archie's law for multiple phases (Archie 1950; Glover 2010; Campanya et al. 2015). We used shale-corrected estimation of porosity based on the sonic velocity values from NGE1 borehole logs. Two geological lithologies are presented as an example with different resulting cementation exponents: estimations for unit 3 represents Conglomeratic Limestones in the shallow depths from 0.450 to 0.650 km and unit five characterizes Quartz Rich Carbonate Grainstones in depths of 0.85–1.01 km. These layers were carefully selected based on the borehole imagery logs as lithologies without the visible presence of fractures. We used the KCGL approach (Glover et al. 2006, Campanya et al. 2015), based on porosity and cementation exponents, to derive permeability estimations (Fig. 12). Figure 12. Open in new tabDownload slide Porosity modelling: (a) Borehole lithology with total dilation space of fractures; (b) Cross-plots of bulk conductivity, fluid conductivity and porosity from borehole NGE1 measurements; (c) Archie's law estimations plots; (d) Probability distribution of cementation exponent and (e) Probability distribution of permeability based on KCGL approach Figure 12. Open in new tabDownload slide Porosity modelling: (a) Borehole lithology with total dilation space of fractures; (b) Cross-plots of bulk conductivity, fluid conductivity and porosity from borehole NGE1 measurements; (c) Archie's law estimations plots; (d) Probability distribution of cementation exponent and (e) Probability distribution of permeability based on KCGL approach The results of cementation exponent modelling for different lithologies within the borehole vary with depth as they are affected by fissures and fractures (secondary porosity) in the studied rocks. Therefore, the reliability of permeability estimation based on our generalized Archie's law and assumptions of cementation exponent is not high. The observed intrinsic rock porosity within the structures from borehole information is low, and the dominant sources of fluid flow are large faults within the lithology. Both the fluid flow and electrical conduction rely strongly on the distribution of pore space as well as the total volume (Brown 1989), therefore relating resistivity and permeability directly, as presented in, for example Kirkby et al. (2016), could be more relevant to evaluate the fluid flow properties of fractures and faults. We focused our permeability analysis on the bottommost part of the borehole geology (GT Energy 2009) with highest fracture density (in depths of 1.015–1.075 km) and highest fluids inflow (in depths of 1.325–1.345 km). Following synthetic studies by Kirkby et al. (2016), we set the primary input for the ratio m of the rock matrix (∼10 000 Ωm) to fluid resistivity (∼1 Ωm) as approximately 10 000, and ratio of the matrix to fracture resistivity (∼200 Ωm) M ∼ 50. The minimum matrix permeability can be set according to the estimation shown in Fig. 12 for non-fractured zones, which is about 0.1 mD (10−16 m2). The number of fractures and faults registered by borehole imagery logs (GT Energy 2009) is restricted by resolution of the method and this number is 34 events (from 364 registered for whole borehole) over a 60-m-thick layer (1015–1075 m). We used synthetic modelling results for fractures with one fault in the strata analysis. If we visualize approximate m and M values following Fig. 6 from Kirkby et al. (2016), the permeability interval based on curves for single fracture estimations in Fig. 13(a) is placed in a permeability threshold zone, where the permeability changes by two magnitudes for very small changes of M value. A rock with system of faults is more realistic scenario than a rock with intrinsic porosity/permeability. The faulted rock is simulated with variable width (perpendicular to the fault plane), containing a single fault, and changing the width of the modelled rock is equivalent to changing the spacing between faults (Fig. 13b based on fig. 10 from Kirkby et al. 2016). For the second case the permeability interval is increased by more than an order of magnitude to many orders of magnitude, from ∼10−14 × 10−12 m2 to ∼10−13–10−10 m2 (100 mD–100 D). According to Bear (2018), the range of estimated permeability in this study can be classified from pervious to semi-pervious; values herein are above the permeability of most geothermal reservoirs (1–100 mD, Björnsson & Bodvarsson 1990). The resulting permeability estimated by the fault/fissure/fractures approach is higher than the permeability estimated for layers without fractures by Archie's law approach. Note that Kirkby et al.’s (2016) approach is based on synthetic modelling at 10 × 10 cm scale, which can result in underestimation for large fault structures with unknown aperture. Based on the observed high fluid inflow and existence transport properties of the rocks, higher permeability is expected. Figure 13. Open in new tabDownload slide (a) The selected ranges of single fracture permeabilities (green lines) based on our estimated conditions in Newcastle area (red lines) and (b) Estimation for different spacing between faults, which are estimated based on Kirkby et al. (2016) (modified Figs 6 and 10). Figure 13. Open in new tabDownload slide (a) The selected ranges of single fracture permeabilities (green lines) based on our estimated conditions in Newcastle area (red lines) and (b) Estimation for different spacing between faults, which are estimated based on Kirkby et al. (2016) (modified Figs 6 and 10). 10 CONCLUSIONS The combined 2-D and 3-D MT modelling has imaged a dextral shear zone with the BNF as a highly fractured and partly conductive zone. The MT data are strongly affected by artificial noise due to nearby industry and railway systems, nevertheless, because of resistive structures, the penetration depth of the good quality data extends down to 10 km, and we discern that multidimensional MT modelling yields reliable geological interpretations to 5 km depth. We have carried out stability and reliability studies of the prepared geoelectrical models with comparison to borehole information in the study area. The investigated area below Newcastle can be divided to two different layered zones. The first zone, up to 1–2 km deep, is dominated by NE–SW oriented conductors that can be spatially connected with shallow faults probably filled with saline waters. These conductors cross the surface trace of the BNF. The second zone can be identified from depths of 2–4 km, where structures are along the BNF and the observed conductivity is lower. The deeper conductive layers are interpreted as water- or geothermal-fluid-bearing rocks, and the porosity and permeability estimations from the lithological borehole logs indicate the geothermal potential of the bedrock to deliver warm geothermal waters to the surface. The BNF is visible in the models as a conductive feature in the second zone and is interpreted to be a highly fractured fault system infilled by saline waters. Generally, the southwestern part of the area is more resistive and compact with a horizontal conductive layer at approximately 1 km depth, with a very thin sedimentary layer on top. The structures north of the BNF are more heterogeneous, with deeper conductive layers (2–3 km) and thicker (several hundred metres) sedimentary layers above. The conductive zones are interpreted as saline water bearing dilation zones rather than conductive meta-sedimentary rocks. The observed fluid transport properties in boreholes and permeability estimation based on Archie's law and fractures indicate a lower permeable zone in the area SSE from NGE1. SUPPORTING INFORMATION Figure S1. Time-series data set example of the electromagnetic field components (Ex, Ey, Hx, Hy and Hz) and its spectrogram from site BUS13. The green line marks data used for processing with less EM artificial noise. Figure S2. Example of processing results for three nights in a row by Smirnov code for site NCW12 without any post-regularization. Figure S3. Example of good sounding curves for site NCW12 and bad sounding curves for BUS13 from three different processing tools (Phoenix commercial processing tool, EMTF package, Maxim Smirnov processing tool). Figure S4. Phase tensor dimensional analysis for 23 most robust points from area. Left-hand panel: all periods with skew angle beta. Right-hand panel: horizontal plot for same periods. Figure S5. Misfits for NC11 2-D model. Differences (RLS) are represented by reduced chi-square values. Figure S6. Misfits for NC12 2-D model. Differences (RLS) are represented by reduced chi-square values. Figure S7. Misfits for NC_All 2-D model. Differences (RLS) are represented by reduced chi-square values. Figure S8. Misfits for 3-D model. Differences (RLS) are represented by reduced chi-square values. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. ACKNOWLEDGEMENTS We would like to thank the Science Foundation Ireland (SFI) for their financial support (Grant No. 10/IN.1/I3022 for IRETHERM to AGJ), Geoserv and GT Energy. The research has been also financed by the program SASPRO (1497/03/01-b) and was cofunded by the People Programme (Marie Curie Actions 7FP, grant agreement REA no. 609427-SK) and cofinanced by the Slovak Academy of Sciences. The work was accomplished with partial support from the grants by the APVV agency by means of projects 16-0482 and 16-0146 and from the Slovak Grant Agency VEGA by means of project 2/0006/19. Sarah Blake, Robert Delhaye and Thomas Farrell helped with fieldwork and data collection. The authors would like to thank Alison Kirkby and an anonymous referee for their very helpful suggestions and comments that improve final manuscript. We thank Volker Rath for comments and suggestions to improve the manuscript. Also authors would like to thanks the Dublin Institute for Advanced Studies and ICHEC for the providing their software and computational resources. Most figures were plotted using the Generic Mapping Tools (GMT, Wessel & Luis 2017). REFERENCES Adam A. , Szarka L. , Vero J. , Wallner A. , Gutdeutsch R. , 1986 . Magnetotellurics (MT) in mountains—noise, topographic and crustal inhomogeneity effects , Phys. Earth planet. Inter. , 42 , 165 – 177 . 10.1016/0031-9201(86)90089-0 WorldCat Crossref Archie G.E. , 1950 . Introduction to petrophysics of reservoir rocks , Am. Assoc. Petrol. Geol. Bull. , 34 ( 5 ), 943 – 961 . WorldCat Asaue H. , Koike K. , Yoshinaga T. , Takakura S. , 2006 . Magnetotelluric resistivity modeling for 3D characterization of geothermal reservoirs in the Western side of Mt. Aso, SW Japan , J. appl. Geophys. , 58 ( 4 ), 296 – 312 . 10.1016/j.jappgeo.2005.05.006 WorldCat Crossref Avdeev D.B. , Avdeeva A. , 2009 . 3D magnetotelluric inversion using a limited-memory quasi-Newton optimization , Geophysics , 74 ( 3 ), F45 – F57 . 10.1190/1.3114023 WorldCat Crossref Avdeeva A. , Moorkamp M. , Avdeev D.B. , Jegen M.D. , Miensopust M.P. , 2015 . Three-dimensional joint inversion of magnetotelluric impedance tensor data and full distortion matrix , Geophys. J. Int. , 202 ( 1 ), 464 – 481 . 10.1093/gji/ggv144 WorldCat Crossref Baba K. , Chave A.D. , Evans R.L. , Hirth G. , Mackie R.L. , 2006 . Mantle dynamics beneath the East Pacific Rise at 17° S: insights from the Mantle Electromagnetic and Tomography (MELT) experiment , J. geophys. Res. , 111 ( B2 ), 1 – 18 . 10.1029/2004JB003598 WorldCat Crossref Bear J. , 2018 . Modeling Phenomena of Flow and Transport in Porous Media , Springer International Publishing AG , p. 742 . Google Preview WorldCat COPAC Berdichevsky M.N. , Dmitriev V.I. , 1976 . Basic principles of interpretation of magnetotelluric sounding curves , in Geoelectric and Geothermal Studies , ed. Adam A. , Akad. Kiad , pp. 165 – 221 . Google Preview WorldCat COPAC Berdichevsky M.N. , Dmitriev V.I. , Pozdnjakova E.E. , 1998 . On two-dimensional interpretation of magnetotelluric soundings , Geophys. J. Int. , 133 , 585 – 606 . 10.1046/j.1365-246X.1998.01333.x WorldCat Crossref Björnsson G. , Bodvarsson G. , 1990 . A survey of geothermal reservoir properties , Geothermics , 19 ( 1 ), 17 – 27 . 10.1016/0375-6505(90)90063-H WorldCat Crossref Blake S. et al. . , 2016 . Understanding hydrothermal circulation patterns at a low-enthalpy thermal spring using audio-magnetotelluric data: a case study from Ireland , J. appl. Geophys. , 132 , 1 – 16 . 10.1016/j.jappgeo.2016.06.007 WorldCat Crossref Boerner D. , Kurtz R. , Jones A.G. , 1993 . Orthogonality in CSAMT and MT measurements , Geophysics , 58 ( 7 ), 924 – 934 . 10.1190/1.1443483 WorldCat Crossref Booker J.R. , 2014 . Magnetotelluric phase tensor evolution , Surv. Geophys. , 35 , 7 – 40 . 10.1007/s10712-013-9234-2 WorldCat Crossref Brown S.R. , 1989 . Transport of fluid and electric current through a single fracture , J. geophys. Res. , 94 ( B7 ), 9429 – 9438 . 10.1029/JB094iB07p09429 WorldCat Crossref Brück P.M. , Potter T.L. , Downie C. , 1974 . The Lower Palaeozoic stratigraphy of the northern part of the the Leinster Massif , Proc. R. Ir. Acad , 74B , 75 – 84 . WorldCat Cagniard L. , 1953 . Basic theory of the magneto‐telluric method of geophysical prospecting , Geophysics , 18 , 497 – 724 . 10.1190/1.1437915 WorldCat Crossref Caldwell T.G. , Bibby H.M. , Brown C. , 2004 . The magnetotelluric phase tensor , Geophys. J. Int. , 158 ( 2 ), 457 – 469 . 10.1111/j.1365-246X.2004.02281.x WorldCat Crossref Campanya J. , Jones A.G. , Vozar J. , Rath V. , Blake S. , Delhaye R. , Farrell T. , 2015 . Porosity and permeability constraints from electrical resistivity models : examples using magnetotelluric data , in Proceedings of the World Geotherm. Congr. 2015 , 19–25 April 2015, Melbourne, Australia . WorldCat Chave A.D. , Jones A.G. , 2012 . The Magnetotelluric Method: Theory and Practice , Cambridge Univ. Press . Google Preview WorldCat COPAC Chew D.M. , Strachan R.A. , 2013 . The Laurentian Caledonides of Scotland and Ireland , Geol. Soc. Lond., Spec. Publ. , 390 ( 1 ), 45 – 91 . 10.1144/SP390.16 WorldCat Crossref DeGroot-Hedlin C. , Constable S.C. , 1990 . Occam's inversion to generate smooth, two-dimensional models from magnetotelluric data , Geophysics , 55 ( 12 ), 1613 – 1624 . 10.1190/1.1442813 WorldCat Crossref Domra Kana J. , Djongyang N. , Raïdandi D. , Njandjock Nouck P. , Dadjé A. , 2015 . A review of geophysical methods for geothermal exploration , Renew. Sustain. Ener. Rev. , 44 , 87 – 95 . 10.1016/j.rser.2014.12.026 WorldCat Crossref Egbert G.D. , 1997 . Robust multiple-station magnetotelluric data processing , Geophys. J. Int. , 130 ( 2 ), 475 – 496 . 10.1111/j.1365-246X.1997.tb05663.x WorldCat Crossref Egbert G.D. , Kelbert A. , 2012 . Computational recipes for electromagnetic inverse problems , Geophys. J. Int. , 189 ( 1 ), 251 – 267 . 10.1111/j.1365-246X.2011.05347.x WorldCat Crossref Fritschle T. , Daly J.S. , Whitehouse M.J. , McConnell B. , Buhre S. , 2018 . Multiple intrusive phases in the Leinster Batholith, Ireland: geochronology, isotope geochemistry and constraints on the deformation history , J. Geol. Soc. Lond. , 175 ( 2 ), 229 – 246 . 10.1144/jgs2017-034 WorldCat Crossref Gamble T.D. , Goubau W.M. , Clarke J. , 1979 . Magnetotellurics with a remote magnetic reference , Geophysics , 44 ( I ), 53 – 68 . 10.1190/1.1440923 WorldCat Crossref Garcia X. , Jones A.G. , 2002 . Atmospheric sources for audio-magnetotelluric (AMT) sounding , Geophysics , 67 ( 2 ), 448 , doi:10.1190/1.1468604 . 10.1190/1.1468604 WorldCat Crossref Glover P.W.J. , 2010 . A generalized Archie's law for phases , Geophysics , 75 ( 6 ), E247 – E265 . 10.1190/1.3509781 WorldCat Crossref Glover P.W.J. , Zadjali I.I. , Frew K.A. , 2006 . Permeability prediction from MICP and NMR data using an electrokinetic approach , Geophysics , 71 ( 4 ), F49 – F60 . 10.1190/1.2216930 WorldCat Crossref Goodman R. , Jones G.L. , Kelly J. , Slowey E. , O'Neill N. , 2004 , Geothermal Energy Resource Map of Ireland - Final Report . COPAC Groom R.W.R.W. , Bailey R.C. , 1989 . Decomposition of magnetotelluric impedance tensors in the presence of local three-dimensional galvanic distortion , J. geophys. Res. , 94 ( B2 ), 1913 – 1925 . 10.1029/JB094iB02p01913 WorldCat Crossref GT Energy, ltd , 2009 . Deep Geothermal Exploration , Newcastle, South County Dublin . Google Preview WorldCat COPAC Hitzman M.W. , 1999 . Extensional faults that localize Irish syndiagenetic Zn-Pb deposits and their reactivation during Variscan compression , Geol. Soc. Lond., Spec. Publ. , 155 ( 1 ), 233 – 245 . 10.1144/GSL.SP.1999.155.01.17 WorldCat Crossref Holland C.H. , Sanders I.S. , 2009 . The Geology of Ireland , Dunedin Academic Press Ltd . Google Preview WorldCat COPAC Jones A.G. , 1983 . On the equivalence of the “Niblett” and “Bostick” transformations in the magnetotelluric method , J. Geophys. , 53 , 72 – 73 . WorldCat Jones A.G. , 2011 . Three-dimensional galvanic distortion of three-dimensional regional conductivity structures. Comment on “Three-dimensional joint inversion for magnetotelluric resistivity and static shift distributions in complex media” by Yutaka Sasaki and Max A. Meju , J. geophys. Res. , 116 ( B12 ), 1 – 5 . 10.1029/2011JB008665 WorldCat Crossref Jones A.G. , 2012a . Distortion decomposition of the magnetotelluric impedance tensors from a one-dimensional anisotropic Earth , Geophys. J. Int. , 189 ( 1 ), 268 – 284 . 10.1111/j.1365-246X.2012.05362.x WorldCat Crossref Jones A.G. , 2012b . Distortion of magnetotelluric data: its identification and removal , in The Magnetotelluric Method: Theory and Practice , pp. 219 – 302 ., eds Chave A.D. , Jones A.G. , Cambridge Univ. Press . Google Preview WorldCat COPAC Jones A.G. et al. . , 2015 . IRETHERM : developing a strategic and holistic understanding of Ireland’s geothermal energy potential through integrated modelling of new and existing geophysical, geochemical and geological data , in Proceedings World Geothermal Congress 2015 , Melbourne, Australia, 19 – 25 April 2015 ., pp. 1 – 6 . Google Preview WorldCat COPAC Jones A.G. , Chave A.D. , Egbert G.D. , Auld D. , 1989 . A comparison of techniques for magnetotelluric response function estimation , J. geophys. Res. , 94 , 14 201 – 14 213 . 10.1029/JB094iB10p14201 WorldCat Crossref Jones A.G. , Jodicke H. , 1984 . Magnetotelluric transfer function estimation improvement by a coherence-based rejection technique , in Proceedings of the SEG Technical Program Expanded Abstracts 1984 , pp. 51 – 55 . WorldCat Jones G.L.L. , Somerville I.D. , Strogen P. , 1988 . The lower carboniferous (dinantian) of the Swords area: sedimentation and tectonics in the Dublin Basin, Ireland , Geol. J. , 23 ( 3 ), 221 – 248 . 10.1002/gj.3350230304 WorldCat Crossref Junge A. , 1996 . Characterization of and correction for cultural noise , Surv. Geophys. , 17 , 361 – 391 . 10.1007/BF01901639 WorldCat Crossref Kalscheuer T. , de los Ángeles García Juanatey M. , Meqbel N. , Pedersen L.B. , 2010 . Non-linear model error and resolution properties from two-dimensional single and joint inversions of direct current resistivity and radiomagnetotelluric data , Geophys. J. Int. , 182 ( 3 ), 1174 – 1188 . 10.1111/j.1365-246X.2010.04686.x WorldCat Crossref Kelbert A. , Meqbel N. , Egbert G.D. , Tandon K. , 2014 . ModEM: a modular system for inversion of electromagnetic geophysical data , Comput. Geosci. , 66 , 40 – 53 . 10.1016/j.cageo.2014.01.010 WorldCat Crossref Kirkby A. , Heinson G. , Krieger L. , 2016 . Relating permeability and electrical resistivity in fractures using random resistor network models , J. geophys. Res. , 121 ( 3 ), 1546 – 1564 . 10.1002/2015JB012541 WorldCat Crossref Ledo J. , 2005 . 2-D versus 3-D magnetotelluric data interpretation , Surv. Geophys. , 26 ( 5 ), 511 – 543 . 10.1007/s10712-005-1757-8 WorldCat Crossref Licciardi A. , Agostinetti N.P. , 2017 . Sedimentary basin exploration with receiver functions: seismic structure and anisotropy of the Dublin Basin (Ireland) , Geophysics , 82 ( 4 ), KS41 – KS55 . 10.1190/geo2016-0471.1 WorldCat Crossref Lowes F.J. , 2009 . DC railways and the magnetic fields they produce-the geomagnetic context , Earth, Planets Space , 61 , 1 – 15 . WorldCat Mardia K.V. , 1972 . Statistics of Directional Data , Academic Press . Google Preview WorldCat COPAC McConnell B. , Philcox M.E. , 1994 . Geology of Kildare - Wicklow , Geological Survey of Ireland . Google Preview WorldCat COPAC McNeice G.W. , Jones A.G. , 2001 . Multisite, multifrequency tensor decomposition of magnetotelluric data , Geophysics , 66 ( 1 ), 158 – 173 . 10.1190/1.1444891 WorldCat Crossref Menke W. , 2012 . Geophysical Data Analysis: Discrete Inverse Theory , Academic Press, Inc . Google Preview WorldCat COPAC Mottaghy D. , Pechnig R. , Vogt C. , 2011 . The geothermal project Den Haag: 3D numerical models for temperature prediction and reservoir simulation , Geothermics , 40 ( 3 ), 199 – 210 . WorldCat Muffler P. , Cataldi R. , 1977 , Methods for regional assessment of geothermal resources , 7 ( 2–4 ), 53 – 89 . COPAC Muñoz G. , 2014 . Exploring for geothermal resources with electromagnetic methods , Surv. Geophys. , 35 ( 1 ), 101 – 122 . 10.1007/s10712-013-9236-0 WorldCat Crossref Muñoz G. , Rath V. , 2006 . Beyond smooth inversion: the use of nullspace projection for the exploration of non-uniqueness in MT , Geophys. J. Int. , 164 ( 2 ), 301 – 311 . 10.1111/j.1365-246X.2005.02825.x WorldCat Crossref Newman G.A. , Gasperikova E. , Hoversten G.M. , Wannamaker P.E. , 2008 . Three-dimensional magnetotelluric characterization of the Coso geothermal field , Geothermics , 37 ( 4 ), 369 – 399 . 10.1016/j.geothermics.2008.02.006 WorldCat Crossref O'Neill N. , Pasquali R. , 2005 , Deep Geothermal Site Charac terisation, Interim Report to Sustainable Energy Ireland . Deep Geothermal Site Charac terisation, Interim Report to Sustainable Energy Ireland, CSA Group Ltd . COPAC Oettinger G. , Haak V. , Larsen J.C. , 2001 . Noise reduction in magnetotelluric time-series with a new signal-noise separation method and its application to a field experiment in the Saxonian Granulite Massif , Geophys. J. Int. , 146 ( 3 ), 659 – 669 . 10.1046/j.1365-246X.2001.00473.x WorldCat Crossref Olasolo P. , Juárez M.C. , Morales M.P. , Damico S. , Liarte I.A. , 2016 . Enhanced geothermal systems (EGS): a review , Renew. Sustain. Ener. Rev. , 56 , 133 – 144 . 10.1016/j.rser.2015.11.031 WorldCat Crossref Park S.K. , Mackie R.L. , 2000 . Resistive (dry?) lower crust in an active orogen, Nanga Parbat, northern Pakistan , Tectonophysics , 316 ( 3–4 ), 359 – 380 . 10.1016/S0040-1951(99)00264-4 WorldCat Crossref Parker R.L. , Booker J.R. , 1996 . Optimal one-dimensional inversion and bounding of magnetotelluric apparent resistivity and phase measurements , Phys. Earth planet. Inter. , 98 ( 3–4 ), 269 – 282 . 10.1016/S0031-9201(96)03191-3 WorldCat Crossref Phillips W.E.A. , Rowlands A. , Coller D.W. , Carter J. , Vaughan A. , 1988 . Structural studies and multidata correlation of mineralization in Central Ireland , in Mineral Deposits within the European Community , pp. 353 – 377 ., eds Boissonnas J. , Omenetto P. , Springer-Verlag . Google Preview WorldCat COPAC Readman P.W. , O'Reilly B.M. , Murphy T. , 1997 . Gravity gradients and upper-crustal tectonic fabrics, Ireland , J. Geol. Soc. Lond. , 154 ( 5 ), 817 – 828 . 10.1144/gsjgs.154.5.0817 WorldCat Crossref Rodi W. , Mackie R.L. , 2001 . Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion , Geophysics , 66 ( 1 ), 174 – 187 . 10.1190/1.1444893 WorldCat Crossref Rothery B.A. , Phillips W.E.A. , 1983 , A report on the structures of the Clane Block . COPAC Schnaidt S. , Heinson G. , 2015 . Bootstrap resampling as a tool for uncertainty analysis in 2-D magnetotelluric inversion modelling , Geophys. J. Int. , 203 ( 1 ), 92 – 106 . 10.1093/gji/ggv264 WorldCat Crossref Siripunvaraporn W. , Egbert G.D. , Lenbury Y. , Uyeshima M. , 2005 . Three-dimensional magnetotelluric inversion: data-space method , Phys. Earth planet. Inter. , 150 ( 1–3 ), 3 – 14 . 10.1016/j.pepi.2004.08.023 WorldCat Crossref Smirnov M.Y. , 2003 . Magnetotelluric data processing with a robust statistical procedure having a high breakdown point , Geophys. J. Int. , 152 ( 1 ), 1 – 7 . 10.1046/j.1365-246X.2003.01733.x WorldCat Crossref Spratt J.E. , Jones A.G. , Nelson D. , Unsworth M.J. , the I. M. Team , 2005 . Crustal structure of the India-Asia collision zone, southern Tibet, from INDEPTH MT investigations , Phys. Earth planet. Inter. , 150 ( 1–3 ), 227 – 237 . 10.1016/j.pepi.2004.08.035 WorldCat Crossref Strogen P. , Somerville I.D. , 1984 . The stratigraphy of the Upper Palaeozoic Rocks of the Lyons Hill Area, Count Kildare , Irish J. Earth Sci. , 6 ( 2 ), 155 – 173 . WorldCat Szarka L. , 1987 . Geophysical aspects of man-made electromagnetic noise in the earth—a review , Surv. Geophys. , 9 ( 3–4 ), 287 – 318 . WorldCat Tikhonov A. , 1950 . On determining electrical characteristics of the deep layers of the earth's crust , Sov. Math. Dokl , 73 ( 2 ), 295 – 297 . WorldCat Weckmann U. , Magunia A. , Ritter O. , 2005 . Effective noise separation for magnetotelluric single site data processing using a frequency domain selection scheme , Geophys. J. Int. , 161 ( 3 ), 635 – 652 . 10.1111/j.1365-246X.2005.02621.x WorldCat Crossref Wessel P. , Luis J.F. , 2017 . The GMT/MATLAB Toolbox , Geochem., Geophys. Geosyst. , 18 ( 2 ), 811 – 823 . 10.1002/2016GC006723 WorldCat Crossref Author notes Now at: Earth Science Institute of the SAS, Bratislava, Slovakia. Now at: Complete MT Solutions Inc., Ottawa, Canada. Now at: Independent Consultant, Cambridge, UK. © The Author(s) 2019. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - A geothermal aquifer in the dilation zones on the southern margin of the Dublin Basin JO - Geophysical Journal International DO - 10.1093/gji/ggz530 DA - 2020-03-01 UR - https://www.deepdyve.com/lp/oxford-university-press/a-geothermal-aquifer-in-the-dilation-zones-on-the-southern-margin-of-a5a9z7E3FD SP - 1717 VL - 220 IS - 3 DP - DeepDyve ER -