Flow Division Dynamics in the Mekong Delta: Application of a 1D-2D Coupled Model

Abstract

The Mekong Delta constitutes a complicated multi-channel estuarine system, exchanging water with a delta-wide irrigation system. A 1D–2DH coupled numerical domain is calibrated and validated for water level and discharge during the dry season. This approach benefits from the simplicity of a 1D network within the estuarine and irrigation systems, while maintaining the interaction with the spatial tidal dynamics of the 2DH coastal domain. First, the role of the irrigation system on tidal dynamics is quantified; then, tidal propagation, freshwater budget, and the effect of offshore subtidal water level on discharge division are investigated. The results show that the complex irrigation system, in a friction-like manner, reduces the tidal amplitude up to 25%. The channels aggregate to 1% of the total water volume in the delta, while accommodating up to 10% of the tidal prism. Tidal amplitude reduces upstream, while subtidal water level is highly sensitive to upstream discharge, spring–neap cycles, and wind-generated offshore surge. Although cumulative discharge division within the estuarine network is consistent, temporal discharge division can be significantly sensitive to offshore wind-surge. During the dry season, it can reverse the expected subtidal discharge division within the time-scale of a few days and potentially influence salt intrusion.

1. Introduction

The Vietnamese Mekong Delta (VMD) covers approximately 60,000 km² of land [1], is home to 18 million people, contributes to 26% of Vietnam’s national GDP [2], and provides half of the nation’s food supply [3,4]. The delta drains the Mekong catchment of approximately 760,000 km², along 4910 km of the Mekong River [5,6]. The Mekong River within the VMD transforms to a multi-channel estuarine system draining to the coastal waters. A distinct feature of the Mekong estuarine system is the presence of thousands of man-made small irrigation/navigation channels, exchanging water with the main rivers. Globally, the VMD is one of the lowest-lying deltas [7], making it disproportionately vulnerable to various environmental and anthropogenic forces (see digital elevation map in Figure 1a). Salt intrusion [8], along with land subsidence [9], flooding [10,11,12], and coastal erosion [1], are identified as hazardous processes risking nation’s food safety and strategic development of the delta [2,13,14]. Different development plans and investments that facilitate the rapidly urbanizing delta are being implemented [8], and numerous modeling efforts are being carried out to support various phases of these investments. Despite many field and numerical studies, there is not yet a comprehensive large scale understanding of physical processes governing tidal and subtidal (averaged over a tidal cycle) flow division dynamics of the multi-channel estuarine system of the VMD, while it has been shown that this is of importance to understand for example the freshwater distribution in the system.

The Mekong River, in terms of river discharge, shows significant seasonal behavior. Discharge in Kratie, Cambodia, varies by as much as 40,000 m³/s to 3000 m³/s from the wet season (SW monsoon season, June to November) to the dry season (NE monsoon season, December to May). The source of freshwater to the VMD alters between wet and dry seasons. The Tonle Sap lake (TSL), the largest body of water in the mainland of South East Asia, plays a unique and a slightly complicated role in freshwater supply to the VMD [6,16]. During the wet season, Tonle Sap river that links TSL and the Mekong River drains the Mekong River, fills the lake, and floods the surrounding lowlands when the wet surface area of the lake rises to as large as 14,500 km². Therefore, in the wet season, the Mekong River is the main source of freshwater to the VMD. During the dry season, when the Mekong River level is dramatically lower, the Tonle Sap river reverses and drains the lake into the VMD, during which the area of the lake shrinks down to 2500–3500 km² [6,16].

In response to discharge variations, tidal propagation and salt water intrusion change seasonally. Tidal propagation varies between nearly 200 km in the wet season to almost 400 km in the dry season [17]. Saltwater intrusion changes from a few kilometers in the wet season to tens of kilometers in the dry season [18,19]. The internal dynamics of the multi-channel estuarine system, its interaction with nearshore hydrodynamic processes, and lack of knowledge on the role of the irrigation system and water demand make it utterly difficult to have a comprehensive image of the freshwater budget. This becomes increasingly challenging in the dry season when limited freshwater resources are vital to agriculture, and tidal dynamics play a crucial role on saline water dynamics and water quality. Figure 1b–e shows the monthly averaged tidal difference, measured at four stations (averaged over 2008 to 2013). In Chau Doc and Tan Chau, tidal difference is reduced nearly 80%–90% during the wet season, emphasizing the significance of tidal forces in the dry season.

Within the VMD, the irrigation system has a hierarchy level of primary and secondary to tertiary and quaternary levels of channels, where even the largest irrigation channels have cross-sectional areas hundreds of times smaller than the main rivers. The majority of numerical studies in the VMD were carried out on either 1D domains [8,20,21,22], where the complicated network of primary and secondary channels were included in the model domain, or in 2DH (averaged over depth) or 3D e.g., [23,24] where the main rivers benefit from higher resolution (in horizontal and/or vertical) and the coastal waters were included but the irrigation system was neglected. General consensus among Vietnamese experts is that the irrigation system (see Figure 2) plays a crucial role in hydrodynamics of the delta, but we were not able to find a study that systematically demonstrates the effect of this system on tidal and freshwater dynamics of the VMD.

There is limited agreement among the existing literature on freshwater division within the estuarine network of the VMD and the physical processes influencing that. Various domestic studies estimated freshwater division during the 1990s [25]. Nguyen et al. [25] estimated discharge division of seven main downstream distributaries, using a salt intrusion analytical model [26], matching those derived by some of the older domestic numerical models. More recently, freshwater distribution of lower Song Hau distributary channels was estimated based on a 25-hour field measurement [19], and a Delft3D (2DH) model [23]. Since the two later studies were focused on the Song Hau channel (see Figure 1), a summary of freshwater distribution estimates between Tran De (TD) and Dinh An (DA) are presented in

Table 1. Summary of recent studies addressing freshwater division between the downstream tributaries of the Song Hau channel during the dry season.

Study	Fraction of Song Hau Discharge		Approach
	Dinh An (DA)	Tran De (TD)
Nguyen et al. [25]	70%	30%	Analytical salt intrusion model [26]
Nowacki et al. [19]	40%	60%	25-h field measurements (calculating percentage of DA+TD discharge)
Xing et al. [23]	94%	6%	2DH numerical modeling of lower Song Hau distributary channel

, demonstrating inconsistency between various estimates. In addition, none of these studies addressed the temporal changes in freshwater distribution over the dry season, potentially in response to oceanic forcing or spring–neap cycles that are common in estuaries (e.g., see References [27,28,29]).

In earlier studies, freshwater division in rivers was related to hydraulic roughness and the downstream channel geometries [30]. Bathymetrical asymmetries upstream of the river junction, such as channel direction or bends can also impact flow division between the branches [31,32,33]. Furthermore, it has been shown that subtidal discharge division is sensitive to asymmetry in downstream channel properties [34] (e.g., cross-sectional area, friction), and bifurcation angle [29].

Downstream processes can also influence discharge division. During low discharge regimes, tides can considerably influence river flow variations, and the low-frequency subtidal water level variations can disturb the flow regime beyond tidal penetration [27]. Fortnightly oscillation of discharge division can be driven by imbalance in Stokes transport and return flow [29]. Furthermore, subtidal waves generated in the ocean can attenuate or amplify and travel hundreds of kilometers upstream [28]. Since river flow can fluctuate in an abrupt and event-like manner [35], it is verified [36] that tides can modulate discharge division by oscillatory and subtidal water level variations, while changes in river flow can in turn obscure the oscillatory water level gradients. However, the effect of subtidal offshore water level variations, generated independent of tides (e.g., as in Reference [28]), on flow division is not yet addressed in the literature.

This research is set out to investigate tidal propagation and freshwater distribution as well as underlying barotropic processes during the dry season in the VMD. To that end, we (a) study the tidal dynamics of the VMD and examine the role of primary and secondary irrigation channels in tidal dynamics of the estuarine system, (b) calculate the freshwater budget/distribution in the VMD, and (c) study the effect of offshore subtidal water level on temporal variations in freshwater distribution during the dry season.

2. Materials and Methods

2.1. Approach

A 1D–2DH coupled model in DFlow-FM (Figure 3) was calibrated against a widespread water level and discharge dataset of the dry season of 2010 and validated against the dry season of 2012. The overall model performance, as well as tidal and subtidal dynamics of water level and discharge were scrutinized to ensure an acceptable reproduction of various physical processes. Then, to address the processes mentioned in the last section, the following steps were taken:

(1)

Studying the role of the irrigation channels:

(a)

For sensitivity analysis, the irrigation channels were excluded from the 1D domain to see the immediate effect on tidal propagation and prism.

(b)

By increasing roughness along different channels, the model was re-calibrated to quantify the irrigation system contribution to tidal propagation and prism.

(2)

The complete domain was then used to study cumulative freshwater distribution during the dry season.

(3)

Temporal variation of subtidal water level and discharge was then related to other physical processes such as wind and surge level in the sea.

(4)

To isolate the role of different physical processes, we exclude all the sink terms (water demand and evaporation) and zoomed in on the Hau River and its two distributaries (see Table 1 and Figure 1a). Following a method of factors separation (e.g., see Reference [29,34]), various forcings in the model were one-by-one included or excluded to examine their individual impact and their interaction with other forcings.

2.2. Model Description

DFlow-FM solves the shallow water equations, under Boussinesq assumptions, with a finite volume method on staggered unstructured grids in 1D, 2DH, 2DV, or 3D schematizations [37]. The model provides the possibility of benefiting from a domain represented by 1D network connected to polygon-shaped cells of arbitrary degrees (such as triangles, quadrilaterals, pentagons, or hexagons). The model solves the depth averaged homogenous shallow water equations [37]:

$$\frac{\partial \mathrm{H}}{\partial \mathrm{t}}+\nabla \left(\mathrm{H}\overrightarrow{\mathrm{u}}\right)=\mathrm{q}$$

(1)

$$\frac{\partial \overrightarrow{\mathrm{u}}}{\partial \mathrm{t}}+\mathrm{adv}\left(\overrightarrow{\mathrm{u}}\right)+\mathrm{g}\overrightarrow{\nabla }\mathsf{\zeta }+\frac{1}{\mathrm{H}}\frac{\mathrm{g}}{{\mathrm{C}}^2}\overrightarrow{\mathrm{u}}\parallel \overrightarrow{\mathrm{u}}\parallel +2\mathsf{\Omega }\times \overrightarrow{\mathrm{u}}=\mathrm{d}.$$

(2)

Here, $\mathrm{H}$ denotes total water depth, $\mathsf{\zeta }$ is water level, $\overrightarrow{\mathrm{u}}$ is depth-averaged horizontal velocity, $\overrightarrow{\nabla }$ is the horizontal gradient operator, $\mathsf{\Omega }$ is the earth rotation vector, $\mathrm{adv}\left(\overrightarrow{\mathrm{u}}\right)$ is the advection term, q is a source/sink term, d is external forcing, C is Chézy coefficient, and g is gravitational acceleration. Time integration in both 1D and 2DH grids are done implicitly, except for the advection term, $\mathrm{adv}\left(\overrightarrow{\mathrm{u}}\right)$ in Equation (2). The advection term is integrated explicitly in time in a momentum conservative fashion in a finite volume sense [38]. This results in a symmetric system for the water levels that is efficiently solved by Gaussian elimination, conjugate gradient, and back substitution (following Reference [38]). The advantage of this procedure is that the time-step restriction for numerical stability applies to flow speed instead of gravity wave celerity, which would impose more severe restrictions, especially in deeper water. The scheme uses an identical approach in 1D and 2DH grids. The 1D and 2DH grid cells are solved in the same system of equations to ensure a smooth 1D–2DH coupling. Apart from neglecting cross-channel flows, the difference between 1D and 2DH is that in 2DH, cross-sectional area is based on X-, Y-, and Z-coordinates of the network, while in 1D, these values are based on profile definitions, interpolated in between the user-defined cross-sectional profiles. We refer to Reference [37] for further details on model description.

2.3. Model Set Up

Table 2. Summary of input data of measurement stations within the VMD (Chau Doc, Tan Chau, Vam Nao, Can Tho, My Thuan) and Cambodian stations (Prek Kdam, Koh Khel, Phnom Penh, Neak Luong).

Input data	Source	Description
River network, channel geometry and profile definitions	SIWRP	Mainly from 2008 and updated when possible
Lower Hau River bathymetry	from Nowacki et al. [19]	multiple cross-sections along Hau River and DA channel
2DH model grid and bathymetry of TSL and the continental shelf	Vo et al. [24]	published in 2018
Water level data within VMD	SRHMC	Hourly time-series
Discharge data within VMD	SRHMC	Hourly time-series
Water level at Cambodian stations	SIWRP	Once a day at 7:00AM
Discharge at Kratie	MRC	Daily averaged discharge
Offshore wind	Climate Forecast System Reanalysis (CSFR), originally published by NCEP NOAA downloaded from online DHI metocean data portal	The hourly wind parameters at 10 m
Evaporation	Literature [16,42]	Daily [mm/day], in consultation with SIWRP
Water demand	SIWRP	Per region [m3/s]
Offshore Tidal constituents	TOPEX/Poseidon global inverse tide model (TPXO 8.0)	Amplitude and phase of tidal constituents (M2, S2, N2, K2, K1, O1, P1, Q1, MF, MM, M4, MS4, MN4)

summarizes and describes the input data and their sources. Figure 3 shows the 1D–2DH domain and bathymetry of the VMD with reference to Hondau datum (official Vietnamese benchmark system, identical to mean sea level). From upstream, the model is bound to Kratie and the TSL (see Figure 3). In downstream direction, the river distributaries drain into the sea and the model boundary is approximately 70 km offshore. The model includes the network of primary and secondary channels, neglecting tertiary and quaternary levels. All the rivers and channels are modeled in 1D, while the TSL and the coastal sea are in 2DH. The simulations were run from 1 January to 25 April in 2010 and 2012, with the first 15 days considered as spin-up time.

From upstream, the model is mainly forced by measured daily river discharge. During the dry season, freshwater input to VMD is controlled by discharge in the Mekong River and contribution from the TSL [39,40]. TSL water level at the end of the wet season varies between years and it is core to the absolute amount of freshwater that gradually drains into the VMD during the dry season. To accommodate for the behavior of the TSL, initial water level of TSL was set to the observed daily water level of Prek Kdam (see Figure 2) on 1 January of the simulation year. In addition, evaporation was gradually increasing from 4 to 7 mm/day from January to April. Note that evaporation is only from the open surface waters and transpiration was not incorporated. Furthermore, extracted groundwater perhaps partly evaporates and partly adds to the surface water system. However, the average of 2.5 million m³/day groundwater pumping [9] compared to almost 85 million m³/day of water demand is negligible.

From downstream, the model is forced by spatially varying astronomic components at the alongshore seaward boundary (20 boundary points). At the two cross-shore boundaries, Neumann conditions were applied. Apart from the tidal constituents, time series of subtidal water level, derived from Godin-filtered [41] hourly water level at Vam Kenh (See Figure 2), after correction for phase lag, were superposed at the offshore boundary. This approach proved able to effectively include the wind-induced surge level in the model (see calibration and validation results in Section 2.5).

Figure 3 shows five colored polygons each presenting a sub-region for water demand (modeled as sink). The demand discharge is uniformly sunk from the open surface of the model within the corresponding polygon (dashed line within each area). This is an intuitive measure to extract water gradually from the system instead of lateral discharges at multiple points that may abrupt local water levels. A total of 52 dams, primarily containing the irrigation system from saline water intrusion, were included in the model that remain closed throughout the dry season (see Figure 2).

2.4. Time Series Analysis

Hourly time series of water level and discharge can be decomposed to tidal and subtidal (tidally averaged) components. To derive the subtidal component, hourly data were de-tided with the low-pass Godin filter [41]. To separate tidal components, harmonic analysis was carried out by applying a least square fit of tidally varying time series (subtidal component subtracted) to a tidal function including independent tidal constituents (satisfying Rayleigh criterion) [43]. Cumulative discharge was calculated by integrating subtidal hourly discharge between 15 January and 25 April of the target year. When addressing discharge division within the VMD, cumulative discharge is presented as percentage of total freshwater inflow (sum of cumulative discharge of Tonle Sap and Kratie between 15 January and 25 April). Daily discharge contribution of Tonle Sap lake could be approximated by subtracting the observed Mekong (i.e., Kratie) discharge from daily sum of Chau Doc and Tan Chau.

We decompose the subtidal discharge time-series into Eulerian discharge and Stokes transport (see References [29,34,44,45]). The cross-sectional area $a$ and cross-sectionally averaged along-channel velocity v were decomposed into Godin-filtered subtidal (denoted with〈 〉) and tidally varying (denoted by ’) components:

$$\mathrm{a}=\langle \mathrm{a}\rangle +{\mathrm{a}}^{\prime }\&\mathrm{v}=\langle \mathrm{v}\rangle +{\mathrm{v}}^{\prime }.$$

(3)

The subtidal discharge $\langle \mathrm{q}\rangle$ can then be decomposed:

$$\langle q\rangle = \langle av\rangle = \underset{A1}{\underbrace{\langle a\rangle \langle v\rangle }} + \underset{A2}{\underbrace{\langle a\rangle \langle v^{\prime }\rangle }} + \underset{A3}{\underbrace{\langle v\rangle \langle a^{\prime }\rangle }} + \underset{A4}{\underbrace{\langle a^{\prime }{v}^{\prime }\rangle }} \simeq \underset{Q_{Residual}}{\underbrace{\langle a\rangle \langle v\rangle }} + \underset{Q_{Stokes}}{\underbrace{\langle a^{\prime }{v}^{\prime }\rangle }}.$$

(4)

A2 and A3 tend to be insignificant (the average of tidally varying components are very small, but not the average of their product). A1 is the (Eulerian) residual/subtidal discharge, which includes river discharge and Stokes return flow. A4 is the Stokes transport.

To demonstrate M₂ tidal amplitude variation in different simulations relative to the reference simulation, the amplitude change reads as:

$${\mathrm{M}}_2 \mathrm{amplitude} \mathrm{change} = \frac{{\mathrm{M}}_{2,\mathrm{Sim}.}\left(\mathrm{x}\right) - {\mathrm{M}}_{2,\mathrm{Ref}.}\left(\mathrm{x}\right)}{{\mathrm{M}}_{2,\mathrm{Ref}.}\left(\mathrm{x}\right)} \times 100$$

(5)

where ${\mathrm{M}}_{2,\mathrm{Sim}.}\left(\mathrm{x}\right)$ and ${\mathrm{M}}_{2,\mathrm{Ref}.}\left(\mathrm{x}\right)$ are the ${\mathrm{M}}_2$ tidal amplitude of the modified and reference simulations, respectively, as a function of distance from the estuary mouth.

The variation of tidal prism in simulations in comparison to the reference simulation are calculated as:

$$\begin{array}{ll}\mathrm{Variation} \mathrm{of} \mathrm{tidal} \mathrm{prism}& = \frac{\mathrm{T}{P}_{\mathrm{Sim}.} - T{P}_{\mathrm{Ref}.}}{T{P}_{\mathrm{Ref}.}} \times 100\\ & = \frac{{{\displaystyle \int }}_{0h}^{24h}{Q}_{\mathrm{flood}, \mathrm{Sim}} - {{\displaystyle \int }}_{0h}^{24h}{Q}_{\mathrm{flood}, \mathrm{Ref}.}}{{{\displaystyle \int }}_{0h}^{24h}{Q}_{\mathrm{flood}, \mathrm{Ref}.}} \times 100\end{array}$$

(6)

where, ${\mathrm{TP}}_{\mathrm{Sim}.}$ and ${\mathrm{TP}}_{\mathrm{Ref}.}$ are daily tidal prism of the modified and reference simulations, respectively. Tidal prism of reference and modified simulations are calculated as the sum of flood discharge (${\mathrm{Q}}_{\mathrm{flood},\mathrm{Ref}.}$ and ${\mathrm{Q}}_{\mathrm{flood},\mathrm{Sim}.}$) at the mouth of seven estuaries.

2.5. Mekong Model Calibration and Validation

Figure 2 shows the Mekong model network and the observation points used for calibration and validation. Both tidal and subtidal water level and discharge representation throughout the model domain were evaluated qualitatively and quantitatively, by applying model performance indicators. The Pearson correlation coefficient R² was used to represent the linear relation between observations and the model results. The Nash and Sutcliffe coefficient of efficiency (NSE, [46]) was used as an indicator of goodness of fit, representing the ratio of residual variance (between model and observation) and observation standard deviation, see Equation (7):

$$\mathrm{NSE} = 1 - \frac{{{\displaystyle \sum }}_{i=1}^N{\left(O_i-M_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^N{\left(O_i-\overline{\mathrm{O}}\right)}^2}$$

(7)

where $\mathrm{N}$ is the number of observations, ${\mathrm{O}}_{\mathrm{i}}$ and ${\mathrm{M}}_{\mathrm{i}}$ are observed and modeled values, respectively, and $\overline{\mathrm{O}}$ is the mean of the observations.

The main calibration variable was spatially varying bed roughness, represented by Manning coefficient, as it was also adopted by previous studies [20,24,47]. It spatially varied between 0.018 and 0.04 [48], based on recent bed material distributions [17]. In the coastal waters and over the ebb-tidal deltas, the Manning coefficient is between 0.018–0.02, increasing to 0.022 towards Can Tho and My Thuan, with the exception of Ham Luong, Tieu, and Dai having a slightly higher roughness of 0.021. Upstream from Can Tho and My Thuan, roughness increases up to 0.035 in Phnom Penh and remains constant for the rest of the upstream domain. Given the high traffic (navigation), population, vegetation, and residential construction density along the banks, we expected higher roughness in the primary and secondary channels and adopted a constant Manning coefficient of 0.04.

Figure 4 shows a qualitative comparison of modeled hourly and subtidal water level against daily-averaged values at four Cambodian stations. In 2010 and 2012, Tonle Sap lake contributed to, respectively, 33% and 38% of the total Mekong Delta freshwater inflow. Modeled hourly discharge of the Tonle Sap river at Prek Kdam (and its subtidal component) shows a satisfactory comparison to daily discharge of Prek Kdam (Figure 4b). Note that from mid-March onwards, there is hardly any freshwater inflow from the lake to the VMD.

Figure 5 shows the time-series of modelled and observed total and subtidal water level in the left panels at seven different stations. Note that subtidal water level variations (despite slight deformation) propagate as far as 200 km upstream; e.g., compare water level at Vam Kenh (estuary mouth) and Tan Chau (Cambodian border). A nearly perfect comparison of subtidal modeled and observed water levels show that this behavior is very well captured in the model. In the right panels, the amplitude (left axes) and phase (right axes) of tidal components of observed and modeled water levels are distinguished. Overall, the most important diurnal, semi-diurnal, and their higher harmonics are simulated well by the model.

In Figure 6, hourly tidal and subtidal discharge are illustrated in the left panels and the tidal constituents are shown in the right panels. The subtidal discharge is in good agreement with observations and decreases over time during the dry season. The subtidal discharge in Can Tho and My Thuan is nearly five times smaller than the hourly tidal signal. The tidal discharge amplitudes and phases of the major semi-diurnal (M₂ and its higher harmonics) and diurnal (K₁ and O₁) tides are well-developed across the model. At Chau Doc, where tidal discharge is reduced by a factor 5 compared to Can Tho, the model overestimates the M₂ component, while at Tan Chau and Vam Nao, with stronger tidal signals, this is not the case.

Table 3. R² and NSE model performance indicators for calibration (2010) and validation (2012) models.

Obs. Stations	Correlation Coefficient (R2)				Nash-Sutcliffe Efficiency (NSE)
	Water Level		Discharge		Water Level		Discharge
	2010	2012	2010	2012	2010	2012	2010	2012
Ba Nho	0.87	-	0.91	-	0.86	-	0.83	-
Vam Kenh	0.97	0.97	-	-	0.97	0.97	-	-
Tan An	0.82	-	-	-	0.82	-	-	-
Chau Doc	0.77	0.78	0.9	0.88	0.66	0.65	0.35	0.4
Can Tho	0.96	0.96	0.92	0.93	0.96	0.95	0.92	0.92
My Hoa	0.97	-	-	-	0.96	-	-	-
An Thuan	0.97	0.98	-	-	0.97	0.97	-	-
My Thuan	0.92	0.93	0.92	0.93	0.92	0.93	0.92	0.92
Tan Chau	0.92	-	0.9	0.91	0.88	-	0.9	0.9
Ben Trai	0.97	0.81	-	-	0.97	0.8	-	-
My Tho	0.91	-	-	-	0.91	-	-	-
Vam Nao	0.96	-	0.88	0.89	0.64	-	0.87	0.88

shows R² and NSE of the model versus observations. In almost all stations, for calibration (2010) and validation (2012), both parameters are well above 0.85 (good or very good [49]), indicating a reliable model performance for vertical and horizontal tidal and subtidal barotropic tides. Only at Chau Doc with weak tidal signal, because of a mismatch in average water level, NSE is low while R² is still acceptable. However, due to limited contribution of Chau Doc to the overall model performance, we expect negligible impact from this discrepancy on addressed physical processes.

3. Results

3.1. Tidal Propagation and the Role of the Irrigation/Navigation System

Figure 7 shows spatial distribution of M₂ tidal amplitude and the co-phase lines over the main rivers and along the eastern coast of the VMD. Tidal wave travels and amplifies southward and propagates through the estuarine system with a speed of approximately 35–40 km/hour. Co-phase lines show differences in the order of a few minutes between the downstream distributaries (diminishing towards the junctions). Nevertheless, the general tidal wave propagation is similar through the Tien and Hau branches, classifying the system as hyposynchronous [50] and dominated by friction.

Figure 8a shows tidal amplitude along the Hau River and its two distributaries, starting at nearly 0.9 m (TD) or 0.8 m (DA) at estuary mouth decreasing to 0.45 m at 150 km upstream. Only near the junction, where tidal amplitude between the two branches converge, there is some tidal amplification in the DA. The Tien system, due to larger number of channels, shows a more complicated pattern. Through Co Chien–Cung Hau, the wave has a shorter distance to reach the upstream junction compared to Ham Luong–Tien–Dai through Chol Lach; therefore, the amplitude and phase lines do not meet on the same distance (x-axis) but meet on the same amplitude/phase (y-axis). In Co Chien–Cung Hau channels, the bifurcation and confluence, approximately 30 km upstream, forces an abrupt tidal amplitude attenuation at Cung Hau. This is mainly because of the interaction between the two channels, and the rapid change in Cung Hau geometry. Tidal amplitude gradually reduces upstream through Ham Luong–Cho Lach. Tieu–Dai system experiences a stronger decay of tidal amplitude compared to Ham Luong, but that leads to amplification through My Tho channel in order to match the water levels to Ham Luong–Cho Lach.

Following the approach explained in Section 2.1, Figure 9 shows that in immediate response to removal of the primary and secondary channels, M₂ tidal water level amplitude along various channels (Equation (5)) can increase up to 25% (left panel). The difference grows to 150 km upstream, where the tidal signal diminishes. Nevertheless, tidal prism variation (Equation (6)) remains within 5% of the reference simulation (right panel).

The re-calibrated model for water levels generally did not deviate more than 5% from the reference simulation except after 180 km from the sea, where the absolute tidal signal was too small to analyze relative tidal amplitude accuracy (see the dashed lines in Figure 9 left panel). However, in the simplified domain, although overall volume of primary and secondary channels (6 Gm³) is hardly 1% of the total water bodies (643 Gm³), tidal prism is then underestimated by approximately 10%. Hence, without the primary and secondary channels, it is nearly impossible to have an accurate simultaneous representation of both water level and discharge in the model.

3.2. Cumulative Freshwater Distribution

Cumulative freshwater distribution in the multi-channel estuarine system was calculated and compared with the observed data (see Figure 10). Tan Chau observation station shows that 84% of total Mekong freshwater enters through the Tien River, while only 16% enters through Chau Doc, but approximately 32.1% to 36.5% of the total freshwater inflow diverts from the Tien River to the Hau River through Vam Nao channel. Water demand is often overlooked in modeling practices within the VMD, but the data show a significant consumption, varying in time and space (as high as 2000 m³/s in April and November). Between the Cambodian border and the lower Mekong Delta (e.g., Can Tho and My Thuan), 15% of the freshwater inflow is consumed through sink terms of water demand and evaporation, while 40% of the total freshwater entering the Mekong Delta drains through Can Tho (Hau River) and 44% through My Thuan (Tien River). In total, 27.4% of the freshwater inflow of the VMD, does not reach the estuary mouth. Compared to the previous discharge division estimates, My Tho system shows the largest difference with Nguyen et al. [25], who estimated a total of 31.7% of the total discharge passing through the My Tho channel (nearly 5 times larger). Fundamental sources of difference could be in geometry of their reference hydraulic models (from 1998) versus the more recent geometry of 2008 used in this study (see Reference [51]) and inclusion of water demand and evaporation.

3.3. Temporal Freshwater Distribution

Tidal and non-tidal (e.g., wind surge and river discharge) physical processes govern subtidal water level variations [27,28,29]. Figure 11a–f demonstrates the correlation between subtidal water level variations and temporal variation of absolute and relative discharge in lower VMD distributaries. Figure 11a,b show observed hourly water levels in Can Tho and My Thuan, respectively, as well as their subtidal water levels. On the right axes, time series of daily-averaged offshore wind speed and direction (pointing towards) are plotted in grey. Wind speed has relatively large fluctuations, but wind direction predominantly blows from E/NE towards W-SW. A strong correlation can be detected between significant wind events (e.g., last week of January, third week of February, second, third, and fourth weeks of March) and subtidal water level peaks, even at stations nearly 80 km inland, highlighting the importance of wind-generated ocean surge.

Figure 11c,d show the modeled subtidal discharge (solid lines, positive downstream) at different tributaries of Hau and Tien systems, respectively; line colors are consistent with branch colors in Figure 10. As an indication of spring and neap cycles, the daily variation of semi-diurnal tidal amplitude in Can Tho and My Thuan are also plotted in Figure 11c,d. Discharge in Can Tho and My Thuan descends from nearly 3500 m³/s to below 1000 m³/s in late March and April, with large fluctuations around the overall descending trends. In the same panels, the spring–neap oscillation is most pronounced in upstream Stokes transport (dashed lines), demonstrating a maximum during spring tide. In shallower channels towards estuary mouth, due to higher tidal velocities and larger channel width, Stokes transport is larger in the lower tributary channels compared to further upstream. While subtidal discharge decreases over time, the proportion of Stokes transport relative to the overall subtidal discharge increases in time.

Figure 11e,f highlight the relevance of surge events to freshwater division by showing the relative discharge (relative to their immediate upstream main discharge station) of lower distributaries for Hau and Tien systems, respectively This shows strong modulations especially during the second and fourth weeks of March. As an example, in Hau system, downstream of Can Tho, discharge division is on average 60% versus 27% in favor of DA (the wider and deeper channel). During the second week of March, this balance collapses when DA is partially importing water from the ocean, while TD is still draining with lower than average discharge. During the fourth week of March, TD channel is importing nearly as much water into the estuary as is being discharged downstream in Can Tho. This results in storage of water in the estuary contributing to upstream subtidal wave propagation.

Comparison of absolute (Figure 11c,d) and relative (Figure 11e,f) discharge with surge effects on subtidal water level (Figure 11a,b) indicates that the system seems to first react with remarkable drops in total discharge (e.g., in Can Tho and My Thuan) followed by dramatic alteration in subtidal discharge division in response to the surge events. However, although the spring–neap cycle can be observed in Stokes transport, the exact effect on discharge division cannot be directly understood. Whether the Stokes return flow remains in the same channel or diverts at the junction in a neighboring channel (during ebb) is a function of geometry, tidal stage, and temporal roughness variations and is difficult to identify; there can be significant exchange between channels in terms of Stokes return flow [29,34], which is further addressed in Section 3.4.

3.4. Physical Processes and Sensitivity Analysis

To discriminate the effects of river discharge, tide, ocean surge level, and their interactions, a set of six sensitivity analysis simulations were run (following References [29,34,52]).

Table 4. List of simulations decomposing the impact of various forces on discharge division.

Sim. Name	Physics Included	Outcome
Reference	R + T + S	N.A.
No Surge	R + T	Surge drives large discharge variations.
No River Discharge	T + S	River discharge moderates the effect of event-like wind-driven ocean surge levels
No Tide	R + S	Tides, along with surge, influence subtidal discharge and its division, but less significant compared to surge
Only Discharge	R	In absence of downstream processes (e.g., tide and surge), flow division, independent of river discharge, only depends on distributary channel properties
Only Surge	S	In absence of varying upstream discharge and downstream tidal forcing, discharge division depends on channel geometry (similar to Only Discharge)
Only Tide	T	development of a fortnightly tidal wave due to fortnightly oscillation of frictional forces

shows simulations carried out to decompose various contributions to subtidal flow division. The ‘Physics’ column of the table indicates what processes were included per simulation, with R, T, and S representing river discharge, and tide and ocean surge. The ‘Outcome’ column summarizes the conclusion drawn per simulation. Figure 12 compares discharge division in the Hau system, per sensitivity analysis simulation (dashed lines) to the reference simulation (solid lines).

The No Surge and No Discharge simulations show that surge is responsible for large discharge variations, but event-like variations can be moderated by river discharge. Hence, the event of the last week of January in the No Discharge simulation is not significant in the reference simulation. The No Tide simulation exaggerates the importance of surge on peak flow division modulations, but also shows that tides, along with surge, influence subtidal discharge and its division, although less significant than surge. The Only Discharge simulation shows that river discharge defines the overall descending trend of subtidal discharge. It also shows that in absence of downstream processes (tide and surge), flow division is independent of upstream discharge magnitude [53] and mainly depends on downstream distributary channel properties [30,33]. The Only Surge simulation shows that during peak events, water is diverted through the channel with larger cross-sectional area. Similar to Discharge Only simulation, DA seems to be importing and exporting approximately 50% more discharge than TD. The Only Tide simulation shows a clear fortnightly oscillation of discharge and its division, with a zero-sum when integrating over longer period; this highlights development of a fortnightly tidal wave due to fortnightly oscillation of frictional forces [54].

To accentuate the importance of ocean surge and tide interaction with discharge and their relevance to flow division, in two No Surge and No Tide simulations (including water demand and evaporation), we compare the relative discharge (relative to Can Tho) to the reference simulation. Figure 13 shows that by excluding surge level or tide, during the low discharge regime, we may make errors of generally in the order of 15%–20% and occasionally up to 100% in temporal discharge division estimation. Note that during the surge events, due to storage of water in the system, the sum of relative discharge in two channels does not add up to 100% of the upstream discharge.

4. Discussion

To demonstrate the importance of various physical processes and properties, a 1D–2DH barotropic numerical domain was calibrated against extensive observation data during the dry season. Particular attention was given to the low flow regime, since due to lack of upstream freshwater inflow, (1) tidal and subtidal dynamics become significant, (2) relative importance of water demand peaks, and (3) water quality, freshwater, and salt intrusion become imperative to the VMD’s livelihood. Spatial and temporal variability of tidal amplitude, and the interaction of upstream discharge and downstream subtidal water level emphasizes the importance of a coastal model domain to capture the interactive nature of the multi-channel estuarine system with the coastal sea.

4.1. Effect of Primary and Secondary Channels

The Mekong estuarine system can be classified as Mesotidal Hyposynchronus [50], where friction dominates convergence and tides decay upstream. The importance of the system of primary and secondary channels on tidal propagation in the VMD estuarine system was quantified. The inclusion of lateral channels in an estuary is analogous to acoustic waves travelling in a pipe with side branches [55,56,57]. Studying the effect of a random lateral channel in a semi-enclosed basin has concluded that depending on the geometry and location, a lateral channel can result in damping or enhancement of incoming tide in the main channel [57,58]. However, to the best of our knowledge, so far, no study has addressed the effect of numerous small channels on the overall tidal behavior along a multi-channel estuarine system. The role of the network of primary and secondary channels in the VMD is shown to be closest yet not identical to an additional friction in the system. When eliminated from the model domain, M₂ tidal amplitude can increase up to approximately 25%. Despite their negligible volume (1% of the total volume), these channels still contribute to nearly 10% of the tidal prism of the VMD. In a modeling exercise, where the model was calibrated for water level, and the domain excluded the primary and secondary channels, the artificial suppression of water level through friction reduced the tidal prism by approximately 10%. Similar systems can be seen in other highly populated deltas such as Yangtze and Pearl River where agriculture dominates regional productivity. It is perhaps good practice to evaluate the relevance of these irrigation systems in other deltas in large scale models. Note that during the wet season, due to the smaller tidal intrusion length, the relative effect of the irrigation system is smaller.

4.2. Temporal and Cumulative and Temporal Discharge Division

Offshore subtidal water level of the VMD, in response to variations in north-easterly monsoon, is subject to modulations in the order of 50–60 cm, with time scales in the order of 2–3 days. By introducing a phase lag to measured subtidal water levels at the estuary mouth, and superposing it to the offshore boundary, we were able to simulate the subtidal water level dynamics throughout the model domain. This implies that subtidal water level is not dominated by nonlinear shallow water effects but is initiated further offshore. Furthermore, since subtidal water level variations along different estuaries show a strong correlation with offshore wind speed, it can thus be concluded that the main source of temporal variations in subtidal water level is wind-induced surge.

We have shown that subtidal water level variations can influence and occasionally reverse temporal discharge division within the multi-channel estuarine system, especially during the low flow regime (e.g., the events of the second and fourth weeks of March 2010). These modulations are sensitive to upstream discharge, and higher discharge could moderate their effect. Therefore, the chance of freshwater division deviating radically from its average increases towards the end of the dry season. In addition, tides, in timescales of spring–neap variations, introduce an alternating import–export (respectively) cycle that can magnify or damp the surge-induced net discharge. Furthermore, during the low-flow regime, exclusion of offshore subtidal water level variations can result in relatively large errors in temporal discharge division estimates. Note that salt intrusion is highly sensitive to discharge, and it is expected that these variations in discharge division have a direct impact on salt intrusion in different channels. During the wet season, discharge division is more dependent on morphology and geometry of the distributary system and independent of upstream traveling surges.

Despite the effect of subtidal water level on temporal discharge division, it has limited impact on cumulative (integrated over a longer period) discharge division. When integrated over a year, the effect of water demand dwindles in the wet season. Freshwater proportions of Hau and Tien Rivers increase from 40% and 44% in the dry season to approximately 48% and 52%, respectively [59].

Table 5. Discharge division estimates in some of the lower distributaries from different studies.

Study	Dinh An	Tran De	Tieu	Dai
Nguyen et al. [25]	70% of Can Tho, 27.4% of total VMD freshwater inflow	30% of Can Tho, 11.5% of total VMD freshwater inflow	10.7% of total VMD freshwater inflow	20.8% of total VMD freshwater inflow
Nowacki et al. [19]	40% of the sum of DA and TD	60% of the sum of DA and TD	N.A.	N.A.
Xing et al. [23]	94% of Can Tho	6% of Can Tho	N.A.	N.A.
Present study	60% of Can Tho, 23.8% of total VMD freshwater inflow	28% of Can Tho, 10.8% of total VMD freshwater inflow	2.1% of total VMD freshwater inflow	3.1% of total VMD freshwater inflow

highlights the major discrepancies in discharge divisions of downstream distributaries between this study and the existing literature. To estimate discharge per channel, Nguyen et al. [25] calibrated Savenije’s analytical salt intrusion model [26], and compared it to estimates obtained with numerical models with geometries from the mid-1990s. This may partly explain the differences compared to this study. The fact that lower Hau estuarine system is morphologically stable [60], may explain why we observe less difference in Hau and larger differences in the Tien tributaries.

The study of Xing et al. [23] calculated discharge division in March 2015, using a 2DH model of the lower Song Hau system (downstream of Can Tho), and arrived at a radical freshwater division with TD receiving only 6% of Can Tho discharge. This, in practice, implies a highly saline TD channel, which is not aligned with previous reports/measurements [17,18,19]. Nowacki et al. [19] measured discharge in DA (25–26 April 2013) and upstream of the junction (27–28 April 2013) concluding that DA only transports 40% of the sum of DA and TD. From synthetic studies (e.g., Reference [34]), it is expected that the deeper and wider channel (DA) carries larger discharge, while Nowacki et al. [19] arrived at a contrasting division. We have shown that considering large temporal discharge variations in the Hau system, without simultaneous discharge measurement in TD over a significantly longer period, estimating freshwater division from one-day measurements of different days at different channels is prone to significant errors.

4.3. Limitations and Uncertainties

Geometry is perhaps the most vital input in modeling the VMD and is likely to have changed in the past 10 years, except the Hau River, downstream of Can Tho, where the most recent measurements were used and the overall morphological changes are negligible [60]. Note that due to the quality of observed data in Cambodia and the old morphology of the Cambodian channels, model evaluation in the Cambodian territories becomes a qualitative exercise. However, this is still the best and most recent available bathymetric data. Although this could be source of an error we cannot quantify, we see a good model performance within the VMD.

The routinely carried out water level measurements are most likely accurate in terms of variation with respect to the mean (not always accurate in referencing), but available discharge observations may carry certain limitations. Hourly discharge in five major stations of the VMD are translation of water level and near-bank point velocity measurements to discharge through rating curves. The rating curves are updated per quarter, following 10–15-day-long transect measurements (in communication with the Southern Regional Hydro-Meteorological Centre), which is prone to errors and uncertainties (e.g., see Reference [35]). Therefore, in assessing discharge, we did not solely refer to the measurements or the simulations, but we treated them as best available indications. Lastly, discharge division is sensitive to spatial distribution of water demand. It is always a delicate matter how the substantial amounts of water demand are introduced as sources/sinks to the numerical model, which perhaps also varies in years in response to immediate local demand. Moreover, although it is expected that during the dry season most of the extracted water from the surface water system is evaporated within the agriculture and aquaculture fields, the exact fate of this water is not clear.

5. Conclusions

The multi-channel estuarine system of the VMD, with the countless number of irrigation channels, constitutes a complicated aquatic network of natural rivers and artificial channels. Using a coupled 1D–2DH model, the irrigation network was shown to play an important role in damping tidal amplitude along the estuary (up to 25%). Despite their limited volume (1% of the total volume), the presence of the irrigation system contributes to up to 10% of tidal prism. In absence of the irrigation system, accurate model calibration for both water level and discharge cannot be obtained; i.e., if the model is calibrated for water levels, tidal discharge will be underestimated. Therefore, although it is perhaps best practice to include the system of primary and secondary channels in any modeling effort of the VMD, our research clarifies on what compromises are made in case of neglecting them in the model domain.

The model results showed that the Mesotidal estuarine network is hyposynchronous, and over the majority of the system, friction dominates convergence and tidal amplitude gradually declines upstream. An updated view of cumulative discharge division in the Mekong Delta is presented. It was shown that, in absence of ocean forces (tide and surge), regardless of magnitude, freshwater division remains relatively uniform in time. However, surge level in the ocean can impose significant temporal alterations to freshwater distribution during low flow regimes. While surge level in the ocean is event-like, tide influences freshwater division by creating subtidal waves at the spring–neap time scale. When discharge is large, the subtidal flow division in the delta is relatively constant in time, meaning that the role of subtidal water level (ocean surge and spring–neap cycle) increases by reduced discharge towards end of the dry season, when salt intrusion and freshwater division become critically important and responsive to freshwater inflow. The interactive nature of the above processes shows the importance of modeling marine and fluvial processes within the same domain.