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

Understanding, quantifying, and forecasting water flow and its behavior in environment is made possible by the use of computational hydraulics in con-junction 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 con-cerns.

Besides that, any partial space derivative of a function C in a centered space difference can be represented as Equation (4) 1 1 2 n n j j where, the temporal and spatial discretizations are t ∆ and x ∆ , respectively.
When space discretizations are performed on the present time level, the scheme 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: 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)).
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. C 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 solu-

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].

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].
where x and y are CarCS coordinates, and w denotes a weight factor defined by The condition of non-linear orthogonality is the boundary condition for this system, namely Additionally,

( )
, 0 f x y = denotes that x and y are points on a particular curve. The system described previously generates an orthogonal grid with streamlines ( ζ lines) and potential lines (η 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 ( 1 w = ), the system is said to be conformal; consequently, a river is the ideal example of a conformal mapping constraint, that can be specified as respectively, and x L and y L 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 x and p y . 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.

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)).
where, C is Chezy roughness coefficient, p, q Mass fluxes in the sand 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: where, ( ) , P Q are fluxes described in a CarCS and ( ) , p q 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 s i may be defined as Equation (24): where u and s R 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: where, s ∂ is the angle between bed shear stress and depth averaged shear stress or flow [2]. The parameter β can be defined as 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 ( sf λ ).
This results in the equation for advection-dispersion, which may be numerically solved using the MIKE 21C. where, 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.

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): where, C = concentration of suspended sediment, ε = is the turbulent diffusion Velocities are represented in the following way . 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].
where, c θ = modified critical shields parameter, 0 C θ = Critical Shields parameter in uniform shear flow and b z = bed level. If a model does not assume zero bed load transfer at a critical shear stress, the following bed load modifica- where, bl S = bed load as calculated from sediment transport formula and α = 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.
where, G = transverse bed slope factor (calibration coefficient) and tan s δ = bed shear direction change due to helical flow strength. In a fixed ( ) , x y coordinate system, bed slopes are calculated as Hence, the bed slope in the transverse direction will be where, φ = is the angle of stream line compared to fixed ( )

Morphology
Hydrodynamic and sediment transport models are combined in a morphological ( ) where, e S = lateral sediment supply from bank erosion, n = bed porosity, x S = total sediment transport in x direction and y S = 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 where, b E = bank erosion rate, S = near bank sediment transport, h = local water depth, z = local bed level and , , α β γ calibration coefficients specified in the model. Therefore, this model also mimic the lateral processes in morphological models that are entitled as planform module.

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.

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  Finally model generated cross-sections for intake location also suggested that intake is located at safe reach (Figures 2-7).

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      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 3 /s and 139 m 3 /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. Open Journal of Modelling and Simulation      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 Open Journal of Modelling and Simulation       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 de- 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    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 Open Journal of Modelling and Simulation 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). Open Journal of Modelling and Simulation 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.

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.