Electrical Resistivity Anomalies Offshore a Carbonate Coastline: Evidence for Freshened Groundwater?

Carbonate lithologies host considerable quantities of the Earth's freshwater resources and partially supply a quarter of the global population with drinkable water. In addition, carbonates constitute substantial amounts of the global coastlines, yet it is not known if and how they can sustain freshened groundwater offshore. Here, we use controlled‐source electromagnetic, seismic reflection, and core sample data to derive a lithological model for the eastern margin of the Maltese Islands and identify four distinct resistivity anomalies within the Upper Coralline and Globigerina Limestone formations. The anomalies hosted in the former are likely associated with low porosities, whereas the anomaly within the latter is indicative of pore fluid freshening. Hydrogeological modeling suggests that freshened pore fluids, emplaced during sea‐level lowstands and preserved in low permeability units, are potentially still found within carbonate shelves. However, resource potential is low due to its relict nature and low permeability host environment.

• Geophysical data and hydrogeological modeling are applied to detect offshore freshened groundwater in a semi-arid carbonate setting • Globigerina Limestone and Blue Clay located offshore SE Malta likely host a disconnected offshore freshened groundwater body • The resistive anomalies within the Upper Coralline Limestone are interpreted as localized porosity variations Supporting Information: Supporting Information may be found in the online version of this article.
characteristics of OFG are poorly constrained, and questions regarding their geometry and distribution, as well as the geological controls and timing of emplacement remain unanswered (Micallef et al., 2021).
The bulk electrical resistivity is sensitive to both changes in pore fluid salinity and sediment matrix porosity (Archie, 1942).Therefore, integrative analyses of CSEM with multi-channel seismic (MCS) reflection data (e.g., Bertoni et al., 2020) have proven effective in differentiating lithological units from pore fluid salinity anomalies in siliciclastic margins (e.g., Gustafson et al., 2019).In comparison, freshened groundwater within carbonate margins has rarely been investigated using CSEM and MCS (one of the few studies was presented by Evans and Lizarralde [2011]) and their hydrogeology remains poorly constrained.This is an important gap in knowledge, considering that carbonates represent major aquifers (Laugié et al., 2019), comprise approximately 16% of the global coastline excluding Antarctica (Goldscheider et al., 2020), and outcrop along >15% of the Mediterranean catchment area (Bakalowicz, 2015;Margat, 1998).
In this study, we combine CSEM, MCS, and lithological data with numerical hydrogeological modeling to map freshened groundwater offshore the Maltese Islands and to identify the factors governing its emplacement and understand its resource potential.This archipelago has a semi-arid climate and one of the lowest ratios of water resources per inhabitant globally (FAO, 2006), motivating the exploration for unconventional water resources using a combined geophysical, geological, and hydrogeological approach.

Stratigraphic Framework
The Maltese Islands are located in the central part of the carbonate Pelagian Platform (Figure 1).An Oligo-Miocene sedimentary succession of four formations is exposed across the archipelago (Pedley et al., 1976): Lower Coralline Limestone (LCL), Globigerina Limestone (GL), Blue Clay (BC), and Upper Coralline Limestone (UCL); the limestone formations are further subdivided into members (Oil Exploration Directorate, 1993).
Three types of meteorically recharged groundwater bodies are hosted within the carbonate lithologies onshore.The largest groundwater body is hosted within the LCL and GL formations and comprises a freshwater Ghyben-Herzberg lens floating above saltwater (Bakalowicz & Mangion, 2003).It provides the primary source of potable water in Malta and has a residence time of 15-40 years (Stuart et al., 2010).The other two groundwater bodies have been identified from onshore drilling within the UCL and are either perched over the impermeable BC formation above sea-level, or in valleys below sea-level (ERA, 2015).Onshore recharge of all groundwater bodies is equivalent to ∼25% of the annual precipitation of 550 mm (FAO, 2006;Galdies, 2013).
Our study area is located offshore the eastern coast of the Maltese Islands (Figure 1).The seafloor predominantly consists of gently sloping terrain and represents a submerged paleo-landscape (Micallef et al., 2013), the boundary of which is delineated by fault escarpments, paleo-shorelines, and shore platforms.

Physical Property Measurements
Seven core samples were collected at two sites on Malta and Gozo in 2018 (blue markers in Figure 1).They represent lithological members that most commonly outcrop across the islands (cf., Table 1).The samples were analyzed in terms of porosity, density, hydraulic conductivity, and electrical resistivity to define plausible ranges for interpreting geophysical and hydrogeological models.Additionally, Table 1 lists reported literature data along with P-wave velocities and unit thicknesses measured at borehole BH3 (pink marker in Figure 1a).
Resistivities measured on samples saturated with seawater (0.2 Ω m) range between 1.96 and 4.43 Ω m and tended to increase with stratigraphic age of the corresponding unit.Measured porosities lie within the broad ranges found in literature.Higher P-wave velocities are measured for the UCL and LCL formations compared to the Globigerina members.This stands in agreement with a positive density correlation, but departs from the inverse correlation of the formation porosities.

Geophysical Methods and Results
Five 2D MCS lines (labeled 2, 4, 5, 6, and 9 in Figure 1) were acquired in October 2018 using a single mini-GI gun source with a Geometrics Geo-eel solid-state streamer consisting of four 12.5 m long sections.Data processing included streamer geometry correction, seismic trace balancing, bandpass-filtering, NMO-correction, and 2D-Stolt-migration using a seismic velocity of 1,500 m/s.
During the same cruise, a seafloor-towed CSEM experiment was carried out using a modular, electric dipole-dipole system consisting of a 100-m-long transmitter dipole followed by three inline receivers at offsets of 150, 280, and 510 m, respectively (Figure S4 and Schwalenberg et al., 2017).Data were acquired along five profiles labeled lines 2, 5, 6, 8, and 9 in Figures 1 and 2. Following the shallow-water sensitivity considerations of Haroon et al. (2020), only step-on data were processed at stationary waypoints marked by black triangles in Figure 2. Data errors are quantified by stacking repeated step functions at each waypoint resulting in a time-dependent noise model that varied each day of the survey and also for each receiver (presumably due to electrode noise).Overall, the source-normalized absolute noise floor is quantifiable at 5e −12 V/m at 1 s.Additionally, a 5% minimum relative error floor is assigned to account for systematic distortions due to imprecise navigation data of the seafloor-towed system.The processed CSEM data were interpreted using isotropic models in the time-domain using MARE2DEM (Key, 2016;Haroon, Hölz, et al., 2018).Data fits for each profile are presented in Figures S5-S10.Note that vertical transverse isotropic models lead to consistent results (see Figures S11 and S12) and, therefore, anisotropy was not considered in this study.
HAROON ET AL.
10.1029/2020GL091909 3 of 10 Due to the strong reflectivity of the carbonate seafloor, only limited seismic energy penetrated the sub-seafloor, restricting lithological inferences based on seismic data alone (Figure S14).Thus, seismic interpretations along lines 2, 5, 6, and 8 were supported by data from borehole BH3, morphological observations, CSEM resistivity models, and in-situ data from remotely operated vehicle (ROV) imagery.Lithological inferences are traced above the electrical resistivity models in Figure 2. Along line 9, seismic data had higher quality and was correlated with Gatt (2012), who shows a >60 m thick GL formation with strata tilted 4° overlying LCL.
The derived geological models are jointly interpreted with the electrical resistivity models and core samples to identify anomalous features indicating lithological boundaries, local porosity anomalies, or OFG occurrences.The inferred offshore geology displayed in Figure 2 shows an overburden layer consisting of marl, gravel, and sand, which generally comprises a mixture of hemipelagic and lowstand deposits.Below, a sequence of UCL, BC, GL, and LCL of varying thicknesses exists.ROV imagery and samples collected at the base of the escarpment show outcrops of LCL (Figure 2b).Box canyons, known to occur only within the UCL formation onshore, are visible at the top of the escarpment (Figure 1), suggesting that UCL comprises the upper part of the escarpment.
All 2D resistivity models in Figure 2 achieve a target misfit of RMS = 1 and overlay the interpreted lithology.
The corresponding error models and data fits are found in Figures S5-S10.The overburden is represented by a conductor of ∼1 Ω m (labeled C1 where corresponding to channel infill).Along lines 2, 5, 6, and 8 (Figures 2b-2e), a prominent conductive feature (∼1 Ω m) is also located at a depth of ∼100 m, which corresponds to the BC formation.Above the BC, UCL appears as a more resistive unit (2-4 Ω m) that contains local resistivity anomalies (>6 Ω m).On line 8, a ∼50 m thick resistor of 14 Ω m (R2) is located within UCL toward the SE of the paleo-channel (labeled C3, Figure 2d).On line 6 (Figure 2e), an anomalous high resistor R1 (>15 Ω m) is observed within UCL near the seafloor at 3,800 m.Line 5 (Figure 2c) shows a dipping resistive structure R3 (∼6 Ω m) within UCL between 0 and 2,000 m below an NW-SE oriented graben (Illies, 1981).Across the northern half of the study area (lines 2, 5, 6, and 8), GL and LCL feature below the BC as a deep-lying resistor (>8 Ω m).However, the interface between BC and this resistor is poorly defined due to the limited vertical resolution of unconstrained CSEM inversion and a lack of adequate seismic reflection data.Along line 9 (Figure 2a), tilted strata appear as a resistive structure (6-10 Ω m) underlying a conductive fine sand to muddy sediment layer.A local resistivity anomaly of ∼10 Ω m (labeled R4) lies within GL HAROON ET AL.   (2015).

Table 1 Physical Properties of the Outcropping Lithologies on the Maltese Islands
at a depth of 40-50 m below the seafloor and is bounded below and toward the coast by conductor C2 (1-2 Ω m).At depths >300 m, the deep-lying resistor is associated with LCL.

Hydrogeological Modeling
A shore-normal cross-sectional hydrological model of groundwater flow and solute transport was constructed along line 9 based on the above-mentioned observations (Figure 3b).A scenario in which sea-level is gradually increased from −130 m to 0 over the past 20 ka was considered to assess how OFG emplaced during the Last Glacial Maximum responds to salinization.The aquifer is considered as unconfined and each geological formation is modeled as homogeneous.No layers of impermeable material are located both above and below the aquifer.The water table is at atmospheric pressure, and thus is able to rise and fall depending on external forcing.Three different hydraulic-conductivity anisotropy factors (K h /K v of 1, 10, and 100) are investigated to assess their influence on OFG evolution.The present-day result from the model simulation that is most comparable to the geophysical inversion model is shown in Figure 3b and corresponds to an anisotropy factor of 100.Such an anisotropy is justified for the geological setting and scale we are modeling; anisotropies of numerical models lumping multiple geological layers into a single hydrostratigraphy over km scales can go up to 1,000-10,000, especially for sedimentary rocks (Freeze & Cherry, 1979).Detailed descriptions regarding the hydrogeological modeling are listed in Section S6 and Figures S15-S19.
HAROON ET AL.In this scenario, OFG occurs as a 1 km extension of the terrestrial groundwater system within LCL (−2,000 to 0 m), and as an isolated groundwater body (2,500-5,000 m), predominantly hosted in the low-permeability BC and GL formations.The evolution from the initial freshwater conditions to this final result entails a progressive salinization with the isolated OFG preserved in the lower permeability units (Figure S19).The salinization is slower than the sea-level rise.A decrease in anisotropy results in an increase in the overall salinity of the groundwater and a reduction in its offshore extent (Figures S17-S19).

Discussion
CSEM resistivity models reveal complex internal architectures in the UCL and GL formations, composed of local resistivity anomalies within a background resistivity range between 2 and 4 Ω m.Assuming Archie's relationship (Archie, 1942) for fully saturated carbonate rock (Glover, 2016), a distinct increase of bulk resistivity (ρ) is either caused by a porosity (Φ) decrease, a cementation exponent (m) increase, or a porewater resistivity (ρ w ) increase.The latter is proportional to a salinity decrease (Figure 3a).Note that a cementation increase alone is rather unlikely to cause the localized resistivity variations we observe for R1-R4, whereas variations in porosity and salinity can easily account for the observed resistivity variations (cf., panels of Figure 3a).
Using the reported UCL and GL porosities (cf., Table 1), and assuming a cementation exponent of m = 1.8 (Figure S3), we can estimate that seawater-saturated (38 PSU) GL resistivities lie between 1 and 3 Ω m, and between 1.5 and several hundred Ω m for UCL.Note that this upper limit is due to the extremely low porosity of 0.02 reported in one published study (Cassar, 2010).
The background resistivities for UCL observed along lines 2, 5, 6, and 8 indicate a porous seawater-saturated carbonate matrix, with local resistivity variations that can range between 6 and 15 Ω m (i.e., anomalies R1-R3).R1 (>15 Ω m) is located at the seafloor within UCL along line 6 (Figure 2e).R3 on line 5 (Figure 2c) is interpreted as dipping UCL strata at an inferred graben structure (Illies, 1981).The observed resistivities are 2-5 times larger than the background seawater-saturated UCL resistivities along lines 5 and 6.Both R1 and R3 extend to the seafloor, yet no freshened groundwater discharge was registered by the CTD (conductivity, temperature, depth) sensor (mounted on the CSEM system) and pore water sampling from gravity cores (Berndt et al., 2021).These observations, together with a lack of a confining unit above the resistivity anomalies, suggest that the latter are unlikely to be associated with OFG.Along line 8 (Figure 2d), R2 is embedded in the UCL formation at several tens of meters depth below seafloor to the east of Gozo.Here, an OFG system with local pore fluid freshening is a more plausible scenario, with the permeable UCL interbedded between less permeable confining units above and below.However, similar to anomalies R1 and R3, R2 can also be explained by a local decrease in the porosity, which in UCL can be attributed to the occurrence of UCL-T.
Resistivity anomaly R4 on line 9 (Figure 2a) is found within GL, which exhibits lower variability in porosity compared to UCL.Thus, R4 is indicative of a salinity anomaly rather than a localized porosity decrease.It extends <1 km along the profile and is overlain by BC and muddy sands, which act as confining units.
The porosities of GL and BC lie above 0.24 and 0.35, respectively, implying that local pore water salinity values could be between 3.5 and 10 PSU (converted with the porosity range in Table 1 and Figure 3a).The numerical simulations demonstrate that the geological setting along line 9 is conducive to the emplacement of OFG (Figure 3b), and that the model-derived salinity pattern qualitatively compares with the CSEM inversion model (Figure 3c).The present-day OFG in the hydrogeological simulation occurs as a disconnected body, although it also includes a smaller 1 km extension of the terrestrial groundwater system that was not covered by the CSEM survey.The model simulations suggest that the OFG was emplaced during previous sea-level lowstands, when the offshore geological formations were charged with freshwater by enhanced topographically driven flow and infiltration.As sea-level rose, the freshwater was better preserved in GL and BC due to their low permeability (Figure S19).The disconnection between the OFG body and the terrestrial groundwater system at present is due to salinization of the more permeable intervening units and the thin GL formation (cf., Figure S19).The modeling results suggest that active meteoric recharge from onshore does not occur at present.The preservation of the near-shore OFG body was thus likely enhanced by its shallow occurrence and, accordingly, a shorter duration of the salinization process in comparison to deeper areas.The sub-seafloor geology thus plays a key role in OFG evolution.The preservation of OFG in low permeability units is similar to what has been reported by Lofi et al. (2013) using borehole information from a siliciclastic margin.
At depths >300 m below sea-level, all resistivity models show a resistive formation (5-30 Ω m) labeled as either GL or LCL.CSEM inversion can detect this resistivity increase with depth but sensitivity is low.Yet, it is seen in all models and is considered a robust feature showing resistivities that exceed the core sample resistivity analysis by as much as one order of magnitude.However, if this resistivity increase is also indicative of pore fluid freshening at depth remains unanswered.LCL porosities as low as 0.02 are reported in literature, which provide an alternative explanation for these anomalously high resistivities.
The implications of our study are twofold.First, it highlights the difficulties of mapping OFG systems with geophysical techniques in shallow lithified carbonate settings in comparison to unconsolidated siliciclastic sediments.Our MCS reflection data had poor quality and penetration, hindering the development of a robust geological model.The high spatial variability in porosity within some of the limestone formations, on the other hand, can lead to resistivity anomalies that are difficult to distinguish from pore water freshening.
Constraining the origin of each resistivity anomaly would require high-quality seismic data, borehole drilling, or other ground-truthing data (e.g., evidence of freshened groundwater seepage).Second, our study demonstrates that near-coastal freshened groundwater can exist offshore semi-arid, carbonate coastlines.This is an important outcome in view of the widespread occurrence of similar settings in the Mediterranean Sea and beyond.However, the OFG resource potential in such lithified carbonates should be questioned if the latter has developed in conditions similar to those of the Maltese Islands.Both extraction rates from low permeability units and sustainability of exploiting a relict groundwater system can be considered limiting factors from a resource perspective.

Conclusions
We present a marine CSEM study from a carbonate shelf targeting OFG.By combining 2D resistivity models with seismic reflection profiles, laboratory core log measurements, and borehole data, we derived a lithological model for the seafloor to the NE of the Maltese Islands.Localized resistivity anomalies were found within the Upper Coralline and Globigerina Limestone formations.In view of the reported variability in porosity, the interpretation of geophysical data from UCL is characterized by a high degree of interpretation ambiguity.This, in combination with the occurrence of some of the anomalies at the seafloor and the absence of corresponding groundwater discharge, leads us to attribute the anomalies to a decrease in porosity.
The resistive anomaly along line 9 is found within the more homogeneous Globigerina formation and is more likely indicative of OFG.Numerical hydrogeological modeling suggests that the geological setting along line 9 is conducive to the formation of a disconnected and relict OFG body by emplacement during sea-level lowstands and its preservation in low permeability units during ensuing sea-level rise.OFG can potentially occur offshore the Maltese Islands, as well as in similar settings with carbonate shelves and semi-arid climates.However, in view of the relict nature of the OFG and its occurrence in low permeability units, its potential as an unconventional source of potable water is likely low.This study has demonstrated that CSEM is a suitable geophysical tool for quantitively characterizing potential OFG bodies hosted within carbonate lithologies, although difficulties related to poorly constrained offshore geological models need to be addressed.

Figure 1 .
Figure 1.(a) Geological map of the Maltese Islands with the main formations highlighted by the colors shown in the legend (OED, 1993).Main faults are illustrated by black lines, whereas yellow and orange lines located to the east of the archipelago depict the controlled-source electromagnetic and seismic reflection profiles collected in October 2018, respectively.The white contour marks the paleo shoreline at 20 ka before present.Locations of where core samples were retrieved and the borehole BH3 are illustrated by blue and pink circles, respectively.(b) Close-up image of line 8 located east of Gozo.(c) Close-up image of lines 2, 5, and 6.The red shaded areas mark box canyons.Island topography is displayed by the corresponding color scale in panels (b) and (c).

Figure 2 .
Figure 2. Controlled-source electromagnetic (CSEM) resistivity overlays with geological units derived using seismic reflection (cf., Figure S14) and borehole (BH3) data.Lithological units are illustrated and labeled by the pale color scheme shown in the legend, whereas resistivity models adhere to the color scale.Location of each CSEM line is schematically illustrated on the top right.(a-e) Resistivity models overlying the lithological interpretation for each CSEM line.Triangles mark the stationary waypoints of the transmitter dipole.(b) Lithological classification of borehole BH3 is displayed as an extension of line 2 toward the southwest.ROV imagery was taken at the base of the escarpment along the extension of line 2 toward the northeast.The inserted picture shows the retrieved LCL sample.

Figure 3 .
Figure 3. Illustrations to aid the interpretation of the resistivity models shown in Figure 2. (a) Resistivity relationship for varying pore fluid salinity and matrix porosity based on Archie's equations(Archie, 1942).The white lines depict specific resistivity contours.Lower panels show resistivity as a function of porosity and cementation exponent (m) for (left) seawater (38 PSU) and (right) freshened water (10 PSU).(b) Salinity model in PSU at present state derived by hydrogeological modeling for an anisotropy factor of 100 using offshore geology and past sea-level lowstands.Note that the salinity was converted from g/l to PSU assuming a density of 1,000 g/l.Light colors at >3,000 m depict the possible OFG body hosted within GL and BC.(c) Resistivity model along line 9 for qualitative comparison to the salinity model displayed in (b).Note that the x-axis of (c) was aligned so that 0 m represents the present-day coastline and is therefore shifted in comparison to Figure2a.
Schembri, Debattista, and Theuma [2017]andSapiano, Schembri, and Debattista [2017]; f: ERA[2015]) and combined with measured data obtained from core samples (porosity, hydraulic conductivity, density, and resistivity) and P-wave velocities measured within borehole BH3 using a three-component geophone and hydrophone streamer with an offset of 10.6 m.UCL-T and LCL-X are incomplete due to insufficient saturation during laboratory measurements, indicating low permeability.Computed porosities derived from the ratio of saturated water volume to sample volume coincide well with literature values, whereas measured hydraulic conductivities are several orders of magnitude smaller.Literature values were largely derived from pump tests in the field, whereas our values are from small-scale lab measurements, which may explain the above-mentioned differences.Data descriptions are displayed in FiguresS1-S3.