Mechanisms of Non‐Fresh Groundwater Presence at Water Tables in Highly Permeable Coastal Aquifers

Coastal aquifers with high hydraulic conductivities on the order of 10−2 m s−1 have unconventional salinity distributions with the presence of non‐fresh groundwater at the water table over a wide swath near the coast. This study aims to unravel the mechanisms underlying the phenomenon via numerical simulations for variably saturated, density‐driven flow and solute transport in porous media. The simulation results indicate that the existence of non‐fresh groundwater at the water table is attributed to the upward mass flux in the saturated zone near the coast, which transports solute from deeper groundwater toward the water table. With high hydraulic conductivity, the upward mass flux becomes prominent at shallower elevations because of the high Darcy flux and the shallow saline groundwater. The upward mass flux has two main drivers, upward advection by the upward flow component and transverse dispersion by the seaward flow component. The advective mass flux dominates over the transverse dispersion in the deep part of the saturated zone where only groundwater with sea water salinity exists. In contrast, the transverse dispersion becomes more pronounced than the upward advection in the shallow saturated zone just beneath the water table and in the unsaturated zone immediately above the water table. Our findings help interpret the unconventional salinity distributions observed and elucidate the unique dynamics of groundwater flow and solute transport in highly permeable coastal aquifers.


Introduction
Coastal aquifers are characterized by the occurrence of interactions between fresh groundwater derived from terrestrial recharge and saline groundwater originating from the sea (Cooper et al. 1964;Jiao and Post 2019).Groundwater at the water table is assumed to be fresh as described by the Ghijben-Herzberg analytical model (Drabbe and Badon Ghijben 1889;Herzberg 1901;Fetter 1972).In contrast, observations on islands with highly permeable aquifers indicate that non-fresh groundwater exists at the water table (Ghassemi et al. 1996;Ishida et al. 2011; The Japanese Institute of Irrigation and Drainage [JIID] 2014).Field measurements on Tarama Island, Japan, showed that non-fresh groundwater with electrical conductivity higher than 2 mS cm −1 exists at the water table at approximately 1.1 km from the coast (Ishida et al. 2011;JIID 2014) (Appendix S1).Ghassemi et al. (1996) estimated the subsurface salinity distribution of Nauru Island, the Republic of Nauru, by electrical conductivity measurements of groundwater in boreholes and geoelectrical soundings.Their results indicate that non-fresh groundwater (electrical conductivity >2.2 mS cm −1 ) exists at the water table at approximately 200 m from the coast.The aquifers of these two islands are characterized by their high hydraulic conductivity (K ) values on the order of 10 −2 m s −1 (Ghassemi et al. 1996;Shirahata et al. 2018), one to two orders of magnitude higher than those reported for Holocene aquifers of atolls (Underwood et al. 1992;Ketabchi et al. 2014;Werner et al. 2017).Then, what physical processes contribute to the phenomenon observed in the high-K coastal aquifers that non-fresh groundwater exists at the water table?
Dynamic environmental processes can at least partly explain this phenomenon.For example, tidal fluctuations result in the presence of non-fresh groundwater at the water table in unconfined coastal aquifers (Werner and Lockington 2006).Badaruddin et al. (2015) indicated that water table salinization occurs with a decline in water table elevation by a landward hydraulic head gradient, called "active sea water intrusion" as defined by Werner (2017).In addition, dynamic processes associated with infrequent events such as storms and tsunamis can contaminate coastal zones with salt spray because of strong winds (Du and Hesp 2020), wave overwash owing to storm surge (e.g., Terry and Falkland 2010;Bailey et al. 2014), and groundwater flooding accompanying storm surge (Tajima et al. 2023).Following these processes during storms, the land surface is inundated with sea water, and the puddled sea water percolates into the subsurface, inducing the long-term degradation of fresh groundwater resources (e.g., Chui and Terry 2012;Bailey and Jenson 2014;Holding and Allen 2015).Tsunamis also result in a similar sea water inundation of the land surface (e.g., Alsumaiei and Bailey 2018;Liu and Tokunaga 2019).These findings imply that infrequent events may also contribute to the presence of non-fresh groundwater at the water table.
However, numerical simulations by Tajima et al. (2022) imply that, even in a quasi-steady state without the transient forces described above, non-fresh groundwater can exist at the water table if the aquifer is highly permeable.Badaruddin et al. (2015) also indicated that high-K aquifers may have non-fresh groundwater at the water table in a quasi-steady state.These authors attribute the observed characteristic salinity distributions to dispersive mass transport processes (Badaruddin et al. 2015;Tajima et al. 2022), yet quantitative analyses are required to further understand the dynamics of groundwater and solute in the subsurface.This study hence aims to elucidate the quantitative mechanisms underlying the existence of non-fresh groundwater at water tables in high-K coastal aquifers.To this end, we run numerical simulations for a conceptual highly permeable coastal aquifer with hydrogeological settings imitating those of Tarama Island.In the simulations, we test various scenarios of aquifer hydraulic conductivity values.We then quantitatively analyze the mass transport processes in the simulation results by calculating various components of the vertical mass flux.

Method
We consider variably saturated, density-dependent groundwater flow and mass transport in a two-dimensional axisymmetric domain.FEFLOW (Diersch 2014), a numerical simulation code for finite-element analysis, was used for the numerical simulations.The governing equations in FEFLOW are summarized by Diersch (2014).Information on the unsaturated flow models applied in these simulations is provided in Appendix S2.

Model Geometry
The model domain (Figure 1) extended 2.55 km from the island's center, including both the land and sea domains.The island's radius was 2.5 km.A constant land surface elevation of 10 m above mean sea level (a.m.s.l.) was assumed from the island's center to a radial distance of 2.4 km.Between 2.4 and 2.55 km from the center, the land surface had a constant slope of 1:10.The bottom of the model domain was set to 100 m below sea level.The model domain was built upon the spatial scale of Tarama Island, which can be approximated as a circular island (Appendix S1).Note that the objective of this study is to investigate general salinity distributions in high-K aquifers, rather than reproducing the site-specific subsurface salinity distribution in the aquifer of Tarama Island.The entire aquifer was assumed to be unconfined with homogeneous physical properties imitating the Pleistocene limestone aquifer of Tarama Island (Appendix S1), which has a high K of 10 −2 m s −1 (Shirahata et al. 2018) because of karstification (Yoshimoto et al. 2016).

Initial and Boundary Conditions
A constant sea water head of 0 m was set along the seafloor (Figure 1).In FEFLOW, the sea water head h s is converted to an equivalent freshwater head.Along the seafloor, a mass concentration of 35 kg m −3 was applied when sea water flowed into the model domain.
Here, we also considered the effect of outflow from the model domain into the sea without constraints on the mass concentration.To the land surface above 1 m a.m.s.l., we assigned a constant recharge (w ) and a mass concentration of 0 kg m −3 .A seepage face was assumed along the land surface extending from sea level to 1 m a.m.s.l.Along the seepage face, water was discharged out of the model domain with its mass concentration, and no inflow into the model domain was considered.We tested various values of seepage face height (1, 2, and 4 m), but observed no notable differences in the simulation results.A no-flow boundary condition was applied to the bottom and the axis of symmetry.
Initial conditions were h = 0 m a.m.s.l. and c = 0 kg m −3 for the domain above the Ghijben-Herzberg interface (Fetter 1972).Below the interface, sea water head h s = 0 m a.m.s.l. and c = 35 kg m −3 were taken as initial conditions.

Spatial and Temporal Discretization
The model domain was spatially discretized into 43,866 nodes and 86,506 triangular elements.The maximum mesh diameter is 3.83 m, which satisfies an upper limit in the mesh Péclet number of 2 to prevent numerical oscillations (Daus et al. 1985).
The simulations employed the predictor-corrector method with the forward/backward Euler scheme for time integration (Diersch 2014).The initial time step length is 0.864 s (10 −5 days), and the maximum length is 172,800 s (2 days).We restricted the maximum ratio between the two consecutive time steps to 1.2.The simulations were run for 50 years (18,250 days).We considered that the output of the final time step had reached a quasi-steady state by confirming that temporal changes in the head and salinity distributions became smaller than 0.1% (e.g., Ajami et al. 2014).

Parameters
We selected parameter values for the baseline simulation (Table 1), which imitate the hydrogeological conditions of Tarama Island (see Appendix S1 for details).Ishida et al. (2013) estimated the ratio of the annual recharge to the rainfall to be approximately 42%.As the average annual precipitation is 1987 mm a −1 (Japan Meteorological Agency 2021), a constant recharge rate (w ) of 2.6 × 10 −8 m s −1 (835 mm a −1 ) was used for each simulation case.We also ran additional simulations with w = 1.3 × 10 −8 and 5.2 × 10 −8 m s −1 because not only K but also w controls the water table elevation, which may influence the presence of non-fresh water at the water table.K was assumed to be anisotropic: . Hereafter, the "horizontal hydraulic conductivity" is simply referred to as the "hydraulic conductivity" (K ) unless explicitly mentioned.A hydraulic conductivity of 2.0 × 10 −2 m s −1 was adopted as the baseline value following the estimation by tidal analysis (K = 1.1 × 10 −2 to 2.8 × 10 −2 m s −1 ) (Shirahata et al. 2018).We designed various scenarios with K = 8.0 × 10 −4 , 1.8 × 10 −3 , 4.0 × 10 −3 , 8.9 × 10 −3 , 2.0 × 10 −2 (baseline), and 4.5 × 10 −2 m s −1 .These values fall within the range typical for atolls (Werner et al. 2017).
In addition, we performed sensitivity analyses for the K anisotropy and dispersivity (Appendix S3).The results indicate that the K anisotropy and longitudinal dispersivity (α L ) have minimal impacts on the presence or absence of non-fresh groundwater at the water table, whereas the transverse dispersivity (α T ) can affect it.The latter corresponds to our finding that the transverse dispersion is a key driver of the upward mass flux, which results in the existence of non-fresh water at the water table (cf.Section 3).

Vertical Mass Flux Analysis
To investigate mass transport processes, we calculate the vertical mass flux in the simulation results.In a twodimensional axisymmetric coordinate system, the vertical mass flux (I ) (M L −2 T −1 ) is defined as where c is the mass concentration (M L −3 ).The dispersion-diffusion tensor (D) in Equation 1 is defined as follows (Eeman et al. 2011): where q r and q z are the radial and vertical components of the Darcy flux (L T −1 ), respectively, α L and α T are the longitudinal and transverse dispersivities (L), respectively, and ε is the porosity.
In the transition zone between the fresh and saline water, the horizontal gradient of the mass concentration is negligibly small compared with the vertical gradient (Eeman et al. 2011).Thus, ∂c ∂r = 0 is assumed, and the NGWA.org S.  following equation is formulated from Equation 2: where I z (M L −2 T −1 ) is the vertical component of mass flux, referred to as the "vertical mass flux" (positive is upward).Note that "vertical" here means "along the z -direction," not "orthogonal to the Ghijben-Herzberg interface." The vertical mass flux is divided into four components: advection (I zA ; M L −2 T −1 ; Equation 4), longitudinal dispersion with the vertical flow component (I zL ; M L −2 T −1 ; Equation 5), transverse dispersion with the horizontal flow component (I zT ; M L −2 T −1 ; Equation 6), and molecular diffusion (I zM ; M L −2 T −1 ; Equation 7).I zM is ignored in this study because it is negligibly small compared with other components.
The vertical gradient of mass concentration is approximated by where c ij +1 , and c ij −1 (M L −3 ) are the mass concentrations of the upper (ij + 1) and lower (ij − 1) adjacent nodes of the node ij , respectively, and z (L) is the distance between the nodes ij + 1 and ij − 1.

Quasi-Steady-State Salinity Distributions
The baseline simulation result with K = 2.0 × 10 −2 m s −1 shows that non-fresh groundwater exists at the water table within approximately 1.3 km from the coast (Figures 2a and 2b).Non-fresh water also fills the lower parts of the unsaturated zone near the coast.In the case with K = 8.0 × 10 −4 m s −1 , only fresh groundwater exists at the water table across the entire aquifer (Figures 2c and 2d). Figure 3 demonstrates that the lower the ratio w /K , the larger the horizontal extent of non-fresh groundwater at the water table (L Sal /R).Here, L Sal /R is defined as the width of the water table where groundwater with a mass concentration higher than 1.36 kg m −3 is observed (Ishida et al. 2011;JIID 2014), relative to the island's radius R = 2.5 km.
Tidal fluctuations can increase the horizontal extent of non-fresh groundwater because they enhance the mixing between fresh and saline water (Oberdorfer et al. 1990) and cause water table salinization in unconfined aquifers (Werner and Lockington 2006).However, our results indicate that even without tidal fluctuations, saline groundwater exists at the water table if the hydraulic conductivity is high (Figures 2 and 3).Therefore, not only tidal fluctuations but also the high hydraulic conductivity of the aquifer explains the presence  of non-fresh groundwater at the water table observed on Tarama and Nauru Islands (Ghassemi et al. 1996;Ishida et al. 2011;JIID 2014).

Mechanisms behind the Existence of Non-fresh Groundwater at the Water Table
In this subsection, we first compare the overall trends in the vertical mass flux, calculated from the baseline simulation result, at r = 2.3 km with those at r = 0.5 and 1.0 km.Recall that non-fresh groundwater exists at the water table at r = 2.3 km, whereas only freshwater exists at the water table at r = 0.5 and 1.0 km (Figure 2a).In addition, we discuss the difference in the dominant mass transport process that contributes to the vertical mass flux at different depths at r = 2.3 km.Then, we compare the vertical profiles of the vertical mass flux at r = 2.3 km with the baseline K to that with a lower K of 8.0 × 10 −4 m s −1 .
At r = 2.3 km, the vertical mass flux by advection (I zA ) and transverse dispersion (I zT ) are higher than those at r = 0.5 and 1.0 km, with peak values approximately six times larger (Figures 4a, 4f, and 4k).The high I zA at r = 2.3 km is because the vertical Darcy flux (q z ) shows larger absolute values at r = 2.3 km than those at r = 0.5 and 1.0 km, and q z is always upward (positive) in the saturated zone at r = 2.3 km, whereas it becomes downward (negative) at some elevations at r = 0.5 and 1.0 km (Figures 4e, 4j, and 4o).In addition, the high I zA can also be explained by the higher mass concentration (c) at r = 2.3 km than that at r = 0.5 and 1.0 km, especially above approximately 5 m a.m.s.l.(Figures 4b, 4g, and 4l).The high I zT at r = 2.3 km is attributed to the higher horizontal Darcy flux (q r ) at r = 2.3 km compared with that at r = 0.5 and 1.0 km (Figures 4d, 4i, and 4n).
From these results, we attribute the presence of non-fresh groundwater at the water table under high-K conditions (Figures 2a and 4a) to the high upward mass flux in the saturated zone near the coast, mainly driven by advection with the upward flow component and transverse dispersion with the seaward flow component.
From here, we focus on the vertical profile of the vertical mass flux at r = 2.3 km (Figures 4 and 5).The vertical mass flux by advection (I zA ) dominates over the transverse dispersion (I zT ) in the deep part of the saturated zone (below approximately −15 m a.m.s.l) (Figure 4k) that is entirely filled with groundwater with the same salinity as sea water (Figure 4l).In contrast, I zT becomes larger than I zA in the shallow saturated zone just beneath the water table and the unsaturated zone immediately above the water table (below approximately −1 m a.m.s.l.) (Figure 5a).The high I zA in the deep saturated zone is attributed to the high mass concentration (c) and vertical upward Darcy flux (q z ) (Figures 4l and 4o).The dominance of I zT in the shallow saturated zone is due to the high q r (Figure 5d), and the peak in ∂c/∂z observed in the unsaturated zone just above the water table (Figure 5c) explains the dominance of I zT in this zone.Above approximately −1 m a.m.s.l., q z becomes negative (Figure 5e) and hence I zA turns downward (Figure 5a).However, the net vertical mass flux remains upward because I zT has larger absolute values than I zA .
Figure 4 indicates high positive (upward) I zA at r = 2.3 km in the high-K case, which might appear counterintuitive against the simulation results showing notable downward flow due to recharge (Figure 4) and the flowlines that approximately parallel the isochlors (Figure 2b).However, the results displayed in Figure 4 hold reasonable from two aspects.First, in contrast to those in the lower-K case, the flowlines in the high-K case are almost horizontal immediately below the water table near the coast.Second, the flow along the isochlors has a non-zero upward component, although the horizontal component is dominant.
Compared with the case with K = 8.0 × 10 −4 m s −1 , I zA and I zT in the baseline case (K = 2.0 × 10 −2 m s −1 ) peak at shallower depths, and their maximum values are more than three times larger (Figures 4k and 6a).This can be explained by two mechanisms.First, in the baseline case with K = 2.0 × 10 −2 m s −1 , q z and q r are higher than those in the case with a lower K of 8.0 × 10 −4 m s −1 (Figures 4n and 4o; 6d and 6e), resulting in the higher I zA and I zT in the baseline case.Second, the peaks of q z , q r , I zA , and I zT in the baseline case are observed at shallower elevations than those in the low-K case .This observation is corroborated by the finding that the Ghijben-Herzberg interface becomes shallower as the ratio w /K decreases (Fetter 1972)  at shallow elevations because of the high Darcy flux and shallow saline groundwater, which explains the presence of non-fresh groundwater at the water table near the coast observed from our baseline simulation (Figures 2a and  2b) and the field measurements (Ghassemi et al. 1996;Ishida et al. 2011;JIID 2014).

Implications for Submarine Groundwater Discharge
Figure 7 shows that the mass concentration of discharged groundwater at the coast at z = 0 m a.m.s.l.increases with a decrease in w /K .In the baseline simulation (K = 2.0 × 10 −2 m s −1 ), all the terrestrially derived freshwater is mixed with the saline groundwater and discharged into the sea as non-fresh water (Figures 2b  and 7).In contrast, the quintessential submarine groundwater discharge includes fresh groundwater discharge from the seafloor (Taniguchi et al. 2002;Burnett et al. 2003) as observed in our simulation result with K = 8.0 × 10 −4 m s −1 (Figures 2d and 7).These results imply that no submarine fresh groundwater discharge occurs under high-K or low-w settings.

Model Validity and Limitation
Although the hydraulic conductivity was set extremely high in our simulations, Darcy's Law is still valid because the Reynolds number is sufficiently small.Assuming that the characteristic length is L = μK ερ f g (Collins 1976), we write the Reynolds number Re as From the maximum Darcy flux ||q|| < 10 −5 m s −1 obtained from the simulations (Figure 4), we obtain Re < 1.1 × 10 −3 .This value is sufficiently smaller than the limit of the Reynolds number (1-10) for applying Darcy's Law (Bear 1972).
A limitation lies in the Fickian solute transport model used in this study.If groundwater flows toward a high-concentration zone, the Fickian model calculates dispersive mass flux to be opposite to the flow direction, and the solute spreads upstream (Konikow 2011).This phenomenon called "upstream dispersion" may be unrealistic (Irvine et al. 2021) because no macroscale flow occurs in the upstream direction (Konikow 2011).Upstream dispersion may have occurred in our baseline simulation in the unsaturated zone, where upward longitudinal dispersion is observed (Figure 5a).Estimating and mitigating the impacts of upstream dispersion on the simulated salinity distributions may facilitate an accurate discussion of the salinity distribution in the unsaturated zone.

Conclusions
This study aimed to unravel the mechanisms underlying the presence of non-fresh groundwater at the water table near the coast in highly permeable coastal aquifers.We ran numerical simulations considering various scenarios for hydrological parameters and analyzed the mass transport processes.
The results show that, under high-K conditions, prominent upward mass flux occurs at shallow elevations.This is because high Darcy flux occurs with the high K , and saline groundwater exists at shallow depths due to the low w /K .The upward mass flux transports the solute toward the water table and contributes to the existence of non-fresh groundwater at the water table.The dominant processes of the upward vertical mass flux are advection with the upward flow component and transverse dispersion with the seaward flow component.In the deep part of the saturated zone where only groundwater with a sea water salinity exists, the vertical mass flux by advection contributes greater than the transverse dispersion.In contrast, the transverse dispersion becomes more pronounced than the advection in the shallow saturated zone just beneath the water table and the unsaturated zone immediately above the water table.These findings help interpret the unconventional salinity distributions observed and elucidate the unique dynamics of groundwater flow and solute transport in highly permeable coastal aquifers.
Future research could estimate and mitigate the impacts of upstream dispersion, inherent to the Fickian model for solute transport, on the simulated salinity distributions.

Figure 1 .
Figure 1.Model domain with boundary conditions.The part figure (b) depicts enlarged view of coastal area in (a).

Figure 3 .
Figure 3. Relationship between w /K and L Sal /R (width of water table where non-fresh groundwater exists, relative to island's radius, R = 2.5 km).

Figure 5 .
Figure 5. Enlarged views of Figures 4k-4o (r = 2.3 km) around the coast.(a-f) depict vertical mass flux by each component (I z A , I z L , and I z T ), c, ∂c/∂z , q r , q z , and degree of saturation (s), respectively.Recharge rate is w = 2.6 × 10 −8 m s −1 (baseline).Blue dashed lines indicate water table.Positive mass flux means upward.

Figure 6 .
Figure 6.Vertical profiles of parameters analyzed with K = 8.0 × 10 −4 m s −1 at r = 2.3 km.(a-d) depict vertical mass flux by each component (I z A , I z L , and I z T ), c, ∂c/∂z , q r , q z , respectively.Recharge rate is w = 2.6 × 10 −8 m s −1 (baseline).Blue dashed lines indicate water table.Positive mass flux means upward.

Table 1 Simulation Parameters and Their Tested Ranges Parameter Baseline Value Tested Range Unit Source
. Our results suggest that, if the hydraulic conductivity is high, prominent upward mass flux occurs