Assessing methods for estimating roughness coefficient in a vegetated marsh area using Delft3D

The following symbols are used in this paper:

a

the frontal area of vegetation stem per unit volume (m⁻¹);

B

surface width (m);

B^V, B^X

volumetric and cross-section blockage factors, respectively (−);

C

Chezy's coefficient (m^1/2/s);

C_b

Chezy's coefficient for non-vegetated channel (m^1/2/s);

C_D

drag coefficient (−);

C*

coefficient for the shear stress at the interface between vegetated and non-vegetated zones (−);

g

gravity acceleration (m/s²);

H

flow depth (m);

h_v

height of vegetation (m);

m

vegetation density (stem/m²);

n

Manning's coefficient (⁠⁠);

n_o, n₁

the maximum Manning's coefficient for non-vegetated channels (⁠⁠);

n_v n₂

the maximum Manning's coefficient for vegetated channels (⁠⁠);

N

total number of observed WSE;

O_i

ith value of observed WSE;

average value of the observed WSE;

S_i

ith value of the simulated WSE;

U_vo

average velocity in the vegetation zone (m/s);

V

averaged velocity across the flow depth (m/s);

w

vegetation frontal width (m);

z₀

length scale of bed roughness (m);

coefficient (−);

β

vegetation resistance parameter (⁠⁠)

0.4 (for sparse, low density vegetation H > 0.3 m);

1.6 (for moderately dense vegetation H = 0.3 m);

6.4 (for very dense vegetation H < 0.3 m); and

η

turbulence length scale (m);

von Karman constant 0.41 (−);

parameter relating to vegetation resistance;

density of water (kg/m³).

Early researchers (e.g. Ree & Palmer 1949; Cowan 1956; Chow 1959; Petryk & Bosmajian 1975) formulated many empirical relations for Manning's roughness coefficient correlating with flow and vegetation properties. However, these relations are only applicable to limited vegetation types and experimental flow conditions (e.g. flow depth) (Klopstra et al. 1997; Wu et al. 1999; Abood et al. 2006). SED2D model, a depth averaged two-dimensional hydrodynamic model, has included an empirical formula to calculate Manning's n for channels with vegetation. The n value in SED2D model is a function of the maximum n value for channel free of vegetation, vegetation height, flow depth, and the maximum n value of vegetated water (Letter et al. 2011). Through measuring turbulence around vegetation, researchers (e.g. Nepf et al. 1997; Nepf 1999) found bed shear stress affected near-bed turbulence production, while far away from the bed, the turbulence wake generated by stems becomes dominated (Nepf et al. 1997). The vortex from the stem wake extracts energy from the mean flow, and feeds it into turbulence kinetic energy (Nepf 1999). Consequently, turbulence intensity initially increases with stem density, and then reaches a peak at a given vegetation density, but eventually subsides as the density becomes very high (Nepf et al. 1997).

However, to date, there is no consensus on which method is the most feasible for calculating Manning's n in freshwater marshes. Vegetation-induced resistance is highly variable in the field due to complex vegetation types, density, and height. To quantify these varieties, we selected the Delft3D-Flow, a widely used open source three-dimensional hydrodynamic model, to study vegetation impact on tidal flow hydrodynamics.

The Delft3D-FLOW open source program (http://oss.deltares.nl/web/delft3d/source-code) is a three-dimensional (3D) hydrodynamic and sediment transport model capable of simulating unsteady incompressible flow (Deltares 2011). The Reynolds Averaged Navier–Stokes equations, under shallow water (hydrostatic pressure) assumption, were solved using a finite difference scheme on a structured staggered curvilinear grid (Stelling & van Kester 1994). Delft3D-FLOW allows users to choose from two vertical grid systems: σ-grid and Z-grid, and four turbulence closure models: constant eddy viscosity coefficient, algebraic eddy viscosity model, k-L model, and k-ɛ model.

Delft3D-FLOW has a function, Trachytopes, for users to define bed and flow resistance on each sub-grid (Deltares 2011). Three classes are available in the Trachytopes function: area, line, and point class. The area class has three types: the first is a constant coefficient for bed roughness, such as White Colebrook, Chezy, and Manning's coefficients, the second accounts for form resistance resulting from sand dunes, and the third is for calculating roughness coefficient in vegetated channels. In a simulation run, the roughness coefficient often remains a constant with time if choosing the first type of area class, while for the second type, it is determined by dune height, and the third by vegetation properties. The line class of Trachytopes function is used to approximate flow resistance for elements with hedge, bridge piers, and other structures. The point class is used to represent a set of point flow resistance elements, such as groups of individual trees or small plant.

Delft3D-Flow has incorporated vegetation effect through an adjusted bed roughness. For example, an implementation based on Klopstra et al. (1996) was added to the Klopstra et al. (1997) equation for calculating the C value of emergent vegetation (Deltares 2011) (Table 2). Additionally, an artificial term was added to account for extra momentum loss due to vegetation using Baptist's (2005) equation (Table 2). Furthermore, the momentum equations and the k-ε turbulent closure model in Delft3D-FLOW have been modified to include vegetation-induced momentum loss as well as influences of vegetation on turbulence generation and dissipation.

Temmerman et al. (2005) studied the vegetation impact on flow hydrodynamic and sedimentation processes in a tidal creek within the Paulina salt marsh in the Scheldt estuary, located at the southwest of Netherland, by modifying the momentum and k-ɛ equations in Delft3D-FLOW. Lately, Horstman et al. (2013) used Delft3D-FLOW for modeling the tidal dynamics in the mangrove forest in Trang Province, Thailand. Both cases are in saltwater marshes where tidal flow is dominant, and vegetation types are distinct from ones in freshwater marshes. None of them have compared the model's performances with results using other methods not embedded in Deflt3D to improve their modelling results.

The objective of this study is to evaluate the accuracy of different methods for calculating vegetation-induced roughness. The methods are not limited to the ones programmed in Delft3D, but all the methods available in literature. Therefore, the first and third types of area class in the Trachytopes function were adopted for incorporating other methods to calculate vegetation-induced roughness. Besides, we selected the orthogonal curvilinear coordinate system, and the vertical plane is σ-grid. The k-ɛ model was chosen to determine the coefficient of eddy viscosity. The sensitivity of modeling results to the selection of vertical grid and key parameters (e.g. a value) was analyzed.

The computational grid of the study area was constructed using available geometric and bed elevation data. The bathymetry of the study area was obtained from USACE, and the elevation is referenced to North American Vertical Datum of 1988 (NAVD88). The computational mesh is a structured grid having 299 points in the main flow direction (M), and 53 points in the direction normal to main flow (N), and 15 layers in the vertical (Figure 4). A finer mesh was used in the inlet canal to capture detailed bathymetry because of rapidly varied bed elevation in this region. The maximum grid sizes in main flow and normal to main flow directions are 177.84 and 462.5 m, respectively, and the minimum are 7.13 and 19.19 m, respectively. The grid aspect ratio (N-grid size/M-grid size) ranges from 1.0 to 20.55. The regions of high aspect ratio are located along the border of the study area, where flow is predominately along the N-direction. Because the velocity gradient close to the bed is very high, the distance between σ-layers is smallest at the bottom, and increases toward free surface. The water depth has reached 6.5 m in the inlet canal, and 2.5 m in the other parts during the flood period from 1:00 pm (CDT) on November 30 to 7:00 pm (CDT) on December 3, 2003.

Eleven stream gauges are available in the study area (Figure 2). The time step is one minute for both approaches to achieve numerically stable solutions. Table 3 is a summary of vegetation roughness parameters for each subarea using the roughness equations in SED2D model (McAlpin et al. 2008). In the equations proposed by Luhar & Nepf (2013), Klopstra et al. (1997) and Baptist (2005), the average frontal width of an individual plant is calculated as 2.25 times the leaf width (approximately equal to 1.25 cm) plus the mean stem diameter (d = 0.4 cm). Since the calculated frontal width, symbolized by w, is 3.2 cm, equal to 8d, this study assumes 8d as the vegetation frontal width. Since the α value is the product of vegetation density and its frontal width, the larger the a value, the larger the vegetation-induced resistance. In a typical Panicum hemitomon grass dominated freshwater marsh in Louisiana, vegetation density is 255 stem/m², and the average stem diameter is 0.4 cm, thus the a-value is 8.160 m^–1. This value is within the range, 1–10 m^–1, suggested by Luhar et al. (2008), Lightbody & Nepf (2006), and Leonard & Luther (1995) for marsh grasses. The drag coefficient for vegetated channel, C_D, is equal to 1.0. Shear stress at the interface between vegetated and unvegetated regions is quantified by another drag coefficient, denoted as C_*, and is set as 0.10 (Luhar & Nepf 2013).

Table 3

Roughness parameter values for all sub-areas used in SED2D model equation

Sub-area .	no .	hv .	nv .	.
1	0.06	0.3048	0.10	0.05
2	0.06	0.381	0.40	0.20
3	0.06	0.3048	0.40	0.20
4	0.12	0.6096	0.82	0.41
5	0.12	0.762	0.82	0.41
6	0.06	0.381	0.40	0.20
7	0.12	0.6856	0.82	0.41

Sub-area .	no .	hv .	nv .	.
1	0.06	0.3048	0.10	0.05
2	0.06	0.381	0.40	0.20
3	0.06	0.3048	0.40	0.20
4	0.12	0.6096	0.82	0.41
5	0.12	0.762	0.82	0.41
6	0.06	0.381	0.40	0.20
7	0.12	0.6856	0.82	0.41

View Large

In the first approach, Figure 8 showed that the range of RMSE for the simulated WSE using the constant Manning roughness for unvegetated channels (n₁) is smaller than the results using Manning's roughness for vegetated channels (n₂). The NSE value is also closer to 1.0 for the results using n₁ (Figure 9). It can be seen from Figure 7 that WSE were underestimated when using the constant n₁, while they were overestimated when using the constant n₂. Although the results using n₁ are slightly better than the ones using n₂, none of them matched the results well at all gauges. Therefore, several n values between n₁ and n₂ values were selected to re-run the model. However, for each value, the WSEs matched the observation data at one gauge in one time interval but not at other times. Also, WSEs matched the observation data better at some gauges but got worse at the rest. Regardless of the n value being used, there is no general improvement of modeling results. This excludes the feasibility of using a constant roughness for simulating flow hydrodynamics in a freshwater marsh.

The simulation was run on a PC with Intel (R) Xeon (R) CPU X5550 (two processors) and 32 GB RAM. The CPU times using the first approach were 74,218.84 and 73,815.20 s for the simulation using constant n₁ and n₂, respectively, whereas they are 72,790.89 and 77,252.20 s using Reed et al. (1995) and Klopstra et al. (1997) equations, respectively.

Using the equations in Luhar & Nepf (2013), Klopstra et al. (1997) and Baptist (2005), results showed the best WSE matches with the observation. However, the results can be sensitive to key parameters (e.g. a value). To quantify the sensitivity of modeling results to the variation of a values, three a values were selected based on the average frontal width of an individual grass, equal to d, 11d, and 22d. Since the average vegetation density is 255 stem/m² and the stem diameter is d = 0.4 cm, the calculated a values are 1.020, 11.220, and 22.440 m⁻¹, respectively.

Despite significant efforts in quantifying vegetation resistance in estuarine marshes, there is no consensus on which equation is the best. Because of the varieties of vegetation, their complex stem and leaf structures, and their elastic properties in flows, vegetation-induced flow resistance still needs to be determined by using empirical relations in Table 2 or others in the literature. Although this study showed the simulated WSEs best matched the observations using the equations in Luhar & Nepf (2013), Klopstra et al. (1997) and Baptist (2005), this conclusion may not be valid for another marsh because vegetation species and ages can alter vegetation's mechanic properties, and thus their influences on flow resistance. Nevertheless, this study found the product of vegetation density and its frontal width is a reliable quantitative measure of vegetation property. Vegetation resistance is dependent on its submergence, and therefore varies with flow depth. Engineers need to validate a selected vegetation resistance equation using observed data of WSE, flow depth, and velocity before adopting it to any field study.

This study applied the Delft3D-FLOW model to evaluate the accuracy of various methods for calculating vegetation-induced resistance in a maidencane marsh in the Mississippi River delta. Besides two equations embedded in Delft3D model, four other equations were programmed into the Delft3D-FLOW model. The simulated results of WSE were compared with the corresponding observation at eleven USGS gauges. Results showed that vegetation-induced roughness is a spatial and temporal variable that changes with submergence and vegetation blockage index. If treating flow roughness in each sub-area as constant, the results of WSE deviated considerably from the measurements. On the other hand, when using the equations by Luhar & Nepf (2013), Baptist (2005), and Klopstra et al. (1997), reasonable matches with the observed WSE were obtained because these equations have taken the degree of submergence and the vegetation blockage index into account. The best results of WSE were obtained by using a constant vegetation density for the Panicum hemitomon vegetation type, and the resulting a values ranged from 8.160 to 11.220 m^–1. However, this does not warrant the universal applicability of these three equations. Further research is needed to better quantify vegetation property and understand their interactions with flow field.

This project is partially funded by an Iraq Governmental Graduate Student Training Scholarship. In-kind technical and data support from the Coastal Hydraulics Laboratory in Engineering Research and Development Center, the Corps of Engineers, is greatly appreciated.