3D Characterization of a Coastal Freshwater Aquifer in SE Malta (Mediterranean Sea) by Time-Domain Electromagnetics

Abstract

Electromagnetic (EM) geophysical methods are well equipped to distinguish electrical resistivity contrasts between freshwater-saturated and seawater-saturated formations. Beneath the semi-arid, rapidly urbanizing island of Malta, offshore groundwater is an important potential resource but it is not known whether the regional mean sea-level aquifer (MSLA) extends offshore. To address this uncertainty, land-based alongshore and across-shore time-domain electromagnetic (TDEM) responses were acquired with the G-TEM instrument (Geonics Ltd., Mississauga, ON, Canada) and used to map the onshore structure of the aquifer. 1-D inversion results suggest that the onshore freshwater aquifer resides at 4–24 m depth, underlain by seawater-saturated formations. The freshwater aquifer thickens with distance from the coastline. We present 2D and 3D electromagnetic forward modeling based on finite-element (FE) analysis to further constrain the subsurface geometry of the onshore freshwater body. We interpret the high resistivity zones that as brackish water-saturated bodies are associated with the mean sea-level aquifer. Generally, time-domain electromagnetic (TDEM) results provide valuable onshore hydrogeological information, which can be augmented with marine and coastal transition-zone measurements to assess potential hydraulic continuity of terrestrial aquifers extending offshore.

1. Introduction

Groundwater resources in many coastal regions worldwide are currently under stress because of increasing population, agricultural demands, tourism and economic growth. Fresh groundwater in coastal regions may be a resource that can help to mitigate the water scarcity experienced by coastal communities [1]. However, several first-order questions need to be addressed before the fresh groundwater can be used sustainably. There is a lack of understanding regarding the location, nature, and geometry of coastal aquifer systems and their offshore connectivity. Large-scale desalination of seawater is a technologically viable solution, but there are important energy and environmental impacts that must be considered [2]. Terrestrial time-domain electromagnetic (TDEM) methods of geophysical exploration employing a loop or grounded dipole source can be used to explore the onshore component of coastal aquifers that may extend offshore. TDEM methods are useful because they have good depth penetration in saline environments compared to other geophysical techniques, such as ground-penetrating radar (GPR) and electrical resistivity tomography (ERT). TDEM methods are sensitive to electrical resistivity which, in turn, is diagnostic of important aquifer parameters such as porosity, water saturation and salinity [3]. Transient electromagnetic (EM) methods respond to the interaction between an applied time-varying magnetic flux and the geoelectrical structure beneath the transmitter.

In order to interpret transient electromagnetic responses, the key physical mechanism is the induction process governed by Faraday’s law, which is equivalent to diffusion into a conducting medium of an image of the transmitter (TX) loop current. A fundamental overview of the physical principles underlying the EM geophysical method is given elsewhere [4,5], and there are many reviews related to near-surface applications of EM techniques, e.g., [6,7]. The TDEM method has been widely used for groundwater studies [8,9,10,11,12,13,14], including coastal aquifer characterization. The EM methods distinguish electrical conductivity contrasts between freshwater-saturated or brackish-water saturated formations (resistivity ~10–100 $\mathsf{\Omega }\mathrm{m}$, or more) and seawater-saturated formations (resistivity ~1–10 $\mathsf{\Omega }\mathrm{m}$). A hydraulic connection implies a continuous hydraulic pathway beneath the coastal zone, such that an offshore aquifer could be recharged by its onshore counterpart, or depleted by pumping of its onshore counterpart. Terrestrial EM geophysical surveys can provide valuable information about the existence of such connections. This knowledge is important if we are to ensure the long-term sustainability of groundwater resources and it can serve as a valuable constraint for hydrogeological modeling studies. Onshore-offshore connectivity also has implications for possible onshore land subsidence due to offshore drilling and extraction [15,16]. It should also be noted that the electromagnetic method remains largely undeveloped for data acquisition in the important coastal shallow-water transition zone.

This study is part of a multi-disciplinary project that aims to investigate potential onshore-offshore groundwater aquifer connections based on terrestrial and offshore TDEM geophysical data from SE Malta (Mediterranean Sea). Here we utilize 2D and 3D electromagnetic forward modeling based on finite-element (FE) analysis to constrain the 3D geometry of the onshore freshwater body, in this case, the mean sea-level aquifer in SE Malta. To accomplish this objective, TDEM responses were acquired with the Geonics G-TEM instrument. The analysis of the TDEM data generated geoelectrical models that are used to better understand the variable coastal hydrogeology along a short segment of the Maltese coastline. This helps to characterize the potential groundwater resources of the semi-arid, rapidly urbanizing island of Malta. Future work, currently in the planning stages, will involve conducting additional EM measurements both offshore and within the coastal transition zone.

2. Characterization of Mean Sea-Level Aquifer and Study Area

The Maltese Islands, comprising of Malta, Gozo and Comino, are composed of marine sedimentary rocks deposited between the Late Oligocene and Late Miocene epochs [17]. The five sedimentary formations outcropping across the Maltese Islands include, from top to bottom: Upper Coralline Limestone (162 m), Greensand (11 m), Blue Clay (75 m), Globigerina Limestone (207 m), and Lower Coralline Limestone (up to 1 km, of which the top ~140 m is exposed) (numbers in brackets denote maximum thickness) (see Figure 1a) [18,19]. This succession contains a range of lithologies and facies, but overall it is dominated by marine carbonates of shallow water origin. The rock formations exhibit a gentle regional flexure and normal faulting is widespread [20]. The older (Early Miocene) and most widespread system of faults is oriented SW–NE and includes the Great Fault or Victoria Fault, which is ~11 km long and traverses the entire width of the island. A younger system of faults (Late Miocene-Early Pliocene) is present along the southern coastline and often cross-cuts pre-existing faults. The longest of the younger faults is the NW–SE striking Maghlaq Fault. The climate of Malta is semi-arid Mediterranean characterized by a hot, dry summer and a mild, humid winter. The mean annual precipitation is 550 mm, which mainly falls between September and April [21,22].

The Maltese Islands obtain ~55% of their potable water supply from groundwater, while the rest comes from seawater desalination [22]. Aquifers are the primary source of portable water as there is no appreciable surface water streamflow. The mean sea-level groundwater body lies within the pores and fissures of Lower Coralline Limestone (LCL) in the interval where the formation sub-crops at sea-level south of the Victoria fault [23]. The LCL formation is predominantly composed of an algal fossiliferous limestone with sparse corals. The rocks exhibit moderate, irregular or channel-like permeability [24]. The primary porosity of LCL ranges between 7 and 20%, whereas its intrinsic permeability is low (10⁻⁷–10⁻⁹ m/s). Effective porosity and secondary permeability, both of which are dependent on fissures and weathering, have values of 10–15% and 10⁻⁶ m/s [23,25]. The mean sea-level groundwater body is in lateral and vertical contact with seawater. A body of fresh water in the form of a ‘lens’ floats on saline water due to its lower density [26,27]. The thicker part of the lens is situated in the central part of Malta, with its height decreasing towards the coastline. The mean-sea level groundwater is not at rest, but flows horizontally outward from the thickest part. The aquifer is recharged by the infiltration of rainwater in every winter, and groundwater is either discharged offshore at the coastline or else removed by abstraction (pumpage) for agricultural purposes. The mean sea-level aquifer (MSLA) has a mean thickness of 67.5 m and covers an area of >200 km² [28]. This water is mainly abstracted for potable supply and agricultural use. A number of discontinuous perched aquifers with a limited saturated thickness occur north of the Great Fault in the Upper Coralline Limestone above the impermeable Blue Clay, and they are exclusively used for agricultural purposes.

The study site is situated on the SE coast of the island of Malta ~6 km SE of Valletta, the capital city. The elevation of the study site is ~10 m above the sea surface (see Figure 1a). The rocks exposed at the study site consist of the lower member of the Globigerina Limestone formation, which overlies the upper members of the LCL formation. There is a geological well at 3.5 km west of the study site that shows 35 m of Lower Globigerina above 34 m of Lower Coralline and the elevation of the well is 35 m [29]. There is another well drilled for hydrological purposes, located 2.5 km west of the study site, where the top of the water table is 1 m above sea level [26].

3. Methods

This study utilizes the near-surface TDEM geophysical method to determine the geometry and characteristics of the onshore MSLA along the coast at the survey site in SE Malta. The TDEM measurements were carried out using the Geonics G-TEM system consisting of a portable battery-operated transmitter-receiver (TX-RX) console, a TX antenna deployed as 4 turns of a 10 × 10 m square loop of wire laid on the ground, a 0.6 m diameter RX rigid coil with pre-amplifier, and the supporting cables. In field operations, the equipment was deployed as shown in Figure 2. In this study, all soundings were acquired in the 20 time-gate mode, corresponding to investigation depths of 60–100 m. The 30-gate mode with longer acquisition time allowing for deeper exploration was not used. The depth of investigation also depends on the TX power, which is a product of the loop size, current, and its number of turns, in addition to the subsurface conductivity and the RX sensitivity [30].

The operating principle of TDEM is based on the EM induction process. An abrupt shut-off of a steady value of TX current in the wire loop, according to Faraday’s law, generates an impulsive electromotive force (emf) that drives eddy current flow in the conductive earth. After the shut-off, the emf vanishes and the eddy currents start to decay. A weak secondary magnetic field is produced in proportion to the deceasing amplitude of the eddy currents. The multi-turn receiver coil located at the ground surface measures the time rate of change of the decaying secondary vertical magnetic field, the decay rate being diagnostic of the subsurface electrical resistivity.

3.1. Data Acquisition

The geophysical survey procedure was as follows. At each sounding location, the wire loop was laid out on the ground. Then the RX coil with its pre-amplifier was set up in the center of the wire loop, to achieve a central-loop sounding. The portable TX-RX console was set up immediately outside the TX loop for convenience. A ramp-off current was passed through the wire loop using the signal generator in the TX console. The resultant signal received by the RX coil was recorded by the RX console, averaged over several thousand repetitions to improve signal-to-noise ratio. The overall time to acquire each sounding response was ~5 min. Then the TX loop and RX coil were picked up, along with the TX-RX console, and moved forward to the next sounding location. The center of the RX coil represents the location of each sounding, the latter recorded by handheld GPS. The operating frequency, i.e., the repetition rate of the TX on/off cycle, is in the kHz range (i.e., well outside the main power supply at 50 Hz and cell phones at ~1 GHz.) The TDEM method is non-invasive, and no significant environmental disturbance is made to natural flora, wildlife, or agriculture.

The two orthogonal transects acquired in SE Malta comprising a total of 23 TDEM central-loop soundings in July 2018 are marked as black-and-white symbols in Figure 1b. Profile A is oriented from SE to NW along the shoreline with a total length of 150 m, while profile B is aligned from NE to SW with a length of 60 m (Figure 3). The two profiles cross each other at stations A7 and B4, respectively. An additional dataset of 31 soundings (cyan and purple symbols in Figure 1b; see also Figure 3) was added to this area from a second field survey conducted during June and July 2019. One of the 2019 soundings (station 431) was performed at the crossing point of the two previous transects (A, B) to check signal repeatability. The 2019 survey was performed in order to expand the coverage of the survey from the previous year. The rationale for adding more soundings is that the denser spatial distribution of data enables us to better construct a fully 3D geoelectrical model.

3.2. Data Analysis

3.2.1. D Inversion

The 1D inversion of G-TEM transient EM sounding curves is performed using the IXG-TEM software from Interpex Limited (Golden, Co, USA). After importing a data file containing a measured sounding curve, the software generates a consistent 1D smooth model of electrical resistivity vs. depth, based on the iterative Occam regularization method [31]. The user is required to define the minimum and maximum depths, and also the starting resistivity for initiating the model iterates. We seek the 1D inverted model that gives a satisfactory fit to the TDEM data with minimal variation in electrical resistivity between adjacent layers. Such a “smooth” model generally provides a preferable representation of subsurface geoelectrical structures compared to a “rough” model that may fit the data better but contains unrealistically large variations in resistivity between adjacent layers. An example of an inversion to 100 m depth of G-TEM sounding from station A7 is shown in Figure 4. The data points on the left indicate the Earth-response signal recorded by the RX coil as a function of time (in ms) after current is shut off in the TX loop. The red line on the right indicates an initial guess of Earth resistivity, which in this case is a uniform 10 $\mathsf{\Omega }\mathrm{m}$ half-space. The dark green line on the right indicates the calculated smooth depth profile of Earth resistivity, the predicted response of which (continuous dark green curve passing through the data points on the left) best fits the observed response, subject to the smoothness constraint. At this location, a high-resistivity zone of ~80–100 $\mathsf{\Omega }\mathrm{m}$ appears at ~4–22 m depth, underlain by a uniform low-resistivity zone of ~2–3 $\mathsf{\Omega }\mathrm{m}$.

3.2.2. D Forward Model Context

In this study, we used a well-tested, in-house FORTRAN program to compute 1D transient responses based on a finite-radius, inductively-coupled loop source deployed over a layered earth. For such a 1D model, the resistivity changes only in a vertical direction. A series of 1D responses at different frequencies is computed using the well-known frequency-domain analytic solution [32,33]. The transient response is then obtained by taking an inverse Fourier transform of the frequency-domain responses using a Padé summation method [34].

3.2.3. D and 3D Forward Model Context

This study also utilizes 2D and 3D forward modeling of transient EM responses to further constrain the geometry of the onshore geoelectrical structure of the SE Malta aquifer system. The computation of time-harmonic EM responses of aquifer geoelectrical models is performed using a finite-element (FE) analysis of the governing Maxwell equations in the magnetoquasistatic regime. The FE algorithm [35,36] generates a rectangular mesh that is used to discretize buried 1D, 2D and 3D structures by defining rectangular prisms, or slabs, and assigning them certain dimensions, locations, and electrical conductivities (the inverse of resistivity).

In our simulations, the G-TEM transmitter (TX) in ‘vertical dipole’ mode is approximated by 4 turns of a circular current loop with 5.64 m effective radius (equivalent to the in-field-survey of a 10 × 10 m square loop) lying on the air-earth interface at the origin of the computational grid. A single receiver position is assigned to the center of the TX loop to simulate the central-loop configuration. The resistivity model is discretized using 100 × 100 × 100 nodes of a uniform rectilinear mesh with cell-size 0.8 × 0.8 × 0.8 m. The modeling-domain limits are {−40 m, 40 m}, {−40 m, 40 m} and {−20 m, 60 m} in the x, y and z directions, respectively. A typical mesh contributes roughly 4 million degrees of freedom to the finite-element system of equations since there are four complex degrees of freedom associated with each interior mesh node in the $A$, $\Psi$ formulation, described below. The CPU time required to compute a single controlled-source electromagnetic (CSEM) response for a model of this size at one frequency on the Aspen Systems Texas A&M cluster is ~20 min.

The FE formulation is cast in terms of two Coulomb-gauged electromagnetic potentials, namely a magnetic vector potential $A$ and a scalar electric potential $\Psi$. The Coulomb gauge condition is applied, $\nabla ·A=0$. A set of known primary potentials $(A_p,{\Psi }_p)$ is specified, which consists of the analytic expression for electromagnetic induction in a homogeneous formation with ${\sigma }_p$ constant (see [35,36]. Secondary potentials $A_s$ and ${\Psi }_s$ are then defined according to $A=A_p+A_s$ and $\Psi ={\Psi }_P+{\Psi }_s$, in which case the governing equations become

$${\nabla }^2{A}_s-i\omega {\mu }_0\sigma \left(A_s+\nabla {\Psi }_s\right)=i\omega {\mu }_0{\sigma }_s\left(A_p+\nabla {\Psi }_p\right),$$

(1)

$$\nabla ·\left[-i\omega {\mu }_0\sigma \left(A_s+\nabla {\Psi }_s\right)\right]=\nabla ·\left[i\omega {\mu }_0{\sigma }_s\left(A_p+\nabla {\Psi }_p\right)\right],$$

(2)

where ${\sigma }_s=\sigma -{\sigma }_p$ is the difference between the conductivity distribution $\sigma \left(r\right)$ whose response is required and the background value ${\sigma }_p$ whose response is known. The value of electric field E and the induction field B are derived, after calculation of the Coulomb-gauged electromagnetic potentials, according to

$$E=-i\omega \left(A+\nabla \Psi \right),$$

(3)

$$B=\nabla \times A$$

(4)

The spatial derivatives in the above equations are performed numerically in the post-processing stage of the algorithm.

To summarize, Maxwell’s equations are formulated in terms of frequency-domain magnetic vector and electric scalar secondary potentials. The primary potentials are set by the aforementioned analytic solution and added to the calculated secondary potentials in order to obtain the total response at the prescribed frequency. At a given receiver location, such as the center of the TX loop, the total vertical magnetic field component is computed by numerical differentiation of the computed potentials. This procedure is repeated for a number of frequencies spanning several decades, building up the frequency-domain response. For this study, responses are evaluated at 43 logarithmically-spaced frequencies, at 6 frequencies per decade over the range 10¹–10⁸ Hz. After its inverse Fourier transform into the time-domain, the resulting computed transient responses may be directly compared to the G-TEM sounding curves measured in the field.

4. Results

4.1. D Scenario

First we analyze the transient EM soundings from the two orthogonal G-TEM transects of July 2018 comprising 23 locations along and across-shore SE Malta. The field dataset is divided into two transects, labeled A and B. All soundings are plotted in terms of Earth-response voltage as a function of time on a single log-log display for each transect (Figure 5a,b). This format illustrates the variability, or scatter, in the temporal decay of the signals following shut-off in the TX current. A definition of time gate is provided in Appendix A. At station A3, a distinctive and unusual decay curve is observed, which is thought to be caused by effects of localized 3D subsurface structures of unknown origin. This curve, plotted as blue dots in Figure 5a, is clearly distinguished from the other curves and it cannot be fit by the response of a 1-D model. At A3, the unusual response—perhaps from inductive or IP coupling to steel infrastructure—exhibits a sign reversal (from solid to open circles) after gate 13 of the transient and it is not considered for further analysis. The central-loop response of a 1-D layered model cannot generate such a sign change.

After examining the remaining 21 sounding curves comprising the alongshore profile, the responses from the southernmost stations A2, and A4–A9 may be classified as one group since they exhibit very similar decay patterns. A separate 1D inversion was performed for each of these soundings. The resulting 1D resistivity models from each station were used as initial resistivity distributions in an attempt to find a single 1D model that could fit these southernmost soundings. After many iterations of computation and model adjustment using the 1-D analytic forward code, a simple 3-layer 1D model was found to be the most consistent with the field responses (Figure 6, right). This resistivity model for the southern section of profile A consists of a three layered-earth of 5.5 $\mathsf{\Omega }\mathrm{m}$ and 25 $\mathsf{\Omega }\mathrm{m}$ resistivity with 4 and 15 m thicknesses, respectively, and including a basal resistivity of 1.8 $\mathsf{\Omega }\mathrm{m}$. The fit of this model to the sounding curves A2, A4–A9 is shown in Figure 6, left.

We used the same procedure to analyze the sounding curves from all 21 stations comprising the reduced SE Malta 2018 dataset and the additional 31 soundings from the field survey conducted during June and July 2019. Another unusual decay curve is found at station 432 (black squares, Figure 7). At the late-time of this sounding, the observed signal decays significantly slower compared to neighboring stations, i.e., station 433, which is located only 15 m to the north. The anomalous response may be due to the effects of a localized highly-conductive body; this will be discussed later.

After investigating all of the decay curves, we observe certain systematics in the spatial variability in the measured responses. Over distances of a few tens of meters, for example, it is shown below that the lateral changes in subsurface resistivities in the across-shore direction are much stronger than those in the alongshore direction. As regards the locations of soundings and similarity of decay patterns, many of the more recently acquired soundings are similar to those of the earlier acquired transects A and B. For example, consider soundings 430 and 429, which are situated 15 and 25 m south of station B3, respectively (see Figure 3). The responses from these three stations, along with that of station 428, can be sorted as one group due to their similar decay pattern. The best-fitting 1D model of these four stations, whose response is illustrated as the thin black line in Figure 8, consists of a three layered-earth of 5.5 $\mathsf{\Omega }\mathrm{m}$ and 18.2 $\mathsf{\Omega }\mathrm{m}$ resistivity with 4 and 12 m thicknesses, respectively; with the basal resistivity of 1.8 $\mathsf{\Omega }\mathrm{m}$. This model is displayed as the column beneath station B3 in Figure 9b.

From the 1-D forward modeling results, two pseudo-2D resistivity models have been constructed and they are depicted in Figure 9. These models are obtained by merging, or “stitching”, the 1D forward model results from groups of adjacent stations. With regards to the alongshore profile shown at the top, the upper-layer resistivity value is constant and it exhibits no variation in thickness observed along the 150 m transect. In the very near surface, from the surface to 3–4 m depth, the uppermost layer represents a spatial average over a heterogeneous region and we do not attempt to interpret this layer. The second layer spans the depth range 4–19 m in the SE part of the profile, but the layer becomes thinner and slightly more conductive in the NW part. A huge contrast in vertical resistivity variations of maximum 0.1 $\mathsf{\Omega }$m beneath the sounding A1 compared with a neighboring sounding A2 is suggestive of structure with very low resistivity at depth, such as steel infrastructure. The lateral variations in resistivity of geological origin are much stronger in the across-shore transect, shown at the bottom. The top layer of this profile becomes slightly less resistive and thicker towards the coast. In contrast, the underlying resistive zone becomes thinner as the sounding location is located closer to the sea.

In Figure 10, another set of pseudo-2D resistivity models, corresponding to TDEM profiles C and D (see Figure 3), are obtained by combining 1D forward model results from the 2018 and 2019 datasets. These models enable visualization of the resistivity structure in the western and northern parts of the study area. The model from profile C shown at the top (Figure 10a) is located ~30 m west of Profile A. Profile C runs NW–SE alongshore and intersects profile B at station B7 (Figure 3). The resistivity model of this transect appears to be similar to that of Profile A. However, the resistivity values of the second and basal layers are higher in Profile C due to its greater distance inland, i.e., away from the seawater. The second layer of Profile C is also thicker compared to that of Profile A. Lateral heterogeneity of resistivity at depth of ~17.5 m to 60 m can be observed in the SE part of this transect similar to that of Profile A beneath soundings 432 and 433. The across-shore resistivity distribution in the northern part of the study area, labeled Profile D, is shown in Figure 10b. Profile D is located 90 m northward from the intersection of transects A and B. There is no significant change in either the thickness or resistivity of the uppermost structure of this profile as compared to Profile B. With respect to the middle, resistive layer along Profile D, the shape is comparable to the resistive zones found in Profile B, except they are less resistive and somewhat thinner.

4.2. D and 3D Scenarios

In the previous section, we used the analytic solution of the TDEM forward problem to determine stitched 1D resistivity depth-profiles across the SE Malta study area. In this section, we use the FE analysis to compute frequency-domain responses of 2D and 3D models. The time-domain response is obtained by splining the frequency-domain response evaluated at each of the designated discrete frequencies. Subsequently, the set of time-domain responses are used as the input from which we develop a series of 2D and 3D forward model iterative adjustments. The best 2D and 3D models that result from this analysis are then further evaluated and interpreted. The adjustments are made by trial and error since insufficient computational resources are available to achieve an automated inversion process. Since a single forward run takes ~14 h of CPU time on our computational platform, and an automated inversion would require many thousands of forward runs, even with a highly efficient algorithm it is envisioned that both coarse-grained and fine-grained massive parallelization are a prerequisite for a fully 3D inversion. Such algorithmic development is beyond the scope of this study, but is definitely recommended for future work.

Figure 11 shows FE-calculated responses at two stations based on the fully 2-D model constructed from the stitched 1-D resistivity models shown in Figure 9. The calculated response at station A7 obtains from the alongshore 2D resistivity model in Figure 9a. This model allows spatial variations in resistivity only in the SE-NW direction. That criterion is kept for all soundings along Profile A. Similarly, the 2D lateral resistivity distribution used to compare with the observed soundings at each station along Profile B is based on the across-shore transect shown in Figure 9b. The yellow dots in Figure 11, left, represent the field response actually measured at station A7. The modeling result of the alongshore 2D structure, computed using the 3D FE code, is marked as the solid line. Another sounding at the intersection of the two transects, namely station B4, is displayed in green diamonds in Figure 11, right, with the corresponding across-shore 2-D model response shown by the solid line.

At the bottom of Figure 12, the step-off voltage response as a result of 3D forward modeling was computed at station locations A7 and B4. The 3D model shown at the top of Figure 12 is constructed by combining the 2D models from the three transects, namely Profiles A, B (in Figure 9) and C (Figure 10a). For the sake of better visualization, only a local portion of the complete 3D model that is indicative of the structure beneath station A7 is illustrated within the modeling-domain limits of {−40 m, 40 m}, {−40 m, 40 m} and {0 m, 60 m} in the x, y and z directions, respectively. Part of the model that is above ground surface up to 20 m high is also excluded for better visualization; the size of 10 × 10 m square TX loop is shown for scaling. The complete 3D model representing the subsurface structure beneath SE Malta, covering a surface area of 16,500 m², is depicted in Figure 13. Some of the sounding points are included to better indicate the location and orientation of the 3D model with respect to the G-TEM survey layout. A brief sensitivity analysis is provided in Appendix B.

The computed misfit at each sounding location is plotted in terms of relative error, visualized using various circle sizes, for the 1D, 2D and 3D models. These misfit circles are shown in black, blue, and red, respectively (Figure 14). The misfits of the responses at gate 1 and 2 for all soundings are excluded from the display since the amplitude of the early-time responses is very large. The relative errors of the 2D model are shown only at the 21 stations of the reduced SE Malta 2018 dataset located along transects A and B. The reader should note that the misfit of the preferred 3D model at a given station may exceed the misfit of the 1D model at that station. The important point is that a single 3D model has been found that fits all the observations reasonably well, sometimes at the cost of locally increasing the misfit compared to a 1D model that strictly applies only to an individual station. The actual geoelectrical structure of the Earth is 3D rather than locally 1D beneath the G-TEM measurement stations.

5. Discussion

The results presented in this study suggest that the structure of the mean sea-level groundwater aquifer near the shore in SE Malta exhibits a lenticular shape, with decreasing thickness towards the coast. The combination of TDEM models derived from the summer 2018 and 2019 datasets shows distinct high-resistivity zones. These are interpreted as the signature of a brackish water-saturated geological medium, in this case corresponding to the LCL formation hosting the mean sea-level groundwater body.

Based on the preferred 3D TDEM model in Figure 13, the top layer up to 5 m deep is interpreted as the overlaying Globigerina Limestone and there is no geophysical indication of freshwater in this low-resistivity formation. We also find that the depth to the top of LCL and water table in the study area is 4–5 m. Our study is in good agreement with a regional groundwater modeling study (MARSOL, 2015) of the South Malta region which points out that the elevation of the top of the formation holding fresh groundwater is in the range +20 to −20 m with reference to mean sea level [37]. The zones of high resistivity below the depth of 4–5 m in the 2D and 3D models are indicative of a (moderately brackish) freshwater-bearing formation with resistivity in the range $\rho$ ~10–100 $\mathsf{\Omega }$m (i.e., the purple-blue-cyan colors in Figure 13). The steeper base compared to the gentler top of the groundwater body is consistent with the geometry of a Ghyben-Herzberg-type lens. Below the freshwater region, from depths ~13–22 m down to the TDEM depth of investigation at 60 m, the underlying rocks situated beneath sea level are much less resistive, attaining values $\rho$ ~1.25–2.5 $\mathsf{\Omega }$m (i.e., the green-yellowish green colors in Figure 13). These low resistivity (conductive) zones are indicative of a seawater-saturated formation. The shallow resistive freshwater lens sits on top of a more conductive formation, the latter being indicative of lateral landward movement of saltwater, i.e., intrusion. The boundary between the zones of high and low resistivity indicates the presence of the interface or transition zone along the two across-shore transects. Within the areas closest to the shoreline a mixing zone of freshwater and seawater appears to be present. Zones of moderate resistivity $\mathsf{\rho }$ ~5 $\mathsf{\Omega }$m are observed along the northeastern parts of the across-shore transects by the coast and this could be indicative of brackish groundwater.

In order to assess the groundwater quality implications of our model, we calculate the bulk resistivity of the fluid-saturated rock using Archie’s law [38] for various porosities of limestone assuming that all pore spaces are filled with freshwater with resistivity of 2 $\mathsf{\Omega }$m. This latter value is equal to the water resistivity found in a well located 2.5 km inland. For porosities of 10% and 15% we find 126.2 and 60.8 $\mathsf{\Omega }$m, respectively, as the formation bulk resistivity. These values of estimated resistivity are slightly higher than the range of those in the 3-D TDEM model ($\mathsf{\rho }$ ~10–100 $\mathsf{\Omega }$m). In addition, we have estimated bulk resistivity for saturated limestone filled with seawater of 0.2 $\mathsf{\Omega }$m with 10%, and 15% porosities. These are 12.6 and 6.1 $\mathsf{\Omega }$m, respectively, which is also slightly higher than our model’s prediction of low resistivity seawater-saturated formations of 1.25–2.5 $\mathsf{\Omega }$m. As we move inland, our values are consistent with the borehole’s fluid resistivity saturating a formation of 10 to 15% porosity. Closer to the shoreline, the TDEM bulk resistivity is lower, reflecting more brackish water. Thus, the groundwater freshens as we move inland. Of course, Archie’s law is not a perfect petrophysical model for the fractured limestone lithology, since the law was founded on lab measurements made on clean sandstone cores, but an Archie-type calculation should be approximately correct.

To assess confidence in the spatial structure of our model, we also consider which of the model slabs indicative of the freshwater-bearing formation are best resolved based on the sensitivity analysis (see Appendix B for details). At the lower frequency of 100 Hz, the best-resolved slab is slab 3; whereas the responses from slabs 5, 6 and 7 are more sensitive to perturbations in their resistivity than the responses of slabs 1 and 2. At the high frequency 1 MHz, the misfit-change distribution indicates that changing a slab’s resistivity affects only the sounding that is situated directly over that slab. Slab 6 seems to be the most well-resolved slab at the intermediate frequency 31.6 kHz. Slabs located further inland appear to be not as well resolved as those closer to the sea. The latter are thin and more conductive relative to the thicker, more resistive inland slabs and it is well known that terrestrial TDEM better resolves thin conductive layers.

We do not have sufficient data coverage to infer a possible offshore extension of the freshwater aquifer at the SE Malta study site. The landward encroachment of seawater decreases the resistivities of the near-coastline region and possibly interacts with the fresh groundwater of the MSLA. There are unpublished ground-penetrating radar (GPR) data that appear to show infiltrated meteoric water trapped in fractures above water table in some areas of the study site [39]. More across-shore measurements throughout Malta are recommended in order to investigate the lateral subsurface geoelectrical variation in the direction perpendicular to the coast.

6. Conclusions

This study demonstrates the utility of the TDEM geophysical method along with 1D, 2D and 3D forward modeling as a means to study coastal freshwater aquifers in water-scarce regions. Here we image the geometry of the onshore aquifer within the permeable Lower Coralline Limestone formation along the SE Malta coast. Our results show 2D and 3D resistivity models found by iterative adjustments of FE forward modeling. The final preferred 3D model provides information to depth of 60 m, covering an area of ~16,500 m² and shows diagnostic spatial variations in subsurface electrical resistivity. The geophysical modeling provides a basis for determining important characteristics of the MSLA that fit our observations, namely the decreasing thickness of fresh groundwater bodies towards the coastline. Zones of fresh groundwater have been identified, but these are located preferentially inland from the coast. Thus, there is no indication from the electromagnetic data of a robust offshore extension of the MSLA at this location. However, it is argued that method that we used can be applied across the entire Maltese archipelago to better constrain the geometry, dimensions and distribution of terrestrial and coastal aquifers providing valuable information for future water management of the stressed groundwater reserves of Malta.