Geochemistry and evolution of groundwater resources in the context of salinization and freshening in the southernmost Mekong Delta, Vietnam Journal of Hydrology: Regional Studies

Study Ca Mau Province (CMP), Mekong Delta (MD), Vietnam. Study focus: Groundwater from deep aquifers is the most reliable source of freshwater in the MD but extensive overexploitation in the last decades led to the drop of hydraulic heads and negative environmental impacts. Therefore, a comprehensive groundwater investigation was conducted to evaluate its composition in the context of Quaternary marine transgression and regression cycles, geochemical processes as well as groundwater extraction. New hydrological insights for the region: The abundance of groundwater of Na-HCO 3 type and distinct ion ratios, such as Na + /Cl - , indicate extensive freshwater intrusion in an initially saline hydrogeological system, with decreasing intensity from upper Pleistocene to deeper Miocene aquifers, most likely during the last marine regression phase 60 – 12 ka BP. Deviations from the conservative mixing line between the two endmembers seawater and freshwater are attributed to ion-exchange processes on mineral surfaces, making ion ratios in combination with a customized water type analysis a useful tool to distinguish between salinization and freshening processes. Elevated salinity in some areas is attributed to HCO 3- generation by organic matter decomposition in marine sediments rather than to seawater intrusion. Nevertheless, a few randomly distributed locations show strong evidence of recent salinization in an early stage, which may be caused by the downwards migration of saline Holocene groundwater through natural and anthropogenic pathways into deep aquifers.


Introduction
Over the past decade much attention has been paid to environmental challenges in the world's most densely populated mega deltas (Bucx et al., 2014). Population growth, urbanization, climate change and the fast development of industry and agriculture exponentially increased the demand for freshwater and put the already limited resources under pressure (Bucx et al., 2014;Fares, 2016;Lam et al., 2018;Minderhoud et al., 2017). In many deltas, surface water is saline and polluted (Nguyen et al., 2019;Rahman et al., 2019), and rainwater is only available during the rainy season. Therefore, groundwater is the major and most reliable source of freshwater, which consequently leads to heavy overexploitation accompanied by various negative environmental impacts (Erban et al., 2014;Ferguson and Gleeson, 2012;Minderhoud et al., 2020Minderhoud et al., , 2017Werner et al., 2013). This also applies for the Vietnamese Mekong Delta (MD), the third largest delta on earth, which has around 17.3 million inhabitants (GSOV, 2019) and relies heavily on fresh groundwater Gunnink et al., 2021;Ha et al., 2018). The most prominent negative environmental impact associated with groundwater overexploitation in the MD is the decrease of hydraulic heads and the potential of seawater intrusion, which would substantially degrade groundwater quality, limit its use and affect human health (Erban et al., 2014;Xiao et al., 2021;Rahman et al., 2019;Sohrabi et al., 2021). In addition, severe local land subsidence of several centimeters a year, which is induced by decreased pore pressure due to pumping, results in irreversible compaction of unconsolidated sandy aquifers and highly compressible clay layers (Erban et al., 2014(Erban et al., , 2013Karlsrud et al., 2020;Minderhoud et al., 2015;Parker, 2020;Thoang and Giao, 2015). In combination with additive processes, the low-lying MD is threatened in its very existence. Those include natural compaction of the upper Holocene sediment layers (Zoccarato et al., 2018), a decline of fluvial sediment supply transported by the Mekong River due to upstream dam constructions (Allison et al., 2017;Karlsrud et al., 2017;Kummu and Varis, 2007), salinization of surface water (Ferguson and Gleeson, 2012;Gugliotta et al., 2017) as well as sea-level rise due to global climate change (Smajgl et al., 2015;Woodroffe et al., 2007). Some parts of the MD are additionally facing heavy coastal and riverbank erosion (Anthony et al., 2015;Karlsrud et al., 2017;Liu et al., 2017;Luijendijk et al., 2018;Tamura et al., 2020). Weather extremes and natural hazards, including floods and prolonging dry seasons leading to severe droughts, like in early 2020, further threaten the area and potentially increase the demand for fresh groundwater (Moser et al., 2012;Phan et al., 2020). In coastal areas the consequences of groundwater extraction often surpass the impact of sea-level rise (Erban et al., 2014;Ferguson and Gleeson, 2012). Overall, the problems can be summarized as a progressive "loss of land and freshwater", which will eventually lead to inundation of large parts of the delta (Minderhoud et al., 2020). As this was acknowledged in the last years by researchers and decision makers, initiatives like the Mekong Delta Plan and Basin Development Strategy Plan by the Mekong Delta Commission aim to solve abovementioned challenges, reduce groundwater overexploitation, and work towards a sustainable development of the MD (MDP, 2013;MRC, 2016).
In order to forecast the future development of groundwater resources in the MD, it is essential to understand the hydrogeological system and the main drivers leading to present groundwater geochemistry. Pham et al. (2019) showed in a regional fresh-saline groundwater distribution model over the whole MD that freshwater present in deep aquifers probably infiltrated during the marine regression phase during the last glacial maximum (LGM) 26-12 ka BP. However, in this previous approach, only total dissolved solids (TDS) and groundwater age were considered and no further insights on groundwater geochemistry were available. Detailed analysis, such as ion ratios and the assessment of water types provide the opportunity to describe the state of the aquifer and to distinguish between geochemical processes (Giménez-Forcada, 2010;Russak and Sivan, 2010;Amiri et al., 2016). Therefore, this study focuses on groundwater evolution and dynamics from a geochemical perspective. A comprehensive groundwater investigation was conducted in CMP, which was chosen as a representative study area for the southern MD to study salinity patterns in a spatial and temporal context. CMP is the southernmost province in Vietnam located on the Ca Mau Peninsula in the MD. It is surrounded by three sides to the sea and is therefore prone to seawater intrusion. It is important to identify the status of the aquifers and to delineate geochemical processes. For example, direct seawater intrusion would lead to accelerated salinization rates, while mixing of fresh and saline groundwater is much slower (Vengosh, 2013). Therefore, this study has direct implications for the whole MD region and contributes to a better understanding of the groundwater system and offers new insights to overcome groundwater related challenges in the region.

Materials and methods
2.1. Study area 2.1.1. Overview CMP is the southernmost province in Vietnam and covers an area of 5221 km 2 (GSOV, 2019). It is part of the Ca Mau Peninsula in the MD, which includes neighboring provinces Kien Giang in the north and Bac Lieu in the northeast. Fig. 1 depicts the location of the study area in South-East-Asia and its main features. CMP is surrounded by the West Sea (Gulf of Thailand) to the west and by the East Sea to the south and east. While the delta has an extremely low average elevation level below 2.0 m above mean sea level (MSL), the average elevation level of CMP is estimated to be less than 1.0 m (Minderhoud et al., 2019;Zoccarato et al., 2018). In the tropical monsoon climate of the MD, the rain season in CMP lasts from May to November, which is followed by a pronounced dry season from December to April. CMP receives the highest precipitation in the MD due to the southwest monsoon, which can reach 2500 mm/year. The temperature ranges between 23 • C and 32 • C and relative humidity is constantly high at 76-86%. (Cosslett and Cosslett, 2014;GSOV, 2019).
The region inhabits around 1.2 million people, whereas around 300,000 live in the provincial capital Ca Mau City (CMC) (GSOV, 2019). A large proportion of the population lives in rural areas outside the smaller district capitals along the dense network of rivers and canals and are working in the agricultural sector. CMP has a special role in the agricultural sector of Vietnam, especially for shrimp farming (Cosslett and Cosslett, 2014;GSOV, 2019;Thong et al., 2010). Therefore, land use is dominated by brackish saltwater aquaculture and one-crop rice farming, sometimes alternating depending on the season and water availability, and rarely two-crop rice cultivation (Le et al., 2018;Nguyen et al., 2020).

Hydrogeology
In the late Neogene, comprising Miocene and Pliocene epochs, the former Mekong basin, a northwest-trending tectonic trough, started to subside due to the uplift of the Himalaya. High erosion rates resulted in a high sediment supply, which was deposited by fluvial transport and accumulated in the former Mekong basin (Nguyen et al., 2000). Within the Miocene to Holocene deposits Table 1 Overview of aquifers in CMP. Depth and thickness of layers are derived from 27 borehole logs (Fig. 2) specifically for CMP. Average thickness calculated by interpolating down to the Miocene aquifer n 1 3 (chapter 2.2.4). Each aquifer is overlain by the corresponding aquitard (not included here). Ages after Cohen et al. (2013). The lower boundary from n 1 3 is estimated after Wagner et al. (2012).
overlying the Mesozoic basement, a complex sequence of sedimentary layers of marine, alluvial and estuarial origin developed. Seven confined aquifers and aquitards were identified (Table 1), which differ in appearance and extend over the whole delta (Anderson, 1978;Nguyen et al., 2000;Wagner et al., 2012). The aquifer system in detail is described in Wagner et al. (2012). Briefly, aquifers consist of alluvial, marine-alluvial, and marine fine to coarse sands as well as gravel, aquitards of marine fine sands, silt and clay (Doan et al., 2018;Pechstein et al., 2018;Wagner et al., 2012). Saline and fresh groundwater is heterogeneously distributed across the aquifers (Anderson, 1978;Pham et al., 2019;Wagner et al., 2012). Intensive Quaternary marine transgression and regression phases had and have a significant influence on delta development worldwide as they define the depositional environment as well as intensive periods of erosion (Custodio and Bruggerman, 1987;Pham et al., 2019;Schimanski and Stattegger, 2005). During the last glacial maximum (LGM) 26-18 ka BP sea level was as low as − 130 m below present MSL, which correspond to around − 68 m below the paleo surface of the MD at that time, leading to erosion of sedimentary deposits of several tens of meters (Pham et al., 2019;Ta et al., 2002). Based on a fresh-saline groundwater model, Pham et al. (2019) assume, that freshwater primary infiltrated during this phase. While sedimentary top layers with a high permeability allowed extensive freshwater infiltration, rising sea levels after the LGM, which reached − 60 m at the end of the Pleistocene (12 ka BP) and 2-4 m above present MSL in the following Holocene (12-6.5-4.5 ka BP), promoted rapid sedimentation forming large-scale clayey mud facies sealing of underlying layers from saltwater intrusion (Pham et al., 2019). At that time, the entire MD up to the Cambodian border was inundated by the sea. Since then, while sea levels declined to present MSL, until recent sea level rise, deltaic sedimentation continued, and the MD propagated rapidly over 250 km southeastward of the Cambodian border to the East Sea (Nguyen et al., 2000;Ta et al., 2002). The dominant propagation characteristics changed from tidal to tidal-wave influenced, leading to a more southward-directed sediment transport and an open-bay system, which eventually created the Ca Mau Peninsula (Liu et al., 2017;Nguyen et al., 2010;Ta et al., 2002). It is estimated that the Ca Mau Peninsula was inundated until 2.0-1.0 ka BP (Liu et al., 2017).
Many households obtain their freshwater, mostly for daily needs and small businesses, from unlicensed private groundwater wells tapping the Pleistocene aquifers between 70 and 200 m. Over half a million extraction wells are estimated in the MD (Danh and Khai, 2015). Shallow groundwater is often brackish and potentially polluted (Anderson, 1978;Pechstein et al., 2018). Groundwater is furthermore utilized for commercial drinking water production, smaller businesses and aquaculture related industry. In the last few years, much effort was made to implement a tap water infrastructure, where groundwater from Pliocene wells deeper 170 m is subtracted, treated and further distributed. So far, however, tap water is only available in more urbanized areas.
Due to this development, hydraulic pressure heads in the main production aquifers qp 2-3 , qp 1 and n 2 2 (Table 1) are now steadily declining, especially in populated areas like CMC, where hydraulic pressure heads dropped to − 25 m (Duy et al., 2021;Erban et al., 2014;Pechstein et al., 2018). Subsequently, a large drawdown cone around CMC developed (Erban et al., 2014), which disturbs the natural groundwater flow. A daily or seasonal fluctuation of groundwater levels is attributed to alternating extraction rates determined by water demand, rather than changes in natural recharge (Pechstein et al., 2018). Natural groundwater recharge from outside the MD is highly limited, if there is any, as the low hydraulic gradient as well as artesian conditions before 1990 reduced groundwater velocity significantly, leading to high residence times in the order of several thousand or ten thousand years (Anderson, 1978;Benner et al., 2008;Pham et al., 1994). Therefore, groundwater might have even been stagnant before large-scale groundwater extraction began in recent decades or generally in case of pristine deeper aquifers. Since the early study of Anderson (1978), the lack of reliable well data to deduce the recharge area and hydraulic conditions remains. Only twelve observation wells exist in the Ca Mau Peninsula, whereas seven were constructed after 2017. Anderson (1978) postulated that groundwater might originate in Cambodia where surface water and shallow groundwater is fresh. Based on 14 C and Deuterium/ 18 O analysis, Ho et al. (1991) determined the recharge area in the northeast extension of the delta and formation outcrops in Cambodia at higher altitudes. Furthermore, the flow directions were determined by 14 C isochrones. They consist of a northwest to southeast and a northeast to southwest component with transit times of 5.6 m/a for Pleistocene and around 2.5 m/a for Pliocene aquifers (Ho et al., 1991). Hydraulic conductivities in the aquifers are in the range of 1.0 * 10 − 3 m/s and 3.0 * 10 − 5 m/s and porosity is around 33% (Pechstein et al., 2018;Shrestha et al., 2016). Due to an impermeable Quaternary top layer a direct connection or recharge by surface water or precipitation in the study area is currently expected to be negligible (Anderson, 1978;Pham et al., 2019). Previously, it was assumed that the aquifers in the MD were not directly hydraulically connected, but Pham et al. (2019) assume that aquitards are heterogeneous and disconnected. This assumption is supported by the observation of similar geochemical and isotopic patterns as well as hydraulic pressure heads throughout the aquifers at some locations in the MD (Doan et al., 2018;Hanh and Roland, 2017;Pechstein et al., 2018). In addition, the moderate hydraulic conductivities of the aquitards do not necessarily preclude the interconnection of aquifers, especially since some aquitards comprise of considerable amounts of fine sand (Pechstein et al., 2018;Wagner et al., 2012). Furthermore, clear identification of aquifer and aquitard delineations in the area is difficult and not always possible. Drill logs conducted in the whole MD reveal that discontinues clay lenses are heterogeneously embedded within aquifers (Erban et al., 2013). Age determination using 14 C stable isotope techniques revealed, that the deep groundwater around U Minh (see Fig. 1) is older than 25 ka (Pechstein et al., 2018) and between 15 ka and > 40 ka in other parts of the delta, tending to increase from Pleistocene to Miocene (Doan et al., 2018;Ho et al., 1991;Hoang and Bäumle, 2019;Nguyen et al., 2004). Other studies showed inconclusive results with an apparent 14 C age ranging from recent to older than 40 ka BP. However, data validity is sometimes not given because sampling-and measurement techniques are not available or insufficiently documented (Pechstein et al., 2018). Groundwater in the aquifer qp 2-3 and deeper is assumed to be at least semi-fossil (Pechstein et al., 2018). The present geochemical composition of groundwater is considered to be influenced by factors such as Tertiary geomorphology, Quaternary transgression and regression phases during the past 60 ka BP, and water-rock interactions, rather than by recent recharge or seawater intrusion (Doan et al., 2018;Pham et al., 1994Pham et al., , 2019. A baseline study was conducted in 2017 summarizing some hydrological information specifically for CMP but also pointing out that much more research is required in the region (Jenn et al., 2017).

Geochemical investigation 2.2.1. Sample locations and salinity analysis
The distribution of fresh (TDS< 1 g/l) and saline groundwater (Rusydi, 2018) in the MD is highly heterogeneous and seems not to be correlated to depth or distance to the coast (Anderson, 1978;Pham et al., 2019). Saline groundwater prevails in the upper Holocene and late Pleistocene aquifer while fresh groundwater is unevenly distributed in Middle and Early Pleistocene, Pliocene, and Miocene aquifers (Wagner et al., 2012). In order to determine the salinity levels and distribution in CMP, a total of 190 groundwater samples were taken in five field campaigns between February 2017 and December 2019. As recent groundwater recharge in deep aquifers is negligible due to several less permeable aquitards and the low hydraulic gradient (Anderson, 1978;Pham et al., 2019) it was estimated that the time or season of sampling does not influence groundwater properties. The focus was on private household wells, but also samples from school wells, drinking water production companies, administrative buildings and smaller businesses were taken. Sample locations were randomly selected with the intention to cover the whole province. In addition, to cover the main production aquifers in CMP, wells of varying depth between 60 m and 310 m were chosen to represent the Middle Pleistocene qp 2-3 , Lower Pleistocene qp 1 , Middle Pliocene n 2 2 and Lower Pliocene n 2 1 aquifer (see Table 1). Since all wells are constructed in a closed system including a tap and a pump, construction details like well depth, diameter and year of construction rely on information given by the owner. Almost all private household wells are made of PVC with a diameter of 49 mm and are screened a few meters in the target aquifer. Prior to sampling, standing water in the pipe was exchanged and the equipment was cleaned thoroughly with fresh groundwater at each location. Physiochemical parameters including temperature, pH and electrical conductivity (EC) were determined on site using a multi-parameter portable meter (WTW Multi 3630 IDS).

Groundwater composition
Groundwater in coastal aquifers is often the result of mixing between seawater, both recent and paleo, and freshwater (Giménez-Forcada, 2010; Jones et al., 1999). The deviation from the conservative mixing line between the two endmembers seawater and freshwater due to cation exchange reactions between groundwater and aquifer sediments is one of the most remarkable phenomena in coastal groundwater geochemistry (Jones et al., 1999;Russak and Sivan, 2010). In order to gain deeper insights on the geochemical composition and to derive these processes, a comprehensive analysis for cations and anions was conducted for each sample. For both analysis, 25 ml of sample groundwater was separately filtered through a 0.45 µm cellulose-acetate filter (Satorius Stedim Biotech GmbH). 50 µL of high-purity nitric acid was added to prevent cation precipitation (APHA, 2017) and microbiological processes were inhibited for accurate anion analysis by adding 50 µL sodium azide (Vanderford et al., 2011). The groundwater samples were transported to Germany and analyzed with ICP-MS (Thermo Fisher X-Series 2) for cations and Ion Chromatography (Dionex ICS-1000 with column IonPac As14 Supressor ERS 500 and Methrom Compact IC 930 with column Metrosep A Supp 5-150. Eluent: NaHCO 3 /Na(CO 3 ) 2 ) for anions in the laboratory of the Institute of Applied Geosciences at the Karlsruhe Institute of Technology, Germany.
Additionally, total alkalinity was determined by measuring the acid capacity to pH 8.2 and pH 4.3 with a titration test kit (Merck KGaA). Concentration of HCO 3 and CO 3 2were directly derived from the alkalinity results. However, especially in redox-sensitive marine and coastal settings, the difference between total alkalinity and carbonate alkalinity may not be neglected due to the impact of other buffering bases such as B(OH) 4 -, NH 3 , Si(OH) 3 O -, or F -(Lukawska-Matuszewska, 2016).
The conversion factor for EC to TDS is highly depended on ionic water composition and may differ from the instrument value, which is typically calibrated by KCl solution (Vengosh, 2013). Therefore, TDS was calculated for each sample using the software Geochemist`s Workbench (v. 12.0, Aqueous Solutions LLC) with the whole elemental composition as input data. Furthermore, the direct relationship between all parameters was calculated using pairwise Spearman correlation. The calcium carbonate system in coastal groundwater settings is complex, and many processes actually lead to carbonate dissolution with subsequent increase in HCO 3 and Ca 2+ (Custodio and Bruggerman, 1987). In saline aquifers, mixing of calcium carbonate saturated fresh groundwater with seawater, which is typically oversaturated, produces brackish water which is undersaturated with regard to calcium carbonate due to the reduction in ionic strength, leading to calcite dissolution (Custodio and Bruggerman, 1987 Diagram (Piper, 1944). The water type can be derived from the upper diamond shape (Back and Hanshaw, 1965). Water type analysis can help to distinguish geochemical processes and to derive the origin of solutes. However, divisions and percentages used for classification should be selected as to best fit the characteristics of the investigated water (Freeze and Cherry, 1979). Therefore, to achieve a more detailed water type specification and to get more insights on mixed water types, the typical description was customized in order to more accurately display different states of salinization with an emphasis on cation-exchange processes. Proportions of major ions were calculated and if only one cation or anion was above 33.3%, it is written alone, like Na-Cl water type, while if two ions were above 33.3%, the ion with the higher proportion was written first, followed by the second ion in brackets, for example Na-Cl-(HCO 3 ). A mixed water type is thus only defined if all three anions or cations have an equal share of around 33.3%, e.g. Na-Ca-Mg-Cl is written as Mix-Cl.

Salinization and freshening
The evaluation of molar ion ratios is an appropriate tool to distinguish between seawater and freshwater intrusion in coastal aquifers (Appelo, 1994;Jones et al., 1999;Russak et al., 2016). Therefore, Na + /Cl -, Ca 2+ /Cland Mg 2+ /Clmolar ratios were calculated and compared to the conservative mixing line. Generally, salinization of a freshwater aquifer by seawater intrusion proceeds in the following steps (Appelo and Postma, 2004;Custodio and Bruggerman, 1987;Jones et al., 1999;Vengosh, 2013): In the early stages of salinization, fresh Ca 2+ -dominated groundwater is enriched in Na + and Mg 2+ . In general, divalent cations are exchanged on the mineral surface preferentially to monovalent cations (Carroll, 1959;Steudel and Emmerich, 2013). However, at a certain intensity of salinization ion selectivity is masked by the proportional increase in sorption with increasing concentration, especially of the monovalent ion Na + (Sawhney, 1972). In conclusion, during further salinization Na + exchanges divalent cations Ca 2+ and Mg 2+ from sedimentary rocks forming Ca-Cl and Na-Mg-Cl waters. Simultaneously, overall salinity increases, while also the proportion of divalent cations increases.
Conversely, freshwater inflow into saline coastal aquifers, also referred to as freshening, reverses the cation-exchange reaction contrary to the salinization process (Appelo, 1994). When freshwater of Ca-HCO 3 type intrudes the saline aquifer exchangeable Na + on mineral surfaces, especially clay minerals, is substituted by Ca 2+ and Mg 2+ , due to a higher sorption selectivity, and the water type changes to the characteristic Na-HCO 3 . The Na + /Clratio increases above the level of seawater, while Ca 2+ as well as Mg 2+ become depleted (Appelo and Postma, 2004;Custodio and Bruggerman, 1987;Jones et al., 1999).
In the natural environment, Cloften acts as a conservative tracer and ion ratios are commonly calculated relative to Cl - (Jones et al., 1999). Sea water has a Na + /Clmolar ratio of around 0.86. Due to the above-mentioned cation exchange reactions, Na + /Cldecreases below the conservative mixing line during salinization and increases during freshening.

Hydrogeological model
To fully evaluate the fresh-saline groundwater patterns in the context of stratigraphy and depth, samples were assigned to their respective aquifer based on available lithological information. As owners of private wells do not know their target aquifer and because no sufficient subsurface model was available for the region, a simplified hydrogeological model was constructed using lithological data from 27 borehole logs in CMP (Vuong, 2014), provided by the National Center for Water Resources Planning and Investigation Vietnam (NAWAPI), whereas six of them are part of the national monitoring network. The locations of the boreholes are shown in Fig. 2. The model serves two purposes: (i) assignment of samples to an aquifer and hence (ii) visualization of samples in depth to delineate aquifer dependent groundwater properties. A layer-based approach was chosen, where surfaces were created separately and Fig. 2. Location of 27 boreholes (Vuong, 2014) and extent of the interpolation area (purple) used for the calculation of the hydrogeological model. then stacked (Turner, 2006). This approach is common in sedimentary environments and is widely used for groundwater mapping (Natali et al., 2012). Subsurface modelling was conducted using MATLAB 2020b (MathWorks Inc.). Top and bottom boundaries from 14 layers recorded in 27 boreholes were used as 3-D scattered input data.
While almost all boreholes reach aquifer qp 1 , except one, the number of boreholes reaching deeper layers are gradually decreasing with only seven reaching aquifer n 2 1 and only one the aquitard beneath. As this is a typical phenomenon in lithological data acquisition (Turner, 2006), all boreholes were vertically interpolated down to the Miocene aquifer n 1 3 based on mean layer thickness from the available data. No information from outside the study area was used because depth and thickness of stratigraphic units significantly vary throughout the delta (Nguyen et al., 2000) and might have therefore distorted the calculation. For the model, a 2-D grid was created with a 100 m resolution in latitude and longitude direction using the meshgrid command, and for each layer a 3-D data interpolation was done using triangulation-based natural neighbor method with griddata, which was then visualized using the mesh command. It should be noted that the quality of the provided borehole logs could not be further evaluated here and thus no extensive model validation was conducted. Based on the straightforward model approach, the model does not consider discontinues clay lenses within aquifers (Erban et al., 2013) and other potential small-scale characteristics, such as paleo channels as well as areas of non-deposition and areas of erosion (Turner, 2006). Groundwater samples were plotted in this subsurface model to enable a sample to aquifer assignment. In order to account for a thin aquitard between aquifer qp 2-3 and qp 1 and to account for possible aforementioned uncertainties in the interpolation, it was only distinguished between Pleistocene and Pliocene samples. Samples outside of the convex hull were correlated with the nearest borehole and assigned to the respective aquifer. In addition, a 71.5 km long geological cross section was drawn through the model area. Samples, larger cities and monitoring wells within a 10 km perpendicular distance to the profile were included. Location of monitoring wells and the course of the cross section are displayed in Fig. 2. Most monitoring wells accumulate around CMC. Extrapolation outside the convex hull of the available borehole data was avoided intentionally to assure a more valid model within the verifiable boundaries. Therefore, the cross-section was not extrapolated to the sea.

Water geochemistry
Major geochemical water parameters are listed in Table 2. EC values range between 571 µS/cm and 5180 µS/cm. 21 samples have EC between 2500 and 3850 µS/cm while two samples in the northwest have EC around 5000 µS/cm (Fig. 5c). A mean conversion factor for EC to TDS of 0.81 was calculated. Note that this factor deviates from 0.67, which is commonly used by NAWAPI. As this conversion factor is highly depended on ionic water composition the conversion factor presented here is more accurate for the study area and should be considered for future water use guidelines. Na + is the dominant cation in the majority of samples in CMP. The Piper Diagram in Fig. 4 shows that Na-Cl, Na-HCO 3 and mixed water types occur in Pleistocene and Pliocene aquifers. Pliocene samples are depleted in Ca 2+ + Mg 2+ and contain over 70% Na + , while some Pleistocene samples contain more divalent cations, especially those with high EC and a higher proportion of Cl -. Sulfate only plays a minor role in the water composition, is mostly below 10% and never exceeds 25%. The spatial distribution of the detailed water type (see chapter 2.2.2) is presented in Fig. 3. In this classification the dominant water type is Na-HCO 3 and Na-HCO 3 -(Cl). Only eleven samples are solely dominated by Cl -. It is striking that these samples of Na-Cl type, shown in red, are randomly distributed in CMP. Two samples in the north of CMP are of Ca-HCO 3 and Ca-(Na)-HCO 3 type, colored in blue.
The spatial distribution of EC in CMP (Fig. 5c) shows that EC is generally lower in the northeast and east. In the west and south, EC is mostly between 1000 µS/cm and 2000 µS/cm with some scattered hotspots of high EC above 2500 µS/cm. The concentration of HCO 3 shows a rather sharp spatial pattern and is more than tripling from north to south from around 250 mg/L to 950 mg/L (Fig. 5a).
By contrast, the concentration of Clis more evenly distributed ranging mostly between 25.0 and 250 mg/L (n = 121) and 250-750 mg/L (n = 40). Only seven samples exceed 750 mg/L and 22 samples contain less than 25 mg/L Cl -, mainly accumulating northeast of CMC (Fig. 5b). Ca 2+ and Mg 2+ show a high correlation (p = 0.94) ( Table A1 in appendix A) and high concentrations are linked to samples with highest Clcontents ( Fig. A2 and A3 in appendix A). Furthermore, except a few samples, SO 4 2concentrations are relatively low (Table 2) and are more evenly distributed (Fig. A4 in appendix A). Temperature ranges between 25.7 • C and 37.2 • C and generally increases with depth, following a normal geothermal gradient (Morgan, 1984). No significant correlation between depth

Hydrogeological model
The 71.5 km long hydrogeological cross section shown in Fig. 6 includes samples within a 10 km distance to the profile, which are colored according to their EC levels (compare Fig. 5c). Aquifers are labeled according to Table 1 and are illustrated in blue colors. Based on the 27 available borehole logs in the project region, the model approach and interpolation method could adequately construct the 14 distinct layers that overly each other smoothly without significant irregularities. Therefore, layered sequences seem to be dominant in the study area.
While Holocene layers are flat, deeper layers appear to dip down southwestward. 167 samples were assigned to the Pleistocene aquifers qp 2-3 and qp 1 and 22 samples to Pliocene aquifers n 2 2 and n 2 1 . One sample might originate from the Miocene aquifer n 1 3 from a depth of 310 m, but is added to the Pliocene samples for further interpretation for simplification. The profile illustrates that EC is generally below 1000 µS/cm north of CMC and increases southwestward to 1000-2000 µS/cm. Along the cross section, groundwater samples with EC > 2000 µS/cm are found solely and randomly distributed southwest of CMC, primarily in qp 2-3 but also in deeper aquifers towards Nam Can.

Groundwater composition
The shift of the fresh-saline water boundary further inland is a major concern for water management in CMP, as an increase of salinity would significantly degrade groundwater resources and limit their use. 40 out of 190 samples exceed the Vietnamese limit for total dissolved solids (TDS) of 1500 mg/L (MONRE, 2015), which is equivalent to an EC value of 1850 µS/cm (see chapter 3.1). 47 samples exceed the threshold for Clof 250 mg/L. Thus, some groundwater is already unsuitable for direct drinking purpose. Data collected in this study, however, also show that locations with elevated EC are randomly distributed across CMP and are most often represented by more than one sample at each location from both Pleistocene and Pliocene aquifers. This leads to the conclusion that locations of elevated EC are not a single well phenomenon but a local spatial pattern. This is also true for other ions and ion ratios, as shown in Figs. 5 and 7. This indicates that aquifers are to some extend connected, at least partly, either through natural discontinuities in the stratigraphy, as also recently proposed by other authors (Doan et al., 2018;Pechstein et al., 2018;Pham et al., 2019), or hydrological short circuits caused by well leakage or well failure (Doan et al., 2018), which are rarely considered. In a regional perspective, no significant correlation between geochemical patterns and depth or aquifer was found, as shown by the correlation analysis (Table A1 in appendix A) and by the geological cross section (Fig. 6). In addition, locations with higher EC are not correlated to the distance to the coast.
As a simple and straightforward measure of salinity, EC is a useful parameter to characterize and classify groundwater resources in coastal regions potentially affected by seawater intrusion, especially when detailed analysis is not available. Seawater of Na-Cl type is largely uniform with a TDS of around 35 g/L, while freshwater is often defined with a TDS concentration below 1 g/L (Vengosh, 2013), meaning that even slight seawater intrusion is detectable. However, groundwater composition is often more complex and other parameters than Na + and Clare affecting EC as well (Jones et al., 1999). This also applies to data presented in this study. As shown in Fig. 5a-c, elevated levels of EC in the southern part of the study region cannot solely be attributed to Classociated with seawater but rather to the presence of high HCO 3 concentrations. In fact, around Nam Can groundwater with EC > 2000 µS/cm is of Na-HCO 3 and Na-HCO 3 -(Cl) type (Fig. 3) and HCO 3 concentrations are lowest inland in the northeastern area, where recent freshwater inflow is estimated to occur, even though it might be limited (Fig. 5a). This indicates that the source of HCO 3 is not freshwater intrusion. Almost all samples are oversaturated to calcite (n = 158), especially in the south and central south between CMC and Nam Can (Fig. A5 in appendix A), which should lead to calcite precipitation. Only a few samples in the northern part of the study area are slightly undersaturated to calcite (n = 5) or in equilibrium (n = 27). One reason for the oversaturation of calcite might be the addition of Ca 2+ due to cation exchange processes or the dissolution of gypsum (Plummer et al., 1990), which has a SI between − 1.56 and − 6.21. This contradicts typical hydrochemical facies evolution in coastal areas, where SI calcite tends to decrease towards the sea (Custodio and Bruggerman, 1987). Ca 2+ and HCO 3 have a negative correlation and while HCO 3 concentrations increase, Ca 2+ concentrations are relatively low, probably because Ca 2+ is involved in complex ion exchange processes (Jones et al., 1999). It is therefore concluded that observed high concentrations of HCO 3 in the south of the study area cannot be explained by simple mixing and calcite dissolution.
Another possible source of HCO 3 in coastal aquifers is dissolved organic carbon derived from organic matter, which is oxidized to dissolved inorganic carbon. It was shown in other studies that the concentration of HCO 3 originating from organic matter oxidation can surpass the concentration in fresh groundwater at recharge areas multiple times (Chapelle and McMahon, 1991;Orem and Gaudette, 1984). In addition, organic matter reduction also shifts the calcium carbonate equilibria due to an increase in P CO2 as well as changes in pH, and causes dissolution (Custodio and Bruggerman, 1987;Jones et al., 1999). Due to their marine and estuarial origin, high amounts of organic matter are consistently found across all stratigraphic layers in CMP (Doan et al., 2018;Nguyen et al., 2010;Pechstein et al., 2018) and could therefore account for high HCO 3 concentrations. Furthermore, sulfate reduction directly produces HCO 3 and leads to sulfate depletion, which is typical for freshening aquifers rich in organic matter and might therefore also apply to Fig. 8. Na + /Clplot of groundwater samples including water type markers associated with salinization of freshwater aquifers.
CMP where as a consequence sulfate concentrations are rather low (Kim et al., 2017;Chapelle and McMahon, 1991). Overall, this is an important finding, demonstrating that the increase of EC values in coastal areas are not necessarily an evidence of recent seawater intrusion, but can also be attributed to the generation of additional HCO 3 -.

Salinization and freshening processes
Coastal aquifers are influenced by marine transgression and regression cycles leading to seawater or freshwater intrusion, respectively (Custodio and Bruggerman, 1987;Han et al., 2020). Fig. 8 reveals that the Na + /Clmolar ratio of most samples in the study area is well above the seawater ratio of 0.86. Only a few samples show ratios lower than seawater. It is remarkable that the water type of those samples below the seawater ratio line is changing with increasing Clconcentration from Mix-Cl-(HCO 3 ) and Ca-Cl-(HCO 3 ) over Na-Mg-Cl and finally to Na-Cl, while almost all samples above the seawater ratio are of Na-HCO 3 and Na-HCO 3 -(Cl) type. This observation in surveyed data convincingly demonstrates the impact of cation exchange processes during salinization and freshening. The customized water type is well suitable to reflect different stages and the geochemical processes associated with salinization. It can be observed that Na-Mg-Cl seem to occur at higher Clconcentrations compared to Ca-Cl-(HCO 3 ) (see Fig. 8), probably because seawater contains more Mg 2+ over Ca 2+ with a Mg 2+ /Ca 2+ ratio of 4.5-5.2 (Vengosh, 2013). The Mg 2+ /Clratio is shown in Fig. 9 and highest ratios correlate with those of lowest Na + /Clratios and salinization water types, again demonstrating the cation-exchange processes. The same pattern can be observed for Ca 2+ /Cl - (Fig. 10) but with a wider spread and deviation from the mixing line, as Ca 2+ concentrations are more influenced by other factors such as calcite precipitation and dissolution, which in turn depend on HCO 3 concentrations.
In conclusion, freshwater intrusion is, or was, the predominant process leading to the abundant Na-HCO 3 and Na-HCO 3 -(Cl) type groundwater in CMP. The rarely occurring groundwater of Na-Cl type above the seawater line might be in a less advanced stage of freshening, where Clis still the dominant anion. Further freshening produces waters of Na-HCO 3 type and Clbelow 25 mg/L, as found in the northeast of the study area. Due to the general hydrogeological conditions, it can also be assumed that the freshening process has progressed furthest in this region. This assumption is supported by Na + /Clmolar ratios greater than 10.0 in the north-east and Na + /Clbetween 5.00 and 10.0 in the central east (Fig. 7).
It can be expected that both processes, cation exchange reactions and the release of organic carbon by organic matter decomposition, contribute to observed Na-HCO 3 groundwater in CMP. Detailed water type analysis in combination with the evaluation of Na + / Clmolar ratios and the consideration of cation exchange processes proved to be an appropriate tool to identify signs of salinization, even in early stages at low Clconcentrations. Focusing on total salinity via EC analysis alone, however, is not a suitable method. For example, three samples in the central north (Fig. 5c) have EC values below 1000 µS/cm (Fig. 5c). However, based on their Na + /Clratio, which is below the seawater line (Fig. 7), and their water type (Fig. 3), particular these samples seem to be in an early stage of salinization.

Groundwater evolution and development
Marine transgression and regression phases in the Quaternary (Fig. A1 in appendix A) caused numerous cycles of alternating saltwater and freshwater intrusions in the MD and coastal regions worldwide (Delsman et al., 2014;Pham et al., 2019;Spratt and Lisiecki, 2016;van Engelen et al., 2021). In this period, phases of sedimentation during marine transgression were alternating with phases of extensive erosion during marine regression. Derived and edited after Pham et al. (2019), Fig. 11a-c illustrate three characteristic stages of past hydrogeological situations in the lower MD. The formation of Na-HCO 3 type groundwater dominant in CMP at moderate EC levels with high Na + /Clratios supports the assumption of Pham et al. (2019) that extensive freshwater intrusion occurred in the aquifer system of the MD. Previously, this assumption was only covered by the distribution of TDS. It furthermore Fig. 9. Mg 2+ /Clplot of groundwater samples including water type markers associated with salinization of freshwater aquifers. proofs that the system was saline before freshwater intrusion occurred through permeable top layers. These findings are evident in both shallow and deep aquifers concurrently. The freshening event presumably happened in the last marine regression phase 60-12 ka BP, especially during the LGM when sea levels were extremely low (Fig. 11a), which would explain groundwater ages in that range (Pham et al., 2019).
The formation of large-scale impermeable Holocene aquitards, particularly QP3, between 8 and 6 ka BP prevented extensive seawater intrusion, which otherwise would have been expected due to rapidly rising sea levels in the Holocene transgression phase (Fig. 11b) and would have also led to the respective salinization signatures, namely Na + /Cl -< 0.86 and associated water types. The phenomenon of Pleistocene freshwater conservation can also be observed in other coastal areas worldwide (Cohen et al., 2010;Delsman et al., 2014;Han et al., 2020;Post and Kooi, 2003;van Engelen et al., 2021). However, such salinization signatures are present at a few locations. Due to the high groundwater exploitation rates as well as the high hydraulic connectivity within and apparently across the aquifers it is rather unlikely that these samples reflect remnants of older saline groundwater unaffected by the latest freshening event. Furthermore, these samples plot below the conservative mixing line at low to moderate EC levels, which indicates that they developed from freshwater and are in an early stage of salinization (Fig. 8) rather than in a less advanced stage of freshening, as it is the case for same samples of Na-Cl type above the mixing line. Therefore, it is more likely that former freshwater mixed with saline water more recently. However, the origin of this recent salinization is not further examined in this study, but it should be noted that the MD was completely inundated at around 5.5 ka BP and CMP until estimated 2 ka BP (Fig. 11c) (Nguyen et al., 2000;Pham et al., 2019). The current stratification with high-density Holocene saline aquifers on top of low-density deep freshwater aquifers reflects an unstable aquifer situation, which might lead to a downward migration of the younger saline groundwater along natural and anthropogenic pathways (Custodio and Bruggerman, 1987;Post and Kooi, 2003). It must be noted that the subsurface in the study region is virtually perforated with abandoned as well as active wells of different construction quality and reliability.
The Piper diagram (Fig. 4) shows that in average the proportion of Na + is notably higher in Pliocene than in Pleistocene samples, which is illustrated in Fig. 11c. This finding might indicate the heterogeneous downwards migration of saline water from top to bottom, which leads to the release of divalent cations in Pleistocene aquifers, which therefore show a wide range of Na + (%). In addition, differences in geomorphology and discontinuities in the stratigraphy probably resulted in non-uniform downward migration  of freshwater, which would explain highly heterogeneous geochemical patterns and partly inconsistent groundwater ages observed over the whole MD (Anderson, 1978;Pham et al., 2019;Wagner et al., 2012). Although information on Miocene groundwater in CMP is rare, one sample probably originating from the Miocene aquifer n 1 3 (Fig. 6). With EC of 3090, Na + over 86%, Na + /Clof 1.5 and Na-Cl water type, the sample is also in a less advanced stage of freshening.
In addition, this sample is consistent with samples from the n 1 3 aquifer in Soc Trang Province bordering the southern Mekong River mouth. They have similar geochemistry and are described as slightly brackish but generally more uniform. Pleistocene samples in Soc Trang do also show a wide range in Na (%) (Hoang and Bäumle, 2019;Tran et al., 2020), which demonstrates the applicability of the concept to other coastal provinces in the MD, even outside the Ca Mau Peninsula. Considering that direct vertical freshwater intrusion must have passed through several layers of various permeability, intensity of freshening is estimated decrease with depth (Doan et al., 2018;Pham et al., 2019). In addition to the depth-depended observations the spatial component in surveyed data can be further examined. Freshening signatures are more pronounced in the northern part inland of CMP with Na-HCO 3 water type, Cl -< 25 mg/L and high Na + /Clratios. This cannot be explained by seaward-directed, naturally occurring lateral freshwater recharge northeast to southwest parallel to the delta propagation axis, which is controlled by the minor marine regression in the last 5.5 ka, until recently. Due to the confined aquifer situation and the extremely low hydraulic gradient in the flat low-lying MD, the natural recharge is estimated to be negligible. Therefore, the almost completed freshening may be the result of accelerated freshwater inflow from outside the study area due to the significant drop of hydraulic heads in recent decades, mainly around CMC. Decreasing hydraulic heads lead to steeper hydraulic gradients and depression cones promoting more recharge at higher velocities from outside the pumping area. Furthermore, the random distribution of groundwater with high Clconcentration and salinization signatures highlights that recent large-scale seawater intrusion into deep aquifers is not occurring and these samples are rather the result of the aforementioned local downwards migration of shallow saline water. There is no correlation between these samples and the distance to the coast or any other consistent pattern, which are regularly observed in other coastal areas and are mostly related to one or more discrete saline sources, such as the upwelling of old saline brines or anthropogenic sources (Marie and Vengosh, 2001;Vinson et al., 2011). A clear fresh-saline boundary was not identified. As pointed out be Delsman et al. (2014), coastal groundwater is only rarely in equilibrium with present boundary conditions and some processes are slow. However, aquifer system salinization can be a fast process under some circumstances (Post and Kooi, 2003;Russak et al., 2016). The direct intrusion of seawater due to pumping would result in accelerated groundwater salinization (Vengosh, 2013) while the downwards migration of shallow saline groundwater might be much slower. Observed geochemical patterns might be more related to a disturbed hydrological system due to pumping rather than direct lateral seawater intrusion. In a short term, appropriate counter measures should be considered at locations where a significant increase in salinity is evident. In a longer term, however, it should be expected that salinization of deep aquifers will arise or continue and that the quantity of freshwater will decline if current practices are not reconsidered.

Conclusions
Groundwater data and detailed evaluation of geochemical signatures, including water types, salinity and element distributions, as well as element ratios, allowed for the development of a comprehensive model for the evolution of groundwater geochemistry in CMP. Based on the predominance of Na-HCO 3 water types, salinity patterns, and the distribution of Na + /Cl -, the groundwater evolution in the study area can be explained by the downward dilution of infiltrating freshwater into an initial saline system. The freshening event happened during the last marine regression phase 60-12 ka BP and especially during the LGM, 26 ka BP to 18 ka BP, and shows slightly decreasing intensity from Pleistocene to Miocene aquifers. Element ratios of Na + /Cl -, Ca 2+ /Cl -, as well as Mg 2+ /Cland their deviation from the seawater line were used to determine cation exchange processes between groundwater and aquifer sediments, and thus to characterize salinization and freshening processes. As demonstrated, evaluation of characteristic element rations is a powerful tool to determine groundwater state and evolution in coastal environments. This is particularly valid in contrast to the sole analysis of salinity via EC measurements, which has not proven to be a suitable indicator of seawater intrusion.
Moreover, hydrogeological layers in CMP are estimated to be connected, at least partly, either through natural discontinuities or anthropogenic influences. Even though the formation of impermeable aquitards in the Holocene marine transgression phase prevented extensive seawater intrusion, some samples show distinct salinization signatures of unknown origin. Salinization due to seawater intrusion would result in accelerated groundwater degradation. Therefore, it is important to distinguish direct seawater intrusion from mixing with other saline water sources, like the slow downwards migration due to unstable aquifer conditions, to establish a sustainable water management plan necessary to provide clean freshwater to the people. These findings have a direct implication to the whole MD region.
More work about the current state of salinization and freshening in the MD will be conducted using geochemical fingerprinting methods like ion ratios, trace elements, stable isotope techniques and geochemical models. By considering the new insights on natural drivers leading to the current groundwater composition in CMP presented in this study, the effect of pumping induced changes in hydraulic parameters and geochemistry can be appropriately evaluated in a future study, also considering pumping induced land subsidence. Lastly, concepts for alternative water resources substituting groundwater should be considered and developed to mitigate negative environmental impacts.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.