The term “hydraulics” is concerned with the conveyance of water that can consist of very simple processes to complex physical processes, such as flow in open rivers, flow in pipes, the flow of nutrients/sediments, the flow of groun d- water to sea waves. The study of hydraulics is pr imarily a mixture of theory and experiments. Computational hydraulics is very helpful to quantify and predict flow nature and behavior. The mathematical model is the backbone of the computational hydraulics that consist s of simple to complex mathemat i- cal e quations with linear and/or non -linear terms and ordinary or partial di f- ferential equations. Analytical solution to these mathematical equations is not feasible in the majority of cases. In these consequences , mathematical models are solved using different numerical techniques and associated schemes. In this manuscript , we aim to review hydraulic principles along with their m a- thematical equations. Then we aim to learn some commonly used numerical technique s to solve different types of differential equations related to hydra u- lics. Among them, the Finite Difference Method (FDM), Finite Element M e- thod (FEM) and Finite Volume Method (FVM) have been discussed along with their use in real -life ap plications in the context of water resources eng i- neering. Keywords Computational Hydraulics, Finite Difference Method, Finite Element Method, Finite Volume Method 1. Introduction Hydraulics is a subfield of fluid mechanics that studies the behavior of water at rest and in motion. This term is frequently used in connection with hydrody nam- ics to refer to the forces acting between water and its boundaries and their effects How to cite this paper: Sarker , S. (2022) A Short Review on Computational Hydraulics in the Context of Water Resources Enginee r- ing. Open Journal of Modelling and Simul a- tion, 10, 1-31. https://doi.org/10.4236/ojmsi.2022.101001 Received: October 5, 2021 Accepted: December 6, 2021 Published: December 9, 2021 Copyright © 20 22 by author (s) and Scientific Research Publishing Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). http://creativecommons.org/licenses/by/4.0/ Open Access S. Sarker DOI: 10.4236/ojmsi.2022.101001 2 Open Journal of Modelling and Simulation on subsequent flow patterns. To measure and predict flow behavior, it is critical to apply physical conservation rules, which will be discussed i n detail in the su b- sequent section. Although they attempted to solve numerous fluid mo tion pro b- lems strictly theoretically, assuming an ideal fluid (frictionless or non-viscous), their results were still extremely effective in specific scenarios and condit ions. On the other hand, hydraulics developed mathematical formulas to handle real -world fluid flow problems using experimental data from numerous laboratory exper i- ments and field observations. Ludwig Prandtl’s essential idea of boundary layer combined hyd rodynamics and hydraulics. Additionally, the illu strious research of Reynolds, Froude, Prandtl, and von Karman persuades us that studying fluids requires a combination of theory and experiment. Here, we can trace the evolu-tion of hydraulics all the way bac k to 200 BCE [1] ( Figure 1 ). Hydraulics originated around the time of Archimedes, who established the now -famous rules of floatin g bodies. Following Galileo (1564 -1642), Gastelli (1577- 1644), Torricelli (1600- 1647), and Guglielmini (1655- 1710) developed mathematical models for barometers, flow from containers, steady flow, and r e- sistance in rivers. When Newton (1642 -1727) published his famed laws of m o- tion, they expanded our quantitative understanding of fluid resistance in terms of velocity gradient and drag on spheres. Following Newton, Daniel Bernoulli, Leonhard Euler, Clairaut, d’Alembert, Lagrange (1736 -1813), Laplace (1749- 1827 ), and Gerstner (1756 -1832) contributed significant concepts on fluid flow and surface waves to hydrodynamics. Additionally, De Pitot invented a tube for me a- suring speed; Chezy devised a formula for open channel resistance; Borda co n- ducted orifice -related experiments; Bossut installed a towing tank; and Venturi Figure 1. Contributors in hydrodynamics and hydraulics ( sources : Google image). S. Sarker DOI: 10.4236/ojmsi.2022.101001 3 Open Journal of Modelling and Simulation experimented with the flow in varying cross -sections, all of which contributed to the field of hydraulics’ enrichme nt [1]. Coulomb (1736- 1806) conducted flow resistance studies concurrently with Ernst Brothers (1795 -1878) and Weber (1804- 1891), all of whom made signif i- cant contributions to the field of hydraulics. As a result, the Buron (1790- 1873), Fourneyman (1802- 1867), Coriolis (1792- 1843), Francis (1815- 1892), and Russel (1808- 1882) all made significant contributions to engineering. The Hagen (1797- 1889), Poiseuille (1799- 1869), and Sazon Weisbach (1806 -1871) all co n- tribute significantly to pipe flow, whereas Saint -Venant (1797- 1886) contributes to open channel flow. On the other hand, Dupuit (1804 -1866), Bresse (1822- 1883), Basin (1829- 1917), Darcy (1803- 1858), and Manning (1816- 1897) conducted r e- search on bo th open channel hydraulics and pipe flow. William Froude (1810- 1879) and Robert Froude (1846- 1924) in particular provide a valuable criterion for clas- sifying the flow based on their ship model tests. Additionally, Osborne Reynolds did a successful experiment based on stream line analysis to discern between the laminar and turbulent flow. Navier (1785- 1876), Cauchy (1789- 1857), Poisson (1781 -1840), Saint -Venant (1797- 1886), Boussinesq (1842 -1929), Stokes (1819 - 1903), Lord Rayleigh (1842- 1919), Lamb (1849- 1934), Helmholtz (1821- 1894), and Kirchoff (1824 -1887) all made significant contributions to the development of theoretical and applied hydrodynamics during the nineteenth century. Additio n- ally, Euler’s equation of motion for an ideal (non -viscous) fluid had advanced sig- nificantly. However, it does not ac count for certain signifi cant observed results, such as the decrease in pipe pressure. As a result, engineers in practice created their own empirical hydraulics. Recent computational hy draulics is one of several do- mains of science in which computer technologies enable an intermediate mode of operation between pure theory and experiment. This discipline, which is fre quently referred to as computational hydraulics, is more of a synthesis of hydraulics and hydrodynamics than an independent development. The objective of comput a- tional hydraulics is to use computers to mimic various physical processes invol v- ing seas, estuaries, rivers, canals, and reservoirs. We will exp lore funda mental h y- draulic processes and their quantification in this article, focusing on water re- source applications such as open channels, pipes, groundwater, and coastal waves [1]. 2. Hydraulics of Open Channels Open channel flows deal with the flow of water in open channels, where the pres-sure at the surface of the water is ambient and the pressure at each portion is pro- portional to the depth of the water. Pressure head is the ratio of pressure and the specific weight of water ( Pγ). Elevation head ( Z) is the height of the section under consideration above a datum and Velocity head (22Vg ) is attributed to the average velocity of flow in that vertical section. Hence, the total head can be expressed by the Equation (1). 2 2PVZgγ=++ (1) S. Sarker DOI: 10.4236/ojmsi.2022.101001 4 Open Journal of Modelling and Simulation In an open channel, water flow is primarily due to the gradient of the head and gravity. The Open Cha nnels are primarily used for irrigation supply, indu s- trial and domestic water supply [2]. 2.1. Conservation Laws Mass Conservation, Momentum Conservation and Energy Conservation are the main conservation laws used in open Channel. A control volume commonly co n- sidered under the Con servation of Mass, the overall mass change in the control volume due to inflow and outflow, is equal to the net rate of weight change in the control volume. This results to the classical continuity equation balancing the inflow, out flow and the storage ch ange in the control volume. Since we only find water that is regarded as incompressible, that is, it is possible to disregard the impact on density [2] . AV Q= (2) The rate of change of momentum in the control volume under the Conserv a- tion of Momentum is equal to the net forc es acting on the control volume. Since the water under consideration is flowing, external forces are acting upon it. This leads ultimately to the Newton ’s second law of motion [2]. 2QyA ConstantgA+= (3) Energy conservation states that neither energy can be generated nor destroyed, it only changes its form. The energy would be in the form of potential energy and kinetic energy, primarily in open channels. The elevation of the water parcel is due to potenti al energy, while the kinetic energy is due to the movement of the parcel. In the context of open channel flow the total energy due these factors b e- tween any two sections is conserved, that leads to the classical Bernoulli’s equ a- tion. 2 2VZ y Constantg++ = (4) This equation must account for the energy loss between the two sections when used between two sections, which is due to the resistance of bed shear etc. to the flow [2]. 2.2. Types of Open Channel Flows The flow in an open channel is classified as Sub critical flow, Super Critical flow, and Critical flow depending on the Froude number (rF), where the Froude number can be described as 5: rVF gy= (5) Types of open channel flow can be further categorized on the basis of time and space parameters (see Figure 2(a) ). Flow is said to be steady when discharge S. Sarker DOI: 10.4236/ojmsi.2022.101001 5 Open Journal of Modelling and Simulation Figure 2. (a) Types of open channel flows ; (b) Specific energy curve for a given discharge. does not change along the course of the channel flow. On the other hand, as the discharge varies over time, flow is said to be unsteady. When both the depth and discharge are the same at any two sections of the channel, flow is said to be un i- form. Furthermore, flow is s aid to be gradually varied whenever the depth changes gradually along the channel, whenever the flow d epth changes rapidly along the channel the flow is termed rapidly varied flow. Whenever the flow depth grad u- ally varies due to the change in discharge, the flow is called a spa tially varied flow [2]. Steady uniform flow (kinematic wave), steady non-uniform flow (diffusion wave) and unsteady non-uniform flow (dynamic wave) can be seen as types of po - ssible flow based on the above classifi cation. 2.3. Specific Energy and Specific Force Specific energy is defined as the energy acquired by a section of water due to its depth and the velocity with which it is flowing, specific Energy E is given by, 2 2VEyg= + (6) where y is the depth of flow at that section and v is the average velocity of flow. Specific energy is minimum at critical condition (see point C in Figure 2 ), and the corresponding depth at point C is termed as critical depth yc. This yc can be measured for a specific section under steady flow [2]. Specific force is defined as the sum of the momentum of the flow passing through the channel section per unit time per unit weight of water and the force per unit weight of water [2]. 2QF yAgA= + (7) The specific forces of two sections are equal provided that the external forces and the weight effect of water in the reach between the two sections can be i g- S. Sarker DOI: 10.4236/ojmsi.2022.101001 6 Open Journal of Modelling and Simulation nored. The term specific force is very useful to quantify the property of rapidly varied flow such as hydraulic jump [2] . At the critical state of flow the specific force is a minimum for the given discharge. When the specific energy is minimal, flow is critical. Also, the flow has to go through critical conditions if the flow changes from sub critical to super critical or vice versa. Sub -critical flow -the depth of the flow is higher, although the velocity is lower and the super -critical flow-the depth of the flow is lower, but the v elocity is higher. Critical flow is the flow over a free over- fall. For Specific energy to be a m inimum d0 dE y=: 2 3dd1ddE QA yy gA= − (8) However, ddAT y= ∗ , where T is the width of the channel at the water su r- face, then applying d0 dE y=, will result in following: 22 2 321; ;c c cc cc ccQT A A V Q T Tg gA gA= = = (9) For a rectangular channel c c cAyT= . Following the derivation for a rectangular channel, 1c r cVF gy= = (10) The same notion applies to trapezoidal and other cross- sections [2]. The cri t- ical flow state describes an unique relationship between depth and discharge that is very useful in the design of flow measurement structures. 2.4. Uniform Flows This is one of the most important principles of open channel flow. The most important uniform flow equation is Manning’s equation stated as : 23 12 1V RSn= ∗∗ (11) where R = the hydraulic radius = A/p and p = wetted perimeter = ()0,f yS , y = depth of the channel bed, S0 = bed slope (same as the energy slope, Sf) and n = the Manning’s dimensional empirical constant. Uniform Flow concept is used in most of the open channel flow design. The uniform flow implies that there is no acceleration of the flow due to the weight portion of the flow being balanced by the resistance of the shear of the bed. In terms of discharge the Manning’s Equ a- tion (11) is given by: 23 12 1Q AR Sn=∗∗ ∗ (12) This is a non linear equation in y the depth of flow for which most of the computations will be made. Derivation of uniform flow equation is given below, S. Sarker DOI: 10.4236/ojmsi.2022.101001 7 Open Journal of Modelling and Simulation where sinWθ = weight component of the fluid mass in the direction of flow, 0τ = bed shear stress and Px∆ = surface area of the channel. The force ba l- ance equation can be written as: 0 00 sin 0 sin 0 sinAW px Ax pxpθτ γ θτ τ γ θ− ∆=→ ∆ − ∆=→ = (13) Now A/p is the hydraulic radius, R, and sinθ is the slope of the channel So. The shear stress can be expressed as: 2 02fVCτρ= (14) where fc is resistance coefficient, V is the mean velocity ρ is the mass de n- sity. Therefore the Equation (14) can be written as: 22 2fo o o fVgC RS V RS V RSCργ= →= →=  (15) where C is Chezy’s constant. For Manning’s equation, 1 61.49Rn= (16) 2.5. Gradually Varied Flow Flow is said to be gradually varied whenever the depth of flow changed gradually. The governing equation for gradually varied flow is given by: 2d d 1of rSSy x F−= − (17) where the variation of depth y with the channel distance x is shown to be a fun c- tion of bed slope oS, Friction Slope fS and the flow Froude number rF. This is a non linear equation with the depth v arying as a non linear function (Figure 3). Figure 3. (a) Steady uniform flow in an open channel ; (b) Total head at a channel section. S. Sarker DOI: 10.4236/ojmsi.2022.101001 8 Open Journal of Modelling and Simulation Gradually varied flow can be derived from the conservation of energy at two sections of a reach of length x∆, can be written as: 22 12 1222ofVVy Sx y Sxgg++ ∆ =++ ∆ (18) Now, let 21 yy y∆= − and 22 2 21 d 2 2 d2VV vxg g xg−= ∆  , Then the above equation becomes: 2d d2ofvy SxS x xxg∆= ∆− ∆− ∆   (19) Dividing through x∆ and taking the limit as x∆, approaches zero gives us: 2dd d d2ofyvSSx xg+= −  (20) After simplification, can be done in terms of Froude number and the general differential equation can be written as Equations (21) to (23): 2d d d1d2ofSS y x V yg− = +  (21) 2 22 2 2 23 2d d 2d 1 d2 d d 22rV Q QA QFyg y y D gA gA gA   −− = = = =−          (22) 2d d 1of rSSy x F−= − (23) Numerical integration of the gradually varied flow equation will give the water surface profile along the channel. Depending on the depth of flow where it lies when compared with the normal depth and the critical depth along with the bed slope compared wit h the friction slope different types of profiles are formed such as M (mild), C (critical), S (steep) profiles. M (mild) -If the slope is so small that the normal depth (Uniform flow depth) is greater than critical depth for the given discharge, then the sl ope of the channel is mild. C (critical)- if the slope if the slope ’s normal depth equals its critical depth, then we call it a Critical slope, de- noted by C. S (steep) -if the channel slope is so steep that a normal depth less than critical is produced, then the channel is Steep, and water surface profile desi g- nated as S (see [1] [3]). 2.6. Rapidly Varied Flow This flow has very pronounced curvature of the streamlines. It has pressure di s- tribution that cannot be assumed to be hydrostatic. The rapid variation in flow regime often takes place in short span. When rapidly varied flow occurs in a sudden -transition structure, the physical characteristics of the flow are basically fixed by the boundary geometry of the structure as well as by the state of the flow. Channel expansion and channe l cont raction, Sharp crested weirs and Broad S. Sarker DOI: 10.4236/ojmsi.2022.101001 9 Open Journal of Modelling and Simulation crested weirs can be seen as examples. Sp ecific force before and after the flow r e- gime can be considered same i.e., 12FF= [2]. 2.7. Unsteady Flows When the flow conditions vary with respect to time, we call it unsteady flows. Some terminologies such as wave is defined as a temporal or spatial variation of flow depth and rate of discharge. Wave length is the distance between two adj a- cent wave crests or trough. Amplitude is the height between the maximum water level and the still water level [4] . Wave celerity ( ) is the relative velocity of a wave with respect to fluid in which it is flowing with V. Absolute wave velocity (wV) is the velocity with respect to fixed refere nce as below: wVV= ± (24) Plus sign indicates the wave is traveling in the flow direction and vice versa. For shallow water waves 0gy= , where y0 = undisturbed flow depth. Un- steady flows occur due to following reasons: surges in power canals or tunnels, surges in upstream or downstream channels produced by starting or stopping of pumps and opening and closing of control gates. Furthermore, navigation looks can also generate waves in the navigation channel. Flood waves in streams, rivers, and drainage channels due to rainstorms and snow -melt, tides in estuaries, bays and inlets. Unste ady flow commonly encountered in an open channels and deals with translatory waves. A translatory wave is a gravity wave that propagates in an open channel and results in appreciable displacement of the water particles in a direction parallel to the flow. For purpose of analytical discussion, unsteady flow is classified into two types, namely, gradually varied and rapidly varied u n- steady flow. In gradual ly varied flow the curvature of the wave profile is mild, and the change in depth is gradual. In the rapidly va ried flow the curvature of the wave profile is very large and so the surface of the profile may become vi r- tually discontinuous. Continuity equati on for unsteady flow [5] in an open cha n- nel becomes: 0fV yyDVx xt∂ ∂∂+ +=∂ ∂∂ (25) For a rectangular channel of infinite width, the Equation ( 25) may be written as: 0qy xt∂∂+=∂∂ (26) When the channel is to feed laterally with a supplementary discharge of q′ per unit length, for instance, into an area that is being flooded over a dike, then the equation becomes: 0qyqxt∂∂′++=∂∂ (27) The general dynamic equation for gradually varied unsteady flow is given by (see [5] [6]): S. Sarker DOI: 10.4236/ojmsi.2022.101001 10 Open Journal of Modelling and Simulation 10y VV V x g x gtα ′ ∂ ∂∂+ +=∂ ∂∂ (28) 3. Hydraulics of Pipe Flows Pipe flows are mainly due to pressure difference between two sections. The total head here also consists of the pressure head, elevation head and velocity head. The principle of continuity, energy, momentum is also used in this type of flow. For example, w hen designing a pipe, we use continuity and energy equations to obtain the appropriate pipe diameter. Then, applying the momentum equation [5] [6], for a given discharge, we get the forces that act on bends. The key factors in the design and operation of a pipeline are head losses, pressures and stresses acting on the material of the pipe, and discharge. Head loss for a given discharge relates to flow effici ency; i.e., an optimum size of pipe will yield the least overall cost of installa tion and operation for the desired discharge. Choosing a small pipe results in low initial cost s, but due to high energy costs from significant head losses, subsequent costs may be excessively high. The design of conduit should be such that it needs least cost for a given discharge. The hydraulic aspect of the problem requires applying the o ne dimensional steady flow form of the energy equation. Energy equation can be written as: 22 1 11 2 2 2 1222p tLPV P VZh Z hhggαα γγ+ ++ = + + ++ (29) where, P γ = pressure head, 2 2V gα = velocity head, Z = elevation head, hp = head supplied by a pump, ht = head supplied by a turbine and hL = head loss between 1 and 2. The kinetic energy correction factor is given by α, and it is defines as, where u = velocity at any point in the section. ()3 3duA VAα=∫ (30) α has minimum value of unity when the velocity is uniform across the se c- tion. α has values greater than unity depending on the degree of velocity vari- ation across a sectio n. For laminar flow in a pipe, velocity distribution is par a- bolic across the section of the pipe, and α has value of 2.0. However, if the flow is turbulent, as is the usual case for water flow through the large conduits, the ve-locity is fairly uniform over most of the conduit section, and α has value near unity (typically: 1.04 1.06α<< ). Therefore, in hydraulic engineering for ease of application in pipe flow, the value of α is usually assumed to be unity, and the velocity head is then simply 2 2V g. The head supplied by a pump is directly related to the power supplied to the flow as ppP Qhγ= , similarly, if head is supplied to turbine, the power supplied to the turbine will be ttP Qhγ= . The head loss term Lh accounts for the co n- version of mechanical energy to internal energy (heat), when this conversion oc curs, S. Sarker DOI: 10.4236/ojmsi.2022.101001 11 Open Journal of Modelling and Simulation the internal energy is not readily converted back to useful mechanical energy, there- fore it is called head loss. Head loss results from viscous resistance to flow (fri c- tion) at the conduit wall or from the viscous dissipation of turbulen ce usually o c- curring with separated flow, such as in bends, fittings or outlet works. Head loss is due to friction between the fluid and the pipe wall and turbulence within the fluid. The rate of head loss depends on roug hness element size apart from velo c- ity and pipe diameter. Further the head loss also depends on whether the pipe is hydraulically smooth, rough or somewhere in between. In water distribution sy s- tem, head loss is also due to bends, valves and changes in pipe diameter. Head loss for steady fl ow through a straight pipe can be written as below: 0wrA PAτ=∆ (31) 04PDτ∆= (32) 2 0 8 fVτρ= (33) 2 2PVhfDgγ∆= = (34) This is known as Darcy -Weisbach equation. hS=L , is slope of the hydra u- lic and energy grade lines for a pipe of constant diameter. Head loss in laminar flow can be calculated using Hagen- Poiseuille equation, which gives: 232VS Dgµ ρ= (35) Combining above equations with Darcy -Weisbach equation provides : 64f VDµ ρ= (36) That can be further written in terms of Reynolds number: 64 efR= (37) This relation is valid for 1000eR< . In turbulent flow, the friction factor is a function of both Reynolds number and pipe roughness. As the roughness size or the velocity increases, flow is wholly rough and f depends on the relative roug h- ness. Where graphical determination of the friction factor is acceptable, it is possible to use a Moody diagram. This diagram gives the friction factor over a wide range of Reynolds numbers for laminar flow and smooth, transition, and rough turbulent flow. The quantities shown in Moody Diagram (see details in [2]) are dimensionless so they can be used with any system of units. Minor losses caused by valves, bends and changes in pipe diameter. This is smaller than fric tion losses in straight sections of pipe and for all practical purposes ignored. Minor losses are significant in valves and fittings, which creates turbulence in excess of that pr o- duced in a straight pipe. A minor loss coefficient K may be used to give head loss as a function of velocity head. 2 2mVhKg= (38) S. Sarker DOI: 10.4236/ojmsi.2022.101001 12 Open Journal of Modelling and Simulation A flow coefficient Cv which gives a flow that will pass through the valve at a pressure drop of 1psi may be specified. Given the flow coefficient the head loss can be calculated as: 42 6 218.5 10 2vDVh gC= × (39) The flow coefficient can be related to the minor loss coefficient by: 2 6 218.5 10 vDK C= × (40) Major losses are due to friction between the moving fluid and the inside walls of the duct. The Darcy -Weisbach formula is generally considered to calculate ma- jor losses in pipe s. This method is gen erally considered more accurate than the Hazen -Williams method. Additionally, the Darcy -Weisbach method is valid for any liquid or gas. Moody Friction Factor f can be calculated as below: ( ) ( )8 62 0.964for 2100 i.e. laminar flow 1.325for 5000 10 i.e. turbulent flow i.e. 10 10 5.74ln3.7ee e e eVDRRR e fRD e D Rν −−= ≤→ =≤≤ → →→ ≤≤ + (41) Major loss in pipes can also be calculated Using Hazen- Williams friction loss equation : 0.63 0.54where and or circular pipe 4f HWh DV kC R S S R f= = =  (42) Hazen -Williams is only valid for water at ordinary temperatures (40˚ - 75˚F). The Hazen -Williams method is very popular, especially among civil eng ineers, since its friction coefficient ( CHW) is not a function of velocity or duct (pipe) diameter. Hazen -Williams is simpler than Darcy -Weisbach for calculations where one can solve for flow -rate, velocity, or diameter. When the flow conditions are change d from one steady state to another, the intermediate stage flow is referred to as transient flow. This occurs due to design or operating errors or equip ment malfunction. This t ransient state pressure causes lots of damage to the network system [4] [7]. Pressure rise in a close conduit caus ed by an instantaneous change in flow velocity. If the flow velocity at a point does vary with time, the flow is unsteady. The terms fluid transients and hydraulic transients are used in practice [4]. Consider a pipe length of length L Water is flowing from a constant level upstream reservoir to a valve at downstream. Assume valve is instantaneously closed at time 0tt= from the full open position to half open position. This r e- duces the flow v elocity through the valve, thereby increasing the pressure at the valve. The increased pressure will produce a pressure wave that will travel back and forth in the pipeline until it is dissipated because of friction and flow cond i- tions have become steady again. This time when the flow conditions have b e- come steady again, let us call it 1t. So the flow regimes can be categorized into : 1) S. Sarker DOI: 10.4236/ojmsi.2022.101001 13 Open Journal of Modelling and Simulation Steady flow for 0tt<; 2) Transient flow for 01t tt<< ; and 3) Steady flow for 1tt>. Transient Transient -state pressures are sometimes reduced to the vapor pressure of a liquid that results in separating the liquid colu mn at that section; this is referred to as liquid -column separation. If the flow conditions are r e- peated after a fixed time interval, the flow is called periodic flow, and the time interval at which the conditions are repeated is called period. The analysi s of transient state conditions in closed conduits may be classified into two categories: lumped may be classified into two categories: lumped -system approach and di s- tributed system approach. In the In the lumped system lumped system approach the conduit walls are assumed rigid and the liquid in the conduit is assumed i n- compressible, so that it behaves like a rigid mass, other way flow variables are functions of time only. In the distributed system approach the liquid is assumed slightly compressible. Ther efore flow velocity varies along the length of the co n- duit in addition to the variation in time. The 1D− form of momentum equa- tion for a control volume that is fixed in space and does not change shape may be written as: ()()22 dddout inAV x AV AVtρρ ρ= +−∑∫ (43) If the liquid is considered to be incompressible and the pipe is rigid, then at any instant the velocity along the pipe will be same, i.e., ()()22 out inAV AVρρ= . Substituting for all the forces acting on the control volume we get: ()0dsindA A D AVtγ ατ ρ+ −= π    (44) where, 2 2Vhgγ= −  , α = pipe slope, D = pipe diameter,  = pipe length, γ = specific weight of fluid and 0τ = shear stress at the pipe wall. Frictional force is replaced by fhAγ , and 0 sin Hh α = + and fh from Dar cy-Weisbach friction equation. The resulting equation yields: 22 0d 22 dfV V VHDg g gt− −= (45) When the flow is fully established d0dV t=. The final velocity 0V will be: 2 0 012V fHDg= + (46) We use the above relationship to get the time for flow to establish: 22 02ddDVtDf VV=+−  (47) If the flow changes are rapid, it is important to take into account fluid co m- pressibility [4]. Changes are not instantaneous in the system, because in the piping system, pressure waves travel back and forth. The walls of the pipe must be rigid and the liquid slightly compressible [2] [4]. Assume that the flo w velo c- S. Sarker DOI: 10.4236/ojmsi.2022.101001 14 Open Journal of Modelling and Simulation ity at the downstream end is changed from V to VV+∆ , thereby changing the pressure from P to PP+∆ . The change in pressure will produce a pressure wave that will propagate in the upstream direction. By assuming that the velocity re f- erence system travels with the pressure wave, the unsteady flow condition can be translated into steady flow. Using momentum equation with control volume a p- proach to solve for p∆ (see details in [2] [4]). The system is now steady, the mo- mentum equation now yield: ()( )()( )()() P A P P A Va V Va V AVa Va A ρρ ρ − +∆ = + +∆ +∆ + +∆ − + + (48) 4. Hydraulics of Ground Water Flow Velocity, discharge and head are the essential quantities used to characterize the flow of groundwater. They are followed by the Darcy’s law and some other m e- trics related to soil physical characteristics. The specific discharge is defined as the volume of water flowing through a unit area of soil per unit time, which has the unit of velocity. Specific discharge usually denoted by q and is defined as b e- low: Qq A= (49) where Q is the flow rate of water through a cross sectional area A. The average velocity of fluid at a certain point of the porous medium is called Seepage velo c- ity. The seepage velocity is represented by v, and is defined as: eeqQvn An= = (50) where ne is the effective porosity of the soil. The effective porosity, ne is defined as ratio of the volume of continuous pore spaces (open to groundwater flow) of a soil sample to the total volume of the sample. For ground water flo w, the total head denoted by h is defined as the sum of the pressure head, and the elevation head, z as follows: pPh zh zgρ= += + (51) where P is the pressure at any point, ρ is the density of water at the prevailing temperature. If streamlines are not perpendicular to the vertical piezometer, then hp may include the contribution of the vertical component of the velocity head for most groundwater flows, the velocity head const itutes a very small part of the tot al head and is usually neglected. The important relationship is that the specific discharge at a certain point is proportional hydraulic gradient at that point and the discharge occurs in the direction of decreasing gradi ent, which is known as Darcy’s law : d dxhqKx= (52) where qx is the specific discharge along the x direction, d dh x is the gradient of S. Sarker DOI: 10.4236/ojmsi.2022.101001 15 Open Journal of Modelling and Simulation head causing the flow along the x direction, and K is known as hydraulic co n- ductivity. The hydraulic conductivity can be further expressed as: gKρ µ= (53) If the hydraulic conductivity K is independent of position within a geologic formation, the formation is homogeneous. In this case, K is a constant. On the other hand, if K is dependent of position then it is called heterogeneous form a- tion. If the hydraulic conductivity K is independent of the direction of me a- surement at a point in the geologic formation, the formation is isotropic at that point. On the other hand, if K varies with the direction of measurement at a point, the formation is anisotropic at that point. In an aquifer with layered het e- rogeneity, flow may occur parallel to the layers or across the layers or both. To simplify the flow calculations, it is possible to use estimates of equivalent K. If the layers are horizontal then the equivalent K for a flow parallel to the layers, Kx, and the same for a flow across the layers, Kz, are given by the following: 1n ii x iKDK ==∑ (54) 1izni iKD K== ∑ (55) where 1in iD==∑ ; iD and iK are thickness and hydraulic conductivity of the ith homogeneous layer. The Darcy law is for a homogeneous and isotropic medium. Assuming that the principal axes of K are aligned with the coordinate axes then the governing equation for a three dimensional anisotropic confined aquifer can be written as: x y zsh h hhK K KSx xy yz z t ∂∂∂∂∂∂ ∂ ++= +  ∂ ∂∂ ∂∂ ∂ ∂   (56) Here,  is source or sink term. If the medium is homogeneous and isotro p- ic, xyzKKKK= = = . In addition, if the aquifer is horizontal and has a co n- stant thickness b, then Kb= , and S S bS= , in two dimensions the equation becomes, 22 22h h Sh t xy∂∂ ∂+=∂ ∂∂ (57) where  is called the transmissivity of the aquifer, and S is called the storage coefficient. The above equation called the diffusion equation. For the steady state condition this equation becomes the well known Laplace equation: 22 220hh xy∂∂+= ∂∂ (58) The two dimensional flow equation for an unconfined aquifer is given by: x yyh hhKh Kh Sx xy y t ∂ ∂∂ ∂ ∂+= + ∂ ∂∂ ∂ ∂  (59) S. Sarker DOI: 10.4236/ojmsi.2022.101001 16 Open Journal of Modelling and Simulation Here, yS is called the specific yield, which is the quantity of water that can be drained out of a saturated volume of porous medium due to unit lowering of the water table. The above equation is also known as the Boussinesq equation. It is a nonlinear, second ord er partial differential equation but linear in 2h. 5. Water Wave Dynamics Small -amplitude or linear wave theory is the most elementary wave theory d e- veloped by Airy (1845). With a large range of wave parameters, this theory is easie r to apply, offering a reasonable approximation of wave characteristics. It is necessary to assert the key assumptions of this theory as below. Homogeneous, ideal and incompressible fluid, irrotational flow, surface tension and Coriolis e f- fect is neglected , pressure at the free surface is uniform and constant, small a m- plitude and bed is a horizontal, fixed, impermeable boundary. Two dimensional irrotationality can be stated by Equations (60) and (61): uxφ∂=∂ (60) and: wzφ∂=∂ (61) where, φ is known as velocity potential is a scaler fu nction [8]. In addition, in- compressible flow implies that there is another function termed as the stream function ψ (62-63). xzφψ∂∂=∂∂ (62) and: zxφψ∂∂=−∂∂ (63) ψ is orthogonal to the potential function φ. The lines of constant values of the potential function ( i.e., equipotential lines) are mutually perpendicular to the lines of constant values of the stream function. Both φ and ψ satisfy the Laplace equation which governs the flow of an ideal water (see Equations (64 ) and (65) waves ). 22 220 xzφφ∂∂+= ∂∂ (64) and: 22 220 xzψψ∂∂+= ∂∂ (65) The speed at which a wave form propagates is termed the phase velocity or wave celerity C LT= . Where, the length L is the horizontal distance between corresponding points on two successive waves and the period T is the time for two successive crests to pass a given point. The most commonly used terms to describe water wave dynamics are depicted in Figure 4 , where, a progressive S. Sarker DOI: 10.4236/ojmsi.2022.101001 17 Open Journal of Modelling and Simulation Figure 4. Definition of elementary terms used to define a progressive sinusoidal w ave [8]. wave represented by the spatial variable x, temporal variable t and η denotes the displacement of the water surface relative to the SWL . As a result their phase can be def ined as kx tθω= − , where k and ω are described as wave number (2kL=π) and angular frequency (2 Tω=π) respectively. In this consequence, the wave profile can be writt en as fo llowing Equation (55): ()22cos cos 2x tHa kx t LTηωππ = = −− (66) Therefore, the amplitude of wave can be depicted as @Ha= , where, H is the wave height. Wave motion can be often defined in terms of some dimensionless parameters such as wave steepness (H L), relative depth (d L) and relative height (H d). It is important to point out that sometimes relative depth and relative height quantified as kd and ka since they differ only by a constant factor of 2π. Water waves usually classified into three following categories based on relative depth: shallow water wave (1020d L<< i.e., 010kd<π<∗ ), transitional water wave (11 20 2d L<< i.e., 10kdπ<<π) and deep water wave (1 2d L< <∞ i.e., kdπ< <∞ ). In general, wave phase speed can be expressed as: () tanh2gTkd =π (67) Therefore, the measurement of wave length can be: ()2 tanh2gTL T kd =π= (68) In addition to different wave parameter, the local flow velocities and acceler a- tions during the passage of a water wave must often be calculated using Equ a- S. Sarker DOI: 10.4236/ojmsi.2022.101001 18 Open Journal of Modelling and Simulation tions (69) -(72). () ()()cosh cos2 coshkz d H gTu kx tL kdω+= − (69) () ()()sinh sin2 coshkz d H gTw kx tL kdω+= − (70) () ()()cosh sincoshxkz d u gHkx tt L kdαω+ ∂ = = −∂π (71) () ()()sinh coscoshzkz d w gHkx tt L kdαω+ ∂ ==−−∂π (72) Water particle displacements from mean position for shallow/transitional -water (i.e., 1 2d L<) and deep -water waves ( i.e., 1 2d L>) can be tracked by following equations: () ()()cosh sin2 sinhkz d Hkx tkdξω+=−− (73) () ()()sinh cos2 sinhkz d Hkx tkdζω+= − (74) Assuming, () ()cosh 2 sinhkz d H kd+= and () ()sinh 2 sinhkz d H kd+= , we can r e- write Equation (73) and 74 in the form of 22 221ξζ+= . Therefore, water particle displacements follows elliptical orbit for shallow/transitional -water while circ u- lar orbit for deep -water waves. In other words, for shallow/transitional- water 22HL d=π and 2Hz d d+= ; for deep -water waves 2exp2Hz L= =π . Water pressure waves are classified into three main components: dynamic co m- ponent due to acceleration, static component of pressure and atmospheric pre s- sure. Mathematically expressed using Equation (75). () ()()cosh cos2 coshakz d HP g kx t gz Pkdρ ωρ+′= −− + (75) Therefore, relative pressure will be a PPP′= − . If pressure response factor is denoted by () ()cosh coshzkz d Kkd+= then the Equation (75) becomes ()z P gK zρη′= − . Besides, pressure force, it is possible to divide water wave energy into traditional potential (2 1 16pE gH Lρ= ) and the kinetic (2 1 16kE gH Lρ= ) form of energy. Therefore, the total energy per wave length and unit width can be written as 2 1 8E gH Lρ= . Furthermore, the rate of energy S. Sarker DOI: 10.4236/ojmsi.2022.101001 19 Open Journal of Modelling and Simulation transfer by water wave is known as energy flux, that c an be quantified using the Equation (76). ()2 11 21 8 2 sinh 2kd LP gH kd Tρ = +  (76) The Equation (76) can simply written as the form of g P En EC= = , where gC is commonly known as group velocity. Therefore, group velocity can be quantified using Equation (77). ()121 2 sinh 2gkd LCn kd T = = +  (77) Coastal hydraulics usually considers problems near the shoreline and the key processes that can affect a w ave as it propagates from deep into shallow water includes: shoaling, refraction, reflection, diffraction, breaking and damping. The transformation of the wave form due to interaction with bathymetry is known as wave shoaling. Wave refraction is the change s in the direction of wave propag a- tion due to differences in wave velocity along the crest, which are generally ill u- strated by lines drawn perpendicular to the wave crest in the direction of wave propagation known as wave rays. The angle of the wave along a ray follows Snell’s law. Wave diffraction is the bending of wave crests (changes in direction) due to along crest gradients in wave height. The point at which wave form becomes u n- stable and breaks is known as breaking point and the phenomenon is known as wave breaking. Wave breaking occurs when water particles at the crest travel much faster/farther than water particles in the trough. Wave breaking may be classified in four types as spilling, plunging, collapsing, and surging [8]. 6. Differential Equations for Water Movement Water movement is the key process in Water Resources Engineering. Most m a- thematical models describing water movement are consisting of different ial eq u- ations deriving from the fundamental principles of the conservation of mass, energy and momentum. The continuity equation for one- dimensional unsteady gradually -varied flow, the Belanger equation, and unsteady flow equation are some example of diffe rential Equations. The highest derivative is called the order of the differential equation and the power of its derivative of the highest order is called the degree of the differential equation. In g eneral, the differential equ a- tions for water movement are second -order partial differential equations. The con- servation laws in physics can be described by a PDE of the general form (Equ a- tion (78)): 222 220C C C CCa b c d e fC gxy x y xy∂ ∂ ∂ ∂∂+ + + + + +=∂∂ ∂ ∂ ∂∂ (78) where a, b, c, d, e, f and g may be constants or functions of x and y or of C. If a, b, c, d, e, f and g are constants or functions of x and y, the equation is linear. If ()240b ac−> , the above equation is hyperbolic. For example wave equation is a S. Sarker DOI: 10.4236/ojmsi.2022.101001 20 Open Journal of Modelling and Simulation hyperbolic equation. If ()240b ac−= , the above equation is parabolic. For e x- ample diffusion equation is a parabolic equation. If ()240b ac−< , the above eq- uation is elliptic. The Laplace equation is an elliptic equation. Besides general form of PDE the conservation law can be seen as advection -diffusion process which is discussed in the following section. 7. Modeling Advection -Diffusion Processes Advection -Diffus ion process can be expressed by an equation known as adve c- tion-diffusion equation which comes from the conservation of contam inant mass. Advection refers to the transport of contaminants by the moving water where contaminants travel at the same velocity as the mass of water in which they are dissolved. Diffusion refers to the spreading of the contaminants in all direction due to the turbulence and non -uniform velocity distribution. Here, molecular diffusion is much smaller than that of turbulent diffusion. In addition to being advected and diffused, the amount of a contaminant dissolved in water may increase or decrease in time d ue to chemical reactions with other agents present in the environment. When such reactions are absent, the co ntaminant process is called a conservative process meaning that the total mass of the di s- solved substance remains constant through the transport pr ocess. When such reaction processes are present, the contaminant process is called non -conservative indicating that the mass of the contaminant is growing or de caying at a certain rate. 2 2 fCC Cu D rC tx x∂∂∂= −+ − ∂∂ ∂ (79) Wave equation, diffusion equation and Laplace equation is also example of Advection -Diffusion process. Analytical solutions for advection -diffusion -reaction processes can be found only for simplified problems. For complicated boundary conditions or initial conditions, one needs numerical methods to get solution for the problem. 8. Transport of Particles in Water Movement In addition to water movement, transports of sediment, salt and pollutants lead to differential equations. Advection -diffusion process is a common process in trans - ports of sediment, salt and pollutants. Advection refers to the transport of co n- taminants by the moving water where contaminants travel at the same velo city as the mass of water in which they are dissolved. Diffusion refers to the spreading of the contaminants in all direction due to the turbulence and non- uniform v eloci- ty distribution. In an open channel flow, molecular diffusion is much lower co m- pared to turbulent diffusion. The general differential form of conservation of con- tami nant mass equation can be expressed as: y x zq q q Cstxyz∂∂ ∂ ∂+++=∂∂∂∂ (80) S. Sarker DOI: 10.4236/ojmsi.2022.101001 21 Open Journal of Modelling and Simulation where, C is the mass concentration per control volume. In addition, Fick’s law define the flux ( ,xyqq and zq) due to an advection and a diffusion process, that can be written as: xxCq uC Dx∂= −∂ (81) yyCq vC Dy∂= −∂ (82) zzCq wC Dz∂= −∂ (83) Apart from the general type of contaminants transporting particle, some other dissolved gas transport along with their reaction rate can also be part of hydra u- lics, which is very important for the supply and treatment of water. 9. Hydraulics Used in Water Supply and Treatment The concen tration of a gas (such as O 2) that can be present in solution is g o- verned by: 1) the solubility of the gas as defined by Henry’s Law ; 2) temperature; 3) presence of impurities (such as salinity); and 4) the partial pressure of the gas in the atmosphere. The saturation concentration of a gas dissolved in a liquid is a function of the partial pressure of the gas. Henry’s law states that g g TxPHP= , where, Pg and xg is the mole frac tion of the gas in air and in the liquid (water) respectively, PT is the total pressure and H is the Henry’s constant that depends on the type of gas. For example, dissolved oxygen ( DO) is measured in standard solution units such as mmol/L, mg/L, mL/L, or ppt. O 2 saturation is calculated as the percent of DO relative to a theoretical maximum concentration given the temperature, pressure, and salinity of the water. Well- aerated water will usually be 100% saturated. In general, lower temperatures, lower salin ity, and higher atmospheric pressures lead to higher values of dissolved O 2. Higher barometric pressures lead to higher values of saturated DO. The correction for barometric pressure may be written as 0760PuDO DOu−′=−, where, P is the barometric pressure and u is the vapor pressure in mm (Hg). Microorganisms in water di- gest organic material as food in the presence of O 2 that lead to introduce a new term commonly known as the biochemical oxygen demand ( BOD ). The BOD on a water b ody depends on the microbial population dynamics which includes a lag phase, a constant growth phase, a stationary phase, and a decay phase. A c- cording to the popular mathematical model for the growth of BOD (i.e., the ex- ponential model where the BOD grows asymptotically to the so -called ultimate BOD of the water), BOD can be written as () 1ekt tu BOD BOD−= − , where k is the reaction rate constant (day−1). The rate constant at T temperatures is given by ()20 201.047T Tkk−= . When a stream of untreated or partially treated sewage is discharged into a large body of water, the wastewater is diluted and their ult i- mate BOD (i.e., oL), DO and temperature (oT) of the river -wastewater imm e- S. Sarker DOI: 10.4236/ojmsi.2022.101001 22 Open Journal of Modelling and Simulation diately after m ixing can be given by ( )ww rrQ L QL Q+ , ( )ww rrQ DO Q DO Q+ and ( )ww rrQ T QT Q+ respectively, where, wr QQ Q= + , wQ = waste water flow rate and rQ = river flow rate. Initial oxygen deficit after mixing is given by 0 O satD DO DO= − . The Streeter- Phelps equations (see 84) describe the DO “sag curve” obtained as a result of simultaneous deoxygenation and reoxygenation that occurs when wastewater is discharged into receiving stream. According to the Streeter- Phelps model 5, the oxygen deficit after time t (0t= at the instant of mixing) is given by : ()0 0 ee ed rr kt kt kt d t rdkLDDkk− −−= −+− (84) where, kd and kr are deoxygenation (day−1) and reoxygenation rate (day−1) re- spectively. The time at which the dissolved oxygen is critical is given by : 1lndo r o d o r c r d do dkL kD kD ktk k kL k−+= − (85) At this time, the critical (maximum) O 2 deficit is given by : edckt do c rkLDk−= (86) The distance from the point of mixing to the critical location is given by ccx vt= , where, v = average stream velocity. In the treatment of wastewater, operations in which some transformation is effected through the use of chemical reactions are termed chemical unit oper a- tions. The rate at which a chemical reaction progresses depends on the nature of the reaction and hydraulics used to describe that nature is called reactor hydra u- lics. Chemical reactions usually classified as zero order, first order, second order, and saturation reactions. Zero -order reactions may be expressed as d dCkt=− that leading to 0 tC C kt= − . Similarly, for a first -order reaction follows d dCkCt=− that leading to 0ekt tCC−= and for a second -order reaction follows 2 d dCkCt=− that leading to 011 tktCC= + . A reaction is classified as a saturation reaction if the ra te of reaction saturates as the reaction progresses. This may be expressed as d dC kC t aC=+ that leading to 0 0 lnt tCa C C ktC+−= . Required hydraulic detention time for various orders of reaction also depicted in Figure 5, which is very useful for water treatment reactor. There are several further co m- putational hydraulics applications that address water quality concerns [9]. 10. Hydraulics of Sediment Transport Sediment transport function is very crucial to predict complex river [10] hy- drodynamics and morphodynamics. For decades, various methods have been used to establish sediment transport functions or formulas. They are ofte n differ S. Sarker DOI: 10.4236/ojmsi.2022.101001 23 Open Journal of Modelling and Simulation Figure 5. (a) Dissolved oxy gen sag curve (Streeter -Phelps) ; (b) Hydraulic detention time for various or ders of reaction. drastically from each other and from observations in the field. These formulas have been used to solve engineering and environmental problems. The sediment transport mechanics for cohesive and non- cohesive materials are different. Most of the sediment transport functions are proposed under the assumption of non- cohesive sediment. In this section, we will review of the basic concepts and a p- proaches used in the derivation of incipient motion criteria and sediment tran s- port functions. Due to the stochastic nature of sediment movement, inc ipient motion is critical for the study of sediment transport, channel degradation and stable channel design. In other words, it is difficult to define precisely at what flow condition a sediment particle will begin to move, hence it depends more or less on definition of incipient motion. Significant progress has been made on the study of incipient motion, both theoretically and experimentally. Out of all ran - ges, the emphasis here is more on Shear Stress Approach and Velocity A pp- roach. A sediment particle is at a state of incipient motion when one of the following conditions is satisfied: LSFW= , DRFF= and ORMM= , where, overturning moment due to FD and FL, and resisting moment due to FL and WS (see Figure 6). One of the most prominent and commonly used incipient motion criteria is the Shields diagram (1936) based on shear stress approach, where, where the central premise is the shear stress τ, the difference in density between sedime nt (Sρ) and water (wρ), the diameter of the particle ( ds), the kinematic viscosity and the gravitational acceleration can be grouped into two dimensionless quantities that provides seminal plot to determine i ncipient motion (see Figure 6). *12 c w sdUdτ ρ νν  = (87) and: ()1cc sS w S s wdg dττ ρρ ργρ=− −  (88) S. Sarker DOI: 10.4236/ojmsi.2022.101001 24 Open Journal of Modelling and Simulation Figure 6. Sediment transport modes and forces acting on a single spherical sediment particle along with Shields and HJulstrom diagram [11]. where, Cτ = critical shear stress at initial motion and *U = shear velocity. Rouse (1939), White (1940), Brooks (1955), Liu (1958), Yang (1973), Govers (1987), Yalin and Karahan (1979) provides useful changes on the Shields dia gram for practical use. While Lane (1953) developed stable channel design curves for tra- pezoidal channel s with different typical side slopes, however, recently USBR (1987) synthesized stable channel design criteria based on the critical shear stress required to move sediment particles in channels under different flow and sedi- ment conditions [11]. On the other hand, Fortier and Scobey (1926) carried out an extensive field survey of the maximum permissible mean velocity values in the canals of diffe r- ent m aterials. Hjulstrom (1935) carried out a detailed study of the movement of uniform materials at the bottom of the channels, which provides a convenient diagram (see Figure 6 ) based on the velocity approach. The relationship between sediment size and average flow velocity for erosion, transport and sedimentation can be illustrated by the Hjulstrom curve. Yang (1973) and Vanoni (1977) fur- ther improved the velocity approach that Govers (1987) and Talapatra and Ghosh (1983) tested experi mentally [11]. Yang used theoretical fluid mechanics to calculate critical velocity by providing FL, FD, FR, and WS equations, which i n- volve the fall velocity of the sediment particles. Sediment particle fall velocity is one of the crucial parameters used in most sedi ment transport functions or fo r- mulas. Different methods have been developed for the computation of sediment particle fall velocity. Rubey’s (1933) introduce fall velocity formula ()1SsF dg Gω= − , where, F is a function of particles diameter, specific gravity S. Sarker DOI: 10.4236/ojmsi.2022.101001 25 Open Journal of Modelling and Simulation of sediment, kinematic viscosit y of water and g, however for particle greater than 1mm, F can be considered as 0.79 [11] . Yang and Simoes (2002) use a shape factor cSF ab= for natural sand, where, a, b, c the length of the longest, the intermediate, and the shortest mutually perpendicular axes of the particle, re- spectively. In water treatment, fall velocity commonly calculated based on the type of particle settling, that can be classified into four types: discrete particle settling, flocculent settling, hinde red or zone settling, and compression settling. In addition, discrete particle fall velocity depends on the laminar, transition or turbulent region of flow that is defined in earlier section. Regime, regression, probabilistic, and deterministic approaches are the basic approaches used in the derivation of sediment transport functions or formulas. The concept of “regime” is similar to the concepts of “dynamic equilibrium” and “hydraulic geometry.” Various collections of regime equations have been formulated by various inve s- tigators, such as Blench (1969), Kennedy (1895) and Lacy (1929). Lacy’s (1929) regime equation describing the relationships among channel slope S, water di s- charge Q, and silt factor Sf for sediment transport can be written as 53 160.0005423SfS Q= . Leopold and Maddock’s (1953) hydraulic geometry rela- tionships commonly known as bW aQ= , fy cQ= and mV kQ= , where, W = channel width, y = channel depth, V = average flow velocity, Q = water di s- charge and a, b, c, j, k, m are local constants. Yang et. al. (1981) applied the unit stream power theory for sediment transport, and the hydraulic geometry rela-tionships between Q and S as jS iQ= , where, i, j, constants. Shen and Hung (1972), Karim and Kennedy (1990) proposed sediment transport function adopting regression approach by considering flow velocity, sediment discharge, bed-form geometry, and friction factor. The regression approach may provide fairly accurate results due to the fact that the sediment transport is such a co m- plex phenomenon and no single hydraulic parameter or combination of par a- meters can be found to describe sediment transport rate under all conditions. Einstein (1950) w as a pioneer in sediment transport studies focused on a prob a- bilistic approach, which is under the assumption of the beginning and ceasing of sediment motion can be expressed in terms of probability and the movement of bedload is a series of steps followed by rest periods. In addition to that, Einstein used the hiding correction factor and the lifting correction factor to get the best fitted theoretical result with the observed experimental data. Although Einstein’s bedload function is not common in enginee ring applications due to complex computational procedures, it has however been used as a theoretical foundation for the formulation of other transport functions. For example, Colby and Hembree (1955) used modified Einstein methods for the computation of to tal bed-material load. On the other hand, the deterministic approach is the exi s- tence of one- to-one relationship between independent and dependent variables. Commonly used independent variables are water discharge, average flow veloc i- ty, shear stress, and energy or water surface slope. More recently, the use of S. Sarker DOI: 10.4236/ojmsi.2022.101001 26 Open Journal of Modelling and Simulation stream power and unit stream power have gained increasing attention as impo r- tant parameters to calculate sediment concentration. Bagnold (1966) introduced the stream power concept based on classical physics, while, Engelund and Ha n- sen (1972), and Ackers and White (1973) later used the concept as the theoret i- cal basis for developing their sediment transport functions. Yang (1972) defines unit stream power as the rate of energy per unit weight of water available for transporting water and sediment in an open channel with reach length x and total drop of y can be written as d dd d ddy xyVSt tx= = , where, V = average flow velocity, and S = energy or water surface slope. Besides that, Velikanov (1954) derived his transport function from the gravitational power theory and Pache- co-Ceballos (1989) derived a sediment transport function based on power ba l- ance between total power available and t otal power expenditure in a stream. This principle is also applicable to understanding the evolution of the natural channel network [12] [13] [14] [15] [16]. In this section, we comprehensively reviews basic approaches and theories used in the determination of noncohesive sediment concentration. In the fo l- lowing section we will discuss some basic finite differences numerical techniques to solve governing equations. Finite element and finite volume are another two other numerical technique which will be covered under subsequent sections. 11. Finite Difference Method In Finite difference methods ( FDM ), partial derivatives are replaced with finite differences approximations, which results in a set of algebraic equations. These algebraic equations can be solved explicitly or implicitly. FDM uses a structured grid and thus requires grid transformation in multidimensional flow case. A typical finite difference grid can be shown in Figure 5 . Finite difference method is popular because of its simplicity in formulation and implementation. Based on different ways of discretization, differen t numerical schemes are developed. Some of them are shown as bellows: In a forward time difference, any partial time derivative of a function C can be expressed as Equation (89) : 1nn jjCC C tt+− ∂=∂∆ (89) In a forward space difference, any partial space derivative of a function C can be expressed as Equation (90): 1nn jjCC C xx+− ∂=∂∆ (90) In a backward space difference, any partial space derivative of a function C can be express ed as Equation (91): 1nn jjCC C xx−−∂=∂∆ (91) S. Sarker DOI: 10.4236/ojmsi.2022.101001 27 Open Journal of Modelling and Simulation In a centered space difference, any partial space derivative of a function C can be expressed as Equation (92): 11 2nn jjCC C xx+−− ∂=∂∆ (92) where t∆ and x∆ is the time and space discretization respectively. When space discretizations are done on the current time level, then the scheme is known as an explicit scheme. For example Equations (90) to (92) are the examples of e x- plicit space discretized scheme. On the other hand, when space discretizations are done on the next time level or using both the current and next time levels, the generated scheme is known as an implicit scheme. For example, if Equation (91) is discretized as follows : 11 1nn jjCC C xx++ −− ∂=∂∆ (93) or: 11 11 1 2nn n n jj j jCC C C C xx x++ −−−− ∂= + ∂∆∆ (94) then both will generate implicit schemes. An explicit scheme will not generate any systems of equation. Therefore, an explicit scheme can be solved directly. That means, unknowns at any nodal value can be calculated without solving other nodal values. On the contrary, an impl icit scheme will generate a system of eq u- ations, and therefore, unknowns at any nodal value can not be solved indepe n- dently. Apart from that anot her popular implicit scheme is Box Finite Difference Scheme, which is also known as four -point implicit scheme. According to this scheme, time derivative can be expressed as : 11 11 1 2nn nn jj jjCCCC C t tt++ ++−− ∂= + ∂∆∆ (95) and space derivative can be expressed as : 11 11 1 2n n nn j j jjCC CC C x xx++ ++−− ∂= + ∂∆∆ (96) and C can be expresses as : ( )11 111 4nn n n jj j j C CC C C++ ++ = +++ (97) In order to improve accuracy, further schemes such as the Crank -Nickelson Scheme, Box Finite Difference Scheme, Lax -Wendroff Scheme, MacCormack Scheme were proposed to solve advection -diffusion processes. Prior to solving the advection -diffusion equation or any governing equation, the initial condition of the system has to be specified. Solution of the governing equations also r e- quire boundary conditions. The most common type of boundary conditions are Dirichlet, Neumann, and Cauchy boun dary conditions. Dirichlet condition o c- curs when a portion of the boundary is at a prescribed concentration level ( i.e., C S. Sarker DOI: 10.4236/ojmsi.2022.101001 28 Open Journal of Modelling and Simulation specified), Neumann condition occurs when a portion of the boundary has sp e- cified flow which is crossing the boundary curve normally (i.e., C n∂ ∂ specified), and Cauchy condition is a mixed boundary condition which represents some combination of the Direchlet and Neumann conditions on the boundary curve. It is important to mention that numerical method used to ref ormulate governing equations in an approximate manner where continuous operations replaced by discrete operations that causes truncation error. In addition, size of the discr e- tized solution domain ( i.e., cells) causes round off error ( Figure 7). Model di s- cretization and numerical solution should be consistent ( i.e., correct), conve r- gent ( i.e., 00xε ∆→→ → ), and stable ( i.e., error should not increase). 12. Finite Element Method Another numerical technique used to obtain an approximate solution to a g o- verning differential equation is the Finite Element M ethod. This is relatively new to the field of computational hydraulics and is supported by rigorous mathe- matical theory. FEM based on variational form of PDE, which derived u sing the method of weighted residuals. In this method, under a given boundary condition, the PDE (i.e., strong form) is used and the residuals are determined for the a s- sumed solution due to the fact that the assumed solution is not exact. Then push the residual to zero at various points or intervals using different weighted resi- dual method. There are three commonly used weighted residual methods ( i.e., Collocation, subdomain, and Galerkins method) in FEM. Among them Gale r- kins method is popular, where the integral of residual equation R, multiplied by a weighting function, w, is pushed to be zero. The weighting functions are cho- sen to be the same form as each part of the approximate so lution. If the residual, ()Ru , for a given PDE , should equal 0 for the true solution. In other words, for approximate solution, non -zero residual measures accuracy huu≈. Here, we require that the weighted residual equal zero: Figure 7. (a) Schematic 1D space and time discretization ; (b) Computational error as a function of x∆. S. Sarker DOI: 10.4236/ojmsi.2022.101001 29 Open Journal of Modelling and Simulation ()d 0; wR u x w W Ω= ∀∈∫ (98) and hu is a linear combination of basic functions, can be ex pressed as ()()() ,h jj ju xt u t x φ =∑ (99) Eventually, transforms governing equation into a system of linear equations (Ax b=), which can be solved using matrix operations. The FEM can be explain by three phases prepossessing, processing/solution and post -processing. In pr e- possessing phases, model geometry is discretized or subdivided into a finite nu m- ber of smaller pieces, each piece is called an element, and all intersection points are called nodes. The process of dividing geometry is called meshing, which is one of the key steps of the FEM. After that, one element is considered at a time and assumes a shape function that depicts the physical behavior of the element. The stiffness matrix of each element is then developed and assembled them into a global matrix system for the entire geometry. Apply the required boundary condition and initial condition to the global matrix system. In solution phase, primary unknown is computed from the global system of linear equati on. Finally, in the post -processing phase, other derived variables are determined on the basis of the nodal value of the primary unknowns. For example, SLIM is a numerical model. The full form of SLIM is Second -ge- neration Louvain- la-Neuve Ice- ocean Model . It is a hydrodynamic model based on finite element method ( FEM ). The special advantage of finite element m e- thod is that, user can use this technique for unstructured grids. One can refined computational grid arbitrarily in the areas of interest that focu sing the comput a- tional power where it is needed. In this model no need of nested grids. A single model is able to resolve both the large- scale features, such as in the open sea and also small -scale phenomena in rivers, coasts and estuaries. SLIM consists o f a 1D river model, a 2D depth averaged model and 3D model. For 1D river model co n- sists of linear river where variable river width and cross -section, 2D model the domain is divided into triangular elements allowing accurate representation of complex geomet ry and 3D model uses triangular prismatic elements that are formed by extending the 2D mesh vertically. 13. Finite Volume Method Besides computational hydraulics, majority of computational fluid dynamics (CFD) codes based on Finite Volume Method (FVM ). Thi s method is based on integral form of PDE , where the governing equation integrated over control v o- lume. For example, the integral form of 1D advection -diffusion conservation law can be written as: dd dfvvvCC Cv u vD vt x xx∂ ∂ ∂∂= −+∂ ∂ ∂∂∫∫∫ (100) The integral form is based on a control volume v. This v volume can be cell - centered or vertex -centered as shown in Figure 8 . Typically, the special model domain divides into discrete control volumes. Then an assumption is made on S. Sarker DOI: 10.4236/ojmsi.2022.101001 30 Open Journal of Modelling and Simulation Figure 8. (a) Cell- centered and (b) Vertex- centered finite volume mesh along with (c) Midpoint approximation rule. how the value of integral changes within the control volume v, usually average value is taken at the center of the control volume v. ()()1,d iiV iC t C xt xV=∫ (101) The suitable time integration scheme is selected to solve the integral of Equ a- tion (101). Midpoint rule, Trapezoidal rule and Simpson’s rule are the most common approximations of integrals. In this method, all fluxes usually co n- verted to surface integrals. This is an exact procedure due to the fact that the v o- lumes are defined as polygons or polyhedra . 14. Concluding Remarks Hydraulics is a prerequisite for those interested in a career in Water Resources Engineering. The fundamental principles of hydraulics, the computational a s- pects of hydraulic system analysis and design, and various engineering applic a- tions of the concepts are discussed in this manuscript. This manuscript can assist in developing an app ropriate learning environment for comprehending and conducting pipe and open channel flow research. Acknowledgements The author wishes to express his gratitude to the editor and anonymous revie w- ers for their constructive suggestions and comments, which substantially i m- proved the manuscript’s presentation and interpretation. S. Sarker DOI: 10.4236/ojmsi.2022.101001 31 Open Journal of Modelling and Simulation Conflicts of Interest The author declares no conflicts of interest regarding the publication of this p a- per.