Understanding Saltwater Origins and Mechanisms in the Coastal Aquifers of Da Nang Area (Central Vietnam)

1. Introduction

Currently, the world’s population is grappling with the pressing issue of water scarcity and saltwater intrusion. The availability of water has become a major concern for many regions, especially those that are arid or semiarid. In the past, the management of water resources has typically treated surface water and groundwater as distinct and separate from each other [1], even though they are closely interconnected components of the hydrological cycle [2]. In Vietnam, most of the water resources management is focused on surface water [3,4], while issues concerning groundwater are neglected and not managed properly. In many studies, modelling has been identified to be the most effective approach for estimating both surface and groundwater availability [5,6,7,8].

Saltwater origin and mechanisms in coastal aquifers in Vietnam have been studied extensively and the main causes identified are groundwater exploitation [9]; inter-aquifer leakage due to diffusion and density flow or high hydraulic gradient [10].

A major urban center located in the central region of Vietnam is Da Nang with over 70 km of coastline and is renowned for its pivotal role in tourism, industry, and agriculture. Despite this, the city is facing multiple challenges such as high population density, water scarcity, environmental pollution, and saltwater intrusion. Due to its vast coastline, the aquifers in this region are greatly impacted by the saline intrusion, but no study has yet analyzed this phenomenon using numerical models. A survey was conducted to create a hydrogeological map of the city, and it was revealed that about 30% of the aquifer area in this region is affected by salinity [11]. Thus, it is essential to evaluate the current condition of saltwater intrusion and identify an effective solution to this problem.

In this study, geophysical, hydrogeochemical and isotopic methods were applied to investigate the distribution of saltwater and its origin and mechanism in the coastal aquifers. The status of saltwater intrusion in the Da Nang area was assessed using a coupled SWAT-MODFLOW numerical model. The model was calibrated and employed to forecast potential scenarios of saltwater intrusion, considering the influence of economic and urban development through changes in groundwater pumping and recharge rates.

2. Study Area

The focus of this study is the coastal region of Da Nang, as highlighted in Figure 1a. The area spans a total of 360 km². The Han River (Figure 1a) acts as the boundary between the coastal plain and the southern parts of the city, while sand dunes can be found in the Ngu Hanh Son and Nam O areas, largely due to wind activity. The geology and hydrogeology of the study region drawn for a cross-section A-B and C-D is shown in Figure 1b. and comprises three primary hydrogeological formations: aquitards, porous aquifers, and fracture aquifers.

The study area is predominantly covered by porous sediments from the Holocene and Pleistocene ages (a blue color in Figure 1). Holocene deposits in the study region are found with thickness varies from 5 m to 29 m and contain sand, silty sand and gravel (light blue, Figure 1). The area of Holocene aquifer (qh) is around 160 km², with an outcrop area of approximately 90 km². The groundwater levels in this aquifer vary seasonally, and it is in the range of 0.45 to 3.0m from the surface in Lien Chieu Hai Chau area, 1.8 to 11.2m in the Son Tra – Ngu Hanh Son area.

Following Holocene, there are fine to coarse sand, as well as sand with gravel and pebbles of Pleistocene age with a thickness up to 34m, average is about 25m. The sediment in Pleistocene deposit contains gravel sand to silty sand (dark blue, Figure 1). The area of the qp aquifer is around 160 km².

An observation for the water level in different boreholes and geological structures revealed that qh and qp aquifers are the same hydraulic system in the study area.

The fissure aquifer in the study zone possesses an area of less than 40 km² (green color in Figure 1), majority of those areas encompassed by younger sediments, predominantly composed of sericite schist, calcareous limestone, quartz schist, and shale, the aquifer is considered to have a moderate to poor capacity for water storage.

In light of the hydrogeological characteristics, the Holocene and Pleistocene porous aquifers need particular attention because water from these aquifers is currently exploited as the clean water source of the local population. Both, however, are under the vulnerability of saline intrusion.

3. Methods

3.1. Sampling Procedure and Field Measurement

Water from 35 boreholes was sampled during two sampling campaigns in Dec/2019 (dry season) and May/2020 (wet season) for major cations and anions, TDS and isotopic compositions of deuterium and oxygen 18 in water determination. The positions of boreholes from those samples were taken were selected around the saltwater boundaries from hydrogeological mapping in 2012 (Figure 1). The isotopic compositions of δ¹⁸O  δ²H in water samples are expressed in delta (δ) notation by the Equation (1) and Equation (2):

$$\delta {}^{2}H=\left(\frac{{}^2{R}_{sample}}{{}^2{R}_{std}}-1\right)x1000$$

(1)

$$\delta {}^{2}O=\left(\frac{{}^2{R}_{sample}}{{}^2{R}_{satd}}-1\right)x1000$$

(2)

Where ²R_sample, ²R_std, ¹⁸R_sample and ¹⁸R_std are, respectively, the ratios of abundances (in Mole) of the heavier to those of the lighter isotopes in water samples subject to analysis and in standard. The standard for water isotopic composition determination is VSMOW2 (Vienna Standard of Mean Ocean Water batch 2), supplied by the International Atomic Energy Agency (IAEA) based in Vienna, Austria.

Before taking groundwater samples, the groundwater levels in boreholes were recorded, then the stagnant water in each sampling borehole was pumped out till the pH, temperature, EC in water became unchanged. The pH, temperature and EC of water were measured using a WTW Multi 197i device equipped with appropriate probes.

Groundwater samples were taken from the boreholes and passed through polycarbonate membranes with a pore size of 0.45 µm to eliminate suspended solids. For major cations analysis, 50 mL of the filtered samples were acidified by HNO₃ acid (65%, PA grade, Merck, Germany) to pH 1-2. The samples were kept in HDPE (high-density polyethylene) bottles for transport to the laboratory. For major anions determination, 50 mL of the filtered samples were stored without acidification in HDPE bottles for transport to the laboratory. For water isotopic composition analysis 50 mL of the filtered water was taken and stored in HDPE bottles with double caps to avoid evaporation. The samples were subjected to analyses within 2 weeks of the sampling day.

Rainwater was collected using a device made according to the IAEA’s recommendation [12]. The device was installed on the roof of an office building located within the study area. Sea water was sampled at a distance of 20 m from the shore, at a depth of 0.5 m from the surface for δ¹⁸O and δ²H analyses.

3.2. Analytical Procedure

The major ionic contents of water samples were analyzed at the Centre for Excellence in Analysis and Experiment (CEAE) of the Hanoi University of Mining and Geology (HUMG) by using ion chromatography (IC) method on a DIONEX 600.

The water isotopes compositions (δ¹⁸O  δ²H) were analyzed by laser spectrometry method on a Picaro’s cavity ring down spectrometer CRDS L2130-I at the Laboratory of Water Resources of the Ha Noi University of Natural Resources and Environment (HUNRE). The results of δ¹⁸O  δ²H were expressed in ‰ relative to the VSMOW2.

3.3. Geophysical Method

Electrical resistivity tomography (ERT) was used in this study to validate the developed saltwater numerical model. The ERT technique provides a continuous characterization of subsurface electrical conductivity (by the means of electrical resistivity distribution) [13,14]. Besides, the ERT method provided high-quality electrical resistivity data obtained [15] with substantial spatial precision and the capability to survey underground areas continuously in both 2D and 3D dimensions [16] that were used to calibrate the numerical model in this study.

In the research area, 17 ERT profiles across the salt-fresh water interface were measured (Figure 1a) by using a multi-electrode SuperSting R1/IP equipment, with inter-electrode spacing of 10m.

3.4. Numerical Model

MODFLOW is a comprehensive groundwater modeling system that employs finite-difference technology to simulate three-dimensional flow. It can accurately model steady and variable flows in saturated environments and has the ability to represent aquifers as confined, unconfined, or a combination of both [2]. The model considers different boundary conditions such as unchanging pressure head, recharge, and consistent or varying flows. The discretized domain of the aquifer system is created using a combination of nodes and finite difference cells [17,18]. The governing equation for MODFLOW is based on Darcy’s Law, which is represented by a partial differential equation, and it is described in equation (2). With its powerful capabilities and advanced technology, MODFLOW is a highly effective tool for analyzing groundwater systems and their behavior in various environments.

$$\frac{\partial }{\partial x}\left(K_{xx}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yy}\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_{zz}\frac{\partial h}{\partial z}\right)-W=S_s\frac{\partial h}{\partial t}$$

(3)

The behaviour of groundwater is governed by a fundamental equation expressed in terms of partial differential equations. The Equation (3) takes into account several variables such as hydraulic conductivities in the x, y, and z axis (K_xx, K_yy, and K_zz), piezometric head (h), volumetric flux per unit volume (W), specific storage of the porous medium (S_s), and time (t). The hydraulic conductivities in the major axes, which are parallel to the x, y, and z axis, are used to calculate the rate of flow in each direction. The height of the water column above a reference point, known as the piezometric head, is a crucial factor for understanding the pressure distribution in the aquifer. Sources or sinks of water that may affect the flow direction and rate are represented by the volumetric flux per unit volume. The specific storage of the porous medium indicates the quantity of water stored in a unit volume of the medium under unit hydraulic gradient change. Time plays a crucial role in modelling the behaviour of groundwater over time. Thus, Equation (3) provides a comprehensive framework for understanding the complex behaviour of groundwater in various geological settings.

3.4. SWAT-MODFLOW Integration

The combination of surface water and groundwater models can provide valuable insights into the spatial-temporal replenishment of groundwater and the interaction between river water and groundwater. SWAT-MODFLOW is a powerful tool that can simulate these complex interactions [6,19]. Figure 2a illustrates that SWAT is used to model the top layers such as the root zone, vadose zone, and shallow aquifer, while MODFLOW is used to model the lower layer, which is the deep aquifer.

In this study, both the SWAT and MODFLOW models were employed to evaluate recharge rates, first separately and then in combination. The SWAT model was initially used to evaluate recharge rates and calculate groundwater recharge values at the HRU level. In order to combine the SWAT and MODFLOW models, the HRU recharge rates had to be transferred and utilized as input data for MODFLOW, as illustrated in Figure 2b. Since SWAT is semi-distributed, it was uncertain where each HRU was located within the sub-basins. To overcome this issue, a single HRU was created for each sub-basin based on the dominant soil, land use, and slope option [2]. By combining the two models and accounting for the limitations of each, it was possible to obtain a more accurate and comprehensive understanding of the groundwater replenishment process.

4. Results and Discussion

4.1. Distribution of Salinity

The classification of water types according to TDS (total dissolved solids) is an important step in water resource management. Zhou [20] proposed a classification system based on TDS levels, which includes four categories: brine (TDS > 50 g/l), saline water (TDS = 3–50 g/l), brackish water (TDS = 1–3 g/l), and freshwater (TDS < 1 g/l). This system allows for a clear understanding of the quality of water resources and helps to guide decision-making related to their use and management. An understanding of the TDS levels present in various water sources can help in assessing the suitability of water for different applications, such as drinking, irrigation, or industrial use. Data from 35 sampling locations were used to generate a limit between fresh and salt water. Among points that were collected during 2020 (in dry and rainy seasons), TDS varies from 0.1 to 3.6 g/l. Five samples with TDS higher than 1 g/l are located along the rivers (Figure 1.a). It was revealed that a sharp contrast of electrical resistivity along Profile 1, 2 (Figure 3). Below surface up to 15 m depth, resistivity varies from about 15 Ohm.m close to the river and up to more than 50 Ohm.m is observed corresponding to superficial Holocene sediments made of sandy clay or unsaturated zone. Resistivity readings ranging from 1 to 10 Ohm.m at a depth of 15-30 m suggest the presence of salt and brackish water from the nearby river. Aquifer with freshwater is reflected by a high resistivity (> 40 Ohm.m) below 30 m depth due to the presence of gravel. The low resistivity zone in the deep part close to the river is demonstrated as salty water affected by surface water from the river.

4.2. The Chemical Composition of the Groundwater in the Area under Investigation

In the study of groundwater hydrochemistry, it is essential to examine the variations between rainy and dry seasons. As shown in Table 1, the hydrochemical and stable isotopic results of groundwater in the rainy season were analyzed. Similarly, Table 2 displays the corresponding results for the dry season. These findings provide valuable insights into the dynamics of groundwater quality and its potential impact on the environment and human health.

The chemical characteristics of groundwater in the study area were presented in a Piper diagram shown in Figure 4.

It was found that the groundwater in the area is dominated by Na-Ca-Cl-SO₄ type. However, some boreholes were found to have chloride concentrations that exceeded the national standard for groundwater quality set by the Ministry of Natural Resources and Environment [21] by up to 40 times. The high levels of chloride in groundwater may be attributed to the ongoing salt intrusion occurring in the study area.

4.3. Isotopic Compositions in Water as Indicators for the Genesis and Recharge Area of Groundwater in the Study Region

In this study, the isotopic composition of groundwater sampled from the study area is displayed alongside the local meteoric water line (LMWL) for the Da Nang city area in Figure 5. The LMWL in this region follows a model with a slope of 8.98 and an intercept of 10.71, as depicted in Figure 5. These values are slightly different from the Global Meteoric Water Line (GMWL) of 8 and 10, respectively, as established by Craig [22]. The variation in values could be attributed to the unique weather conditions in Da Nang city, which is characterized by high atmospheric temperatures and humidity throughout the year.

In the rainy season isotopic compositions in groundwater (open dots in Figure 5) from the study area positioned close to the LMWL indicate that groundwater in this area is recharged from the local precipitation. The only point positioned below the LMWL, possibly is being affected by seawater making the isotopic signatures in water from that borehole to be enriched.

For the dry season, a quite different trend of the δ²H vs δ¹⁸O relationship for groundwater was observed (solid dots in Figure 5), oxygen 18 in the water was depleted but deuterium was enriched. The reason for this could be the recharge area was at a higher altitude for which isotopic compositions in precipitation were depleted, and the flow direction of the recharge would be through and beneath a landfill. In fact, a large landfill is currently located around 12-15 km towards Northwest of Da Nang for collecting and processing waste from the city. The leachate from the landfill containing sulphide (H₂S) and methane (CH₄) as the products of domestic waste bio-decomposition leak into aquifers. This leakage would be so intensive during the dry season because of not enough pressure from recharge to protect from the infiltration of the leachate to aquifers. It was reported that hydrogen in H₂S and CH₄ generated from domestic waste is enriched in deuterium [23,24] and the enriched deuterium in sulphide and methane participated in the isotopes exchange with hydrogen in water in aquifers making deuterium in water to be enriched. However, oxygen 18 has no chance to take part in isotopes exchange, therefore its signature was unchanged as it was from the recharge area.

4.4. Saltwater Intrusion in the Study Area

Saltwater intrusion modelling was validated and calibrated by using EC measurement data from boreholes in two campaigns 2019-2020 and electrical resistivity tomography (ERT) measurement. The results from the calibrated numerical model in Holocene and Pleistocene aquifers in the rainy season (05/2020) show that the salinity area with TDS higher than 1 g/l in the Holocene aquifer decreases faster than the Pleistocene aquifer compared to which in the rainy season of the year 2012. The salinity decreases because of the precipitation recharge to these aquifers. The replenishment of the aquifer due to rain and the pumping process both impact the aquifer, causing it to pale. In areas located near the Han and Cu De River, and the coastline, the salinity level rises with the passing of time, a result of the high influence of rivers and the sea on the surface aquifer. Figure 7 indicates areas with TDS greater than 1g/l of Holocene and Pleistocene aquifers during wet season (4/2020). Accordingly, the total area with salinities higher than 1 g/l in Holocene aquifer is 41.8km² and in Pleistocene aquifer is 53.4km².

Figure 6. Simulated and observed groundwater head values at observation borehole QT3A (a) and QT3B (b).

Figure 6. Simulated and observed groundwater head values at observation borehole QT3A (a) and QT3B (b).

4.5. Saltwater Intrusion due to Recharge and Pumping in Coastal Aquifers

The potential control of saltwater intrusion by reducing groundwater pumping or by increasing surface water infiltration was assessed using the calibrated model. To predict the saltwater intrusion in 2020, four distinct scenarios were simulated:

• status quo: with the current groundwater pumping rate and recharge from rain/surface water
• stop pumping: use 100% surface water resource for water supply
• double pumping rates: increase GW pumping rate to fit the demand of drinking water of the city
• no recharge occurs: no recharge area due to urban development

Figure 6 displays a comparison between the simulated water heads in scenarios 1 and 2 and the water heads that were acquired from the existing pumping rates. The changes in the water heads ranged from 0.2 to 1 m, and they depended on the position of the observation borehole. Those boreholes located in proximity to pumping wells showed a higher degree of change. Based on the results of the numerical modeling for saltwater intrusion in the Holocene and Pleistocene aquifers in April 2020, it was observed that the area of the Holocene aquifer with initial high salinity decreased at a faster rate than that of the Pleistocene aquifer. This difference may be due to the higher levels of precipitation recharge in the Pleistocene aquifer, as shown in Figure 7.

Predicted saltwater areas in different scenarios for both aquifers Holocene and Pleistocene in the years 2035 and 2050 are shown in Table 3. Most of the pumping wells in this area penetrated to Pleistocene aquifer, the simulated model shows the saltwater area in both aquifers decrease quickly when all pumping wells stop, the saltwater area decreases 23.4% and 28.9% in Pleistocene and Holocene aquifer respectively in 2050. The simulated salt intrusion model from the scenario of groundwater of pumping rate doubled to fit the water demand in future times is relatively similar to the scenario of no recharge to groundwater due to urban development.

5. Conclusions

In Da Nang area, data from groundwater analysis and ERT measurements shows the distribution of saltwater (TDS>1g/l) in Holocene with about 60.3km² and in the Pleistocene aquifer with 41.3km2. In some boreholes, water contains chloride at a concentration exceeding 40 times compared to the National Standard for groundwater quality. Hydrochemical and stable isotopic data of 35 points from both aquifers demonstrated that in the rain season, groundwater is recharged from the local precipitation, in the dry season groundwater is recharged from a higher altitude area and the flow direction of the recharge would be through and beneath a landfill. The replenishment of groundwater from surface sources is a valuable resource that can reduce the extent of saltwater intrusion in many areas. Along the river downstream up to the hydraulic dams, saltwater in groundwater mainly recharges by surface water from rivers.

The study area was simulated for surface and groundwater, including saltwater intrusion, through the development of an integrated SWAT and MODFLOW model. Thanks to groundwater pumping management, with freshwater recharge from the surface, the area of saline groundwater in both aquifers is decreasing. Double pumping rate scenario and no recharge scenario give the same salt intrusion effect to qh and qp aquifers. In this area, groundwater is an effective resource for water supply to the city, but detailed studies are needed to control saltwater intrusion into aquifers.