Dynamics of Saltwater Intrusion in a Heterogeneous Coastal Environment: Experimental, DC Resistivity, and Numerical Modeling Approaches

Abstract

Saltwater intrusion (SWI) is a critical concern affecting coastal groundwater sources due to natural and anthropogenic activities. The health of coastal aquifers is deteriorated by excessive SWI, mainly caused by the disturbance of the freshwater–saltwater equilibrium due to the escalating population, climate change, and the rising demand for freshwater resources for human activities. Therefore, gaining insight into the dynamics of SWI is crucial, particularly concerning the various factors that influence the intrusion mechanism. The present study focuses on the experimental simulation of saltwater in freshwater aquifers, considering boundary conditions and density-dependent effects. Two geological scenarios within coastal environments were investigated: First, a uniform, homogeneous case consisting of only sand, and second, a heterogeneous case in which layers of sand, clay, and sand mixed with pebbles are used. During the experiment, DC resistivity sounding data, as part of a widely recognized geophysical method, were collected and subsequently inverted to determine the depth of the freshwater–saltwater interface (FSWI). A finite element analysis was employed to generate numerical models based on experimental feedback. Further, for validation purposes, electrical resistivity tomography (ERT) data were collected from two distinct locations: near the seacoast and an aquaculture area. The ERT results show the presence of salinity intrusion in the study area, attributed mainly to groundwater overpumping and fish farming practices. The experimental findings indicate that the advancement of saltwater is affected by the geological properties of the media they traverse. The porosity (ϕ) and permeability (k) of the geological layer play a crucial role during the passage of saltwater flux into freshwater aquifers. The FSWI deviated along the clay boundary and hindered the easy passage of saltwater into surrounding layers. The alignment of experimental, numerical, and geophysical data suggests that this integrated approach could be valuable for studying SWI and can be applied in different geological settings, including tidal flats and alluvial plains.

1. Introduction

In recent years, an exponential increase in population has created enormous pressure on coastal groundwater resources. According to the high end forecast scenario, the global population residing in low-elevation coastal zones (LECZs) could increase by more than 50% between the base year of 2000 and 2030 [1]. The exponential growth of coastal populations imposes a significant strain on groundwater resources. Coupled with this, the elevated living standards prevalent in these hotspot areas exacerbate groundwater depletion. The extensive extraction of groundwater intensifies the threat of saltwater intrusion, particularly in coastal regions, emerging as a critical environmental concern. Researchers worldwide have reported salinity problems in their respective areas [2,3,4,5,6,7]. Anthropogenic and catastrophic events can also affect the dynamics of coastal aquifers; therefore, research must focus on analyzing and understanding the behavior of saltwater intrusion along coastal margins. Numerous efforts have been made to understand the mechanism of saltwater intrusion, with Henry’s problem emerging as the benchmark for density-dependent SWI models [8]. Various experimental data sets were used to validate the analytical solution for groundwater flow, thereby establishing correlations with the freshwater–saltwater mixing conditions. Many simulation models have been created to address various aspects, including the impact of water levels, tidal effects, seawater intrusion concentrations, the freshwater–saltwater interface’s migration rate, and the three-dimensional variable-density advection–dispersion model. Several laboratory experiments have been conducted to replicate the behavior of seawater mixing under controlled conditions [9,10,11,12]. Most of these studies concentrate solely on the progression of the saltwater wedge as a result of salinity contour distribution through a homogeneous geology. Thus, it is advantageous to trace the flow path of saltwater in various geological settings, utilizing geophysical methods that have been proven to be the most effective tools. Among these techniques, DC resistivity is particularly efficient in delineating saline boundaries in coastal aquifer environments [13,14,15,16,17]. Since the migration of saline water depends upon the geological characteristics of the host region, resistivity methods can help characterize various geological layers, such as clay, sand, silt, and shale. These layers exhibit different responses in terms of porosity and permeability based on their composition. Therefore, it is crucial to map the depth and thickness of different subsurface layers accurately. Resistivity methods have been previously used to characterize the subsurface geology for mapping saline zones [18,19,20].

The primary focus of laboratory simulations is to replicate the conditions under which salt diffusion arises due to the concentration disparity between fresh and saline water. However, flow velocity is also crucial for the advancement of the saline contours [21]. Several attributes have been incorporated into these experiments to address the dynamics of the freshwater–saltwater interaction, including an uncertainty analysis in fractured aquifers [22], saltwater up-coning [23], the pumping effect [24,25], the impact of beach face slope variation on SWI [26], subsurface dams to protect aquifers from SWI [27], the effect of an inclined boundary on SWI [28], and the impact of cutoff walls [29]. These experiments typically focus on a homogeneous geological environment, often represented by a sand-filled sandbox. However, the Earth consists of anisotropic behavior corresponding to various geological layers. Until simulation conditions closely resemble actual ground conditions, it becomes challenging to comprehend the interaction of salinity intrusion in complex geological environments. Therefore, accounting for the heterogeneity factor when conducting these laboratory experiments is crucial. There is always potential for research into the migration of the freshwater–saltwater interface during SWI experiments. This interface migrates when its natural state is disturbed by any external factor, and the extent of its intrusion depends on the rate at which saltwater flux intrudes. External factors such as groundwater overpumping result in deeper saltwater penetration into freshwater aquifers. Additionally, the migration of the interface depends on the geological heterogeneity of the media. The primary focus of the present study is to understand this migration behavior of the FSWI, particularly concerning the geological heterogeneity of the media.

The current study explores the boundary between freshwater and saltwater through a laboratory experiment using DC resistivity and numerical simulations. Initially, we conducted a sandbox experiment that involved two distinct scenarios: one characterized by a homogeneous setup with sand as the background material and the other featuring a heterogeneous environment incorporating sand, clay, and pebbles. Once a stable state was reached at the interface between freshwater and saltwater, DC resistivity measurements are obtained along the center of the profile (in a vertical cross-section). A forward model was generated to assess the sensitivity of the array utilized for resistivity data collection, followed by the creation of a one-dimensional inversion of the sounding model, which provides information about the depth of the interface. High-resolution imagery was captured throughout the experiment, and numerical solutions for both scenarios were obtained using the initial Henry parameters. Our main goal was to investigate the behavior of advancing saltwater wedges in various geological settings using a numerical model and the direct current resistivity method, focusing on the influence of clay layers in coastal environments. The present study emphasizes the effect of geological heterogeneity on the mechanisms of saltwater intrusion. In earlier research, uniform background materials like sand or artificial silica beads have typically been used to study intrusion behavior. For example, a homogeneous sandbox has been utilized in a study on subsurface dams to control SWI [30]. Few researchers have considered heterogeneity in their experiments by using artificial beads of different sizes [31,32]. However, the most significant aspect of this experimental study is the use of actual heterogeneous layers of sand, clay, and pebbles rather than artificial silica beads, making the study more realistic and applicable in terms of geological nature. Additionally, the study incorporates real-time geophysical data acquisition using DC resistivity to determine the depth of the freshwater–saltwater interface, providing a better understanding of subsurface contaminant flow. Considering the above parameters, this study suggests a more realistic observation and behavior of saltwater intrusion when it encounters different layers.

The present study holds significant utility as it offers valuable insights into the behavior of geological layers to SWI. The experimental findings were validated by an electrical resistivity tomography (ERT) investigation in the coastal area of West Bengal, India. The ERT study observed that anthropogenic activities, such as groundwater extraction for paddy crop cultivation and aquaculture practices, significantly impact groundwater quality. Consequently, it is essential to implement effective management strategies and policies, such as sustainable pond practices for aquaculture and crop rotation, to reduce dependency on water-intensive paddy crops, ensuring the sustainable use of groundwater in these regions.

Our study employs an integrated method that combines experimental data, numerical modelling, and geophysical approaches. This comprehensive approach provides valuable insights that can directly inform policymaking and the development of sustainable practices for managing coastal water. This research establishes a scientific foundation for developing targeted groundwater extraction regulations, implementing efficient land-use planning and designing physical barriers to protect freshwater supplies by demonstrating how various geological conditions influence saltwater intrusion dynamics.

2. Materials and Methods

2.1. Experimental Setup

The study conducted an experimental examination of saltwater intrusion, considering a conceptual model of a coastal scenario. This model accounts for the potential contribution of numerous human activities to the incursion of salinity in coastal aquifers (Figure 1). The experiment was carried out in a flow tank made up of acrylic glass with dimensions of 50 cm (length) × 30 cm (height) × 5 cm (width). The cross-section represents various coastal environments (Figure 2). The flow tank was divided into three chambers; the left chamber represents the freshwater reservoir, the middle chamber is filled with porous background material (as per two different scenarios), and the right chamber is the saltwater reservoir. Both freshwater and saltwater chambers are kept at a constant head for a continuous flow supply maintained by the outflow valves. Porous media are separated in these chambers using a US #18 fine mesh screen at a distance of 5.5 cm. Instead of artificial silica beads as a porous material, we used natural sand grains with a diameter ranging from 0.6 mm to 1.2 mm. We kept a head difference (∆h = 1 cm) between freshwater and saltwater reservoirs.

The saltwater was prepared in a 50-litre barrel using commercial salt. To maintain the level of salinity (close to that of seawater), we dissolved 35 g of salt in 1 L of tap water. Saltwater was dyed with a carmine color in a 1 g/L solution concentration to distinguish it from freshwater. Instead of ordinary dye for mixing with saltwater, we used carmine color, as it does not show any adsorption effect on aquifer medium, and carmine can migrate at the same rate as cl^- ions [34]. The saltwater density was maintained at 1.025 kg/L and measured with a WKM hydrometer. The grid marking was performed at the base of the flow tank to monitor the saltwater flow. Two models were adopted for experimentation to simulate natural subsurface conditions. Model A is homogenous, in which sand (fine-to-medium grain) was used as a porous filled material, and Model B is heterogeneous, in which different layers of sand, clay, and sand are mixed with pebbles at various depths. Sand grains were delicately compressed to avoid air-filled voids and ensure the homogeneity of the media. For experimental purposes, natural sand was used instead of glass beads or silicon balls to improve the accuracy of replicating the geology of the real-field aquifer. A Terrascience instrument acquired the DC resistivity measurement using an external 90 V battery supply. Thin stainless-steel electrodes were used as two current electrodes (C1, C2) and two potential electrodes (P1, P2) in a dipole–dipole array. A high-resolution DSLR camera was used throughout the experiment to monitor the saltwater movement.

2.2. Experimental Procedures

The experimental procedure resembles that of prior research, although most previous studies only focused on homogeneous cases. This study conducted experiments for two scenarios to achieve optimal responses under controlled conditions. The middle chamber of the flow tank was filled with porous material (sand) (homogeneous case) and clay, sand, and pebbles (heterogeneous case) in multiple horizontal layers, as shown in Figure 3. Before the actual saltwater experiment started, the system was set to allow freshwater flow from left to right (to the saltwater chamber) under gradient conditions (∆h = 1 cm). Excess freshwater was allowed through the valves fixed at different heights.

After achieving a steady equilibrium of freshwater flow from left to right, the intrusion experiment was initiated by opening the valve from constant saltwater head tank B. It was observed that saltwater rapidly flushed out the freshwater from the right chamber and began to invade the porous medium. Due to concentration differences, saltwater slowly migrated towards the freshwater chamber. As it was a density-dependent progression, the complete experiment was recorded, and high-resolution time-lapse images were captured to delineate freshwater and saltwater. During the investigation, considerable mixing of freshwater–saltwater flow was observed until the system achieved a steady-state condition. After reaching a stable state, there were no further observational changes in the location of the saltwater wedge.

2.3. Numerical Modeling

In the realm of mathematical formulations regarding the dynamics of coastal salinity, the Henry problem stands out as a widely recognized and accepted tool by numerous researchers [8]. This problem aims to streamline the experimental behavior by focusing on a vertical SWI near coastal aquifers. In this context, a balance is maintained between the inland flow of freshwater and the intrusion of seawater from the coast, until it is disturbed by an external factor. The aquifer is presumed to exhibit homogeneity and isotropy. The SWI is influenced by various factors, including anisotropy [35], heterogeneity [36,37], hydraulic conductivity [38], dispersivity [39], and the impact of the inland boundary conditions on SWI [40]. The Henry problem (HP)’s formulation is based on the concept of density-dependent flow, which involves the integration of variable-density flow equations, an advection–dispersion equation, and mixture density as a function of saltwater concentration. The flow system is governed by Darcy’s law in the following manner:

$$v=-K(\nabla H+\frac{{\rho }_{mix}-{\rho }_f}{{\rho }_f}{g}_z)$$

(1)

where $v$ is Darcy’s velocity (m/s), K is the hydraulic conductivity tensor (m/s), H is the hydraulic head (m), ${\rho }_{mix}$ is the mixed fluid density (kg/m³) which depends on concentration c, ${\rho }_f$ is the freshwater density (kg/m³), and $g_z$ is the unit vector corresponding to the direction of gravity.

$${\rho }_{mix}={\rho }_f(1+\alpha \frac{c}{c_s})$$

(2)

where $c_s$ is the concentration of seawater and $\alpha =({\rho }_s-{\rho }_f)/{\rho }_f$

The medium’s pressure gradient, fluid viscosity, and porosity structure influence the Darcy velocity field. The simulation modules provide the facility to solve Darcy's law, where diluted species are governed by diffusion and convection processes. The species are assumed to be diluted in such a way that density and viscosity are consistent in the mixture. The mass balance equation used for such a system is as follows:

$$V.\nabla c=\nabla .(D_F\nabla c)$$

(3)

where V is the flow velocity (m/s) obtained using Darcy’s law, c is the species concentration (mol/m³), and $D_F$ is the diffusion-dispersion tensor (m²/s).

The initial model parameters used in the present study for generating the SWI flow model are given in

Table 1. Summary of initial simulation model parameters used for this study.

Input Parameters	Value
Porosity, ϕ	0.35
$\mathrm{Freshwater}\mathrm{density}(\mathrm{kg}/{\mathrm{m}}^3){\rho }_f$	1000
$\mathrm{Saltwater}\mathrm{density}(\mathrm{kg}/{\mathrm{m}}^3),{\rho }_s$	1025
Saltwater concentration (mol/m3), c	1
$\mathrm{Inflow}\mathrm{velocity}(\mathrm{m}/\mathrm{s}),V_{Inflow}$	3.3 × 10−5
Pressure (Pa), P	0
Permeability (m2), k	1.02 × 10−9
Fluid diffusion coefficient (m2/s), D	1.886 × 10−6
Dynamic viscosity (kg/m.s), µ	0.001
Gravity (m/s2), g	9.8
Model Dimension
Length, x (m)	0.39
Height, y (m)	0.27

Note: For each input parameter, symbols are included.

and Equations (1)–(3) were implemented for numerical simulation are inspired by simple aquifer conditions [41]. The “subsurface flow” module environment of COMSOL version 4.4 was used, with boundary conditions set to no flow for both upper and lower boundary faces. The numerical solutions for SWI problems have previously been compared with semianalytical solutions and are well documented [42]. For a better approximation of the solution, the system was discretized and prepared to solve the problem using a finer mesh, as depicted in Figure 4. Finite element analysis (FEA) in COMSOL involves the use of the finite element method (FEM) to solve such fluid flow problems. FEM and appropriate mesh selection can effectively be used in simulations and analyses with high accuracy and computational efficiency. The numerical modeling approach offers a more comprehensive understanding of subsurface density-dependent flow conditions, particularly in addressing challenges such as salinity intrusion in coastal regions [43,44,45,46,47,48,49,50,51,52,53].

2.4. DC Resistivity Sounding

The resistivity data were collected once the freshwater–saltwater interface reached a steady state. The direct current (DC) resistivity method was used to determine the depth of the freshwater–saltwater interface. This method works on principle so that two current electrodes (C1, C2) are used to inject current (I) into the subsurface, and two potential electrodes (P1, P2) measure the potential difference (∆V) generated due to the interaction of current lines with different geological layers. The ratio of potential difference and current provides the resistance value (R), which is multiplied by the geometric factor (G) to calculate the apparent resistivity (${\rho }_{\mathrm{a}}$) value.

The equation used to calculate apparent resistivity (${\rho }_{\mathrm{a}}$) is as follows:

$${\rho }_{\mathrm{a}}=\mathrm{G}\left(\frac{\Delta \mathrm{V}}{\mathrm{I}}\right)$$

(4)

where the geometrical factor (G) is the linear arrangement of electrodes for the dipole–dipole array; G = πn (n + 1) (n + 2) a, where n is the dipole separation factor that varies from n = 1, 2, 3…; and a is electrode spacing.

Under the DC resistivity method, the Schlumberger array is conventionally used for vertical electrical sounding (VES) purposes to achieve a greater depth of investigation. However, it is less sensitive in identifying inclined subsurface bodies. Therefore, we have used a dipole–dipole array to investigate the saltwater wedge in this experiment (Figure 5). This array has greater sensitivity to detect the lateral resistivity variation and can detect vertical/inclined subsurface features with greater accuracy. The inclined features are progressive salinity contours associated with the freshwater–saltwater interface. The data were collected by readings in both forward and reverse modes to enhance the accuracy of resistivity measurements. After obtaining the average of these readings, the measured data were obtained. The current signal strength was also improved during the experiment by connecting the DC power supply batteries in series. This adjustment was necessary because the current was significantly attenuated due to the highly conductive nature of the host medium (water-saturated). A forward model response was initially generated using Res2dmod ver. 3.03 (Geotomosoft Solutions, Malaysia). For the dipole–dipole array, data points were measured along the center of the profile (vertical cross-section). The observed resistivity data were processed, and smooth models were obtained using Occam’s inversion method [54,55].

Some researchers have combined VES with other methods, such as time-domain electromagnetic (TDEM) techniques, to enhance our understanding of salinity intrusion in coastal aquifers [56,57]. In the present study, the multielectrode ERT field data collected for validation purposes using ABEM Terrameter with a Wenner array (Location 1) and Wenner–Schlumberger array (Location 2). After data acquisition, processing and interpretation were conducted utilizing Res2Dinv ver. 3.71 (Geotomo Software, Malaysia). The field data were inverted using a smoothness-constrained least-squares method. The choice of a suitable inversion scheme is crucial in obtaining high-resolution inverted subsurface images [58,59]. The forward model successfully detected a lateral change in resistivity distribution, enabling it to identify the inclined interface between freshwater and saltwater. The forward response generated is shown as apparent resistivity pseudosection (X vs. Ps.Z), which can later be used in an inversion engine to obtain an actual subsurface model (Figure 6).

3. Results and Discussion

3.1. Qualitative Observations

The saltwater intrusion dynamic was examined for both Model A and B; the comparative movement of the saltwater wedge is shown in Figure 7. This approach was designed to enhance our understanding of the flow pattern and solute transport associated with the system. It was generated using a head difference of ∆h = 1 cm between the left (freshwater) and right (saltwater) sides. When the experiment began, the time window was initiated and continuously recorded as the saltwater wedge (SW) advanced throughout the experiment. For Model A, the SW flow was smooth, and after 5 min, the initial height (y) of the SW was measured at 0.05 m, while its lateral extent (x) was found to be 0.4 m, as shown in Figure 7a. The FSWI SS-1 crossed the centre of the experimental box after 30 min, with a height of 0.23 m and a lateral extent of 0.14 m, as shown in Figure 7b. Model B began with a smooth progression of the SW, although minor shape alterations were attributed to the sand and pebble layer at the bottom of the container, as shown in Figure 7c. When the SW reached a height of 0.13 m and encountered a clay layer, the interface was disturbed, and the SW rate became slow compared to that of Model A. SS-2 developed after 40 min, showing a deviation from the expected smooth pattern.

The deviation was likely triggered by a tightly packed clay layer with a low porosity exhibiting a nearly impermeable behaviour. A small interface developed on the upper-right side of the box, as shown in Figure 7d, possibly due to the slightly porous clay layer that allowed its development. The comparative height and lateral positions of the SW were observed for both models, and a subsequent analysis revealed distinct placements once a stable interface was established. At a lateral distance (x) of 0.4 m, the SW for Model A reached a height of 0.17 m, while for Model B it was just 0.10 m. This was probably due to the existence of a clay layer that hindered the development of the contact. According to our observations, the clay layer behaves like a natural barrier, which does not allow the saline water easy passage into the freshwater zone, depending on the degree of saturation within the layer.

After reaching a steady condition, the SW was observed to have different heights for Models A and B, with δ representing the relative difference in their height. The development of saline contours with respect to length (m) and their respective heights (m) was predominantly influenced by subsurface heterogeneity, notably the presence of a clay layer (Figure 8). The advancement of the SW depends on the medium’s geological composition, specifically whether it moves smoothly or encounters any obstacles. The key factors include the porosity (ϕ) and permeability (k) of the geological layer. A tightly packed layer with less pore space will not be easily penetrated by the saline water, whereas a loosely packed layer with more pore space will allow the saline water to penetrate it more easily.

The experimental behaviour exhibited a two-dimensional salinity incursion pattern across various geological conditions. A stable equilibrium was achieved between the freshwater and saltwater under steady-state conditions despite variations in the shape of the interface. However, it is important to understand that the occurrence of three-dimensional heterogeneity in coastal aquifers may not behave in the same way as two-dimensional heterogeneity. An effort was undertaken to investigate the response of various geological layers to SWI.

3.2. Numerical Model Based on Experimental Feedback

In this section, the experimental results obtained were compared to the numerical model to establish the relationship between the two. The numerical model was created for the homogeneous and heterogeneous scenarios using the initial model parameters (Table 1) shown in Figure 9. For Model A, the flow of the SW was smooth as it advanced in the homogenous media (Figure 9a). The 17% contour line reached its lateral position (x) at 0.14 m after 30 min (Figure 9c). For Model B, the 17% contour line crossed the lateral position (x) of 0.10 m after 40 min but with a deviated shape (Figure 9d). For Model B, a numerical layer (L1) was created to have the similar properties of a geological clay layer, being both porous and impermeable (Figure 9b). To justify these properties, we assigned the upper and lower boundaries of the L1 layer as impermeable while keeping the left and right sides as open boundaries and allowing freshwater and saltwater flux to enter the layer from the left and right sides, respectively.

However, as shown in Figure 3b, the clay layer used during the experiment is tightly packed (sandwiched between an upper sand and a lower sand mixed with pebbles), reducing its porosity. The surrounding compactness affects a layer’s porosity (ϕ). Here, the clay layer behaves as an impermeable layer with a lower porosity. Meanwhile, there is a constraint with the numerical layer (L1) as its behaviour slightly differs from that of the experimental clay layer. This discrepancy arises from the boundary conditions (BCs) chosen for layer L1—impermeable on the upper and lower boundaries and open on the left and right sides. In contrast, an actual geological clay layer is porous and impermeable on all surface boundaries, making it difficult to assign the same boundary conditions in a numerical model. If we were to assign such boundary conditions to the numerical layer (L1), we would not be able to solve the 2D freshwater–saline water flow problem. Therefore, while a numerical layer cannot fully replicate the actual behaviour of a clay layer, an attempt can be made to observe the behaviour of the FSWI across impermeable boundaries.

The isochlor contour lines were plotted for both models with concentrations of 83%, 50%, and 17% (Figure 9c,d). These contour lines represent the decreasing concentration levels, expressed as a percentage, of the saline water flux (initially at concentration c = 1) as it permeates the medium against the freshwater flux.

The computational model exhibits a strong correlation with the experimental images in the case of Model A, where the FSWI crossed the centre of the experimental box after 30 min at a lateral extent (x) of around 0.14 m. In both the experimental and the model images, the SWs had relatively comparable heights at the centre of the vertical cross-section. For Model B, both the experimental and model images validate the deviation of the interface between freshwater and saltwater when it comes into contact with the clay layer (referred to as layer L1 in the model). However, there was a notable disparity in the height and the lateral extent of the SW between the experimental and model images, particularly evident at the centre of the profile. This significant variation in the SW can be attributed to an impermeable clay layer in the middle of the box. The behaviour of this layer differs between the experimental setup and the numerical model boundaries. For Model B, the numerical image depicts a greater SW height at the centre of the experimental box compared to that of the experimental image. The alignment between the experimental and numerical model results is illustrated in Figure 10. For Model A, E1, E2, and E3 denote experimental images taken at various x (m) positions, corresponding to numerical models N1, N2, and N3, respectively. Similarly, for Model B, E4, E5, and E6 represent experimental images captured at different x (m) positions, corresponding to numerical models N4, N5, and N6. These observations lead to two significant findings. First, it was observed that when advancing salt contours encounter impermeable layers like clay, the interface deviates from its initial trajectory. Additionally, the porous characteristics of the layer (such as clay) were observed to influence the height of the saltwater wedge. Specifically, when the layer is densely packed, the height of the saltwater wedge decreases at the same lateral distance. In contrast, a clay layer with a high level of porosity can lead to a greater height of the saltwater wedge, allowing for a smooth build-up of the interface.

3.3. Vertical Electrical Sounding (VES)

The DC resistivity data were obtained using the VES method once a steady-state equilibrium was achieved in both Model A and Model B. Data were collected along the center of the profile using a dipole–dipole array, with each reading noted in both the direct and reverse modes of the current (I) direction. The average value was plotted for each particular data point for better accuracy. Each observed data point for both Models A and B in Figure 11a,c represents the apparent resistivity value corresponding to the N-spacing (m) along the center of the experimental box. Figure 11b,d depict the one dimensional inverted response obtained from those data points, illustrating the variation in apparent resistivity values along the center of the profile with corresponding depth for experimental Models A and B.

For Models A and B, the inverted response is plotted using a smooth-layer model rather than a layered model. The noise-free data were acquired for model A due to having a consistent sand layer throughout its volume. In contrast, model B exhibited various layers, such as sand, clay, and sand mixed with pebbles, which can lead to an abrupt change in apparent resistivity values and thus may lead to noisy data points. If a layered model had been used, we would have encountered a high fitting error for the noisy points and would not have accurately represented the actual variation in resistivity values with respect to depth. Therefore, smooth layers were used for both models to obtain a better fit for such noisy points. Also, for a better comparison between models A and B, we maintained both models as smooth-layer models.

Model A shows a consistent and gradual change in the resistivity data due to the homogeneous geological medium (sand) within the box. The resistivity value obtained for freshwater-saturated sand was approximately 54 ohm-m at a depth of 3 cm below the measurement surface of the box. The depth of each layer was measured from the point at which the resistivity data were initially measured using the potential electrodes. It was observed that beyond the 17 cm mark, the resistivity begins to decrease upon encountering the FSWI. Beyond this depth, the measured resistivity was 0.72 ohm-m due to the saline sand, which indicates that the saline sand exhibited a higher level of conductivity (less resistivity) than that of the freshwater-saturated sand. For Model B, initially, the data points exhibited smooth characteristics. However, sparse data points were observed at a certain depth, indicating noisy data excluded from the one-dimensional inversion process. The first layer consisted of freshwater-saturated sand, with a measured resistivity of 26.40 ohm-m at a depth of 4 cm. In the second layer, the resistivity decreased to 6.90 ohm-m, corresponding to a clay layer situated at an 8 cm depth. A significant decrease in the resistivity data was observed at a depth of 20 cm, indicating the presence of the FSWI. Below this depth, the measured resistivity was 1.77 ohm-m, which was attributed to the presence of saline sand mixed with pebbles.

The observed resistivity data vary depending on the medium. In the case of Model A, the saline sand exhibited a resistivity value of 0.72 ohm-m. In contrast, for Model B, the resistivity was measured as 1.77 ohm-m, which is slightly higher due to the surrounding effect of the sand mixed with pebbles. The sounding results demonstrate a strong correlation with the experimental data. In the case of Model A, the interface depth at the center of the profile was approximately 20 cm according to the experimental image and 17 cm according to the depth derived from the one-dimensional inversion model. For Model B, the interface depth was approximately 22 cm according to the experimental image and 20 cm based on the one-dimensional inverted model. Slight discrepancies were noted with depth in the interface observations, as VES is an indirect geophysical measurement, and model error can arise from factors like data noise, signal attenuation, and errors in one-dimensional inversion, such as equivalence and suppression errors. However, our findings still fall within an acceptable range, and we successfully identified the FSWI utilizing the dipole–dipole array.

3.4. Validation with ERT Field Data for Location 1: Near Sea Coast

The ERT field data were used to validate the experimental findings. The correlation between the experimental and numerical results was evaluated by analysing the ERT subsurface image. The two-dimensional profile was acquired along the coast of the Mandarmani area of West Bengal in India, as illustrated in Figure 12a. This region is a tourist hotspot well known for its salt industries, aquaculture practices, and agricultural lands. Over the past few years, there has been a significant increase in water demand in this region.

Additionally, the overpumping of groundwater has been identified as a concern. Salinity intrusion presents a considerable challenge for the region, arising from natural phenomena and human interventions. The study area is characterized by alluvial deposits comprising sand, silt, and clay. Furthermore, an analysis of the borehole and lithology data in the area reveals that surface sediment deposits predominantly comprise medium-to-coarse sand interspersed with clay and small patches of dune sand [20]. The region is abundant in borewells utilized for agricultural practices, primarily for cultivating paddy crops and household water needs.

The ERT data were collected over a 160 m profile using a Wenner array, with electrodes positioned at 2 m intervals. This choice was made deliberately to maintain a fine electrode spacing to achieve high-resolution images. The inversion routine is based on the smoothness-constrained least-squares method with the L2 norm [60]. The two-dimensional inverted image reveals a prominent high-conductive zone with resistivity values ranging from 0.8 to 1.8 ohm-m, indicating the presence of a saline clay layer marked as a zone of SWI at a depth of 6 to 16 m below ground level. The sand layer at a depth of 20 to 26 m below ground level has a resistivity value of 6 to 7 ohm-m, as shown in Figure 12b. The area has resistivity values for a sandy freshwater formation that vary from 20 to 60 ohm-m. However, in ERT location 1, the shallow layers are completely saline (with resistivity values less than 1.8 ohm-m) for 6 to 16 m depths below ground level.

Saltwater intrusion naturally occurs near coastal areas. However, during the field survey, local people reported a high level of groundwater salinity in this coastal region. They have drilled borewells deeper than 20 m to access freshwater, as the shallow groundwater is completely saline. This issue arises from the extensive use of hand pumps and borewells in the confined coastal area. Overpumping lowers the water table, allowing saltwater from the sea to exert pressure and infiltrate deeper into freshwater aquifers. The ERT results confirm SWI in the coastal region, primarily caused by the disruption of the natural balance between saline and freshwater interfaces due to overpumping by residents.

Interestingly, it was observed that saltwater is trapped in unconfined shallow clay layers, preventing the further deterioration of groundwater in deeper sand layers. This study suggests that these clay layers can act as natural barriers, trapping saline water due to their porous but impermeable nature. However, the clay layer can only hold saline water up to a certain threshold. Excessive overpumping will increase salinity intrusion, causing the clay layer to fail to retain saline water and potentially impacting the surrounding layers. Similarly, the experimental results show that the tightly packed clay layer prevented the saltwater wedge from easily passing through, resulting in a deviated interface shape due to its impervious nature.

3.5. Validation with ERT Field Data for Location 2: Aquaculture Area

The ERT data were acquired in the aquacultural ponds with an area of 2.25 km² near the Mandarmani-Contai region of West Bengal in India, where they are artificially formed with seawater collected through a network of interconnected canals (Figure 13a). The area is committed to aquaculture practices focused on a unique breed of fish that thrives exclusively in saltwater environments. The geological area resembles fluvial deposits with clay, silt, and sand layers. The fish farming practices in this region serve as a crucial economic aspect of the local community’s livelihood. Nonetheless, the extensive nature of such production methods also exerts a detrimental impact on environmental health.

The ERT data were collected over a 800 m profile using a Wenner–Schlumberger array with 10 m electrode spacing. The data were inverted using the smoothness-constrained least-squares method with the L2 norm. The two-dimensional inverted section reveals the presence of a saline zone (with levels of resistivity ranging from 0.9 to 2 ohm-m) extending to a depth of 40 m below ground level (Figure 13b), which correspond to silty clay layers. The deeper zone is indicated by a sand layer with resistivity values of more than 12 ohm-m. The resistivity values for the sandy freshwater formation in the area range from 20 to 60 ohm-m. It can be observed from the ERT result that deeper layers (more than 50 m below ground level) are not affected by saltwater intrusion as resistivity values are greater than 12 ohm-m. However, considering shallow layers, the likely cause of the heightened salinity in this area is attributed to aquacultural practices, particularly fish farming. In such practices, seawater is transported to artificial ponds through connected canals from coastal regions. A considerable portion of this saltwater infiltrates the subsurface. The two-dimensional section illustrates that saline water permeates silty clay layers. However, no further penetration of salt water into deeper sand layers was observed, likely due to the silty clay layer acting as an impermeable layer.

Researchers worldwide have reported that SWI has contributed to the shrinkage of lakes [61], degrading water quality and threatening freshwater coastal resources [62]. The present study can be useful as its originality lies in integrating experimental data with numerical modeling and geophysical techniques, offering a comprehensive analysis of SWI dynamics in homogeneous and heterogeneous models. This approach enhances our understanding SWI mechanisms, informing coastal groundwater management by predicting the impacts of groundwater extraction and aquaculture on SWI. The findings can be applied to other coastal regions with similar geological characteristics to the study area in West Bengal, India, such as the Nile Delta [63], the Gulf of Mexico [64], and the Mekong Delta [65], which face similar challenges of saltwater intrusion due to intensive groundwater extraction and agricultural practices. By adopting this integrated methodology, tailored groundwater extraction policies and land-use planning can be developed, improving the sustainability of coastal groundwater resources. This approach can inspire future research and the development of sustainable groundwater management practices in coastal areas globally.

3.6. Limitations of This Experimental Study and Future Directions

The present study introduces new insights by incorporating heterogeneous background materials in the laboratory simulation of SWI. In this study, the alignment of experimental and numerical model results was quite good, though some uncertainties remain due to the following factors: (1) The numerical models could better align with the experimental findings if the same type of filler material, such as clay, had been used in the numerical simulation. It can provide better results as, in our case, the layer (L1) with constrained boundary conditions was not able to fully justify the actual behavior of a clay layer; (2) Instead of using VES for the experimental study, which provides a one-dimensional inverted model, a two-dimensional ERT method could be adopted. This method offers a better resolution and accuracy in predicting the depth of the FSWI. Some researchers have used a multielectrode ERT setup for laboratory experiments [58], yielding promising results as they can observe the high-resolution two-dimensional resistivity distribution rather than the one-dimensional resistivity change. These one-dimensional inverted models can be affected by acquisition errors, cultural noise, and inversion errors, leading to false interpretations and inaccurate estimates of the actual depth of layers. This discrepancy is evident in our study, in which the FSWI depth is slightly mismatched between the experimental and one-dimensional results. These issues may need to be addressed to reduce uncertainties associated with such experimental studies.

4. The Significance of Research Outcomes in Informing Policies for Coastal Water Management and Sustainable Aquaculture Practices

The current research findings are crucial for shaping policies related to coastal water management and promoting sustainable practices in aquaculture ponds and lakes. In the study area of West Bengal, India, groundwater overpumping and aquaculture ponds have emerged as significant concerns due to elevated salinity levels in the region. By providing scientific evidence and insights into the environmental dynamics of coastal areas, research helps policymakers make informed decisions. Understanding factors such as saltwater intrusion, pollution levels, ecosystem health, and the impact of human activities allows for the development of effective management strategies. Moreover, research contributes to identifying best practices for sustainable aquaculture management, including aquaculture techniques, specifically fish farming, that minimize environmental degradation and maximize productivity. By integrating research outcomes into policy formulation, governments and organizations can work towards safeguarding coastal waters and preserving the ecological balance of ponds and lakes for future generations. The flowchart below illustrates the connection between sustainable water management (SWM) and sustainable aquaculture practices, which is crucial for the overall sustainability of the environment (Figure 14).

This study highlights the negative effects of artificial ponds and lakes used for aquaculture, especially in causing groundwater salinization. Effective governance of aquaculture practices is crucial in implementing measures that balance economic benefits with environmental preservation, ultimately contributing to sustainable water management efforts. This study employs an experimental approach to investigate salinity intrusion in diverse coastal environments. The ERT results validate the impact of SWI due to various anthropogenic activities. Additionally, such studies can aid in developing decision-making tools for policymakers to maintain coastal aquifer sustainability [66].

The present study is also relevant to hydrogeology and environmental engineering, particularly for coastal groundwater management and mitigating saltwater intrusion. The integrated approach combines experimental data, numerical modeling, and geophysical techniques and offers valuable insights that can directly inform policymaking and sustainable coastal water management practices. By demonstrating how different geological conditions influence saltwater intrusion dynamics, this research provides a scientific basis for developing targeted groundwater extraction policies, implementing effective land-use planning, and designing physical barriers to protect freshwater resources. Emphasizing these potential applications can help policymakers and practitioners adopt more effective strategies for managing coastal groundwater resources and mitigating the adverse impacts of saltwater intrusion.

However, while the present study provides a better understanding of SWI behavior in a heterogeneous coastal environment, some mitigation measures still need to be physically implemented for the better prevention of SWI. Researchers have focused and aligned their experimental studies on physical barriers to minimize SWI [27,29]. Some physical and hydraulic management approaches can be used to mitigate SWI, such as abstraction barriers, cutoff walls, recharge wells, and tidal regulators. When exploring suitable locations for recharge wells, ERT combined with TDEM surveys can provide a better picture of the subsurface. Therefore, more geophysical studies are needed to complement laboratory experiments.

5. Conclusions

The present study employs an integrated approach, combining experimental, geophysical, and numerical modeling methodologies to gain insight into the dynamics of salinity intrusion in coastal aquifers. The experimental setup outlined in this study simulates various geological coastal scenarios by incorporating saltwater boundary conditions. A good agreement was found between the experimental results and models obtained from the numerical simulation and the DC resistivity method. The two critical findings indicate that the progression of the saline contours depends on the geological composition through which they propagate. First, homogeneous formations, such as sand, facilitate the unhindered infiltration of saline water into freshwater aquifers. Second, heterogeneous media featuring layers of sand, silt, and clay, especially impervious clay layers, act as natural barriers, impeding the advancement of saline water by trapping it within their structure.

This study offers valuable information to help us understand saltwater wedges’ interaction and solute transport mechanisms under diverse geological conditions. In addition, this study serves as a valuable benchmark for validating density-coupled flow and transport models, particularly those that incorporate flux-type boundary conditions using different background materials. The ERT data also confirm SWI in coastal areas resulting from anthropogenic activities, such as groundwater overpumping (location 1) and aquacultural activities, leading to saline water infiltration into subsurface aquifers (location 2). This study also examines the behavior of clay layers, which intriguingly can function as natural barriers to salinity intrusion issues. However, it cannot be viewed as a complete solution to saltwater intrusion, as clay layers have saturation thresholds. Beyond these thresholds, they allow saline water infiltration into the surrounding layers. External physical barriers are necessary to prevent saline water intrusion into freshwater aquifers. Furthermore, to ensure groundwater sustainability, a continuous subsurface investigation is needed for monitoring the health of coastal aquifers. A comprehensive approach that can be adopted involves integrating a three-dimensional time-lapse ERT survey with a geochemical analysis.

The experimental findings hold validity for geological environments characterized by clay, silt, and sand layers and are applicable to coastal and river depositional settings, such as the coastline margin area of India, where the biggest challenge is saltwater intrusion. The present experimental model incorporating a heterogeneous geology can serve as input for future investigations. It can be enhanced by introducing physical barriers, pumping wells, recharge wells, and inclined slopes, facilitating multiple attributes for detailed experimental studies of SWI. The findings can be applied to other coastal regions with similar geological characteristics, such as delta regions worldwide, which face similar challenges of saltwater intrusion due to intensive groundwater extraction and agricultural practices; further remedial measures can be planned accordingly.