Essence of MIKE 21C (FDM Numerical Scheme): Application on the River Morphology of Bangladesh

Shiblu Sarker

Bureau of Watershed Management and Modeling, St. Johns River Water Management District, Florida, USA.

DOI: 10.4236/ojmsi.2022.102006                

Abstract

Understanding, quantifying, and forecasting water flow and its behavior in environment is made possible by the use of computational hydraulics in conjunction with numerical models, which is one of the most powerful tools currently available. It is made up of simple to complex mathematical equations having linear and/or nonlinear elements, as well as ordinary and partial differential equations, and it is used to solve problems in many areas. In the vast majority of cases, it is not useful to reach analytical solutions to these mathematical equations using conventional methods. In these settings, mathematical models are solved by employing a variety of numerical algorithms and associated schemes. As a result, in this manuscript, we will cover the most fundamental numerical approach, the Finite Difference Method (FDM), in order to reformulate the governing equations for water and sediment flow from a system of partial differential equations to a system of linear equations. As part of our analysis into the inner workings of a computer program known as MIKE 21C, we will attempt to gain a better understanding of the hydrodynamic processes that take place in major rivers in Bangladesh. In addition to that, we will go over some of the most commonly used morphological studies that have been conducted on Bangladesh’s major rivers, including morphological solutions that have been developed in response to water supply concerns.

Keywords

Curvilinear Model, Finite Difference Method, MIKE 21C

1. Introduction: An Overview on Finite Difference Method

Finite-difference methods (FDM) are a class of numerical frameworks for solving differential equations using finite differences to approximate the derivatives. This technique’s popularity stems from its ease of formulation and application. It can be used to solve any governing equation using numerical analysis, which generates a collection of algebraic equations by substituting finite differences for partial derivatives [1]. Thus, these algebraic equations can be solved explicitly or implicitly, where FDM employs a structured grid and thus requires grid transformation in the case of multidimensionality [1]. A typical finite difference grid is depicted in Figure 1(a). Numerous numerical schemes are established on the basis of various discretization techniques [1], some of which are presented as Equations (1)-(9). Any partial time derivative of a function C can be written as Equation (1) in a forward time difference (1), forward space difference (2) and backward space difference (3).

$\frac{\partial C}{\partial t}=\frac{C_j^{n+1}-C_j^n}{\partial t}$ (1)

$\frac{\partial C}{\partial x}=\frac{C_{j+1}^n-C_j^n}{\partial x}$ (2)

$\frac{\partial C}{\partial x}=\frac{C_j^n-C_{j-1}^n}{\partial x}$ (3)

Besides that, any partial space derivative of a function C in a centered space difference can be represented as Equation (4)

$\frac{\partial C}{\partial x}=\frac{C_{j+1}^n-C_{j-1}^n}{2\partial x}$ (4)

where, the temporal and spatial discretizations are $\Delta t$ and $\Delta x$, respectively. When space discretizations are performed on the present time level, the scheme

Figure 1. (a) Schematic 1D space and time discretization (taken from [1] ), (b) Flow parameter P and Q, flow depth H in a Curvilinear coordinate system (MIKE 21C).

is referred to as an explicit scheme. For instance, Equations (2) to (4) illustrate an explicit space discretized scheme. When spatial discretizations are performed on the next time level or using both the current and next time levels, the resulting scheme is referred to as an implicit scheme [1]. For instance, if the following Equation (3) can be discretized:

$\frac{\partial C}{\partial x}=\frac{C_j^{n+1}-C_{j-1}^{n+1}}{\partial x}$ (5)

or

$\frac{\partial C}{\partial x}=\frac{1}{2}\left(\frac{C_j^n-C_{j-1}^n}{\partial x}+\frac{C_j^{n+1}-C_{j-1}^{n+1}}{\partial x}\right)$ (6)

Afterwards, both will produce implicit schemes. Due to the fact that an explicit scheme generates no systems of equations, it can be solved directly [1]. This means that unknowns at any nodal value can be estimated independently of the solution of other nodal values. On the other hand, an implicit scheme generates a system of equations, implying that unknowns at any nodal value cannot be solved independently [1]. Additionally, another popular implicit scheme is the Box Finite Difference Scheme, also referred to as the four-point implicit scheme [1]. Using this scheme, the time derivative can be expressed as Equation (7), space derivative can be expressed as Equation (8), and C can be expressed as Equation (9) (depicted in Figure 1(a)).

$\frac{\partial C}{\partial t}=\frac{1}{2}\left(\frac{C_{j+1}^{n+1}-C_{j+1}^n}{\partial t}+\frac{C_j^{n+1}-C_j^n}{\partial t}\right)$ (7)

$\frac{\partial C}{\partial x}=\frac{1}{2}\left(\frac{C_{j+1}^{n+1}-C_j^{n+1}}{\partial x}+\frac{C_{j+1}^n-C_j^n}{\partial x}\right)$ (8)

$C=\frac{1}{4}\left(C_j^n+C_{j+1}^n+C_j^{n+1}+C_{j+1}^{n+1}\right)$ (9)

Additional schemes for solving governing equations, such as the Crank-Nickelson Scheme, the Box Finite Difference Scheme, the Lax-Wendroff Scheme, and the MacCormack Scheme, have been proposed to improve accuracy [1]. Prior to solving the governing equations, the initial and boundary conditions must be specified. The most frequently encountered boundary conditions are Dirichlet, Neumann, and Cauchy. Dirichlet condition exists when a portion of the boundary is at the specified C level, Neumann condition exists when a portion of the boundary has specified C and is crossing the boundary curve normally (i.e. $\frac{\partial C}{\partial n}$ ), and Cauchy condition exists when a portion of the boundary has specified C and is crossing the boundary curve normally [1]. Numerical method for approximate reformulation of governing equations in which continuous operations are replaced by discrete operations, resulting in truncation error [1]. On the other hand, the discretized solution domain’s size results in round-off error (see Figure 1(b) [1] ). Model discretization and numerical solution ideally consistent (i.e. correct), convergent (i.e. $\partial x\to 0\to \epsilon \to 0$ ), and stable (i.e. error should not increase) [1].

2. Methods: MIKE 21C Curvilinear Model

MIKE 21C uses curvilinear finite difference grids (Curve FDG) to predict river bed and channel planform development in two dimensions. Curve FDG provides major advantages in simulating the scenario of construction, dredging, seasonal flow fluctuations, and other events that contribute to bank erosion, scouring, and shoaling [2]. This section aims to describe and review mathematical and numerical aspects in plain language [2].

2.1. Curvilinear Grid Generator

This section introduces the Curvilinear Grid (CG) Generator using a basic mathematical and numerical framework. These include discretisation, solution, initial conditions, smoothing methods, and residual evaluation. MIKE 21C’s grid generator composed two primary PDE [2].

$\frac{\partial }{\partial \zeta }\left(w\frac{\partial x}{\partial \zeta }\right)+\frac{\partial }{\partial \eta }\left(w\frac{\partial x}{\partial \eta }\right)=0$ (10)

$\frac{\partial }{\partial \zeta }\left(w\frac{\partial y}{\partial \zeta }\right)+\frac{\partial }{\partial \eta }\left(w\frac{\partial y}{\partial \eta }\right)=0$ (11)

where x and y are CarCS coordinates, and w denotes a weight factor defined by Equation (12)

$w=\sqrt{\frac{x_{\zeta }^2+y_{\zeta }^2}{x_{\eta }^2+y_{\eta }^2}}$ (12)

The condition of non-linear orthogonality is the boundary condition for this system, namely $x_{\zeta }{x}_{\eta }+x_{\zeta }{x}_{\eta }=0$. Additionally, $f\left(x\mathrm{,}y\right)=0$ denotes that x and y are points on a particular curve. The system described previously generates an orthogonal grid with streamlines ( $\zeta$ lines) and potential lines ( $\eta$ line). This boundary condition is equivalent to the kinematic boundary condition, in which streamlines are parallel to a boundary. When the grid weight of the system is equal to one ( $w=1$ ), the system is said to be conformal; consequently, a river is the ideal example of a conformal mapping constraint, that can be specified as $\frac{N}{M}=\frac{L_x}{L_y}$, where N and M denote the number of grid points in each direction, respectively, and $L_x$ and $L_y$ denote the geometry’s extent. This conformal mapping formulation is invariant with regard to changes in coordinate systems, which helps the overall aspect ratio of the complex geometry governing the weight function. To solve the Equations (10) and (11) Newton-Raphson method for the boundary conditions along with Stone’s implicit elliptic solution method are adopted [2]. Hence, the discretization of the grid can be stated as

$\begin{array}{l}{x}_p=C_N{x}_N+C_E{x}_E+C_S{x}_S+C_W{x}_W\\ y_p=C_N{y}_N+C_E{y}_E+C_S{y}_S+C_W{y}_W\end{array}$ (13)

where, N, E, S and W are indicates four directional grid around $x_p$ and $y_p$. Implicit conditions can be used to calculate the coefficients in 13 equations [2]. The boundary conditions are regarded as evolving Dirichlet conditions throughout the solution process. The implicit process is also utilized in iteratively solving equations using elliptic solution methods, which are strongly related to the notion of fast short wave error decay [2]. The initial conditions were constructed using a widely utilized transfinite interpolation approach and an algebraic grid generator. The algebraic mapping methodology is substantially faster than elliptic approaches; however, the grids formed by elliptic grid generators are often not as smooth, which is why the grid weights are smoothed prior to solving the elliptic system. Grid weights have been placed at cell corners, where the x and y coordinates can also be obtained, and the grid weight function for the current single block version of the grid generator has been generated using central differences in the interior and one-sided differences at the boundaries [2]. The adaptive filter determines the smoothness of boundary lines based on their relative curvature; if a specified criterion is not met, the boundary point is subjected to the running average; this procedure is repeated until all boundary points have a relative curvature less than the limit; the corner points are fixed during the smoothing process [2]. The CG Generator uses splines to represent the boundary lines; the equation and orthogonality condition are non-linear and are solved using a modified Newton-Raphson technique, and the iteration is ended when a given time step requirement is fulfilled [2]. The residual field is as Equation (14), where L is the scale of linear cell length on the geometry, i.e.

$L=\frac{\sqrt{Area}}{Cel{l}_{No}}$.

$r=\frac{1}{L}\left(\begin{array}{c}{x}_p^{new}-x_p^{old}\\ y_p^{new}-y_p^{old}\end{array}\right)$ (14)

2.2. 2D Hydrodynamic Model

Despite the intricacy of rivers’ three-dimensional flow patterns, conducting long-term simulations to examine river morphology would require unfeasible computational efforts. As a result, the Navier-Stokes equations can be simplified and reduced to two-dimensional mass and momentum conservation equations in two horizontal directions, while the depth-averaged model can retain secondary flow by introducing a separate helical flow component and assuming vertical distribution of flow velocities [2]. Two parallel horizontal axes can be coupled to provide a curvilinear orthogonal coordinate system, with one axis following the river’s bank lines, which provides a more detailed description of the flow field near the boundary, which is particularly important when computing bank erosion [2]. The Equations (15)-(18) define the transformation from Cartesian coordinate systems (CarCS) to curvilinear coordinate systems (CurCS) (depicted in Figure 1(b)).

$h=H\text{and}\frac{\partial h}{\partial x}=\frac{\partial H}{\partial s}\text{and}\frac{\partial h}{\partial y}=\frac{\partial H}{\partial n}$ (15)

$u=U\text{and}v=V$ (16)

$\frac{\partial u}{\partial x}=\frac{\partial U}{\partial s}-\frac{V}{R_s}\text{and}\frac{\partial u}{\partial y}=\frac{\partial U}{\partial n}-\frac{V}{R_n}$ (17)

$\frac{\partial v}{\partial x}=\frac{\partial V}{\partial s}+\frac{U}{R_s}\text{and}\frac{\partial v}{\partial y}=\frac{\partial V}{\partial n}+\frac{U}{R_n}$ (18)

where, s and n axes are two horizontal axes intersecting at right angle, h and H are depth in CarCS and CurCS respectively, similarly, $\left(u\mathrm{,}v\right)$ and $\left(U\mathrm{,}V\right)$ are velocity components in the CarCS and CurCS respectively, $R_s$ and $R_n$ denote radius of curvature of the s-lines and the n-lines, respectively. Under three major approximations, the hydrodynamic model solves the vertically integrated Saint Venant equations in two directions [2]. Approximation of shallow water, hydrostatic pressure distribution, and rigid lid approximation. Lateral momentum exchange due to fluid friction, as well as wall effects along riverbanks, are ignored in the shallow water approximation. Under the hydrostatic pressure distribution assumption, the gradients of the vertical velocity component are ignored. For the water surface condition, the rigid lid approximation [2] proposes that the water surface is a rigid, impermeable, shear-free plate that encounters only normal stresses. In other words, the flow model applies to shallow, gently moving topography as well as mildly curved and vast river channels with low Froude numbers. Thus, river equations incorporate flow acceleration, convection, and cross-momentum, as well as water surface slopes, bed shear stress, momentum dispersion, Coriolis forces, wind forces, flow curvature, and helical flow. As a result, the governing equations for the curvilinear hydrodynamic model can be represented as Equations (19)-(21):

$\frac{\partial H}{\partial t}+\frac{\partial p}{\partial s}+\frac{\partial q}{\partial n}-\frac{q}{R_s}+\frac{p}{R_n}$ (19)

$\frac{\partial p}{\partial t}+\frac{\partial }{\partial s}\left(\frac{p^2}{h}\right)+\frac{\partial }{\partial n}\left(\frac{pq}{h}\right)-2\frac{pq}{h{R}_n}+\frac{p^2-q^2}{h{R}_s}+gh\frac{\partial H}{\partial s}+\frac{g}{C^2}\frac{p\sqrt{p^2+q^2}}{h^2}=RS$ (20)

$\frac{\partial q}{\partial t}+\frac{\partial }{\partial s}\left(\frac{pq}{h}\right)+\frac{\partial }{\partial n}\left(\frac{q^2}{h}\right)+2\frac{pq}{h{R}_s}-\frac{q^2-p^2}{h{R}_n}+gh\frac{\partial H}{\partial n}+\frac{g}{C^2}\frac{q\sqrt{p^2+q^2}}{h^2}=RS$ (21)

where, C is Chezy roughness coefficient, p, q Mass fluxes in the s- and n-direction, respectively and RS is the force balance terms such as Reynolds stresses, Coriolis force and atmospheric pressure. Reynolds stresses can be described in a CG for the p and q-direction as:

$\begin{array}{l}\frac{\partial }{\partial x}\left(E\frac{\partial P}{\partial x}\right)+\frac{\partial }{\partial y}\left(E\frac{\partial P}{\partial y}\right)\\ =\frac{\partial }{\partial s}\left(E\frac{\partial p}{\partial s}\right)+\frac{\partial }{\partial n}\left(E\frac{\partial p}{\partial n}\right)-\frac{2E}{R_s}\frac{\partial q}{\partial s}-\frac{\partial E}{\partial s}\frac{q}{R_s}-\frac{2E}{R_n}\frac{\partial q}{\partial n}-\frac{\partial E}{\partial n}\frac{q}{R_n}\end{array}$ (22)

$\begin{array}{l}\frac{\partial }{\partial x}\left(E\frac{\partial Q}{\partial x}\right)+\frac{\partial }{\partial y}\left(E\frac{\partial Q}{\partial y}\right)\\ =\frac{\partial }{\partial s}\left(E\frac{\partial q}{\partial s}\right)+\frac{\partial }{\partial n}\left(E\frac{\partial q}{\partial n}\right)+\frac{2E}{R_s}\frac{\partial p}{\partial s}+\frac{\partial E}{\partial s}\frac{p}{R_s}+\frac{2E}{R_n}\frac{\partial p}{\partial n}+\frac{\partial E}{\partial n}\frac{p}{R_n}\end{array}$ (23)

where, $\left(P\mathrm{,}Q\right)$ are fluxes described in a CarCS and $\left(p\mathrm{,}q\right)$ are fluxes described in CurCS [2]. The curvature of the coordinate lines introduces additional terms into the flow PDE, which is implicitly solved on a computational grid consisting of variables P and Q in two horizontal directions and a water depth H, with grid coordinates referring to the flow variables defined in all grid locations. Due to the curved grid, the space steps also vary throughout the domain [2]. When water approaches a river bend, it causes a centrifugal force imbalance, which results in outward motion near the free surface and inward motion near the bed [2]. Because the main stream velocities are greater in the upper part of the flow than in the lower part, water particles in the upper part of the water column must take a more curved path than water particles in the lower part to maintain a nearly constant centripetal force throughout the depth [2]. This classical analytical solution predicts the formation of a single helical vortex that transports fluid downstream in spiral paths while simultaneously creating a lateral free surface slope to maintain equilibrium between the lateral pressure force, centripetal force, and lateral shear force generated by friction along the bed [2]. This helical flow pattern is made up of a longitudinal component and a circulation in a plane perpendicular to the principal flow direction. Helical flow intensity $i_s$ may be defined as Equation (24):

$i_s=u\frac{h}{R_s}$ (24)

where u and $R_s$ denote the longitudinal flow and curvature radius of streamlines, respectively. Due to the critical nature of the direction of bed shear stress in a curved flow field in a bed topography model for river bends, logarithmic technique (see [2] ) offers the bed shear stress direction:

$\tan {\partial }_s=-\beta \frac{h}{R_s}$ (25)

where, ${\partial }_s$ is the angle between bed shear stress and depth averaged shear stress or flow [2]. The parameter $\beta$ can be defined as $\beta =\alpha \frac{2}{{\kappa }^2}\left(1-\frac{\sqrt{g}}{\kappa C}\right)$, where, $\kappa =0.4$ is commonly known as Von Kárman’s constant and α = calibration constant [2]. In numerical testing, secondary flow profile adaptation happens substantially faster approaching the bed, complicating secondary flow adaptation. As a result, the adaption length depends on water depth and friction [2], and MIKE 21 C uses the differential length scale ( ${\lambda }_{sf}$ ).

${\lambda }_{sf}=\frac{1.2hC}{\sqrt{g}}$ (26)

As a result, under steady flow conditions with constantly varying curvature, the direction of bed shear stress can be approximated as Equation (27):

${\lambda }_{sf}\frac{\partial \left(\tan {\partial }_s\right)}{\partial S_s}+\tan {\partial }_s=-\beta \frac{h}{R_s}$ (27)

Equations (28) are used to convert Equation (27) into a fixed $\left(s\mathrm{,}n\right)$ coordinate system:

$\frac{\partial s}{\partial S_s}\frac{\partial \left(\tan {\partial }_s\right)}{\partial s}+\frac{\partial n}{\partial S_s}\frac{\partial \left(\tan {\partial }_s\right)}{\partial n}+\frac{\tan {\partial }_s+\beta \frac{h}{R_s}}{{\lambda }_{sf}}=0$ (28)

This results in the equation for advection-dispersion, which may be numerically solved using the MIKE 21C.

$u\frac{\partial \left(\tan {\partial }_s\right)}{\partial s}+v\frac{\partial \left(\tan {\partial }_s\right)}{\partial n}+\frac{\sqrt{p^2+q^2}}{h{\lambda }_{sf}}\left(\tan {\partial }_s+\beta \frac{h}{R_s}\right)=0$ (29)

where, $\frac{\partial s}{\partial S_s}=\frac{p}{\sqrt{p^2+q^2}}$ and $\frac{\partial n}{\partial S_s}=\frac{q}{\sqrt{p^2+q^2}}$.

While the hydrodynamic model is based on depth-averaged flow equations, the morphological model uses vertical velocity profiles to determine bed shear stress and suspended sediment movement. As a result, the Reynolds stress concept and the Prandtl mixing length hypothesis may be used to determine shear stresses in a fluid [1], with the assumption that viscous (laminar [3] ) friction is significantly less than turbulent friction.

2.3. Sediment Transport

Sediment movement is often classified into three modes: bed load, suspended load, and wash load. However, the latter defines bed load transport as the movement of bed material along the bed; suspended load transport as the movement of sediment suspended in the fluid for an extended period of time; suspended sediment can be classified as both bed material and wash load; and wash load transport as the movement of material finer than bed material. MIKE 21 C studies just bed material transport and the portion of suspended load originating from bed material. This is because only bed material transport is relevant for the morphological evolution of alluvial rivers as a result of the interaction of bed bathymetry and hydrodynamics. Suspended load acts substantially differently from bed load, which can be incorporated into sediment movement modeling. The model for suspended sediment transport is based on the governing Equation (30):

$\begin{array}{l}\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}+w\frac{\partial C}{\partial z}\\ ={\omega }_s\frac{\partial C}{\partial z}+\frac{\partial }{\partial x}\left(\epsilon \frac{\partial C}{\partial x}\right)+\frac{\partial }{\partial y}\left(\epsilon \frac{\partial C}{\partial y}\right)+\frac{\partial }{\partial z}\left(\epsilon \frac{\partial C}{\partial z}\right)\end{array}$ (30)

where, C = concentration of suspended sediment, $\epsilon$ = is the turbulent diffusion coefficient and ${\omega }_s$ = fall velocity of sediment particles in suspension. If all diffusion terms are ignored except the vertical diffusion term, the equation along a streamline is

$\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial s}+v\frac{\partial C}{\partial n}+w\frac{\partial C}{\partial z}={\omega }_s\frac{\partial C}{\partial z}+\frac{\partial }{\partial z}\left(\epsilon \frac{\partial C}{\partial z}\right)$ (31)

Velocities are represented in the following way

$\begin{array}{l}u\left(z\right)=U_h{p}_1\left(\eta \right)\\ v\left(z\right)=\frac{U_{h}H}{R_s}{p}_2\left(\eta \right)\end{array}$ (32)

where, $U_h$ = depth-averaged flow velocity, $\eta$ = non-dimensional vertical coordinate $\frac{z}{H}$, $p_1$ = longitudinal velocity profile and $p_2$ = transverse velocity profile. While the surface boundary condition ensures that no silt flows past the border, the bottom boundary condition might specify a bottom concentration. A special asymptotic approximation technique can be utilized to provide information about the concentration profile in circumstances when the vertical diffusion coefficient and the fall velocity term are dominant terms, or when e is very tiny, as in $e=\frac{h}{L}\frac{U}{{\omega }_s}$. For a slowly variable flow, the general solution is based on the notion that the concentration profile is the sum of ascending order profiles. Equation (31) is used to determine the concentration of various orders. The equilibrium between suspended sediment settling and vertical diffusion is examined to identify the 0^th order concentration, and this procedure is repeated for the first and second order concentrations [1]. The interplay between the bed load and the alluvial bed is a fundamental element of a river bend’s morphological behavior. In comparison to the suspended load, it is assumed that the bed load is more responsive to changes in the local hydraulic conditions, necessitating the omission of advection-dispersion models. However, the effect of a sloping river bed must be considered, as must the deviation of the bed shear stress direction from the mean flow direction caused by helical flow. As a result, modeling helical flow must be distinct from modeling bed load. To determine the capacity of a local bed load sediment transport system, it is required to examine only the uniform shear flow solutions that have been proposed for this situation over the years (1975 and 1984). Due to the fact that bed slope has an effect on both the rate and direction of sediment transport, both of which are crucial for morphological modeling, two approaches were used. The first alters the threshold for motion start due to shear stress [2].

${\theta }_c={\theta }_{C0}\left(1+\frac{\partial z_b}{\partial s}\right)$ (33)

where, ${\theta }_c$ = modified critical shields parameter, ${\theta }_{C0}$ = Critical Shields parameter in uniform shear flow and $z_b$ = bed level. If a model does not assume zero bed load transfer at a critical shear stress, the following bed load modification can be applied along the streamline.

$S_s=S_{bl}\left(1-\alpha \frac{\partial z^{\mathrm{<mo>\; *\; </mo>}}}{\partial s}\right)$ (34)

where, $S_{bl}$ = bed load as calculated from sediment transport formula and $\alpha$ = model calibration parameter that has to be specified. River engineers place a premium on the prediction of transverse depth distributions in alluvial channel bends because they are crucial in studies of river bend navigability enhancement and channel bank protection. They recommended Equation (35) as a possible solution.

$S_n=\left(\tan {\delta }_s-G{\theta }^{-\alpha }\frac{\partial z^{*}}{\partial n}\right)S_{bl}$ (35)

where, G = transverse bed slope factor (calibration coefficient) and $\tan {\delta }_s$ = bed shear direction change due to helical flow strength. In a fixed $\left(x\mathrm{,}y\right)$ coordinate system, bed slopes are calculated as

$\frac{\partial z^{*}}{\partial s}=\frac{\partial z}{\partial x}\cos \varphi +\frac{\partial z}{\partial y}\sin \varphi$ (36)

Hence, the bed slope in the transverse direction will be

$\frac{\partial z^{*}}{\partial n}=\frac{\partial z}{\partial y}\cos \varphi -\frac{\partial z}{\partial x}\sin \varphi$ (37)

where, $\varphi$ = is the angle of stream line compared to fixed $\left(x\mathrm{,}y\right)$ coordinate system. Transformation from streamline coordinates to fixed coordinates gives the following expressions

$\begin{array}{l}{S}_x=S_s\cos \varphi +S_n\sin \varphi \\ S_y=S_n\cos \varphi -S_s\sin \varphi \end{array}$ (38)

The transport capacity of sediment in uniform shear flow has been extensively studied over the years; the model used in this modeling system to calculate bed load and suspended load transport capacity (equilibrium concentration at the riverbed) are primarily Engelund and Hansen model, Van Rijn model, Engelund, Fredsϕe and Zyserman model, Meyer-Peter and Müller model, Smart and Jaeggi model, Yang’s model for sand and gravel transport model (see details in [2] [3] ).

2.4. Morphology

Hydrodynamic and sediment transport models are combined in a morphological model. The hydrodynamic flow field is constantly adjusted as the bed bathymetry varies. Historically, linked and uncoupled morphological models have been distinguished. Coupled models integrate the flow and sediment movement governing equations into a single set of equations that may be solved concurrently. The hydrodynamic equations are solved first in uncoupled models, followed by the sediment transport equations. Following that, a new bed level is determined, at which point the hydrodynamic model continues to the next time step. This model adopts the latter strategy. Additionally, additional sub-models may be integrated, such as bank erosion, bank line updates, alluvial bed resistance, bed morphologies, and graded sediment. Equation (39) may be used to compute the change in bed level after computing the sediment transport of bed material (bed load and suspended load).

$\left(1-n\right)\frac{\partial z}{\partial t}+\frac{\partial S_x}{\partial x}+\frac{\partial S_y}{\partial y}=\Delta S_e$ (39)

where, $S_e$ = lateral sediment supply from bank erosion, n = bed porosity, $S_x$ = total sediment transport in x direction and $S_y$ = total sediment transport in y direction. The difference is computed using a space-centered-time-forwarded method. The time step is constrained by the Courant criteria, which requires that the Courant number be smaller than one. The boundary condition may be provided for each border point individually or as a sum across the whole boundary. The sediment transport model can include bank erosion in the continuity Equation (40).

$E_b=-\alpha \frac{\partial z}{\partial t}+\beta \frac{S}{h}+\gamma$ (40)

where, $E_b$ = bank erosion rate, S = near bank sediment transport, h = local water depth, z = local bed level and $\alpha \mathrm{,}\beta \mathrm{,}\gamma$ calibration coefficients specified in the model. Therefore, this model also mimic the lateral processes in morphological models that are entitled as planform module.

3. Application Example of MIKE 21C

In this section, we aim to address the application example of MIKE 21C in the water supply engineering context of morphological study of two major rivers of Bangladesh.

3.1. Morphological Study of Upper Meghna River for Saidabad Water Treatment Facility

The Dhaka Water Supply and Sewerage Authority (DWASA) is responsible for distributing potable water via a pipe network [4] to Dhaka City and its surroundings. DWASA now fulfills 78% of its water needs from groundwater sources, causing a 2 - 3 m annual decline in the ground water table. As a consequence of exploring alternate drinking water sources, DWASA has decided to switch from ground water to surface water. With this study, DWASA intends to increase the third phase capacity of Saidabad water treatment plant from 4.5 × 10⁵ m³/day to 9 × 10⁵ m³/day. The main objective of this study is to assess morphological stability of Haria Intake site. The study further aims to develop two dimensional morphological models to assess stability of the site, to assess river bed scour, river bank erosion, velocity with direction of flow and siltation/char formation near intake point and hydraulic design variables for bank protection/river training works. The upper Meghna River (outlet of Meghna basin [5] [6] ) carries the combined flow of the Surma and Kusiyara Rivers, which originate in the Indian highlands to the northeast of Bangladesh. The Surma River passes through Sylhet, which is gradually sinking. Several Surma tributaries contribute sediment that probably settles in Sylhet. The Meghna joins the Old Brahmaputra near Bhairab Bazar. The Dhaleswari meets the Upper Meghna near Munshiganj on the right bank. Some Tripura hill tributaries join the Upper Meghna on the left bank. Between Bhairab Bazar and the Padma river confluence, flow is increased by left bank spill from the Jamuna and local drainage through the Old Brahmaputra, Lakhya, Arjam, and Dhaleswari rivers. Left bank tributaries like Titas and Gumti contribute.

IWM has built up a huge collection of hydrometric and sediment data from a variety of sources. Several WL and Q stations were identified as important for this analysis; Bhairab, Badyar Bazar, and Kalagachia were identified as relevant for Q and WL. Additionally, Eklaspur and Haria are determined to be significant for data on Suspended and Bed Sediment. Statistical tests such as the Spearman Rank Correlation Trend Test, the t-Test [5] with Welch modification, and the Spearman Rank Serial Correlation Test were used to determine the consistency and homogeneity of the data. Gumbel, Log-Normal, and Log-Pearson tests, as well as Chi-square and Kolmogorov-Smirnov (K-S) tests, were performed on Bhairab Bazar Station data, indicating that average event hydrographs from 2003 can be utilized as model boundary data (Q and WL). No historical discharge peak, on the other hand, corresponds to a 100-year event. The 1988 year event was found to be the largest ever recorded in historical data set, matching the top of the 50-year event. As a result, using a 1D model, a synthetic hydrograph based on the 1988 hydrograph was produced.

Echo-sounder was used to conduct extensive bathymetric and land/char surveying in order to schematize recent river bed, sand bar, and bank elevations in order to generate the present Digital Elevation Model (DEM). A 25 km length was chosen for the study, which included near circuit branching and sand bars.

MIKE 21C was used to build a two-dimensional morphological model for the 25 km river reach to simulate hydrodynamic and morphological characteristics in each computational grid that were produced using the observed bank lines based on the previously stated key criteria. Boundary data has been collected from 1D model. At upstream boundary-discharge data and at downstream boundary-water level data has been considered. Calibration of 2D hydrodynamic module under MIKE 21C model was accomplished comparing 1D model generated water level data with 2D model generated data for location within the model domain by tuning parameters like Chezy’s C value [7] [8] and eddy Viscosity.

Direction of flow and magnitude of velocity is important for the estimation of bank protection measures. More specifically, as seen in the intake location, from the velocity vectors (flow field), near bank velocity differ (0.4 - 0.6 m/s) and maximum velocity differs (1.23 - 1.44 m/s) from peak of average to high flood events respectively. Therefore, it appears that the velocity distributions are nearly uniform at intake location at the peak of high (100 yr) and average (2.33 yr) flood events. The average and high event results show no significant bed level change near the intake site. After a high event, there is some siltation in the interest region, but the bed level stays nearly unchanged. The model findings show that our project reach is nearly stable. The accumulated bank erosion at right bank around intake location has been extracted from the model for average event and high event. No erosion of bank is predicted around the intake location. Finally model generated cross-sections for intake location also suggested that intake is located at safe reach (Figures 2-7).

3.2. Investigation of the Ganges for Suitable Location of Water Intake along with the Left Bank of Ganges

The purpose of this study is to determine the likelihood of channel alteration, flow concentration, and other characteristics occurring at the planned intake locations (Godagari, Harupur under Rajshahi and Yusufpur under Sarda). In this context, a two-dimensional Ganges erosion management model was constructed that spans a 70 km stretch from 8 km upstream of the Ganges-Mohanonda

Figure 2. The study area for two MIKE 21C application examples (Collected from IWM project archive).

Figure 3. (a) Computational CG (inset show zoomed view of Grid at intake site), (b) Model bathymetry showing the bed level.

Figure 4. (a) Water level calibration of 2D model, (b) Sediment calibration results.

Figure 5. Depth average velocity vector of upper Meghna river near Haria intake during the (a) peak of average and (b) peak of high flood events.

Figure 6. Bed profile and velocity profile at intake locations on upper Meghna river at (a) peak of average and (b) peak of high flood events.

Figure 7. Bed level changes for (top panel) average event and (bottom panel) high event (bed level in m).

confluence to 10 km downstream of Sardah. To determine compatibility, a 2D model comprising the anticipated intake areas was examined. The model has been calibrated for use in simulating various alternative scenarios. Primary data from field surveys as well as data from other research in the area were used to support the modeling investigation. Using historical data, including satellite pictures, we investigated and interpreted the results of model simulations under autonomous development settings as well as other choices. Estimate water availability by observing channel movement away from the left bank at proposed intake locations and vertical and lateral scour extent at proposed intake locations.

The Bangladesh Water Development Board has collected historical data on the Ganges River discharge at Hardinge Bridge and water level data at Sardah. Maximum flow data were used for frequency analysis, the findings of which were used as boundary conditions for the model to determine the Intake structure’s sustainability; whilst the lowest flow records were used to determine the water availability at the intake locations. As shown in the Figure, the Ganges has a maximum and minimum flow of 79,059 m³/s and 139 m³/s, respectively, with occurrences in 1998 and 2011. To establish the range of Godagari’s water levels, a series of water levels were used to estimate standard high and low water levels. Additionally, SLWL was utilized to anticipate the lowest water depth for distinct flood events in order to ensure water availability during dry periods. SHWL for Godagari is 21.27 mPWD, meaning that water levels are less than this level for 95% of the year, whereas normal low water level for the same site implies that water levels are larger than 9.35 mPWD for 95% of the year. SHWL for Sardah is 18.29 mPWD, meaning that water levels stay below this level for 95% of the year, whereas standard low water levels for this location imply that water levels remain above 7.51 mPWD for 95% of the year. The SHWL and SLWL were determined using data from Godagari station from 1945 to 1982 and Hardinge Bridge from 1998 to 2014. On the other side, the SHWL and SLWL have been determined using data from the Sardah station from 1929 to 2013.

From 1976 to 2014, it is possible to see that the channel remains common at a distance of 590 m from the bank at the Intake location by superimposing a series of satellite images depicting the channel layout over time. The channel incidence map for the Godagari site (prepared from 1976 to 2014) indicates that the main channel stays within a radius of more than 90% of the time, indicating the channel’s stability, i.e. the channel’s migration tendency is quite low. Similar analyses of the other two study locations, Harupur and Yusufpur, indicate that the channel remains common at 1900 m and 300 m downstream from the bank at the Harupur and Yusufpur Intake locations, respectively. Additionally, channel prevalence rates of 75% and 88% were observed in Harupur and Yusufpur, respectively. Additionally, it is discovered that the Ganges River’s sinuosity varies according to its location. Godagari, in particular, has a higher sinuosity of 1.4 than Harupur, with Yusufpur following closely behind. This increased sinuosity could be linked to a faster rate of land erosion and accretion (Figure 8).

A two-dimensional model has been developed that constitutes computational grid, bathymetry into the grid and boundary data. Based on the digitized bankline [9] data from the image of 2014 and verified with the surveyed bankline data, the grid of the 2D model from Godagari to Sardah has been developed, shown in Figure 9. The resolution of this grid is 775 × 126. Having generated the grid, a bathy of the model domain has been prepared from the survey data of Jan-Feb 2015 (see in Figure 9).

In order to observe the bed scour and its movement, bathymetry contour for different years from some past studies (collected from IWM archive) near the three intake locations have been plotted using developed computational grid in Figures 10-12, which shows the bed contour of the Ganges for different years.

Figure 8. Banklines (left), erosion/accretion map (middle) and channel incidence map (right) of the Ganges at intake locations.

Figure 9. Developed (a) computational grid of the 2D model of the Ganges and (b) bathymetry.

Figure 10. Observed bathymetry for different years near the proposed Intake location at Godagari.

Figure 11. Observed bathymetry for different years near the proposed Intake location at Harupur.

Figure 10 shows that bed scour just in front of the proposed intake location at Godagari, may reach up to −11 mPWD whereas at immediate downstream, it is in the order of −27 mPWD (1700 md/s from the intake location), and these values correspond to the year of 2005. Also the dominating channel flows along the left bank at Godagari side. In the case of Harupur (Figure 11), the dominating Channel remains at Intake location side only during 2015 having bed level of −14 mPWD from 2002 to 2015. Whereas at Yusufpur (Figure 12), the dominating Channel remains at the intake location side only from 2002 to 2015 having bed level of deepest point −21 mPWD during 2015. In terms of bathymetric variation and main channel shifting the Godagari has been found less dynamic in

Figure 12. Observed bathymetry for different years near the proposed Intake location at Yusufpur.

comparision of Yusufpur followed by Harupur. Above bed level contours also reveal the fact of channel changing pattern in the vicinity of the proposed Intake location, i.e. the channel approaches to the proposed locations (Godagari, Harupur and Yusufpur) with different angle; however the channel persists near the proposed Intake location at Godagari and Harupur for all the years. The channel alignment, that is revealed from the bed contour, indicates that the right channel may also get prominent at Godagari site. Such phenomenon is very prominent at Godagari for the period of 2011 and 2015. Also for Yusufpur during 2011.

For better understanding of the morphological behavior and character of the river, several cross-sections of different years near proposed intake locations have been plotted in Figures 13-15. From these figures, the following observations have been made:

According to the information above, the bed level of the near bank channel at Godagari is the deepest, followed by Yusufpur and Harupur. Additionally, Harupur’s thalweg position remains the furthest (from the left bank).

Apart from historical bathymetry, the developed model was calibrated using 2014 observed water level, discharge, and sediment data. Prior to applying the model to bankfull and extreme flood events, it is necessary to calibrate the model to ensure that it exhibits satisfactory agreement between simulated and observed data. To finalize the model’s input parameters, a calibration simulation was run using January-February 2015 bathymetry. A comparison of the Ganges’ hydromorphological behavior in terms of water level, discharge, and sediment transport is shown in Figures 16(a)-(c). Calibration was performed using observed discharge and sediment transport data at Hardinge Bridge, as well as the water level at Godagari. According to these figures, the model performs quite satisfactorily. It has been determined that the trend of simulated sediment discharge

Figure 13. Cross-sections near Godagari for different years.

Figure 14. Cross-sections near Harupur for different years.

Figure 15. Cross-sections near Yusufpur for different years.

Figure 16. (a) Water level, (b) discharge and (c) sediment discharge calibration for 2014 at marked location within model domain.

follows the trend of observed data (see Figure 16(c)). Additionally, the magnitude is within the observed data band, which is statistically acceptable. As a result, the calibration parameters used in the simulation were chosen to optimize the performance of the 2D model.

The 2D model was applied to two different flow conditions in order to determine the Ganges’ response. Bankfull and extreme flood (100 yr) conditions were considered during design and assessment. According to the statistical analysis conducted in 2011 for the Ganges Barrage study Project, the bankfull discharge of the Ganges corresponds to the hydrological year 2005 (Peak WL = 15.70 mPWD and Q = 38,080 cumec), while the extreme flood event was simulated using the flow of the year 1998 (Peak WL = 19.08 mPWD and Q = 83,607 cumec). The 2D model’s boundary conditions (water level and discharge) for these two events are depicted in Figure 17(a). To consider the worst-case scenario, another three-year model simulation was conducted using a combination of an extreme flood (1998) and two consecutive low floods (1993). In this long-term scenario, the minimum water level and discharge are 5.08 mPWD and 290 cumec in March, and the maximum water level and discharge are 19.08 mPWD and 83,607 cumec in August. The boundary conditions for this long-term event (water level and discharge) are depicted in Figure 17(b). For both extreme and bankfull conditions, the simulation period is July to October, as the monsoon season has the highest morphological activity. Additionally, a three-year simulation was conducted to assess the morphological activity following a combination of flood events (one extreme and two low flood events).

After simulation of 2D model for extreme, bankfull and long term flood event, the bathymetry along with cross-section at the end of simulation has been given in (Figures 18-20). This prediction is important in understanding the probable bed movement and channel shifting if a flood of mentioned event occurs.

According to Figure 18, erosion of 9 m at the left channel in 1998 and 3 m in 2005 revealed that a large area of the char along with the mid-channel was eroded during 1998 and in the long run. The right channel was found to be

Figure 17. Water level (WL) and discharge (Q) hydrographs as boundary conditions for (a) extreme (synthetic 1998 yr) and bankfull (2005 yr) events and (b) three-year-long simulation.

Figure 18. Simulated bathymetry and cross-section near Godagari for different flood event.

Figure 19. Simulated bathymetry and cross-section near Harupur for different flood event.

Figure 20. Simulated bathymetry and cross-section near Yusufpur for different flood event.

deposited during 1998 and long event simulations, but it was eroded in 2005, while the mid-channel remained unchanged. The channel near Godagari (left channel) has been scoured; additionally, its width has been reduced, and the left channel thalwegh has been shifted toward the left bank at approximately 100, 400, and 300 meters for 1998, 2005, and long event, respectively. In the case of Harupur, Figure 19 demonstrates that a single left channel splits into two channels with char in between during both flood events and the long run, and left channel thalwegh erosion of 1 m and 4 m has been observed during 1998 and 2005 flood events, respectively. Whereas in the long run, it is deposited to a depth of approximately 1 m, the newly created parallel channel is prominently eroded to a depth of approximately 12 m. This indicates that in the event of a flood of this magnitude, the river may divert from its current path and flow through the parallel channel approximately 2 km from the bank. For 1998, the dominant channel was discovered to be 1.5 km from the left bank. Similar phenomena were observed at Yusufpur, where left channel thalwegh erosion of 2 m was observed in 2005 and over time, and deposition of 4 m was observed during the 1998 flood event. From 1998 and long-run simulations, it has been observed that the river may become more prominent on the right side.

Water depth analysis has been conducted in order to ensure that the criteria for water availability are met. The required depth of the intake pipe has been set at 3.5 m to account for the diameter of the intake pipe, freeboard, and safety factors. According to Figure 21(a) and Figure 21(b), at the end of the monsoon, the lowest water depth of 3.5 m (calculated in relation to the SLWL) was discovered 77 m from the Godagari intake location during long run simulation. However, this bed level may not be the predicted maximum bed level for the duration of the simulation. Thus, in order to determine the worst-case scenario, the maximum bed level was determined statistically over the duration of the simulation. After determining the maximum bed level over the simulation period, it was determined from Figure 21(c) that the 3.5 m depth persists at 125 m from the Godagari intake location in the long run simulation. According to Figure 22(a) and Figure 22(b), the Harupur intake location has the lowest water depth of 3.5 m. However, as shown in the figure, the bed level never reaches 4.5 mPWD in the long run. At the end of the long run, the lowest water depth is 1.5 m at a distance of 20 m from the Harupur intake location. After determining the maximum bed level for the duration of the simulation, it is clear from Figure 22(c) that 3.5 m depth does not occur within a 500 m radius. However, the maximum depth of 2.18 m at a distance of 25 m from the Harupur intake location is only possible during extreme and bankfull flood events. According to Figure 23(a) and Figure 23(b), the lowest water depth of 3.5 m was discovered 32 m from the Yusufpur intake location at the end of the monsoon of bankfull and long run. However, it did not reach the 4 mPWD bed level during the extreme event. Again, after determining the maximum bed level throughout the simulation period, it can be seen in Figure 23(c) that the bankfull event occurred at a depth of 3.5 m at a distance of 40 m from the Yusufpur intake location. However, it never

Figure 21. (a) Depth contour calculated from SLWL and cross-sectional view of (b) bed level with SLWL (c) maximum bed level with SLWL near Godagari Intake location.

Figure 22. (a) Depth contour calculated from SLWL and cross-sectional view of (b) bed level with SLWL (c) maximum bed level with SLWL near Harupur Intake location.

Figure 23. (a) Depth contour calculated from SLWL and cross-sectional view of (b) bed level with SLWL (c) maximum bed level with SLWL near Yusufpur Intake location.

reaches the 4 mPWD bed level during an extreme event or for an extended period of time. To ensure the sustainability and safety of the intake structure, the worst-case scenario must be considered. It is obvious that 3.5 m depth is only possible at Godagari site from a distance of 125 m from the intake location for maximum bed level occurrence.

It is critical to protect these intake structures from damage caused by high velocity. Thus, this section of the report analyzes predicted velocity at study locations for various flood events. Figure 24 (left-top) illustrates the maximum current speed over the entire simulation period near Godagari for the extreme, bank-full flood event and the three-year long-run scenarios. A water depth of 3.5 m has been discovered 125 m from the Godagari intake location (stated in the previous section). 3.07 m/s velocity has been determined at 125 m distance from the Godagari intake location for the extreme flood event case (synthesized 1998) and long-run scenario. And in the case of the bankfull flood event (2005), the velocity was determined to be 2.80 m/s at a distance of 125 m from the intake location. These high velocities must be considered in order to ensure the long-term viability of existing protection works. Additionally, it is ideally considered when designing the water treatment plant’s intake structure. Whereas low velocities can provide information about the fall velocity and particle suspension state. Again, Figure 24 (left-bottom), illustrates the minimum current speed

Figure 24. Maximum and minimum velocity along the cross-section for extreme (synthesized 1998), bankfull (2002) and 3 year-long flood events at the Godagari (left), Harupur (middle) and Yusufpur (right) Intake locations.

throughout the simulation near the Godagari intake location. The minimum current speeds for the extreme, bankfull, and long-run have been determined to be 1.27 m/s, 0.29 m/s, and 0.06 m/s, respectively, at a distance of 125 m from the Godagari intake location. Figure 24 (middle-top), illustrates the maximum current speed over the entire simulation period near Harupur for the extreme, bank-full flood event and the three-year long-run scenarios. Within 500 m of the Harupur intake location, there is no water depth of 3.5 m. However, the 25 m distance can be used to calculate the velocity (stated in the previous section). For the extreme flood event scenario (synthesized in 1998) and the long run, a velocity of 1.98 m/s was observed 25 m away from the Harupur intake location. And in the case of the bankfull flood event (2005), the velocity was determined to be 2.41 m/s at a distance of 25 m from the intake location. The minimum current speed at Harupur (Figure 24 (middle-bottom)) has been determined to be 0.68 m/s and 0.25 m/s, respectively, at a distance of 25 m from the Harupur intake location. However, in the event of a prolonged run during a dry period, the total flow decreases to approximately 290 cumec (boundary condition stated in earlier section). Thus, the velocity of the near bank channel at Harupur decreases to almost zero during this period. Again, 3.5 m depth was discovered 40 m from the Yusufpur intake location (mentioned in the previous section). In the case of Yusufpur (Figure 24 (right-top)), the velocity was found to be 1.58 m/s for the extreme flood event (synthesized 1998); 1.49 m/s for the bankfull flood event (2005); and 2.47 m/s for the long-run case at a distance of 40 m from the Yusufpur intake location. Again, 24 (right-bottom), illustrates the minimum current speed over the entire simulation period near the Yusufpur intake location. For the bankfull flood event simulation, the minimum current speed at 40 m from the Yusufpur intake location was determined to be 0.07 m/s. A similar phenomenon to Harupur has been observed near the Yusufpur Intake. During extreme and long-run simulations, the velocity approaches zero.

The deepest bed level has been determined along the channel near Intake Locations using bathymetry data from different years (2002 to 2015). Additionally, the deepest points were identified from simulated results of an extreme flood event (synthesized in 1998), a bankfull flood event (2005), and a three-year-long run event. In 2002, the deepest point was discovered to be the closest to the bank, while in 2011, the deepest point was discovered to be the furthest from the bank. Whereas, after simulating the long-run event, it was determined that the deepest point is closest to the bank (106 m). While proximity to the bank may jeopardize the left bank’s stability and also existing bank protection work at Godagari, such proximity also ensures flow availability near the bank, which may be required for the proposed intake. After simulating an extreme flood event, a bankfull flood event, and the long run, thalwegs were extracted at the end of the simulation period and superimposed on a satellite image of November 2014, along with the initial thalweg line (extracted from pre-monsoon 2015). At Godagari, it has been determined that the initial thalweg is located 259 m from the intake location. Additionally, it was discovered that the thalweg is located 184 m, 59 m, and 96 m away from the intake location for the extreme flood event, bankfull flood event, and long run, respectively. At Harupur Intake, the deepest point was discovered to be 4 km from the left bank in 2002, and 338 m away from the left bank in 2015. After simulating an extreme flood event, it was determined that the deepest point is closest to the bank (47 m). Again, thalweg has been extracted and overlaid with the initial thalweg for the extreme flood event, bankfull flood event, and long-run at the end of the simulation period. The initial thalweg at Harupur has been determined to be 624 m from the intake location. Additionally, it was discovered that the thalweg is located 46 m, 120 m, and 465 m away from the intake location for the extreme flood event, bankfull flood event, and long run, respectively. At Yusufpur Intake, the deepest point was found to be very close (166 m) to the left bank in 2002, but far away (308 m) from the intake location in 2015. After simulating a three-year flood event, it was determined that the deepest point is located closest to the left bank (81 m). After simulating an extreme flood event, a bankfull flood event, and a three-year run, thalweg was extracted at the end of the simulation period and superimposed on a satellite image from November 2014. Though it was discovered that the main channel flows along the right bank during the simulation of the extreme flood event and the three-year run, the thalweg was extracted from the left channel for the purposes of this study. At Yusufpur, it has been determined that the initial thalweg is located 177 m from the intake location. Additionally, it was discovered that the thalweg is located 339 m, 317 m, and 223 m away from the intake location for the extreme flood event, bankfull flood event, and long run, respectively. According to the analysis above, the predicted thalweg for the bankfull and long-run events at Godagari is the closest. Although it was discovered 184 m away from the intake location during an extreme flood event (more than Harupur that is 46 m from intake location) (Figure 25).

Figure 25. (a) Location map of deepest points along the channel (b) thalweg at initial (pre-monsoon 2015) and end of the simulation period of extreme flood event, bankfull flood event and long-run near Godagari (left), Harupur (middle) and Yusufpur (right) intake locations.

From the aforementioned findings and discussions, it can be concluded that Godagari is the most suitable site for constructing the intake infrastructures; however, there is a large char that has been visually observed to be extremely dynamic in nature.

4. Concluding Remarks

As a result of the above description, it can be inferred that the curvilinear FDM numerical scheme (MIKE 21C) is effective for both understanding River dynamics and solving real-world engineering decision support tools. One notable advantage of this model is that it can generate accurate scenarios for river bends and local erosion and deposition events, although numerous river systems require a unique rectangular grid to address this problem. In this scenario, a combined curvilinear FDM and FEM approach may be investigated, which could be the subject of a future study.

Acknowledgements

The author wishes to express his gratitude to the Institute of Water Modelling (IWM) (https://www.iwmbd.org/) for providing training and computational resources for MIKE 21C simulations, as well as to the project leaders of the two aforementioned projects of the Institute of Water Modelling (IWM), whose suggestions and constructive comments significantly improved the author’s understanding of MIKE 21C. The author also wishes to express his gratitude to the editor and anonymous reviewers for their constructive comments.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Sarker, S. (2022) A Short Review on Computational Hydraulics in the Context of Water Resources Engineering. Open Journal of Modelling and Simulation, 10, 1-31. https://doi.org/10.4236/ojmsi.2022.101001
[2]	DHI (2022) MIKE 21C, Curvilinear Model—Scientific Documentation. Danish Hydraulic Institute, Horsholm, Denmark. https://manuals.mikepoweredbydhi.help/latest/Water_Resources/MIKE21C_Scientific_documentation.pdf
[3]	Sarker, S. (2021) Hydraulics Lab Manual. engrXiv, 1-66. https://doi.org/10.31224/osf.io/mxcvw
[4]	Sarker, S. (2021) Water Distribution (Pipe) Network Analysis with WaterCAD. International Journal of Engineering Development and Research, 9, 149-153. http://www.ijedr.org/papers/IJEDR2104028.pdf
[5]	Sarker, S. (2021) Investigating Topologic and Geometric Properties of Synthetic and Natural River Networks under Changing Climate. UCF STARS.
[6]	Sarker, S., Veremyev, A., Boginski, V. and Singh, A. (2019) Critical Nodes in River Networks. Scientific Reports, Nature Publishing Group, 9, No. 1, 1-11. https://doi.org/10.1038/s41598-019-47292-4
[7]	Reza, A.A. and Sarker, S. and Asha, S.A. (2014) An Application of 1-D Momentum Equation to Calculate Discharge in Tidal River: A Case Study on Kaliganga River. Technical Journal River Research Institute, 12, 77-86.
[8]	Sarker, S. (2021) Separation of Flood Plain Flow: 1-D Momentum Equation Solver. engrXiv, 1-7. https://doi.org/10.31224/osf.io/sjcmv
[9]	Sarker, S. (2021) Understanding the Complexity and Dynamics of Anastomosing River Planform: A Case Study of Brahmaputra River in Bangladesh. Earth and Space Science Open Archive, 1. https://doi.org/10.1002/essoar.10508926.2