Sharp Interface Approach for Regional and Well Scale Modeling of Small Island Freshwater Lens: Tongatapu Island

Abstract

Sustainable management of small island freshwater resources requires an understanding of the extent of freshwater lens and local effects of pumping. In this study, a methodology based on a sharp interface approach is developed for regional and well scale modeling of island freshwater lens. A quasi-three-dimensional finite element model is calibrated with freshwater thickness where the interface is matched to the lower limit of the freshwater lens. Tongatapu Island serves as a case study where saltwater intrusion and well salinization for the current state and six long-term stress scenarios of reduced recharge and increased groundwater pumping are predicted. Though no wells are salinized currently, more than 50% of public wells are salinized for 40% decreased recharge or increased groundwater pumping at 8% of average annual recharge. Risk of salinization for each well depends on the distance from the center of the well field and distance from the lagoon. Saltwater intrusions could occur at less than 50% of the previous estimates of sustainable groundwater pumping where local pumping was not considered. This study demonstrates the application of a sharp interface groundwater model for real-world small islands when dispersion models are challenging to be implemented due to insufficient data or computational resources.

1. Introduction

Groundwater occurring in the form of a thin lens floating on denser seawater is the primary source of freshwater in most small islands [1]. This limited groundwater resource is highly vulnerable to saltwater intrusion due to natural causes such as droughts, storm surges, sea level rise, and anthropogenic activities such as increased groundwater withdrawals [2,3]. Many small islands in the Pacific face groundwater shortages during droughts associated with El Niño—Southern Oscillation (ENSO) events [2] which are predicted to be frequent in the future due to global warming [4]. Public water supply in small raised coral islands of the Pacific including Tongatapu and Niue [5] relies on multiple vertical wells. Large-scale localized and unplanned abstraction can expedite the salinization of the freshwater lens [2]. Over-pumping can severely contaminate the aquifer, and hence, it is essential to modify and manage the pumping rates, especially during droughts. Management of groundwater withdrawal requires knowledge about the extent of freshwater lens and identification of wells that are under high risk of saltwater intrusion [6]. Hence, regional and wellscale modeling of freshwater lens are essential for planning water management strategies in small islands.

Various numerical models are available to assess the freshwater resources in small islands [7]. However, quantitative assessment of freshwater resources and seawater intrusion impacts have been performed in only a small number of real-world islands [8,9]. Numerical models, which consider the dispersion zone between the freshwater lens and underlying seawater, can closely represent typical island groundwater processes [10,11,12]. With fine discretization of space and time and sufficient monitoring data, dispersion models can provide a quantitative assessment of freshwater lens including the local effects of pumping [7,11,13]. Proper monitoring of well salinities is often scarce, if not absent, in many real islands [2] and it is difficult to obtain a sufficient quantity and quality of monitoring data for calibration [14,15]. Hence, due to the high computational demands and lack of sufficient monitoring data for calibration of heads and salinity, field-scale modeling of regional and well scale freshwater lens have been rarely attempted using dispersion models [7,11], and use of empirical models are suggested for well scale modeling [16]. As a result, use of a single numerical model to evaluate both the island-wide extent of freshwater lens and the local effects of pumping remains a challenge when the size of an island is not so small.

Numerical models that assume freshwater and saltwater to be immiscible and separated by a sharp interface have been used for the evaluation of risks of saltwater intrusion to freshwater resources when dispersion models are difficult to be implemented in regional scales [17]. Sharp interface numerical models have lesser data and computational requirements and have been used to predict salinity at pumping wells [18,19]. Hence, a sharp interface approach is used in this study for regional and well scale modeling of the freshwater lens on an island. It ensures a single numerical model is capable of providing a first-hand prediction of the extent of freshwater lens and the local effects of pumping.

Sharp interface two-phase numerical modeling was applied for Nantucket Island [20] to assess the impact of projected pumping on the freshwater lens. The model predicted the general possibility of saltwater upconing into the well field when the interface was within the intake of deepest well, though no attempt was made to quantify individual well salinization. Freshwater heads were used for calibration of the model, and no island-wide comparison of the model predicted freshwater thickness with the field observations were done. The midpoint of freshwater-saltwater mixing zone was matched to the sharp interface and hence the encroachment of water with salinity higher than the freshwater limit (chloride content of 250 mg/L) into well screen could not be identified.

When a numerical model for groundwater flow is developed, the model is generally calibrated with observed groundwater heads [21]. Some issues associated with the monitoring data of groundwater heads in small islands make it practically difficult to be used for calibration. The highest freshwater heads in small islands are typically of the order of 0.2 to 0.5 m above the mean sea level [2,22] due to the high hydraulic conductivity of the aquifer. Since the depth of the interface estimated by Ghyben-Herzberg’s approximation [23,24] is 40 times freshwater heads, highly accurate groundwater heads in the order of centimeter accuracy may be required to predict freshwater thickness. Additionally, due to the insufficiency of salinity monitoring data [2,7], the vertical profile of groundwater density is ambiguous, and so conversion of observed heads to equivalent freshwater heads is also challenging. In some islands like Tongatapu, which is the case study selected here, accurate land surface elevations are not available at monitoring locations.

Freshwater thickness is the most critical parameter to be assessed for groundwater management in small islands [7]. In this study, the freshwater thickness is selected for calibration of the groundwater model to overcome the challenges associated with the use of freshwater heads in calibration. The two-fluid sharp interface approach used in this study [25] represents groundwater flow in terms of freshwater and saltwater heads. Though it is unable to predict the nature of the transition zone, it can predict the response of the interface [26] and pumping wells [18,27] to applied stresses. By adopting the upper limit of freshwater electrical conductivity as the interface, it is possible to predict the extent of freshwater thickness in the island.

The objective of this work is to develop a methodology based on a sharp interface numerical model to conduct regional and well scale modeling of island freshwater lens under long-term stresses. The methodology is applied to Tongatapu, a small Pacific Island. A numerical model is developed and calibrated for the average recharge and the current pumping condition and is then used to assess the effect of long-term stresses of recharge and pumping variations on saltwater intrusion and well salinization. This work expands the capacity of the sharp interface numerical model to regional and well scale modeling of freshwater lens in real-world small islands. The methodology developed can be applied to other similar island systems for the prediction of freshwater lens and to quantify the local effects of pumping. The impact of long-term stresses on saltwater intrusion along with well salinization is evaluated using the groundwater system of Tongatapu and generalized outcomes are relevant to other similar island systems with concentrated pumping. Tongatapu is the largest and the most populated island of the Kingdom of Tonga, a Pacific island nation. There has been an increasing trend of migration to Tongatapu, more specifically, to its capital Nuku’alofa from outer islands [28]. Groundwater extracted from vertical wells meets the public water supply demands for Nuku’alofa due to the lack of surface water sources [29]. With the increase in groundwater demands and expected climate changes in the future, there has been an increase in public concern over groundwater availability. Previous numerical modeling studies on Tongatapu [30,31] used a large grid size (about 1 km) for simulations and ignored the effects of groundwater pumping and future change in recharge patterns. For small islands like Tongatapu where groundwater pumping is concentrated over a small area, the risk of well salinization is high due to saltwater upconing [2]. The results of this study can help to identify the wells with a high risk of salinization, which was rarely attempted in earlier modeling studies and can contribute to planning and management of groundwater development.

2. Materials and Methods

2.1. Study Area

Tongatapu (21°03′–21°16′ S and 175°02′–175°21′ W) is the main island of Kingdom of Tonga (Figure 1). It has a total area of 259 km² and is the location of the national capital Nuku’alofa. The population of Tongatapu was 75,416 (73% of total population of the Kingdom of Tonga), with a density of 290 persons per square kilometers in 2011 [32]. The island is 31.5 km along its Northwest-Southeast directions and 18.5 km along the Northeast-Southwest direction.

Tongatapu is a tilted raised limestone island with a maximum elevation of 65 m [22] above mean sea level. The aquifer consists of karst-type limestone that is generally very porous and has hydraulic conductivity ranging from 500 to 3600 m/day [22,33]. The drill logs indicated a single layer of limestone aquifer [34]. The freshwater stored in the limestone aquifer is generally in the form of a lens floating on denser salt water. The water table of the island of Tongatapu is almost flat with an average hydraulic head of 0.4 m [22], and the groundwater lens with an electrical conductivity less than 2500 µS/cm was measured to have a maximum thickness of about 14 m [22,33,35]. The freshwater thickness at a few locations like Liahona and Fua’amotu (Figure 1) was estimated from monitoring wells. However, it is less understood across the remainder of Tongatapu as only limited observation data is available [34]. Fanga’uta lagoon (Figure 1) opens to the sea and salinity is reported to be approximately 25,700 to 32,900 ppm dissolved solids [36].

The semi-tropical climate of Tonga is highly vulnerable to the effects of El Nino, which frequently coincides with drought incidence [37]. Tongatapu has a mean temperature of about 24 °C [38]. Mean annual rainfall for Nuku’alofa, Tongatapu from the year 1993 to 2016 was 1781 mm. The rainfall pattern is characterized by the contrast between a wet season (November–April) which accounts for over 60% of total annual rainfall and a dry season (May–October).Values ranging from 20–35% of the mean annual rainfall have been reported as annual recharge rates in Tongatapu [22,33,36].

Groundwater extraction in Tongatapu is divided into public/urban and village/rural water supply. The public well field at Mataki’eua-Tongamai (Figure 1) is the main source of water supply to the capital city Nuku’alofa. The well field consists of 39 vertical wells concentrated over an area of 0.57 km² [29]. As per the Ministry of Lands and Natural Resources (MLNR) Tonga, the groundwater extraction from the public well field in 2016 was 10410 m³/day. Rural water supply varies between villages and most of the wells are located at the edges of villages. It is estimated that there are about 200 village wells in Tongatapu [36]. It is observed that the current pumping from village wells is difficult to be estimated, as there are no flow meters installed [39].

There are sixteen salinity monitoring boreholes (SMB) throughout the island, ten around the Mataki’eua-Tongamai well field, and three each around Hihifo and Fua’amotu (Figure 1). Each salinity monitoring borehole (SMB) in Tongatapu consists of a nested group of four to six measuring tubes with 1 m long screen near the bottom of each tube (Figure 2). A screen of each tube either contains fresh or non-fresh water. The freshwater limit adopted in Tongatapu is based on the European Commission (1998) maximum electrical conductivity (EC) of 2500 µS/cm [22]. Figure 2 illustrates a situation where the screens of tubes 3 and 4 contain freshwater, while those of 1 and 2 contain non-freshwater. In this case, the lower boundary of the freshwater lens lies between screens 2 and 3. The freshwater thickness at the monitoring location is estimated as depth to EC value of 2500 µS/cm from the water table. The Tonga government (MLNR, Tonga 2017) adopted linear interpolation between EC values measured below and above 2500 µS/cm and their corresponding depths to obtain the freshwater thickness. This linear interpolation method results in a slightly conservative estimate of freshwater thickness compared to that based on relative salinity used for the estimation of freshwater thickness based on monitoring data in Majuro Atoll [40].

2.2. Sharp Interface Freshwater-Seawater Groundwater Flow Numerical Model

2.2.1. Governing Equations

In this study, a sharp-interface numerical model described by Huyakorn et al. [25] is chosen to simulate fresh and saline groundwater in the island. The sharp interface model is less rigorous than a dispersion model in treating physics relevant to the problem at hand. However, the former model requires fewer data and low-resolution grids so that it can be applied with relative ease to a large area. This model considers the flow dynamics of both freshwater and saltwater. The sharp interface freshwater-seawater model for groundwater flow with free surface consists of two vertically integrated governing equations, one for freshwater flow and the other for saltwater flow:

$$\frac{\partial }{\partial \mathrm{x}}\left(K_{\mathrm{x}}^{\mathrm{f}}{b}^{\mathrm{f}}\frac{\partial h^{\mathrm{f}}}{\partial \mathrm{x}}\right)+\frac{\partial }{\partial \mathrm{y}}\left(K_{\mathrm{y}}^{\mathrm{f}}{b}^{\mathrm{f}}\frac{\partial h^{\mathrm{f}}}{\partial \mathrm{y}}\right)=S_{\mathrm{y}}^{\mathrm{f}}\frac{\partial h^{\mathrm{f}}}{\partial \mathrm{t}}-\theta \frac{\partial \xi }{\partial \mathrm{t}}-Q^{\mathrm{f}}$$

(1)

$$\frac{\partial }{\partial \mathrm{x}}\left(K_{\mathrm{x}}^{\mathrm{s}}{b}^{\mathrm{s}}\frac{\partial h^{\mathrm{s}}}{\partial \mathrm{x}}\right)+\frac{\partial }{\partial \mathrm{y}}\left(K_{\mathrm{y}}^{\mathrm{s}}{b}^{\mathrm{s}}\frac{\partial h^{\mathrm{s}}}{\partial \mathrm{y}}\right)=S_{\mathrm{y}}^{\mathrm{s}}\frac{\partial h^{\mathrm{s}}}{\partial \mathrm{t}}+\theta \frac{\partial \xi }{\partial \mathrm{t}}-Q^{\mathrm{s}}$$

(2)

where the superscripts $\mathrm{f}$ and $\mathrm{s}$ refer to freshwater and saltwater, respectively, and $h$ is the hydraulic head, which is vertically averaged for the freshwater zone and saltwater zone. $b^{\mathrm{f}}$ and $b^{\mathrm{s}}$ are the thickness of freshwater and saltwater zones in the aquifer. $K_{\mathrm{x}}^{\mathrm{f}}$, $K_{\mathrm{y}}^{\mathrm{f}}$ and $K_{\mathrm{x}}^{\mathrm{s}}$, $K_{\mathrm{y}}^{\mathrm{s}}$ are the hydraulic conductivities in $\mathrm{x}$ and $\mathrm{y}$ directions for freshwater and saltwater, respectively. $\theta$ is effective porosity, $\xi$ is the height of saltwater- freshwater interface above the datum. $S_{\mathrm{y}}^{\mathrm{f}}$ and $S_{\mathrm{y}}^{\mathrm{s}}$ are specific yield of the aquifer in freshwater and saltwater zone. The coupled freshwater and saltwater flow equations are solved simultaneously to obtain freshwater heads and saltwater heads. The interface elevation is calculated by equating pressures at the interface [41];

$$\xi =\frac{{\rho }_{\mathrm{f}}}{\left({\rho }_{\mathrm{s}}-{\rho }_{\mathrm{f}}\right) }\left[\left({\rho }_{\mathrm{s}}/{\rho }_{\mathrm{f}}\right)h^{\mathrm{s}}-h^{\mathrm{f}}\right]$$

(3)

where ${\rho }_{\mathrm{f}}$ and ${\rho }_{\mathrm{s}}$ are the freshwater and saltwater densities.

When the freshwater-saltwater interface is within the well screen, it is assumed that both freshwater and saltwater are extracted simultaneously. The total pumping rate $Q^{\mathrm{T}}$ at the well is considered to be the sum of $Q^{\mathrm{f}}$ and $Q^{\mathrm{s}}$, the volumetric fluxes per unit volume of aquifer due to pumping of freshwater and saltwater, respectively. $Q^{\mathrm{f}}$ and $Q^{\mathrm{s}}$ are related to the total pumping rate $Q^{\mathrm{T}}$ as follows.

$$Q^{\mathrm{f}}=\left(\frac{K^{\mathrm{f}}{L}^{\mathrm{f}}}{K^{\mathrm{f}}{L}^{\mathrm{f}}+K^{\mathrm{s}}{L}^{\mathrm{s}}}\right)Q^{\mathrm{T}}$$

(4)

$$Q^{\mathrm{s}}=\left(\frac{K^{\mathrm{s}}{L}^{\mathrm{s}}}{K^{\mathrm{f}}{L}^{\mathrm{f}}+K^{\mathrm{s}}{L}^{\mathrm{s}}}\right)Q^{\mathrm{T}}$$

(5)

where $L$ is the total length of the screen, $L^{\mathrm{f}}$ and $L^{\mathrm{s}}$ are the length of the screen exposed to freshwater and saltwater, respectively. $K^{\mathrm{f}}{L}^{\mathrm{f}}$ and $K^{\mathrm{s}}{L}^{\mathrm{s}}$ represent transmissivities of freshwater and saltwater regions, respectively. The saltwater ratio (SWR) at a pumping well is obtained as the ratio of saltwater pumped to the total pumping at the well.

$$SWR=\left(\frac{Q^{\mathrm{s}}}{Q^{\mathrm{f}}+Q^{\mathrm{s}}}\right)$$

(6)

Despite the simplicity of the above approach, it was demonstrated that saltwater ratios at pumping wells could be modeled with reasonable accuracy [18,27].

A modified Galerkin finite element method was used to approximate governing equations [25]. The numerical model [25] has been improved in various aspects and applied to diverse groundwater problems such as assessment of saltwater intrusion due to groundwater developments [42], evaluations of potential groundwater development [43], coastal groundwater discharge for complex coastlines [44], estimation of saltwater contents in pumping wells [18,27], assessment of impacts of combined fresh and saltwater developments [45], and evaluation of the effectiveness of freshwater injection and/or saltwater pumping to prevent saltwater intrusions [46,47].

2.2.2. Conceptual Model for Groundwater Flow in Tongatapu Island

The numerical model was applied to the entire Tongatapu Island. Rectangular finite elements of 100 m × 100 m are used for areas with less groundwater development, and elements of 25 m × 25 m are applied for the central part of the island where public wells are concentrated. The grid was selected based on a grid convergence test. Consequently, the numerical model consists of 89,200 elements and 89,824 nodes. A single homogeneous isotropic unconfined aquifer with the base elevation at −40 m with reference to mean sea level (MSL) is considered. The effective porosity of the aquifer is assumed to be uniform as 0.3 [22]. Densities of freshwater and seawater are 1000 and 1024 kg/m³, respectively. A no-flow boundary is assumed at the bottom surface of the aquifer. At those nodes along the coastline and lagoon first type and efflux-only third type boundary conditions are specified for the saltwater equation and the freshwater equation, respectively. Recharge is applied with flux boundary conditions at all internal nodes.

Recharge was estimated from a water balance model provided by APEC Climate Center (APCC) based on the soil water balance model for a small island with no runoff [22]. The surface runoff was neglected as in other small islands due to flat terrain, high infiltration, and lack of surface flow observations. The water balance model parameters used in the present study (Table S1) are based on the sensitivity analysis by White et al. [22] for Tongatapu. As the daily rainfall data was unavailable, monthly rainfall data was used as input. The available monthly rainfall from 1993–2016 at Nuku’alofa was used to estimate the average annual recharge as 570 mm/year. This is 32% of the average annual rainfall during the period and falls within the range of recharge estimates by previous studies [22,33,36].

The locations of 39 public wells and 56 village wells are known [22], and hence a total of 95 pumping wells have been considered (Figure 1). Actual depths, lengths of well screens and individual pumping rates are not known. Screens are assumed to extend from −2 m to −5 m (MSL) for public wells and 0 m to −1 m (MSL) for village wells [34]. Groundwater extraction estimates for the public well field in 2016 was 10,410 m³/day (MLNR Tonga, 2017). The above pumping rate is assumed to be equally distributed in the 39 public wells as 267 m³/day/well. Tonga Water Board estimates rural per capita demand to be 80 L per person per day with rainwater availability [48]. Using the rural population of Tongatapu as per 2011 census and assuming 50% losses [22] the rural groundwater extraction is estimated to be 6300 m³/day in this study. This is equally distributed to the 56 village wells included in the model as 112.5 m³/day/well. Thus, the total groundwater extraction for Tongatapu is approximated to be 16,710 m³/day. This value is close to estimates of groundwater demands for Tongatapu by Waterloo and Ijzermans [39].

2.3. Calibration of the Model

Depths to the groundwater table and electrical conductivity (EC) values are measured at salinity monitoring boreholes (SMB) in Tongatapu. However, measured depths cannot be converted to groundwater heads as land surface elevations of the reference point at SMBs from where depth to water table is measured and is not known precisely. In this study, freshwater thickness values are used as the calibration targets. The aquifer parameter of horizontal hydraulic conductivity based on available data was selected and later adjusted to obtain simulated freshwater thickness, which satisfactorily matches observed freshwater thickness in 16 salinity monitoring boreholes (SMB) distributed across the island (Figure 1). The Tonga government (MLNR Tonga) provided an average freshwater thickness in each SMB based on quarterly salinity monitoring data from July 1997 to April 2017. Seven out of the sixteen salinity monitoring boreholes were constructed in 1997, and the rest were added in 2012. On average, 56 observations are available for the older SMBs and 20 for the new SMBs. The average coefficient of variations (standard deviation divided by mean) are 0.18 and 0.13 for old and new SMBs, respectively, indicating a relatively stable freshwater thickness at each SMB. The manual trial and error model calibration, though unable to lead to a unique solution, is constrained based on the data available [26]. The simulated groundwater flow field was based on the steady state assumption despite temporal variation in recharge. Groundwater pumping was known to begin in 1960’s and has been increasing since then. Nevertheless, it is assumed that the present groundwater flow field has reached the equilibrium state to the current pumping rate.

2.4. Long-Term Stresses

2.4.1. Impact of the Current Pumping

Simulations start with an initial condition of the aquifer fully saturated with saltwater below the MSL. Transient simulations are continued till the total saltwater and freshwater volume in the domain reached the steady state. Pumping wells are removed from the calibrated model to estimate the groundwater flow for the predevelopment condition prior to the groundwater development. All other conditions except for the pumping wells were assumed applicable to the predevelopment simulation. The predevelopment model was also run under the steady state assumption.

2.4.2. Impact of Long-Term Stresses

In this study, two major long-term stresses, droughts, and overpumping, which can cause saltwater intrusion and well salinization in Tongatapu are considered. Hydrological drought analysis of past rainfall records of Tongatapu observed that the average duration of droughts was 37 months and they occurred on nearly every 11.5 years with a wide range in frequency and duration [22]. The longest period of zero recharge was 18 consecutive months for past rainfall records. It was also observed that the salinity of water pumped from Mataki’eua/Tongamai well field is most strongly correlated with the amount of rainfall received in the preceding 12 to 18 months [22,49]. Groundwater development from Mataki’eua/Tongamai well field began in 1966 (MLNR, Tonga 2017). The current groundwater withdrawal in the same well field is over eight times the extraction at the beginning of pumping (MLNR, Tonga 2017). Many wells in the well field use diesel pumps that pump at about 270 m³/day. Possible use of electric pumps and the growing water demands can increase the pumping rates (MLNR, Tonga 2017) in the future.

Steady-state simulation results can indicate the potential impact and the worst case of the climate change scenario and have been used in many previous studies [13,50]. Current recharge value was progressively decreased to estimate the recharge rate at which well salinization started. While pumping rate was maintained at current value, well salinization occurred at about 69% of recharge. Similarly, pumping rate was increased from the current value to obtain pumping rate at which saltwater intrusion occurred in public wells. For current recharge values, 160% of present pumping initiated well salinization. Thus, in this study, six different scenarios (

Table 1. Long-term stresses scenarios.

Sl.No	Scenario	Recharge (m/year)	Pumping from Public Wells (m3/day)
1	Base (current)	0.570 (R)	10,410 (P)
2	0.69 R	0.393	10,410
3	0.65 R	0.371	10,410
4	0.60 R	0.342	10,410
5	1.6 P	0.570	16,656
6	1.8 P	0.570	18,738
7	2.0 P	0.570	20,820

) are considered to simulate the possible impacts of changes in recharge and pumping on the freshwater lens and well salinization. Scenarios are selected to cover the increasing impact of drought and groundwater development conditions on well field salinization. The base scenario represents the calibrated model at base recharge and pumping conditions. The first three scenarios represent potential droughts, and the next three represent scenarios of increased groundwater development. Either recharge or pumping is varied from the base case in each scenario so that the relative impact of recharge and pumping can be compared.

3. Results

3.1. Results of Calibration

The hydraulic conductivity value of 3600 m/d was obtained as a result of manual calibration. The value is near the upper end of the values reported in the literature and close to that from pumping tests [22]. Statistical results of freshwater thickness calibration are given in

Table 2. Statistical results of simulated freshwater thickness.

Parameters	Values
Number of observation points	16
Correlation coefficient	0.87
Range in observations	8.23 m
Maximum residual	2.65 m
Minimum residual	−2.84 m
Residual mean	0.27 m
Absolute residual mean	1.15 m
Root mean squared error	1.44 m
Normalized RMS error	13.8%
Scaled RMS error	17.5%

. The simulated and observed freshwater thickness at salinity monitoring boreholes indicate a correlation coefficient of 0.87 with a normalized root mean square error of 13.8%.

Figure 3a shows the comparison of simulated freshwater thickness and observed freshwater thickness at the 16 SMBs. The contours based on average observed freshwater thickness at SMBs along with simulated freshwater thickness across the island are shown in Figure 3b. Spatial distribution of residuals at SMBs are also indicated in Figure 3b. The model captures the spatial distribution of freshwater thickness though there are differences at salinity monitoring boreholes. It is reported that the agricultural fields around SMB08 rely heavily on groundwater from village wells [22]. As all the village wells in this region have not been included in the model due to the lack of well data (MLNR Tonga, 2017), freshwater thickness may have been overestimated at SMB08, indicating a minimum residual of 2.84 m Fua’amotu region where SMB09 and SMB10 are located and receive higher rainfall than the rest of Tongatapu [22]. The non-uniform spatial distribution of recharge could have resulted in an underestimation of freshwater thickness by the model in the above region. The maximum residual of 2.65 m is at SMB10. The assumption of uniform recharge and pumping could have led to differences in freshwater thickness at other locations.

3.2. Current State of Freshwater Lens and Impact of Current Pumping

The current state of freshwater lens thickness is indicated in Figure 3b. The total volume of recharge is 147.6 MCM/year, and the total pumping rate is 6.1 MCM/year, which amounts to roughly 4.1% of the recharge. The maximum thickness is 14.1 m occurring around the Fua’amotu region, and the average thickness of freshwater is 8.4 m for the island. The regions with a freshwater lens thicker than 10 m are Fua’amotu, Liahona, and Kolonga region. Freshwater volume is calculated as the product of porosity and sum of the area of each element multiplied by freshwater thickness in the element. The estimated volume of freshwater contained in regions where the lens thickness is greater than 10 m is 294 MCM considering a porosity value of 0.3. The maximum freshwater head under sustained steady-state pumping is obtained as 0.329 m. At the current rate of pumping used in the study, it is observed that there is no seawater intrusion into both the public and village wells. Detailed statistics of the lens are specified in

Table 3. Impacts of the current pumping condition.

Parameters		Current Pumping	Predevelopment	Percentage Decrease
Maximum freshwater head (m)		0.329	0.347	5.19
Minimum interface elevation(m)		−13.72	−14.46	5.12
Freshwater thickness (m)	Maximum	14.05	14.81	5.13
	Average	8.40	9.02	6.87
Area (km2) where lens thickness is	>5 m	206.5	208.2	0.79
	>7 m	167.3	170.8	2.08
	>10 m	83.1	91.7	9.39
Volume (MCM) of freshwater where lens thickness is	>0 m	636.7	658.4	3.30
	>5 m	581.4	604.7	3.86
	>7 m	510.1	536.9	5.00
	>10 m	294.8	334.6	11.90

.

The impact of the current groundwater development on the predevelopment freshwater lens is investigated. Statistics of the freshwater lens for predevelopment and post-development conditions are presented in Table 3. Average freshwater head, maximum freshwater thickness, and minimum interface elevations decrease by about 5%. However, area and volume of the freshwater lens with thickness greater than 10 m decrease by 9% and 12%, respectively, for the post-development simulation. The impact of concentrated pumping from the public well field at Mataki’eua/Tongamai is evident in Figure 4 where the contours indicate a decrease in predevelopment freshwater thickness due to current pumping. Though the well field covers less than one km², the freshwater lens thickness reduced by over 0.5 m, 1 m, and 2 m in an area of 43 km², 19.3 km², and 12.3 km², respectively. Contours are asymmetric in shape and are concentrated along the coast as drawdown is restricted at the coast where specified head boundary conditions are defined. As the village wells are distributed across the island and have lower pumping rates effect of the village well pumping on the freshwater lens is not significant.

3.3. Predicted Effects of Long-Term Stresses

3.3.1. Impact on Freshwater Lens

The effect of long-term stresses is indicated as a change in the present state of the freshwater volume. The percentage of change in total freshwater volume (>0 m thick), the volume of freshwater where the thickness is more than 5 m, 7 m, and 10 m in comparison to the base case are indicated in Figure 5. Detailed freshwater statistics under scenarios of reduced recharge are presented in Table S2. The total volume of freshwater lens is decreased by 18% for scenario 0.69 R and by 24% for scenario 0.6 R. The deeper freshwater lens indicated by freshwater volume with lens thicker than 10 m decreased by 65% and 86% for scenarios 0.69 R and 0.6 R, respectively, indicating a higher sensitivity of the thicker part of the freshwater lens to applied stresses. The effect of increased pumping on freshwater lens volume change is much less compared to that of reduced recharge as pumping is concentrated over a smaller area. Detailed freshwater statistics under scenarios of increased pumping are presented in Table S3. The total freshwater volume decreased by 2% and 3% for scenario 1.6 P and 2 P, respectively. The part of freshwater lens thicker than 10 m resulted in reductions of 9% and 13% for scenarios 1.6 P and 2 P, respectively.

3.3.2. Well Salinization

Saltwater intrusion into wells is quantified in terms of the number of wells salinized and saltwater ratios in each well. With the best estimates of pumping rates and screen lengths, no well salinization was observed at the village well locations in the model, under the scenarios considered in this study.. Spatial distribution of saltwater intrusion into public wells under various scenarios of reduced recharge rates is shown in Figure 6. A saltwater ratio of 0.001 is considered to be the onset of well salinization. For scenario 0.69 R, two wells located at the center of the well field get salinized. With further reduction in recharge wells around the center and those near the lagoon are intruded by saltwater. Though freshwater thickness is reduced across the well field due to a decrease in recharge, concentrated pumping expedites saltwater intrusion into well field from the lagoon side, indicating lateral saltwater intrusions. Spatial distribution of saltwater intrusion into public wells under increased pumping is shown in Figure 7. As expected, wells located at the center of the well field get salinized first. Increased pumping rates cause saltwater upconing and more wells around the center are salinized.

The percentage of total wells intruded by saltwater, and maximum saltwater ratios for different scenarios are given in

Table 4. Well salinization under different scenarios.

Sl.No	Scenario	Percentage of Wells Intruded by Saltwater	Maximum Saltwater Ratio
1	Base (current)	0.0	0.00
2	0.69 R	5.1	0.004
3	0.65 R	38.5	0.13
4	0.60 R	51.3	0.22
5	1.6 P	12.8	0.05
6	1.8 P	38.5	0.24
7	2.0 P	53.8	0.35

. For scenarios with reduced recharge 0.69 R, 0.65 R, and 0.6 R, saltwater intrusion occurs in two, 15, and 20 of the public wells, respectively. The maximum saltwater ratio was 0.21 for 0.6 R. While increased pumping causes saltwater intrusion in five, 15, and 21 public wells for 1.6 P, 1.8 P, and 2 P, respectively. Maximum saltwater ratios for the worst case pumping scenario (2 P) is 0.35, whereas it is 0.22 for the worst case recharge (0.6 R). Both the above scenarios cause around 50% of wells to be intruded by saltwater.

4. Discussion

The numerical model developed in this study revealed the quantitative behavior of freshwater lens and replicated observed freshwater thickness with a root mean square error of 1.44 m. Island aquifers are characterized by high heterogeneity and lack of data [7]. Previous modeling studies on real islands using dispersion models reported scaled-RMSE of 10% for Bonriki Island model with an area of 9.5 km² [12] and 21% for Kish Island model with an area of 112 km² [10] for modeled versus measured salinity. Scaled RMSE for modeled versus measured freshwater thickness is 17.5% in this study (Table 2). Given the large size of the island (259 km²) and simplifying assumptions, the model-measured difference that is not unacceptable.

The total volume of freshwater lens is estimated to be 637 MCM for an annual recharge of 147.6 MCM. Current groundwater development is estimated at 6.1 MCM/year. The volume of the freshwater lens with the thickness greater than 10 m is less than 300 MCM, and exist in less than one-third of the total area of the island (Table 3). It is seen that the current pumping rate of about 4.1% of annual recharge can be sustained by a freshwater lens with no saltwater intrusion under the assumption of uniform recharge throughout the year. However, its impact on predevelopment groundwater lens appears to be significant. The area and volume of the freshwater lens with thickness greater than 10 m are reduced by 9.4% and 11.9%, respectively, as compared with those for the predevelopment lens. Current pumping reduced freshwater thickness in an area of 12.3 km² around well field by over twometers, as compared to the predevelopment state of freshwater lens. Jakovovic et al. [51] adopted an interface rise of two meters to delineate a region around a pumping well as saltwater upconing zone of influence (SUZI). A limit of two meters is about 25% of the average freshwater thickness of Tongatapu. Considering the above criteria of the saltwater upconing zone of influence, an area over 20 times well field size indicates a very significant reduction in freshwater thickness in Tongatapu.

Scenarios considered in the study investigated the effect of drought and increased pumping. For a scenario with a 40% decrease in average annual recharge, 86% of freshwater lens volume with thickness more than 10 m is lost. For pumping at 8.2% of the average annual recharge, loss of freshwater lens volume with thickness more than 10 m is 13%. Groundwater development in Tongatapu is concentrated over a smaller area which leads to gradient reversals (ocean to land) resulting in saltwater intrusions (Figure 7). Salinization risk of each well is influenced by location from the center of the pumping well field and also the distance from the lagoon (Figure 6 and Figure 7). For scenarios with a reduction in recharge (0.69 R, 0.65 R, and 0.6 R), lateral saltwater intrusions are prominent, whereas saltwater intrusion due to upconing under wells is significant for scenarios with increased pumping rates (1.6 P, 1.8 P, and 2 P). However, wells with the highest risk of salinization are almost common for scenarios of reduced recharge and increased pumping. Identification of high-risk wells helps in the planning and management of groundwater pumping.

Among the scenarios considered in this study, the worst cases of 0.6 R and 2 P cause over 50% of public wells to be intruded by saltwater (Table 4), even though total freshwater volume change is 24% for 0.6 R and 3% for 2 P, as compared to the base case (Figure 5). Hence, evaluation of freshwater volume change alone may not indicate the potential threats due to pumping. In scenario 2 P, though the pumping rate is about 8.2% of the average annual recharge, 50% of wells are intruded by saltwater. Well salinization of over 50% of number of public wells indicates a reduction in freshwater production by 50%. This further emphasizes the need for groundwater management and planning.

To evaluate the relative impact of pumping and recharge, the total volume of water extracted is quantified as the reduction in recharge [22]. This approach can also evaluate the relative impact of pumping distributed across total area and localized pumping. As base pumping rates are about 4% of recharge, scenario 1.8 P represents a 7% decrease in recharge. Well salinization, in terms of the number of wells intruded, under this scenario, is comparable to that of 0.65 R where effective recharge reduction is 39%, including pumping. Similarly, scenario 2 P and 0.6 R are quite comparable in terms of the number of wells salinized. It is seen that pumping impacts on well salinization are about five times that of recharge for Tongatapu Island. The high relative contribution of pumping to well salinization can be attributed to the high concentration of wells in the small area in Tongatapu.

Since seawater intrusion into pumping wells is a local phenomenon, sustainable groundwater development from a limited number of pumping wells is much smaller than the estimated volume of freshwater. White et al. [22] estimated island-wide sustainable groundwater development to be 20% of groundwater recharge and total sustainable groundwater yield considering areal pumping rates as 54,000 to 72,000 m³/d. The numerical model, considering the current distribution of pumping wells, indicates that extractions at about 6.5% of the annual average recharge, and 50% of the lower range of areal pumping rates, as mentioned above, can cause saltwater intrusions.

Even though the precipitation exhibits seasonal and yearly patterns, the steady-state assumption was not too unrealistic as the groundwater system responded quickly due to large hydraulic conductivity. The salinity monitoring wells in Tongatapu use multi-nested tubes resulting in different values of heads in each tube due to the density differences in the tubes. The mean water table elevation for Tongatapu from observations in village wells for 1971–2007 was 0.41 m above mean sea level [22]. Error in deducing head from multi-nested monitoring wells might be large, as compared to the head above the MSL at the location. Hence, careful attention is required to make full use of the water level observations from multi-nested tubes.

Spatial heterogeneity of the aquifer and recharge distribution have not been considered. The tidal effects were not included in the study. The model developed considers the sharp interface at the freshwater limit. Hence, the fluid above the interface is freshwater, but in reality, the fluid below is not completely saltwater. Error of neglecting mixed saltwater below the interface on general flow pattern needs to be investigated further. While all the public wells for urban water supply have been included in the model, only 56 village wells were represented out of the over 200 village wells known to exist [36] as there was no information about the location of other wells. Estimated total groundwater development rate of the island was applied uniformly to all wells in the well field as individual pumping rates were not known. A more accurate representation of wells, in terms of locations, pumping rates, and screen length could result in better estimates of well salinization.

5. Conclusions

Prediction of the extent of freshwater lens and quantification of the local effects of pumping is essential for the design and development of effective water management strategies in small islands. In this study, a steady-state sharp-interface numerical model of fresh and saline groundwater flow is developed to evaluate saltwater intrusion and well salinization in Tongatapu Island. The model has a comparatively fine grid size of 25 m × 25 m in public well field and 100 m × 100 m in other areas, which is able to better represent the freshwater distribution and saltwater intrusion characteristics at pumping wells. Freshwater thickness is used as the calibration target in this study with the sharp interface corresponding to the upper limit of freshwater electrical conductivity adopted in the Kingdom of Tonga. Such an approach has not been reported in previous applications of sharp interface numerical models for islands. Thus, the current study is the first, to the best of our knowledge, to conduct a comprehensive analysis of freshwater lens and well salinization using sharp interface approach in real islands.

Model results indicate that at present conditions, 46% of the total freshwater volume of Tongatapu, estimated around 630 MCM, is from freshwater lens thicker than 10 m and concentrated in two regions towards the center of the island. This may indicate sufficient availability of resources for groundwater development. However, the concentration of pumping wells in a single well field increases the risk of saltwater intrusion. Various scenarios of decreased recharge and increased pumping were considered to evaluate saltwater intrusion and well salinization. More than 50% of public wells are intruded by saltwater for 40% decreased recharge (0.6 R) or groundwater pumping at 8% of average annual recharge (2 P). It is seen that pumping impacts on well salinization are about five times that of recharge for Tongatapu Island. The higher relative contribution of pumping to well salinization can be attributed to the high concentration of wells in the small area. The results suggest that the risk of salinization for individual well depends on the distance from the center of well field and the distance from lagoon or coast. However, the closeness of the well to the center of well field results in higher saltwater ratios and the extent of upconing. Major factor leading to well salinization was saltwater upconing rather than lateral saltwater intrusions. Hence, it is imperative to include pumping wells in a groundwater model, especially when pumping is localized.

While the model is currently limited to steady state analyses, it can help water managers of the island to assess potential impacts on groundwater resources caused by changes in some important hydrogeologic parameters, such as pumping rates and recharge rates. It also helps to identify wells with a high risk of salinization. Currently, an unsteady model is being developed to simulate drought conditions. An enhanced model will be linked with an optimization method to identify the best groundwater management schemes under drought conditions. The study highlights the applicability of a sharp interface numerical model calibrated for freshwater thickness to assess the extent of freshwater lens, saltwater intrusion, and well salinization in real islands. It can be used as the first stage analysis or when sufficient data and computational resources are not available for implementation of dispersion models.