Estimating freshwater lens volume in highly permeable aquifers

: A unique freshwater lens shape observed in Tarama Island, Japan, in which hydraulic conductivity is on the order of 10 −2 m s −1 , has posed a question as to how well we can estimate the fresh groundwater volumes in extremely permeable aquifers. We applied both an analytical model and numerical simulations with various hydraulic conduc‐ tivities, including extremely permeable conditions, and compared their results. The simulation showed that, when the hydraulic conductivity was extremely high, saline groundwater existed near the coast. The analytical model overestimated the freshwater volume compared with those estimated from the numerical simulations, and the discrep‐ ancy became more significant with increasing hydraulic conductivity. These findings imply that, when hydraulic conductivity is extremely high, numerical simulations con‐ sidering density-dependent flow and dispersive mass trans‐ port processes should be applied to better assess the shapes and volumes of freshwater lenses.


INTRODUCTION
Freshwater resources on coral islands depend principally on groundwater because surface water is scarce owing to their propensity for infiltrating precipitation into highly permeable aquifers (Bailey et al., 2010;Robins, 2013). Beneath a small island, fresh groundwater exists as a "freshwater lens" where the fresh groundwater floats upon saline groundwater. Freshwater lenses are vulnerable to natural and anthropogenic disturbances, such as storm-surges and groundwater abstraction (Werner et al., 2017). To provide fundamental information for sustainable utilization of freshwater lenses, it is essential to assess their volumes.
Correspondence to: Satoshi Tajima, Graduate School of Frontier Sciences, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, The analytical model has been applied as a first-stage assessment tool (Kalf and Woolley, 2005) to sites with scarce data (Elshall et al., 2020). Numerical simulations have also been used to assess the volume of freshwater lenses because they can handle dispersive mass transport processes (Underwood et al., 1992), variably saturated flow (Elshall et al., 2020), and transient processes such as stormsurges (e.g. Chui and Terry, 2012). The analytical model suggests that hydraulic conductivity mainly controls the volume of a freshwater lens, i.e. the lens thickens with a decrease in hydraulic conductivity (e.g. Fetter, 1972;Vacher, 1988). Numerical simulations show similar changes in lens volume (Bailey et al., 2010).
Recent studies conducted in Tarama Island, Japan, showed a unique freshwater lens shape in the aquifer in which hydraulic conductivity is on the order of 10 −2 m s −1 (Figure 1d) (Ishida et al., 2011(Ishida et al., , 2013Shirahata et al., 2018). This shape is distinctive from that described by the analytical model and poses a question as to how we can properly assess freshwater resources, especially in such extremely permeable aquifers, either through the analytical model or numerical simulations. Previous studies also mentioned that the analytical model may overestimate the freshwater lens volumes compared with numerical simulations in some cases (Ketabchi et al., 2016;Oberdorfer et al., 1990).
In this study, we applied both the analytical model and numerical simulations to explore how well we can estimate fresh groundwater lens volumes in extremely permeable aquifers and how hydraulic conductivity affects the estimation. For this purpose, Tarama Island was chosen as an example because in-situ observation data are available (Ishida et al., 2011), and the freshwater lens volume was estimated by Ishida et al. (2011).
The aquifer is composed of Pleistocene limestone (the Ryukyu Group), which is underlain by Pliocene sandstone (the Shimajiri Group) (Figures 1b and 1c) (Ishida et al., 2011). The thickness of the Pleistocene limestone aquifer is approximately 50 m, in which the freshwater lens is developed. The underlying sandstone has a hydraulic conductivity of 10 −7 m s −1 , five orders of magnitude lower than that of the limestone aquifer (Ishida et al., 2011). Ishida et al. (2011) installed seven observation wells on the island and measured electric conductivity profiles inside the wells. Their study showed that groundwater with an electrical conductivity higher than 200 mS m −1 was observed near the water table. The study also showed that the bottom of the freshwater lens, defined as an isoline of an electrical conductivity of 200 mS m −1 , showed a concave-down shape near the coast (Figure 1d). These characteristics differ from other cases (e.g. Ketabchi et al., 2014;Underwood et al., 1992) and from those delineated by the analytical model.  Fetter (1972) and Vacher (1988) developed an analytical model to estimate the thickness of a freshwater lens formed in a circular island. Figure 2 shows the schematic shape of a freshwater lens. Here, b(r) [m] is the thickness of the freshwater lens as a function of the distance r [m] from the center of the island, R [m] is the radius of the island, and w [m s −1 ] is the recharge rate.

Analytical model
The lens thickness is expressed by Equation is the seawater density, and ρ f [kg m −3 ] is the freshwater density.

Numerical simulations
Numerical simulations were performed with FEFLOW (Diersch, 2014) for a conceptual axisymmetric model of Tarama Island (Figure 3). FEFLOW is a finite element based numerical simulation code that can handle variably saturated, density-dependent groundwater flow and mass transport processes (Diersch, 2014). In this study, the van Genuchten-Mualem model was employed to represent the unsaturated zone properties (Mualem, 1976;van Genuchten, 1980).
The radius of the model was set to 2.5 km, which is the average radius of Tarama Island. The land surface elevation was set at 10 m a.m.s.l. from the island's center to a radial distance of 2.4 km. The bottom of the model was set at −50 m a.m.s.l., which is equivalent to the bottom of the limestone aquifer.
The boundary conditions applied for the model ( Figure  3) are summarised as below. For fluid flow processes, a constant influx equivalent to a recharge rate (w) was assigned to the land surface above 1 m a.m.s.l., and a constant seawater head of 0 m was applied to the seafloor. For mass transport processes, a constant mass concentration of 0 kg m −3 was assigned to the land surface above 1 m a.m.s.l., and a constant mass concentration of 35 kg m −3 Figure 2. Schematic of a freshwater lens formed in a circular island FRESHWATER LENS VOLUME: HIGH-PERMEABILITY AQUIFERS was applied to the seafloor when seawater flowed into the model domain. When groundwater was discharged to the sea at the seafloor, its mass concentration was not constrained. A seepage face was specified at the land surface below 1 m a.m.s.l., where water was allowed to flow out of the model domain without any constraints on the mass concentration. Different heights of the seepage face were tested (1 to 4 m), and it was confirmed that the simulation results were not sensitive to the changes within the ranges tested. Along the central axis and the bottom boundary, a no-flow boundary condition was assigned.  The initial mass concentrations above and below the freshwater-seawater interface calculated from the analytical model (Equation 1) were set to 0 kg m −3 and 35 kg m −3 , respectively. The domain above the sea level was initially considered unsaturated, whereas below the sea level fully saturated.
To prevent numerical oscillations, the meshes were generated to satisfy the mesh Péclet number less than 2 (Daus et al., 1985). The mesh sizes above the sea level ranged between 0.12 and 0.95 m and below the sea level between 0.73 and 3.64 m. The initial time step length was set to 10 −5 d, and the maximum time step length was 2 d. The maximum growth factor of each time step length was restricted to 1.2. The simulations were run for 100 years (36500 d), and the final outputs were considered to have reached a quasi-steady state as the temporal changes in the salinity distributions became sufficiently small. Table I summarizes the parameters used for the baseline numerical simulation. To investigate how hydraulic conductivity affects the lens shapes, numerical simulations with different hydraulic conductivities (2.0 × 10 −4 to 6.3 × 10 −2 m s −1 ) were conducted. In the case where the hydraulic conductivity was smaller than 3.6 × 10 −3 m s −1 , the model domain was expanded downwards and seawards to prevent the domain boundaries from touching the lens bottoms.

Calculation of freshwater volumes
The volumes of the freshwater lens were calculated both from the analytical solution and the numerical simulations. From the analytical solution, the volume of the freshwater (V [m 3 ]) within r [m] from the island's center was calculated as follows: where ε [-] is the porosity. In the numerical simulations, the bottoms of the freshwater lenses were defined as isochlors (mass concentration Seawater density (ρ s ) 1.024 × 10 3 kg m −3 isolines) of 1.36 kg m −3 . These volumes were compared with those based on the in-situ measurements by Ishida et al. (2011). They defined the bottom of the island's freshwater lens as 200 mS m −1 , which is equivalent to a mass concentration of 1.36 kg m −3 according to the conversion equation (Jiao and Post, 2019): where c is the mass concentration [kg m −3 ] and SC is the specific conductance [mS m −1 ].

Simulated salinity distributions
The baseline simulation showed the existence of a wide and thick transition zone between freshwater and seawater nearby the coast (Figure 4c). The Darcy flux was almost horizontal and increased closer to the coast ( Figure S1). The horizontal width of the transition zone increased with increasing hydraulic conductivity (e.g. Figures 4a and 4c).
When the hydraulic conductivities (K) were 2.0 × 10 −2 m s −1 (the baseline value) and 6.3 × 10 −3 m s −1 , the simulation showed that groundwater with mass concentration higher than 1.36 kg m −3 existed at the water table and in the unsaturated zone near the coast (Figures 4b and 4c), and that the freshwater zones were located in inland positions. When K = 2.0 × 10 −2 m s −1 , a larger portion of the island was subjected to the unsaturated zone salinization compared with the case where K = 6.3 × 10 −3 m s −1 . When K = 6.3 × 10 −4 m s −1 , no saline water was observed at the water table, and the freshwater zone existed near the coast ( Figure  4a). These simulation results suggest that, when hydraulic conductivity is extremely high, no fresh groundwater exists near the coast. These results are consistent with the observation by Ishida et al. (2011) but different from the freshwater distribution estimated from the analytical model (Figure 4). The difference was caused by the assumption in the analytical solution that the entire subsurface domain above the Ghijben-Herzberg interface was filled with freshwater (Fetter, 1972;Vacher, 1988), which resulted in the existence of fresh groundwater at the coast.

Freshwater volumes
When the hydraulic conductivity was lower than 6.3 × 10 −4 m s −1 , the freshwater volumes calculated from the analytical model (Equation 1) were almost identical to those from the numerical simulations. In contrast, when the hydraulic conductivity was higher, the freshwater volumes calculated from the analytical model deviated from those estimated from the numerical simulations, and the differences became more significant when the hydraulic conductivity became larger ( Figure 5). In the baseline case (K = 2.0 × 10 −2 m s −1 ), the analytical model estimated the lens volume around six times as large as the numerical simulation, whereas the freshwater volume calculated from the in-situ measurements (Ishida et al., 2011) was within the precipitation-related uncertainty of that calculated from the numerical simulations ( Figure 5).
When hydraulic conductivity is extremely high, dispersive mass transport becomes significant because of the  In the case with a hydraulic conductivity of 6.3 × 10 −4 m s −1 , the model domain was expanded downwards and seawards to prevent the lens bottom and the transition zone from touching the boundaries of the model domain increase of average linear velocity through the porous medium. In particular, upward mass transport is enhanced by the transverse dispersion within the capillary fringe (Tajima et al., 2021). The analytical model does not consider the dispersive mass transport processes, whereas the numerical simulations can handle these processes. This difference results in the notable difference between the analytical model and the numerical simulations, especially for the extremely permeable aquifers.
Even when hydraulic conductivity is extremely high, numerical simulations can obtain a freshwater volume close to the analytical model by adjusting the dispersivity values to be extremely small. However, such small dispersivities may not be reasonable considering the scales of interest and the possible heterogeneous nature of the limestone aquifer. Our findings imply that numerical simulations considering density-driven flow and dispersive mass transport processes with a proper parameterization of dispersivity should be applied to extremely permeable aquifers to assess the freshwater lenses more accurately.

Future research
The numerical simulations employed in this study included several simplifications, e.g. the homogeneous material properties such as hydraulic conductivity, and the constant sea level as the boundary condition. Future research will investigate these factors to represent actual situations more realistically.

CONCLUSION
This study assessed the volume of a freshwater lens in an extremely permeable aquifer by comparing the in-situ measurements, the analytical model, and the numerical simulations. The key findings of the study are: 1) When the hydraulic conductivity was extremely high, no freshwater resources were available near the coast, which is different from the concept of a freshwater lens represented by the analytical model.
2) The analytical model overestimated the freshwater volumes compared with the numerical simulations when the hydraulic conductivity was extremely high, and the discrepancy became more significant with increasing hydraulic conductivity.
These findings imply that, when hydraulic conductivity is extremely high, dispersive mass transport process should not be ignored. Therefore, it is suggested that numerical simulations considering density-dependent flow and dispersive mass transport processes should be applied to extremely permeable aquifers to assess the freshwater lenses more accurately. As indicated by this study, numerical simulations can provide us a detailed understanding on the dynamics of fresh/saline groundwater and the formation mechanisms of the unique salinity distribution in a coastal aquifer.