**Unsteady 1D Flow Model of Natural Rivers with Vegetated Floodplain – An Application to Analysis of Influence of Land Use on Flood Wave Propagation in the Lower Biebrza Basin**

Dorota Miroslaw-Swiatek *Warsaw University of Live Sciences - SGGW Poland* 

#### **1. Introduction**

144 Water Resources Management and Modeling

Ortiz-Pérez, M.A. and D.L. Lanza-Espino, 2006. *Diferenciación del espacio costero de México: Un* 

Ocampo-Torres, F.J. 1980. Análisis y predicción de velocidad mediante un modelo

Stigebrandt, A. 1977. On the effect of barotropic current fluctuations on the two-layer transport capacity of a constriction. Journal of Physical Oceanography 7, 118-122.

UNAM, México.137p.

California BSc 90 pp.

*inventario regional*. México: Serie deTextos Universitarios, Instituto de Geografía-

unidimesional en Bahía San Quintín, B.C, México. Universidad Autónoma de Baja

The present studies on renaturalization of river valleys, river beds as well as quantitative estimation of water demand of protected hydrogenic habitats exposed to flooding, in many cases require application of flooding flow models. The state of fluviogenic ecosystems is dependent mostly on the conditions of their hydrological alimentation by flooding waters. The possibility of particular swamp vegetation occurrence on those areas is related to the presence of annual floods of certain duration and depth. So there is a significant relationship between the plant communities' pattern and persisting hydrological conditions. The most important hydrological characteristics, conditioning of the growth and the development of swamp vegetation, are, first of all (Oswit, 1991; Hooijer, 1996; Zalewski et al., 1997): inundation surface, mean inundation depth, the frequency and the duration of inundation. That indicates that there is a very strong relationship between a vegetation structure and water conditions. On the other hand, floodplain vegetation significantly affects flood extent in the valley. Hydrodynamic models coupled with GIS techniques make it possible to obtain necessary data for the determination of abovementioned hydrological characteristics. They are also tools which facilitate estimation of the influence of different river valley management methods on hydraulic conditions of water flow. They can be applied in the process of decision-making for: projects, investments and operational activities in the field of flood protection and environmental impact assessment.

The valley of the Biebrza River, in particular the Lower Basin of Biebrza, is an exceptionally valuable natural object with inundation of surface waters occurring annually as a vital siteshaping factor – this builds specific habitat conditions, crucial for the unique character of that area. This area is an extremely valuable wetland site of global significance protected by Ramsar Convention and annual flooding influences formation of a unique character of the study site. A zonal system of various plants reflects water conditions of the Lower Biebrza Basin (Oswit, 1991) and causes spatial division of the resistance to water flow (Swiatek et al., 2008).

The hydrodynamic model, providing a correct description of flooding water flow in the area of the Biebrza Valley, could be used as research tool for executing effective policy of natural values protection within Biebrza National Park. The influence of the different floodplain land use, described by changes in vegetation structure, on flood wave propagation in the Lower Biebrza River Basin (LBRB) was analyzed as an example of developed model applying.

Mathematical modelling of river flow with water levels not exceeding a bank elevation is widespread and described by many authors (Cunge et al., 1980, Szymkiewicz, 2010). Applications of hydraulic models, particularly for lowland river valley flooding with water overflowing the main channel, flooding adjacent areas, and flowing in the floodplains covered with various vegetation, are not very common. The models used usually apply simplified approaches to the influence of land use in river valleys on flow conditions. Flow resistance, in rivers as well as in floodplains covered with vegetation, is usually described with generalized Manning's coefficient values subject to plants vegetative alterations. In another simplified approach widely used, the complexity of water flow on floodplains, which is very hard to be sufficiently represented, is considered in the form of so-called 'dead' areas of the cross-section (storage zones), which do not conserve momentum, but take part only in storing water in the river valley. In this approach, floodplain geometry is accounted for in only one of St Venant's equations – a continuity equation – and the momentum equation reduces that to hydraulic parameters within the main channel geometry (Cunge et al., 1980).

The advantage of one-dimensional models is their computational effectiveness for largescale areas. These models provide fairly good results in those cases, when the assumption of one-dimensional flow is valid, that is for river valleys with no extensive floodplains (Horrit & Bates, 2002). On account of their applicational usefulness, these models are constantly developed and upgraded.

The developed model, as opposed to many existing commercial models, withdraws from a simplified description of flow resistance expressed by the spatially differentiated Manning's coefficient and the use of the Darcy-Weisbach relationship. It also enables introduction of water mass and momentum exchange process between the main channel and floodplains, and parts of a cross section covered with vegetation and those with no vegetation. To this end, additional flow resistance along imaginary vertical boundaries between the main channel and floodplains has been introduced (Pasche & Rouve, 1985). Thus the developed model enables accounting for unsteady flow calculations flow resistance resulting from both vegetation covering a cross section and momentum exchange between the main channel and floodplains, proposed in the Pasche approach (Pasche ,1984). The 1D unsteady river flow model with vegetated floodplain was proposed by Helmio (2002). In this model the flow resistance for the riverbed and the floodplains was calculated according to the method by Nuding (1991). That method is poorly documented in literature, especially in respect of the determination of applied parameters. The method of Pasche was also applied in SPRUNNER model (Laks & Kałuża, 2007), where it is used for the determination of active zone of flow, extended by the interaction of the main channel with floodplains.

In this paper an unsteady 1D flow model was developed for a river with vegetated floodplains. The basic form of the non-linear St Venant equations combined with retention effects of the vegetated areas on flood wave conveyance were used in the model. In this approach friction caused by vegetation and additional resistance caused by interaction between the main channel and vegetated areas (the Pasche method) were taken into account. The model developed in this paper, on the contrary to the other models mentioned above, enables also considering of the flow resistance, resulting from any configuration of vegetation cover in a compound channel, in the calculations of unsteady water flow. The model was applied to the Lower Biebrza River, situated in the north-eastern part of Poland, flowing through the last extensive, fairly undisturbed river-marginal peatland in Europe. The roughness height coefficients of vegetation (reeds, sedges, grasses, bushes as well as trees and scrubs) in floodplains were determined by field monitoring. Water level observations collected by an automatic monitoring system were used in calibration and validation process The model was calibrated and validated for flood events in 2006 and 2007 year. Elaborated model will be used as a tool for assessment of various agricultural practices with regard to the effective management for protection of the unique wetland site of the Biebrza National Park

#### **2. Numerical model development**

#### **2.1 The unsteady flow model**

146 Water Resources Management and Modeling

values protection within Biebrza National Park. The influence of the different floodplain land use, described by changes in vegetation structure, on flood wave propagation in the Lower Biebrza River Basin (LBRB) was analyzed as an example of developed model

Mathematical modelling of river flow with water levels not exceeding a bank elevation is widespread and described by many authors (Cunge et al., 1980, Szymkiewicz, 2010). Applications of hydraulic models, particularly for lowland river valley flooding with water overflowing the main channel, flooding adjacent areas, and flowing in the floodplains covered with various vegetation, are not very common. The models used usually apply simplified approaches to the influence of land use in river valleys on flow conditions. Flow resistance, in rivers as well as in floodplains covered with vegetation, is usually described with generalized Manning's coefficient values subject to plants vegetative alterations. In another simplified approach widely used, the complexity of water flow on floodplains, which is very hard to be sufficiently represented, is considered in the form of so-called 'dead' areas of the cross-section (storage zones), which do not conserve momentum, but take part only in storing water in the river valley. In this approach, floodplain geometry is accounted for in only one of St Venant's equations – a continuity equation – and the momentum equation reduces that to hydraulic parameters within the main channel

The advantage of one-dimensional models is their computational effectiveness for largescale areas. These models provide fairly good results in those cases, when the assumption of one-dimensional flow is valid, that is for river valleys with no extensive floodplains (Horrit & Bates, 2002). On account of their applicational usefulness, these models are constantly

The developed model, as opposed to many existing commercial models, withdraws from a simplified description of flow resistance expressed by the spatially differentiated Manning's coefficient and the use of the Darcy-Weisbach relationship. It also enables introduction of water mass and momentum exchange process between the main channel and floodplains, and parts of a cross section covered with vegetation and those with no vegetation. To this end, additional flow resistance along imaginary vertical boundaries between the main channel and floodplains has been introduced (Pasche & Rouve, 1985). Thus the developed model enables accounting for unsteady flow calculations flow resistance resulting from both vegetation covering a cross section and momentum exchange between the main channel and floodplains, proposed in the Pasche approach (Pasche ,1984). The 1D unsteady river flow model with vegetated floodplain was proposed by Helmio (2002). In this model the flow resistance for the riverbed and the floodplains was calculated according to the method by Nuding (1991). That method is poorly documented in literature, especially in respect of the determination of applied parameters. The method of Pasche was also applied in SPRUNNER model (Laks & Kałuża, 2007), where it is used for the determination of active

zone of flow, extended by the interaction of the main channel with floodplains.

In this paper an unsteady 1D flow model was developed for a river with vegetated floodplains. The basic form of the non-linear St Venant equations combined with retention effects of the vegetated areas on flood wave conveyance were used in the model. In this approach friction caused by vegetation and additional resistance caused by interaction

applying.

geometry (Cunge et al., 1980).

developed and upgraded.

In practice unsteady flow in natural rivers is usually treated as one-dimensional flow, and is based on St Venant equations. St Venant equations consist of the dynamic equation and the continuity equation. When water discharge and water level are dependent variables, these equations are written respectively as:

$$\frac{\partial \mathbb{Q}}{\partial t} + \frac{\partial}{\partial \mathbf{x}} \left( \frac{\partial \mathbb{Q}^2}{A} \right) + gA \frac{\partial \mathbb{h}}{\partial \mathbf{x}} + A \frac{|\mathbb{Q}| \mathbb{Q}}{K^2} + q \frac{\mathbb{Q}}{A} = \mathbf{0},\tag{1a}$$

$$
\partial \frac{\partial \mathbf{h}}{\partial t} + \frac{\partial \mathbf{Q}}{\partial \mathbf{x}} = q \tag{1b}
$$

where: *Q* = discharge [m3/s]; *h* = water level [m]; *x* = distance [m]; *t* = time [s]; *A* = cross area of flow [m2]; *B* = width of water surface [m]; *K* = conveyance factor [m3/s]; g = gravitational acceleration [m/s2]; *q* = lateral inflow [m3/s/m]; = momentum correction factor.

The conveyance factor *K* is expressed as

$$K = A \left(\frac{8\,\mathrm{g}\,\mathrm{R}}{\mathrm{A}}\right)^{1/2} \tag{2}$$

where: *R* = hydraulic radius [m]; = average dimensionless friction factor.

The equation (1) requires also determining of boundary conditions as well as initial conditions. Boundary conditions refer to hydraulic properties at the upstream and downstream ends of a river. The model developed in this paper enables only a subcritical flow description. The upstream boundary condition is determined as a discharge hydrograph *Q(t)*. A water level hydrograph *h(t)*, discharge hydrograph *Q(t)*, rating curve *Q(h)* or friction slope *Sf* can be used as the downstream boundary condition.

The initial condition refers to the state of flow in the river when the simulation starts. In the model presented in this paper a steady flow in the channel is used as the initial condition (Swiatek et al., 2006). The finite element method in the Galerkin formulation was used to solve a pair of Eqs. (1) (Szymkiewicz, 2010). This approach leads to the non-linear ordinary differential equations, in which the time-weighted finite difference method is used in the approximation of the time derivative. This method forms a system of algebraic non-linear equations which are numerically solved by iteration method.

The total conveyance (Eq. (2)) for a compound cross-section is obtained by summing the subdivision conveyances of the channel and floodplains:

$$\mathbf{K} = \mathbf{K}\_{\rm ff} + \mathbf{K}\_{\rm ch} + \mathbf{K}\_{\rm rf} \tag{3}$$

where: *K* = total conveyance; *Klf* = total conveyance of the left floodplain; *Kch* = conveyances of the channel; *Krf* = total conveyance of the right floodplain.

The total floodplains conveyances are calculated according to the vegetation distribution with the Pasches's method used for computing the total Darcy-Weisbach friction factor . The total conveyance *K* was introduced to the St Vetnant Eq. (1) computation. It was calculated for each cross-section and water level in the iterative method of solving a system of algebraic non-linear equations (Swiatek, 2008).

#### **2.2 Determination of friction factors**

Hydraulic calculations of the flow in natural rivers with floodplain require methods which include natural vegetation structure of the river waterside zones and the floodplain. The values of the resistance coefficients of the floodplain vegetation still belong to rareness. Presently the most often used are roughness coefficients for Manning's equation, but the choice of the sufficient one among the tabular values is subjective. Another method is determination of plant characteristics elaborated by Rouvè (1987). The division on high, medium and low vegetation (proposed by Bretschneider and Schulz 1985) is used in such calculations. High vegetation is considered here as higher than water flow depth (trees and shrubs) and in small degree go under hydrodynamic water pressure, medium vegetation as approximately equal to water depth (mostly shrubs) and low vegetation which refers mostly to sedge and grass communities (Fig.1). Established criterion is ambiguous and in fact the same vegetation can be ranked into different types in view of the natural water levels variability.

The basis of hydraulic calculations of river flow including its natural vegetation structure of high and medium vegetation is assumption that water flow resistances are the same as resistances which occur when water overflow regularly distributed vegetation with averaged geometric parameters (DVWK, 1991; Kubrak & Nachlik, 2003).

Parameters which describe vegetation of the floodplain and are used in calculations are an average tree diameter or shrubs branches *dp* and distances between them in the direction of the water flow *ax* and transversal to it *ay* (Fig.2). Named parameters are determined on the basis of field measurements in the area of the water flow.

The resistance of flow caused by the roughness of low vegetation occurring on the scarp, the channel bed and floodplain is calculated from the formula given by Colebrook-White:

$$\frac{1}{\sqrt{\lambda\_{\rm s}}} = -2.03 \log \left( \frac{2.51}{\text{Re}\sqrt{\lambda\_{\rm s}}} + \frac{k\_s}{14.84 \text{ R}} \right) \tag{4}$$

solve a pair of Eqs. (1) (Szymkiewicz, 2010). This approach leads to the non-linear ordinary differential equations, in which the time-weighted finite difference method is used in the approximation of the time derivative. This method forms a system of algebraic non-linear

The total conveyance (Eq. (2)) for a compound cross-section is obtained by summing the

where: *K* = total conveyance; *Klf* = total conveyance of the left floodplain; *Kch* = conveyances

The total floodplains conveyances are calculated according to the vegetation distribution with the Pasches's method used for computing the total Darcy-Weisbach friction factor

The total conveyance *K* was introduced to the St Vetnant Eq. (1) computation. It was calculated for each cross-section and water level in the iterative method of solving a system

Hydraulic calculations of the flow in natural rivers with floodplain require methods which include natural vegetation structure of the river waterside zones and the floodplain. The values of the resistance coefficients of the floodplain vegetation still belong to rareness. Presently the most often used are roughness coefficients for Manning's equation, but the choice of the sufficient one among the tabular values is subjective. Another method is determination of plant characteristics elaborated by Rouvè (1987). The division on high, medium and low vegetation (proposed by Bretschneider and Schulz 1985) is used in such calculations. High vegetation is considered here as higher than water flow depth (trees and shrubs) and in small degree go under hydrodynamic water pressure, medium vegetation as approximately equal to water depth (mostly shrubs) and low vegetation which refers mostly to sedge and grass communities (Fig.1). Established criterion is ambiguous and in fact the same vegetation can be ranked into different types in view of the natural water levels variability.

The basis of hydraulic calculations of river flow including its natural vegetation structure of high and medium vegetation is assumption that water flow resistances are the same as resistances which occur when water overflow regularly distributed vegetation with

Parameters which describe vegetation of the floodplain and are used in calculations are an average tree diameter or shrubs branches *dp* and distances between them in the direction of the water flow *ax* and transversal to it *ay* (Fig.2). Named parameters are determined on the

The resistance of flow caused by the roughness of low vegetation occurring on the scarp, the channel bed and floodplain is calculated from the formula given by Colebrook-White:

1 2.51 2.03log Re 14.84

 

*s s*

*s*

(4)

*k*

*R*

averaged geometric parameters (DVWK, 1991; Kubrak & Nachlik, 2003).

basis of field measurements in the area of the water flow.

*KK K K lf ch rf* (3)

.

equations which are numerically solved by iteration method.

of the channel; *Krf* = total conveyance of the right floodplain.

of algebraic non-linear equations (Swiatek, 2008).

**2.2 Determination of friction factors** 

subdivision conveyances of the channel and floodplains:

Fig. 1. Vegetation classification proposed by Bretschneider and Schulz (1985)

Fig. 2. The high vegetation geometric characteristics

where: *<sup>s</sup>* = friction factor of low vegetation or not overgrown part of the cross- section [-]; *ks* = roughness height of the overgrown (low vegetation) or not overgrown cross-section [m]; *Re* = the Reynolds number for part of the cross-section.

As it results from the Colebrook-White law, the friction factors of flow depend on the Reynolds number and on relative roughness *ks/R*. The influence of the Reynolds number on friction factors decreases along the increase of its value and of relative roughness of channel sides. In natural channels the influence of the Reynolds number for values higher than 25000 may be neglected without any harm to the precision of calculations. Therefore Rickert (1988) recommends application of the equation (4) for practical calculations, in the following form:

$$\frac{1}{\sqrt{\lambda\_s}} = -2.03 \log \left( \frac{k\_s}{14.84 \ R} \right) \tag{5}$$

The friction factor *<sup>v</sup>* for trees and bushes (submerged part of high vegetation) is calculated from the following expression (Pasche & Rouve, 1985):

$$\mathcal{A}\_{\upsilon} = \frac{4h\_z d\_p}{a\_x a\_y} \mathcal{C}\_{\text{\tiny{\text{\tiny{\upsilon}}}}} \tag{6}$$

where: *hz* = height of submerged part of trees [m] (Fig.3); *dp* = trees diameter [m]; *ax*, *ay* = distance between plants along the flow and transversal; C*WR* = dimensionless drag coefficient for submerged part of the trees or bushes.

The drag coefficient *CWR* depends on the ratio of the *Vi* flowing velocity to the average velocity *Vv* of the flow going through tree overgrown areas (relative velocity *Vr*), and is described by the following empirical formula (Rickert, 1986):

$$\mathcal{C}\_{\mathsf{VMR}} = \left( 1.1 + 2.3 \frac{d\_p}{a\_y} \right) \left( V\_r \right)^2 + 2 \left( \frac{1}{1 - d\_p \,/\ a\_y} - 1 \right) \tag{7}$$

$$\left(V\_r^2 = \left(\frac{V\_i}{V\_v}\right)^2 = 0.6 + 0.5\log\left(\frac{a\_x}{a\_y}\right)\tag{8}$$

The expression (7) is valid when 0.05<dp/ay<3 and 0.2<ax/ay<2.

The resistance of flow in parts of cross-sections overgrown by high vegetation depends on both vegetation and bed roughness. The friction factor for this area, according to the concept issued by Einstein and Banks (Indlekofer, 1981), is the following sum:

$$
\lambda = \lambda\_s + \lambda\_v \tag{9}
$$

where: = average friction factor in part of the cross-section [-]; *<sup>s</sup>* = friction factor caused by channel bed roughness or low vegetation[-]; *v* = friction factor for non-submerged and nonflexible vegetation (high vegetation, Fig.3.) [-].

Physically, composite roughness along the wetted perimeter of the compound cross section modifies velocity distribution in the cross section (Fig.3). A detailed examination of the effects of varying wall roughness and cross sectional geometry would require a threedimensional analysis of the flow. Pasche (1984) proposed one-dimensional analysis of steady flow in a compound cross-section of the lowland river based on the Darcy-Weisbach formula. According to the observed velocity distribution a compound river cross section is divided into sections with vertical imaginary walls with roughness height of channel *kT* and friction factor *<sup>T</sup>* between the main channel and neighboring floodplains (Fig.4). The heights of these boundaries (*hT1*, *hT2*) are taken into consideration in calculations of the wetted perimeter of the main channel. Mean velocity in each section (in Fig.4. signed by symbols *T1* and *T2* for vertical imaginary walls and numbers 1, 3, 2 for scarps and channel bed) for a channel is calculated from the Darcy-Weisbach equation, which results from the momentum balance in part *'i'* of the cross-section:

$$w\_i = \sqrt{\frac{8 \text{g} \mathcal{R}\_i S\_f}{\lambda\_i}} \tag{10}$$

4 *z p v WR x y h d C a a*

where: *hz* = height of submerged part of trees [m] (Fig.3); *dp* = trees diameter [m]; *ax*, *ay* = distance between plants along the flow and transversal; C*WR* = dimensionless drag

The drag coefficient *CWR* depends on the ratio of the *Vi* flowing velocity to the average velocity *Vv* of the flow going through tree overgrown areas (relative velocity *Vr*), and is

*p*

2

<sup>2</sup> 0.6 0.5 log *i x <sup>r</sup> v y*

The resistance of flow in parts of cross-sections overgrown by high vegetation depends on both vegetation and bed roughness. The friction factor for this area, according to the concept

Physically, composite roughness along the wetted perimeter of the compound cross section modifies velocity distribution in the cross section (Fig.3). A detailed examination of the effects of varying wall roughness and cross sectional geometry would require a threedimensional analysis of the flow. Pasche (1984) proposed one-dimensional analysis of steady flow in a compound cross-section of the lowland river based on the Darcy-Weisbach formula. According to the observed velocity distribution a compound river cross section is divided into sections with vertical imaginary walls with roughness height of channel *kT* and

of these boundaries (*hT1*, *hT2*) are taken into consideration in calculations of the wetted perimeter of the main channel. Mean velocity in each section (in Fig.4. signed by symbols *T1* and *T2* for vertical imaginary walls and numbers 1, 3, 2 for scarps and channel bed) for a channel is calculated from the Darcy-Weisbach equation, which results from the momentum

8 *<sup>i</sup> <sup>f</sup>*

*i gR S*

*i*

*v*

*<sup>T</sup>* between the main channel and neighboring floodplains (Fig.4). The heights

*V a <sup>V</sup> V a* 

*d*

*WR r*

*C .. V*

<sup>2</sup> <sup>1</sup> 11 23 2 <sup>1</sup>

1

*y p y*

*a d /a* 

*<sup>v</sup>* for trees and bushes (submerged part of high vegetation) is calculated

(6)

*s v λ λ λ* (9)

*v* = friction factor for non-submerged and non-

*<sup>λ</sup>* (10)

(7)

(8)

*<sup>s</sup>* = friction factor caused by

The friction factor

where: 

friction factor

balance in part *'i'* of the cross-section:

from the following expression (Pasche & Rouve, 1985):

coefficient for submerged part of the trees or bushes.

described by the following empirical formula (Rickert, 1986):

The expression (7) is valid when 0.05<dp/ay<3 and 0.2<ax/ay<2.

channel bed roughness or low vegetation[-];

flexible vegetation (high vegetation, Fig.3.) [-].

issued by Einstein and Banks (Indlekofer, 1981), is the following sum:

= average friction factor in part of the cross-section [-];

where: *Sf* = hydraulic slope [-]; *vi* = average velocity in the 'ith' sub-domain of the main channel; *<sup>i</sup>* = friction factor of the *'ith'* sub-domain.

The average friction factor in the whole main channel *<sup>g</sup>* is calculated with consideration of the friction factors in every sub-section of the main channel (Fig.4):

Fig. 3. Scheme of compound channel with floodplain with high vegetation and velocity distribution

Fig. 4. Scheme of main channel sub-sections with different hydraulic parameters

The friction factor *<sup>i</sup>* is calculated from the formula (5), in which hydraulic radius *Ri* and *RT1*, *RT2* related to different roughness factors in cross-sections is determined according to the concept by Einstein (Kubrak, 2003), assuming the equality of average velocity (4) in the main channel *v* and average velocity *vi* in the every sub-domain *i*:

$$\sqrt{\frac{8\overline{\text{g}R\_iS\_f}}{\lambda\_i}} = \sqrt{\frac{8\overline{\text{g}RS\_f}}{\lambda\_g}}\tag{12}$$

When floodplain is covered by low vegetation, the interface of the floodplain from the main channel is treated as the rough wall with roughness height *kT1*, *kT2* , that have the same values like left and right scarp of the main channel. The calculation of the hydraulic radius for different roughness factors in cross-sections is conducted using the iterative method and basing on the formula (12).

The way of calculation of flow intensity in the floodplain with trees is more sophisticated and based on the theory proposed by Passche and Rouve (1985). According to the observed velocity distribution (Fig.3) the compound channel is divided into four sections. Due to the horizontal mass and momentum exchange between areas with steep velocity gradients the flow is reduced in the main channel (sub section IIIa) and accelerated on the floodplain (sub section II). At calculation of the discharge in the main channel, the interface plane of the floodplain from the main channel is treated as the rough wall with roughness of channel *kT* and friction factor *T* (Fig.3). The friction coefficient *<sup>T</sup>* in interface plane is introduced to determine the flow reduction in the main channel along the width bIII, while the discharge on the floodplain is increased along the width bII.

In reality the resistances to flow by the interface are made by intensive cyclical impulses of mass and momentum in transverse direction to the main flow, as well as by associated high turbulence strains and whirls on the surface of water in round-the-tree flows. The roughness height in the interfaces (Fig.3) is calculated from the formula (Pasche, 1985):

$$k\_T = 0.854 \,\, R\_{,T} \,\, \Omega \left(\frac{2 \,\, b\_{\rm II}}{b\_{\rm ch}}\right)^{1.07} \tag{13}$$

where: *RT* = hydraulic radius concerning the resistance by the interface [m]; *bII* = width of the floodplain zone of influence of the area with trees [m]; *bch* = width of the main channel [m].

The width of the area of floodplain influenced on flow in main channel is calculated from the formula:

$$h\_{\rm II} = \frac{h\_T}{\lambda\_v (0.068 \, e^{0.56 \, C\_T} - 0.056)} \tag{14}$$

where: *hT* = depth of interfaces plane between main channel and floodplain [m]; *cT* = dimensionless velocity in interfaces.

The friction factor *<sup>v</sup>* of the high vegetation (trees or bushes) is calculated from expression (6). The dimensionless velocity in interfaces plane is calculated from the formula:

$$C\_T = -\Im \mathcal{Z} \nabla \log \Omega + \text{2.85} \tag{15}$$

where: = vegetation parameter.

152 Water Resources Management and Modeling

*RT2* related to different roughness factors in cross-sections is determined according to the concept by Einstein (Kubrak, 2003), assuming the equality of average velocity (4) in the main

> 8 8 *<sup>i</sup> f f i g gR S gRS*

When floodplain is covered by low vegetation, the interface of the floodplain from the main channel is treated as the rough wall with roughness height *kT1*, *kT2* , that have the same values like left and right scarp of the main channel. The calculation of the hydraulic radius for different roughness factors in cross-sections is conducted using the iterative method and

The way of calculation of flow intensity in the floodplain with trees is more sophisticated and based on the theory proposed by Passche and Rouve (1985). According to the observed velocity distribution (Fig.3) the compound channel is divided into four sections. Due to the horizontal mass and momentum exchange between areas with steep velocity gradients the flow is reduced in the main channel (sub section IIIa) and accelerated on the floodplain (sub section II). At calculation of the discharge in the main channel, the interface plane of the floodplain from the main channel is treated as the rough wall with roughness of channel *kT*

determine the flow reduction in the main channel along the width bIII, while the discharge

In reality the resistances to flow by the interface are made by intensive cyclical impulses of mass and momentum in transverse direction to the main flow, as well as by associated high turbulence strains and whirls on the surface of water in round-the-tree flows. The roughness

> , <sup>2</sup> 0.854 *II*

where: *RT* = hydraulic radius concerning the resistance by the interface [m]; *bII* = width of the floodplain zone of influence of the area with trees [m]; *bch* = width of the main channel

The width of the area of floodplain influenced on flow in main channel is calculated from

where: *hT* = depth of interfaces plane between main channel and floodplain [m]; *cT* =

0.56 (0.068 0.056) *<sup>T</sup> T*

*<sup>v</sup>* of the high vegetation (trees or bushes) is calculated from expression

*<sup>b</sup> k R*

1.07

*<sup>e</sup>* (14)

*ch*

*b*

*T* (Fig.3). The friction coefficient

height in the interfaces (Fig.3) is calculated from the formula (Pasche, 1985):

*T T*

*II C v <sup>h</sup> <sup>b</sup>* 

(6). The dimensionless velocity in interfaces plane is calculated from the formula:

*<sup>i</sup>* is calculated from the formula (5), in which hydraulic radius *Ri* and *RT1*,

*λ λ* (12)

*<sup>T</sup>* in interface plane is introduced to

(13)

The friction factor

basing on the formula (12).

and friction factor

[m].

the formula:

The friction factor

dimensionless velocity in interfaces.

on the floodplain is increased along the width bII.

channel *v* and average velocity *vi* in the every sub-domain *i*:

The coefficient is calculated from the following form:

$$\mathbf{Q} = \left( \mathbf{0}.07 \frac{a\_{\rm NL}}{a\_{\rm x}} \right)^{3.29} + \left( \frac{a\_{\rm NB}}{a\_{\rm y}} \right)^{0.95} \tag{16}$$

$$a\_{\rm NL} = 128.87 \,\text{C}\_{\rm PV} d\_p \left( 1 + \frac{\text{g } a\_{\rm NL} \, S\_f}{V\_T^2 / 2} \right)^{-2.143} \tag{17}$$

$$a\_{\rm N\beta} = 0.24 \,\, a\_{\rm NL}^{0.59} \left( \mathbf{C}\_{\rm \mathcal{W}} d\_p \right)^{0.41} \tag{18}$$

where: *aNL*, *aNB* = the length and the width of the Karmann path formed at a single plant [m]; *CW* = drag coefficient; *VT* = velocity in the interface plane between the main channel and its floodplain.

The drag coefficient *CW* is determined for a single tree at ideal two-dimensional flow. Its variability for different forms of the turbulent flow was given by Lindner (1982), and in natural channels, where Reynolds number for a single plant varies between 8103-1105 has value 1.2.

The velocity *VT* in the interface plane between the main channel and its floodplain is calculated from the following formula:

$$V\_T = C\_T \sqrt{\frac{\lambda\_v}{8}} \upsilon \tag{19}$$

The velocity *VT* calculated from the equation (19) is usually different from the one assumed at the beginning of *aNL* calculation (17), therefore the whole set of calculations should be repeated, accepting the calculated value *VT* according to formula (19) as the next approximation.

The method described above has the following limitations (Schumacher, 1995):


In the elaborated model, when the first or the second of the above conditions is not valid, in the formula (11) in place of *T1* , *T2* the friction factors of the flow left, right scarp *<sup>i</sup>* is used. The third condition is usually fulfilled for natural lowland rivers with floodplain.

The next problem, which can occur in calculation when Colebrook-White equation (5) is used, is related to the value of relative roughness *ks/R*. When the value of *ks/R* is approximately 14.84, the friction factor *<sup>s</sup>* is indeterminate, and in the case of relative roughness a little bit smaller than this number it has a huge non–physical value. Therefore, to avoid these difficulties the limited value of *s* is introduced in the calculation in the model.

### **3. Application: the lower Biebrza river basin**

#### **3.1 Research area**

The Valley of the Biebrza River is situated in the north-eastern part of Poland (Fig.5). The research focused on the Lower Biebrza Basin (LBB), which is located in the southern part of the Biebrza Valley, from the village of Osowiec to the mouth of Biebrza to the Narew River. The length of the lower basin amounts to 30 km, while its width ranges from 12 to 15 km. The largest surface in the area of that basin is occupied by a flood terrace, which includes vast, flat peatlands and gently folded, loamy near-riverbed zone of the width from 1 to 2 km. The riverbed of Biebrza has a winding course in the lower basin, it forms numerous meanders, side reaches and ox-bows, through which water flows during floods. The length of the riverbed varies between nearly 20m to 35m. The riverbed is significantly emerged from the valley, and in the southern part, close to the mouth to the Narew River, the riverbed is sharply incised in the bottom of the valley. The most important tributaries to the Biebrza River in the lower basin are: the right-side tributaries - the Klimaszewnica and the Wissa Rivers, and the left-side tributary: the Kosodka River. The Valley of the Lower Biebrza has a characteristic, zonal pattern of plant communities, which entirely reflects hydrological conditions dominating in the study area.

Fig. 5. Location of the Lower Biebrza Basin and plant characteristics used in this study

According to the previous investigation by Oświt (1991) five vegetation zones are distinguished: reed communities Phragmition, sedge-reed communities Magnocaricion, sedge-moss communities Scheuzchzerio caricetea fuscae, willow and birch shrubs and also the communities of alder carr and birch-alder forest Alnetea glutinosae. The most valuable ecosystems here are not only natural areas of peatlands but also large areas of open semimeadows, which are the result of extensive agricultural use. The hydrogenic dependent habitats in the Biebrza marshes run on a stable ground water inundation or river flooding which occurs regularly and reinforce these representative water ecosystems every year. Changes of water conditions and discontinuance of extensive agriculture causes transformation of meadows and pastures into tall herb vegetation, reed and in the end the succession of shrubs and forest upon non-forest ecosystems of peatlands which can be observed in some places in the Biebrza River valley.

#### **3.2 Hydraulic and topographical data**

154 Water Resources Management and Modeling

The Valley of the Biebrza River is situated in the north-eastern part of Poland (Fig.5). The research focused on the Lower Biebrza Basin (LBB), which is located in the southern part of the Biebrza Valley, from the village of Osowiec to the mouth of Biebrza to the Narew River. The length of the lower basin amounts to 30 km, while its width ranges from 12 to 15 km. The largest surface in the area of that basin is occupied by a flood terrace, which includes vast, flat peatlands and gently folded, loamy near-riverbed zone of the width from 1 to 2 km. The riverbed of Biebrza has a winding course in the lower basin, it forms numerous meanders, side reaches and ox-bows, through which water flows during floods. The length of the riverbed varies between nearly 20m to 35m. The riverbed is significantly emerged from the valley, and in the southern part, close to the mouth to the Narew River, the riverbed is sharply incised in the bottom of the valley. The most important tributaries to the Biebrza River in the lower basin are: the right-side tributaries - the Klimaszewnica and the Wissa Rivers, and the left-side tributary: the Kosodka River. The Valley of the Lower Biebrza has a characteristic, zonal pattern of plant communities, which entirely reflects

Fig. 5. Location of the Lower Biebrza Basin and plant characteristics used in this study

**3. Application: the lower Biebrza river basin** 

hydrological conditions dominating in the study area.

**3.1 Research area** 

The developed hydrodynamic model was used in river and floodplain flow simulations in the LBB. A river water flow model was elaborated for the reach BD1 –Burzyn, where the geometry of the riverbed and the floodplains were described by 47 cross-sections (Fig. 6). The upstream boundary condition was defined in the form of flow hydrograph for BD1 cross-section as the sum of flows for Osowiec gauge at the Biebrza River and for Przechody gauge at the Rudzki Channel. The downstream boundary condition was defined by the rating curve for Burzyn gauge at the Biebrza River (cross-section BD17).The alimentation by the main tributary – the Wissa River – was described as a point lateral inflow, located at the cross-section BD11, applying the discharge data for Czachy gauge. Available hydrological data for gauges: Osowiec, Przechody, Czachy and Burzyn proves that in the period of annual spring floods the outflow from the study area (recorded at Burzyn gauge) is much higher than the sum of inflows from Osowiec, Przechody and Czachy, which indicates a significant alimentation from a sub-catchment (Fig. 7).

In the elaborated hydrodynamic model, on account of insufficient recognition of hydrological alimentation by ground and snowmelt waters, the inflow from a subcatchment was described with the use of the following, simplified formula:

$$Q\_{lateral}(t) = Q\_{out}(t + \Delta t) - Q\_{in}(t) \tag{20}$$

where: *Qlateral(t)* – sub-catchment discharge at time t [m3/s*]; Qout(t)* – discharge at Burzyn at time *t* [m3/s]; *Qin(t)* – inflow discharge (at BD1 and Czachy gauge) at time *t* [m3/s], *Δt* – time lag [s].

Next, the *Qlateral* inflow was uniformly imposed along the river section. The value of *Δt* was adapted as a result of numerical experiment during a calibration process. The developed solution of sub-catchment alimentation is obviously a rough approximation of the actual situation, as it does not account for the spatial variation of the lateral inflow. A more realistic approach to this issue needs to be elaborated in further study.

The analyzed reach of the Biebrza River has a length of 41.5 km. The Rudzki Channel flows into the Biebrza River about 8 km upstream of the opening cross-section BD1, however, the Wissa River, which is the right-side tributary, inflows to Biebrza at the kilometer 14.5 (Fig. 6). The river bed elevations were measured by manual sounding, the neighbourhood of

Fig. 6. Digital Elevation Model of the LBB with location of the cross – sections used in hydrodynamic model

river banks was leveled and points situated in floodplain were calculated from Digital Elevation Model and topographic maps for each cross section. Close to the river, the cross sections were measured in the field by coupling the traditional leveling survey with Differential GPS (DGPS) techniques. The topography of the floodplains was determined by the analysis of the existing Digital Elevation Model (Swiatek and Chormanski, 2007) (Fig. 6). The DEM was elaborated with the use of ArcGis software. The elevation model represents the shape of the valley in the form of a raster map of 25 m resolution. The main sources of data for DEM development were contour lines digitized on a 1:25 000 scale topographical map and also elevation points measured by GPS technique at transects perpendicular to the river. The elaborated DEM is characterized by mean square error not exceeding 0.35 m. (Chormanski, 2003).

Fig. 7. The comparison of inflow with outflow for Biebrza River in the period 01 March 1999 – 17 June 1999 time period (Q– flow [m3/s])

#### **3.3 Floodplain vegetation data**

156 Water Resources Management and Modeling

#\* *D1-Pale*

**Legend**

#\* divers

[m asl]

**DEM**

cross-sections rivers

*PRZECHODY*

*D2-Klimaszewnica* 

Fig. 6. Digital Elevation Model of the LBB with location of the cross – sections used in

river banks was leveled and points situated in floodplain were calculated from Digital Elevation Model and topographic maps for each cross section. Close to the river, the cross sections were measured in the field by coupling the traditional leveling survey with Differential GPS (DGPS) techniques. The topography of the floodplains was determined by the analysis of the existing Digital Elevation Model (Swiatek and Chormanski, 2007) (Fig. 6). The DEM was elaborated with the use of ArcGis software. The elevation model represents the shape of the valley in the form of a raster map of 25 m resolution. The main sources of data for DEM development were contour lines digitized on a 1:25 000 scale topographical map and also elevation points measured by GPS technique at transects perpendicular to the river. The elaborated DEM is characterized by mean square error not exceeding 0.35 m.

#\*

*WISSA*

*CZACHY*

hydrodynamic model

(Chormanski, 2003).

#\* *D4-Chyliny*

*D3-Kosodka* #\*

The vegetation parameters of Biebrza floodplains used in calculations were determined on the basis of vegetation inventory in the Lower Biebrza Basin, vegetation map and PHARE aerial photos. Every part of a cross-section has one bottom roughness parameter, related to the type of plants overgrowing in this site of a floodplain. An example of vegetation distribution in selected BD3 cross-section is shown in Fig.8.

The methodology of plant characteristics determination proposed by German engineers (DVWK, 1991) was applied in this study for high and medium vegetation in every crosssection established in the Lower Biebrza Basin. Water flow resistance in flowing round of high vegetation such as bushes or trees is calculated on the basis of supplementary average vegetation diameter *dp* and supplementary distances between plants in direction of water flow *ax*, and transversal to it *ay* (equation 6). Resistances of totally flooded vegetation are characterized by roughness height of flexible and stiff vegetation *ks*. It was determined on the basis of values given by Ritterbach (1991) for the Wupper River Valley in Germany and field measurements in LBB of particular vegetation type's height (Szporak et al, 2006). The vegetation types and parameters of Biebrza floodplains are presented in Fig. 5.

#### **3.4 Model calibration and validation**

The main role of the model is the reproduction of flood water levels and discharges in the LBB. The model, which is elaborated in a correct manner, should also provide proper reproduction of water flow not exceeding bankful water. For the sake of this, the hydrodynamic model was also calibrated and validated for unsteady state in flood condition and when water level is less than bankful. The values of roughness height coefficients in the riverbed were determined for each cross-section as a result of numerical experiments, calibrating the model for the steady water table elevation, which was measured in the October 2005 (the discharge at BD1 crosssection recorded 22.10 m3/s, while the lateral inflow from Wissa was equal to 2.15 m3/s).

Mean square error (RMSE) for the steady state reached 0.03 m. The calculations for the steady state were performed in order to identify *ks* coefficients at the computational cross-sections of the main riverbed for non-vegetation period. The calibrated values of roughness height coefficient for the riverbed varied from 0.1 to 0.7 m. The resultant roughness coefficients were next used for model validation for unsteady state in the period from 1 November 2005 to 30 November 2005. For that period the available data considered constant measurements of water stages by automatic sensors (divers) in the riverbed. The sensors were located at four following places along the river: D1-Pale at 48+648 km, D2-Klimaszewnica at 40+483 km, D3- Kosodka at 24+195 km and D4-Chyliny at 17+843 km (Fig. 6). The sampling rate of water levels was equal to 6 hours. Maximum value of the RMSE is 0.06 m in this case (Tab.2) and takes place at diver D1. Results of the executed simulation confirmed that calibrated values of *ks* coefficient can be used in river flow modelling for main channel in non-vegetation period.

The elaborated model was next calibrated for the spring flood period: 1 April 2006 to 30 April 2006, based also on constant measurements of water stages collected by automatic sensors. In the process of model calibration the following parameters were estimated: dimensions of active zones on the floodplain, hydraulic features for low and high vegetation in floodplains, time lag in the *Qlateral* assessment (equation 20). Because the analyzed flood event took place during the non-vegetation period, *ks* coefficients for the riverbed were the same as in the previous simulation.

The dimensions of active zones on the floodplain were determined by flow capacity calculation for each cross-section, which is used in the model. An example of a relationship between the discharge and the width of a cross-section, which takes part in water flow, is shown in Fig.9. Calibrated values of hydraulic parameters of Biebrza floodplain vegetation types are presented in Tab.1. The best correspondence between the observed discharge hydrograph and the one produced by the model at the Burzyn gauge was calculated for a time lag *Δt* of 24 hours.

The methodology of plant characteristics determination proposed by German engineers (DVWK, 1991) was applied in this study for high and medium vegetation in every crosssection established in the Lower Biebrza Basin. Water flow resistance in flowing round of high vegetation such as bushes or trees is calculated on the basis of supplementary average vegetation diameter *dp* and supplementary distances between plants in direction of water flow *ax*, and transversal to it *ay* (equation 6). Resistances of totally flooded vegetation are characterized by roughness height of flexible and stiff vegetation *ks*. It was determined on the basis of values given by Ritterbach (1991) for the Wupper River Valley in Germany and field measurements in LBB of particular vegetation type's height (Szporak et al, 2006). The

The main role of the model is the reproduction of flood water levels and discharges in the LBB. The model, which is elaborated in a correct manner, should also provide proper reproduction of water flow not exceeding bankful water. For the sake of this, the hydrodynamic model was also calibrated and validated for unsteady state in flood condition and when water level is less than bankful. The values of roughness height coefficients in the riverbed were determined for each cross-section as a result of numerical experiments, calibrating the model for the steady water table elevation, which was measured in the October 2005 (the discharge at BD1 crosssection recorded 22.10 m3/s, while the lateral inflow from Wissa was equal to 2.15 m3/s).

Mean square error (RMSE) for the steady state reached 0.03 m. The calculations for the steady state were performed in order to identify *ks* coefficients at the computational cross-sections of the main riverbed for non-vegetation period. The calibrated values of roughness height coefficient for the riverbed varied from 0.1 to 0.7 m. The resultant roughness coefficients were next used for model validation for unsteady state in the period from 1 November 2005 to 30 November 2005. For that period the available data considered constant measurements of water stages by automatic sensors (divers) in the riverbed. The sensors were located at four following places along the river: D1-Pale at 48+648 km, D2-Klimaszewnica at 40+483 km, D3- Kosodka at 24+195 km and D4-Chyliny at 17+843 km (Fig. 6). The sampling rate of water levels was equal to 6 hours. Maximum value of the RMSE is 0.06 m in this case (Tab.2) and takes place at diver D1. Results of the executed simulation confirmed that calibrated values of *ks* coefficient can be used in river flow modelling for main channel in non-vegetation period. The elaborated model was next calibrated for the spring flood period: 1 April 2006 to 30 April 2006, based also on constant measurements of water stages collected by automatic sensors. In the process of model calibration the following parameters were estimated: dimensions of active zones on the floodplain, hydraulic features for low and high vegetation in floodplains, time lag in the *Qlateral* assessment (equation 20). Because the analyzed flood event took place during the non-vegetation period, *ks* coefficients for the riverbed were the

The dimensions of active zones on the floodplain were determined by flow capacity calculation for each cross-section, which is used in the model. An example of a relationship between the discharge and the width of a cross-section, which takes part in water flow, is shown in Fig.9. Calibrated values of hydraulic parameters of Biebrza floodplain vegetation types are presented in Tab.1. The best correspondence between the observed discharge hydrograph and the one produced by the model at the Burzyn gauge was calculated for a time lag *Δt* of 24 hours.

vegetation types and parameters of Biebrza floodplains are presented in Fig. 5.

**3.4 Model calibration and validation** 

same as in the previous simulation.


Table 1. The hydraulic plant characteristics calibrated for various vegetation types in the LBB (non-vegetation period)

Fig. 9. Relationship between the river discharge and the active part of the floodplain (BD3 cross-section)

Calculation results for the flood in 2006 were compared with the water level observations at four automatic sensors (Fig. 10-13). The RMSE water stage varies from 0.07 m to 0.14 m (Tab.2) and indicates that the developed model correctly reproduces the range of variability of water levels for flood simulation. On the other hand, values of correlation coefficient (R) (Tab.2, Fig. 12-13), which is a measure of compliance of observed and calculated water stage hydrographs variation, shows that 1D model is unable to reproduce flood dynamics in a vast part of floodplain (Fig.10). Two divers D3 and D4 are located in places in the river, where surrounding floodplain is very wide, and the calculated R coefficients for these sensors have small values: 0.43 and 0.69 respectively. The R coefficient for the two divers located in a narrow part of valley is higher than 0.9, and observed and calculated water stage hydrographs are congruent (Fig. 10-11).

Fig. 10. Observed and calculated water stage at D1-Pale diver (01.04.2006 – 30.04.2006) (H – water level [m asl.])

Fig. 11. Observed and calculated water stage at D2-Klimaszewnica diver (01.04.2006 – 30.04.2006) (H – water level [m asl.]

Fig. 12. Observed and calculated water stage at D3-Kosódka diver (01.04.2006 – 30.04.2006) (H – water level [m asl.])

The elaborated model was validated for the flood in 2007 (1 March – 30 April), which was characterized by a significantly higher water level than the flood in 2006. The maximum water level in Burzyn gauge was 101.76 m a.s.l. in 2006 and amounted to 102.28 m a.s.l. in 2007. In this case the model gives similar results to the ones of water stage observations performed at D1-D4 divers (Tab.2).

Fig. 13. Observed and calculated water stage at D4-Chyliny diver (01.04.2006 – 30.04.2006) (H – water level [m asl.])


Table 2. Variation of the mean square error (RMSE) and correlation coefficient (R) for water stage-divers (model calibration and validation)

#### **3.5 Scenarios determination**

160 Water Resources Management and Modeling

0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 time [days]

0 3 6 9 12 15 18 21 24 27 30 33 36 39 42 45 48 51 54 57 60 time [days]

0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 time [days]

Fig. 12. Observed and calculated water stage at D3-Kosódka diver (01.04.2006 – 30.04.2006)

model D3 - measurement

model D1 - measurement

Fig. 10. Observed and calculated water stage at D1-Pale diver (01.04.2006 – 30.04.2006) (H –

model D2 - measurement

Fig. 11. Observed and calculated water stage at D2-Klimaszewnica diver (01.04.2006 –

106.75 106.80 106.85 106.90 106.95 107.00 107.05 107.10

105.00 105.10 105.20 105.30 105.40 105.50 105.60 105.70

30.04.2006) (H – water level [m asl.]

103.00

103.10

103.20

103.30

H [m asl

(H – water level [m asl.])

103.40

103.50

103.60

H [m asl

H [m asl

water level [m asl.])

Various types of land use were determined for three scenarios and two of them present protective activities carried out nowadays in the National Park, which consist of preservation of open grasslands by mowing, elimination of biomass in selected areas, and mechanical scrub cutting. In the first analyzed scenario (scenario I), wet meadows, sedge plant communities and pastures are areas of extensive agricultural use. In this scenario, these areas are mowed for cattle bedding, used for hay production and as feeding areas. In the second scenario (scenario 2), there is willow scrub and birch tree cutting, in addition to the former scenario. The overgrowing process in many cases seem to be a natural succession so the third scenario (scenario 3) allows the natural succession of willow shrubs and birch forest upon non-forest ecosystems in the valley.

The numerical calculations for the three scenarios, which show changes in the land use, were made for the maximum flood event in the LBRB in 1979. In scenario I, roughness height coefficients for wet meadows, tall sedge plant communities, and pastures were given according to values in Tab.1. In the model scenario II, it was described by *ks* = 0.4 m for the floodplain area with willow scrubs. In the last simulation (scenario 3), roughness coefficients for the study site were established like willow scrubs (Tab.2) because of natural succession in this area.

#### **3.6 Results and discussion**

The results were analyzed in four river cross sections which represent a variability of the hydromorphological shape of the main channel and floodplain in the Lower Biebrza River Basin. The distances from BD17, where the river course ends, to the selected cross sections are 7.62 km, 29.62 km, 33.04 km and 39.96 km. Fig. 14 and Fig. 15 show the calculated hydrographs of water stages and water stage profile of the maximum flood event of 03.04.1979 for the various land use scenarios. The average depth in scenario I and scenario II was decreased by 0.25 m compared to the actual state (scenario 0). In scenario III, the average depth was increased by 0.20 m. Some differences in the water stage at the outlet (near BD17) are related to the downstream boundary condition which is formed by the rating curve, which is the same for all scenarios. Therefore in calculation, the influence of the vegetation structure on the rating curve was excluded which results in lack of influence of a land use type on a calculated water stage at BD17. Tab. 3 shows average water depth in the floodplain, top flood width, and flood duration for the various land use scenarios.

The results for scenario I and scenario II are similar, because of willows scrubs covering a small area in the valley (Fig. 5). The highest changes of flooding time were calculated for P4 localized in the northern part of the Lower Basin, where the width of the river valley is lower than in the middle and southern parts of the research area. Time of flood in scenario I and II, when the resistance of flow is lower, was decreased by about 40 %. When the resistance of flow is higher because of growing vegetation (scenario III), flood duration was increased by 17 % compared to the actual state.

An increasing valley width resulted in decreased flooding for scenarios I and II by 8 to 6 % along the river from P17 to BD14. In scenario III, flood duration was increased by 5, 7, and 9 % for P17, BD6, and BD14 in comparison to the scenario 0. Changes in the average depth do not exceed 10 %. The variation of flood extent is connected with the cross sections geometry. The highest decrease of width occurs at P4 for scenarios I and II and is about 40 %. In these scenarios it is about 10 % at P17 and BD6. The flood extent was increased by several percent in scenario III at the analyzed cross sections. In this scenario birch forest and willow shrubs occurred over large areas of wetland. The flood extent, water depth and flood duration increased which means that the area of plant communities dependent on rich surface water from the river like Phragmition and Magnocaricion will increase at the cost of more valuable plant communities (like sedge-moss meadows). Long lasting inundation also does not support the development of woody vegetation except the alluvial forests. Using the hydrodynamic model it was impossible to predict to what extent changes in flooding characteristic will affect the vegetation, changing in the next turn its structure and location. We can only state here that secondary succession probably does not form climax vegetation on the floodplain.

The numerical calculations for the three scenarios, which show changes in the land use, were made for the maximum flood event in the LBRB in 1979. In scenario I, roughness height coefficients for wet meadows, tall sedge plant communities, and pastures were given according to values in Tab.1. In the model scenario II, it was described by *ks* = 0.4 m for the floodplain area with willow scrubs. In the last simulation (scenario 3), roughness coefficients for the study site were established like willow scrubs (Tab.2) because

The results were analyzed in four river cross sections which represent a variability of the hydromorphological shape of the main channel and floodplain in the Lower Biebrza River Basin. The distances from BD17, where the river course ends, to the selected cross sections are 7.62 km, 29.62 km, 33.04 km and 39.96 km. Fig. 14 and Fig. 15 show the calculated hydrographs of water stages and water stage profile of the maximum flood event of 03.04.1979 for the various land use scenarios. The average depth in scenario I and scenario II was decreased by 0.25 m compared to the actual state (scenario 0). In scenario III, the average depth was increased by 0.20 m. Some differences in the water stage at the outlet (near BD17) are related to the downstream boundary condition which is formed by the rating curve, which is the same for all scenarios. Therefore in calculation, the influence of the vegetation structure on the rating curve was excluded which results in lack of influence of a land use type on a calculated water stage at BD17. Tab. 3 shows average water depth in the

floodplain, top flood width, and flood duration for the various land use scenarios.

The results for scenario I and scenario II are similar, because of willows scrubs covering a small area in the valley (Fig. 5). The highest changes of flooding time were calculated for P4 localized in the northern part of the Lower Basin, where the width of the river valley is lower than in the middle and southern parts of the research area. Time of flood in scenario I and II, when the resistance of flow is lower, was decreased by about 40 %. When the resistance of flow is higher because of growing vegetation (scenario III), flood duration was

An increasing valley width resulted in decreased flooding for scenarios I and II by 8 to 6 % along the river from P17 to BD14. In scenario III, flood duration was increased by 5, 7, and 9 % for P17, BD6, and BD14 in comparison to the scenario 0. Changes in the average depth do not exceed 10 %. The variation of flood extent is connected with the cross sections geometry. The highest decrease of width occurs at P4 for scenarios I and II and is about 40 %. In these scenarios it is about 10 % at P17 and BD6. The flood extent was increased by several percent in scenario III at the analyzed cross sections. In this scenario birch forest and willow shrubs occurred over large areas of wetland. The flood extent, water depth and flood duration increased which means that the area of plant communities dependent on rich surface water from the river like Phragmition and Magnocaricion will increase at the cost of more valuable plant communities (like sedge-moss meadows). Long lasting inundation also does not support the development of woody vegetation except the alluvial forests. Using the hydrodynamic model it was impossible to predict to what extent changes in flooding characteristic will affect the vegetation, changing in the next turn its structure and location. We can only state here that secondary succession probably does not form climax vegetation

of natural succession in this area.

increased by 17 % compared to the actual state.

on the floodplain.

**3.6 Results and discussion** 

Fig. 14. Calculated stage hydrograph at select cross sections for different land use scenarios

Fig. 15. Calculated water stage profile for different land use scenarios (03.04.1979)


Table 3. Variability of the average water depth on the floodplain, top flood width and flood duration for different land use scenarios

### **4. Conclusions**

The application of a flow law in the model, in the form of the Darcy-Weisbach formula, made it possible to consider the influence of vegetation on water flow conditions. It allowed for implementation of equations mass and momentum conservation processes between the main riverbed and the floodplains in 1D St Venant.

In order to explore the potential of the model use for the analysis of unsteady flood flows, simulations of water flow were performed for the Lower Biebrza Basin, which is a natural reference valley. The model was calibrated and validated for floods which took place in 2006-2007 based on the data achieved as a result of executed field research. During model calibration process there were generally two difficulties: the determination of inflow from a sub-catchment and delineation of those areas of a very flat and vast valley that take part in water flow.

The solution to the first problem was applied in the model in the form of uniformly distributed lateral inflow estimated as a result of numerical experiments, and, even if it

Scenario

Table 3. Variability of the average water depth on the floodplain, top flood width and flood

The application of a flow law in the model, in the form of the Darcy-Weisbach formula, made it possible to consider the influence of vegetation on water flow conditions. It allowed for implementation of equations mass and momentum conservation processes between the

In order to explore the potential of the model use for the analysis of unsteady flood flows, simulations of water flow were performed for the Lower Biebrza Basin, which is a natural reference valley. The model was calibrated and validated for floods which took place in 2006-2007 based on the data achieved as a result of executed field research. During model calibration process there were generally two difficulties: the determination of inflow from a sub-catchment and delineation of those areas of a very flat and vast valley that take part in

The solution to the first problem was applied in the model in the form of uniformly distributed lateral inflow estimated as a result of numerical experiments, and, even if it

Average depth [m]

0 0.58 2346 53 I 0.51 1432 32 II 0.51 1332 31 III 0.64 2400 62

0 0.24 3439 86 I 0.19 3010 79 II 0.19 3008 79 III 0.29 3744 90

0 0.32 5103 92 I 0.28 4685 86 II 0.28 4666 85 III 0.35 5367 100

0 0.49 6950 97 I 0.44 5655 91 II 0.44 5592 91 III 0.53 7182 104

Top width [m]

Flood duration [days]

Distance along the river from BD17 [km]

Number

Cross– section name

1 P4 39.96

2 P17 33.04

3 BD6 29.62

4 BD14 7.62

duration for different land use scenarios

main riverbed and the floodplains in 1D St Venant.

**4. Conclusions** 

water flow.

leads to a satisfactory match of flow hydrographs at the closing water gauge, it is a very rough solution. A proper solution should involve the elaboration of a hydrological model for the LBB with the consideration of river water-ground water interaction. The application of flow capacity calculations allowed for delineation of that part of the valley which takes active part in a flood wave propagation. In future this approach should be coupled with a field monitoring.

The measurement technique used for inundation dynamics monitoring, which involves installation of automatic sensors of D-Diver type in the river bed, is considered to be the most reasonable measurement solution in the field of achieving data on water levels in the moment of high water occurrence. That sort of data is indispensable to perform proper identification and validation of a hydrodynamic model. In order to properly investigate the variability of inundation in the Lower Biebrza Valley, that type of sensors should not only be installed at selected points located in the river, but also on the floodplains. The location of measurement points may be chosen on the basis of collected measurement data and GIS analyses against the indications of existing digital elevation model. Such an organized monitoring network should lead to a more complete recognition of spatial and temporal variability of inundation in that area.

The application of aerial photographs, vegetation maps and time consuming field monitoring made it possible to delineate the low and high vegetation communities on the floodplains and to determine the values of roughness of high and mean geometrical parameters for each plant community.

In the case of the Lower Biebrza the biggest differences of water levels occur in those places where the valley is very wide and the simplification of water flow to a one dimensional case results in impossibility to reproduce real dynamics of the floods. The water flow in this part of the valley can be performed by a 2D model.

The elaborated model, which was calibrated and validated for the area of the LBB, allows for recognition the dynamics of inundation process and the range of inundation in that region. It is particularly vital for the implementation of an effective policy of natural values protection in the area of the Biebrza National Park

The obtained results show that vegetation structure on the floodplain influences the flood wave propagation in the LBRB. An intensive land use on the floodplain decreased the resistance of flow, which decreased flood duration, flood extent, and average flood water depth. These flood characteristics are strongly connected with the cross section geometry. The valley in the northern part is narrow, which causes flood duration shortening; for an intensive land use on the floodplain the time was decreased by about 40 % compared to the current state, while in the other part of the Lower Basin - by less than 10 %. The influence of the natural succession on the flood wave is less significant than the influence of intensive agriculture. In this case, flood duration changes by a few percent.

The achieved results were obtained with the use of a one dimensional hydrodynamic model which does not cover water flow across the river. It means that the used model represents simplified flow conditions in the Lower Biebrza River Basin, which is a wide area and for which this aspect is very important. The essential influence also has assigned ineffective flow areas that contain water which is not actively conveyed. Ineffective flow areas are used to describe portions of a cross section in which water ponds, but velocity of that water is close to or equal to zero.

Moreover, the hydrodynamic model is an appropriate tool for assessment of various agricultural practices with regard to the effective management to protect the unique wetland site of the Biebrza National Park.

The developed model makes it possible to find a solution for unsteady flow problems in natural rivers with vegetated floodplains. It may be used as a tool to estimate new water surface level for renaturalized rivers, especially for flood conditions, as well as to ensure suitable conditions for habitat diversity in projects of environmental flood management. It is an appropriate tool for estimation of floodplain vegetation influence on flow conditions, which can be considered a 1D problem.

#### **5. References**


which this aspect is very important. The essential influence also has assigned ineffective flow areas that contain water which is not actively conveyed. Ineffective flow areas are used to describe portions of a cross section in which water ponds, but velocity of that water is

Moreover, the hydrodynamic model is an appropriate tool for assessment of various agricultural practices with regard to the effective management to protect the unique wetland

The developed model makes it possible to find a solution for unsteady flow problems in natural rivers with vegetated floodplains. It may be used as a tool to estimate new water surface level for renaturalized rivers, especially for flood conditions, as well as to ensure suitable conditions for habitat diversity in projects of environmental flood management. It is an appropriate tool for estimation of floodplain vegetation influence on flow conditions,

Bretschneider H., Schulz A. 1985: Anwendung von Fließformeln bei naturnahem

Chormanski J. 2003: Methodology for the flood extent determination in the Lower Biebrza

DVWK – MERKBLÄTTER 220, 1991: Hydraulische Berechnung von Fließgewässern, *DK* 

Cunge, J., F. Holly Jr., and A. Verwey, Practical Aspects of Computational Hydraulics,

Helmio , T. 2002, Unsteady 1D flow model of compound channel with vegetated

Hooijer A. 1996: Floodplain hydrology: An ecologically oriented study of the Shannon

Horritt, M.S. and Bates, P.D. 2002. Evaluation of 1-D and 2-D numerical models for

Indlekofer, H., 1981, Überlagerung von Rauhigkeitseinflüssen beim Abfluß in offenen

Kaiser, W., 1984, Fliewiderstandsverhalten in Gerinnen mit durchströmten

Kubrak J., Nachlik E. (red.), 2003: Hydrauliczne podstawy obliczania przepustowości koryt

Laks I., Kałuża T., 2007, Analiza przejścia fali wezbraniowej z 1997roku na odcinku od

Gerinnen. *Mitt. Institut für Wasserbau und Wasserwirtschaft,* RWTH Aachen, Heft 37,

Ufergehölzzonen. *Thesis presented for the degree of Doctor in Applied Sciences* TH

Jeziorska do Obornik w modelu obliczeniowym, uwzględniającym aktywną część

predicting river flood inundation. Journal of Hydrology, 268, 87-99.

*551.51/54 Fließgewässer, DK 532.543 Hydraulik, DVWK - Merkblätter 220/1991*,

close to or equal to zero.

**5. References** 

s. 105-145.

Darmstadt.

site of the Biebrza National Park.

which can be considered a 1D problem.

Gewässerausbau. *DVWK - Schriften, Heft 72.*

*Pitman Publishing Ltd*., London, UK, 1980.

rzecznych. *Wydawnictwo SGGW,* 317 stron.

przekroju. Nauka Przyr. Technol. 1, 2, #27.

floodplains*, J. Hydrol*., 269(1– 2), 89–99.

Basin. *Ph.D. thesis, Warsaw Agricultural University,* 186.

Callows, Ireland. Thesis, Virje Universiteit Amsterdam

Kommissionsvertrieb Verlag Paul Parey, Hamburg und Berlin.


Zalewski M., Janauer G.A., Jolankai G. 1997: A new paradigm for the sustainable use of aquatic resources. *UNESCO IHP Technical Document in Hydrology* 7. IHP - V Projects 2.3/2.4, UNESCO Paris, Ecohydrology
