2.1 Indicators Used in Analytical Analysis

The mathematical expression of the sharp interface model and related analytical solutions are based on the Ghyben- Herzberg principle. When computing the SWI process, mixing is commonly ignored, and a sharp interface between the saltwater and freshwater zones is assumed. Most analytical models adopt a simplified, sharp interface approximation ^{(Strack, 1976)}. The equation for the sharp interface model is a single-potential analysis and was derived from a study by Custodio and Bruggeman ^{(Custodio and Bruggeman, 1987)}. In early studies, Tamai and Shima ^{(Tamai and Shima, 1967)} developed an approximate solution based on the Duipuit assumption to predict the transient movement of an intrusion wedge toe. The authors compared the mathematically predicted results to the experimental results and observed a good match when the experiment was conducted using coarse sand. ^{Vappicha and Nagaraja (1976)} developed a similar model ^{(Bear and Dagan, 1964)}, which adopted indicators such as the rate of movement of the interface and the length of intrusion. They compared the mathematically predicted results with the experimental results and observed a good match when the experiment was conducted with coarse sand. ^{Chesnaux and Allen (2008)} used closed-form analytical solutions for circular and strip islands, applied the Ghyben-Herzberg relation, and compared the location of the water table and the time of groundwater transport. ^{Werner and Simmons (2009)} used analytical models to classify the extent and nature of SWI in aquifers with immediate SLR according to the type of freshwater inflow. To evaluate the SWI extent, they adapted an equation representing the salt wedge toe, originally developed by ^{Strack (1976)} and ^{Custodio (1987)}. The analytical index, considering the representative parameters of the Nile Delta region, showed results consistent with the numerical solution of ^{Sherif and Singh (1999)}. Nevertheless, the authors recommended using it to manage coastal aquifers impacted by SLR separately rather than proposing a replacement with a numerical solution. ^{Werner et al. (2011)} improved the existing analytical solutions to devise new SWI vulnerability indicators using two partial derivatives: 1) the rate of change in the wedge toe and 2) the rate of change in seawater volume when the toe of the wedge and the volume of water are expressed in the analytical solution. The stressors affecting the vulnerability of the SWI include SLR, changes in recharge due to climate change, and changes in seaward discharge. Based on these SWI vulnerability indicators, a national inventory of SWI vulnerability for Australia was introduced based on 28 coastal aquifers in Australia ^{(Morgan and Werner, 2015)}. As a useful MODFLOW family tool, the SWI package was developed by ^{Bakker et al. (2013)}. This package is based on sharp interface approximations and is more suitable for modeling systems with negligible dispersion effects.

Recently, researchers have begun to consider uncertainties when studying the extent of SWIs. ^{Ferguson and Gleeson (2012)} investigated the extent of the SWI based on the impact of SLR and saltwater inundation. They concluded that saltwater inundation is more important than SLR inundation because of the topographic gradient in low-lying coastal areas. By taking into account the range of one standard deviation of hydraulic conductivity owing to various aquifer types, from coarse-grained classic rocks to coarse-grained unconsolidated sediment, the study included a statistical indicator, the uncertainty of the SWI. Representative analytically modeled SWI studies are summarized in Table 1.

2.2 Indicators Used in Conceptualized Numerical Models

The extent of the SWI and movement of the mixing zone depend on the geography, hydrogeological properties of the aquifer, and meteorological parameters, such as precipitation. Rectangular, two-dimensional, and three-dimensional precise numerical simulations are used for density-dependent groundwater flows. The advantage of more sophisticated numerical models is that mixing effects due to dispersion can be simulated across the interface. Numerical models are useful for modeling field-scale problems because they include topographic complexity and geological variability ^{(Dausman and Langevin, 2005; Oude Essink et al., 2010; Praveenaa et al., 2010; Vandenbohede and Lebbe, 2002)}. A particular set of benchmark problems is solved with tried-and-true analytical, laboratory, and numerical solutions in order to validate a numerical model. The Henry problem ^{(Henry, 1964)} is the only density coupling numerical example related to the saltwater-freshwater interface and SWI for which a semi-analytical solution exists in a two-dimensional flow. This solution defines the location and shape of the freshwater boundary, showing that freshwater flow moves toward the ocean boundary. During the data validation process, the isochore, which includes linearized points corresponding to 50 % concentration of saltwater, was used as the main indicator defining the saltwater-freshwater interface. After being introduced as a benchmark solution, the model design and numerical schemes were comprehensively evaluated and modified by researchers. These studies were conducted by verifying the model parameters and modifying a few according to Henry's solution. ^{Pinder and Cooper (1970)}, ^{Lee and Cheng (1974)}, ^{Segol and Pinder (1976)}, Frind ^{(Frind, 1982)}, ^{Huyakorn et al. (1987)}, and ^{Voss and Souza (1987)} developed several modified solutions. ^{Croucher and O'Sullivan (1995)} concluded that Henry's original interpretation may have been slightly erroneous. ^{Segol (1994)} reexamined the errors related to Henry's solution and suggested a few corrections. The proposed correction is a 183-term truncation scheme. These solutions to the Henry problem have mainly investigated the steady-state isochore distribution. A criterion for determining the extent of SWI is the locations of additional isochores, like 25 %, 75 %, and 50 % isochores, because of the mixing zone in the domain of the numerically simulated Henry problem. ^{Simpson and Clement (2003)} adopted a novel transient indicator for the Henry problem: the transient position of the toe. Simpson and Clement proposed a 293-term truncation scheme ^{(Simpson and Clement, 2004)}. When Segol ^{(Segol, 1994)} recalculated Henry’s solution using additional terms, the new solution differed slightly from the original. With this new solution, ^{Segol (1994)} showed that numerical codes could reproduce the correct solutions for the Henry problem. The coupling of the flow and transport based on the density difference owing to the high concentration was a key factor in the numerical analysis. ^{Simpson and Clement (2003)} compared the isochores of domains for coupled versus uncoupled cases and reported that the coupling process effect of the Henry problem was not evident when the default parameter values were given. Based on these findings, Simpson and Clement ^{(2003; 2004)} discussed a modified Henry problem by reducing the dimensionless freshwater flux among the input parameters. ^{Abarca et al. (2007)} investigated the heterogeneity of the Henry problem by considering a dimensionless input parameter, and the results illustrated a dimensionless saltwater flux for dispersive and diffusive cases. ^{Kerrou and Jaouher (2010)} investigated the two-dimensional and three-dimensional Henry problems by following the dimensionless parameters of ^{Abarca et al. (2007)} to avoid scale-dependent analysis. The dimensionless penetration length of the saltwater wedge and the dimensionless width of the mixing zone described the geometry of the saltwater-freshwater interface. The amount of saltwater entering the aquifer system is described by the dimensionless saltwater inflow flux. The results showed that in two dimensions, anisotropy significantly affects the permeability coefficient; however, in three dimensions, the degree of diffusion through the dispersion coefficient changes, changing the length of the SWI. In addition to the dimensionless salt wedge toe, width of the mixing zone, and saltwater flux, recent studies of the Henry problem considered additional dimensionless indicators ^{(Fahs et al., 2016; Fahs et al., 2018)}.

A laboratory-scale experimental setup was used to analyze the crosscut model in the SWI. The experimental results enabled visualization of the SWI movement under the surface that could not be observed in the field. ^{Oostrom et al. (1992)} considered a laboratory experimental setup that is widely used for SWI. In contrast to the Henry problem, which involves a wide mixing zone, most laboratory experimental results indicate a relatively narrow mixing zone. ^{Goswami and Clement (2007)} provided the first benchmark problem based on experimental results by specifying the boundary conditions of the experimental model. The numerical simulations using 10 %, 50 %, and 90 % isochores for steady-state simulations were compared with the experimental results. In addition, the authors used only 50 % of the isochores in the transient simulations. An experimental dataset of horizontal and vertical distances with respect to time was provided for use as a benchmark dataset. ^{Abarca and Clement (2009)} succeeded in visualizing only the mixing zone, including the freshwater interface, in different colors by setting boundary conditions for freshwater and saltwater with different pH values in a small-scale water tank experiment.

Most benchmark problems, including the Henry problem, exist only as numerical simulations or semi-analytical solutions. ^{Chang and Clement (2012)} extracted and quantified the location of the salt wedge in a laboratory experiment and quantified the velocity of the movement of the freshwater–saltwater interface using an exponential function to describe a time-dependent equation. To quantify the difference between the intruding and receding SWI rate, the term ∆XT(t) was newly defined, which is the distance that remains to be traversed by a migrating saltwater wedge toe to reach the steady-state condition at the end of the simulations. This indicator showed that the wedge rapidly advanced during the early simulation period when the boundary condition was changed and gradually slowed down. ^{Chang and Clement (2013)} investigated the movement of contaminants in freshwater and saltwater regions based on the freshwater interface and reproduced the migration of slugs using numerical simulations. This study was extended to a site scale of decades to simulate the movement of a contaminant or nutrient plume. When comparing water circulation patterns between the freshwater and saltwater regions, the results showed that the movement of contaminants or nutrients in the saltwater area was relatively low. ^{Oz et al. (2014)} obtained experimental and numerical results from laboratory experiments. The results compared the salinity values of the 40 measuring points; these points were divided into three different ranges based on their salinity values. Sensitivity tests were performed using three dimensionless parameters: density, thickness, and flux.

Many studies have used the detailed vertical movement of the SWI in nearshore aquifers using two-dimensional numerical models. ^{Taniguchi et al. (2002)} highlighted the relationship between the extent of SWI inland and the increase or decrease in submarine groundwater discharge (SGD) and the importance of evaluating SGD as a coastal water resource management issue. SGD contributes greatly to the discharge and circulation of nutrients and pollutants near the coast and has attracted the interest of groundwater experts and experts in the marine engineering field. Changes in sea-level oscillations due to tides and waves can influence offshore and large-scale SGD patterns, which can affect the aquifer-ocean interface in terms of the upper saline plume (USP) and saltwater wedge ^{(Li et al., 2008; Li and Jiao, 2003a; 2003b; Michael et al., 2005; Robinson et al., 2007a; 2007b)}. The shape of sloping beaches has been altered by ocean oscillations, and numerical models have been used to study the effects of wave and tidal forcing. These models have confirmed that surface changes in unconfined aquifers can affect beaches ^{(Horn, 2006; Robinson et al., 2006)}. ^{Robinson et al. (2006)} extracted the vertical hydraulic gradient from numerical simulations to indicate upward flow according to tidal cycles. This study further investigated the effects of wave forcing on beach profile evolution using large-scale laboratory experiments with numerical modeling ^{(Bakhtyar et al., 2011)} and groundwater dynamics. The tidal influence simulated by ^{Kuan et al. (2012)} is indicated by the toe and upper edge of the salt wedge, width of the freshwater discharge zone, area of the USP, tide-induced circulation rate, and root mean square error (RMSE). Here, the RMSE was calculated based on the difference between the x coordinates of the interface given by the Glover solution and the simulation results for a given z value. The aquifer response to waves simulated by Robinson et al. ^{(Robinson et al., 2014)} is indicated by 1) the total salt mass per unit width of the aquifer in the USP, 2) the x-coordinate of the centroid of the salt plume, and 3) the z-coordinate of the salt plume centroid in the USP of a subterranean estuary. ^{Boufadel (2000)}, ^{Boufadel et al. (2006)} mesured concentrations of salt water and pressure heads from laboratory-scale experiments on a tidal-induced fresh brine interface.

One of the most intensively studied coastal aquifers in North America is the Biscayne Aquifer in southeastern Florida, USA ^{(Kohout, 1960)}. For example, ^{Dausman and Langevin (2005)} evaluated the relative importance of drought, well-field pumping, SLR, and canal management in SWI. ^{Ranjan et al. (2006)} focused on changes in groundwater recharge and SLR due to the loss of fresh groundwater resources in water- resource-stressed coastal aquifers. ^{Masterson (2004)} and ^{Masterson and Garabedian (2007)} simulated the underground water system of Cape Cod Island, Massachusetts, USA, which is well-known for tourism. Due to the significant abstraction of groundwater from the aquifer, the seawater-freshwater interface continues to move inland. ^{Gingerich and Voss (2005)} used the three-dimensional model SUTRA ^{(Provost and Voss, 2019)} to simulate SWI in Hawaii, USA. Oki ^{(Oki, 2005)} simulated SWI in Hawaii. ^{Langevin (2001)} investigated coastal groundwater problems in Biscayne Bay, Florida, USA. ^{Werner and Simmons (2009)} studied the SWI problem in Pioneer Valley, Australia. Table 2 provides an overview of representative numerically modeled SWI studies.

2.3 Application of Index in Changing Future Environment

Climate change, which induces environmental changes in coastal aquifers, is a critical issue that significantly increases uncertainty in water resource management due to groundwater salinization ^{(Green et al., 2011)}. In recent years, increased anthropogenic activities, extreme droughts, and SLR due to climate change have emerged as issues for coastal groundwater management. Climate change can affect long-term groundwater resources by affecting precipitation patterns and recharge pathways. Long-term recharge changes can affect groundwater resources by inducing drought or excessive water supply during rapidly changing or intensified seasonality. Numerous modeling-driven coastal groundwater studies have investigated the effects of climate change on future water resources. Future anthropogenic and natural changes were incorporated into these scenarios. ^{Nicholls et al. (2008)} defined a scenario for climate change and coastal vulnerability as “a description of potential future conditions developed to inform decision-making under uncertainty” ^{(Parson et al., 2007)}.

The SLR boundaries of SWI models were hypothetically simulated or exaggerated in early studies on the effects of climate change on SWI. The results were primarily analyzed using the salt wedge in the aquifer's crosscut or the salt concentration, measured in terms of salinity or total dissolved solids (TDS). ^{Meisler et al. (1984)} numerically analyzed the effect of eustatic sea-level fluctuations on the transition zone between freshwater and seawater in the Northern Atlantic coastal plain by lowering the relative sea level from 0 to 150 m. ^{Fesker (2007)} simulated an SLR of 0.5 m per 100 years over a simulation time of 250 years by comparing the salt wedge for these three cases. The results showed that a constantly rising sea level resulted in a linear increase in the salt load in the model domain with successive salinization of the marsh area from the coastline towards Geest, which could represent the effect of SLR at the freshwater-saltwater interface. ^{Oude Essink (2001)} established three types of simplified SLR scenarios to perform numerical analysis in aquifers in the Netherlands and concluded that a 0.5 m rise in sea level over 100 years would lead to an increase in salinity in all low-lying areas near the sea, providing maps of the freshwater head and the chloride concentrations. By comparing the measured and simulated water levels and the toe positions of 250 mg/L isochores, ^{Dausman and Langevin (2005)} concluded that the coastal aquifer of Broward County, Florida, could cause chloride contamination if the sea level exceeded 48 cm over 100 years. By comparing the volumetric loss of freshwater due to SLR in Israel, ^{Melloul and Collin (2006)} proved that 77 % of the loss was due to lateral movement and 23 % was due to head change, assuming that the SLR is 0.5 m in the coastal aquifers. ^{Navoy (1991)} graphically compared the positions of sharp interfaces from a crosscut model to study aquifer-estuary interactions based on the impact of several hypothetical SLRs and various combinations of hydraulic conductivities in the coastal aquifer of New Jersey, USA.

Aquifers are dynamically linked to the watershed hydrology and SWI studies have been conducted to integrate detailed hydrological analysis results ^{(Chang et al., 2016;} El-Kadi et al., 2016). This makes it even more important to generate detailed projection of groundwater recharge based on climatic future stresses. The climatic effect on groundwater resources is considered dependent on local climatic conditions and topography ^{(Green et al., 2011; Ranjan et al., 2006)}. ^{Ranjan et al. (2006)} calculated freshwater loss by considering land-use changes and the hydrological behavior of a catchment in Sri Lanka. Ranjan et al. studied the impact of climate change on high- (A2) and low-emission scenarios (B2) and observed that it may result in a decrease in semi-arid rainfall at the local scale, whereas increased temperatures are projected to increase precipitation at the global scale ^{(Ranjan et al., 2006)}. These two studies adopted an aridity index to study fresh groundwater loss. ^{Green and MacQuarrie (2014)} simulated the impact of climate change and groundwater abstraction in Atalatic, Canada to complete a vulnerability assessment and quantified the impact of SLR and changes after recharge on SWI in coastal aquifers.

Due to the effects of climate change, islands are significant research sites that are greatly impacted by changes in the water resources that are available. Several island regions have geological and terrestrial characteristics that make it difficult to use surface water resources, and most groundwater recharge in these regions relies on rainfall. Therefore, compared to coastal aquifers along the inland coast, groundwater resources on islands are more susceptible to SWI, and the islands can be sensitive to natural disasters depending on the island's hydrology, geology, topography, or shape ^{(van der Geest et al., 2020)}. Among the island regions, the most vulnerable to SWI are the atoll islands in the Pacific and Indian Oceans, which are in the lowlands less than a few meters above sea level and simultaneously have a relatively small area ^{(Bailey et al., 2009; Kundzewicz and DÖll, 2009; White and Falkland, 2010)}. In contrast, the impact of groundwater recharge due to rainfall changes will have a greater effect on the freshwater lenses of islands if shoreline erosion by SLR is insignificant ^{(White and Falkland, 2010)}. Studies have focused on SWIs impacted by climate-related stressors. The impact of changes in precipitation and SLR was reported in the Special Report on Emission Scenario (SRES) from the IPCC 4th report ^{(IPCC, 2007)}, which showed the movement of the saltwater wedge and change in volumetric freshwater lenses under both the most and least favorable scenarios applied to Shelter Island, USA ^{(Rozell and Wong, 2010)}. Mapping the water level distribution, impact of SLR, and impact of increased precipitation were projected until 2100 based on the IPCC SRES A2 scenario on freshwater lenses of the North Sea Island of Borkum, Germany ^{(Sulzbacher et al., 2012)}. The calibration method was illustrated by comparing the measured TDS at several monitoring wells to the computed results and the vertical profiles of electrical conductivities. The impact of a 25 % increase in the recharge rate was reported in Malaysia ^{(Praveena and Aris, 2010)}, and the impact of future climate and land-use changes during 2070-2099 under the SRES A2 scenario on the entire Jeju Island of South Korea was studied ^{(El-Kadi et al., 2014)}. ^{Babu et al. (2018)} simulated a reduction in freshwater volume on Tongatapu Island under recharge reduction and increased pumping scenarios. ^{Chang et al. (2016)} examined the effects of anthropogenic activities and recharge scenarios caused by climate change based on the IPCC SRES on Dauphin Island, USA, and indicated a volumetric reduction in freshwater by indicating changes in the freshwater head, salt wedge in the vertical crosscut of the model, and volumetric changes in freshwater resources.

^{Loáiciga et al. (2011)} simulated SLR scenarios for the city of Monterey, California, USA. The study illustrated numerically simulated locations of the vertically averaged 1000 and 10,000 mg/L isosalinity lines in the plan view of the model, where predictions of mean SLR are expected to range from 0.10 to 0.90 m. The simulation analysis included a detailed resolution of SWI, and its causes were quantified via numerical simulation under scenarios of changes in both groundwater extraction and SLR in the 21st century. One of the most recent studies was conducted by ^{Chun et al. (2018)} who simulated SLR and precipitation scenarios under the Representative Concentration Pathway (RCP) 4.5 and 8.5 scenarios based on the 5th IPCC report. The results showed the highest salinity increase of 40 % in the case of SLR and a less favorable freshwater recharge rate under RCP 4.5.

Recently, studies have considered more sophisticated climate change scenarios that consider recharge and SLR based on projected changes in precipitation and air temperature from global climate models and emission scenarios. In addition, interdisciplinary research has been conducted to map the SWI vulnerability in groundwater aquifers that are impacted by climate change ^{(Klassen and Allen, 2017; Werner et al., 2013)} and to inform the status of groundwater resources using indexed equations, various hydrogeological data, and GIS techniques. Various indicators have been developed in field research to assess the severity of SWI. ^{Giambastiani et al. (2007)} conducted a numerical study to investigate SWI in an unconfined coastal aquifer in Ravenna, Italy. By indicating the salt concentration distribution in the aquifer, the results showed that the mixing zone between fresh and saline groundwater shifted 800 m farther inland for an SLR of 0.475 m per century. In addition, the impact of the SLR was illustrated by indicating the seepage rate and salt load. Representative SWI studies related to changing environments are summarized in Table 3.

2.4 Recent Development of New Technology and Perspectives for SWI Indicators

Approaches to SWI research include numerical modeling, analytical solutions, and GIS. Recent modeling includes an uncertainty analysis and assesses the effect of alternative SWI reduction measures. Using new observational data and statistical techniques, uncertainty analysis was performed. ^{Coulon et al. (2021)} considered the uncertainties of an observation dataset and quantitatively reported the standard deviation and 95 % confidence intervals obtained from shallow, deep open, and pumping wells on Magdalen Island, Canada. To obtain an optimized extraction rate for pumping, ^{Rajabi and Ketabchi (2017)} tested five new formulations for uncertainty-based simulation-optimization approaches, which have the objective function of minimizing the sum of energy distances for various reliability levels in the management plan. ^{Roy and Datta (2020)} attempted ensemble-based SWI prediction by combining artificial intelligence-based prediction models and selecting five prediction models: adaptive neuro-fuzzy inference system (ANFIS), Gaussian process regression (GPR), multivariate adaptive regression spline (MARS), support vector regression (SVR), and probabilistic linear regression. The prediction performance was compared to that of the dataset obtained from the five monitoring locations using the following statistical metrics: rRMSE, MAPRE, R, NS, IOA, and KGE. To reduce uncertainty, new observation data or a combination of new and existing data were used in the model. For example, uncertainties can be reduced by applying observation data using new electrical signals such as airborne electromagnetic (AEM) ^{(Meyer et al., 2019)}. The AEM resistivity data in the cross-sectional model was mapped and visually compared to the simulated distribution of TDS. This demonstrated a further application for quantitative comparison with conventional indicators. ^{Goebel et al. (2019)} combined offshore AEM and onshore electrical resistivity tomography (ERT) and discovered that this approach was successful in regions previously considered inaccessible by traditional monitoring methods. ^{Hermans and Paepen (2020)} combined inland and ocean ERT to improve the reliability of observations.

Studies related to SGD and oceanic oscillations have considered the effects of tidal ^{(Veerapaga et al., 2019)} and seasonal changes in confined systems ^{(Qu et al., 2020)}. ^{Qu et al. (2020)} noted that most studies have investigated the effects of oceanic oscillations on unconfined aquifer systems and applied numerical approaches to confined systems using several indicators. To compare the groundwater discharge from the system, seasonally varying fresh submarine groundwater discharge (QFSGD), recirculated saline groundwater discharge (QRSGD), and fluxes flowing out of the aquifer across the seaward boundary (QSGD) were averaged over a seasonal cycle. Comparing freshwater storage capacities showed that an idealized circular system could store much more freshwater than an equivalent strip island. ^{Tang et al. (2020)} estimated the temporally varying cumulative volume of the channel (CVC) under different tidal levels from the coastal aquifer caused by surface and groundwater interactions through the river, and the change in the seawater-occupied channel volume, as a prediction and adaptation plan for SWI.

Compared to experimental results, Memari et al. investigated a circular, two-dimensional axisymmetric model for experiments on circular island systems ^{(Memari et al., 2020)}. ^{Veerapaga et al. (2019)} developed new indicators for the validation of simulation results using (a) water level, (b) salinity interface gradient (SIG), which is defined as the slope of a 10 % isoline, (c) salinity intrusion length for a salinity of 10 (SIL10), and (d) salinity intrusion length for a salinity of 5 (SIL5). Table 4 summarizes representative SWI studies that have recently developed SWI indicators.