Assessing Aquifer Salinization with Multiple Techniques along the Southern Caspian Sea Shore (Iran)

Abstract

This study focuses on the salinization of the coastal aquifer in the Mazandaran Province (Iran) within four different sites. Many factors can lead to declining groundwater quality, but this study focuses on the seawater intrusion area. Therefore, locating the interface between saltwater and freshwater is very important. For this purpose, three characterization methods with different accuracies have been employed: the Verruijt equation, vertical resistivity sounding, and an electromagnetic survey. Vertical resistivity sounding and the electromagnetic survey were performed near existing exploration boreholes and were used to determine the saltwater interface. The results showed that the Verruijt equation provides a reliable localization in two of the sites, but in the other two sites, the determined interface is lower than the observed data. The geoelectrical method showed acceptable results, but often this method cannot distinguish between the saltwater and saline aquitard boundary. The electromagnetic method showed a high accuracy in all the study sites and proved to be the most reliable method compared with the other techniques employed in this study. The results from this study are useful in helping to identify the most suitable technique for locating the freshwater/saltwater interface, especially in those sites where a detailed characterization via multilevel sampling is not feasible for technical or economic reasons.

1. Introduction

In coastal aquifers, because of the high density of marine saltwater, usually, gravity forces the distribution of saltwater below the freshwater lenses [1]. The boundary between saltwater and freshwater is not sharp because of the diffusion processes and mechanical dispersion, but it is rather a volume [2,3] that is called the transition zone [4].

Determining the depth and the extent of the interface between saltwater and freshwater in coastal aquifers is very important. In fact, in many studies, different methods have been used to locate this interface [5,6,7]. Badon-Ghyben [8] and Herzberg [9] were the first scientists that studied the interface between saltwater and freshwater with a simple one-dimensional analytical solution. After that study, Glover [10] provided a two-dimensional approach, which was subsequently improved [11,12,13]. Delineating the mixing zone with this last method is more challenging than using the original approach since it depends more on hydrogeological data collected. In many areas, because of the rough physiographical characteristics and the limited number of observation wells, the application of the two-dimensional approach is not straightforward [14].

The same applies to density-dependent numerical models, which are more flexible than analytical solutions but need even more data to be set up, calibrated, and validated [15,16,17,18].

Therefore, alternative methods (e.g., geophysical methods) can be more cost-effective, provided that their typical drawbacks are carefully addressed, such as the large uncertainties that may occur if the geological setting is not properly constrained by lithostratigraphic data [19]. The most widespread method used for determining the saltwater/freshwater interface is the electrical resistivity tomography. Because of its good performance, it was used in Belgium [20], in the Lekki Peninsula in Nigeria [21], and in Martil-Alila coastal aquifer in Morocco [14], although it showed to be affected by some limitations. For example, it is impossible to clearly distinguish between geological formations of similar resistivity such as saline sand and saline clay [22]. In recent years, the transient electromagnetic method (TEM) has been increasingly used [23], though it is still less frequently employed than electrical resistivity tomography methods. The accuracy of the TEM method has been confirmed in numerous studies [1,14,24].

The Southern Caspian Sea aquifer is a peculiar case with respect to the saltwater intrusion issue in coastal aquifers. In fact, the Southern Caspian Sea is characterized by salinities of approximately 13.0 ± 0.5 g/L [25] and seasonal salinity changes of less than ~0.3 g/L. The mean annual salinity increases from the surface to the bottom of the sea, reaching values of approximately 13.3 g/L [26]. Despite this low salinity boundary, the coastal aquifer system hosts saline groundwater due to the presence of relict connate seawater within the permeable marine-sandy and sandy-loam Holocene sediments of the semi-confined aquifers of the Amol–Ghaemshahr plain [27]. In addition, as far as the authors are aware, in this area, there are no studies that have used multi-level samplers to precisely locate the transition zone, while a few studies have used integrated depth sampling in long screen wells to characterize groundwater salinity and geochemistry [27,28]. Unfortunately, the integrated depth sampling technique is known to produce misleading results when used to interpret the origin of different saline sources [29] due to artificial mixing within the long screen wells. The latter is mainly ascribable to sampling procedures, but can also be due to the head and density gradient effects. Given the lack of a depth-dependent sampling dataset on salinities in this area, well-known geophysical and analytical methods were combined in this study to achieve a robust conceptual model of the saltwater wedge position in different sites along the Southern Caspian Sea. Though the methods here employed are not novel, the proposed approach could be applied to many other coastal aquifers as a cost-effective characterization technique, where a network of multi-level samplers is not always available, except for few well-equipped research sites [30,31].

2. Materials and Methods

2.1. Study Area

In the Mazandaran Province, the freshwater supply for agricultural and domestic needs is mostly covered by groundwater exploitation of the coastal aquifers, either because the water table is easily accessible, being less than 10 m below ground level (BGL) [28], or because the storage capacity is wide, as this alluvial aquifer can be up to 250 m thick. An average of 0.12 Mm³/km² of groundwater is pumped from this aquifer to supply the water demand for cultivation [28], providing also drinking water for 75% of the local population [27]. The drilling and operation of many wells in this area have provoked a progressive lowering of the water table in recent years [27]. Therefore, in the east of the Mazandaran Province, farmers often find saline waters during the midsummer season, which provoke many problems for agriculture and residents.

In addition to pumping, during the summer period, the aquifer also suffers from a lack of recharge from rainfall, surface watercourses, and the numerous coastal ponds that characterize this area. In fact, since the area has a humid subtropical climate (following the Köppen classification), the average annual rainfall ranges from 900 to 1300 mm but is unevenly distributed with a peak in late autumn and a minimum in early summer. The average monthly temperatures instead record a minimum (11 °C) during winter and a maximum in late summer (30 °C), so the average annual evaporation rate is quite elevated and ranges from 800 to 1100 mm [29]. As a consequence, the groundwater table is affected by seasonal fluctuation that on average amounts to 1.45 m, with minimum values (0.45 m) in the western part of the study area and maximum (up to 2.75 m) in the eastern part.

The Southern Caspian Sea area has a well-developed surface drainage network (Figure 1). The rivers are perennial since they are fed by the Elburz Mountains and their waters are used for agricultural purposes. From a hydrogeological point of view, this area consists of four hydrogeological units of unconsolidated sediments. The top one is a highly permeable unit consisting of sand and gravel calcareous materials. This unit constitutes the shallow unconfined aquifer and is limited to the upper alluvial plain. This unconfined aquifer gradually becomes thinner toward the north and disappears at a distance of 10 km from the Caspian Sea (Figure 2c). The second unit is an aquitard, composed of clay and silty-clay sediments, which impede the groundwater flux between the unconfined aquifer and the third hydrogeological unit. The latter is composed of gravel and fine- to medium-sized sands and forms the major aquifer of the area. The fourth hydrogeological unit is made of marine silty and sandy sediments saturated with paleo-brackish waters and behaves like a poorly transmissive hydrogeological unit with respect to the major aquifer of the area. In the study area (central panel of Figure 1), the aquifer is unconfined and formed by the third hydrogeological unit described above. The bedrock below the Quaternary sediments is generally found at a depth of more than 500 m [32].

With regard to the importance of the freshwater/saltwater interface location, different transects perpendicular to the coast were chosen according to the following principles: (i) different types of hydrogeological conditions had to be represented to account for lateral heterogeneities; (ii) a minimum of four monitoring wells had to be present along a flow line to allow for direct observations of the interface; (iii) site accessibility for field investigations and sampling had to be guaranteed.

Each of the four field sites (Nur, Mahmudabad, West-Babolsar, and East-Babolsar), covering an area of approximately 6 km² (3 km along the shoreline × 2 km inland), include 3 transects spaced by a distance of 1500 m. Along each transect, 15 points were selected and investigated with both electrical resistivity sounding and electromagnetic sounding (bottom panel of Figure 1).

2.2. Groundwater Heads and Quality Data

The groundwater heads and quality data (TDS and major ions concentration) in agricultural and domestic wells were provided by the Ministry of Energy and by the Mazandaran Regional Water Company. Groundwater heads were measured in static conditions using an electric contact phreatimeter and converted into equivalent freshwater heads (H_f) to account for the variable density and salinity [33]. Groundwater heads measured in 340 wells were interpolated for May 2015 and October 2015 via ordinary kriging with a linear variogram, an anisotropy ratio of 2, and the nugget effect. A blank file was used to trim away the Caspian Sea from the contour map. A vector map was built on the interpolated grid to show the groundwater flow directions.

Chemical analysis of the water samples was carried out at the Mazandaran Regional Water Company, following an analytical procedure suggested by the American Public Health Association [34]. Specifically, Cl⁻ was estimated using the standard AgNO₃ titration method, Na⁺ was measured via the standard F-AAS method, and TDS was measured gravimetrically after heating at 105 °C for 24 h. The TDS map was interpolated by ordinary kriging with a linear variogram, an anisotropy ratio of 2, and the nugget effect. A blank file was used to trim away the Caspian Sea and the mountain area, where the groundwater quality data were sparse. To assess the accuracy of the TDS interpolation, open borehole logs of specific conductance (electrical conductivity normalized at 25 °C) were used to detect the saltwater interface in the monitoring wells located along the transects at each site. The downhole specific conductance data were acquired using a Solinst Levelogger LTC connected to a nylon wire. The freshwater/saltwater interface was set at a specific conductance of 4000 μS/cm, according to the classification scheme used in the Abu Dhabi Emirate [35]. Given that open borehole logs of specific conductance can be affected by biases with respect to the multi-level samplers’ data, the relative uncertainty of the observed saltwater interface was estimated to be around ±5% of the measured specific conductance as reported for a similar alluvial coastal aquifer [36].

2.3. Analytical Methods

The Ghyben–Herzberg analytical solution was the first to address the problem of the saltwater wedge. After that, many other solutions were presented [10,11,12,13]. The Verruijt analytical solution, widely accepted and used in recent years, has been employed in this study and is provided below:

$$Z=-\sqrt{\left(\frac{q^2}{{\beta }^2{K}^2}\right)\left(\frac{1-\beta }{1+\beta }\right)+\frac{2qx}{\beta K\left(1+\beta \right)}}$$

(1)

$$H_f=\sqrt{\frac{2\beta q}{K\left(1+\beta \right)}x}$$

(2)

$$X_0=\frac{q}{2\beta K}.$$

(3)

In Equation (1), Z is the depth of the freshwater/saltwater interface below the water table (m); in Equation (2), H_f is the equivalent freshwater head (m); in Equation (3), X₀ is the position of the freshwater/saltwater interface at the ground surface (positive values are from the shoreline toward the sea and negative values are inland from the shoreline). The terms of the equations are defined as follows:

$$q=Ti$$

(4)

where q is the groundwater specific discharge (m²/d); T is the aquifer transmissivity (m²/d), and i is the hydraulic gradient (-).

$$\beta =\frac{{\rho }_s-{\rho }_f}{{\rho }_f}$$

(5)

where β is a coefficient (-) accounting for the density contrast between freshwater (ρ_f = 1.0 gr/cm³) and saltwater (ρ_s = 1.013 gr/cm³).

$$K=\frac{T}{b}$$

(6)

where K is the hydraulic conductivity (m/d), and b is the aquifer thickness (m).

x is the distance (m) from the shoreline inland.

High K values lead to a shallow interface, while low K values give a deep interface [5], since the interface’s depth is related to aquifer lithological characteristics [37]. The values of T and K for each transect in the four sites were taken from the pumping tests performed by the Mazandaran Regional Water Company (see Supplementary Information).

2.4. Vertical Electric Sounding

The geo-electrical methods may contribute to an improved understanding of aquifer characteristics such as thickness of alluvium, particle size, bedrock depth, and an estimated depth of the freshwater/saltwater interface. Vertical electric sounding (VES) was performed using the Schlumberger configuration with a maximum current electrode separation of 150 m along each transect. An ABEM SAS 1000 Terrameter was used to obtain the apparent resistivity at each VES point. The data were interpreted by the conventional partial curve matching technique with two-layer master curves, in conjunction with an auxiliary point diagram. The data obtained from the VES were controlled using the software RESIST; in this way, true resistivity values were obtained. Drilling exploration borehole data were used to interpret the information derived from the VES, for all 12 transects.

2.5. Transient Electromagnetic Sounding

In order to determine the position of the freshwater/saltwater interface, the transient electromagnetic (TEM) method was used in the four field sites, in addition to the Verruijt equation and the VES method. In each field site, the TEM vertical profiles were calibrated along the transects using stratigraphic logs. The TEM method has the capability of achieving waves reflected from the underground layers to the ground surface at a very low electromagnetic frequency. This method can be performed using natural ground signals (inactive, such as magnetotelluric signals) or with artificial transmitter (active) in the near square or in the far square, using powerful military transmitters in the very low frequency (VLF) radio spectrum. The Developers Water Saraveyeh Company, combining the magnetotelluric signals and the VLF, created a tool with sensitive receivers that can acquire the low telluric turbulent waves as the minimum and maximum difference phase. The standard configuration was applied using a single-turn, 40 × 40 m² transmitter loop with the receiver coil in the center of the loop. The transmitter of the TEM system had a maximum output of 3 A, resulting in a magnetic moment of 4800 m²A. The groundwater depth was determined using a brass rod. The electromagnetic square mobile receiver tool near this brass rod was able to determine the groundwater depth by taking the distance from the earth surface. The freshwater/saltwater interface was determined following the Fitterman and Stewart approach [38].

2.6. Methods’ Evaluation Using Statistical Operators

The reliability in capturing the position of the freshwater/saltwater interface with respect to the observed transition zone in the monitoring wells was assessed by calculating and comparing the value of the mean absolute error (MAE), the root mean square error (RMSE), the percent of bias (PBIAS), and the Wilcox index of agreement (d), for all the applied methods in the four study sites. PBIAS values were calculated as shown in [39]:

$$PBIAS=\frac{{{\displaystyle \sum }}_{j-1}^n\left(Ob{s}_j-Cal{c}_j\right)\times 100}{{{\displaystyle \sum }}_{j-1}^n\left(Ob{s}_j\right)}$$

(7)

where Calc_j is the calculated value corresponding to the observed one Obs_j. The optimal value of PBIAS is 0.0, with low-magnitude values indicating an accurate model simulation. Positive values indicate model underestimation bias and negative values indicate model overestimation bias [38]. d values were calculated as shown in [40]:

$$d=1-\frac{{{\displaystyle \sum }}_{j-1}^n{\left(Cal{c}_j-Ob{s}_j\right)}^2}{{{\displaystyle \sum }}_{j-1}^n{\left(\left|Cal{c}_j-Ob{s}_j\right|+\left|Ob{s}_j-Obs\right|\right)}^2}$$

(8)

where Obs is the mean value of the observed data. This statistical indicator can vary from 0 to 1; a computed value of 1 indicates a perfect agreement between the measured and predicted values, and 0 indicates no agreement at all [40].

3. Results and Discussion

3.1. Groundwater Flow Directions and Quality

From the spatial analysis of the H_f piezometric maps of May 2015 and October 2015 (Figure 2), a recharge area located in correspondence of the Elburz Mountain is visible in the southern part of the study area. In particular, in the western portion of the Amol–Ghaemshahr plain, the recharge area is not far away from the Caspian Sea and the hydraulic gradient is high, leading to a considerable groundwater discharge toward the sea; moving eastward, this pattern becomes less pronounced. The main rivers of the area supply water to the aquifer, since the bottom of the riverbeds are located always above the water table and the piezometric contours diverge from the main river courses. By comparison of the two maps (May versus October 2015), it is also evident that, at the end of the rainy season (the beginning of May), the groundwater table is higher than it is at the end of the summer season (the beginning of October), being +0.80 m on average with peaks in the mountains of +1.55 m. Given this transient behavior of the groundwater heads, the most exposed period to a possible salinization of groundwater resources is likely to be late summer–early autumn. For this reason, it was decided to perform the sampling campaign to record the TDS and the major ions in October 2015.

Before using TDS as a proxy of aquifer salinization, a scatter diagram correlating TDS versus Cl⁻ and Na⁺ was plotted to demonstrate that these ions are governing the TDS trend, which is only marginally affected by other species (such as nitrates and sulfates), which are attributable to human activities. The high R² (0.91 for Cl⁻ and 0.88 for Na⁺) confirmed this hypothesis (Figure 3).

Figure 4 shows that most of the wells located in the Western portion of the study area are characterized by very fresh (0 < TDS < 1000 mg/L) and fresh waters (1000 < TDS < 2000 mg/L), while low brackish (2000 < TDS < 4000 mg/L) and medium brackish (4000 < TDS < 6000 mg/L) groundwater become more frequent in the eastern portion of the study area [34].

As reported by other authors [27,41], the aquifer system of the Amol–Ghaemshahr plain is characterized, especially in its eastern portion, by residual connate saline waters due to the eustatic variations in the Caspian Sea over the last 10,000 years [42]. Holding this, while Sites A, B, and C may essentially be affected by actual sweater intrusion, Site D may also be affected by autonomous salinization processes, exacerbated by the upconing of saline groundwater due to aquifer overexploitation. In general, defining a correct conceptual model of coastal aquifer salinization is vital not only for the sustainable water management but also for the implementation of saltwater intrusion mitigation strategies via coastal engineering systems [19]. When a dedicated network of multi-level samplers is not available, geophysical techniques can help to constrain the problem, as shown in the following sections.

3.2. Comparison of the Analytical Solution, VES, and TEM Results

In order to confirm the existence of groundwater discharge toward the sea, as suggested by the spatial analysis of the H_f piezometric maps of May 2015 and October 2015 (Figure 2), the location of the freshwater/saltwater interface at the ground surface (X₀) was calculated by means of the Verruijt analytical solution in each transect (

Table 1. The distance of the freshwater flux from the coast toward the sea (X₀) in each section of the four field sites, for May 2015 and October 2015.

Site Name	Section	X0 (May 2015) m from the Coast	X0 (October 2015) m from the Coast
Site A (Nur)	A1	19.84	19.05
	A2	6.73	6.56
	A3	18.77	17.65
Site B (Mahmudabad)	B1	8.85	8.69
	B2	12.37	12.14
	B3	8.42	8.18
Site C (West-Babolsar)	C1	0.15	0.13
	C2	1.93	1.72
	C3	0.94	0.81
Site D (East-Babolsar)	D1	2.95	2.15
	D2	1.04	0.97
	D3	0.95	0.69

). Positive values stand for positions from the coast toward the sea, suggesting that a freshwater discharge is taking place, and negative values correspond to positions inland from the coast, suggesting that saltwater intrusion is occurring. Tidal fluctuations can affect the location of the freshwater/saltwater interface [19]; however, given that in the Southern Caspian Sea the maximum tidal range is just 21 cm with an average value of 10 cm [43], this effect was considered negligible.

All the calculated X₀ values are positive, meaning that there is always a freshwater discharge toward the sea in the two available monitoring campaigns of May 2015 and October 2015 (Figure 2). In addition, X₀ is always lower in October 2015, witnessing a decrease of freshwater discharge at the end of the summer months, as postulated by the analysis of the piezometric maps.

Accordingly, in the western transects (Nur and Mahmudabad sites), X₀ is on average 12.05 m while in the eastern transects (West and East Babolsar) X₀ is on average 1.08 m, testifying that, in the western part of the Amol–Ghaemshahr plain, the groundwater discharge toward the sea is higher (Figure 2).

The Verruijt analytical solution was also used to calculate the freshwater/saltwater interface for each transect of the four field sites. The analytical solution (AS) results were then compared with the VES and TEM interpolations, as well as with the monitoring campaign records of October 2015.

The average interface’s depths calculated with the AS were −215 m above sea level (ASL) in Nur (Site A in Figure 5), −178 m ASL in Mahmudabad (Site B in Figure 6), −104 m ASL in West-Babolsar (Site C in Figure 7), and −84 m ASL in East-Babolsar (Site D in Figure 8).

VES data were interpreted with the help of the boreholes’ log data available for each field site, then the freshwater/saltwater interface’s depth was determined in each transect. The average interface depths in Sites A, B, C and D were −119, −128, −91 and −65 m ASL, respectively (Figure 5, Figure 6, Figure 7 and Figure 8).

The same was done with the TEM data, and the average depths of the freshwater/saltwater interface in Sites A, B, C and D were calculated as −111, −108, −75 and −60 m ASL, respectively (Figure 5, Figure 6, Figure 7 and Figure 8).

More specifically, Figure 5 shows that, in Site A, the best method for locating the freshwater/saltwater interface is TEM, while the AS does not perform well. Likewise, the VES results are poor in A1 and A3, where the location of the interface with respect to the observed freshwater/saltwater interface in the monitoring wells is grossly overestimated in proximity to the shoreline and underestimated inland.

In Site B, the best method is still TEM (Figure 6), while AS provides satisfactory results only for Transect B1, overestimating the location of the freshwater/saltwater interface’s depth in the other two transects. The VES method seems to be less reliable, producing both overestimation and underestimation if compared with the records in the monitoring wells.

In Site C, the best method to capture the elevation of the transition zone is again TEM (Figure 7), while AS provides acceptable results only for Transect C1. The VES method is also not fully reliable, producing unrealistic results, above all moving inland.

When describing the performance of the different methods in Site D, it is worth pointing out that AS cannot deal with the connate saline groundwater since it assumes simplified boundary conditions such as freshwater from inland and saltwater from the sea. For this reason, AS totally failed to correctly locate the saltwater/freshwater interface’s depth in Site D, where connate saline groundwater is most probably present even closer to the coast. The TEM is graphically the best method, although the VES in Site D slightly outperformed the other sites, especially in Transects D2 and D3.

In general, the results of the three methods applied in this study showed the following: (i) The AS method produced fair results in transects where the aquifer is characterized by a considerable thickness. In fact, the closer the saline aquitard is to the ground level (Sites C and D), the less accurate the interpolation. Moreover, the presence of sources of salinity other than actual seawater intrusion generates a high disturbance in the AS accuracy; (ii) In some transects, the VES method seems to be greatly influenced by the presence of the saline aquitard, often reproducing the saline aquitard/aquifer contact at a higher elevation than the actual one, rather than reproducing the saltwater/freshwater interface’s depth, as already raised by other studies [14,22]; (iii) The TEM method has proved to be the most reliable in capturing the freshwater/saltwater interface in all the configurations encountered in the transects. The values of the statistical indicators used to quantify the efficiency of each method in capturing the freshwater/saltwater interface for all the four field sites (A, B, C and D) are reported in

Table 2. Statistical indicators quantifying the efficiency of the three methods AS, VES and TEM in locating the position of the freshwater/saltwater interface, with respect to the elevation of the transition zone measured within the observation wells.

Site Name	Method	MAE (m)	RMSE (m2)	PBIAS (m)	d (-)
Site A (Nur)	AS	47.4	56.8	−36.2	0.49
	VES	30.7	35.6	−12.2	0.72
	TEM	8.4	10.2	−4.4	0.98
Site B (Mahmudabad)	AS	37.3	44.2	−33.1	0.67
	VES	26.7	31.6	−13.8	0.89
	TEM	10.5	12.3	+1.6	0.96
Site C (West-Babolsar)	AS	42.1	57.6	−35.5	0.49
	VES	27.1	41.2	−12.6	0.52
	TEM	5.1	6.7	+0.1	0.97
Site D (East-Babolsar)	AS	58.6	83.9	−96.1	0.13
	VES	5.9	10.5	−7.1	0.96
	TEM	1.6	2.2	−0.5	0.99

.

In general, the AS method displays a very poor efficiency, with an MAE ranging between 37.3 and 58.6 m, and a d between 0.13 and 0.67, depending on the field site. This is due to local heterogeneities of T and K that influence the inland saltwater wedge extension. In fact, the use of K values higher than the mean ones derived from the pumping tests (see Supplementary Information) would considerably increase the saltwater/freshwater interface’s depth. However, to improve the efficiency of the AS method, a K increase of an order of magnitude would be needed. This is in contrast with the observed K variability that is comprised between 1.6 to 5.6 m/d in all sites (see Supplementary Information). Clearly, the choice of the K value is a key factor in the AS method, but the information on the aquifer properties at each site demonstrates that the AS method underestimates the extent of the saltwater wedge at all investigated sites. The VES method is also affected by large errors (5.9 < MAE < 30.7 m and 0.52 < d < 0.96) and cannot be safely employed to determine the freshwater/saltwater interface in complex geological settings, even if it becomes reliable under particular conditions. For instance, VES performance improves when a strong contrast between freshwater and saltwater concentrations exists, such as in Site D, where connate saline groundwater is most likely present. The TEM method outperforms the other two methods in all field sites (1.6 < MAE < 10.5 m and 0.96 < d < 0.99).

Interesting information is given by the PBIAS values that are nearly always negative, indicating a constant overestimation of the freshwater/saltwater interface’s depth, with all the three methods and in every field site. This could be due to a number of causes: local heterogeneities of the aquifer properties (e.g., fine-grained lenses, mineralogical composition of the sedimentary formations, and the degree of compaction), groundwater temperature variation with depth, which may be due to the cold groundwater flowing from the Elburz Mountain into the deepest portion of the aquifer, and the existence of upward fluxes from the deepest part of the aquifer toward the ground surface.

If Table 2 is analyzed with a greater focus on the performance of statistical indicators at the various sites rather than on the methods employed, it can be seen that, in Nur (Site A), the location of the freshwater/saltwater interface’s depth is highly uncertain, while in East-Babolsar (Site D) it is reliably identified. This behavior may be due to both the different hydrological boundary conditions characterizing the two sites and their contrasting salinity levels. In fact, while the Nur site is affected by a large seasonal variation in the groundwater flux, the boundary conditions in the East-Babolsar site are more stable. In addition, in this latter site, the presence of higher salinities may have enhanced the efficiency of the geophysical methods in distinguishing among freshwater and saltwater and thus in locating the interface that separates them.

4. Conclusions

Various methods exist for determining the position of the freshwater/saltwater interference, and each of them has different limitations that must be assessed and estimated. In this study, the Verruijt analytical solution (AS), the vertical electric soundings (VES), and the transient electromagnetic method (TEM) were employed, and their results were compared with the actual position of the transition zone retrieved within the monitoring wells located in a complex geological setting—namely, the one characterizing the Mazandaran Province (Southern Caspian Sea, Iran).

From a methodological point of view, it was demonstrated that the AS method is affected by a low accuracy in reconstructing the interface’s depth, which is in general grossly overestimated, with an increasing bias moving inland. The VES method is strongly affected by the saline aquitard/aquifer boundary’s depth. The TEM method showed the best fit to the observed data. Although the TEM method most suitably and cost-effectively locates the freshwater/saltwater interface in this area, the TEM method should not be employed alone; rather, it should be used in conjunction with other hydro-geophysical methods so as to better cope with subsurface uncertainties and thus ensure a high degree of reliability of the results [44].

From a water management perspective, this study emphasizes that the freshwater/saltwater interface in the area moves progressively closer to the ground level, moving from West to East within the coastal zone of the Amol–Ghaemshahr plain (approximately from −92 to −37 m BGL). This is due either to a large fresh groundwater discharge toward the Caspian Sea in the eastern zone, induced by the proximity to the coast of the Elburz Mountain’s recharge area, or to the upconing of the paleo brines in the western zone [27,41], which adds its effect to the actual seawater intrusion from the Caspian Sea. These natural phenomena are then exacerbated by the excessive withdrawal of groundwater from a dense network of wells used for domestic purposes and for irrigation, which in turn may worsen the picture by adding a polluted returned flow to the system. Given the fragile balance of the water resources in this area, and in view of the possible changes due to the fluctuation of the Caspian Sea level or to a variation of the recharge rates from the Elburz Mountains, it follows that a good management of the same cannot leave aside a detailed monitoring of the freshwater/saltwater interface’s depth. Accordingly, this study contributes in a broader sense to the characterization and monitoring of saltwater intrusion in coastal aquifers, where a network of multi-level samplers is not yet available, providing, for example, insightful information on where to install the monitoring wells, because a sound coastal management for the future will strongly depend on a sound knowledge of the actual salinization status.