An integrated numerical model for vegetated surface-saturated subsurface flow interaction
K.S. Erduran
Department of Civil Engineering, Faculty of Engineering & Architecture,
University of Nigde, Campus, 51245, Nigde, Turkey.
Phone: ++903882252288 Fax: ++903882250112
Abstract: This study involves the development of an integrated numerical model to deal with interactions between vegetated surface and saturated subsurface flows. The numerical model is built up by integrating quasi three dimensional (Q3D) vegetated surface model with two dimensional saturated groundwater flow model. The vegetated surface model is constructed by coupling the explicit finite volume solution of the two dimensional shallow water equations with the implicit finite difference solution of Navier stokes equations for vertical velocity distribution. The subsurface model is based on the explicit finite volume solution of two dimensional saturated groundwater flow equations. Ground and vegetated surface water interaction is achieved by the introduction of source-sink terms into the continuity equations. Two solutions are tightly coupled in a single code. The integrated model has been applied to two test cases and the results are satisfactory.
Key Words: Vegetated surface flow, saturated groundwater flow, flow interactions, tight coupling, finite volume method, finite difference method, flow resistance.
Abbreviations
h : water depth (m)
and : depth-averaged velocity components in x and y directions respectively (m/s)
: flow due to infiltration (m/s)
: excess of water coming from the ground (m/s)
g : acceleration due to gravity (m/s²)
and : bed slope in x and y directions respectively
and : friction terms in x and y directions respectively
and : average drag forces in x and y directions respectively (N/m²).
, and :velocity components in x, y and z directions respectively (m/s)
: density of water (kg/m³)
: vertical shear stresses in x and y directions respectively (N/m²)
: drag forces in x and y directions respectively (N/m²)
: vertical eddy viscosities along x and y directions respectively (m²/s)
: density of vegetation per m²
: drag coefficient (empirical constant)
d : diameter of a reed (m)
: effective reed height (m)
: vertical space step(m)
: specific yield (constant)
H : groundwater head (m)
: hydraulic conductivities in x, y and z directions respectively (m/s)
E : thickness of fully saturated groundwater inside the aquifer (m)
Z : vertical elevation from the bottom (m)
Dt : time step (s)
f : groundwater fluxes normal to the cell interfaces (m²/s)
L : length of a cell side (m)
A : area of the cell (m²)
: updated water depth (m)
: updated groundwater head (m)
: updated depth averaged velocities in x and y directions respectively (m/s)
: updated velocity components in x, y and z directions respectively (m/s)
: next time step values of water depth (m)
: next time step values of groundwater head (m)
: next time step values of depth averaged velocities in x and y directions
respectively (m/s)
: next time step values of velocity components in x, y and z directions
respectively (m/s)
1. INTRODUCTION
With the advent of fast digital computers and numerical methods, it has become possible to solve numerically many engineering problems. In this study, an integrated numerical model, which is constructed for simulation of vegetated surface - saturated subsurface flow interactions, is introduced. Surface and subsurface flow interactions are common in nature and these flow processes are the core of the hydrological cycle. The most frequent example takes place between river and aquifer and it often occurs between overland flow and aquifers. These flow processes get more complicated if the surface flow is obstructed by natural causes or manmade structures. This study concentrates on numerical modeling of such flow processes that occur in areas, where the surface flow is often covered and obstructed by vegetation and there exists interaction between the surface and subsurface flow.
Many integrated groundwater and surface water models have been developed. These models can be distinguished according to their spatial and time dimensions and the type of equations and solutions methods used. These factors help us to understand the complexity of the model. In nature, flow is truly three-dimensional. However, dimensional reduction as well as other simplifications on the governing equations and their solutions are made. Hence, surface water flow processes are generally described by 2D shallow water equations or 1D de Saint-Venant equations. In the case of vegetated flow, such simplifications are limited as the flow velocity generally shows considerable changes in all three spatial directions and also experimental studies prove that the flow shows turbulent characteristic [1]. The main impact of vegetation on flow is that it causes a drag resulting in momentum losses [2]. Vegetative characteristics such as height, diameter, placement and stiffness of the vegetation play important role on the flow components [3, 4, 5]. Flow is also controlled by the condition of the vegetation, i.e. whether the vegetation is submerged or non-submerged. Palmer [6] classifies three flow conditions low flows (vegetation is emergent and no bending occurs), intermediate flows (vegetation is completely submerged and bent and flow resistance shows rapid changes with changes in discharge) and high flows (vegetation is forced to a prone or nearly prone position). As stated by Carollo et al. [7] the most design problems, corresponding to vegetated flows fall within intermediate flows and low flow type. Wu et al. [8] also show that vegetation can cause not only a drag effect but also a blockage effect. All these parameters either reduce or increase the drag effects induced by vegetation; hence, they result in a change in flow components. In addition, vegetation plays significant role on dispersion [9]. These parameters also indicate that the computation of flow through vegetation is not an easy task. Many attempts, most of which are based on experiments, have been made to find an accurate way to compute flow resistance induced by vegetation [1, 4, 7, 10 - 12]. The earliest approach to compute uniform flow through rigid vegetation is to use Manning’s formula. Petryk and Bosmajian III [13] introduce a rather simple approach for one-dimensional steady uniform flow including drag forces caused by vegetation. A more mechanistic method for the computation of flow through flexible vegetation is proposed by Kouwen [14]. Flow through submerged and non-submerged vegetation is studied by Fischenich [2] and Wu et al. [4]. Helmiö [15] studies unsteady one- dimensional flow in a compound channel with vegetated floodplains. The model is later applied to the Rhine River in order to assess the effects of resistance caused by partially vegetated floodplains [16]. The effects of type, placement, and density of vegetation on flow components are experimentally studied by Järvelä [10]. Stone and Shen [17] introduce physically based formulas for computation of flow resistance valid for submerged and emergent rigid vegetation. Among the previously developed numerical approaches the following studies covers more general vegetated flow conditions. Kutija and Hong [18] introduce one-dimensional model, suitable for computation of flow through flexible, rigid, submerged and non-submerged vegetation. Darby [17] develops a 1D flow model, which can be dealt with flexible and rigid vegetation. Vionnet et al. [20] determine the flow resistance as well as eddy viscosity coefficients numerically by so called lateral distribution method for flow through flexible plants. 2D depth averaged model with drag effects is employed for the investigation of the effects of vegetation on flow [21]. Simoes and Wang [22] introduce a quasi three-dimensional (Q3D) turbulence model. Their Q3D model is suitable for simulation of flow through rigid vegetation. Shimizu and Tsujimoto [23] also develop a turbulence model, applicable only for the rigid vegetation. Recently, Fischer-Antze et al. [24] introduce a 3D k- turbulence model, which is suitable for rigid submerged vegetation. Zhang and Su [25] also use 3D LES model for the simulation of flow through rigid, straight and smooth cylindrical vegetation. Erduran and Kutija [3] introduce a Q3D simple turbulence model that is applicable to flexible, rigid, submerged and non-submerged vegetation.
As for the subsurface water flow model, previously developed groundwater models also differ according to the equations and the solutions methods used. Generally speaking, groundwater flow is modelled by solving either Richards equation (saturated, unsaturated and variable-saturated flow can be dealt with), or a rather simple equation, which is based on Darcy law and only suitable for simulation of saturated flow processes [26, 27].
The problem gets more complex when a surface water model needs to be coupled with a
groundwater model. Two types of coupling can be achieved; external (loose) coupling and internal (also known as tight or dynamic) coupling. In external coupling surface and groundwater simulation is done simultaneously (one after another) whereas in internally coupled models coupling is provided within the same time level. The latter requires a single source written for both surface and ground water computations and this coupling is also rather difficult to implement comparing with external coupling [28-30].
The main aim of this study is to construct a numerical model that internally couples surface and subsurface flow solutions as well as that is capable of dealing with flow through flexible or rigid and submerged or non-submerged vegetation. Hence, in this study, much attention is paid to the coupling technique and the model features and limits rather then numerical solution techniques as they are partially presented in details in the previous studies [3, 28, and 31].
In the following sections, the equations and their solutions used in the model are briefly introduced. The developed model is applied to two test cases and the model ability for dealing with vegetated surface-saturated subsurface flow has been demonstrated. Assumptions made to develop the model and the model applicability are discussed. Finally, conclusions are presented.
2. METHODOLOGY
The integrated model is constructed by internally coupling vegetated surface flow solution with the solution of saturated groundwater flow. In order to avoid the repetition, the solutions will not be given in detail but the corresponding references, where the solution in details can be found, will be given.
2.1. Used Equations and Solution Methods
The vegetated surface flow computation is achieved by using two main modules within a program. In the first module, the shallow water equations with drag forces are solved by the finite volume method (FVM). The numerical fluxes are computed by Roe scheme [32] and the upwinded technique [33, 34] is used to deal with the bottom slope. This provides better flux balances in the existence of bottom slope. In this module, the following shallow water equations with drag forces are solved.
(1) (2)
(3)
where h is the water depth, and represent the depth-averaged velocity components in the x and y directions respectively, qi is the infiltration from surface to ground, is the excess of water coming from the ground, g is the acceleration due to gravity, and are the bed slope and friction terms respectively in the x direction and similarly and in the y direction. and are the depth averaged drag forces in the x and y directions due to vegetation.
The reader may refer to the following references [3, 31] for the solution to Equations (1) through (3). In the second module, Navier Stokes equations are solved in the vertical direction using the implicit finite difference method on grids, which lie vertically above the cell centres of the finite volumes in the first module, see Figure 1. The Navier Stokes equations including the drag forces can be given as
(4) (5) (6)
where , and are the velocity components in the x, y and z directions respectively, is the density of water, are the vertical shear stresses in the x and y directions respectively, are the additional drag forces per unit area due to vegetation in the x and y directions respectively. The vertical shear stresses are represented in terms of vertical viscosity and the vertical gradient of horizontal velocities as shown in Equation (7).
, = x, y (7) where are the vertical eddy viscosities along the x and y directions respectively.
For computation of the vertical viscosity values for vegetated flow, Kutija and Hong [18] approach is opted here.
As seen in Equations (4) to (6), the momentum equation in the vertical direction is omitted. Hence, a solution is quasi three-dimensional (Q3D). Drag forces,, in the x and y directions due to vegetation are zero above the vegetation and inside the vegetative watercourse they can be computed as
(8)
where is the density of vegetation, is a drag coefficient, d is diameter of a reed, is effective height of a reed, see Figure 2.
In the solution of equations (4) to (6), different discretisation techniques are used in order to increase the stability as well as ease the computation. The implicit finite difference approximations are employed to the following terms; acceleration, drag forces, and the shear stresses. Remaining terms are treated explicitly to decrease the computational effort. The advective terms are discretisied by so called upwind technique. The horizontal gradients of water depth are approximated using forward difference approximations. Resulting approximated equations produce a system of linear algebraic equations. The number of equations is equal to the number of grid points in the vertical direction. The unknown horizontal velocities are computed using a double-sweep algorithm [35] since the matrix of the system is tri-diagonal and they are computed at all the dicretisation points but the vertical shear stresses are computed halfway between each two grid points.
The depth averaged drag forces given in Equations (2) and (3) can be computed by
, = x, y (9)
Equation (9) shows that the depth averaged drag forces are first computed in every vertical grid point, k in the second module and they are later passed to the first unit. In order to add a solution for flow through flexible vegetation, cantilever beam theory [36] is used.