A Transient ElectroMagnetic (TEM) Method Survey in North-Central Coast of Crete, Greece: Evidence of Seawater Intrusion

Seawater intrusion into near-shore aquifers is one of the main environmental problems that affect Mediterranean islands. Crete is the biggest and most populated island of Greece, characterized by limited surface waters and strong dependence on groundwater sources as the primary source of natural water supply for extensive agricultural activity and human use. Freshwater demand in Crete has increased notably the last decades. The Geropotamos aquifer is located on the north-central coast of Crete and freshwater management is in a delicate balance with saltwater at coastal areas of the aquifer due to the scarce precipitation and high evaporation as well as the intense over exploitation of the groundwater resources. The geological setting of the study area is considered complex and the local tectonic regime is characterized by two sets of faults orientated NW-SE and NE-SW. Investigation of the aquifer using a survey grid of 1179 Transient ElectroMagnetic soundings (TEM) in 372 sites, has resulted in 1D models, and 2D/3D visualization of geoelectric structures, depicting the zones of salination of groundwater in the aquifer. Geological mapping, hydro-lithological data and geochemical analysis of 24 water samples (22 boreholes and 2 springs) are in agreement with results obtained from TEM soundings, supporting our interpretation that the aquifer is degraded by saline intrusion which likely occurs along fractures in a fault zone, emphasising the critical role of fracture pathways in salination problems of coastal aquifers.


Introduction
Recent climate changes caused by both natural processes and human activities highlight the need to secure freshwater supplies in Mediterranean islands, which are affected due to their geographical isolation and climate. A place where seawater intrusion into the aquifer is of great importance is the island of Crete, Greece. Crete has undergone a rapid and substantial economic development over the last years due not only to the long-established tourist industry but also to changes in the pattern of agriculture, placing increasing pressure on the finite and limited water availability on the island. This research investigates the extent of seawater intrusion in the Geropotamos aquifer, on the coast of north-central Crete.

Geological Setting
The study area called Geropotamos, covering approximately 40 km 2 , is located on the north-central coastline of Crete ( Figure 1) where the local water demand for both irrigation and tourism significantly increased in recent years. The study area lies in the northern part of Geropotamos basin (drained by Geropotamos River) and it is characterized by a combination of lowland and semi-mountainous terrain. The elevation range is 0-502 m, average about 106 m. Geropotamos basin consists of 3 sub-basins: Peramatianos, Anogianos and Zonianos sub-basins. The study area is situated at the northern part of the Peramatianos sub-basin ( Figure 2). The geology and complex tectonic setting of Crete is governed by its position in the central fore-arc of the Hellenic Subduction Zone (HSZ). It is characterized by pre-Alpine and Alpine rocks (composing a pile of nappes) and Neogene and Quaternary sediments (Table 1), which formed in various episodes prior, during and after Alpine orogenesis. The geology of Geropotamos is composed mainly of some of the oldest rocks in Crete (Phyllite/Quartzite and Plattenkalk nappes), covered by Neogene and Quaternary deposits [20]. Table 1. The nappe pile of Crete and sediments (data used extracted from [21]). The nappes and sediments that are present in the study area are highlighted and only rock types occurring in the study area are shown in the last column.
At the northern part near the coastline lies the Phyllite/Quartzite nappe (C3-T3). This unit is composed of low-permeability rocks (Phyllites and Quartzites) ( Figure 2). According to the published N-S geological cross-section in the north-central part of Crete, nearby the study area, the thickness of this nappe is approximately 500 m [22]. At the eastern part of the study area the lower slopes of the Talea Ori Mountain are situated, consisting of dolomitic limestones and dolomites in the southern part, and units of marble and dolomitic limestones in the northern part (Plattenkalk nappe). The Phyllite/Quartzite nappe overlies the Plattenkalk nappe with thrust contact, where the Phyllite/Quartzite nappe moved northwards under Alpine compression. Pliocene and Pleistocene marine deposits overlie Neogene (here only Miocene) biogenic limestones, marls, clays and conglomerates. The rocks crop out in the central and the southern part covering the largest portion of the study area. Holocene alluvial deposits are present at the top of the stratigraphic column. Figure 1 depicts the nappes and key sedimentary deposits in simplified form for clarity. The sediments and rock types of every nappe in the study area are listed in Table 1. The local tectonic regime of the study area is mainly characterized by normal faults of NW-SE and NE-SW directions [23,24] (Figure 1).
The hydrogeological setting of Geropotamos has been poorly investigated and some assumptions (especially about the salinization of the groundwater) are still unconfirmed. The northern part of Geropotamos basin was the subject in [25], which forms a small part of the study of all of the geology of Crete. In the northwestern part of the study area a saltwater spring (known as Salty Spring, SL-SPR) emerges within 200 m of a freshwater spring (known as Sweet Spring, SW-SPR) ( Figure 1) and demonstrate the complexity of underground water pathways. Analysis of the saltwater spring led to an interpretation that evaporites (halite) are part of the Neogene sequence [26]. The two reasons for this interpretation are: (a) the assumption that the Phyllite/Quartzite nappe that lies near the coast was considered by these authors [26] to be a barrier to seawater intrusion ( Figure 2); and (b) the presence of gypsum in Pigi-Adele region 12 km west of the study area, that was theorized to influence the groundwater quality. However, these ideas were not supported by evidence of the evaporites influencing the water supply and there has been no further study until the present research [18]. Recharge of the study area's aquifer is expected to come from the mountains, i.e., Talea Ori, in the eastern and southern part of the area, where highly permeable formations (dolomitic limestones) act as groundwater conduits for rainfall in the mountains. Leakages from the dense river network in the basin supplement the groundwater recharge. Two different types of aquifer are expected based on rock types: (a) a shallow aquifer, in porous media (Quaternary and Neogene sediments), and (b) a deep aquifer in fractures in rocks of the Phyllite/Quartzite nappe and karstic rocks of the Plattenkalk.
Existing data from 22 public boreholes were collected ( Figure 1). The wells are divided into two groups: 9 holes for potable water and 13 for irrigation as classified by the municipality of Geropotamos [27]. According to those data, the groundwater of boreholes used only for irrigation is characterized by high salinity and therefore high conductivity. Note that the logs from these available boreholes are the only available geological information for the subsurface, so they are used for calibration and confirmation of the geophysical modelling for the present study.

Geological Cross Sections
Five simple geological cross-sections were constructed for better comprehension of the geological setting of the study area: AB (9 km), CD (4 km), EF (4 km), FG (7 km) and HI (4 km) ( Figure 3). The data sources that were used for their construction were: (a) the 1:50,000 IGME map [18], (b) information from logs, (c) data collected after geological observation where all rock formations of the study area were recognized, located by GPS, and rock samples were collected. A key feature is the contact (thrust fault) between Phyllite/Quartzite and plattenkalk nappe at the north-eastern part of the area ( Figure 1). Concerning the five cross sections, above each borehole is situated its name and its depth. Boreholes marked with "?" continue below the level drawn but the amount of borehole below that level is not known. The boundaries of rock units are the best estimates from the map and borehole data. Solid and dashed black lines indicate certain faults (CF) and inferred faults (IF) respectively, as derived from the IGME map. There is no information about the lengths of the faults (from the geological map or from fieldwork), so they are drawn as short lines. A geological cross-section for AG was not constructed since this section is exclusively within the Phyllite/Quartzite nappe.    Figure 1 (see text for more explanation). A geological cross-section for AG not shown as this section is exclusively within the Phyllite/Quartzite nappe.

The TEM Technique
Transient ElectroMagnetic (TEM) Technique is relatively young in comparison with other EM methods, as it has been developed and refined mainly in the mid-1980s. An analytical presentation of the method can be found in [28]. The TEM technique can be applied in many different configurations, no direct electrical contact with the ground is required and the depth range in which it can be applied is from the top few meters to hundreds of meters. It is a fast and cost effective method, but it does not work properly in highly resistive rocks and is susceptible to interference from nearby (buried or not) metal pipes, cables, fences, vehicles and induced noise from power lines. For the acquisition of the TEM data, the TEM-Fast 48 system from AEMR Ltd. was used [29]. It is a portable, fast and robust geophysical instrument providing efficient operation in several environmental noisy conditions, giving solutions in many environmental and mineral explorations. TEM data were collected over a period of 3 months (summer 2008), a total 1179 TEM soundings were acquired in 372 different locations in a detailed survey grid (about 200-250 m in EW and NS directions) using a single loop 50 × 50 m (size of the antenna) in order to be able to image a possibly complex subsurface ( Figure 1). In addition to the TEM-Fast 48 system and the antennas (2 × 100 m copper cable), for the acquisition were used: (a) a palmtop PC, (b) an external power supply (12 V + 12 V in series), and (c) two sighting compasses in order to install correctly the 50 × 50 m square. The measurements were acquired with "Time parameter" (range of measurement of transient characteristics) set to 5 or 6 (32 or 36 time gates), which correspond to 1024 µs or 2048 µs Max Time (one complete pulse is divided to 48 time gates and corresponds to 16,384 µs). This parameter is essential because high values correspond to distorted results and low values can result in loss of information about the late stages of transients. In addition, "Stack parameter" (digital stacking of a signal) was selected to 5 viz 65 complete cycles (one unit of the Stack parameter corresponds to 13 complete cycles), in order to improve quality of the results of measurements in case of white noise existence. The current in the generating loop was set to 4A and the frequency of the filter was set to 50 Hz.
The large number of the TEM soundings reflects the fact that the measurements were repeated several times at each sounding location for optimisation of the raw data. Some of those acquired data were discarded, as they were affected by resonance of power lines, by strong Super Paramagnetic (SPM) effect, or rarely by aliasing effect from HF noise. Finally, data from 350 out of the 372 locations were used, the remaining 22 locations being unreliable for interpretation. The 350 collected TEM data were filtered and smoothed prior to the final modelling. The output of a TEM instrument consists of an approximately power law voltage decay curve. The way to interpret these data is to convert the field data to apparent resistivities, p(t). TEM-RES software, available with the TEM-Fast 48 instrument, is a Windows integrated software system and was the tool for data processing of TEM data and then for the inverse problem solution (1D modelling). Due to high quality of the data, almost no editing was required (exclusion of points, change position), except the head correction of the files (Name, Coordinates, etc.). The next step before modelling was the smoothing, as strongly suggested by the manufacturer. Two more curves are given as result of the raw data smoothing. The first (upper) continuous curve approximates initial ρ a (t) data and the second (lower) curve is the calculated dependence ρ full (t) (Figure 4).
Once the resistivity data were derived, the next stage was modelling of the data, as the ultimate goal was to deduce the subsurface resistivity distribution (resistivity-depth information) from these sounding curves. TEM soundings are commonly used to define aquifer properties and other subsurface characteristics. Although those kinds of structures are multi-dimensional, 1D inversions are used to determine geological structure for groundwater flow models. TEM-RES software was applied for the inversion procedure. The formulae that TEM-RES software uses for this procedure are not described in this paper and can be found in [29]. The number of layers was up to 5 and the thickness was not fixed. For better interpretation, after manufacturer's recommendation, an upper threshold 1000 Ω·m was set (resistivity limit). The expected depth of investigation of the soundings is about 150-200 m [29]. As already mentioned, the raw data quality was very good and the RMS errors after inversion was 0.5%-5%. Since the 1D modeling is not enough to reconstruct and describe the subsurface, 2D visualization (1D interpolation) is required. TEM-RES software is used for the geoelectrical sections construction. The same positions as for the geological cross-sections were used to allow for direct comparison between the geology and TEM results. 2D geoelectrical sections are the colour presentation of the structure.
The final 1D inverted models from all the available TEM data were finally merged to form slices differentiated with depth (Layers A, B, C and D). The roof and the bottom of each layer [for example: for layer A, roof = 0.0 m (surface) and bottom = 20.0 m (20 m deep from the ground surface)] are defined using the TEM-RES software. Thus, the selected layers are extracted, where each layer includes the information of the average resistivity (in Log Ω·m) of each point. For example, for layer A (0.0-20.0 m depth from the ground surface) for each of the 350 points an average resistivity value was extracted (for these 20 m only). Then, all this information was inserted in a GIS environment for spatial distribution map construction (one map for each layer) using ordinary kriging interpolator.

Geochemical Analyses
As already presented in previous work [19] in which detailed hydrogeochemical investigation was conducted, it was shown that seawater and the extensive human activities are the main sources of groundwater contamination. Along with the TEM geophysical data, more groundwater chemical results are presented in order to check the quality of the waters and to test the hypotheses of evaporite influence or seawater intrusion in the area.
Three periods of water samplings during 2008 were: (a) first week of June, (b) end of July, and (c) 1st of October, that is to say approximately every 2 months. Chemical analyses of the samples are presented in [19]. These water samplings covered the period from the beginning of the dry season until its peak (October), in which the aquifer undergoes the highest stress. It is noted that during the sampling period high pumping rates for water supply take place (irrigation and touristic season). Generally, rainfall stops during May in Crete. October is the beginning of the rain season. Consequently, in the case of the seawater intrusion problem in Geropotamos study area, the worst groundwater quality would be expected to be found in the October samples before the rainfall can diminish the effects of seawater intrusion.

TEM Approach
As already mentioned, 350 good quality apparent resistivity curves were finally extracted. Then processing and inversion of the data followed in order to receive the resistivity-depth information. The 2D visualization is presented in Figure 5 (sections AB, AG, CD, EF, FG, and HI). The sections that were chosen for 2D visualization lie along the geological cross-sections ( Figure 3) for direct comparison between the geology and resistivity results except for geoelectrical section AG (7 km long, parallel to the coastline) where no geological cross-section was constructed since this section is exclusively on the Phyllite/Quartzite nappe. Note that HI geoelectical section is part of the HI cross section as marked in Figure 3 with green arrows. Hot colours (red) represent high resistivity formations, while the cold colours (blue) depict the low resistivity units. The 3× exaggeration of all the sections is the same as that used for the construction of the geological cross-sections (except for the HI and AG where the exaggeration is ×1 and ×4 respectively). For all the profiles the resistivity scale ranges between 1 Ω·m and 1000 Ω·m (resistivity limit set by the analyst).
3D visualization of resistivity data (spatial distribution map of 1D inversion results) is shown in Figure 6 (Layers A, B, C, and D), where geoelectrical interpretation of intervals to a depth of 200 m is presented. Hot (red) colours represent high resistivity formations and the cold (blue) colours depict the high conductivity. For slice A the nappes layer is superimposed, depicting the surficial geology, while for slice C faults are added. Figure 11 illustrates the interpolated maps of TEM with heights obtained by DEM.

Water Quality
As presented in previous work [19] groundwater chemistry data were obtained in order to check the quality of the waters and to test the hypotheses of evaporite domes or seawater intrusion in the study area. Spatial distribution maps were constructed using the results of the determination of the chemical parameters in the samples of the October 2008 sampling period (ordinary kriging method). Representative maps of four parameters (EC, TDS, Na, Cl) are shown in Figure 7. Figure 8 shows maps of water quality distribution and the type of water of each sample location was characterized (based on the dominant ion), as well as Durov diagram (Figure 9) was constructed in order to identify the source of salinization (seawater or halite). The red dotted line represents the marine composition line (mixing line), joining the end-members-fresh water (SW-SP) and Aegean seawater (green and blue star respectively). It is noted that the ion concentrations of the Aegean seawater were taken from [30].
According to the main constituents (dominant ions), different groundwater types can be met classed into different types for the purposes of this paper, according to dominant dissolved components. From the 24 samples collected in October, three types of groundwater were determined: NaCl type (blue colour)-13 samples, Ca(HCO3) 2 type (green colour)-9 samples, and CaCl 2 (red colour)-2 samples. Generally, Ca(HCO 3 ) 2 type of water is considered to be fresh water in coastal areas due to calcite dissolution, while NaCl type of water is considered salty, as in seawater Na + and Cl − are the dominant ions. On the other hand, CaCl 2 type of water is characterized as a transition water type [31]. % start a new page   In the context of the special interpretation of the results of the geochemical analysis, a groundwater quality map was generated. This quality map, known as Groundwater Quality Index (GQI) [32] links seven chemical parameters (Ca 2+ , Mg 2+ , Na + , Cl − , SO 2− 4 , NO 3 − and TDS) to the WHO guidelines [33].
For each used concentration, X , at every point we construct the index C following the World Health Organization, Guidelines for drinking water quality [33] where C is given by the following expression [32]: C = X −X X +X , X corresponds to WHO guideline value. It is obvious that the contamination index values ranging from −1 to 1.The next step is to generate the "rank map" using the following polynomial function: r = 0.5 × C 2 + 4.5 × C + 5 [32], where C represents the contamination index value, and r represents the corresponding rank value. Afterwards, the relative weight (w) of each parameter had to be calculated as the mean value (r: 1-10) of the corresponding rank map and the "mean r + 2" (r ≤ 8) for the parameter that has potential health effects, i.e., nitrate. The GQI map concerning the seven chemical parameters was constructed using the calculation [32]: GQI = 100 − r 1 w 1 + r 2 w 2 + . . . + r 7 w 7 7 Figure 8 comprises the GQI map of the study area. Hot colours indicate the good water quality, while the cold colours indicate the bad (maximum groundwater quality = 100). A scatter plot of the GQI with the resistivity data of Layer D ( Figure 6) is presented in Figure 8 as well. The correlation coefficient R that occurs is 0.3848.
The data presentation is completed by Figure 10, showing correlation between the geological cross sections and geoelectric profiles. Finally Figure 11 illustrates a three-dimensional view of the distribution of resistivity of the rock mass in the study area, highlighting the low resistivity in association with the location of the coastal area.  Figure 11. This image mainly was constructed in order to make clear the physical mechanism of the seawater intrusion. From the physical point of view, seawater intrusion cannot occur above the sea-level, so the NW part of the study area was contaminated firstly through fracture zones (low elevation), while the salinization of the SE part of the study area is shown only in Layer D.

Data Interpretation and Discussion
There are three possible scenarios to explain the pattern of the saline water in the study area: Scenario 1: seawater intrusion, Scenario 2: the presence of evaporites (halites), Scenario 3: combination of seawater intrusion and evaporites.
To assess the relationship between the groundwater salinity in the context of the geomorphology, lithology and geological structure, we examined the potential impact of evaporite deposits in the study area, discussed below, along with the geological and geomorphological aspects. It was previously suggested that buried evaporites are present in the study area [25,26]. The only place that saline water affects the aquifer is the northern part of Geropotamos basin where seawater could be intruding. The lack of evaporites in observed Neogene sediments in the southern part of the basin is evidence that evaporites probably do not exist in the study area. Note, however, that evaporites (gypsum) exist in the neighbouring basin 12 km west of the study area. Overall, this analysis leads to the deduction that there is good evidence of raised salinity in the aquifer due to seawater intrusion and not due to evaporite dissolution and, therefore, intrusion is considered the most likely interpretation.
Hydrolithological classification of the geological units based on permeability was based on previous study by [26]. According to that classification, the northern part of Geropotamos basin (where the Phyllite/Quartzite nappe lies) is characterized as almost impermeable rock (Figure 2), whereas the eastern and southern part (dolomitic limestones) is shown to consist of permeable formations. Previous ideas of an impermeable formation along the coastline (Phyllite/Quartzite nappe) forming a barrier to seawater intrusion, is contrary to the results of this study. Thus the tectonic setting is a very important factor that probably supports the scenario of seawater intrusion, where the faults network acts as the major pathway of groundwater salinity into the aquifer. The sweet-salty spring phenomenon that takes place at the north-western part of the study area demonstrates the complexity of underground water pathways. The situation is not investigated explicitly and it is not part of this research, but it is in the future plans of the authors.
Although the geophysical technique that was applied is considered as capable for deep investigations, it is necessary to combine it with other geoenvironmental data, such as geological/geomorphological.
The examples of Figure 4 show different curves and models corresponding to different contexts. Measurements 334 and 341 (Figure 4a,b) are totally different even though they were acquired very close to each other on the same formation (Plyllites/Quartzites nappe) in the seafront (Profiles AB & AG, Figures 1 and 10). Thus location 334 is interpreted to be affected by seawater intrusion below 40 m depth, while location 341 is not. Concerning measurement 171 (Figure 4c) that was acquired on Sediments near G-17 borehole (measurement is included in profile EF, Figure 5), fresh water can be identified (resistivity values near 20 Ω·m), while this results is confirmed by the geochemical data from G-17 borehole. Chemical spatial distribution maps (Figure 7), GQI map ( Figure 8) and Durov diagram (Figure 9) prove that the ground water of G-17 borehole is fresh and not affected. The depth of the water level can be recognized with some precision too. The freshwater is estimated at 80 m depth from TEM corresponds nicely with the level of the borehole G-17 (85 m) [26]. Very interesting results are extracted by Measurement 246 that was acquired on Plattenkalk nappe near boreholes G-05 and G-80 (measurement is included in profile AB, Figures 5 and 10). Here is it shown that freshwater ranges between 90 m and 120 m with resistivity values more than 10 Ω·m, while the fresh water-salt water interface is seen in 120m depth. The groundwater that is pumped from G-05 and G-80 boreholes is fresh but affected as resulted from Chemical spatial distribution maps (Figure 7), GQI map ( Figure 8) and Durov diagram ( Figure 9). Therefore, we can conclude that resistivities greater than 20 Ω·m result from the presence of fresh water in the subsurface and when resistivity is less than 20 Ω·m indicate the existence of salty/contaminated water (Table 2). For 2D visualization of the 1D models, 6 cross sections (AB, AG, CD, EF, FG and HI) were constructed using 1D interpolation ( Figure 5). The AB cross-section is approximately 9 km long, has NW-SE direction and it comprises the two major units (Plattenkalk and Phyllite/Quartzite nappes) covered by the Neogene and Quaternary sediments. Plattenkalk unit (presented in purple colour) is recognised at the south-eastern part of the profile (Figures 1 and 10) with resistivities greater than 1000 Ω·m (as expected). Phyllite/Quartzite nappe is defined very well in the north-western and central part of the profile with resistivities generally about 100 Ω·m and sometimes in places up to 300 Ω·m, when the nappe is rich in quartzite. Sediments which are composed by marly limestones, marine and alluvial deposits, clays, etc. have resistivities ranging from 100 to 500 Ω·m, and occasionally, when Neogene limestones are more cohesive and stiff, the resistivities are even higher up to 1000 Ω·m. Generally, the sediments are not recognised clearly over Phyllite/Quartzite nappe (Figure 10), since the resistivity values overlap ( Table 2). The thrust fault between Phyllite/Quartzite nappe and Plattenkalk nappe is observed in AB geoelectrical model, which its position and direction are verified by the geological cross-section. Additionally, the transition zone between them is recognised with resistivity values between the values of the two nappes ( Figure 10). The AB geoelectrical profile is in good correlation with the geological AB cross-section concerning the recognition of the nappes (especially Phyllite/Quartzite nappe with Plattenkalk nappe) and confirmation of faults ( Figure 10). In addition, the AB geoelectrical profile gives information about the permeability of the rocks (conclusively Phyllite/Quartzite nappe is not impermeable) and the aquifer characteristics (water table depth and salinity).
As already mentioned, Phyllite/Quartzite nappe is considered almost impermeable, except for faults that could be major potential conduits of water flow inland. The north-western part of AB 2D geoelectrical model ( Figure 10) shows the seawater encroachment and complements the seawater intrusion scenario. It illustrates the low resistivity values (blue colour) situated into the Phyllite/Quartzite nappes 25 m under the surface extending up to 1.5 km inland and it is considered as the strongest evidence for this scenario. Based on rock types, two different types of aquifer are expected. Sediments are expected to host a shallow aquifer (G-63 and G-81 are evidence- Figure 10) which lie above Phyllite/Quartzite unit. The problem is that resistivity values overlap ( Table 2) and it is hard to separate them, which is part of the limitation of geophysics [34]. The deeper aquifer is expected to be found in discontinuities (fracture zone in Phyllite/Quartzite nappe or karst in Plattenkalk unit).
AG geoelectrical profile about 7 km long was constructed along the coastline at the northernmost part of the study area ( Figure 1) presented in Figures 5 and 10. No geological section is appeared in comparison with AG geoelectrical profile as this section is on the same geological formation (Phyllite/Quartzite nappe). If this nappe was impermeable (according to [26], the expected geophysical section should give a homogenous geoelectrical model. Instead of this, low resistivity values (≈1 Ω·m) are illustrated (blue colour) situated of course into the Phyllite/Quartzite nappes 50 m under the surface extending up to 4.8 km to the East. From 4.8 up to 7.2 km (end of the AG profile), no seawater intrusion is observed (resistivity values 30-100 Ω·m) and the formation is considered impermeable. This is in agreement with the geological map ( Figure 1). The geophysical survey confirms the NW-SE inferred fault which crosses the Phyllite/Quartzite nappes ( Figure 1) and it is also marked in Figure 10 with dotted line.
Since AB and AG profiles are sub-parallel (Figure 1), we can conclude that the western part of the study area is governed by salinization that extends to the southeast as shown in profile AB. The fracture zone (1st left black colour ellipse) ( Figure 10) probably acts as a conduit for seawater encroachment into aquifer. The 2nd ellipse in the middle most likely represents the seawater intrusion through the thrust zone. The deeper aquifer (maybe a little influenced by seawater) is recognised in the dolomites with resistivity values greater than 10 Ω·m, while underneath there is brackish water (≈1 Ω·m) probably originating from seawater which intrudes through faults from the north-eastern part of the seashore (3rd ellipsis). Many boreholes through the Plattenkalk nappe with relatively high pumping rate confirm the existence and the amount of groundwater.
The CD, EF and FG geological cross-sections ( Figure 2) cross mainly above sediments, superimposed to Phyllite/Quartzite nappe, for which their thickness was approximately estimated by the available data (geological map and borehole logs). CD and EF are 4 km long, while FG cross-section is 7 km long having NW-SE, NNW-SSE and SSW-NNE directions respectively. Their 2D geoelectrical models (juxtaposed 1D inversion models) (Figure 4) demonstrate that the resistivity values of these sediments vary depending on the rock types from 100 to 500 Ω·m (Neogene limestone, clays, gravels, etc.); however when Neogene limestones are more cohesive and stiff, the resistivity values are even more increased up to 1000 Ω·m. Generally, due to resistivity values overlapping between Phyllite/Quartzite nappe and sediments, it is hard to distinguish them. Boreholes confirm that there is groundwater in sediments, which is confirmed by TEM data, where the recharge zone is defined with resistivity values 20-40 Ω·m. Low resistivity areas at deeper layers are evidence for fracture zones in the Phyllite/Quartzite nappe. These zones are conduits for seawater intrusion, and consequently contamination of the groundwater (low resistivity areas shown with blue colour).
Finally, the HI geoelectrical model is a small part (2 km long) of HI geological cross-section and crosses only the plattenkalk nappe. The high resistivity values indicate the dolomites and dolomitic limestones on the top (up to 100 m depth). Then the aquifer system is clearly identified and it is confirmed once more that TEM technique is very good for identifying conductive zones underneath resistive layers. The water table (80 m-100 m from the surface) and the interface between fresh and brackish waters (beginning from 110 m from surface) are depicted in section HI of Figure 5.
Four 3D geoelectrical resistivity layers were constructed in order to give some extra information about the surficial geology (rock types and general tectonic regime) and aquifer characteristics as well ( Figure 6). Therefore, the upper Layer A can be correlated with the surficial geological setting, thus a layer with the geology is overlain. At the eastern part, the high resistivity Plattenkalk nappe is shown, while sediments and Phyllite/Quartzite nappe have almost the same resistivity values. In that part of the study area there is no evidence of seawater intrusion, but in contrary the phreatic aquifer is shown in the central and southern part (Quaternary deposits and Neogene limestones, respectively).
As already mentioned, two aquifers are expected in the area (shallow and deep); where from log information the pumping level ranges mostly between 20 m and 100 m depth. At both B and C Layers plattenkalk nappe was recognised (high resistivity values) at the eastern part. In the north-western part (especially for Layer C) 50 m below surface the most reasonable interpretation is that seawater intrudes inland. Except for affected areas (where the resistivity measurements are not reliable), the same resistivity values exist for this nappe, as described earlier at 2D TEM interpretation which overlap with sediment resistivity values.
Layer D gives information about the size and geometry of salinity ( Figure 6). Specifically, the south-western part of the 3D model seems unaffected by seawater. The north-eastern part of the area seems not contaminated as well, but (as mentioned above) no TEM measurements are available, and maybe seawater encroaches and contaminates plattenkalk nappe at the eastern part. Finally, a graben with NW-SE direction crosses the study area and it is envisaged here to be a major conduit for seawater.
Moreover, as already studied in previous work [19] in order to confirm the origin of the high salinity, a preliminary geochemical study was performed, too. This approach is important because groundwater chemical results are ground truth data that validate the geophysical survey. GQI Quality map as already discussed showed great relationship between geophysics and geochemical analyses as well (Figure 8). This approach of groundwater quality map based on WHO guidelines, locate the northern-central part of the area near the seashore as contaminated, which is also in agreement with scenario 1. Moreover, apart from the visual comparison with the TEM 3D layers (Figure 6), a scatter plot of the GQI values with the resistivity data of Layer D is presented in Figure 8. The correlation coefficient R that occurs is 0.3848 which consider as weak-medium. Although the spatial distribution of the chemical sample location is not good (Figures 1 and 8), there is still a correlation between resistivity data (Layer D) and water analysis (GQI). Of course a detailed dataset is required that will further support any possible strong existing correlation and it is in the future plans of the authors.

Conclusions
This work emphasizes that the application of TEM technique is a valuable tool in groundwater resources assessment. The reliability of the TEM data, the portable nature of the equipment, the fast data gathering and the minimal participants needed for its operation (3 people maximum) make this electromagnetic technique an excellent tool for hydrogeology. Also TEM in correlation with chemical analyses as used within the frames of this work allows independent comparison. Although different possible explanations exist, regarding Geropotamos case, the evidence presented is best interpreted that the saline water is due to seawater intrusion. During the winter period precipitation (mainly snow and rain), falling mostly in the mountains (Talea Ori Mountain), either runs off surficially (creating Geropotamos River), or soaks into the ground as infiltration (replenishment of deep aquifer-Plattenkalk nappe, or phreatic aquifer-Neogene sediments). During summer season, seawater is here interpreted to intrude inland along faults with NE-SW direction, contaminating the freshwater.
The application of TEM technique provides a crude model of the aquifer (qualitative characteristics) since geophysical measurements are on multi-face materials (aquifer and hosted rocks). In complex environments, like Geropotamos area, the interpretation of geophysical results does not always lead to single interpretations and needs experience, due to resistivity values overlapping (Table 2). Nevertheless, the method is robust in that brackish waters can be easily distinguished by TEM, but the range of variation within brackish waters is not fully determinable. In addition, time-lapse TEM measurements are recommended as future work for monitoring the movement of the fresh/saline water interface over time To sum up, geophysical interpretation showed the most probable source of contamination, as well as the qualitative and quantitative characteristics of the groundwaters. So, the results of this work aid in search for solutions for water management in order to achieve and preserve a long-term protection of the available groundwater resources of the area.