A semi coupled shallow water model for vertical velocity distribution in an open channel with submerged exible vegetation

Flow-vegetation interactions modify the instream roughness and flow characteristics in the river and estuaries. This study proposes a new quasi three-dimensional hydrodynamic framework to compute the vertical velocity profile in an open channel having submerged flexible vegetation. A modified form of two-dimensional depth-averaged shallow water equations coupled with vegetal drag forces is derived and applied in the simulation. The explicit second-order accurate TVD McCormack predictor-corrector finite difference method with operator splitting technique is used to solve the governing equations in MATLAB. The TVD approach is robust and gives accurate results free from numerical oscillations. The bending profile of the flexible stems under various flow events is calculated from the cantilever beam theory. The vertical velocity profile in the vegetation layer and the free water layer is estimated from Reynold's stress equation and Shannon's entropy theory. The present model is used to replicate some popular experimental test cases. Results indicate a conservative and robust model performance under different flow conditions and patch density. Quantitative analysis of the predicted results is carried out using two statistical indices and found satisfactory.


Introduction
The floodplains and the riparian vegetations modify the land use-land cover pattern, channel roughness, turbulence intensity, shear stress characteristics, and water quality of a river system (McBride et al. 2007). These alterations, in turn, affect the flow velocity, solute, and sediment transport processes (Salama and Bakry 1992;Wang et al. 2018). The presence of vegetation also impacts the habitat richness and biodiversity in the river ecosystem (Verschoren et al. 2016). Thus, a precise understanding of flow dynamics under the vegetated environment is essential for implementing habitat restoration and conservation strategies in any natural domain (Jalonen and Jarvela 2014;Verschoren et al. 2016).
Based on stem flexibility, vegetations are classified as flexible and rigid vegetation (Aberle and Järvelä 2013). Flexible canopies are ubiquitous and occupy a more extensive section of the flow area in a vegetated terrain than the rigid stems (Armanini et al. 2005). In recent years, flexible vegetation has gained a lot of research interest because of its potential impact on the flow and velocity structure (Huai et al. 2019;van Veelen et al. 2020). Unlike the rigid vegetation, the stem of the flexible vegetations bends to various degrees depending upon the velocity of flow, height, and flexural rigidity (Wilson 2007;She et al. 2014). The different degrees of bending also have a distinctive impact on the flow resistance (Luhar and Nepf 2011). For instance, with the increase in velocity, the flexible vegetation's flow resistance reduces, leading to further reconfiguration of the plants and the vegetation drag (Wu 2008;Shields Jr et al. 2017). Thus, the flow characteristics with flexible vegetations are more complex to deal with as compared to the case of rigid vegetation.
The presence of flexible vegetation, especially under submerged conditions, modifies the flow domain's velocity profile (Nepf 2012). For the non-vegetation channels, the logarithmic expressions are well accepted in defining the vertical velocity profile (Tang 2019). However, with submerged vegetation, the velocity profile becomes more complex (Nepf 2012). In those conditions,two-layer or three-layer methods are adopted for velocity distribution above and below the vegetation layer (Keramaris et al. 2015). These methods divided the total water depth into a vegetated region and a free water region and distributed the velocity profiles using different approximations methods (Afzalimehr et al. 2019). For instance, the study carried out by Pu et al. (2019) divided the flow region into multiple layers and predicts the velocity profile in those layers from mixing length models and eddy viscosity models. Klopstra et al. (1997) proposed a two-layer analytical model to distribute the velocity profile in a vegetated open channel. The model calculates the turbulent stresses in the vegetation region from the Boussinesq hypothesis. Huai et al. (2013) divide the flow regime into a free surface layer, an outer layer, and an inner vegetated layer. Yang and Choi (2010) distributed the velocity profile using a two-layer model in an open channel with immersed vegetation. A large number of studies that attempted to predict the vertical velocity distributions using similar approaches can be found in the literature (Huai et al. 2014;Wang et al. 2015). Though these methods are simple to apply, they have their shortcomings. The eddy viscosity and mixing length models mentioned above fail to capture the anomaly in the flow structure. For instance, near the flow separation regions in a river, the mixing length models fail to predict the velocity dip (Chin and Murray 1992), which is often the case, especially in natural flow domains. These models also require seven to eight parameters, viz shear velocity, energy slope, hydraulic radius, etc., for velocity computation. The difficulties associated with measuring these parameters are another shortcoming of the analytical formulations for field application (Chen and Kao 2011).
The probabilistic approach can tackle these difficulties to obtain the vertical velocity distribution in open channels. Following Chiu and Said (1995), considering the time-average velocity as a probabilistic variable, the velocity distribution may be obtained by maximizing the Shanon's entropy. The velocity distribution obtained from this method has an advantage as it required few input parameters and able to describe the velocity distribution in all circumstances. The entropy model distributes the vertical velocity profile at a particular location primarily from the depth-averaged velocity information in a channel cross-section (Kubrak et al. 2008). With this method, the location of maximum velocity can also be obtained. A dimensionless entropy parameter M can be derived from the mean and maximum velocity information at a channel crossection (Baruah and Sarma 2020). This parameter M represents the channel cross-section, slope, bedform, geometric shape and can characterize the velocity distribution in an open channel flow (Kumbhakar and Ghoshal 2017). The methodology section of this paper provides further details about this method.
The measurement of depth-averaged velocity, particularly in a natural domain, is a tedious task. However, two-dimensional shallow water models are excellent numerical tools to predict the streamwise and transverse velocities even under worst flow events (Blanckaert 2001). The ease of getting the depth-averaged velocity from the shallow water models makes the authors proposing an integrated formulation by linking the entropy formulation with the 2D hydrodynamic models. A full quasi three-dimensional numerical framework can be developed for vertical profile estimation in open channels evidently under a vegetated environment by coupling these two approaches. As per the author's knowledge, no such study has been reported by linking the entropy theory with the shallow water model to compute the vegetated channel flow's vertical velocity profile.
In this work, the authors developed a semi-coupled approach to compute the three-

2.1.Derivation of the governing equation and numerical scheme
In this work, full Saint Venant equations (continuity and momentum equation) are solved in boundary-fitted coordinates. A modified form of the shallow water equation is derived and applied to compute the flow domain's hydrodynamic parameters. The deviatoric form of twodimensional depth-averaged shallow water equation presented by Liang and Borthwick (2009) Where hu and hv are the momentum fluxes in x and y direction, and are the bed slope in x and y direction, η is the water surface elevation,S fx and S fy are friction slope in x and y direction, is the bed level elevation.
Under the steady-state condition, this formulation is further modified, hereby balancing the flux gradient and the source term in the x-direction momentum equation (equation-1) The above simplification modifies the momentum equation and also the source terms. The Where η is the water surface elevation (m), hu and hv are the unit discharges in x and y direction, and water surface gradients, h is the flow depth (m) measured from the bed up to the free surface, and are the bed friction components in the longitudinal and transverse direction, is the density of water and and are the drag force due to the vegetation. U and V transformed velocity in the computational domain in , coordinate system and expressed as The drag force from the vegetation in the x and y direction is expressed as A splitting algorithm is employed in the solution by dividing the two-dimensional equations into two one-dimensional equations and solved at each time step. Using the operator splitting techniques proposed by Strang(1968), equation (3),(4),(5) can be expressed in matrix form as The sweep in direction is presented below. A similar process is adopted in the ƞ direction for the next sweep by changing the subscript to j. Predictor- Corrector- The function B () in equation (15) is defined as B(x) =0.5*C*[1-ϕ(x)] and super bee flux limiter function is given as Where C = f ( ), is the courant number given as = * (( * ( 2 + 2 )+√ * ℎ)) (20)

2.2Deflection of the flexible vegetation
The Imposedaverage load P on flexible stems due to the upstream water head is expressed as (Huai et al. 2013) Where m is the vegetation density per unit area, s f is the energy slope, ρ is the density of water(KN/ 3 ), H is the total flow depth (m), θ is the angle of rotation after bending,EI is the flexural rigidity (N. m 2 )and g is the acceleration due to gravity. The curve length(S)of the element after bending is found from (Huai et al. 2013) Where h v = height of the vegetation (m) and z=any depth along the vertical(m)

Velocity distribution in the vegetated layer
In

Shannon's entropy theory for velocity distribution in the free water layer
The average and maximum velocity in a channel cross-section can be related to the entropy model as Chiu C L (1989) applied the principle of maximum entropy to the Shannon entropy (Shannon 1948) to obtain the following expression for vertical velocity distribution.
Where is the velocity at any point in a vertical, is the maximum velocity occurring in that vertical, and is the entropic parameter. For the convenience of representing the flow, a special coordinate system was defined such that the velocity is dependent on . The velocity is assumed to increase from zero ( = ) to a maximum value ( = ).
Therefore, is a function of .

Initial and Boundary condition
The initial and boundary conditions are essential for the smooth functioning of the numerical model. In the hydrodynamic simulation, a discharge value is provided as the initial condition.
The solid walls are simulated as a no-slip boundary (Anderson and Wendt 1995) viz both the streamwise and transverse velocity components are assigned with a zero value. At the upstream boundary, constant discharge and the flow depth corresponding to that discharge are set at the downstream boundary.

Semi coupled approach
The present study used the Semi-Coupled approach to link the vertical velocity distribution  (25) and (26) to compute the velocity profile within the vegetation zone and at the crown of the vegetation.
Further, the calculated velocity at the crown is used in the free layer's velocity distribution.
Whereas for the top free water layer, the entropy parameter for the channel cross section is calculated from the depth-averaged velocity and the maximum velocity. The depth-averaged velocity is obtained from the shallow water model and used in equation (27)

3) Application of the numerical model
The performance of the present model is assessed by setting up the model with series of experimental test cases carried out by Amreeva and Kurbak (2007). The model computed output is then compared with the experimental outputs reported in their laboratory flume experiment. Amreeva and Kurbak (2007) conducted the experiments in the hydraulic laboratory containing vegetation without foliage in the department of Hydraulic structure, Warsaw Agricultural University. The experiment was performed in a glass-walled flume of 16m length, 0.58m wide, and 0.6m depth. In their experiment, the cylindrical stems of elliptical cross-sections having diameter d1=.00095m and d2=.0007m were placed in a removal plate of 3m length, .58 m wide made of PVC. The longitudinal and transverse velocity profiles were measured using a programmable electromagnetic liquid velocity meter.
The schematic diagram of the experimental channel is available in (Amreeva and Kurbak 2007). The details of the different test runs of the experimental flume are enlisted in Table-1.
In the present study, for numerical simulation, the entire domain is discretized into 191 finite-difference grid points. The grid spacing in the channel is taken as Δx=0.5m and Δy=0.1m. However, in the vegetated region, the stems are placed at a more refined grid spacing. The Mannings or Chezy's roughness coefficient is an essential parameter for shallow water simulations. These coefficients are generally estimated from the mean diameter of the bed material. Amreeva and Kurbak (2007) The results obtained based on these measures are shown in Table

Conflict of interest: None
Code availability: The code used in this study is developed in MATLAB

Data availability statement
The data that support the finding of this study can be found in https://doi.org/10.1623/hysj.53.4.905 Author contribution:All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by AnupalBaruah. The first draft of the manuscript was written by Anupal Baruah and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.   Exp.no-4.1.1