Sustainable groundwater development in the coastal Tra Vinh province in Vietnam under saltwater intrusion and climate change

Three-dimensional transient groundwater flow and saltwater transport models were constructed to assess the impacts of groundwater abstraction and climate change on the coastal aquifer of Tra Vinh province (Vietnam). The groundwater flow model was calibrated with groundwater levels (2007–2016) measured in 13 observation wells. The saltwater transport model was compared with the spatial distribution of total dissolved solids. Model performance was evaluated by compar-ing observed and simulated groundwater levels. The projected rainfalls from two climate models (MIROC5 and CRISO Mk3.6) were subsequently used to simulate possible effects of climate changes. The simulation revealed that groundwater is currently depleted due to overabstraction. Towards the future, groundwater storage will continue to be depleted with the current abstraction regime, further worsening in the north due to saltwater intrusion from inland trapped saltwater and on the coast due to seawater intrusion. Notwithstanding, the impact from climate change may be limited, with the computed groundwater recharge from the two climate models revealing no significant change from 2017 to 2066. Three feasible mitigation scenarios were analyzed: (1) reduced groundwater abstraction by 25, 35 and 50%, (2) increased groundwater recharge by 1.5 and 2 times in the sand dunes through managed aquifer recharge (reduced abstraction will stop groundwater-level decline, while increased recharge will restore depleted storage), and (3) combining 50% abstraction reduction and 1.5 times recharge increase in sand dune areas. The results show that combined interventions of reducing abstraction and increasing recharge are necessary for sustainable groundwater resources development in Tra Vinh province.


Introduction
Groundwater is an important natural resource with high economic value and social significance. It is estimated that 2.5 billion people worldwide depend on groundwater resources to satisfy their daily needs for water and hundreds of millions of farmers rely on groundwater to sustain their livelihoods and contribute to food security (Mechlem 2012). Groundwater provides valuable services to the Mekong Delta in Vietnam. These include the supply of drinking water to millions and the prevention of saltwater intrusion (Phan et al. 2002). About 4.5 million people depend on groundwater for drinking water in the region. Seawater intrusion is among the main sources of groundwater contamination in coastal aquifers where saline water displaces or mixes with freshwater. Groundwater in coastal aquifers is often saline (Wagner et al. 2012) when the abstraction of fresh groundwater from coastal aquifers exceeds the natural recharge from rainfall or leakage from streams (Custodio 2002) when there is upconing by upward flow in the confining layers (Louw et al. 2011), and because of past marine transgressions in low-lying deltas (Delsman et al. 2017). Saltwater intrusion in coastal areas is one of the most important processes of water quality degradation by raising salinity in exceedance of drinking water standards. This is an urgent issue that must be considered and researched to offer solutions to adapt and minimize the impacts on socio-economic development and food security.
Vietnam is considered one of the most vulnerable countries with regard to climate change impact (IPCC 2014), and groundwater in Mekong Delta is sensitive to future climate change. Shrestha et al. (2016) used the projected future climate changes from two representative concentration pathway scenarios (RCP4.5 and RCP8.5) as inputs to a WETSPASS model to estimate future groundwater recharges, which were then used in a MODFLOW model to predict impacts on groundwater levels and storage. They concluded that groundwater levels and storage are projected to decline in the future due to projected decline of groundwater recharge. Silva (2018) downscaled statistically the projected daily precipitation and temperature data for Tra Vinh province from two global circulation models (CSIRO Mk3.6 and MIROC5) based on RCP4.5 emission scenario. A soil-water-balance model was used to estimate future groundwater recharge in different sandy areas. Silva (2018) concluded that groundwater recharge is projected to increase due to the increase in rainfall intensity and duration which is contradictory to the findings of Shrestha et al. (2016).
There are many modeling studies of saltwater intrusion in coastal aquifers-for example, Lin et al. (2009) studied the variable-density groundwater flow and miscible salt transport to assess the extent of seawater intrusion in the Gulf Coast aquifer of Alabama, USA. The model was calibrated with 35 observed groundwater head data in 1996 and used to predict seawater intrusion in 40 years. Sherif et al. (2012) presented the simulation of seawater intrusion in the coastal aquifer of Wadi Ham, UAE. The transient groundwater flow model was calibrated with groundwater level series from four observation wells. Possible seawater intrusion was simulated with a contact density transport model. A three-dimensional (3D) seawater intrusion model was developed to investigate the current condition of salinization by examining seawater intrusion in the Silifke-Goksu Deltaic plain, and was used to assess the effects of increasing and decreasing groundwater abstractions on the extent of seawater intrusion (Cobaner et al. 2012). The model was calibrated with one-time measurements of groundwater levels and salinity from 21 wells. Scenario analysis indicated that seawater may intrude 700 m inland after 50 years of increased abstraction. Rasmussen et al. (2013) assessed the impacts of climate change, sea level rise, and drainage canals on saltwater intrusion to a coastal aquifer with a saltwater transport model. They calibrated the model with measured heads and validated the simulated fresh/saltwater interfaces with geophysical measurements. Scenario analysis indicated that the fresh/saltwater interfaces were sensitive to groundwater recharge, sea level rise and water levels in the drainage canals. Gopinath et al. (2016) developed a numerical model to investigate the seawater intrusion in the Nagapattinam coastal aquifer of Tamilnadu, India, and found that the cone of depression may induce saltwater up-coning for a 50-year prediction. Mansour et al. (2017) developed the model to study seawater intrusion in a coastal aquifer of Karaburun Peninsula, western Turkey. The model was calibrated with one-time measurements of groundwater levels and total dissolved solids (TDS) from nine wells and used for scenario analysis. The predicted seawater intrusion may extend 420 m in 10 years. Other case studies on seawater intrusion were also conducted in different parts of the world such as in Red River flood plain aquifers, Vietnam ( With the growing population and rapid economic development over the past decades, the water demand in Tra Vinh province has increased and resulted in increasing groundwater extraction. Groundwater abstraction has increased considerably to meet water demand, which has led to the depletion of the fresh groundwater reserves. The rate of groundwater abstraction in Tra Vinh is more than 264,710 m 3 /day, and the number of abstraction wells reaches about 89,000 (Tra Vinh DONRE 2018), resulting in a persistent continuous decline in groundwater levels in almost all aquifers. Although there were modelling studies on saltwater intrusion into rivers in the Mekong Delta (Doung et al. 2018;Eslami et al. 2020), there are no published studies on seawater intrusion into coastal aquifers in Mekong Delta. This study used the integrated groundwater flow and density-dependent saltwater transport models to investigate saltwater intrusion under human activities and climate change and assessed mitigation measures for controlling groundwater abstraction and increasing (managed) aquifer recharge for the coastal Tra Vinh province. The methodology is applicable to assess seawater intrusion in other coastal aquifers.
The main objective of this study is to assess the impact of groundwater extraction and climate change on groundwater depletion in the upper aquifers of Tra Vinh province. More specifically, the study aims to construct and calibrate 3D transient groundwater flow and saltwater transport models in order to analyze the importance of groundwater abstraction vs. climate change, as well as the importance of salinization from inland (trapped) saltwater vs. modern seawater intrusion in the area. In addition, the study looks at adaptation options for controlling groundwater depletion and saltwater intrusion, through a reduction of abstraction and an increase in (managed) aquifer recharge. The findings from the modelling studies provide guidelines for sustainable groundwater resources development in Tra Vinh province.

Locations and geography
Tra Vinh province is located in the lower Mekong Delta, between two large rivers, Hau River and Co Chien River (Fig. 1). Tra Vinh province covers an area of 2,358 km 2 and has a population of 1,040,500 (General Statistics Office of Vietnam 2018). The terrain is relatively complex because it is partitioned by an entangled network of sand dunes, roads, and waterways. In general, the terrain is mainly flat with an elevation of ~0.4-1 m above mean sea level, accounting for 66% of the total area. The terrain is flatter in the north than in the south.

Climate and hydrology
Tra Vinh province is located in the tropical monsoon climate zone. Two distinct seasons are observed-the dry season lasts from December to April the following year, while the rainy season, from May to November, has an average of 115 rainy days. The annual average precipitation in the northeast is about 1,500 mm/year, while it rises to 1,700 mm/year in the southwest (Mai et al. 2014). Approximately 90% of the annual rainfall is concentrated in the rainy months (Fig. 2). Severe droughts have impacted all provinces of the Mekong Delta since the end of 2015 (UNPC Viet Nam 2016). Temperatures are relatively higher in the spring and lower in the winter with an average of 27 °C. The monthly evaporation is almost constant with an average of 80 mm/month. The difference between rainfall and evaporation is runoff and recharge, but recharge is relatively limited, so most of the rainfall generates runoff, which causes flooding problems in the wet season.
The hydrological regime of the Tra Vinh river system is complicated and controlled by the Hau River in the southwest, the Co Chien River to the northeast, the East Sea tides and canals (Fig. 3). Tra Vinh province has a dense network of canals connected to the two large rivers, Co Chien and Hau, which provide irrigation water to the fields. The main canals are distributed evenly in the province, with a density between 4 and 10 m/ha. The density of interior field channels is 45 m/ha at Tieu Can District and 18-28 m/ha at Tra Cu, Cau Ngang district (Fig. 3).

Hydrogeology
There are seven aquifers in the Tra Vinh province (Nguyen HD, unpublished report (2004), see Table 1; Phan et al. 2002;Bui et al. 2016), namely (from young to old, from shallow to deep): Holocene (qh), Upper Pleistocene 3 ) aquifers (Fig. 4). Generally, the lithology of each aquifer consists of fine to coarse sand, gravel and pebbles, separated by silts and clays. The crosssections over the area presented in Fig. 4 provide an overview of the spatial distribution and interconnection of the aquifers and aquitards. The characteristics of the aquifers are summarized in Table 2. Natural recharge by rainfall is limited to the relatively small areas with more permeable soils, where the Holocene aquifer reaches the surface, in the sand dune areas (Fig. 4).

Groundwater abstractions
Groundwater in the study area is mostly abstracted for water supply and industrial purposes. The total abstraction rate in Tra Vinh Province was estimated as 187,685 m 3 / day in 2007 (Sanh 2010) and the quantity of groundwater exploited for domestic activities increased with the yearly increase in population and water use per capita (Pham and Koontanakulvong 2019). In 2016, total groundwater abstraction in Tra Vinh province had increased by 41% to about 264,710 m 3 /day, mainly from the top three aquifers (Table 2), and with a clear dominance (67% of

Groundwater modelling framework
The general framework of this modelling study is shown in Fig. 5. First, a conceptual hydrogeological model was constructed with the groundwater modeling system (GMS) software. Then, a 3D transient groundwater flow model was constructed and calibrated with the MOD-FLOW model (McDonald and Harbaugh 1988). In a next step, a SEAWAT (Langevin et al. 2008) transport model was constructed and calibrated based on TDS measurements. Subsequently, groundwater management and climate change scenarios were formulated and implemented in the MODFLOW and SEAWAT models. Finally, the results of the scenarios were analyzed and compared in order to come up with a strategy for sustainable groundwater resources development in the Tra Vinh Province.

Construction of conceptual hydrogeological model
A conceptual modelling approach in GMS was used to create a digital conceptual hydrogeological model for Tra Vinh province. Conceptual model coverages were created using GIS utilities to represent physical structures, parameters, stresses, and boundary conditions of the aquifer systems.

Model boundary coverage
The model boundary was delineated with arcs following the topographic map of the province. The enclosed model area is about 2,247 km 2 . Figure 6 shows the model boundary, boundary conditions, and major canals simulated in the model.

Boundary conditions coverage
Referring to the boundary conditions, river boundaries were assigned to the Hau and Co Chien River and their main branches. The depth of the riverbed ranges from 5 to 10 m below the surface, cutting into the Holocene aquifer, and therefore river boundaries were assigned to the first layer of the model. General head boundary conditions were assigned below the Hau and Co Chien River for the other model  layers, at the northern administrative boundary for the upper five model layers and at the coastline also for the upper five model layers. The coastline was defined as the general head boundary for all aquifers assuming the head equal to sea level (zero) and conductance depending on the occurrence of clay layers. The aquifers might extend beneath the sea floor, but there are no data to determine the exact extension. The lateral boundaries for model layers 6-13 were assumed to be no-flow boundaries since there were no data to identify the nature of the boundaries; thus, a no-flow boundary is then considered a conservative and adequate assumption.

Model layers and elevations
The hydrogeological formations were classified into seven aquifers (2) separated by six aquitards; therefore, a total of 13 model layers were used to represent the system. The top elevation of the first aquifer was interpolated from the digital elevation model with a resolution of 90 m. The bottom elevations of the 13 layers were interpolated from 56 boreholes in the model domain. Figure 7 shows a diagram of the model structure.

Hydrogeological parameter coverages
The properties of horizontal and vertical hydraulic conductivities (Kh, Kv), specific storage (Ss), and specific yield (Sy) of the model layers were computed from pumping tests conducted between 2002 and 2016 (Dang et al. 2019). Details of pumping test wells and computed parameter values are provided in the electronic supplementary material (ESM1). Parameter zones were delineated according to similar lithologies and are shown in Fig. 8 for the upper three aquifers and Table 3 provides ranges of parameter values for the upper three aquifers.

Areal recharge coverage
The Holocene aquifer was divided into two zones with different values of groundwater recharge. The high recharge zones consist of sand and silty sand areas. The low recharge zones are clay and silty clay areas. The distribution boundaries of sand dunes are taken from the hydrogeological map.
In the study area, there are two monitoring wells in the qh aquifer, which is used to calculate the recharge for the sandy area. For the clay area, the recharge is very small, and the value was taken as 1/100 of the recharge of the sandy area. The water-table fluctuation method was used to estimate recharge the rate (Scanlon et al. 2002): where I is the infiltration rate (m/day), Sy is the specific yield, and Δϕ mod is the water-table rise at time period Δt after a rain event. In the sandy area of the Holocene aquifer, specific yield is estimated at 0.13. The measured annual water-table difference between the peak and valley was used to compute the annual recharge. The computed recharge from 2006 to 2015 varies from the lowest 120 mm/year in 2013 to the highest 197 mm/year in 2007, which corresponds to about 8 and 11% of the annual precipitation, respectively.

Abstraction coverages
There are currently about 88,930 abstraction wells with a total abstraction rate of ~264,710 m 3 /day. In particular, 117 large-capacity wells with a total extraction rate of 55,490 m 3 /day are located in the Upper-Middle Pleistocene aquifer. The remaining 88,813 wells abstract groundwater from the Holocene, Upper Pleistocene and Upper-Middle Pleistocene aquifers. The distribution of extraction wells per aquifer is shown in Fig. 9.

Observation coverages
There are 13 groundwater monitoring wells located in the aquifers qh, qp 3, qp 2-3, n 2 2 , n 2 1 and n 1 3 (Fig. 10). Most observation wells have long-term groundwater level data from late 1992 to present. The manual measurements are taken every 5 days. Some wells are installed with automatic dataloggers with hourly readings.

Dispersivity parameters
The longitudinal dispersivity was assigned a value of 50 m as one-tenth of the grid size for all aquifers and 10 m for all aquitards. The ratio of horizontal transverse dispersion to longitudinal dispersion (TRPT) was assigned a value of 0.1 and the ratio of vertical transverse dispersion to longitudinal dispersion was assigned a value of 0.01.

Seawater boundary and salinity measurements
Saltwater is represented by TDS. For the coastline, the TDS of 35,000 mg/L was specified as the concentration of seawater for all aquifers. The first estimates of TDS concentrations were inferred from geophysical surveys where saline groundwater was found in inland areas . The TDS concentrations of aquifers were subsequently determined in more detail through three electrical conductivity (EC) survey campaigns for the three upper aquifers in March 2016: 140 measuring points in the Holocene aquifer, 357 measuring points in the Upper Pleistocene aquifer and 531 measuring points in the Upper-Middle Pleistocene aquifer. The EC values were converted into TDS concentrations through a linear regression equation. The contour maps of TDS concentrations of the three aquifers (qh, qp 3 and qp 2-3 ) are shown in Fig. 11. Fresh groundwater areas have TDS values from 500 to 1,500 mg/L, the saline groundwater areas have TDS values from 1,500 to 5,000 mg/L, and seawater has a TDS of 35,000 mg/L.

Set-up of the groundwater flow model
MODFLOW was used to construct the transient groundwater flow model for Tra Vinh province. The modeling area was divided into 135 rows and 151 columns with a uniform grid size of 500 m × 500 m. In the vertical direction, 13 model layers were defined, layers 1, 3, 5, 7, 9, 11, and 13 represent aquifers, and layers 2, 4, 6, 8, 10, and 12 represent aquitards. The top semipermeable layer is combined with aquifer qh to form model layer 1.
The transient simulation covered a period of 10 years, from January 2007 to December 2016. A monthly stress period was used to simulate seasonal changes in groundwater dynamics. There are in total 120 stress periods in which all time-dependent inputs were provided.  The MODFLOW River (RIV) package was used to simulate two large rivers and major canals (Fig. 6) assigned in the top model layer. MODFLOW General Head Boundary (GHB) package was used to simulate the corresponding general head boundaries. MODFLOW Recharge (RCH) package was used to simulate recharge from the precipitation infiltration, applied to the cells of layer 1. MODFLOW Well (WEL) package was used to simulate groundwater abstractions, evenly distributed in model cells for each district, but increasing over the 10-year period as described in the preceding. MODFLOW Layer Property Flow (LPF) package was chosen to compute groundwater flow through model cell faces and changes of storage. Inputs for all these packages were mapped and transferred from the conceptual model to MODFLOW packages. Elevations of the tops and bottoms of model layers were entered from interpolated values. The groundwater flow model was calibrated by comparing computed hydraulic heads with observed ones at 13 observation wells distributed in aquifers in the Tra Vinh province. The groundwater model was calibrated by using a trial-and-error method in adjusting model parameters until simulated groundwater level corresponds reasonably with observed ones (Osei-Bonsu et al. 2004). The parameters were calibrated including hydraulic conductivity (horizontal and vertical), specific yield and specific storage. The mean residual error and root mean squared residual error (RMSE) were used to check the calibration accuracy.

Set-up of saltwater transport model
SEAWAT was used to simulate saltwater transport in Tra Vinh province. The processes include advective transport, dispersion, source/sink mixing, and density dependency. SEAWAT is coupled with MT3DMS (Langevin et al. 2008) to simulate these processes.
The MT3DMS advection package (ADV) was used to simulate advective transport. A flow interface with MODFLOW transfers groundwater velocities for advective transport simulation. The third-order total variation diminishing (TVD) scheme (ULTIMATE) was chosen to solve the advective transport. The MT3DMS dispersion package (DSP) was used to simulate dispersion caused by heterogeneity; dispersivity values were defined in the conceptual model and transferred to DSP package. The SEAWAT variable density flow (VDF) package was used to simulate density-dependent flow. The fresh groundwater has a density of 1,000 kg/m 3 ; whereas the seawater boundary, defined as the specified concentration boundary, has a TDS of 35,000 mg/L and a density of 1,025 kg/m 3 . The density/concentration slope was assigned a value of 0.000735.
The SEAWAT model cannot really be calibrated since systematic monitoring of salinity did not exist in Tra Vinh province until the implementation of the DUPC2 SALIN-PROVE project. The EC measurements from monitoring campaigns provide a reference for checking the model results. The computed TDS in March 2016 from the model should resemble the contour maps in Fig. 11.

Formulation of climate change scenarios
Future climate change may further stress the groundwater availability in the province. A warmer and drier climate is projected from climate change research in Mekong Delta (Bui et al. 2016;Shrestha et al. 2016). In this study, the scenarios of climate change were selected from a global climate change scenario RCP4.5 using downscaling techniques by two models, MIROC5 (Watanabe et al. 2010) and CSIRO Mk3.6 (Gordon et al. 2010). The rationale for choosing such models was the following: (1) they are widely used by different climate change impact studies, (2) the countries in which these models were developed have close proximity to the study area and, (3) daily predictor variables are already available for these models (Silva 2018).
In this study, only the effects of projected changes in precipitation and temperature-based potential evapotranspiration (Thornthwaite 1948) were simulated since groundwater recharge comes mainly from precipitation infiltration. Projections of precipitation and temperature from two climate change models were downscaled and validated, these are: • Scenario C 1 : Climate change model CSIRO Mk3.6 • Scenario C 2 : Climate change model MIROC5 The scenario simulation was for 50 years from January 2017 to December 2066. A monthly stress period was used. The daily groundwater recharge series from 2017 to 2066 was computed using a soil-water-balance model in Excel (Nonner and Stigter, unpublished report (2016), see Table 1). The inputs to the model include daily precipitation, daily potential evapotranspiration, soil type, crop type, root zone depth, extinction depth, and moisture content. In five representative locations, fieldwork was conducted to determine these parameters (Silva 2018). The computed daily recharge was aggregated to monthly recharge for climate change impact simulations.  Table 1). With the current increasing demand for water use and projects for securing drinking-water supply with installing new pumping wells in rural areas, groundwater abstraction will most likely increase. For better planning and management of groundwater resources in the province, the consequence of groundwater overabstraction must be analyzed along with possible mitigation measures. Therefore, based on the results of a business-as-usual scenario, three mitigation scenarios were proposed and simulated. In total, the following scenarios were simulated:

Formulation of groundwater management scenarios
• Scenario 0: Business as usual (BAU) • Scenario 1: Reduction of abstraction by 25, 35 and 50% • Scenario 2: Increase of groundwater recharge in sand dune areas by 1.5 and 2 times • Scenario 3: Combination of most feasible scenarios 1 and 2 The scenario simulation period is from 2017 to 2030, continuing from the calibration model. The objective of the BAU scenario (S 0 ) is to assess the impact of current groundwater exploitation on groundwater depletion and saline intrusion. In this scenario the amount of groundwater extraction was fixed at 264,710 m 3 /day during the simulation period, which is still a relatively conservative estimate. For scenario 1 the amount of groundwater abstraction was gradually decreased, first by 25%, then 30%, and finally 50% until to effectively reduce the groundwater level decline and reach an equilibrium between groundwater recharge and discharge. Groundwater recharge data for the period 2017-2030 for scenarios 0 and 1 were obtained from the climate change impact simulation models. The objective of increasing groundwater recharge in scenario 2 is to investigate the potential of managed aquifer recharge in the sand dune areas through methods such as contour farming and infiltration galleries (Dillon et al. 2009). The widespread occurrence of sand dunes in the southeastern part of the area provides possibilities for managed aquifer recharge of the top Holocene aquifer. Possible water sources are rainfall harvesting, storm runoff and surface water from the canals. Other managed aquifer recharge (MAR) methods can be applied in the urban areas including city rain gardens, tree boxes, street infiltration trenches and roof-top rainwater harvesting. Most feasible mitigation measures from scenarios 1 and 2 form scenario 3, which can be considered as a sustainable development plan.

Model calibration
The comparison of groundwater levels as shown in four observation wells from four aquifers in Fig. 12 indicates that the model is able to capture major variations of groundwater levels. The computed mean errors are 0.10, 0.27, 0.36, 0.19 m; and RMSEs are 0.12, 0.47, 0.56, and 0.46 m, respectively, for aquifers qh, qp 3 , qp 2-3 and n 2 2 . The strong groundwater level decline in qp 3 (Fig. 12c, d) despite limited abstractions from these aquifers reveals the good hydraulic connection with qp 2-3 as also shown in the cross-sections of Fig. 4. In contrast, the quite distinct behavior of the Holocene aquifer (Fig. 12a), without any significant trend, reveals poor hydraulic connections with the aquifers below, confirmed by the presence of a thick aquitard in large parts of the area (Fig. 4). The fittings of computed groundwater levels to observed levels in all observation wells can be found in ESM2.
The computed water budget components are plotted in Fig. 13, which shows that groundwater abstraction was increasing and was exceeding groundwater recharge. The change of groundwater storage was negative most of the time and led to a total storage depletion of about 215 × 10 6 m 3 in 10 years.
There were no historical measurements of TDS, so a formal calibration of SEAWAT model was not possible; nevertheless, the computed spatial distribution of TDS from the SEAWAT model should resemble contour maps of TDS from the sample campaigns. The contour maps of TDS values computed from the SEAWAT model in March 2016 are shown in Fig. 14. These contour maps are similar to the contours of the TDS measurements in Fig. 11, indicating that the saltwater transport model can be used for scenario analysis. Figure 15 presents the groundwater recharge and cumulative storage change time series for the two climate models. There are no significant differences in recharge and storage changes between the two climate models. The projected future recharge is stationary (no clear trend of decrease or increase) with large seasonal fluctuations. The cumulative storage change indicates rapid depletion of groundwater

BAU: business as usual (S 0 )
The flow and saltwater transport models were run from 2007 to 2030 for 24 years. The simulated results were checked by plotting the trend of groundwater level declines, expansion of the cone of depression, changes of the water budget, and saltwater intrusion. The water budget components are plotted in Fig. 16. Groundwater abstraction was kept at the same level as in 2016 with seasonal variations, and groundwater recharge for the simulation period was extended from 2017 with the recharge estimations from the climate model CSIRO Mk3.6. The continuous overabstraction caused a slight increase in river leakage and boundary inflow; however, groundwater storage was under continuous depletion. The cumulative groundwater storage depletion in 24 years was estimated as 531 × 10 6 m 3 , which is equivalent to 5.5 years of total groundwater abstraction. Inflow from the coastal boundary indicated seawater intrusion.

Reduction of abstraction (S 1 , S 2 , S 3 )
Three scenarios were simulated: S 1 , S 2 , S 3 with 25, 35, and 50% reduction of the abstraction from S 0 , respectively. The simulation results show that when the abstraction is reduced by 25%, the trend of groundwater level decline is effectively reversed (Fig. 17a, S 1 ), and the depletion of groundwater storage is stopped (Fig. 17b, S 1 ). However, inflows are not sufficient to recover the depleted storage that occurred in the first 13 years (3 years of the scenario period were kept at initial pumping rates, given the fact that these have not changed until 2020). A further Fig. 15 Time series of the a groundwater recharge and b cumulative storage change, from climate scenarios CSIRO Mk3.6 and MIROC 5 reduction of abstraction by 35% leads to a slight increase in groundwater storage (Fig. 17b, S 2 ). A large reduction of abstraction by 50% is required to recover the depleted storage by 2030 (Fig. 17b, S 3 ), when groundwater level will have risen to the level of 2007 (Fig. 17a, S 3 ).

Increase of groundwater recharge (S 4 , S 5 )
Two scenarios were simulated: the recharge rate in the sand dunes was increased by 1.5 (S 4 ) and 2 (S 5 ) times, respectively, from 2020 onwards as a result of proposed managed aquifer recharge. The current groundwater abstraction was kept. The results show that with 1.5 times increase of groundwater recharge in the sand dunes, groundwater storage is still being depleted (Fig. 17b, S 4 ). For a further increase of recharge of 2× in the sand dunes, both groundwater storage depletion (Fig. 17b, S 5 ) and groundwater level decline are stopped (Fig. 17a, S 5 ). It is clear that the increase of groundwater recharge by means of managed aquifer recharge alone will not recover the depleted storage by the year 2030.

Combined scenario (S 6 )
In this scenario, groundwater abstraction was reduced by 50% and groundwater recharge in the sand dunes was increased by 1.5 times. The results show a significant recovery of groundwater levels (Fig. 17a, S 6 ) and depleted storage (Fig. 17b, S 6 ). The sustainable development of groundwater resources can be achieved with these combined mitigating measures. Table 4 presents the estimated groundwater budget in 2016. Groundwater abstraction exceeds groundwater recharge: 28% of abstracted water comes from groundwater storage depletion. As this causes continuous decline of groundwater levels and triggers saltwater intrusion, current groundwater development, therefore, is not sustainable. This was already shown by research in coastal regions of North America (Barlow and Reichard 2010) and observed in the Indus Basin, Pakistan (Khan et al. 2008). One aspect that results in overabstraction is the limited recharge in a significant part of the aquifer, due to the extensive occurrence of clay materials. This causes a large portion of the rainfall to run off and to be drained by the canals, and is also found to occur in other deltas such as in Bangladesh (Nowreen et al. 2020). Unlike the Bengal Basin, indirect leakage from river channels to the Holocene aquifer is not significant, as incision is not deep enough to cut through the Holocene clays and create a direct hydraulic connection with the Upper Pleistocene aquifer (qp 3 ).

Comparison of scenarios
The scenario simulation and analysis indicate that both the reduction of abstraction and increase of recharge will stop the further decrease of groundwater levels; however, the reduction of abstraction is more effective to reverse the trend of groundwater level decline. A new equilibrium state will be reached when abstraction is reduced or recharge is increased. It might take longer for the aquifer system to reach a new equilibrium state so that storage change becomes zero. Allocation of groundwater abstractions from the Holocene aquifer may enhance both direct and indirect recharge, which was also shown by Nowreen et al. (2020). The low permeability is not the only factor influencing the low recharge in the area; very shallow water tables and the high drainage capacity of the system through its canals, further reduce the potential of deep groundwater recharge, even where the top sediment has moderate permeability to allow groundwater recharge.
The water authority in Tra Vinh province may opt for two strategies, with the first being to stop the continuous decline of groundwater levels and depletion of groundwater storage. For this purpose, scenarios S 2 , or S 5 , or a combination of S 1 and S 4 can be considered. The second strategy is to recover depleted groundwater storage for a long-term sustainable development. A combination of reducing abstraction (>35%) and increasing recharge (at least 1.5 times in the sand dunes are required.

Effect of climate changes
The effect of climate change was only considered in groundwater recharge from projected precipitation. The model simulation results indicated a smaller effect of climate change on groundwater depletion in comparison to groundwater abstraction. However, due to high evaporation and indirect stress on groundwater demand, the projected temperature increase might have direct consequence on groundwater recharge (Bui et al. 2016); furthermore, the possible impact from sea level rise or land subsidence on seawater intrusion was not simulated. The impact of extraction-induced land subsidence was shown by Minderhoud et al. (2020) to be very significant in the delta. A more complete study of climate change impact on groundwater resources should therefore be considered in any future study.

Effect of saltwater intrusion
Saltwater transport is a slow process and its effect can only be observed in the long term. The effect of saltwater transport can be seen from the CRISO Mk3.6 model simulation with SEAWAT. The predicted salinity distribution in 2066 shows that saline groundwater in the Holocene aquifer is freshening due to groundwater recharge, whereas salinity in the Upper Pleistocene and Upper-Middle Pleistocene aquifers increases. Figure 18 shows that saline groundwater in the northeast area is spreading and seawater is intruding to coastal zones in the Upper-Middle Pleistocene aquifer; therefore, saltwater intrusion is a long-term risk to fresh groundwater availability in deeper aquifers. About 67% of total groundwater abstraction takes place in the Upper-Middle Pleistocene aquifer, the reduction of which should be the first priority.

Conclusions
Three-dimensional transient groundwater flow and saltwater transport models were developed to assess current groundwater exploitations and to predict future changes under various interventions and climate changes in Tra Vinh province. The results revealed that groundwater has been overexploited, resulting in a continuous decline of groundwater levels, especially in Pleistocene aquifers. The present groundwater abstraction is estimated to be about 264,710 m 3 /day, while average total groundwater recharge is about 190,088 m 3 /day. A daily depletion of groundwater storage is at 74,622 m 3 /day to maintain the abstraction. To achieve sustainable groundwater development, the current groundwater abstraction must be reduced by 35%, to be about 172,726 m 3 /day. However, a better option is to increase groundwater recharge in sand dune areas. In dune areas, the natural groundwater recharge is only 10% of annual precipitation. Groundwater recharge can be increased to 25% of annual precipitation with managed aquifer recharge methods. The possible artificial recharge technologies include infiltration ponds and contour farming in sand dune areas, supplemented with street infiltration trenches, roadside infiltration trenches, rain gardens and tree boxes in urban areas. Saline groundwater is widely distributed in Tra Vinh province. Fresh groundwater is only available in the Upper and Middle Pleistocene aquifers. The results of simulation of saline groundwater distribution show that the area of saline groundwater distribution in the qp 3 and qp 2-3 aquifers increases. The increase of the saline groundwater distribution in qp 2-3 aquifer is largest, and the reason is that the amount of groundwater extraction in this aquifer amounts to 67.4% of the total abstraction. Thus, saline intrusion of the aquifers is mainly caused by groundwater abstraction. The area of saltwater intrusion occurs mainly in the northern part of the study area in the three aquifers where saltwater was trapped in the geological formations.
The projections of monthly rainfall from two climate models (MIROC5 and CRISO Mk3.6) provide similar rainfall in the future. There may be no major changes of natural groundwater recharge from precipitation. The model simulation until 2066 shows small effects from precipitation changes; however, groundwater storage depletion will continue when current groundwater abstraction is not controlled. Therefore, a reduction of groundwater abstraction is necessary in order to mitigate saltwater intrusion and to avoid land subsidence.