Abstract

This study examines the level set method with reinitialization (LSMR) for monitoring the interface in saltwater intrusion problems. The governing equations consist of a parabolic equation modeling groundwater flow and a hyperbolic equation (the level set equation) that tracks the evolution of the freshwater/saltwater interface. The level set equation used for interface tracking is based on the sharp-interface approach, wherein a regular scalar function is defined across the interface using the level set method. The weighted essentially non-oscillatory scheme and a third-order total variation diminishing Runge-Kutta scheme are used to solve the nonlinear partial differential equations governing flow, minimizing mass loss. These schemes are implemented in MATLAB. Numerical simulations for Glover’s problem, the Madras aquifer, and the Kawatani problem demonstrate that the results of the LSMR model closely align with those reported by other researchers. The study highlights that significant saltwater intrusion remains a major concern for coastal aquifers, which play a crucial role in supplying water for domestic, industrial, and agricultural needs. Accurate tracking of the sharp interface is essential to mitigate potentially harmful consequences.

Keywords

Seawater Intrusion, Level Set Method, Sharp Interface, Glover Problem, Madras Aquifer, Kawatani Problem

Share and Cite:

1. Introduction

Aquifers provide drinking water for half the world’s population, and one in five is overexploited [1]. In regions such as northern China, where aquifers are heavily used for intensive irrigation, the water table has fallen by more than 40 m in just a few years. Excessive extraction not only increases the risk of landslides but also promotes saltwater intrusion, eventually rendering the water unsuitable for drinking [2].

Saltwater intrusion occurs in virtually all coastal aquifers, as they are hydraulically connected to the sea. Fresh groundwater from inland areas gradually drains into the ocean, while saltwater simultaneously enters the aquifer, establishing a dynamic equilibrium between the two fluids. This equilibrium forms a boundary or interface between freshwater and seawater. Frequent displacement of this interface from its natural position is largely due to the overexploitation of aquifers, which is the primary driver of saltwater intrusion. Additional contributing factors include the unregulated extraction of surface water, various forms of pollution, untreated wastewater discharge, and rising sea levels. These factors destabilize the natural equilibrium of the freshwater/saltwater interface and significantly accelerate the rate of saltwater intrusion [3].

Under natural conditions, hydraulic exchanges between freshwater and saltwater are typically slow and may be approximated by a quasi-equilibrium between the two fluid phases. This equilibrium is represented by the Ghyben-Herzberg approximation, which has led to several analytical models based on the relationship between interface height and sea-level height. However, these analytical solutions are limited to simple geometries and to conditions where the Ghyben-Herzberg assumption holds [4] [5]. These solutions are primarily used as test cases to validate numerical codes. In more extreme situations caused by meteorological events or human activities such as intensive freshwater pumping, the assumption of a static saltwater zone becomes inappropriate. These events cause the water table to drop, facilitating seawater intrusion into the aquifer [6].

Freshwater and saltwater differ in their physical properties, particularly in density; this difference causes the lighter freshwater to float atop the denser saltwater, forming a mixing region known as the transition zone. The geometry and thickness of this transition zone are influenced by factors such as the aquifer’s hydraulic properties and geometry, transport parameters (e.g., diffusivity), recharge conditions, and other hydrodynamic parameters. The interface’s shape and evolution also depend on the degree of heterogeneity within the aquifer [3]. Saltwater intrusion has been studied extensively since 1888, with numerous analytical, experimental, and numerical models proposed. Two main modeling approaches are commonly used: the sharp-interface approach and the variable-density (or density-dependent) approach.

The sharp-interface approach assumes that the fluids are immiscible and models a distinct interface between them, maintaining hydrostatic pressure continuity at the interface [7] [8]. This model assumes that the thickness of the transition zone resulting from salt diffusion into freshwater is small relative to the aquifer’s dimensions, making the approximation reasonable. This assumption has proven robust and reliable in the absence of significant external forcing. Multiple studies have used this approach [9]-[11]. In contrast, the variable-density approach accounts for a transition zone where freshwater and saltwater mix [12] [13]. It involves solving a convection-dispersion transport equation for salt within the aquifer.

The variable-density approach is considered more realistic because freshwater and saltwater are miscible. However, this approach leads to a system of strongly coupled, nonlinear parabolic partial differential equations, making both analytical and numerical solutions more complex and computationally expensive [14]-[16]. Choquet et al. [17] proposed a new sharp-diffuse-interface model for seawater intrusion in unconfined aquifers, combining the computational efficiency of the sharp-interface method with the realistic representation of the diffuse-interface method. Abudawia et al. [18] modeled a sharp-diffuse interface for the seawater intrusion problem in unconfined aquifers and analyzed a partially implicit time discretization scheme using a P_k (k ≥ 1) Lagrange finite-element approximation. They confirmed that the method achieved first-order accuracy in time and k^th-order accuracy in space. The authors also proposed a finite-volume technique for a structured grid and compared the results from these two approaches.

The sharp-interface approach requires explicit consideration of interface conditions. The flow equations modeling saltwater intrusion may be solved separately in the freshwater and saltwater zones, or in the freshwater zone alone under the assumption that the saltwater is hydrostatic. The interface between freshwater and saltwater is then derived from pressure continuity between the two phases. Although sharp-interface models do not resolve the transition zone, they provide critical information regarding the dynamics of the interface tip and toe, or the salt bevel, defined as the intersection of the interface with the bedrock and the top of the aquifer.

In [19], a two-dimensional saltwater intrusion model based on the sharp interface approach and the level set method was introduced. This formulation couples a parabolic equation for groundwater flow with a hyperbolic equation (the level set equation) to track interface evolution. Validation was initially performed using standard and modified Henry problems, which rely on modeling variable density in confined and homogeneous aquifers. The present study extends this work by applying the LSMR model to additional benchmark problems. First, the Glover problem, one of the most widely used for validating sharp interface models, allows for the evaluation of the toe and tip positions of the interface without considering dispersion. Second, the Madras aquifer problem incorporates longitudinal and transverse dispersion coefficients, not included in the Henry formulation, thus testing the model under more realistic transport conditions. Third, the Kawatani problem examines anisotropy in unconfined aquifers. Indeed, studies by Voss and Souza [20] and Croucher and O’Sullivan [21] highlight the considerable influence of parameters such as dispersion and porosity on the position of the freshwater-saltwater interface.

While numerous studies have modeled saltwater intrusion considering various influencing factors, the level set method with reinitialization (LSMR) has received relatively limited attention. Therefore, this study focuses on monitoring the interface between freshwater and saltwater using LSMR by unifying both saltwater and freshwater phases within a single flow equation. The method introduces a regular scalar function, called a level set function, which solves a first-order hyperbolic transport equation. The zero level set of this function defines the freshwater/saltwater interface. The method is applied to the Glover problem (1959) [22] as presented by Polo and Ramis [23], the Madras aquifer case (1990) in India, and the Kawatani problem (1980) [24] [25]. The LSMR model has already demonstrated its effectiveness in the context of the Henry problem (1964) [19].

2. Presentation of the LSMR Model

The LSMR model is a two-dimensional model developed by Meya Loua-Bouayi et al. [19] for monitoring freshwater/saltwater interfaces in coastal aquifers. This method assumes that the two fluids are immiscible. With its various features, the level set method overcomes the classic limitations of the sharp-interface approach, namely, the lack of information on the freshwater/saltwater interface geometry, instability in evolving states, and the difficulty of tracking the interface when solving equations in different phases with a moving mesh.

The level set method characterizes the interface Γ as the zero-level curve of a regular function or level set function $\varphi \left(x,t\right)$ : $\Omega \times ℝ\to ℝ$ , where Ω is an open subset of ${ℝ}^2$ . That is, at a time $t$ ,

$\Gamma \left(t\right)=\left\{x\in \Omega /\varphi \left(x,t\right)=0\right\}.$ (1)

The interface is thus implicitly defined by the level set function $\varphi$ as follows:

$\{\begin{array}{l}\varphi \left(x,t\right)<0\text{if}x\in {\Omega }_f\\ \varphi \left(x,t\right)>0\text{if}x\in {\Omega }_s\\ \varphi \left(x,t\right)=0\text{if}x\in \Gamma \end{array}$ ,(2)

with $\Omega ={\Omega }_f\cup {\Omega }_s$ , where ${\Omega }_f={\Omega }^{-}$ is the freshwater domain, and ${\Omega }_s={\Omega }^{+}$ is the saltwater domain.

The movement of the interface is obtained by determining the time evolution of the level set function. Instead of tracking the motion of the interface points Γ directly, the method evolves the global surface represented by the level set function $\varphi$ . The interface is then determined at each time step by locating the set $\Gamma \left(t\right)$ , where the level set function $\varphi$ equals zero.

From a numerical standpoint, this global formulation is advantageous because topological changes such as breakups or reconnections are naturally captured without requiring complex treatment. Because velocity drives the motion of the interface, its displacement is computed by transporting the level set function $\varphi$ through the velocity field $v$ , obtained from Darcy’s law. This is implemented using the level set equation:

$\{\begin{array}{l}\frac{\partial \varphi }{\partial t}+ \nabla \cdot \left(v\varphi \right) = 0\text{in}\Omega \times \left]0,T\right[\\ \varphi \left(x,t=0\right)={\varphi }_0\left(x\right)\text{in}\Omega \underset{}{}\\ \nabla \varphi \cdot n=0\text{in}\partial \Omega \times \left]0,T\right[\end{array}$ (3)

The transport Equation (3) is hyperbolic, and its solution proves to be difficult. Classical numerical methods, such as Upwind or Lax-Friedrichs, while simple to implement, suffer from significant numerical diffusion and limited accuracy near discontinuities. These limitations lead to oscillations and mass loss, thus compromising the precise localization of the interface. To overcome this, the 5th-order WENO (Weighted Essentially Non-Oscillatory) scheme was adopted for spatial discretization, known for its ability to reduce numerical diffusion and maintain high accuracy near sharp interfaces. Furthermore, the 3rd-order Runge-Kutta TVD (Total Variation Diminishing) scheme was used for temporal discretization, ensuring stability and mitigating oscillations. The combination of these two schemes offers superior robustness and reliability, as widely demonstrated in the literature [3] [15], and represents a more sensible choice than the simpler methods cited above.

In a discretized framework, gradients of the level set function $\varphi$ must be calculated carefully. To ensure this, the level set function $\varphi$ is initialized as a signed distance function. Numerically, signed distance functions are well suited for differentiation owing to their regularity, allowing high-order precision.

Following the work of Sussman et al. [26] [27], efforts have been directed toward improving the mass conservation properties of the level set method, which are generally weak.

• The first approach to address this problem involves implementing robust numerical methods to solve the transport equation of the level set function $\varphi$ . High-order discretization schemes are required. Building on Osher’s work, the weighted essentially non-oscillatory (WENO) scheme is used for spatial discretization, and the total variation diminishing Runge-Kutta (TVD-RK) scheme is used for temporal discretization.

• The second drawback lies in the spreading or tightening of the contour lines over time, which leads to a loss of stability and alters the qualitative properties of the solution. To limit these distortions in the level set equation, the reinitialization algorithm proposed by Sussman et al. [26] is used. The reinitialization procedure is frequently used to restore a proper signed distance function by solving the following evolution equation iteratively:

$\{\begin{array}{l}\frac{\partial d\left(x,t\right)}{\sigma \tau }=\text{sign}\left(\varphi \left(x,t\right)\right)\left(1-\Vert \nabla d\left(x,\tau \right)\Vert \right)\\ d\left(x,\tau =0\right)=\varphi \left(x,t\right)\end{array}$ (4)

Here, $\tau$ is a fictitious reinitialization time, $d$ is the level set function at the fictitious time $\tau$ , and sign is a smoothed sign function defined as

$\text{sign}\left(\varphi \right)=\{\begin{array}{ll}-1\hfill & \text{if}\varphi <-\Delta x\hfill \\ \frac{\varphi }{\sqrt{{\varphi }^2+\min {\left(\Delta x,\Delta y\right)}^2}}\hfill & \text{if}\left|\varphi \right|\le \Delta x\hfill \\ 1\hfill & \text{if}\varphi >\Delta x\hfill \end{array}$ .(5)

Equation (4) is introduced to construct, from an initial ${\varphi }_0$ defining the interface $\Gamma$ (with $\Vert \nabla \varphi \Vert \ne 1$ ), a new level set function $\varphi$ that captures the same interface but satisfies $\Vert \nabla \varphi \Vert =1$ .

Discretization errors from numerical schemes can significantly distort the interface position. To reduce this, low-diffusive, high-order schemes are employed. In this work, ENO1 (first-order essentially non-oscillatory) and ENO2 (second-order) schemes are used for spatial discretization, along with the explicit Euler scheme for time discretization. Stability is maintained using Courant-Friedrichs-Lewy conditions. This combination was selected to optimize computational time [28].

Solving the level set equation with a signed distance function has been shown to improve result accuracy, particularly regarding mass conservation, and enhance the precision in computing normal vectors and interface curvature [28]-[30].

ENO schemes use polynomial interpolation to compute fluxes at cell boundaries by choosing the smoothest possible stencil, achieving high-order accuracy even near discontinuities and avoiding spurious oscillations. These interpolations are based on sets of neighboring points (sub-stencils), which form the overall stencil for flux approximation.

The WENO scheme builds upon this principle by combining multiple stencils into a weighted average, with weights based on the local smoothness of the solution [31].

The governing equations for the flow are based on the continuity equation combined with Darcy’s law [3] [32]-[34], expressed as

$\{\begin{array}{l}{S}_{\alpha }\frac{\partial H_{\alpha }}{\partial t}+\nabla \cdot \left({\rho }_{\alpha }{v}_{d_a}\right)={\rho }_{{\alpha }^Q}\text{with}\alpha =f,s\\ v_{d_{\alpha }}=-K_{\alpha }\nabla H_{\alpha }\end{array}$ (6)

where $S_a$ , $H_{\alpha }$ , $K_{\alpha }$ , $Q$ , and $v_{d_{\alpha }}$ represent the specific storage coefficient, hydraulic head, hydraulic conductivity, source term, and Darcy velocity for each phase $\alpha$ (with $\alpha =f$ for freshwater and $\alpha =s$ for saltwater).

The model assumes the fluids are immiscible, which introduces an interface containing discontinuous parameters. Additional modeling conditions are therefore required. The continuum surface force (CSF) method is applied to represent the surface tension $\overrightarrow{F}$ between the two fluids [35]:

$\overrightarrow{F}=\sigma Cn\delta =\sigma C\left(\varphi \right)\delta \left(\varphi \right)\nabla \varphi$ ,(7)

where $\sigma$ is the surface tension coefficient, $n$ is the unit normal vector on $\Gamma$ directed from the saltwater domain ${\Omega }_s$ to the freshwater domain ${\Omega }_f$ , $C$ is the curvature of the interface, and $\delta$ is the Dirac delta function.

The CSF method does not directly handle discontinuities; instead, it regularizes them using Heaviside and Dirac functions [19].

By coupling the flow equations for both phases (Equation (6)) with the additional interface condition (Equation (7)), the following equation is obtained for modeling saltwater intrusion in porous media:

$S_{ϵ}\left(\varphi \right){\rho }_{ϵ}\left(\varphi \right)\frac{\partial H}{\partial t}=\nabla \left({\rho }_{\epsilon }\left(\varphi \right)k\nabla H+\frac{k\sigma C\left(\varphi \right){\delta }_{ϵ}\left(\varphi \right)\nabla \varphi }{g}\right)$ .(8)

Here, the density ${\rho }_{ϵ}\left(\varphi \right)$ and storage coefficient $s_{ϵ}\left(\varphi \right)$ are physical parameters expressed as functions of the level set function $\varphi$ . The parameter $ϵ$ defines the fictitious thickness of the smoothed interface.

The finite-volume method is used for the spatial discretization of Equation (8), and the explicit Euler scheme is used for time discretization. For the hyperbolic equation (Equation (3)), the ENO and WENO schemes are applied in space, while a third-order TVD-RK scheme is used in time. Further numerical details are provided in the study by Meya Loua-Bouayi et al. [19].

3. Numerical Results and Discussion

This study tests the LSMR model on three well-known benchmark problems previously solved using two different modeling approaches:

1) The Glover problem (1959), using the sharp-interface approach (the solutions obtained by the LSMR model are compared with those of [23] and the analytical solution of Glover (1959) [22]);

2) The Madras aquifer problem, examined by [24] [25] [36]; and

3) A theoretical case studied by Kawatani [37].

The LSMR model is a two-dimensional framework developed to track the freshwater-saltwater interface using the level set method. Its formulation is based on two fundamental equations: a parabolic equation for groundwater flow and a hyperbolic equation for interface tracking. Unlike classical sharp interface models, which require solving coupled and nonlinear systems separately at each phase, the LSMR handles the entire domain with a single flow equation coupled to the transport equation. To assess the model’s robustness, the fundamental equations were adapted to the boundary conditions and physical assumptions specific to several benchmark problems. The Glover problem focuses on a sharp interface without dispersion, allowing for the evaluation of the toe and tip positions of the interface. The Madras problem incorporates longitudinal and transverse dispersion coefficients, thus testing the model under more realistic transport conditions. The Kawatani problem examines anisotropy in unconfined aquifers, further testing the model’s applicability.

These tests are chosen to analyze the flexibility of the LSMR model and demonstrate its ability to reproduce various saltwater intrusion scenarios.

The latter two cases were solved using the variable-density approach. The three tests were implemented using MATLAB.

The LSMR model has already demonstrated its effectiveness in addressing the Henry problem (1964) [38], which models saltwater intrusion in a confined aquifer using the variable-density approach [19].

3.1. Glover’s Problem (1959)

For over a century, the problem of saltwater intrusion into coastal aquifers has been studied both experimentally and theoretically, using the governing equations of the freshwater/saltwater interface. These studies are now used to predict the impacts of excessive freshwater pumping and to guide water resource management strategies in coastal aquifers.

The Glover problem, as presented by Polo and Ramis [23], is one of the most widely used test cases for validating saltwater intrusion models under the sharp-interface assumption. The assumptions used by Polo and Ramis [23] for solving this problem include the following: the aquifer is confined, isotropic, and homogeneous; the flow is two-dimensional; the initial interface and hydraulic head conditions are arbitrarily prescribed; and the boundary conditions at the domain edges (except the coastal boundary) are constant.

The parameters and boundary conditions used in this problem are shown in Table 1 and Figure 1 [23].

Table 1. Parameters used for the Glover problem (1959) [23].

In general, the boundary and edge conditions in the Glover problem are chosen arbitrarily. Figure 1 shows the conditions used by Polo and Ramis [23] and adopted in this study. These are constant conditions (for both hydraulic head and level set function), except at the coastal boundary, where inflow is defined. The initial interface is represented by the zero-level curve of the function.

${\varphi }_0\left(x,z\right)=-z+\left[\frac{-60}{\sqrt{575}}\sqrt{x}+60\right]$ .(9)

Figure 1. Initial interface condition used for Glover’s problem (1959) [23].

The analytical solution for Glover’s problem, as given by Polo and Ramis [23], is

$\xi \left(x\right)=z^b-\sqrt{\frac{2\delta Ux}{K}+{\left(\frac{\beta \delta U}{K}\right)}^2}$ , (10)

where $\delta =\frac{{\rho }_f}{{\rho }_s-{\rho }_f}$ .

At the interface, $\xi \left(x\right)=z$ , and according to the LSMR model, the freshwater/saltwater interface corresponds to the zero-level set function, i.e., $\varphi \left(x,z,t\right)=0$ . This allows the analytical solution (Equation (10)) presented by Polo and Ramis [23] to be rewritten as

$\varphi \left(x,z,t\right)={\left(z\left(k\right)-z^b\right)}^2-\left[\frac{2\delta Ux}{K}+{\left(\frac{\beta \delta U}{K}\right)}^2\right]$ .(11)

Figure 2 shows the initial interface, Glover’s analytical solution as given by Polo and Ramis [23], and the numerical result from the LSMR model at t = 240 d. Figure 3 shows the comparison between analytical and numerical solutions obtained by Polo and Ramis [23]. The LSMR model demonstrates satisfactory alignment with the interface tip and toe behavior. However, the interface generated by LSMR shows a minor deviation from the analytical solution and the numerical result shown in Figure 3, particularly in the form of reduced convexity in the interface curve.

Figure 2. Comparison between the interface obtained using the LSMR model (ENO1 scheme) and the analytical solution of Glover’s problem (1959), at t = 240 d.

Figure 3. Analytical and numerical solutions of Glover’s problem (1959) obtained by Polo and Ramis at t = 240 d. The interface moves from right to left [23].

3.2. Madras Aquifer Problem

The Madras aquifer is located near the Bay of Bengal in India. Madras is one of India’s most populous cities, with a population of approximately 11.8 million in 2023 (United Nations), and many industries are concentrated around the city. The Madras aquifer is subject to high water demand for domestic use, industry, and agriculture, serving nearby cities such as Minjur, Panjetty, and others (Figure 4).

This aquifer has been the subject of several theoretical, practical, and numerical studies; prior studies [24] [36] [40] have used the finite-element method to predict the evolution of the freshwater/saltwater interface. The Madras aquifer is confined, and in the study by Rouve and Stoessinger [36], it is considered homogeneous and isotropic for model simplification purposes.

To compare existing solutions with those obtained using the LSMR model, the parameter values used in this study are consistent with those in previous studies [24] [25] [36] [39]. The parameters and boundary conditions for this problem are presented in Table 2.

Figure 4. Location of the city of Madras [40].

Table 2. Parameters used for the Madras aquifer modeling.

In Figure 5, the interface is represented by the zero-level curve of the initial state function, defined as ${\varphi }_0\left(x,z\right)=x-200$ . The salt concentration is initialized with $C_s=35\text{mg}/\text{L}$ on the seawater side and $C_f=0\text{mg}/\text{L}$ on the freshwater side. The interface is thus located where the concentration equals 17.5 mg/L. Using the level set function, this interface is always represented by the zero-level set: $\varphi \left(x,z,t\right)=0$ .

Figure 6 shows various solutions obtained by previous studies [24] [25] [36]. The isochlors for 35, 17.5, 10, 5, and 2 mg/L are shown from saltwater to freshwater, i.e., from right to left. Given that the initial saltwater concentration is 35 mg/L and the freshwater concentration is 0 mg/L, the interface location corresponds to $C=17.5\text{mg}/\text{L}$ , with the curve positioned to the left of the right boundary of the domain.

Figure 5. Illustration of the Madras problem with boundary conditions [24] [40].

Figure 6. Results obtained for the Madras aquifer problem by Abd-Elhamid [25] (SUFT), compared with Rouve and Stoessinger [36] and Sherif et al. [24]. Isochlors for 35, 17.5, 10, 5, and 2 mg/L are shown from saltwater to freshwater (from right to left).

The solutions provided by Rouve and Stoessinger [36] and Sherif et al. [24] [40] were obtained using the finite-element method for a steady-state model, while the solution by Abd-Elhamid (saturated/unsaturated fluid flow and solute transport (SUFT)) [25] was derived from a transient-state model with velocity-dependent dispersion coefficients.

Figure 7 compares the interface obtained using the LSMR model (with the ENO1 scheme for the reinitialization procedure) with that reported by Sherif et al. [24]. The figure shows that the LSMR solution closely matches the results of [24] [36]. The LSMR model accurately reproduces the interface tip and toe movements.

Figure 7. Comparison between the interface obtained using the LSMR model (ENO1 scheme for reinitialization) and that obtained by Sherif et al. [24].

3.3. Kawatani Problem (1980)

The developed LSMR model is applied to a problem originally studied by Kawatani in 1980. The aquifer is considered unconfined, with the height of its free surface on the land side being 0.35 m above sea level, resulting in a hydraulic gradient of 7 × 10⁻⁴. Several studies, including Kawatani [37] and Sherif et al. [40], have addressed this problem using the finite-element method. In the present study, the same domain, hydraulic parameters, and boundary conditions as those chosen by prior studies [25] [37] [40] are adopted to facilitate comparison.

The parameters and boundary conditions used for this problem are presented in Table 3 and Figure 8.

In Kawatani’s study, the salt concentration is initialized with $C_s=35\text{mg}/\text{L}$ on the seawater side and $C_f=0\text{mg}/\text{L}$ on the freshwater side. Therefore, the interface is located at $C=17.5\text{mg}/\text{L}$ .

Table 3. Parameters used for the Kawatani problem [24] [37].

Figure 8. Illustration of the boundary conditions for Kawatani’s problem (1980) [37].

Figure 9 shows the various solutions obtained by Kawatani [37], Sherif et al. [40], and Abd-Elhamid [25]. The isochlors shown for 1, 0.5, 0.25, 0.05, and 0.0005, respectively, extend from saltwater to freshwater (i.e., from right to left). Abd-Elhamid [25] normalized the saltwater concentration in the initial state to 1 and the freshwater concentration to 0, while Sherif et al. [40] used 35 and 0 mg/L as saltwater and freshwater concentrations, respectively. Thus, the interface position corresponds to $C=0.5$ in the study by Abd-Elhamid [25] and $C=17.5\text{mg}/\text{L}$ in the studies by Kawatani [37] and Sherif et al. [40].

Comparing the results in Figure 9 with those obtained using the LSMR model in Figure 10 (in blue) reveals that the LSMR model with the ENO1 scheme reproduces the interface behavior, specifically the toe and tip of the salt wedge, accurately.

Figure 9. Results obtained by Abd-Elhamid (SUFT), compared with other studies [24] [25] [37]. For Abd-Elhamid [25], the isochlors are shown for 1, 0.5, 0.25, 0.05, and 0.0005, from saltwater to freshwater.

Figure 10. Comparison between interfaces obtained for the Kawatani problem using LSMR with ENO1 (blue) and ENO2 (red) at the final time t = 250 d.

The same numerical experiment was repeated (Figure 10) by modifying the numerical scheme used in the reinitialization algorithm, replacing the ENO1 spatial scheme with the ENO2 scheme. The resulting solution is shown in red. This solution exhibits a minor deviation from the LSMR result with the ENO1 scheme (in blue). This deviation can likely be attributed to differences in the numerical schemes used during reinitialization.

4. Conclusions

This study focused on modeling saltwater intrusion into coastal aquifers, a phenomenon that poses a significant risk to the sustainability of groundwater resources. By extending the LSMR model beyond the classical Henry problem to benchmark problems such as Glover, Madras, and Kawatani, its flexibility and robustness under diverse hydrogeological conditions were demonstrated. The reinitialization algorithm, combined with high-order schemes, improves numerical accuracy, although the computational cost remains a challenge in complex scenarios. This suggests the need for hybrid approaches in future work. Importantly, reliable monitoring of the freshwater-saltwater interface provides valuable information for sustainable groundwater management, aiding in decision-making regarding pumping strategies and well placement.

Validation against three reference problems—the Glover problem (1959) presented by Polo and Ramis [23], the Madras aquifer problem, and the Kawatani problem (1980)—shows that the LSMR model satisfactorily approximates existing variable-density solutions for saltwater intrusion. In particular, the movement of the freshwater-saltwater interface is well reproduced. The model’s ability to integrate pumping and recharge effects opens promising avenues for future research.

The crucial role of the reinitialization algorithm in improving the accuracy and stability of level set method solutions has been highlighted. Although mass losses or gains can still occur, the algorithm, combined with TVD Runge-Kutta and WENO schemes, significantly improves interface tracking. Nevertheless, in complex scenarios, the computational cost of repeated resets can become substantial. Future work could explore hybrid approaches, such as coupling the level set method with radial basis functions and neural networks (RBFNNs), to optimize the efficiency-accuracy ratio.

Beyond these methodological contributions, the model has important practical implications. Accurate monitoring of the freshwater-saltwater interface can inform sustainable pumping strategies, guide the optimal placement of coastal wells, and help mitigate the risks of saltwater intrusion. By providing reliable predictions under various hydrogeological conditions, the proposed framework strengthens decision-making tools for water resource managers and contributes to the long-term preservation of coastal groundwater quality.

Acknowledgements

The authors would like to thank the African Union Scholarship Scheme for sponsoring the PhD studies of the first of them under grant Mwalimu Nyerere (Ref: HRST/ED/606).

Authors’ Contribution

All authors contributed to the study’s conception, design, and writing. All authors read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References