**1. Introduction**

The prospective for water resources management in the Vietnamese Mekong Delta (VMD) is very challenging in the 21st century due to multiple issues including increasing water use, climate change impacts and sea level rise [1,2]. These changes will likely have considerable impacts on economic development and the livelihoods of people living in the delta [3]. In this context, the Mekong Delta

is considered to be one of the most vulnerable catchment areas throughout the basin because of the cumulative impacts of upstream development activities. The Mekong Delta has long been identified as a hotspot in terms of vulnerability to sea level rise and climate change, as well as to hydropower dam construction in the upper basin [4,5]. As a result of these factors, the Mekong Delta experiences increasing salinity intrusion, which represents a critical challenge for water resources management and agriculture production.

The VMD plays a crucial role in the economy of Vietnam. The delta contributes approximately 27% of the country's GDP, supports 16 million inhabitants (nearly 22% of Vietnam's total population), and provides about 50% of the annual rice production [6,7]. Vietnam is also one of the most vulnerable countries to climate change and sea-level rise [8]. Currently, the VMD is experiencing rapid development in terms of population growth and infrastructure development [9,10], which create substantial pressure on water resources [11,12]. Moreover, hydropower development and water withdrawal in upstream countries cause changes in the hydraulic regime and salinization in the VMD [5,13]. All these changes pose critical challenges for sustaining water resources for economic development. In view of these challenges, it is important to assess future hydrological changes, especially salinization, in order to support decision making and planning for sustainable delta management.

Salinity intrusion is a natural phenomenon occurring in the lands, estuaries, and aquifers adjacent to the sea. Salinization in deltas and estuaries varies considerably depending on the tidal regimes, river flows and topography [14]. Tidal dynamics cause turbulent mixing via the transportation of saltwater, as well as tidal trapping and tidal transport via the interaction between tide and terrain [15]. The tide is a powerful source of mixed fresh and saltwater, that also plays an important role in saltwater intrusion. Recent studies also identified a wide range of factors affecting salinity intrusion, including river flows, topography, morphology, river bed slope, wind velocity and direction, and water temperature [16,17].

Salinization in the Mekong Delta under climate change and sea level rise has been addressed in a number of recent studies. Several studies have shown that a combination of climate change, sea level rise and upstream inflow changes will increase salinization, resulting in important consequences for water supplies and agriculture production [2,8,18–20]. For instance, Wassmann et al., Le et al., and Dinh et al. [21–24] investigated the effects of sea level rise on water levels in the VMD. Khang et al. [20] demonstrated the relationship between salinity intrusion and sea level rise and river flow change in the VMD during the dry season. Most studies covered the whole VMD using the 1D-MIKE 11 hydrodynamic model and, therefore, did not pay particular attention to salinity intrusion and hydraulic regime in the estuaries.

Distributed hydrological models are frequently applied in water resources assessments at the Mekong basin level, but little attention has been paid to the hydraulic regimes in the estuaries due to highly complex river–estuary–ocean interactions. The existing modelling frameworks are not well-designed to cope with the complex hydrodynamic relationship between river flow, tides and human activities. To address this challenge, we propose a combination of hydrodynamic models in order to capture the hydrodynamics of the Mekong Delta's estuary sufficiently. Therefore, the objective of this study is to combine 1D-MIKE 11 and 2D-MIKE 21 hydrodynamic models to investigate the hydraulic regimes and salinity intrusion in a specific river section, namely the Hau River, stretching between Can Tho to Tran De and Dinh An estuaries. To achieve this objective, a 1D model was first applied to simulate the entire hydrodynamic of the VMD with different scenarios of upstream discharge and downscaled precipitation from five general circulation models (GCMs) with two representative pathway scenarios (RCP) from Coupled Model Intercomparison Project Phase 5 (CMIP5). We employed statistical downscaling and bias-correction methods to prepare climate change input data, as shown in detail in Duong et al. [25]. The outputs of downscaled daily precipitation for 2036–2065 were used to run a rainfall-runoff model to simulate runoff for 1D hydrodynamic simulations. Finally, we used the results at Can Tho and branch discharge along the Hau River from MIKE 11 simulation as a boundary condition to force the 2D hydrodynamic model to simulate river–estuary–ocean interactions. As such, our modelling framework allows evaluation of the combined impacts of changing upstream discharges, the variabilities of precipitation and sea level rise from climate change on salinity intrusion in Hau River, VMD.

#### **2. Study Area**

The Vietnamese Mekong Delta (VMD) is located at the end of the Lower Mekong Basin. It is a large, relatively flat and low-lying area, with an elevation of approximately 0.5–1.0 m above mean sea level. The delta's hydraulic regimes are complicated, with two major distributaries—the Tien (Mekong River) and Hau River (Bassac River)—which drain out to the South China Sea through eight estuaries (Figure 1). The structures comprise 7000 km of main canals, 4000 km of secondary on-farm canal systems, 193 spills, 409 reservoirs, 528 junctions, 29 sluices, 749 floodplains and more than 20,000 km of dykes to prevent early floods [26,27]. The mean annual flow is 15,000 m3/s, the maximum discharge is about 60,000 m3/s in the flood season, and the minimum discharge is around 2000 m3/s in the dry season [28,29]. The climate of the VMD has tropical monsoon characteristics, with two separate seasons per year. The rainy season normally lasts from May to October, whereas the dry season lasts from December to March [30]. The mean annual precipitation is 1500 mm for the entire VMD and the range varies from 1600 mm to 2400 mm/year. The total precipitation in the rainy season contributes to nearly 90% of the annual precipitation [31]. The average discharge in the dry season fluctuates greatly, from 1700 m3/s to 6000 m3/s between January and May, and leads to water shortages for irrigating about 1.5 million hectares of irrigated crops [27]. The Mekong Delta experiences two types of tidal regime, namely the semi-diurnal and diurnal tides, affecting the hydraulic condition in estuaries of the East and West seas, respectively [32]. The salinity intrusion in the VMD—and in all the river networks—is substantial in recent years, particularly during the dry season, during which about 2.1 million hectares in the Mekong Delta suffered from salinity intrusion [21,27].

**Figure 1.** *Cont.*

**Figure 1.** The study area for the whole Vietnamese Mekong Delta (VMD) (top panel) and 2D study domain (bottom panel).

#### **3. Methodology and Model Setup**

#### *3.1. Hydrological and Hydraulic Data*

The hydrological input data for the MIKE 11 model includes river discharges at the delta inlet at Kratie (Cambodia), outlet water levels, precipitation, evapotranspiration and water demands for agriculture, industry and domestic sectors in the entire VMD. The hourly Kratie discharge and hourly water levels at 10 major stations, including Vung Tau, Vam Kenh, Binh Dai, An Thuan, Ben Trai, My Thanh, Ganh Hao, Song Doc, Rach Gia and Xeo Ro, were measured from 2009 to 2011 and all input data is summarized in Table 1. More details about these data can be found in Duong et al. [25].

The daily precipitation scenarios for boundary conditions in the MIKE 11 model were generated from the CMIP5 using five GCMs: ACCESS 1.0, CCSM4, CSIRO-Mk 3.6, HadGEM and MPI-ESM-LR, with two RCP4.5 and RCP8.5 scenarios (Table 2). The bilinear interpolation was applied to downscale climate data to 0.5◦ × 0.5◦ grid before applying bias correction. Future precipitation change scenarios were generated using three bias correction methods, including linear scaling [33], local intensity scaling [34] and distribution mapping [35,36]. The motivations for selecting the GCMs and RCP scenarios were explained in detail in [25,37].

The hydraulic data for the MIKE 21 model consists of hourly discharge at Can Tho station, offshore tidal levels and salinity concentration at Tran De and Dai Ngai stations, in the years 2011 and 2010 for calibrating and validating the models, respectively. The tidal levels in the years 2010 and 2011 were derived from the Global Tidal Model in the MIKE Zero Toolbox and were calibrated with tidal levels data at My Thanh station, provided by the Institute of Coastal and Offshore Engineering (http://www.icoe.org.vn/index.php?pid=551). Furthermore, the branch discharges along the Hau River were simulated from the MIKE 11 model, for five main intakes, including Mang Thit, Rach Mop, Cau Quan, Nga Bay and Dai Ngai. The future discharges at Can Tho were simulated from MIKE 11 with four scenarios of upstream discharges at Kratie and precipitation in the VMD.

Regarding the wind data for the MIKE 21 simulation, there are two dominant wind directions; the wind blows south-west in the wet season and northeast in the dry season, with stronger speeds as a result of the monsoon in the South China Sea. The average north-east wind speed varied between

3 and 6 m/s, and the maximum speed was approximately 15–20 m/s towards the east, blowing from the sea. The predominant wind directions were calculated and are shown in Table 3.


**Table 1.** Main input data for the model simulations.

SRHMC: Southern Regional Hydro-Meteorological Center, SIWRR: Southern Institute of Water Resources Research, CMIP5: Coupled Model Intercomparison Project Phase 5, DHI: Danish Hydraulic Institute, GCM: General Circulation Model, MRC: Mekong River Commission.



**Table 3.** Direction of wind in Vung Tau station, Southern Vietnam (SIWRR, 2011).


Note: E = east; ENE = east-north-east; NE = north-east; NW = north-west; SE = south-east; and SW = south-west.

Salinity data was used for calibrating and validating the MIKE 21 model for the years 2011 and 2010, respectively. The salinity concentration in the Hau River changed significantly during the dry season and reached its highest value between March and May (see Figure 2). The Dinh An and Tran De estuaries are characterized by an irregular semi-diurnal tide, which has a strong effect due to tidal oscillations [38,39]. Similar to salinity, the discharge at Can Tho shows considerable fluctuation within one day, and the diurnal differences are evident. The salinity and tidal levels showed phase differences, with the maximum salinity at the Dai Ngai station occurring during the spring tide. The maximum salinity in 2011 at the Tran De station was about 23.0 PSU (practical salinity unit; 1 PSU equals 1‰) and the minimum value was 0.5 PSU [40]. The salinity concentration at Dai Ngai station was smaller than Tran De because this station was located about 31 km further inland. In 2011, the maximum value of salinity level in Dai Ngai reached 11.0 PSU, and minimum values varied between 0.0 and 0.5 PSU.

**Figure 2.** The salinity concentration at Tran De (**a**) and Dai Ngai (**b**) stations in the dry season of 2011.

#### *3.2. Model Schematization*

The detailed bathymetry in numerical modelling plays a critical role in achieving accurate hydrodynamic simulations. In the present study, the bottom topography data in the coastal area and the Hau River estuary were obtained from the DHI Vietnam (Danish Hydraulic Institute), the Department of Transportation in Ho Chi Minh City, and the Southern Institute of Water Resources Research (SIWRR). The modelling domain in the horizontal plane covers 100 × 110 km, including the main river stream and the coastal sea. The greatest depth in the study area is 20 m below mean sea level, near the south-east corner of the modelling domain in the coastal sea. It is important to consider the model stability by satisfying the Courant–Friedrich–Levy number: the selection of time step (Δt) and grid step (Δx) is crucial in order to balance the trade-off between numerical stability and computational time. The simulation grid cell selection for the study area is flexible mesh (FM) or unstructured mesh, where triangular cells of bathymetry are used to optimize the simulation, with small sizes in river domain and larger sizes in offshore settings. The triangular element sizes are about 80–100 m in the river and 800–1000 m offshore, with the total triangular being 35,000 elements and 28,000 nodes. Bathymetry was constructed with the MIKE zero tool and obtained from an updated digitizing map in 2011 from the Southern Institute for Water Resource Planning and DHI Vietnam. The bathymetry and grid resolution is presented in Figure 3.

**Figure 3.** Bathymetry (**left** panel) and computational flexible grid mesh (**right** panel).

### *3.3. MIKE 11 Hydrodynamic Model*

The MIKE 11 is an unsteady 1-dimensional hydrodynamic model, which is based on one-dimensional equations and solves the vertically integrated equations regarding conservation of continuity and momentum. The solution of continuity and momentum equations are employed as an implicit finite difference scheme with a 6-point Abbott scheme [41]. The main governing equations are Saint-Venant equations [42]. In order to effectively simulate floodplain and field areas, and to connect main channels and rivers, we employed the quasi-2D modelling approach developed by Dung et al. [43]. The rice fields and the floodplain are considered to be artificial channels or "virtual channels", having wide cross sections and were extracted from the digital elevation model (DEM). The calibration process of these channels and crest of dikes was explained in detail in [43,44].

#### *3.4. MIKE 21 Hydrodynamic Model*

The MIKE 21 is a dynamic modelling system applicable for coastal and estuarine environments. The MIKE 21 comprises several modules, including the hydrodynamic module, advection–dispersion module, spectral wave module, and transport module. This model is based on the numerical solution of the two dimensional incompressible Reynolds-averaged Navier–Stokes equations with the assumptions of Boussinesq and hydrostatic pressure. It comprises continuity, momentum, temperature, salinity and density equations [45].

This study applies MIKE 21 with two modules, namely hydrodynamic (HD) and transport modules (TR), using a flexible mesh. The mutual interaction between the flow, wind, and velocity and sea tide is considerable. The wind is specified as a spatially constant value for the entire domain and temporally variable values. The schematic presentation of the modelling framework is shown in Figure 4.

**Figure 4.** Schematic of the modelling framework for this study.

The MIKE 21 FM is based on a flexible mesh approach. The spatial discretization of the primitive equation is performed using a cell-centered finite volume method. The spatial domain is discretized by subdivision of the continuum into a non-overlapping cell. The mesh is divided by triangles or quadrilateral elements, as described in Figure 3.

#### *3.5. Model Parameterization*

The model parameters to be defined for the HD model are the roughness coefficient (n), or Manning number, and for the advection-diffusion model, the horizontal dispersion coefficient (Dh), eddy viscosity, Courant-Friedrich-Levy (CFL) and Smagorinsky coefficient. In order to realistically reproduce the physical phenomena of river mechanics, the Manning number is defined as a function of water depth (h) and bed river type, and can be calculated based on depth and drag coefficient [46,47]. In the MIKE 21 manual, the Manning number (M) was defined via the drag coefficient in Equation (1) with g as gravity constant.

$$\mathbf{C\_d} = \frac{\mathbf{g}}{\left(\mathbf{M}\mathbf{h}^{\frac{1}{d}}\right)^2} \tag{1}$$

The Cd can be defined through total water depth (h) and empirical drag coefficient (at 1 m above the bed) for different bottom types (from mud to gravel) by Equation (2):

$$\mathbf{C\_d} = \left(\frac{1}{0.32\mathbf{h}}\right)^{1/7} \cdot \mathbf{C\_{100}}\tag{2}$$

C100 ranges between 0.022 ÷ 0.0047, and a table of full empirical drag coefficients can be found in Soulsby [47]. The Manning coefficients were available for each grid and varied from estuary to upstream river in the entire domain. The range of the Manning values varied from 20 to 40 m1/3/s, and depend on grainsize (rippled sand-gravels) and water depth (shallow to deep water) during model calibration. We adjusted the Manning coefficient through the revision of C100 by defining the type of riverbed and seabed. The eddy viscosity was defined in a Smagorinsky formulation and was adjusted during calibration with a range 0.25–0.27, with an initial value of 0.28.

#### *3.6. Boundary Conditions*

For the delta-wide modelling with MIKE 11, we proposed four scenarios to cover future changes of upstream discharges, sea level rise and in-delta precipitation changes, and used these scenarios as boundary conditions for our modelling exercises (Table 4). More detailed information about rationales and designs of the scenarios can be found in Duong et al. [25]. Results from the Mike 11 modelling, particularly the discharge data at Can Tho, were then used as boundary conditions for the salinity intrusion simulation using MIKE 21. The changes in upstream discharges at Kratie are selected with a range of +10% to −20% relative to baseline discharge in 2011, based on literature review on projected future flow changes of Lauri et al. [48] and Hoang et al. [37]. Notably, both studies [37,48] did not consider the other anthropogenic factors such as irrigated land expansion, urbanization, and inter-basin water transfer. In addition, other studies [5,49] show different ranges of future hydrological changes due to differences in GCM outputs and climate change scenarios selection. Therefore, we extended the range of inflow changes at Kratie to capture possible future hydrological alterations in the upper Mekong Delta. For the climate change scenarios, we selected two RCPs for precipitation, namely RCP4.5 and RCP8.5.

Our sea level rise scenarios were obtained from the scenarios provided by the Ministry of Natural Resources and Environment. We selected an average of the predicted sea level rise, resulting in an increasing 23 cm between 2030 and 2040 and 35 cm between 2050 and 2065 [50]. For the MIKE 21 simulation, the hourly discharges at the Can Tho station were taken from the MIKE 11 simulation with four distinct time series of flows. The predicted sea levels were calculated linearly based on the averaged multi-annual tidal magnitudes from the global tidal model [51]. The seawater density was assumed to be constant and the salinity level at sea was predicted to remain at 35.0 PSU in the future condition. There are two approaches to obtaining the sea level rise condition. One is to impose an instantaneous elevation of sea level at coastal sea boundaries by applying linear adjustment of the mean sea level from the long observed sea level; the other is to run the hydrodynamic model over a long period to obtain predicted sea level. In this study, we obtained the sea level rise projection from the Vietnamese Ministry of Natural Resources and Environment (MONRE) and applied linear adjustment for predicted sea level rise during the 2036–2065 period.

The boundary condition for the MIKE 21 modelling includes the upstream hourly discharge at Can Tho station and tidal magnitude at downstream boundaries in the coast (Figure 1). To obtain salinity boundary conditions under future sea level rise and to reduce the influence of the boundary conditions on the in-delta modelling, a larger offshore domain that covers Hau estuary and adjacent shelf area was used. Secondly, five river branches connecting to the Hau River include Mang Thit, Rach Mop, Cau Quan, Nga Bay and Dai Ngai, were simulated as the sources and sinks. Thirdly, wind speeds and directions data were collected at the Vung Tau hydrological station and applied for the entire domain. Finally, time series of water level at My Thanh, Dai Ngai and Tran De stations were measured during 2010–2011 for model calibration and validation (Table 2). For salinity at downstream boundaries, we set up the salinity concentration as 35.0 PSU and assumed the salinity in upstream to be 0.1 PSU.

**Table 4.** Four scenarios of changing of discharges, sea level rise and precipitation.


#### *3.7. Computation of Flushing Time*

The flushing time can be determined by the freshwater fraction approach [14,52], which can be determined from the salinity distributions. This technique provides an estimation of the time scale over which contaminants and/or other material (saltwater in this study) released in the estuary are removed from the system. Using the freshwater fraction method, the flushing time (*Tf*) in an estuary can be expressed as Equation (3):

$$T\_f = \frac{F}{Q} = \frac{\int f \cdot d(V)}{Q} \tag{3}$$

where *F* is the accumulated freshwater volume in the estuary, which can be calculated by integrating the freshwater volume *d*(*V*) in all the sub-divided model grids over a period of time. In estuaries with unsteady river flow and tidal variations, *F* and *Q* are the approximate average freshwater volume and average freshwater input, respectively, over several tidal cycles for a period of time. The term "*f* " is the freshwater content or the freshwater fraction, which is described by Equation (4):

$$f = \frac{S\_o - S}{S\_o} \tag{4}$$

where *So* is the salinity in the ocean and *S* is the salinity at the study location.

#### **4. Model Calibration and Validation**

#### *4.1. Model Calibration*

The MIKE 21 FM-HD and TR models were first calibrated with the observed averaged water level and salinity concentration at the middle cross section of the Hau River's estuaries. The calibration period was from January to December 2011 for hydrodynamic simulation, and six months from January to June 2011 for salinity simulation. The water levels measured at My Thanh and Dai Ngai were used to calibrate the model. In the inlet of the model at Can Tho station, initial salinity values were set at 0.1 PSU (‰). In the outlet of model (offshore), the salinity concentration was set at 35 PSU.

Several parameters of the HD and TR modules of MIKE 21 were adjusted to improve model performance. The calibration process aims to match simulated results and observed data, including water level and salinity concentration in different locations, by changing the Manning number and CFL in MIKE 21 HD and TR. In the simulations we applied a value of 0.8 for CFL-number and Smagorinsky coefficient within the range 0.25 and 0.27. Model calibration results have achieved the realistic results and indicated that the selected parameters were reasonable, as shown in Figure 5. To compare the observed and simulated water levels and salinity concentrations, the latter were taken at the center point of the cross section, as mean values for comparing with observed values. The performance indices include mean absolute error (MAE), root mean square error (RMSE) and correlation coefficients (R). The comparison between modelling results and observed salinity data also reveals a strong correlation and an excellent prediction, with the R-value higher than 0.8 (Table 5).

Figure 6 compares the observed and simulated salinity concentrations at the Tran De and Dai Ngai stations. The simulated salinity time series compared favorably with the discrete salinity measurements at the two aforementioned stations. Overall, the model reflected the large dynamic variation of salinity between 0 and 30 PSU over a tidal cycle, with decreasing values of mean salinity as freshwater discharge increased.

**Table 5.** Statistical performance of calibration and validation at My Thanh and Dai Ngai stations.


**Figure 6.** Calibration of salinity concentration at Tran De and Dai Ngai stations.

#### *4.2. Model Validation*

The MIKE 21 model was further validated with the calibrated parameters, focusing on water levels and salinity concentration. The time series data of hourly discharge at the Can Tho station in 2010 was used as upstream boundary conditions to drive the model simulations. The hourly tidal values taken from a global tide model were adopted as a forcing function at the coastal sea boundaries. The hourly tidal data and daily freshwater discharges were collected from SIWRR and the Southern Regional Hydro-Meteorological Center (SRHMC). The water level and salinity concentrations from Tran De and Dai Ngai stations were employed to evaluate the model. The comparison of observed data and simulated results for water level and salinity concentration is verified to check the model's performance.

Figure 7 compares the simulated water levels and observed data with the time series at My Thanh and Dai Ngai during the period between 29 April and 11 May 2010. In general, the modelling results show realistically simulated water level variations. The comparison demonstrates the model's capability to reproduce the water levels, even under large variations of daily freshwater influx from the upstream Can Tho station. To compare the observed and simulated water levels and salinity concentrations, the simulated water level and salinity were taken at the center point of the cross section, as mean values for comparing with observed values. Overall, the model satisfactorily simulated the water level at My Thanh and Dai Ngai on the Hau River. The calibrated model parameters were, therefore, adopted for our modelling exercises and scenario analyses.

**Figure 7.** Validation of water level at My Thanh (**a**) and Dai Ngai (**b**) stations.

#### **5. Results and Discussion**

#### *5.1. Changes in Salinity Intrusion*

Simulating the spatial variations of salinity in the estuary and further upstream of the Hau River shows the detailed changes in the salinity dynamics under different discharges, sea level rise (SLR) and rainfall scenarios. The discharge at Can Tho varied significantly throughout the year, with a maximum discharge less than 20,000 m3/s in the wet season and a minimum discharge of −15,000 m3/s (Figure 8). The inverse flow direction in the dry season is caused by the tidal flow. Discharges at Can Tho were influenced not only by the upstream flow from Tan Chau, Vam Nao, but also by tidal regimes from the East Sea. The tidal feature in South China Sea is semidiurnal asymmetry; the peak spring tide reaches 3.0 m between December and January every year and reaches its lowest levels from June to August with a variability of around 0.5 m.

**Figure 8.** River discharges at Can Tho simulated from the MIKE 11 model for four scenarios.

In order to quantify the spatial variations in the salinity concentrations along the river, the minimum and maximum difference in salinity concentrations between the baseline and four scenarios was calculated and is presented in Table 6. Figures 9 and 10 presents the spatial distributions of maximum and average salinity levels along the Hau River estuary under the four scenarios, respectively. The result shows radical changes in the salinity levels across the entire estuary. Compared to the reference isohaline of 4.0 PSU, the salinity profile with 4.0 PSU moves farther upstream by 48.55, 49.13, 49.16 and 49.18 km from Scenario 1 to Scenario 4, respectively (Table 7). The relative changes compared to the baseline are 3.29, 3.87, 3.90, and 3.92 km for the four scenarios, respectively (Table 8). Scenario 4 shows the farthest salinity intrusion, which is explained by strong upstream discharge reduction and substantial sea level rise of 35 cm (Figure 11). Similarly, Scenario 1 shows weaker salinity intrusion from the sea compared to other scenarios. This is explained by slight increases in the upstream discharge and the sea level rise of 25 cm. Figure 12 presents the salinity distribution under four scenarios at different points along the river: from the river mouth (My Thanh station) to upstream (Can Tho station) in spring tide (Figure 12a) and neap tide simulations (Figure 12b). The mean values are taken at the middle point of the cross-section. Figure 9 shows the difference in salinity between the four scenarios and the distance of salinity intrusion from the river mouth.

Here we discuss our results for salinization modelling in view of the relevant studies on this topic. Smajgl et al. [53] studied the effects of a wide range of driving factors on salinization, including land use changes, sea level rise of 30 cm, development of all proposed upstream reservoir and irrigation, and an increasing number of dry years. The results stated that the isohaline of 4.0 PSU is relocated throughout the Hau River, with distance from the river mouth reaching approximately 70 km. Considerable differences in the 4.0 PSU isohaline between this study and Smajgl et al. [53] can be attributed to differences between the scenario setups and boundary conditions in the two studies. Nguyen and Savenije [32] tested analytical solutions of salinity simulation in the Mekong estuaries based on observed data in 2005 for the Co Chien, Cung Hau and Hau estuaries. The results showed that salinity intrusion distances from the river mouth were 41 km and 23 km for spring tide and neap tide during the dry season in 2005, respectively. Furthermore, the finding by Trieu and Phong [54] concluded that the salinity intrusion of 1.0 PSU in Hau River was approximately 55–60 km for the dry season of 2010. This projection is higher than our results for all scenarios. All in all, results from this study show a strong dependency of salinization on upstream inflow and the changing tidal dynamics

under sea level rise. Such a mechanism is also commonly reported in large river deltas of the world, including the Bangladesh Delta [55,56], and the Dutch Delta [57].

**Figure 9.** Spatial distribution of maximum salinity level across the modelling domain under Scenario 1 (**a**), Scenario 2 (**b**), Scenario 3 (**c**) and Scenario 4 (**d**).


**Table 6.** Maximum and minimum salinity concentration at Dai Ngai and Tran De stations.

**Table 7.** Salinity intrusion length for thresholds of 1.0 PSU and 4.0 PSU for four scenarios, under spring tide and neap tide.



**Table 8.** The relative change of salinity intrusion length for four scenarios compared to baseline.

**Figure 10.** Spatial distribution of average salinity concentration across the modelling domain under Scenario 1 (**a**), Scenario 2 (**b**), Scenario 3 (**c**) and Scenario 4 (**d**).

#### *5.2. Changes in Salinity Intrusion Length*

We further assessed the changes in salinity intrusion length under the four future scenarios. The salinity intrusion length is defined as the distance from the estuary mouth to the location upstream with mean salinity at the cross section. There are three types of intrusion length: intrusion length at low water slack, intrusion length at high water slack, and tidal average intrusion length, which is considered to be an average of low and high [58–60]. The salinity intrusion length is defined in this study as the distance from the Hau River mouth to the upstream limit location where the bottom salinity level drops down to a certain threshold, e.g., 1 PSU or 4 PSU. Figure 12 presents salinity distribution along the river for four scenarios, under spring tide and neap tide.

**Figure 11.** Differences of maximum salinity concentration between future scenarios and baseline, for Scenario 1 (**a**), Scenario 2 (**b**), Scenario 3 (**c**) and Scenario 4 (**d**).

**Figure 12.** Salinity distribution at center line of cross section along the river from the mouth to upstream during spring tide (**a**) and neap tide (**b**).

Table 7 shows that the salinity intrusion lengths under spring tide are larger than those under neap tide. By comparing the baseline to the four scenarios, it can be seen that the intrusion lengths for all scenarios are consistently farther than that of the baseline. In Scenario 1, river discharge increases by 10% at Kratie, and therefore the inflow at Can Tho is also slightly higher than the other scenarios. This explains the difference in salinity intrusion length between the baseline and Scenario 1. For Scenarios 2, 3 and 4, (with sea level rises of 23 cm, 35 cm and 35 cm, respectively) the differences in salinity intrusion lengths between the baseline and these three scenarios are relatively small. In essence, the results show that upstream discharge changes do not substantially affect salinity intrusion lengths, in contrast to the more dominant impacts of tidal regime.

Under sea level rise scenarios, salinity intrusion length in the Hau River would increase between 3.29 and 3.92 km (for 1.0 PSU salinity level) compared to the baseline in the spring tide condition. During neap tide, this length ranges between 4.36 and 4.65 km for the 1.0 PSU salinity level (Table 8). Scenario 4 shows the longest salinity intrusion length for spring and neap tide (Table 7) for both 1.0 PSU and 4.0 PSU isohaline. Such variation in the magnitude of salinity intrusion length is caused by changes in the stratification of the estuary. As the sea level rises, the river depth increases and the horizontal gradient of salinity changes, resulting in increased estuarine circulation. The model simulations indicate that the location of 4 PSU isohaline would migrate up from 43.58 to 44.27 km if sea level rises from 23 to 35 cm and the upstream discharges change from +10% to −20% at Kratie. Our findings about increasing salinity intrusion length are in line with observed data, which confirms that when river discharge reduces in the Hau River, saline water moves upstream to Can Tho [54,61]. Any increase in the magnitude or duration of salinity intrusion at the location of the irrigation water intake would likely affect crop yield and aquaculture in Mekong Delta.

#### *5.3. Changes in Flushing Time*

To examine the effect of river discharge changes on the temporal dynamics of salinity intrusion in the Hau River, we calculate flushing time with varying discharge conditions for baseline and four future scenarios. Flushing time is defined as the time required to replace the existing saltwater in the estuary using inflows at certain discharge levels [14,62]. To calculate the flushing time in the estuary, different freshwater discharges were selected based on maximum and minimum discharges from the dry season in 2011, with a range of 400 to 4500 m3/s. The simulated flushing times for the baseline and the four scenarios were calculated and are presented for different discharge values in Figure 13.

**Figure 13.** Estimated flushing time under varying discharge levels (Q, m3/s) for four scenarios.

The flushing time varies considerably under all considered scenarios, ranging between 10 h (at discharge 4500 m3/s) and 109 h (at discharge 400 m3/s). For the most extreme scenario (i.e., Scenario 4), flushing time varies within a range from 20.7 h to 182 h. The results also indicate that, for high flow, the flushing time under Scenario 1 is lower than under baseline conditions, while Scenarios 2, 3 and

4 consistently exhibit longer flushing time than the baseline. Longer flushing time caused by reducing discharges means that it would require a higher volume of freshwater to push saltwater back to the river mouth in the future. Reductions in river discharge also cause the salinity boundary to move further upstream (Figure 11). Salinity deceases substantially when river discharge increases, and as a result the salinity intrusion boundary migrates farther downstream. As freshwater discharge increases from 400 to 4500 m3/s with different sea levels, the flushing time reduces significantly from 87.72 h to 10 h for Scenario 1, and 150 h to 20.7 h for Scenario 4.

#### *5.4. Uncertainties, Limitations and Further Research*

The hydraulic models, including 1D-MIKE 11 and 2D-MIKE 21 HD and TR used in this study, entail limitations and assumptions. These limitations and assumptions exist in both the data and modelling approach. The major data used in this study were predicted river discharge, rainfall scenarios under climate change, and sea level rise. The climate change and sea level rise scenarios used in this study were obtained from different sources, such as GCMs outputs, SRHMC and SIWRR, and therefore our results might be subjected to errors caused by data inconsistencies. However, we implemented quality checks on our input data and controlled for inconsistencies. Additionally, we acknowledge that climate change is a non-stationary and dynamic process; while the projected salinity intrusion from upstream discharge changes, rainfall and sea level rise alterations were treated in this study as part of a steady-state system. That means we assume that the relationship between these driving factors are also valid in the future and did not include their interrelations in the simulations. This assumption could lead to biases in future discharge and rainfall estimates, which in turn result in uncertainties in the simulated salinity intrusion. However, we think that the hydraulic and salinity dynamics in the Hau River of the Mekong Delta are primarily driven by changes in the considered driving factors and their interrelations play marginal roles. Therefore, the developed modelling approaches are useful and applicable for our modelling purposes.

Furthermore, we employed a two-step modelling procedure (see Section 3—Methodology and Model setup) for assessing changes in salinity intrusion in the Mekong Delta. A 1D model was applied to simulate delta-wide hydrodynamic conditions (i.e., discharge and water level), thereby providing boundary conditions for the more detailed hydrodynamic and salinity simulation with the 2D-MIKE 21 model. While the current data availability and computational capacity does not allow for detailed 2D simulation for the whole Mekong Delta domain, we suggest gradually increasing the domain size to cover larger areas in future studies. Last but not least, while our calibrated parameter set yields realistic simulation results for the current state of the delta system (i.e., 2010–2011), we did not consider changing delta settings and its hydraulic characteristics in the future. Therefore, we suggest focusing on this topic in future studies.

#### **6. Conclusions**

This study combined 1D-MIKE 11 and 2D-MIKE 21 hydrodynamic (HD) and salt transport models (TR) to simulate the hydrodynamics and salinity distributions in the Hau (Bassac) River estuary of the Mekong Delta, southern Vietnam. The model was calibrated and verified using observational water level, salinity distribution and tidal data from 2010 to 2011. The model simulation results agree well with the observed data during calibration and validation periods. The calibrated model was used to simulate future salinity intrusion, focusing on the potential impacts of upstream discharges, rainfall variabilities and sea level rise on salinity intrusion length and saltwater flushing time. The simulation results indicate that a combination of upstream discharge reductions, rainfall changes and rising sea level will substantially exacerbate salinity intrusion in the Mekong Delta. By the 2036–2065 period, salinity intrusion moves farther upstream by between 4855 and 4918 m, depending on the scenario considered. As a result, the flushing time required to re-establish freshwater conditions in Hau River will increase in the future. While flushing time under baseline ranges between 125 h and 13.5 h, it tends to increase in the future, ranging between 180 h and 10 h. Increasing salinization in the future will

make it more difficult to re-establish the freshwater condition in the estuary. In particular, the flushing time required to replace saltwater with freshwater at the estuaries tends to increase between 7.27 h for maximum discharge of 4500 m3/s and 58.95 h for discharge of 400 m3/s under the most extreme scenario. Increasing salinity intrusion under upstream discharge changes, rainfall alterations and sea level rise, will likely have serious consequences for crop production, freshwater supplies and freshwater ecosystems, therefore requiring timely adaptation responses.

**Author Contributions:** D.T.A. designed the study, processed and analyzed the data, developed the models, interpreted the results and wrote the paper. L.P.H. provided data, assisted in the data analyses and drafting the manuscript. The study has been carried out under the supervision of M.D.B. and P.R., who contributed to the model development stage with theoretical considerations and practical guidance, assisted in the interpretations and integration of the results and helped in preparation of this paper with proof reading and corrections.

**Funding:** This research was funded by the CUOMO Foundation and the IPCC Scholarship Program, grant number IPCC2013-2017.

**Acknowledgments:** The contents of this paper are solely the liability of the authors and under no circumstances may be considered as a reflection of the position of the Cuomo Foundation and/or the IPCC. The German Research Foundation (DFG) and the Technical University of Munich (TUM) supported this work in the framework of the Open Access Publishing Program. We greatly appreciate their help. Wholehearted thanks are given to DHI Vietnam, Hanoi water resource University, base 2, Southern Regional Hydro-meteorological Center, National meteorological center for providing all necessary data. Special thanks to Dang Quang Thanh and Nguyen Van Song on for valuable support and comments to complete this paper. We thank the editor and three anonymous reviewers for their constructive comments, which helped us to improve the manuscript.

**Conflicts of Interest:** The authors declare no potential conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
