Transient Evolution of Inland Freshwater Lenses: Comparison of Numerical and Physical Experiments

Abstract

Brackish to saline groundwater in arid environments encourages the development and sustainability of inland freshwater lenses (IFLs). While these freshwater resources supply much-needed drinking water throughout the Arabian Peninsula and other drylands, little is understood about their sustainability. This study presents a numerical model using the SEAWAT programming code (i.e., MODFLOW and the Modular Three-Dimensional Multispecies Transport Model (MT3DMS)) to simulate IFL transient evolution. The numerical model is based on a physical laboratory model and calibrated using results from simulations conducted in a previous study of the Raudhatain IFL in northern Kuwait. Data from three previously conducted physical model simulations were evaluated against the corresponding numerical model simulations. The hydraulic conductivities in the horizontal and vertical directions were successfully optimized to minimize the objective function of the numerical model simulations. The numerical model matched observed IFL water levels at four locations through time, as well as IFL thicknesses and lengths (R² = 0.89, 0.94, 0.85). Predicted lens degradation times corresponded to the observed lenses, which demonstrated the utility of numerical models and physical models to assess IFL geometry and position. Improved understanding of IFL dynamics provides water-resource exploration and development opportunities in drylands throughout the Arabian Peninsula and elsewhere with similar environmental settings.

1. Introduction

While saline groundwater systems in semi-arid and arid environments are generally considered a problem for water resource managers, the occurrence of shallow (<100 m) brackish to saline groundwater promotes the development of localized accumulations of subsurface freshwater referred to as inland freshwater lenses (IFLs) [1,2,3]. IFLs serve as alternative freshwater resources for drinking water, livestock management, and micro-oasis agriculture in several arid and semi-arid regions including the Middle East (e.g., Kuwait, Oman, United Arab Emirates (UAE), Saudi Arabia), Australia (e.g., Murray Basin, Queensland), Central South America (e.g., Paraguay), Central Asia (e.g., Turkmenistan), Africa (e.g., Zambia, Namibia), and elsewhere with similar environmental settings [2,3,4,5,6,7,8,9,10].

IFLs form and are sustained by runoff from infrequent, high-intensity rainfall events (>20 mm/h), which occur between every one to three years based on a correlation between the El Niño-Southern Oscillation and rainfall anomalies [11,12]. Runoff collects in low elevation geomorphic features (e.g., depressions, gullies, interdune swales) and infiltrates vertically and laterally through highly permeable sediments towards the saline groundwater [13,14]. The resultant variable-density condition forms a boundary surface or freshwater-saltwater interface, above which a freshwater lens forms [15] (Figure 1). Unlike coastal and oceanic island freshwater lenses, which are stabilized by bounding seawater on one or more sides, IFLs migrate atop the saline groundwater in the direction of the regional gradient changing in geometry and position through time [1,16].

Investigations into IFLs are limited likely due to their sporadic nature and smaller scale relative to the larger-scale freshwater aquifers commonly exploited in populous regions. However, Bedouin communities likely once utilized these resources along with other ancient irrigation techniques (i.e., runoff harvesting), but abandoned these systems due to changes in the governing authority, modern agricultural practices, and inadequate drainage [17]. Recent studies have identified more than 100 potential locations throughout the Arabian Peninsula with conditions favorable (e.g., climatic, geomorphic, hydrogeologic) for the development and sustainability of IFLs [18]. Investigations into the formation and transient evolution of IFLs offer freshwater resources sustainability and development opportunities throughout the Arabian Peninsula and elsewhere with similar environmental settings. Novel solutions for the management of freshwater resources in these regions are critical as the overreliance on desalinization and non-renewable fossil water has created risks for environmental, political, economic, and social conditions [19,20,21,22].

Freshwater lens studies typically occur in coastal settings where the lenses are surrounded by seawater [23,24,25,26]. Physical models [27], analytical solutions [23,24,28], and numerical models have demonstrated the relationships between recharge rate, hydraulic conductivity, differences in fluid densities, and coastal freshwater lens geometry. IFL investigations using geophysical techniques [5,7], analytical methods [8,14,29], and physical laboratory models [1] highlight IFL relationships between the physical and chemical properties (i.e., salinity, age) and geometry. However, none of these studies address the transient nature of IFLs over varying temporal and spatial scales. A three-dimensional, numerical model of the Raudhatain and Umm Al-Aish IFLs in Kuwait by Al-Weshah and Yihdego [30] used MODLFOW-SURFACT code to predict water table elevations over 22 years but does not account for the freshwater extent in the vertical or horizontal direction as it migrates laterally through space and time. Research into the geometry of IFLs is needed to adequately predict the position and sustainability of exploitable freshwater resources. The purpose of this research is to presents a variable-density, finite-difference numerical model that considers geometry (i.e., thickness, length) calibrated by and compared with shared results from physical model observations by Rotz and Milewski [1].

2. Materials and Methods

The goal of this study is to investigate the development and transient evolution of a topographically induced IFL over varying recharge rates by comparing three simulations from a newly developed numerical model with three corresponding simulations from a previously conducted physical model in Rotz and Milewski [1]. The physical model simulations conducted in Rotz and Milewski [1] represent the formation, migration, and degradation of IFLs using an acrylic sandbox model inspired by the Raudhatain IFL in northern Kuwait (Figure 2). The previous study investigated the effect of recharge on IFL development and compared results to the Dupuit–Ghyben–Herzberg analytical solution by Vacher [24]. The statistical analyses in Rotz and Milewski [1] showed the analytical solution overestimated IFL thickness (D) and had no correlation with the recharge rate (R) (ΔD = 278.57R + 0.0872, R² = 0.34, p = 0.1304). The study concluded that the differences in lens geometry between the physically modeled lenses and the analytical solution were due to the transient nature of IFLs, subsequently motivating the need for numerical studies that could address IFL changes in geometry through time.

For this study, the numerical model was designed by using the physical model dimensions, porous medium, boundary conditions, and initial conditions from Rotz and Milewski [1]. In addition, numerical model calibration and comparative analyses used the observational data from three physical model simulations. The observational data (i.e., water table elevation, thickness, length) was obtained from still photographic images generated during the previous study (Figure 3). Uranine tracer-dye added to the freshwater during the physical model simulations allowed for the digital measurement of IFL geometry (i.e., thickness, length) using ImageJ software. Results from the numerical model were then compared to the results from the previously conducted physical model simulations to determine if the numerical model adequately simulated the observed IFLs.

The numerical model was developed using the SEAWAT programming code to simulate saturated, variable-density groundwater flow and solute transport. The SEAWAT program is a density-dependent groundwater flow model utilized for the modeling of groundwater in aquifers with fresh and saline groundwater interaction [31]. SEAWAT couples MODFLOW-2000 [32] and the Modular Three-Dimensional Multispecies Transport Model (MT3DMS) to solve the three-dimensional groundwater flow and solute transport numerically with finite difference approximations [32,33,34,35]. The shell program is Visual Modflow. The numerical model was designed in a two-dimensional manner to represent a cross-sectional slice of the hydrogeological system of a topographically induced IFL simulated with the physical model.

SEAWAT solves the fluid flow equation in terms of the hydraulic head by combining a general form of Darcy’s law for variable density conditions [36] with a hydraulic conductivity tensor [37], and vertical flow components [38,39] (Equation (1)):

$$\nabla ·\left[\rho K\left(\nabla h+\frac{\rho -{\rho }_{\mathrm{f}}}{{\rho }_{\mathrm{f}}}\nabla z\right)\right]=\rho S_{\mathrm{s}}\frac{\partial h}{\partial t}+\theta \frac{\partial \rho }{\partial t}-{\rho }_{\mathrm{s}}{q}_{\mathrm{s}}$$

(1)

where $\nabla ·$ is the divergence operator, $\rho$ is fluid density (ML⁻³), K is hydraulic conductivity tensor (LT⁻¹), h is the freshwater head (L), ρ_f is the freshwater density (ML⁻³), z is elevation, S_s is the freshwater specific storage (L⁻¹), θ is porosity (1), and q_s is a source or sink (T⁻¹) of fluid with density ρ_s.

SEAWAT uses the MT3DMS model to solve the general form of the solute transport equation [34] (Equation (2)):

$$\frac{\partial C}{\partial t}=\nabla ·\left(D\nabla C-vC\right)-\frac{q_s}{\theta }{C}_s$$

(2)

where C is the solute concentration (ML⁻³), D is the solute dispersion tensor, v is the advective velocity (LT⁻¹), and C_s is the source or sink concentration (ML⁻³). The linking of Equations (1) and (2) requires a function to relate concentration and fluid density.

2.1. Numerical Model Design and Parameters

The transient flow model was developed with SEAWAT to assess IFL geometry and transience for 200 h. The domain of the numerical model encompasses a space of 2.0 m (L) × 0.5 m (H) × 0.10 m (W) (Figure 4). The numerical model did not include the upper 0.5 m of the unsaturated zone of the physical model, as these were considered inactive cells within SEAWAT. The bottom left side elevation of the domain is at 0 m. The grid consists of 100 columns to represent the length (x-axis), 25 rows to represent the height (z-axis), and one layer to represent the width (y-axis). The parameters of hydraulic conductivity (K_x, K_y, K_z), were used from measurements taken during the previously conducted physical laboratory model study with a constant head permeameter during the laboratory experiments from sediments comprised of medium to fine sand (d₅₀ = 0.55 mm) to represent the highly permeable gravel sediments in the upper layer of the Kuwait Group Aquifer [36]. The porosity ($\theta$) was calculated during the previous study using the sum of solid and pore volumes. Additional representative properties of variable-density values through porous media were selected after Konikow, et al. [40] and others, as shown in

Table 1. Initial model parameters.

Input Parameter	Value	Description
Kx, Ky, Kz	0.0015 m/s	Hydraulic conductivity
θ	0.39	Porosity
Ss	1 × 10−5	Specific storage
Sy	0.34	Specific yield
Dm	1 × 10−10 m2/s	Molecular diffusion coefficient
Δρ/ΔC	0.7143	Slope of fluid density to solute concentration
δL	1 × 10−5 m	Longitudinal dispersivity
δv/Δl	0.02	Vertical/longitudinal dispersivity

.

2.2. Boundary and Initial Conditions

The boundary conditions of the numerical model corresponded to the physical laboratory model conditions (Figure 4). The base of the domain included a no-flux boundary. The left side of the domain included a constant head boundary of 0.41 m (h₁) with a constant concentration of 35,700 mg/L for the full duration of the model (200 h). The right side of the domain included a general-head boundary with a general-head elevation at 0.405 m (h₂) and boundary-head elevation at 0.40 m for a distance of 2 m, which generated a slope of 0.25% to represent the gentle gradient. The top layer of the domain included a recharge boundary with a width (w_r) of 0.42 m. The initial condition of the freshwater head (h_f) was calculated to 0.4075 m based on the change in water table elevation between the left and right boundaries of the domain. The initial saltwater concentration (C_s) of 35,700 mg/L was selected for the brackish to saline water from the Raudhatain Depression (7000–50,000 mg/L) and used with the physical model simulations. Evaporation, similar to the previously conducted physical laboratory simulations, was considered negligible. A summary of the boundary conditions is shown in

Table 2. Boundary condition and concentration values.

Boundary Condition	Position (x,y)	Values
Constant head	0, 0.41	35,700 mg/L
General head	2, 0.405	calculated
Recharge	0.79, 1.21	0 mg/L
Evaporation	Entire domain	0 m/s

.

2.3. Observation Wells and Geometry Data

Observation wells for water table elevations were used within the numerical model. Measurements from the physical model simulations at four static locations along the cross-section provided data for the wells. Water table elevation measurements were obtained from the physical laboratory model simulations from the previous study when the lens reached a maximum thickness and then every hour until ten hours, which was the time observed for the physical simulations where the tank wall did not interfere with the lens geometry. Wells 1–3 were located beneath the recharge zone, and Well 4 was assigned on the right side or down gradient end of the model (Figure 4). The elevation of the well screens within the numerical model was set directly below the initial water table elevation of 0.40 m. IFL thickness and length measurements from the physical model were taken to compare with the numerical model.

2.4. Numerical Settings

The simulation scenario selections for the SEAWAT numerical engine are described here. The explicitly coupled or “one-time step lag” option was selected where the flow equation is calculated using the fluid density from the previous transport time step, as solute concentrations are not expected to change rapidly. Due to the transient nature of IFLs, the upstream-weighting algorithm was selected to calculate the density value based on the flow direction of a particular iteration. For the viscosity calculations, a single species was used with no temperature dependence.

The Preconditioned Conjugate Gradient (PCG2) solver package was selected for the flow calculations with the maximum inner and outer iterations set to 100 each. The head change criterion was set to 1 × 10⁻⁵ m because of the small range of measurement values from the physical laboratory model. The flow solver used the incomplete Cholesky preconditioning method.

The implicit finite-difference method was used with the Generalized Conjugate Gradient (GCG) package to solve the transport equation. An upstream finite difference method was used to solve the advection term. The maximum number of outer and inner iterations was set to 1 and 50, respectively. A modified incomplete Cholesky preconditioner was also used for the transport solver. The relative convergence criterion was set to 1 × 10⁻⁸. The lengths of transport time steps were calculated using a Courant number of 0.75. The time steps for the flow calculations for stress periods 1, 2, and 3 are 10, 50, and 50, respectively. The maximum number of transport steps was set at 3000 with a maximum step size set at 0.125 h. Simulation results were saved at every time step for both flow and transport calculations to use for analyses.

2.5. Simulations

Three physical laboratory simulations at varying recharge rates were selected to reproduce three numerical simulations (

Table 3. The initial recharge rate (R_i), time of recharge flux (T_q), time of observed maximum thickness (T_max), adjusted recharge (R_a) rates, and stress period times for the numerical model simulations.

Simulation	Ri (m/s)	Tq (s)	Tmax (s)	Ra (m/s)	Stress Periods Starts (1, 2, 3) (s)
1	4.72 × 10−5	1620	6480	3.33 × 10−5	1620, 6480, 720,000
2	6.48 × 10−5	1296	6120	4.72 × 10−5	1296, 6120, 720,000
3	7.86 × 10−5	1008	5400	6.67 × 10−5	1008, 5400, 720,000

). An example of Simulation 2 from the previously conducted physical laboratory simulations is shown in Figure 5. For the physical model simulations, recharge was applied at varying rates for 1 hour. For the numerical model, the recharge rates were adjusted to account for the infiltration travel time through the unsaturated zone (Table 3). Observations made of the time passed for the freshwater to initially flux across the water table (T_q) were recorded and subtracted from the time observed for each simulated lens to reach a maximum thickness (T_max). The initial recharge rates (R_i) used in physical experiments were divided by the calculated time differences to equal the adjusted recharge rates (R_a). The beginning of the second stress period was set to the observed time of initial flux, and the end of the second stress period was set at the time of the observed maximum thickness to match the maximum thickness times of the physical model to the numerical model. Three stress periods were assigned to each simulation to apply the recharge and represent (1) a period of no recharge, (2) a period of freshwater recharge applied across the recharge boundary, and (3) a period of no recharge for a total duration of 200 h. The duration of each stress period was set accordingly (Table 3). Freshwater was infused with uranine, which is a tracer dye that allowed for the measurements of IFL thicknesses and lengths for ten hours for each simulation before the tank wall interfered with the lens geometry, which was compared to the numerical model using concentration data exported from the software. Statistical model evaluation techniques were used to determine the accuracy of the simulated data [41].

3. Results

3.1. Model Calibration and Results

The numerical model successfully converged using the observational data from the three previously conducted physical model simulations and demonstrated IFL transient evolution (Figure 6). The numerical model calibration for each simulation was achieved through a parameter estimation tool (PEST) [42] and a trial and error approach improving estimates of the hydraulic conductivity in the x and z directions while minimizing the error in head values. This value is referred to as the objective function (Φ). PEST adjusted the model parameters until the fit between the physical observations and numerical outputs was optimized to minimize the weighted least squares. The optimization procedure successfully estimated the hydraulic conductivities in the x and z directions (

Table 4. The objective function (Φ), hydraulic conductivity estimates (K_x, K_z), and confidence intervals for the numerical simulations using the parameter estimation tool (PEST) tool.

Run	Φ (m)	Kx (m/s)	Kz (m/s)	Confidence Interval for Kx (m/s)	Confidence Interval for Kz (m/s)
1	4.8 × 10−4	7.6 × 10−4	7.0 × 10−4	7.6 × 10−4–1.2 × 10−3	2.5 × 10−5–1.9 × 10−2
2	5.5 × 10−3	3.6 × 10−4	4.3 × 10−4	2.5 × 10−4–5.3 × 10−4	2.1 × 10−5–8.5 × 10−3
3	1.5 × 10−3	4.3 × 10−4	9.0 × 10−4	5.3 × 10−4–1.0 × 10−3	5.8 × 10−5–4.7 × 10−3

).

3.2. Water Table Elevation Results

The comparison between the observed and simulated water table elevations is presented in Figure 7. Overall, the statistical values indicate accurate model simulation [41]. The slope and y-intercept of the best fit regression line indicate a slight overprediction of simulated water table elevations and a lag effect (m = 1.0477, b = −0.02). The coefficient of determination describes a strong degree of collinearity and a positive linear relationship between the simulated and measured data (R² = 0.89). The Nash–Sutcliffe efficiency (NSE) indicates a good fit between the observed and simulated data (NSE = 0.84). The root mean square error (RMSE) indicates a small error value appropriate for model evaluation (RMSE = 0.007). The percent bias (PBIAS) value suggests nearly zero model overestimation bias of water table elevation (PBIAS = −0.7).

3.3. Lens Geometry Results

The comparison between the observed and simulated IFL thicknesses and lengths is presented in Figure 8 and Figure 9. Overall, the statistical values indicate accurate model simulation for both geometric parameters [41]. For the thicknesses and lengths, the slope and y-intercept of the best fit regression line indicate a slight overprediction of simulated thicknesses and a slight to moderate lead effect (m = 1.0856, b = 0.0133) and (m = 1.0134, b = −0.0094), respectively. The coefficient of determination describes a strong degree of collinearity and a positive linear relationship between the simulated and measured data (R² = 0.94) and (R² = 0.85). The NSE indicated a good to moderate fit between the observed and simulated data (NSE = 0.65) and (NSE = 0.82). The RMSE indicates an error value appropriate for model evaluation (RMSE = 0.03) and (RMSE = 0.06). The PBIAS value suggests a slight to moderate model overestimation bias of IFL thickness (PBIAS = −16.3) and (PBIAS = −0.4). For lengths, Simulations 2 and 3 correlates more linearly; however, Simulation 1 deviates from the linear trend likely from a decrease in velocity during the physical model simulation from hydraulic conductivity variations due to the packing of the sand (Figure 8). A summary of these results is presented in

Table 5. Summary of model evaluation statistics for IFL water table elevation, thickness, and length for Simulations 1, 2, and 3. The evaluation statistics are defined as a measure of fit (R²), the Nash Sutcliffe Efficiency (NSE), root mean square error (RMSE), and percent bias (PBIAS).

Run	Slope-Intercept Form	R2	NSE	RMSE	PBIAS
Water table elevation	y = 1.0477x – 0.0244	0.89	0.84	0.007	−0.7
Thickness	y = 1.0856x + 0.0133	0.94	0.65	0.03	−16.3
Length	y = 1.0134x – 0.0094	0.85	0.82	0.06	−0.4

.

3.4. Velocity Results

The average velocities for each simulation were calculated along a transect on the x-axis (v_x) from the center layer and column to the last column of the downgradient or right side of the domain (x = 1.0–2.0 m, z = 0.25 m) for each simulation before, during, and after the recharge stress period as shown in

Table 6. Average groundwater velocity (v_x) measurements of the numerical model simulations before, during, and after the recharge stress event.

Run	vx Before Recharge Stress Period (m/s)	vx During Recharge Stress Period (m/s)	vx After Recharge Stress Period (m/s)
1	6.88 × 10−6	2.40 × 10−5	7.04 × 10−6
2	3.86 × 10−6	6.99 × 10−5	3.85 × 10−6
3	5.64 × 10−6	6.48 × 10−5	4.84 × 10−6

.

3.5. Degradation Results

The decrease of IFL thicknesses, referred to as IFL degradation, was compared between the physical and numerical simulations. As noted in the lens geometry results, the numerical model over-predicted the thicknesses of all three solutions for the duration of the simulation but consistently predicted the rate of decrease in thickness (Figure 10). The IFL was delineated by locations during the simulation, where the concertation values were less than 1000 mg/L. The changes in IFL thicknesses were compared between models with best fit exponential regression lines. IFL degradation was predicted for both the physical and numerical models by calculating the time required for the thicknesses to equal 0.01 m. A summary of the results is outlined in

Table 7. Lens degradation results between the physical model and numerical model simulations. Units have been converted from seconds to hours.

Run	Physical Model Regression	Degradation Time (h)	Numerical Model Regression	Degradation Time (h)
1	y = 1.699 × 10−1 e−0.069x	41	y = 2.023 × 10−1 e−0.072x	42
2	y = 2.697 × 10−1 e−0.047x	70	y = 3.017 × 10−1 e−0.05x	68
3	y = 2.597 × 10−1 e−0.1x	33	y = 3.185 × 10−1 e−0.087x	40

.

4. Discussion

The numerical model using the SEAWAT programming code successfully modeled the development and transient evolution of three IFLs produced with a physical laboratory model [1]. Although SEAWAT does not model variable-saturated conditions, the adjusted recharge rates were sufficient to adequately simulate water table elevations, thicknesses, and lengths based on the physical laboratory model observations. In all cases, the model evaluation statistics indicated that the numerical model over-predicted thicknesses and lengths under the assumption that the error variance occurs within the simulated values.

The numerical model parameters of hydraulic conductivity were refined with the PEST calibration tool, which estimated hydraulic conductivities within an acceptable range for fine to medium unconsolidated, heterogeneous sand. The initial laboratory measurement of hydraulic conductivity at 0.0015 m/s likely did not consider the heterogeneity and anisotropy present in the physical model sand due to the systematic, horizontal packing of sand during the model setup. However, the statistical evaluations of the water table elevations, thicknesses, and lengths show that the numerical model adequately reproduced the magnitudes of the measured data.

The average velocities estimated by the numerical model aligned with velocities measured from the physical model, as well as the average velocities reported by Al-Weshah and Yihdego [30] of 6.34 × 10⁻⁷–2.85 × 10⁻⁶ m/s. The numerical model showed an increase in average velocity during the recharge stress period by an order of magnitude (10⁻⁵ m/s), which aligns with the range of recharge rates used in the physical and numerical model simulations 3.33 × 10⁻⁶–6.67 × 10⁻⁶ m/s. Based on these values, freshwater recharge that entered the center of the Raudhatain depression is estimated to travel 4000 meters towards the periphery between 20 and 200 years, depending on the hydraulic conductivity. This estimate aligns with those by Yihdego and Al-Weshah [43] and Parsons Engineering and Construction Corporation [44] who reported an inferred flow velocity between 11 and 245 m/year based on ¹⁴C and ³H age dating, as well as an estimated groundwater velocity of 20–90 m/year based on a hydraulic conductivity range of 40–80 m/day. Estimates such as these inform geochemical studies that aim to approximate IFL water age, such as the study by Kuldzhayev [45], which determined the onset of freshwater accumulation of an IFL in the Karakum desert to be 3400 years. Investigations into the formation, geometry, and extent of freshwater lenses over long temporal scales are needed for inland and coastal environments to quantify the supply of drinking water and the effects of external changes such as precipitation and sea-level rise [13,30], but also geomorphological [46] and anthropogenic impacts [47,48].

Errors from the measurements taken from the physical model simulations were likely produced from multiple sources. The water table elevation or top of the IFL was challenging to estimate due to the capillary fringe. Preferential flow paths or heterogeneities of hydraulic conductivity likely developed from the packing of the sand during the setup of the physical model simulations, which required careful consideration when determining IFL length. The numerical model did not account for these heterogeneities or water retained in the pore space of the unsaturated area, which likely further contributes to errors in the numerical simulations.

Physical and numerical models are helpful tools to understand the stochastic relationships for subsurface flow and transport in variable-density aquifers in naturally and artificially recharged aquifers [49,50]. Numerical models such as SEAWAT and others (e.g., MODFLOW-SURFACT, SUTRA (Saturated-Unsaturated Transport)) allow for the incorporation of variable-density conditions and parameters such as recharge volume, recharge location, concentration, temperature, multiple species, and other variables in terms of solute transport and subsurface flow. Although freshwater lens geometry, including transience, can be determined with these models, additional tools and optimization approaches to automate the monitoring of lens geometry are needed to improve the understanding of IFL sustainability. Models should be developed to consider the freshwater extent in the vertical and horizontal directions to predict volume, location, and residence time. Numerical simulations that include concentration data may use PEST to improve model calibration, which in turn improves the position and shape of the lens. These advancements will be useful for studying the dependence of periodic recharge on IFL degradation for water resources management. Figure 11 shows a conceptual example of IFL sustainability after a second recharge pulse is applied during IFL degradation.

5. Conclusions

The purpose of this study was to create a conceptual framework for topographically induced IFL transient evolution using a numerical model and demonstrate the utility of IFL geometric analyses for resource management. We employed the SEAWAT programming code with Visual Modflow to simulate variable-density groundwater flow and solute transport within an inland setting based on observations made during physical laboratory model simulations [1]. The numerical model was calibrated against the observed water table elevations measured in the laboratory by adjusting the values of hydraulic conductivity. The numerical model simulated IFL transient evolution and demonstrated the usefulness of predicting IFL geometry over time for water resource management.

As the demand for water resources in arid environments intensifies in conjunction with climate change, novel approaches for water resources development are needed [51,52,53,54,55]. Desert communities have historically used inland and coastal freshwater lenses combined with ancient irrigation techniques (e.g., runoff harvesting, subsurface canals) for access to freshwater because the methods were uncomplicated and technologically inexpensive [56,57,58]. Similar systems have been identified within the deserts of the Middle East, Asia, Australia, South America, and Africa but have declined in use. A return to these systems with a focus on the interplay between meteoric recharge and the underlying saline groundwater offers development opportunities through modern and traditional methods to meet the growing demand for water resources in arid environments [18,57].

This study provides a useful example for investigating IFL dynamics for the discovery of new resources [18], adjustment of pumping schemes [47], and artificial recharge applications [50]. The SEAWAT programming code serves as an effective tool for investigating the complexity of IFL transient evolution over long temporal scales for the prediction of available water resources in drylands.