Paleo-hydrogeological reconstruction of the fresh-saline groundwater distribution in the Vietnamese Mekong Delta since the late Pleistocene

Study Mekong Delta, Vietnam (MKD). Study focus: Present-day distribution of fresh-saline groundwater is highly heterogeneous in the MKD. Close to the coastline, fresh groundwater is found in aquifers up to 500 m below ground surface. To gain better insight into the fresh-saline groundwater evolution since the late Pleistocene, we simulated long-term groundwater flow and salt transport in a two-dimensional NW-SW cross-section over the MKD. To fully consider the regression and transgression phases of sea-level changes over the past 60ka, variable-density groundwater flow and salt transport was simulated with SEAWAT to reproduce characteristics of the present-day distribution of fresh-saline groundwater and its age. We simulated nine scenarios to evaluate the most important factors controlling freshening and salinization processes of the MKD groundwater system. We compared the final model stage with present day observations of groundwater salinity and age. New hydrological insights for the region: The sedimentation and erosion processes, related to sea-level changes over the last 60ka, were important drivers of the fresh-saline distribution in the present MKD. The two-dimensional model indicates that most fresh groundwater in the MKD was recharged 60–12ka before present, when the sea-level was at its lowest and the top sedimentary layers had a relatively high permeability. Due to deposition of a clayey top layer during the Holocene, at present, groundwater recharge of the deeper MKD groundwater system is very limited.

In recent decades, fresh groundwater resources in deltas worldwide have reduced due to increasing groundwater extraction caused by increasing demands from domestic, industrial and agricultural sectors (Wada et al., 2010;Gleeson et al., 2012;Kuenzer and Renaud, 2012). This increased extraction rate has led to significant challenges in the sustainable management of fresh groundwater resources, especially in effectively preventing the salinization of coastal aquifers Werner et al., 2013). Both sealevel rise and greater demands on fresh groundwater resources in the future will lead to further salt water intrusion. Thus, fresh groundwater resources in coastal aquifers are now being an increasingly rapid rate of threatening. To assess the vulnerability of the groundwater resources, it is essential to have sufficient information on the salinization processes.
Vietnam is one of the countries where groundwater is an important source of water for human activities Vu et al., 2014). The Mekong Delta, Vietnam (MKD) is rapidly changing due to urbanization, land use transformation, and intensification of economic activities (Minderhoud et al., 2017), as illustrated in Fig. 1. The associated increase of freshwater demand has led to large-scale extraction of fresh groundwater with rates depleting the existing groundwater level, from 0.2 to 0.3 m per year (Bui et al., 2013;Erban et al., 2014); salinization of groundwater and surface water resources (Bui et al., 2013); and land subsidence and increased flood risk, resulting in land loss as well as damage to buildings and infrastructure (Minderhoud et al., 2017). The MKD has been highly dependent on groundwater resources since the late 1990s. Over the past few decades, groundwater extraction has increased dramatically. The total quantity of groundwater extracted in 2015 was estimated at approximately 2.5 million m 3 /day from some 500,000 wells (Wagner et al., 2012;Bui et al., 2013;Minderhoud et al., 2017).
To sustainably fresh groundwater resources, we need better insight into the functioning of the (ground)water system, the water balance terms, and the factors controlling groundwater flow and salt transport. Until now, there are no comprehensive accounts of the infiltration of fresh and/or saline groundwater in the MKD: the times that salinization and freshening have taken place, the extent that the groundwater system is now being recharged by surface water, and whether the fresh groundwater resources are being depleted. To understand the present-day fresh-saline groundwater dynamics in the delta, we first need a comprehensive representation of the MKD in its paleo-hydrogeological context. Modeling the MKD in a conceptual way will indicate where additional information is needed to yield more reliable estimates of the present groundwater situation and constructs scenarios about the impact of future climate change (Shrestha et al., 2016;Schmitt et al., 2017). To gain reliable insight into the possible processes and mechanisms that control the (ground)water balance, flow rates, and salinization processes in the MKD groundwater system, we constructed, simulated and analyzed a 2D 'Mekong-like' cross-section.
A similar approach was taken for a case study in the Netherlands (Delsman et al., 2014). A paleo-hydrogeological model was built to study the processes controlling the Holocene evolution of groundwater salinity in a Dutch deltaic coastal aquifer. We applied the same approach to the MKD. Our aims were: (i) to make an inventory of the most likely explanations for the present-day distribution of Changes (e.g. urbanization, sea-level rise) and consequences (e.g. subsidence, salinization) in the Mekong Delta, Vietnam (after Delsman, 2015). fresh-saline groundwater; (ii) to determine to what extent fresh groundwater is being recharged at present; and (iii) to identify the sensitive factors controlling the salinization and freshening of the MKD groundwater system.
In this research, a 2D paleo-hydrogeological model was applied in a representative MKD coastal aquifer system to consider the processes controlling the evolution of groundwater salinity over the entire Holocene and late Pleistocene. We aimed to understand sedimentation and erosion processes related to the deposition of the clayey top layer. We also aimed to understand how the sea-level changes affect the fresh-saline groundwater distribution in the MKD.

Hydrogeology of the Mekong Delta, Vietnam
The MKD is one of the largest deltas in the world, lying at the mouth of the complex Mekong River network (Tanabe et al., 2003). It covers a total land area of about 39,700 km² and is located in the southern part of Vietnam, bounded by the Gulf of Thailand to the southwest, the East Sea to the south and southeast, and Cambodia to the north (Fig. 2a) (Bui et al., 2013;Shrestha et al., 2016). The delta surface is relatively flat, with low elevation ranging from 0.3 to 2.0 m above mean sea-level (MSL), except for a few hills in the southwest in An Giang and Kien Giang provinces (Minderhoud et al., 2017).
Presently, the MKD is classified as a wave-influenced and tide-dominated delta (Ta et al., 2005). Sediments were deposited during repeated transgression and regression events from the late Neogene (Miocene, Pliocene) up to the present Holocene period, which has resulted in a highly heterogeneous stratigraphy of the MKD basin fill. The Cenozoic rifting phases, related to the opening of the East Sea of Vietnam approximately 32-15.4 million years (Ma) ago, resulted in a NW-SE tectonic trough (Nguyen et al., 2004;Bui et al., 2013). This trough is filled with sediments ranging in age from the early Miocene (approximately 23.8-16.4 Ma ago) to the present-  (Nguyen et al., 2004).

Table 1
Lithological descriptions and the hydraulic properties of the studied area. The measurement permeability of aquifers varies, these data were derived from (1) aquifer tests, taken from Bui et al. (2013) and (2) the simulated mean value of both aquifers and aquitards, taken from Minderhoud et al. (2017). day (Anderson, 1978;Nguyen et al., 2004). The bulk of the sediments in the proximal part of the basin are of fluvial and estuarine origin, but a substantial part of the delta, especially the Miocene and Pliocene sediments, is of marine origin. These deposits reach a total present-day thickness of up to 1000 m in the SE part of the delta in the Tra Cu district (Tra Vinh province) situated between the Bassac and Mekong rivers ( Fig. 2a) (Nguyen et al., 2004;Bui et al., 2013;Doan et al., 2016). Tectonic movements, in particular the NW-SE pattern of normal faults (e.g. the fault line stretching along the Bassac river valley), have likely played an important role in the formation of the fluvial patterns in the delta (Nguyen et al., 2004). The hydrogeological cross-sections ( Fig. 2b and c) and data from borehole logs (Nguyen et al., 2004;Kuenzer and Renaud, 2012;Bui et al., 2013) provide an overview of the spatial distribution and interconnection of the age stratification systems in the delta's subsurface. The sediments can be roughly subdivided into seven formations ( Fig. 2b and c) (Nguyen et al., 2004;Bui et al., 2013;Doan et al., 2016). Each formation comprises two parts: the aquitards, the upper part, which is composed of low permeability silt, clay, or silty clay; and the aquifers, the underlying part, which is more permeable and consists of fine to coarse sand. The seven aquitards and seven aquifers, from surface to deeper level, are Holocene (Holo), Late Pleistocene (LPlei), Middle-Late Pleistocene (MLPlei), Early Pleistocene (EPlei), Middle Pliocene (MPlio), Early Pliocene (EPlio), and Late-Miocene (LMio) (Fig. 2b and c) (Nguyen et al., 2004;Bui et al., 2013;Doan et al., 2016;Shrestha et al., 2016;Minderhoud et al., 2017). The lithological nature of the substratum and its hydraulic properties are shown in detail in Table 1. The Vietnam Petrol Corporation determined the offshore stratigraphy and found Neogene-Q formations from 550 m to 650 m thickness that are composed of sand, clay and siltstone interleaved. The Upper Miocene level from 450 m to 820 m thick is composed of sand and clay interleaved. Boehmer and Kootstra (2000) presented groundwater levels in some aquifers of the MKD. The groundwater level data were derived from the average pre-1990 groundwater level elevations in each aquifer from the clusters of monitoring wells. All groundwater level maps show similar groundwater level patterns and similar groundwater flow patterns. Groundwater levels in all aquifers in the coastal low-lying areas were between 0 m and 2 m above MSL pre-1990. In the southern area of the MKD, they were usually less than 1 m above MSL and only occasionally below mean sea-level. This means that before 1990, regional groundwater flow was pretty limited on a regional scale, while groundwater extractions were absent (Bui et al., 2013).

The age of groundwater and salinity data in the MKD
The age of groundwater (apparent age) (Suckow, 2014) and salinity data in the MKD were provided by the Division of Water Resources Planning and Investigation for the South of Vietnam (Louvat and Ho, 1989;Nguyen et al., 2004Nguyen et al., , 2005. Stable isotope and tritium analyses were performed in the IAEA Isotope Hydrology Laboratory and the 14 C content of the TDIC was measured at the Center of Nuclear Techniques in Ho Chi Minh City (Louvat and Ho, 1989). Suckow (2014) suggested that we take care when using the term groundwater age. He presented three basic definitions of the age of groundwater, and we use his definition of apparent age. Fig. 3c shows the locations where the age of groundwater and salinity data were available. Fig. 3c also shows the measurement positions relative to our 2D cross-section (red line in Fig. 2a) and the age of groundwater and salinity data points for the MKD. We projected all the observed data to the 2D cross-section (red line in Fig. 2a and c) approximately perpendicular to the present coastline. Some of the projected samples are some distance from the cross-section, but they still belong to the Mekong Delta, Vietnam basin, and therefore are valuable for our cross-section study.
The distribution of fresh-saline groundwater is very heterogeneous due to sea-level regression and transgression phases over the last few millennia in the MKD and is not clearly related to depth (Fig. 3a). The horizontal axis in Fig. 3a and b is the red line in Fig. 3c, with point X II , Y II on the left. Close to the present coastline, fresh groundwater was found at more than 450 m below ground surface level. All the groundwater samples that measure age and salinity, see Table 2 and Fig. 3, were derived from observation wells belonging to the national groundwater monitoring network in Vietnam (from the Division for Water Resource Planning and Investigation for the South of Vietnam). The age of groundwater observations (Fig. 3b) were collected from Vietnamese reports on the age of groundwater 14 C dating (Louvat and Ho, 1989). These 14 C dated groundwater samples were tested at the Vietnam Atomic Energy Institute. Analytical errors of the age of groundwater were determined to be from 1.0 to 8.2 ka (Nguyen et al., 2004). Well depths varied from 42 to 462 m below the surface at six different aquifers ( Table 2). Most of the observed ages varied from 5 to 18 ka before present (BP) and lay at less than 300 m below MSL. At depths below 300 m, most of the observed ages were relatively old, even more than 25 ka BP.

Paleo-hydrogeology of the Mekong Delta, Vietnam
Analysis of the salinization and freshening processes in such a large region requires a paleo-hydrogeological description over a long timescale (Meisler et al., 1984;Gossel et al., 2010). The present large deltaic system was formed over the past ten million years (Nguyen et al., 2004;Ta et al., 2002;Truong et al., 2011). Groundwater flow and salt transport must be modeled over a period long enough to capture the phenomena that led to the present-day distribution of fresh and saline groundwater (e.g. Delsman et al., 2014). Using the age of groundwater data and research indicating that salt transport in porous media on a delta scale is a very slow process (Meisler et al., 1984;Post, 2004;Larsen et al., 2017), we simulated the salinization and freshening processes over a period of 60 ka.
In this long period, three main phenomena must be taken into account. The first is the changing hydraulic conditions and, in particular, sea-level changes (Fig. 4a). The second is the groundwater interaction with the surface water system, and the third is the Table 2 Overview of the observed groundwater ages (in ka) and salinity (in TDS g/L) and their aquifer, lithology of the layer surrounding the screen, depth below ground surface, depth of the screen, water level (in meters relative to present mean sea level) in the Mekong Delta, Vietnam (Nguyen et al., 2004). See text for abbreviations of Late Pleistocene (LPlei), etc. sedimentation and erosion that has taken place during the last 60 ka (Fig. 4b). We approximate the continuous changes by a number of subsequent stress periods. The stress periods were chosen to capture both the changes in the sedimentation/erosion process (in particular during the Holocene and a part of the Late Pleistocene) as well as the changes in the hydraulic boundary conditions. With respect to the sea-level changes, we used insights from Waelbroeck et al. (2002). Their method is based on a transfer function between sea-level and benthic Δᵹ 18 O, estimations. We divided the past 60 ka into eight distinct stress periods (SP, from SP1 to SP8, see Fig. 4a and Table 3). SP1-SP3 are phases of sea-level decline; SP4-SP6 are phases of sea-level rise; and SP7-SP8 are phases of slight sea-level decline. During SP1, from 60 to 40 ka BP, sea-level was estimated to be 45-80 m below present MSL. During SP2, from 40 to 26 ka BP, sea-level declined from 80 to 130 m below present MSL. SP3, from 26 to 18 ka BP, covers the Last Glacial Maximum (LGM) (Hanebuth et al., 2000;Waelbroeck et al., 2002). The LGM led to a deep valley (where the deepest incision is 68 m below present MSL) being eroded in the subsurface of the eastern branches of the present Mekong river. Erosion during the LGM likely removed a substantial amount of the near-surface deposits, estimated at tens of meters (Ta et al., 2005(Ta et al., , 2002Truong et al., 2011) (Fig. 4b).
During SP1 and SP2 we presume that a kind of braid plain existed that deposited fluvial sands on a part of the MKD. Interfluves may have contained clay and silty sands. During SP3, the LGM phase, deep erosion occurred that removed part of the sediments of the delta surface deposits. During SP4, from 18 to 12 ka BP, the valley was filled up with estuarine and tidal river sands and these were covered by fine silty sands and marsh deposits (Nguyen et al., 2000;Ta et al., 2005). The sea-level at that time rose from 130 m to 60 m below the present MSL (Hanebuth et al., 2000). During SP5, from 12 to 6.5 ka BP, the sea-level rose from 60 m to the present MSL. The estuarine sediments in the deep valleys grade into open-bay mud facies at depths of 35-20 m below the present MSL, with ages around 8-6 ka (Nguyen et al., 2000). These deposits covered a wide area stretching over the inland fluvial plain and consisted of clay, silt and muddy sediments (Nguyen et al., 2000). During SP6, from 6.5 to 4.5 ka BP, the sea-level reached its highest position (4 m above present MSL): the area up to the border with Cambodia was then fully submerged by the sea, and the delta advanced over 200 km towards the northwestern part of the MKD. During SP7, from 4.5 to 2.5 ka BP, and SP8, from 2.5 BP to the present-day, the sea-level declined stepwise to the present mean sea-level (Pirazzoli, 1996). The river floodplain, supratidal and intertidal areas, gradually became silted up with fine grained (clay and silt) clastic material including abundant organic matter in the inland delta plain, resulting in the Holocene aquitard which is characterized by a high resistance ( Fig. 2b and c).

Location cross-section
The location of the 2D model, being a representative 'Mekong-like' cross-section, is visualized by the red line in Fig. 2a. It is fairly similar to the representative NE-SW cross-section A-A', which also has an orientation perpendicular to the present coastline. We used the computer code SEAWAT (Langevin et al., 2008) to simulate a variable-density groundwater flow and coupled salt transport across the 2D model. We further modeled the age of groundwater directly (Goode, 1996) by including a first-order irreversible reaction as an additional factor (Zheng, 2009(Zheng, , 2010. This was set to zero when entering the model domain, and then one year was added for every year of flow inside the model domain. In addition, data on the age of groundwater suggested that the modeling period should be long enough to include the late Pleistocene and the entire Holocene. The simulation time-step was 10 years. This was sufficient to meet the stability criteria as required within the MT3DMS/SEAWAT solute transport modeling (Langevin et al., 2008;Zheng and Wang, 1999), and provided accurate results. Each stress period (out of eight) consisted of a different model simulation, with a different geological setting (Fig. 6), but the model state (concentration and age) at the end of each stress period was used as the starting state for the  Waelbroeck et al., 2002). The eight socalled stress periods (from SP1 to SP8) are shown by red lines (length of line represents the duration of a stress period and the height of the sea level during that period). (b) Erosion and sedimentation processes estimated during the last 60 ka BP in the Mekong Delta. The paleo-land surfaces in erosional and sedimentation settings were reconstructed based on average altitude of the list of 14 C ages from cores taken from the Mekong River Delta (Nguyen et al., 2000;Ta et al., 2002Ta et al., , 2001aTruong et al., 2011). See also Table 3.
Hung Van Pham, et al. Journal of Hydrology: Regional Studies 23 (2019) 100594 subsequent period. The 2D steady-state flow model domain has a total length of 530 km, including 30 km in Cambodia, 200 km of the MKD, and 300 km extending into the Vietnam East Sea (Fig. 2b). The horizontal grid size was 2000 m. We used 50 vertical layers with a thickness of 10 m each, providing 500 m of total thickness for the vertical domain. The surface elevation, which is the upper boundary of the flow domain, has changed over the past 60 ka due to the process of erosion and sedimentation (Table 3) (Ta et al., 2005;Truong et al., 2011).
SEAWAT (Langevin et al., 2008) was used to solve the following partial differential equation for variable density groundwater flow: where ρ is the density of the groundwater (M L −3 ); K 0 is the hydraulic conductivity tensor (L T −1 ); h f is the freshwater head (L); z is the vertical coordinate (L); ρ f is the density of fresh groundwater (M L −3 ); S s is the specific storage coefficient (L −1 ); t is the time (T); ϕ is the effective porosity (-); C is the concentration (M L −3 ); ρ ss is the density of the sink or source (T −1 ); and q ss is the sink and source term (T −1 ). It was assumed that K 0 is isotropic, that the density of fresh groundwater ρ is 1000 kg m −3 and that the maximum density of the water ρ ss is 1025 kg m −3 . A linear equation of state was used to calculate the density of the groundwater from concentration, where it was assumed that only salt concentration influences the density: Relative salt concentrations were used, ranging from 0 (concentration of fresh groundwater) to 35 (concentration of sea water); hence, C equals 0.7143. The governing equation for salt transport in the simulations was: Table 3 List of 14 C ages (right part in the column Period length) from the cores taken from the Mekong Delta, Vietnam. Collected average altitude is included to adapt the elevation of the top layers accordingly. VLM1, BT2, BT3, VL1, TV1, TC1, DT1, CD1 are core locations taken from the Mekong Delta, Vietnam (Ta et al., 2001a(Ta et al., ,b, 2005Truong et al., 2011). These numbers present the altitudes of the cores which were taken from the boreholes at the field in the Mekong Delta, Vietnam. The depth of the cores is adapted to mean sea level.
where D is the hydrodynamic dispersion tensor (L 2 T −1 ); q is the specific discharge vector (L T −1 ) and C ss is the source and sink concentration (M L −3 ). The groundwater flow equation and the dispersion term of the salt transport equation were formulated using a cell-centered finite difference scheme. The groundwater flow equation and the non-advective part of the salt transport equation were solved using the Pre-Conditioned Conjugate Gradient (PCG) and Generalized Conjugate Gradient (CGC) packages, respectively (Langevin et al., 2008;Zheng and Wang, 1999). The advection term was solved using third-order Total-Variation Diminishing TVD method, based on the ULTIMATE algorithm (Universal Limiter for Transient Interpolation Modeling of the Advective Transport Equations) (Zheng and Wang, 1999). Using this method, and splitting up the advection-dispersion equation into two components, we envisage that the solution is mass conservative, without excessive numerical dispersion, and essentially oscillation-free (Zheng and Wang, 1999). We also tested on our modeling case the performance of two other often used solvers within the SEAWAT solver suite, viz. the Finite-Difference (FD) and the Method of Characteristics (MOC) methods. The solvers HMOC and MMOC were not considered as they usually underperform. For MOC, we used two parameter sets: in the accurate case we used NPL = 16, NPH = 45, NPMIN = 45, NPMAX = 100 instead of 4, 15, 15, 75, respectively; see Zheng and Wang (1999). The results of all the solvers appear to be all similar, and thus, the TVD method was used throughout the entire modeling.

Geological dynamics
In addition to these processes, compaction results in a decrease of the permeability of top layers. We adapted the elevation and permeability of the top layers accordingly (Ta et al., 2002;Truong et al., 2011) (Tables 3 and 4). The horizontal hydraulic conductivity for each stress period, in each aquifer and aquitard, is presented in Table 4. The Holocene sediments were not present in SP1-SP5; the sediments were deposited in the last 6.5 ka (Table 4) (Ta et al., 2001a,b;Truong et al., 2011). It is likely that the late Pleistocene layers mostly disappeared due to erosion during SP1 to SP3, but after SP4 the rising sea-level led to further sedimentation ( Table 3).
The estimates of aquifers and aquitards properties rely on geological processes. Clay compaction in particular will reduce the hydraulic conductivity. Furthermore, other relevant variables (e.g. recharge, effective porosity) have to be estimated over the whole simulation period. The lack of data indicates that the conceptual geological model and the values of the hydraulic properties and other variables are disputable. We do not claim to deduce the only correct geology and parameter values, but by analyzing several options that change the hydraulic conductivities of aquitard layers, we were able to establish a hydrogeological model that we can use to perform a plausible analysis of the most sensitive factors. The starting points are the two geological cross-sections given in Fig. 2b and c. For groundwater transport phenomena, the discontinuities and heterogeneity of the layers (in particular the aquitards) at a smaller scale are important. To analyze the discontinuity and heterogeneity, we simulated four 2D steady-state flow model representations with aquifers and aquitards including different hydraulic conductivities. This was done with the aim to find the reference model that fitted best observed salinization data. Four geology models were used. Firstly, geology (i) has more disconnected strata and thinner aquitards (Fig. 5a) with hydraulic properties as shown in Table 4. Secondly, geology (ii) has more connected strata and thicker aquitards (Fig. 5b) with hydraulic properties as shown in Table 4. Thirdly, geology (iii) has more disconnected strata and thinner aquitards (being labeled as 'clay', 'silt', and 'sandy clay', in Fig. 5a) with all horizontal hydraulic conductivities of aquitards 10 times lower than geology (i). And finally, geology (iv) has more disconnected strata and thinner aquitards ( Fig. 5a) with all horizontal hydraulic conductivities of the aquitards equal to 10 −9 m/s. Note that all four options have the same hydraulic properties for very fine sand, fine sand, and coarse sand aquifers (Table 4).
In the Mekong Delta, the hydrogeological conceptual model consists of seven aquifers separated by seven aquitards (Bui et al., 2013;Minderhoud et al., 2017). The aquifers and aquitards are heterogeneous and anisotropic. The hydraulic conductivities are divided into parameter zones. The top aquifer is unconfined and the other aquifers are confined. The impermeable basement consists of intrusive and extrusive formations.
The hydraulic properties (Table 4) are based mainly on recent research from the Division for Water Resource Planning and Investigation for the South of Viet Nam (DWRPIS). These data were derived from many aquifer tests conducted by DWRPIS (Bui et al., 2013). From 1981 to 2015, there were more than 234 aquifer tests performed at the MKD (Bui et al., 2013). The values of the horizontal conductivities were determined for seven aquifers and resulted in a range of 2 × 10 −6 -8 × 10 −4 m/s (Bui et al., 2013;Ha et al., 2015;Minderhoud et al., 2017). There are no tests currently available for the vertical conductivities of the aquifers or for either the horizontal or vertical conductivities of the aquitards. The vertical hydraulic properties for the eight stress periods were assumed to be 1/3 and 1/1 of the horizontal hydraulic properties of aquifers and aquitards, respectively (Bagarello et al., 2009;Minderhoud et al., 2017;Wang and Qiu, 2017). For anisotropy, we used the same values as in Minderhoud et al. (2017), in which a 3D model of the Mekong delta was calibrated over 25 years. The effective porosity was assumed to be 0.25, representing both sand and clay layers for all eight periods.

Boundary conditions
We assumed Dirichlet (constant head) boundary conditions on the southeastern side and the northwestern side of the 2D model. The constant head boundary at the Cambodian border was assigned the same values of the land surface geometry at that time. The southeastern boundary at the coastline was based on the estimated average MSL at that time (Figs. 4 and 6).
The river head boundary in the model was specified relative to the surface elevation at that time. This was 5.0 m (SP1-SP4), 3.0 m   (SP5-SP7) and 0.2 m (SP8) lower than the surface elevation at the Cambodian border at that time and from there, it dropped smoothly over the 2D cross-section towards the Vietnamese East Sea at the corresponding time (Fig. 6). The model had only a limited river head boundary for the three stress periods SP5 to SP7 due to full inundation from a rising sea-level. The formation depths are based on the 14 C ages from the cores taken from the Mekong Delta, Vietnam (Table 3) (Ta et al., 2005(Ta et al., , 2002Truong et al., 2011) which provided the average altitude in the modeling in each stress period.

Salinity conditions
We used Total Dissolved Solids (TDS, in g/L) in our modeling to represent salinity. As mean sea-level at 60 ka BP was approximately 70 m below MSL, the initial fresh-saline groundwater distribution was assumed to be saline (35 TDS g/L, similar to present TDS of sea water) for all model layers more than 70 m below the present MSL. Layers 70 m above MSL were assumed to be filled with fresh groundwater with 0.1 TDS g/L. The TDS concentration at the southeastern boundary was kept constant at 35 TDS g/L during the past 60 ka, while at the northwestern boundary it was kept at a constant level of 0.1 TDS g/L. The TDS of the river and recharge remained constant with a sink/source concentration value of 0.2 TDS g/L and 0.1 TDS g/L, respectively, for the whole simulation period (Fig. 6). For simplicity, we assumed that there was no salt water intrusion through the river system, although salt water intrusion in deltaic areas was possibly important.
The molecular diffusion coefficient (D m ) for porous media was assumed to be 8.64 × 10 −5 m 2 d −1 (e.g. Delsman et al., 2014or Oude Essink et al., 2010. The longitudinal dispersion coefficient α L was set to 2.0 m, and the transversal ratio to longitudinal dispersion coefficient α L /α TH and α L /α TV was set to 0.1.

Groundwater recharge
The present-day water balance for the MKD that links the inflow components to the aquifer system includes infiltration of rainfall and irrigation; downward leakage through the semi-permeable aquitards; seepage from rivers, streams and lakes; downward flow from the overlying aquifers through the separate aquitards; and subsurface inflow over the boundaries of the domain. The discharge components include evapo-transpiration of shallow groundwater; well extraction from the aquifers; flow from the MKD into the sea; flow from the aquifers into the rivers, canals, streams, and lakes; and flow from the aquifers into the extensive network of surface drains (Bui et al., 2013;Boehmer and Kootstra, 2000).
There is no long-term precipitation record available for the MKD over the past 60 ka, however Chabangborn et al. (2014) estimated that in the Asian monsoon area precipitation minus evaporation was in the range of 400-800 mm/year in the Last Glacial Maximum (viz. 23-19 ka BP). We applied a constant, uniform groundwater recharge of 1.5 mm/day, which is equal to about 25% of the current long-term precipitation average, for the first four stress periods (SP1-SP4). There is only a limited recharge boundary for the next three stress periods (SP5-SP7) due to the sea water inundation of the whole MKD. Due to the increase of the resistance of the top aquitards in SP8, natural groundwater recharge was expected to be fairly small. Groundwater recharge was chosen to be 0.2 mm/ day only. In addition, the river conductance changes over time. During the first four stress periods erosion occurs, whereas during the last four periods sedimentation dominates. As a result, the river conductance in the first four stress periods (0.1-1.0 m 2 /day) is smaller than in SP8 (0.4-1.0 m 2 /day) (Table 4). Hung Van Pham, et al. Journal of Hydrology: Regional Studies 23 (2019) 100594

Setting up the simulations
The purpose of the 2D model was to gain insight into the most likely explanations for the present-day fresh-saline groundwater distribution and to identify the factors controlling the salinization and freshening processes. The first step was to set up the reference case, based on our best knowledge of the most likely parameter values. We adjusted the model in such a way that the distribution of fresh-saline groundwater, as well as the age of the groundwater after computing 60 ka, compared reasonably well with present-day observations. Simultaneously, groundwater heads, groundwater velocity, hydraulic parameters, recharge, and discharge had to lie within a plausible range. During the development of the reference case, we were able to identify the factors that were most sensitive Hung Van Pham, et al. Journal of Hydrology: Regional Studies 23 (2019) 100594 to the model results.
Based on the knowledge gained from the reference case, we defined nine scenarios with a minimum and maximum value for each of the identified factors. We changed only one parameter at a time so that the effects of the change could be observed individually. The results of the simulations of the nine scenarios were evaluated and compared to the current distribution of saline and fresh groundwater. Table 5 summarizes the reference case and the nine scenarios. The most sensitive parameters were expected to be hydraulic conductivities, surface water level, recharge, river conductance, longitudinal dispersion, and initial TDS concentration. Scenarios 1-4 tested the sensitivity of both horizontal and vertical hydraulic conductivities, focusing specifically on each set of aquifers, aquitards, top aquitard layers, and the percentage of vertical sandy shortcuts, all for geology (i) (Fig. 5a and Table 4). For scenarios 1-3 the scenarios with the maximum parameter values are labeled "a"-scenarios; the scenarios with the minimum parameter values are labeled "b"-scenarios. The hydraulic conductivities in the maximum "a"-scenarios were 10 times larger than the reference case, while in the minimum "b"-scenarios they were 10 times smaller than in the reference case (except in scenario 1 where they were 4 times larger or smaller). The top layers in scenarios 3a and 3b were continuous. In scenario 4a, around 25% of sandy, permeable shortcuts were present in the top layer, while in scenario 4b, 35% of permeable shortcuts were present. In scenario 5, in comparison to the reference case, the river level rose by +5 m (scenario 5a) and dropped by -5 m (scenario 5b). Furthermore, in scenarios 6a and 6b, the recharges were 20% higher and 20% lower, respectively, than in the reference case. Scenarios 7 and 8 focused on the conductance of the river and the longitudinal dispersion, respectively. They were multiplied ("a"-scenarios) or divided ("b"-scenarios) by a factor of 10, as in scenarios 2 and 3. In scenario 9, we changed the initial concentration to completely fresh groundwater in order to compare with the reference case and evaluate the effect caused by a different initial concentration.

Reference case
The modeled distribution of fresh-saline groundwater is presented in Fig. 7a-i at 60 ,40,26,18,12,6.5,4.5,2.5 ka BP and at present (i.e. the end of stress periods), respectively. The age distribution of the groundwater simulated at the present time is given in Fig. 8. For comparison, the projection of observed salinity and the age of groundwater are given in Fig. 7i (present-day) and Fig. 8, respectively.
Four main phases can be distinguished (Fig. 7a-i): during the first phase, from 60 ka to 18 ka BP, there was freshening of the groundwater, driven by the sea-level regression and low mean sea-level. Fresh groundwater replaced saline groundwater up to 500 m below the surface. During the second phase, the largest volume of fresh groundwater was present between 18 and 12 ka BP.
During the third phase from 12 ka to 2.5 ka BP, a serious salinization of the groundwater occurred because a higher mean sealevel led to sea water infiltrating the groundwater. Kooi et al. (2000) explained that if the rate of transgression is sufficient fast, the salinization of the subsurface cannot keep up by means of a lateral 'Henry Case' salt water intrusion wedge development, and saline water will move on top of the (fresh) groundwater system, causing salt water to infiltrate. Here in our case, the conditions are such that even salt water plumes (also often called seawater fingers, salt fingers or density fingers) are bringing relatively rapidly salt groundwater into the system. This fingering process is a much more rapid salinization process than diffusion is (see Simmons et al., 1999;Zimmermann et al., 2006;Wooding et al., 1997a,b). The driving force of the groundwater flow and salt transport system is the transition of a relatively "high-gradient system" into a "low-gradient system" within a relatively short time span at this period due to the rapid rise in the sea water level. However, we have to remember that our model area of interest, the horizontal scale (a few hundreds of kilometers) is quite a few orders of magnitude larger than the vertical scale (some tens of meters), making the "highgradient system" not really a topography-driven groundwater flow system. During this third phase, saline groundwater was present in most parts of the 2D model.
During the last phase of the last 2.5 ka, sea-level dropped slightly, thus a freshening process is to be expected, as in the early stress  Fig. 7. Groundwater salinity distribution (in TDS g/L) modeled: (a-i) for the collected reference case geology (i) at 60, 40, 26, 18, 12, 6.5, 4.5, 2.5 ka BP and for the present-day; (j-m) for geology (ii) at 40, 12, 4.5 ka and for the present-day; (n-q) for geology (iii) at 40, 12, 4.5 ka and for the present-day; (r-u) for geology (iv) at 40, 12, 4.5 ka and for the present-day. The salinity scale is red for saline (35 TDS g/L) and blue for fresh groundwater. The dots in Fig. 7i are locations where salinity was observed in the MKD (Nguyen et al., 2004. periods. However, due to the high resistance of the top sediment layer, formed during the period of high mean sea-level, there was no significant freshening of the groundwater. To evaluate the feasibility of our reference case in detail, we compared the model results at the end of SP8 with the observed salinity data (Table 2 and Fig. 3). We compared the groundwater salinity and age simulated in SP8 with the salinity and age observed at different depths in the MKD (Figs. 7i and 8). Furthermore, we ensured that the groundwater heads lay within a feasible range of values during the entire calculation period. For geology (ii) (more connected strata and thicker aquitards (Fig. 5b) with hydraulic properties as shown in Table 4), geology (iii) (more disconnected strata and thinner aquitards with all horizontal hydraulic conductivities of aquitards 10 times lower than geology (i), and geology (iv) (more disconnected strata and thinner aquitards with all horizontal hydraulic conductivities of the aquitards equal to 10 −9 m/s), the freshwater from the top system (recharge source) could not penetrate deeper than about −350 m ( Fig. 7j-u) below present-day surface level after a period of 60 ka. This does not match with observed groundwater salinity because, presently, fresh groundwater is found in aquifers up to 500 m below the surface (Fig. 3a). In addition, at the present-day (viz. Fig. 7m, q, u, representing the geologies (ii), (iii) and (iv), respectively) sea water could not have infiltrated into the fresh groundwater system to form the present heterogeneous distribution of fresh-saline groundwater in the study area.
The effect of aquitards in the top system on the development of the salt water plumes depends on their occurrences and their level of permeability per geology model (Figs. 5,7 and Table 4). When the geology model contains clay layers with hydraulic conductivities that are assumed to be normal in this MKD (viz. geology (i)), then wide salt water plumes occur (Figs. 5a, 7f-i). When the geology model contains thick clay layers with (very) low hydraulic conductivities (e.g. geology (ii)), Fig. 5b), salt water plumes moving through the medium were slowed down; hence, it took more time to get into the deeper parts of the groundwater system while less plumes are developed anyway (Fig. 7m) (Kooi et al., 2000). Even when the geology set-up remains the same as in the reference case (Fig. 5a), but the existing hydraulic conductivities of the aquitards are (much) lower (e.g. geology (iii), (iv)), the salt water plumes were also slowed down and the deeper parts of the (onshore) groundwater system would not contain saline groundwater (Fig. 7q, and u). In both geology models with very low permeability of the aquitards, the total amount of salt water plumes was large but the number of deep plumes was pretty limited. For further analysis, we used the geology (i) option as the most likely geology for the reference case ( Fig. 7a-i).
Because we compared a 2D modeled cross-section to 3D real-world observations, we cannot expect absolute equivalence. Nevertheless, the spatial distribution and heterogeneity seen in Figs. 7a-i and 8 appear to be reasonable. In particular, in the final stage, fresh groundwater could be modeled in the deepest aquifers, up to 500 m below the surface and at a short distance from the present coastline. This situation was observed in Tra Vinh and Soc Trang provinces (Bui et al., 2013). Furthermore, the modeled age of groundwater was also consistent with the observed data (Figs. 3b and 8). As shown in Fig. 8, groundwater older than 40 ka was only present in the offshore section of the 2D model. For the inland section, the groundwater was not older than 40 ka BP, except for small volumes in the deepest aquifers close to the bedrock, where it is trapped as fairly stagnant groundwater. This indicates that most of the groundwater in our 2D model was completely replaced during the past 40-60 ka (Fig. 7a-i). Likely, older (> 60ka) saline groundwater is trapped in areas with low permeability in the vertical and lateral direction, and has not been flushed out as suggested by the results of our 2D reference model. In the 3D reality, flushing could more easily happen and this saline water could be better drained by slow mixing or even just only hydrodynamic diffusion. This phenomenon is quite common in old saline groundwater systems and could also explain the heterogeneous distribution of groundwater ages.  (Nguyen et al., 2004).

Scenarios
Scenarios 1 (change in aquifers), 2 (change in aquitards), 3 (change in top aquitards), 7 (change in river conductance), 8 (change in longitudinal dispersion) and 9 (change in initial concentration) appear to be the most sensitive to groundwater salinity and age. The results of the other scenarios appear to be less affected. The present-day distributions of groundwater salinity for the nine scenarios are given in Fig. 9. The present-day distributions of the age of groundwater for these scenarios are given in Fig. 10. Fig. 11 shows the deepest location of fresh groundwater, as a function of time for all scenarios (except for scenario 9). Fig. 12 shows the total volume of fresh groundwater (in million m 3 per stretched meter), as a function of time for all scenarios. Each scenario is shown separately with the reference case and the minimum and maximum values (except for scenario 4). Fig. 13 shows the comparison between the volume of fresh groundwater modeled for the reference case with mean sea-level relative to its fluctuations over the last 60 ka; note that we inverted the y-axis for this volume of fresh groundwater on purpose to better show the correlation between fresh groundwater volume and mean sea-level.
For the reference case, the deepest level of fresh groundwater as well as the volume of fresh groundwater changed dramatically with the sea-level over the last 60 ka (Figs. 11 and 12). Even when the sea-level rose, the deepest level of fresh groundwater kept Fig. 9. Modeled distribution of groundwater salinity for the present-day reference case and for nine scenarios. The salinity scale is red = for saline (35 TDS g/L) and blue for fresh groundwater.
Hung Van Pham, et al. Journal of Hydrology: Regional Studies 23 (2019) 100594 getting deeper and the volume of fresh groundwater remained stable for some millennia before it fell following the mean sea-level (Fig. 13). Fig. 7d and e show where freshening continues between 18ka and 12ka in the center of the study area while saltwater transgresses from the Southeast boundary. The delay is mainly connected to the permeability of the low-permeable layers (Kooi et al., 2000;Post and Kooi, 2003). The volume of fresh groundwater has continued falling, even to the present-day, thus there is apparently no steady-state situation.  Our first four scenarios focus on the conductivities of the sedimentary layers. In scenarios 1a and 2a the groundwater velocity was so high that saline groundwater was quickly replaced by fresh groundwater recharge in the entire model within 60-30 ka BP. Thus, the conductivities of both aquitards and aquifers -as high as 10 and 4 times the reference case, respectively -were probably not truly representative for these scenarios. In contrast, the groundwater in scenarios 1b and 2b (Fig. 12) moved a bit too slowly. This observation was also supported by the age of groundwater modeling: while the ages of groundwater were very low in scenarios 1a and 2a (Fig. 10a and c), they were too high in scenarios 1b and 2b ( Fig. 10b and d). In scenarios 1b and 2b, fresh groundwater could  . Our 2D modeling approach shows that the fresh groundwater volume follows the pattern of Mean Sea Level but thousands of years delayed. Therefore an autonomous further decrease in fresh groundwater volume is to be expected. not penetrate deeper than 350 m below the present-day surface level after a period of 60 ka BP. The high vertical resistance that was set in scenario 2b meant that the infiltration of fresh or saline groundwater from the top layer to deeper aquifers was also too slow (Fig. 9d). The changes in top aquitards (scenario 3) basically only affect the volume and depth of fresh groundwater during the last thousand years. Scenarios 4 (more sandy shortcuts in top layer), 5 (river water level) and 6 (groundwater recharge) appear to be less sensitive to the factors we identified than the others. Scenario 7a (increase the river conductance by a factor of 10) shows that most of the fresh groundwater came from upstream Cambodia and the river system in the higher northwestern parts, causing the simulated age of groundwater to shift from young to old in a northwest to southeast direction (Fig. 10g); this was not consistent with the age data observed. In general, the age of groundwater in the Mekong Delta, Vietnam increases into the direction of groundwater flow. The age of groundwater grows from young groundwater to older groundwater as we move from the upper to the lower aquifers (Figs. 3b and 8). Scenario 7b, in which groundwater was rapidly and completely replaced between 60 and 40 ka BP, is also not plausible. In scenario 8a, the dispersion coefficient was ten times higher than the reference case, resulting in too much mixing of the groundwater, while in scenario 8b fresh groundwater could not penetrate deeper than 460 m below the present-day surface level. Over time in scenarios 8a and 8b, the volume of fresh groundwater was notably different, as observed in Figs. 11h and 12h. The reference case and scenario 9 represent the most extreme initial concentrations (almost completely saline and completely fresh, respectively).

Discussion
Our 2D model of the Mekong delta, Vietnam, simulates variable-density groundwater flow and coupled salt transport over the past 60 ka. The final model stage, after 60 ka of simulation and representing the present situation, shows acceptable agreement between modeled and the observed present-day distribution of both groundwater salinity and age. Analyses of all nine scenarios give strong indications that the groundwater present at 60 ka BP in the entire system has been almost completely replaced within 60 ka. We therefore assume that our simulation period of 60 ka is long enough to avoid any strong or undesired effects from the unknown initial conditions. Moreover, scenarios 1, 2 and 3 revealed our first insight into the rather limited value of the conductivity of the entire groundwater system: low vertical conductivity of the (top) aquitards (scenario 2b) would have led to the deeper parts of the MKD having a substantially larger fresh groundwater volume. This would have been due to the infiltration of fresh groundwater taking place over a sufficiently long period; the vertical conductivity would not have been a serious limiting factor over this time scale.
For the last ten thousand years, a low vertical conductivity would have safeguarded the fresh groundwater resources. Unfortunately, the most realistic vertical conductivity (as used in our reference case) could not protect the fresh groundwater resources from being rapid infiltrated by saline seawater over time. The analysis of the geology (ii) model indicates that connected and homogeneous aquitards are most unlikely. The freshwater from the top system could not penetrate deeper than about −350 m (Fig. 7k) below present-day surface level. Nevertheless, saline water could also not infiltrate into the groundwater system in order to form the very heterogeneous present-day fresh-saline distribution (Fig. 7m). In order to reproduce the present-day age and salinity distribution, water from the upper boundary has to infiltrate into the deeper subsoil. This appears to be possible only with disconnected and/or heterogeneous aquitards.
There is a significant change in fresh groundwater volume during the 60 ka of simulation (see Fig. 13). On changing fresh groundwater volumes and sea-level elevation, it appears the delay to the sea-level rise by volume is on average several thousands of years. The delay is mainly linked to the permeability of the low-permeable layers, though the shallower the system we are focusing on, the shorter the 'memory' of the groundwater system. In the deepest parts of the MKD, the memory is many tens of thousands of years. Additionally, Fig. 13 shows that during the period 60-12 ka BP, most freshwater infiltrated the groundwater system when the sea-level was very low. Substantial parts of the top layer had high hydraulic conductivity, thereby leading to a strong freshwater infiltration.
It appears from our scenarios that the most sensitive parameters are the hydraulic conductivity of the aquifers and aquitards (scenarios 1, 2 and 3), the interaction with the riverbed (scenario 7), which represents the interaction between the surface water system, the longitudinal dispersion (scenario 8) and the initial concentration (scenarios 0 and 9). As can be seen in scenarios 1, 2, 3, 7, 8 and 9 (Figs. 9-12), these factors strongly control the salinization and freshening processes in the 2D model of the MKD. Although the depth of the deepest fresh groundwater and the volume of fresh groundwater change significantly between the minimum and maximum values modeled, we observed that only the reference case can realistically reproduce both the observed data on salinity and the age of groundwater. Furthermore, the longitudinal dispersion coefficient was observed to be important for the degree of mixing, but its effect on the overall availability of freshwater is not substantial (Fig. 12h).
Sedimentation took place during the last transgression period (the last 12 ka). This resulted in top clay layers that have a high resistance to infiltration by fresh groundwater. Currently, the present top layer acts more or less as a seal on the subsurface reservoir. The deeper aquifers are basically separated from surface hydrology due to the low permeability of the top aquitard, but the deeper aquifers still receive some groundwater influx via the upper aquifers. Nevertheless, our 2D model result suggested that the natural groundwater recharge of fresh water in the present-day is limited. As a result, the recharge of fresh groundwater is very limited at present and the volume of fresh groundwater was seen to decrease significantly over the last 12 ka in all the scenarios (Fig. 13). Despite the slightly falling sea-level, almost all the scenarios show a persistent decrease in the volume of fresh groundwater during the last 2.5 ka. This is mainly due to the limited groundwater recharge and over exploitation of groundwater resources is serious problem. It should be noted that we did not include groundwater extraction in the model, which would lead to a much faster decrease in the fresh groundwater resources. Because the extractions are only relevant for the last decades, they are not relevant for our scenarios.
We performed minimum-maximum scenarios for several parameters (hydraulic conductivity, river conductance, groundwater recharge) but did not attempt to quantify the uncertainties in our 2D model. As such, the results from the scenarios form a conceptual analysis rather than a real-world simulation or quantification of the associated uncertainties. Nevertheless, the salinity and age in the reference case resemble the observed present-day distribution of both groundwater salinity and age reasonably well and give us confidence that our model has captured the most important groundwater processes in the MKD. Informed by this research, our future work includes compiling a 3D groundwater model, using 3D geology, groundwater-surface water interaction, salinity distribution and groundwater extractions. The present research contributes to this future 3D modeling tool that can be used to inform (ground)water strategies regarding vulnerability of the system.

Conclusions
We have simulated the evolution of the fresh-saline groundwater distribution and its age through the late Pleistocene and Holocene periods along a 2D cross-section with a variable-density groundwater flow and coupled salt transport model based on the subsurface and geographical characteristics of the Mekong Delta, Vietnam. This model enabled us to analyze the fresh-saline distribution, volume of fresh groundwater and age dynamics of the groundwater reservoir. The erosion and sedimentation processes over the past 60 ka, as well as mean sea-level changes, still appear to have a dominant influence on the present-day distribution of fresh-saline groundwater and its age. Most fresh groundwater infiltrated the subsurface during a period of low sea-level, at least more than 20 ka ago. Sedimentation during the last 20 ka has increased the resistance of the top layer and this currently prevents fresh groundwater from being recharged to a great extent. From the last period of lowstand of mean sea-level, our simulations strongly suggest that the volume of fresh groundwater was consistently declining (under natural conditions). Extraction of groundwater -not included in the model over paleo time scales -will amplify this decline. Despite the limitations of our 2D model, the most sensitive factors controlling salinization and freshening processes were identified. Our analysis is the first contribution towards a full 3D groundwater model that can take variable-density groundwater flow and coupled salt transport into account.