Next Article in Journal
Biosorption of Water Pollutants by Fungal Pellets
Next Article in Special Issue
Microplastic and Fibre Contamination in a Remote Mountain Lake in Switzerland
Previous Article in Journal
Effect of Intermittent Aeration in a Hybrid Vertical Anaerobic Biofilm Reactor (HyVAB) for Reject Water Treatment
Previous Article in Special Issue
Coupled Multifield Response to Coordinate Mining of Coal and Uranium: A Case Study
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Geology, University of Georgia, 210 Field Street, Athens, GA 30602, USA
2
Warnell School of Forestry and Natural Resources, University of Georgia, 180 East Green Street, Athens, GA 30602, USA
*
Author to whom correspondence should be addressed.
Water 2020, 12(4), 1154; https://doi.org/10.3390/w12041154
Submission received: 31 December 2019 / Revised: 23 March 2020 / Accepted: 9 April 2020 / Published: 17 April 2020
(This article belongs to the Special Issue Advances in Hydrogeology: Trend, Model, Methodology and Concepts)

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 (R2 = 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.

Graphical Abstract

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, R2 = 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)):
· [ ρ K ( h + ρ ρ f ρ f z ) ] = ρ S s   h t + θ ρ t ρ s q s
where · is the divergence operator, ρ is fluid density (ML−3), K is hydraulic conductivity tensor (LT−1), h is the freshwater head (L), ρf is the freshwater density (ML−3), z is elevation, Ss is the freshwater specific storage (L−1), θ is porosity (1), and qs is a source or sink (T−1) of fluid with density ρs.
SEAWAT uses the MT3DMS model to solve the general form of the solute transport equation [34] (Equation (2)):
C t = · ( D C v C ) q s θ C s
where C is the solute concentration (ML−3), D is the solute dispersion tensor, v is the advective velocity (LT−1), and Cs is the source or sink concentration (ML−3). 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 (Kx, Ky, Kz), 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 (d50 = 0.55 mm) to represent the highly permeable gravel sediments in the upper layer of the Kuwait Group Aquifer [36]. The porosity ( θ ) 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.

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 (h1) 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 (h2) 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 (wr) of 0.42 m. The initial condition of the freshwater head (hf) 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 (Cs) 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.

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−5 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−8. 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). 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 (Tq) were recorded and subtracted from the time observed for each simulated lens to reach a maximum thickness (Tmax). The initial recharge rates (Ri) used in physical experiments were divided by the calculated time differences to equal the adjusted recharge rates (Ra). 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).

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 (R2 = 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 (R2 = 0.94) and (R2 = 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.

3.4. Velocity Results

The average velocities for each simulation were calculated along a transect on the x-axis (vx) 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.

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.

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−7–2.85 × 10−6 m/s. The numerical model showed an increase in average velocity during the recharge stress period by an order of magnitude (10−5 m/s), which aligns with the range of recharge rates used in the physical and numerical model simulations 3.33 × 10−6–6.67 × 10−6 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 14C and 3H 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.

Author Contributions

R.R., A.M., and T.C.R. conceptualized and designed this project. R.R. and A.M. acquired the financial support for this project leading to the manuscript. R.R. and A.M. developed the methodology. R.R. conducted the experiments, performed the formal analysis, and prepared the initial draft of the manuscript. A.M. provided the software and computing resources. A.M. and T.C.R. provided oversight and mentorship for the research activity. R.R., A.M., and T.C.R. critically reviewed and prepared the final draft of this manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Society of Exploration Geophysicists Lynn and Gary Servos Groundwater Exploration Scholarship and the University of Georgia Department of Geology Gilles Allard Research Grant.

Acknowledgments

The authors would like to thank the University of Georgia Water Resources and Remote Sensing Laboratory and John Dowd for assistance with the research, as well as the reviewers who invaluably helped with the improvement of this manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Rotz, R.R.; Milewski, A.M. Physical modeling of inland freshwater lens formation and evolution in drylands. Hydrogeol. J. 2019, 27, 1597–1610. [Google Scholar] [CrossRef]
  2. Kunin, V. Local waters in deserts and problems in their utilization. Moscow Izd. Akad. Nauk SSSR 1959, 2, 213–231. [Google Scholar] [CrossRef]
  3. Bergstrom, R.E.; Aten, R.E. Natural recharge and localization of fresh ground water in Kuwait. J. Hydrol. 1965, 2, 213–231. [Google Scholar] [CrossRef]
  4. Shevchenko, N. Large fresh water lenses in the deserts of Turkmenistan. In Fresh Water Lenses in a Desert; Academy of Sciences of the USSR: Moscow, Russia, 1963; pp. 24–94. (In Russian) [Google Scholar]
  5. Barrett, B.; Heinson, G.; Hatch, M.; Telfer, A. Geophysical methods in saline groundwater studies: Locating perched water tables and fresh-water lenses. Explor. Geophys. 2002, 33, 115–121. [Google Scholar] [CrossRef]
  6. Young, M.E.; Macumber, P.G.; Watts, M.D.; Al-Toqy, N. Electromagnetic detection of deep freshwater lenses in a hyper-arid limestone terrain. J. Appl. Geophys. 2004, 57, 43–61. [Google Scholar] [CrossRef]
  7. Houben, G.; Noell, U.; Vassolo, S.; Grissemann, C.; Geyh, M.; Stadler, S.; Dose, E.J.; Vera, S. The freshwater lens of Benjamín Aceval, Chaco, Paraguay: A terrestrial analogue of an oceanic island lens. Hydrogeol. J. 2014, 22, 1935–1952. [Google Scholar] [CrossRef]
  8. Werner, A.D.; Laattoe, T. Terrestrial freshwater lenses in stable riverine settings: Occurrence and controlling factors. Water Resour. Res. 2016, 52, 3654–3662. [Google Scholar] [CrossRef]
  9. Macumber, P.G. Lenses, plumes and wedges in the Sultanate of Oman: A challenge for groundwater management. In Developments in Water Science; Elsevier: Amsterdam, The netherlands, 2003; Volume 50, pp. 349–370. [Google Scholar]
  10. Laattoe, T.; Werner, A.D.; Woods, J.A.; Cartwright, I. Terrestrial freshwater lenses: Unexplored subterranean oases. J. Hydrol. 2017, 553, 501–507. [Google Scholar] [CrossRef] [Green Version]
  11. Cluff, C.B. Kuwait Water Harvesting System: Final Report; Water Resources Research Center, The University of Arizona: Tucson, AZ, USA, 1990. [Google Scholar]
  12. Marcella, M.P.; Eltahir, E.A. The hydroclimatology of Kuwait: Explaining the variability of rainfall at seasonal and interannual time scales. J. Hydrometeorol. 2008, 9, 1095–1105. [Google Scholar] [CrossRef]
  13. Kwarteng, A.Y.; Viswanathan, M.N.; Al-Senafy, M.N.; Rashid, T. Formation of fresh ground-water lenses in northern Kuwait. J. Arid Environ. 2000, 46, 137–155. [Google Scholar] [CrossRef]
  14. Hantush, M.S. Growth and decay of groundwater-mounds in response to uniform percolation. Water Resour. Res. 1967, 3, 227–234. [Google Scholar] [CrossRef] [Green Version]
  15. Henry, H.R. Salt intrusion into fresh-water aquifers. J. Geophys. Res. 1959, 64, 1911–1919. [Google Scholar] [CrossRef]
  16. Fadlelmawla, A.; Hadi, K.; Zouari, K.; Kulkarni, K. Hydrogeochemical investigations of recharge and subsequent salinization processes at Al-Raudhatain depression in Kuwait/Analyses hydrogéochimiques de la recharge et des processus associés de salinisation dans la dépression de Al-Raudhatain au Koweit. Hydrol. Sci. J. 2008, 53, 204–223. [Google Scholar] [CrossRef]
  17. Shanan, L. Runoff, Erosion, and the Sustainability of Ancient Irrigation Systems in the Central Negev Desert; Hebrew University of Jerusalem: Jerusalem, Israel, 2000; pp. 75–106. [Google Scholar]
  18. Milewski, A.; Sultan, M.; Al-Dousari, A.; Yan, E. Geologic and hydrologic settings for development of freshwater lenses in arid lands. Hydrol. Proces. 2014, 28, 3185–3194. [Google Scholar] [CrossRef]
  19. Lezzaik, K.; Milewski, A. A quantitative assessment of groundwater resources in the Middle East and North Africa region. Hydrogeol. J. 2018, 26, 251–266. [Google Scholar] [CrossRef]
  20. Sultan, M.; Sturchio, N.; Al Sefry, S.; Milewski, A.; Becker, R.; Nasr, I.; Sagintayev, Z. Geochemical, isotopic, and remote sensing constraints on the origin and evolution of the Rub Al Khali aquifer system, Arabian Peninsula. J. Hydrol. 2008, 356, 70–83. [Google Scholar] [CrossRef]
  21. Sultan, M.; Ahmed, M.; Sturchio, N.; Yan, E.; Milewski, A.; Becker, R.; Wahr, J.; Becker, D.; Chouinard, K. Assessment of the vulnerabilities of the Nubian Sandstone Fossil Aquifer, North Africa. In Climate Vulnerability: Understanding and Addressing Threats to Essential Resources; Academic Press: Cambridge, MA, USA, 2013; Volume 5, pp. 311–333. [Google Scholar]
  22. Lezzaik, K.; Milewski, A.; Mullen, J. The groundwater risk index: Development and application in the Middle East and North Africa region. Sci. Total Environ. 2018, 628–629, 1149–1164. [Google Scholar] [CrossRef]
  23. Fetter, C. Position of the saline water interface beneath oceanic islands. Water Resour. Res. 1972, 8, 1307–1315. [Google Scholar] [CrossRef]
  24. Vacher, H. Dupuit-Ghyben-Herzberg analysis of strip-island lenses. Geol. Soc. Am. Bull. 1988, 100, 580–591. [Google Scholar] [CrossRef]
  25. Chesnaux, R.; Allen, D. Groundwater travel times for unconfined island aquifers bounded by freshwater or seawater. Hydrogeol. J. 2008, 16, 437–445. [Google Scholar] [CrossRef]
  26. Stuyfzand, P.J.; Bruggeman, G. Analytical approximations for fresh water lenses in coastal dunes. In Proceedings of the 13th Salt Water Intrusion Meeting (SWIM), Cagliari, Italy, 5–10 June 1994; pp. 15–33. [Google Scholar]
  27. Dose, E.J.; Stoeckl, L.; Houben, G.J.; Vacher, H.L.; Vassolo, S.; Dietrich, J.; Himmelsbach, T. Experiments and modeling of freshwater lenses in layered aquifers: Steady state interface geometry. J. Hydrol. 2014, 509, 621–630. [Google Scholar] [CrossRef]
  28. Van Der Veer, P. Analytical solution for steady interface flow in a coastal aquifer involving a phreatic surface with precipitation. J. Hydrol. 1977, 34, 1–11. [Google Scholar] [CrossRef]
  29. Kacimov, A.; Obnosov, Y.V. Analytic solutions for fresh groundwater lenses floating on saline water under desert dunes: The Kunin-Van Der Veer legacy revisited. J. Hydrol. 2019, 574, 733–743. [Google Scholar] [CrossRef]
  30. Al-Weshah, R.A.; Yihdego, Y. Flow modelling of strategically vital freshwater aquifers in Kuwait. Environ. Earth Sci. 2016, 75, 1315. [Google Scholar] [CrossRef]
  31. Guo, W.; Langevin, C.D. User’s Guide to SEAWAT; A Computer Program for Simulation of Three-Dimensional Variable-Density Ground-Water Flow; USGS: Reston, VA, USA, 2002.
  32. Harbaugh, A.W.; Banta, E.R.; Hill, M.C.; McDonald, M.G. MODFLOW-2000, The U. S. Geological Survey Modular Ground-Water Model-User Guide to Modularization Concepts and the Ground-Water Flow Process; Open-File Report; USGS: Reston, VA, USA, 2000; p. 134.
  33. Langevin, C.D.; Shoemaker, W.B.; Guo, W. MODFLOW-2000, the U.S. Geological Survey Modular Ground-Water Model—Documentation of the SEAWAT-2000 Version with the Variable-Density Flow Process (VDF) and the Integrated MT3DMS Transport Process (IMT); USGS: Reston, VA, USA, 2003.
  34. Zheng, C.; Wang, P.P. MT3DMS: A Modular Three-Dimensional Multispecies Transport Model for Simulation of Advection, Dispersion, and Chemical Reactions of Contaminants in Groundwater Systems; Documentation and User’s Guide; US Army Corps of Engineers Engineer Research and Development Center: Washington, DC, USA, 1999.
  35. Langevin, C.D.; Thorne, D.T., Jr.; Dausman, A.M.; Sukop, M.C.; Guo, W. SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport; USGS: Reston, VA, USA, 2008; pp. 2328–7055.
  36. Bear, J. Groundwater Hydraulics; McGraw-Hill: New York, NY, USA, 1979. [Google Scholar]
  37. Senger, R.K.; Fogg, G.E. Stream functions and equivalent freshwater heads for modeling regional flow of variable-density groundwater: 1. Review of theory and verification. Water Resour. Res. 1990, 26, 2089–2096. [Google Scholar] [CrossRef]
  38. Holzbecher, E.O. Modeling Density-Driven flow in Porous Media: Principles, Numerics, Software; Springer Science & Business Media: Berlin/Heidelberg, Germany, 1998. [Google Scholar]
  39. Oberlander, P.L. Fluid density and gravitational variations in deep boreholes and their effect on fluid potential. Groundwater 1989, 27, 341–350. [Google Scholar] [CrossRef]
  40. Konikow, L.F.; Akhavan, M.; Langevin, C.D.; Michael, H.A.; Sawyer, A.H. Seawater circulation in sediments driven by interactions between seabed topography and fluid density. Water Resour. Res. 2013, 49, 1386–1399. [Google Scholar] [CrossRef]
  41. Moriasi, D.N.; Arnold, J.G.; Van Liew, M.W.; Bingner, R.L.; Harmel, R.D.; Veith, T.L. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 2007, 50, 885–900. [Google Scholar] [CrossRef]
  42. Doherty, J. PEST model-independent parameter estimation user manual. Watermark Numer. Comput. Brisb. Aust. 2004, 3338, 3349. [Google Scholar]
  43. Yihdego, Y.; Al-Weshah, R.A. Engineering and environmental remediation scenarios due to leakage from the Gulf War oil spill using 3-D numerical contaminant modellings. Appl. Water Sci. 2017, 7, 3707–3718. [Google Scholar] [CrossRef]
  44. Parsons Engineering and Construction Corporation. Ground-Water Resources of Kuwait; Parsons Engineering and Construction Corporation: Los Angeles, CA, USA, 1961. [Google Scholar]
  45. Kuldzhayev, N.K. Origin of the Yaskhan freshwater lens in Karakumy. Int. Geol. Rev. 1974, 16, 247–254. [Google Scholar] [CrossRef]
  46. Holt, T.; Greskowiak, J.; Seibert, S.L.; Massmann, G. Modeling the evolution of a freshwater lens under highly dynamic conditions on a currently developing Barrier Island. Geofluids 2019, 2019, 15. [Google Scholar] [CrossRef] [Green Version]
  47. Senay, Y. Groundwater Resources and Artificial Recharge in Rawdhatain Water Field; Ministry of Electricity and Water: Kuwait City, Kuwait, 1977; p. 35.
  48. Schneider, J.C.; Kruse, S.E. Assessing selected natural and anthropogenic impacts on freshwater lens morphology on small barrier Islands: Dog Island and St. George Island, Florida, USA. Hydrogeol. J. 2006, 14, 131–145. [Google Scholar] [CrossRef]
  49. Trefry, M.G.; McLaughlin, D.; Lester, D.R.; Metcalfe, G.; Johnston, C.D.; Ord, A. Stochastic relationships for periodic responses in randomly heterogeneous aquifers. Water Resour. Res. 2011, 47. [Google Scholar] [CrossRef]
  50. Van Ginkel, M.; des Tombe, B.; Olsthoorn, T.; Bakker, M. Small-scale ASR between flow barriers in a saline aquifer. Groundwater 2016, 54, 840–850. [Google Scholar] [CrossRef]
  51. Becker, D.; Sultan, M.; Milewski, A.; Becker, R.; Sauck, W.; Soliman, F.; Rashed, M.; Ahmed, M.; Yan, E.; Wagdy, A. Integrated solutions for hydrologic investigations in arid lands. Geosphere 2012, 8, 1588–1605. [Google Scholar] [CrossRef] [Green Version]
  52. Seyoum, W.M.; Milewski, A.M. Monitoring and comparison of terrestrial water storage changes in the northern high plains using GRACE and in-situ based integrated hydrologic model estimates. Adv. Water Resour. 2016, 94, 31–44. [Google Scholar] [CrossRef]
  53. Seyoum, W.M.; Milewski, A.M. Improved methods for estimating local terrestrial water dynamics from GRACE in the Northern High Plains. Adv. Water Resour. 2017, 110, 279–290. [Google Scholar] [CrossRef]
  54. Sagintayev, Z.; Sultan, M.; Khan, S.; Khan, S.; Mahmood, K.; Yan, E.; Milewski, A.; Marsala, P. A remote sensing contribution to hydrologic modelling in arid and inaccessible watersheds, Pishin Lora basin, Pakistan. Hydrol. Proces. 2012, 26, 85–99. [Google Scholar] [CrossRef]
  55. Al-Dousari, A.; Milewski, A.; Din, S.U.; Ahmed, M. Remote sensing inputs to SWAT model for groundwater recharge estimates in Kuwait. Adv. Nat. Appl. Sci. 2010, 4, 71–77. [Google Scholar]
  56. Zonn, I.S. Water resources of Turkmenistan. In The Turkmen Lake Altyn Asyr and Water Resources in Turkmenistan; Springer: Berlin/Heidelberg, Germany, 2012; pp. 59–68. [Google Scholar]
  57. Fleskens, L.; Ataev, A.; Mamedov, B.; Spaan, W. Desert water harvesting from takyr surfaces: Assessing the potential of traditional and experimental technologies in the Karakum. Land Degrad. Dev. 2007, 18, 17–39. [Google Scholar] [CrossRef]
  58. Murray, G.W. Water from the desert: Some ancient Egyptian achievements. Geogr. J. 1955, 121, 171–181. [Google Scholar] [CrossRef]
Figure 1. Conceptual diagram of an inland freshwater lens (IFL). Modified from Rotz and Milewski [1].
Figure 1. Conceptual diagram of an inland freshwater lens (IFL). Modified from Rotz and Milewski [1].
Water 12 01154 g001
Figure 2. Map of Kuwait showing the location of the Raudhatain (top) and Umm Al-Aish (bottom) IFLs highlighted in red.
Figure 2. Map of Kuwait showing the location of the Raudhatain (top) and Umm Al-Aish (bottom) IFLs highlighted in red.
Water 12 01154 g002
Figure 3. Diagram of observational data collected including water table elevation and IFL geometry (i.e., thickness, length) from previously conducted physical laboratory simulation.
Figure 3. Diagram of observational data collected including water table elevation and IFL geometry (i.e., thickness, length) from previously conducted physical laboratory simulation.
Water 12 01154 g003
Figure 4. Conceptual diagram of the numerical model with domain dimensions, boundary conditions, and observation well locations.
Figure 4. Conceptual diagram of the numerical model with domain dimensions, boundary conditions, and observation well locations.
Water 12 01154 g004
Figure 5. Physical model Simulation 2 shows IFL development and transient evolution.
Figure 5. Physical model Simulation 2 shows IFL development and transient evolution.
Water 12 01154 g005
Figure 6. Numerical model Simulation 2 of IFL development and transient evolution. The dimensions of each cross-section match the numerical model domain of 2.0 m (L) × 0.5 m (H) × 0.10 m (W).
Figure 6. Numerical model Simulation 2 of IFL development and transient evolution. The dimensions of each cross-section match the numerical model domain of 2.0 m (L) × 0.5 m (H) × 0.10 m (W).
Water 12 01154 g006
Figure 7. Simulated vs. observed water table elevation values.
Figure 7. Simulated vs. observed water table elevation values.
Water 12 01154 g007
Figure 8. Simulated vs. observed IFL thicknesses.
Figure 8. Simulated vs. observed IFL thicknesses.
Water 12 01154 g008
Figure 9. Simulated vs. observed IFL lengths.
Figure 9. Simulated vs. observed IFL lengths.
Water 12 01154 g009
Figure 10. Observed vs. simulated thicknesses with best fit exponential regression lines to represent IFL degradation for Simulation 1 (red circles), Simulation 2 (green triangles), and Simulation 3 (purple squares).
Figure 10. Observed vs. simulated thicknesses with best fit exponential regression lines to represent IFL degradation for Simulation 1 (red circles), Simulation 2 (green triangles), and Simulation 3 (purple squares).
Water 12 01154 g010
Figure 11. A second recharge pulse is applied to Simulation 2 to show the potential applications for IFL models that consider IFL transient evolution.
Figure 11. A second recharge pulse is applied to Simulation 2 to show the potential applications for IFL models that consider IFL transient evolution.
Water 12 01154 g011
Table 1. Initial model parameters.
Table 1. Initial model parameters.
Input ParameterValueDescription
Kx, Ky, Kz0.0015 m/sHydraulic conductivity
θ0.39Porosity
Ss1 × 10−5Specific storage
Sy0.34Specific yield
Dm1 × 10−10 m2/sMolecular diffusion coefficient
ΔρC0.7143Slope of fluid density to solute concentration
δL1 × 10−5 mLongitudinal dispersivity
δvl0.02Vertical/longitudinal dispersivity
Table 2. Boundary condition and concentration values.
Table 2. Boundary condition and concentration values.
Boundary ConditionPosition (x,y)Values
Constant head0, 0.4135,700 mg/L
General head2, 0.405calculated
Recharge0.79, 1.210 mg/L
EvaporationEntire domain0 m/s
Table 3. The initial recharge rate (Ri), time of recharge flux (Tq), time of observed maximum thickness (Tmax), adjusted recharge (Ra) rates, and stress period times for the numerical model simulations.
Table 3. The initial recharge rate (Ri), time of recharge flux (Tq), time of observed maximum thickness (Tmax), adjusted recharge (Ra) rates, and stress period times for the numerical model simulations.
SimulationRi (m/s)Tq (s)Tmax (s)Ra (m/s)Stress Periods Starts (1, 2, 3) (s)
14.72 × 10−5162064803.33 × 10−51620, 6480, 720,000
26.48 × 10−5129661204.72 × 10−51296, 6120, 720,000
37.86 × 10−5100854006.67 × 10−51008, 5400, 720,000
Table 4. The objective function (Φ), hydraulic conductivity estimates (Kx, Kz), and confidence intervals for the numerical simulations using the parameter estimation tool (PEST) tool.
Table 4. The objective function (Φ), hydraulic conductivity estimates (Kx, Kz), 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)
14.8 × 10−47.6 × 10−47.0 × 10−47.6 × 10−4–1.2 × 10−32.5 × 10−5–1.9 × 10−2
25.5 × 10−33.6 × 10−44.3 × 10−42.5 × 10−4–5.3 × 10−42.1 × 10−5–8.5 × 10−3
31.5 × 10−34.3 × 10−49.0 × 10−45.3 × 10−4–1.0 × 10−35.8 × 10−5–4.7 × 10−3
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 (R2), the Nash Sutcliffe Efficiency (NSE), root mean square error (RMSE), and percent bias (PBIAS).
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 (R2), the Nash Sutcliffe Efficiency (NSE), root mean square error (RMSE), and percent bias (PBIAS).
RunSlope-Intercept FormR2NSERMSEPBIAS
Water table elevationy = 1.0477x – 0.02440.890.840.007−0.7
Thicknessy = 1.0856x + 0.01330.940.650.03−16.3
Lengthy = 1.0134x – 0.00940.850.820.06−0.4
Table 6. Average groundwater velocity (vx) measurements of the numerical model simulations before, during, and after the recharge stress event.
Table 6. Average groundwater velocity (vx) measurements of the numerical model simulations before, during, and after the recharge stress event.
Runvx Before Recharge Stress Period (m/s)vx During Recharge Stress Period (m/s)vx After Recharge Stress Period (m/s)
16.88 × 10−62.40 × 10−57.04 × 10−6
23.86 × 10−66.99 × 10−53.85 × 10−6
35.64 × 10−66.48 × 10−54.84 × 10−6
Table 7. Lens degradation results between the physical model and numerical model simulations. Units have been converted from seconds to hours.
Table 7. Lens degradation results between the physical model and numerical model simulations. Units have been converted from seconds to hours.
RunPhysical Model RegressionDegradation Time (h)Numerical Model RegressionDegradation Time (h)
1y = 1.699 × 10−1 e−0.069x41y = 2.023 × 10−1 e−0.072x42
2y = 2.697 × 10−1 e−0.047x70y = 3.017 × 10−1 e−0.05x68
3y = 2.597 × 10−1 e−0.1x33y = 3.185 × 10−1 e−0.087x40

Share and Cite

MDPI and ACS Style

Rotz, R.; Milewski, A.; Rasmussen, T.C. Transient Evolution of Inland Freshwater Lenses: Comparison of Numerical and Physical Experiments. Water 2020, 12, 1154. https://doi.org/10.3390/w12041154

AMA Style

Rotz R, Milewski A, Rasmussen TC. Transient Evolution of Inland Freshwater Lenses: Comparison of Numerical and Physical Experiments. Water. 2020; 12(4):1154. https://doi.org/10.3390/w12041154

Chicago/Turabian Style

Rotz, Rachel, Adam Milewski, and Todd C Rasmussen. 2020. "Transient Evolution of Inland Freshwater Lenses: Comparison of Numerical and Physical Experiments" Water 12, no. 4: 1154. https://doi.org/10.3390/w12041154

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop