Effects of Hydropower Dam Operation on Riverbank Stability

Abstract

The increasing number of extreme climate events has impacted the operation of reservoirs, resulting in drastic changes in flow releases from reservoirs. Consequently, downstream riverbanks have experienced more rapid and frequent changes of the river water surface elevation (WSE). These changes in the WSE affect pore water pressures in riverbanks, directly influencing slope stability. This study presents an analysis of seepage and slope stability for riverbanks under the influence of steady-state, drawdown, and peaking operations of the Roanoke Rapids Hydropower dam on the lower Roanoke River, North Carolina, USA. Although the riverbanks were found to be stable under all the discharge conditions considered, which indicates that normal operations of the reservoir have no adverse effects on riverbank stability, the factor of safety decreases as the WSE decreases. When the role of fluvial erosion is considered, riverbank stability is found to reduce. Drawdown and fluctuation also decrease the safety factor, though the rate of the decrease depends more on the hydraulic conductivity of the soils rather than the discharge pattern.

1. Introduction

Slope stability is generally expressed in terms of its factor of safety (FS), which is defined as the ratio of the available shear strength of the soil to the equilibrium shear stress that is required to maintain a just-stable slope [1]. Several significant factors contribute to stability: slope geometry, soil weight, soil shear strength, location of the phreatic surface, and degree of saturation. When considering a riverbank, these factors are directly affected by the water surface elevation (WSE).

The WSE acts as an external force on the slope surface, confining the riverbanks and affecting the location of the phreatic surface, which influences pore water pressure and the weight of the soil. The WSE also affects the shear stress applied to the riverbank by flow, thus controlling fluvial erosion. The decrease in riverbank stability due to fluvial erosion and the associated changes in bank geometry are well documented in the literature [2,3,4], and cumulative fluvial erosion at the toe of a slope is known to trigger mass failures [5]. Thus, a critical evaluation of the influence of reservoir releases on the stability of the downstream riverbank can prove a helpful tool for the operation of hydropower dams to reduce riverbank failures.

In recent decades, the global climate has changed dramatically, resulting in an increased frequency of extreme events, such as hurricanes, heavy rains, and snows. Storms have become more intense and prolonged in duration [6], and the proportion of category 4 and 5 hurricanes has doubled since 1970 [7]. Thus, there is a higher likelihood that reservoir releases may need to be adjusted in response to these events. More instantaneous dam operations may be required to maintain the reservoir’s safety and the downstream floodplain. However, immediate releases create flow conditions that are not commonly found during the normal functions of the dams. Instead, extreme flows or drastic changes in the downstream WSE may lead to overbank flow or transient drawdown conditions. As a result, riverbank stability under prolonged high flow conditions or drawdown dam releases needs to be considered in slope stability analysis. In addition, hydropower dams operate with consideration of power generation and factors such as fish spawning, drought/flood controls, and recreational use that are typically considered in other reservoirs. Each operational mode focused on dominant factors usually has different discharge patterns that also affect the WSE.

The designated sites for the study presented in this paper were exposed to these types of releases due to the operation of an upstream hydropower dam, the Roanoke Rapids Dam, which directly affects the downstream WSE. The dam and downstream area experienced extreme conditions due to prolonged heat and drought, followed by high flow releases and more frequent peaking for the increased power demand. Thus, it is essential to understand how the dam release patterns impact the downstream flow condition. Consequently, it is critical to determine how to control reservoir releases to minimize the potential negative impacts on the riverbank.

This study aimed to compare the stability of riverbanks under normal and extreme reservoir operations of an upstream hydropower dam that can cause overbank conditions, frequent changes of the WSE, or transient drawdown. The results can be directly interpreted to determine the optimum dam operational modes to minimize bank retreat. The unique contributions of this study include: (i) riverbank stability analysis, taking the actual reservoir operations and downstream conditions into account; (ii) the consideration of unsaturated soil properties for the transient seepage and slope stability analysis; and (iii) the prediction of riverbank stability, taking fluvial erosion into account using actual erosion rates monitored by the U.S. Geological Survey (USGS).

2. Background

2.1. Site Description and Materials

Five study sites, shown in Figure 1, are located on the lower Roanoke River near Scotland Neck, NC, USA, about 72 km downstream from the Roanoke Rapids Dam. Since the dam’s construction was completed in 1955, the river flow has been regulated, resulting in a decrease in extreme high flow events, as shown in Figure 2. The WSE at the study sites is primarily controlled by flow releases from the dam because of the lack of major tributaries (see Figure 1b), the limited effect of precipitation due to the narrow watershed, and the heavy vegetation and swamps that surround the study sites.

The study sites were selected based on evidence of past bank retreats. Sites 2, 4, and 5 are located on relatively straight reaches, while Site 1 and 3 are situated on the outer and inner banks of meander bends, respectively. These sites were close enough to the dam to be directly affected by its peaking releases. Over the course of three years, eight field trips, for total of 55 days, were made to install and monitor river and ground water elevations, to perform in situ tests, and to collect hydraulic and geotechnical data representing field conditions under reservoir release rates ranging from 57 m³/s to 566 m³/s. Visual observation during the field trips confirmed that the riverbanks are relatively steep, indicating that fluvial erosion likely occurs at these sites. In addition, the presence of trees rooted in the submerged slopes within the river channel indicated that deep-seated mass failure events had already occurred. Small scale local failures near the WSE were also noticeable, but both failure types seemed to be randomly distributed throughout the study reach. Other researchers have noted this. Hupp et al. [8] observed and documented bank erosion and mass failure in their study of the same reach. At the same time, the USGS also monitored bank erosion at several locations along the study reach between 2005 and 2009 [9].

Hydraulic properties, such as field flow velocity, discharge rate, and bathymetry, were obtained by acoustic Doppler current profiler (ADCP) and echosounder. The complete bank geometry required for the numerical modeling was obtained by combining the bathymetric data with ground-based light detection and ranging (LiDAR) data at each site.

2.2. Soil Properties

After selecting the study reach and five study sites, disturbed and undisturbed soil samples were collected during the field trips. The soil samples were then analyzed for basic physical properties, hydraulic conductivity, and shear strength for saturated and unsaturated conditions.

Surface soil in the area is primarily Quaternary alluvium deposits up to 7.6 m depth, with Upper Cretaceous sedimentary materials underlying the alluvium soil layer [10]. A loose layer of silty sand exists at the top of the riverbank, and a firm and thick cohesive soil layer underlies it. The underlying layer is classified as CL, ML, or MH by the Unified Soil Classification System (USCS). Even though the thicknesses of the surface and underlying layers in the field were different by the sites, and the soils are classified as different types, the physical properties seemed to be relatively consistent regardless of the sites. For example, a total of 30 soil samples were tested for the liquid limit and plastic limit (18, 4, 3, 2, and 3 samples from Sites 1 to 5, respectively) and were predominantly CL of MH, as shown in Figure 3.

Shelby tube and block samples were mainly used where undisturbed soil samples were necessary. The shear strength parameters for the saturated and unsaturated soils were determined by conventional and suction-controlled multi-stage direct shear tests (MDST). A conventional direct shear test device was modified to allow suction control and multi-stage loading. Even though limitations of the direct shear test and multi-stage loading exist, the test method and results have been proven to be reliable for both saturated and unsaturated soils [11,12,13,14,15]. Particularly for unsaturated soils of lower permeability that generally require a long equilibrium period under controlled matric suction conditions, the MDST can be a practical choice because of having a shorter test time and, thus, has been used successfully [16,17,18]. A total of 15 MDSTs were successfully performed, with five tests conducted under the suction-controlled condition. All four soil types at Site 1 were tested for both conventional and suction-controlled multi-stage loading conditions, whereas the soils from other sites were tested for the conventional multi-stage loading condition only. For the suction-controlled MDST, each sample was tested with a fixed net normal stress of 43 kPa, except the MH soil, with varying the matric suction from 25 up to 290 kPa. With this test setup, the angle of shearing resistance for matric suction (${\varphi }^b$) was determined. The experimental results showed the expected non-linear relationship between the soil suction and shear strength, which was in agreement with others’ findings [16,19,20]. The non-linear behavior is often simplified with a bi-linear model with ${\varphi }^b$ changing around the air entry value (AEV). In this study, the slope of the model, after the AEV, was selected for the numerical models for a more conservative analysis. More information regarding the MDST and the shear strength parameters are presented in Nam et al. [18].

In addition, the borehole shear test (BST) was also conducted in the field. The BST determines the shear strength parameters in the field and is known to show good agreement with the values measured in the laboratory [21,22,23,24]. A total of 26 BSTs were conducted in the field for Sites 1–5. Except for the upper CL soil, the results from both tests seem to be consistent and the average values of ${\varphi }^{\prime }$ and $c^{\prime }$ were used for the modeling. The results of the MSDS and BST are presented in

Table 1. Shear strength parameters determined by conventional and suction-controlled MSDS and BST.

Site	USCS Soil Type	$$\begin{gathered} {\mathit{\varphi }}^{\prime }{}_{\mathit{M}\mathit{S}\mathit{D}\mathit{S}} \\ (°) \end{gathered}$$	$$\begin{gathered} {\mathit{\varphi }}^{\prime }{}_{\mathit{B}\mathit{S}\mathit{T}} \\ (°) \end{gathered}$$	$$\begin{gathered} {\mathit{c}}^{\prime }{}_{\mathit{M}\mathit{S}\mathit{D}\mathit{S}} \\ \left(\mathbf{kPa}\right) \end{gathered}$$	$$\begin{gathered} {\mathit{c}}^{\prime }{}_{\mathit{B}\mathit{S}\mathit{T}} \\ \left(\mathbf{kPa}\right) \end{gathered}$$	$$\begin{gathered} {\mathit{\varphi }}^{\mathit{b}} \\ (°) \end{gathered}$$
S1	SM	35.0	31.0	4.3	6.4	13.3
	CL	32.4	24.2	5	21.5	10.2
	MH	32.0	32.2	15.8	21.0	13.5
	CL	28.4	27.7	22.5	15.0	9.0
S2	CL	N/A	29.5	N/A	7.5	N/A
	CL	34.1	N/A	12.8	N/A	N/A
S3	ML	N/A	28.8	N/A	6.3	N/A
	CL	35.5	31.8	6.8	10.5	N/A
S4	ML	N/A	28.2	N/A	9.5	N/A
	CL	34.5	N/A	8.9	N/A	N/A
S5	ML	33.7	32.1	11.8	4.5	N/A
	CL	34.6	28.7	7.5	12.5	N/A

.

The hydraulic conductivity was measured by laboratory and field tests, and was determined by comparing the transient seepage analysis and the groundwater table observations [25]. A total of 17 constant head permeability tests, with a constant air pressure, and oedometer tests were conducted in the lab. In addition, 19 auger hole tests and Guelph permeameter tests were conducted in the field. The soil-water characteristic curves (SWCC) were obtained from Tempe cell tests, pressure plate tests, a dewpoint potentiometer, filter paper tests, a vapor equilibrium test, and an osmotic technique test to characterize the relationship between suction and water content using disturbed and undisturbed samples from Site 1, and the air entry values were determined from the SWCC. The detailed results are reported in Nam et al. [26]. Due to the restricted site access and limited availability of samples and testing apparatuses, some of the input parameters for the modeling were assumed.

All values used for the modeling are listed in

Table 2. Soil properties for seepage and riverbank stability modeling.

Site	USCS Soil Type	$$\begin{gathered} {\mathit{\gamma }}_{\mathit{t}} \\ (\mathbf{kN}/{\mathbf{m}}^3) \end{gathered}$$	$$\begin{gathered} {\mathit{\varphi }}^{\prime } \\ (°) \end{gathered}$$	c’ (kPa)	AEV (kPa)	$$\begin{gathered} {\mathit{\varphi }}^{\mathit{b}} \\ (°) \end{gathered}$$	VWC (%)	kLab (m/s)	kAuger (m/s)
S1	SM	16.4	33	5.4	10	13.3	46.1	5.09 × 10−7	1.84 × 10−4
	CL	17.7	28.3	13.3	120	10.2	50.3	7.32 × 10−10	2.64 × 10−5
	MH	18	32.1	18.4	160	13.5	48.3	4.99 × 10−9	1.35 × 10−5
	CL	18.5	28.1	18.8	200	9	47.8	1.02 × 10−9	2.58 × 10−5
S2	CL	18.2	29.5	7.5	120 *	10.2 *	48.5	1.65 × 10−8 *	2.64 × 10−5 *
	CL	19.1	34.1	12.8	180 *	11.3 *	47.7	3.13 × 10−9 *	1.97 × 10−5 *
S3	ML	17.6	28.8	6.3	120 *	10.2 *	52.5	1.65 × 10−8	2.13 × 10−5
	CL	18.7	33.7	8.7	180 *	11.3 *	49.7	8.30 × 10−9	1.97 × 10−5 *
S4	ML	17.9	28.2	9.5	120 *	10.2 *	49.5	1.72 × 10−8	2.13 × 10−5 *
	CL	18.3	34.5	8.9	180 *	11.3 *	48.1	3.13 × 10−9 *	1.97 × 10−5 *
S5	ML	17.4	32.9	8.2	120 *	10.2 *	48.5	5.90 × 10−9	3.15 × 10−5
	CL	17.9	31.6	10	180 *	11.3 *	47.3	5.13 × 10−9	1.97 × 10−5 *
S1–S5	N/A	22 *	37 *	200 *	0 *	0 *	26 *	1.26 × 10−10 *	1.26 × 10−7 *

AEV: air entry value; VWC: volumetric water content; * Assumed values; N/A: soil type by the USCS was not available as the layer was assumed to be hard and impermeable.

, and the bank geometries and WSEs used for the seepage analysis and slope stability analysis are provided in Figure 4. Different flow patterns were considered for the transient analysis.

2.3. Reservoir Release Patterns

As shown in

Table 3. Dam operational modes and range of flow rates [modified from Dominion [27]].

Operation Mode	Drought Control		Fish Spawning		Normal		Flood Control
	From	To	From	To	From	To	From	To
Discharge rate (m3/s)	57 1	N/A	99	388	57	566	566	991 2

¹ 42 m³/s between September and November. ² Greater of 100% of the inflow to the dam or 991 m³/s.

, the Roanoke Rapids Dam has four operational modes: normal, fish spawning, drought control, and flood control modes. During regular operation, which occurs throughout the year except during the fish spawning season, the reservoir release rate can vary from 57 to 566 m³/s and is limited only by the hydraulic capacity of the power station [27]. Thus, the discharge is expected to be dependent on the seasonal inflow rate to the dam. However, more periodic and frequent peaking is anticipated to maximize the efficiency of power generation. During spawning season, between 1 March and 15 June, the downstream WSE is maintained at a relatively uniform flow with only minor variations allowed to minimize WSE fluctuations. If drought or flood conditions are declared, the flow release is governed by the John H. Kerr Dam discharge, which is located further upstream and is controlled by the U.S. Army Corps of Engineers. During drought operations, the flow tends to remain at 57 m³/s, and peaking events are reduced, resulting in a low steady downstream WSE. Flood control does not limit the minimum flow, but can release up to 991 m³/s or more depending on the inflow to the reservoir.

Figure 5 depicts an example of the actual release rate at the dam and the resulting downstream WSE near Scotland Neck, NC, as monitored at the USGS stream gage between October 2005 and November 2009. As the figure shows, although the dam has four distinct operational modes for regulating the release rate, the actual release rate and downstream WSE are not easily categorized in terms of these modes. Thus, several representative release rates and flow events were selected for this study based on historical data.

3. Seepage and Slope Stability Analyses

3.1. Transient Seepage Analysis

Transient seepage analysis was utilized in this study to analyze the impact of fluctuating WSE due, mainly, to the influence of dam operations on riverbank stability. When a riverbank is exposed to rapid changes in WSE, the pore water pressure inside the riverbank becomes an essential factor determining slope stability. Excess pore water pressure develops and dissipates over time, which directly affects the riverbank stability. Depending on the soil properties, bank geometry, and changes in the WSE, the influence of excess pore water pressure may become critical for slope stability. Thus, transient seepage analysis was applied to estimate how the excess pore water pressure changes with time.

Transient seepage analysis is typically performed using numerical methods based on Darcy’s Law. Darcy’s Law was originally derived for water flow through saturated soil, but it can be extended to include flows in unsaturated soil. The governing equation for two-dimensional transient flow through unsaturated soils is as follows:

$$\frac{\partial }{\partial x}\left(k_x\frac{\partial h}{\partial x}\right) + \frac{\partial }{\partial y}\left(k_y\frac{\partial h}{\partial y}\right) + q=\frac{\partial {\theta }_w}{\partial t}$$

(1)

where h = total hydraulic head, k_x and k_y = hydraulic conductivity in the x and y directions, respectively, q = external boundary flux, and θ_w = volumetric water content.

The significant difference between water flow in saturated and unsaturated soils is that the hydraulic conductivity varies with saturation in unsaturated soils. In contrast, for saturated flow, hydraulic conductivity remains constant.

3.2. Slope Stability Analysis

The limit equilibrium method (LEM) is a classical approach used to evaluate the stability of a slope and was developed before the application of computers. More advanced techniques, using numerical methods and more powerful computational resources, have since been developed for slope stability analysis. However, LEM remains the most widely used slope stability analysis method due to its simplicity and ease of use.

LEM presents the stability of a slope in terms of the factor of safety (FS). Duncan and Wright [1] defined FS as the ratio of the available shear strength of soil to the shear stress required to maintain a just-stable slope, expressed as follows:

$$FS=\frac{s}{\tau }$$

(2)

where s = available shear strength and τ= equilibrium shear stress.

Leshchinsky [28] described the calculation of the FS of a slope in terms of a two-step process: (1) assuming a potential slip surface, and (2) assembling and solving the limit equilibrium equations for the soil mass defined by the surface. Both processes have benefitted from the application of numerical methods. A large number of potential slip surfaces, including non-circular types, can be examined, and various methods applying different assumptions when defining equilibrium can be utilized. As a result, LEM provides more accurate results and is used more frequently than advanced numerical methods.

Once the basic concepts of unsaturated soil mechanics were proposed, slope stability analyses could be performed using the shear strength of unsaturated soils. This is logical, as most engineered and natural slopes are not fully saturated. Additional factors, such as precipitation, evaporation, or changes in boundary conditions, alter the degree of saturation in soils. The unsaturated shear strength must be considered when performing a slope stability analysis for more realistic modeling. Thus, the concept of matric suction, which is not considered in classical slope stability analysis, is now being incorporated into slope stability analyses. Matric suction can be related to the degree of saturation or volumetric water content using the SWCC. The modified Mohr-Coulomb failure criterion can describe the shear strength of an unsaturated soil as follows [29]:

$$\tau =c^{\prime }+\left({\sigma }_n-u_a\right)\tan {\varphi }^{\prime }+\left(u_a-u_w\right)\tan {\varphi }^b$$

(3)

where $\tau$ = the shear strength of unsaturated soil, $c\prime$ = the effective cohesion, ${\varphi }^{\prime }$ = the effective friction angle, ${\sigma }_n$ = normal stress, $u_a$ = pore air pressure, $u_w$ = pore water pressure, and ${\varphi }^b$ = the angle of shearing resistance for matric suction.

The calculations for transient seepage analysis and slope stability analysis using LEM have improved with advances in computational resources. Recent studies showed how these analyses can successfully be applied to solve various problems (e.g., [30,31,32,33,34,35,36]).

3.3. Coupling with Fluvial Erosion

One of the unique characteristics of riverbank stability is that fluvial erosion may also affect the bank stability, altering the bank geometry by creating steeper slopes. This often occurs through erosion of the bank toe, which may trigger mass failure. However, fluvial erosion is a complex phenomenon that involves both geotechnical and hydraulic characteristics and is difficult to estimate.

A common approach to estimating the erosion rate of fine-grained cohesive soils is to apply the linear excess shear stress equation as follows [37,38]:

$$\epsilon =k_d\left({\tau }_e-{\tau }_c\right)$$

(4)

where ε = the erosion rate, k_d = the erodibility coefficient, τ_e = the effective shear stress applied by the flow, and τ_c = the critical shear stress of the soil.

Although the threshold concept in Equation (4) may not represent the actual physical phenomenon, the equation has been used to predict the erosion rate of cohesive soils, often in conjunction with the results of submerged jet tests [39,40]. However, it is not straightforward to estimate all the parameters in the excess shear stress equation. In addition, although the effective shear stress in the equation may be simplified, for example, by assuming that the shear stress increases linearly with depth, the actual effective shear stress by the flow is much more complex in three-dimensional flows and bank geometry. The critical shear stress and erodibility coefficient are also challenging to estimate due to the spatial and temporal variability of natural soils and the complex interparticle forces in cohesive soils [41,42]. Thus, an experimentally evaluated erodibility of a point on a riverbank would not be appropriate to predict the erodibility and erosion profile of the riverbank in the future.

Schenk et al. [9] installed and monitored a series of pins to estimate the erosion rate in the lower Roanoke River between 2005 and 2009. The study sites presented here coincided with the locations of the erosion pins, with the exception of Site 4. At each site, 6 to 8 erosion pins were installed as normal onto the local bank slope, in positions ranging from near the WSE for drought conditions (57 m³/s) to near the bank crest (566 m³/s). Each erosion pin is approximately 1 m in length and 1 cm in diameter. If erosion occurs, the exposed portion is measured, and then the pins are pushed back to the slope surface. Pins buried due to deposition are located using a metal detector and their depth beneath the new surface is measured. Details of the techniques involved in measuring erosion pin data are provided in Schenk et al. [9]. The annual average erosion rates were calculated using the erosion pin data in 2006, 2007, 2008, and 2009. Although the actual erosion rates and profiles were different each year, these were averaged to roughly estimate erosion at the five study sites for the next ten years.

Fluvial erosion can be challenging for numerical modeling of slope stability. However, in practice, it is possible to neglect erosion for short times periods. In terms of hours or days, erosion is too small to induce changes in the pore water pressure distribution and bank geometry that are sufficient to affect slope stability calculation. Thus, riverbank stability with fluvial erosion was considered as a separate modeling case for different bank geometries, assuming that the monitored annual mean erosion rate is maintained for the next ten years and creates noticeable changes in the bank geometry.

4. Numerical Modeling

Seven steady-state discharges and three transient operations were considered as representative flow scenarios for the modeling. Riverbank stability with changing water surface elevation (WSE) was analyzed using the SLIDE software package developed by Rocscience [43]. Both transient seepage analysis using the finite element method and slope stability analysis using the limit equilibrium method can be used. The factor of safety is calculated by the Morgenstern-Price methods with half sign interslice force function [44].

For numerical modeling, boundary conditions need to be carefully defined in addition to the soil properties. The far-field boundary was assumed to be the constant head. The initial location of the groundwater table at the boundary was supposed to be located 1.3 m below the surface, which is the annual average location of the groundwater table (GWT), as shown in Figure 6. Solid dots present the depths of GWT below the surface between October 1994 and June 2010 at Roxobel, NC, as monitored by the USGS. Each box-and-whisker plot also shows the upper extreme, upper quartile, median, lower quartile, and lower extreme values for that month. The USGS station (USGS 361420077111407 BE-080) is about 12 km away from the lower Roanoke River and 14 km away from Site 1. It was assumed to represent the initial location of GWT for the modeling adequately.

The boundary condition of the riverbank surface was defined in terms of the given total head with a time function simulating the fluctuation of the WSE. The local WSE at each study site was determined by correlating the USGS gaging station data with the estimated WSE using steady-state HEC-RAS simulation and the WSE measured by differential pressure sensors at Site 1, 2, 4, and 5. The real-time surface water data at the USGS gaging stations were available from the USGS website [45]. The dam discharge was measured at the USGS gaging station near Roanoke Rapids, NC (USGS 02080500), and the comparable WSE data for areas near Scotland Neck, NC (USGS 02081000) and Oak City, NC (USGS 02081028) were also available.

The WSE and its fluctuations during the discharge events were monitored and reviewed. It was concluded that, while the locations of the WSE under the steady state conditions were different by sites, the changes of the WSE by discharge events were very similar. This seemed to be because: (i) the sites were relatively close to each other, (ii) the width of the river did not change significantly within the study reach, and (iii) the WSE at the sites were not affected by the local floodplain, but mainly controlled by the dam discharge.

4.1. Steady-State Condition

The steady-state condition assumes that the downstream WSE remains stable even after peaking occurs at the dam. Thus, steady-state conditions represent seasonal variations in the riverbank stability rather than a single event. Based on the regulations for each operational mode, summarized in Table 3, seven different flow rates representing the full range of fish spawning, flood control, draught control, and regular operational modes, along with an overbank flow representing an extreme flood control scenario, were considered for the steady-state slope stability analysis and are listed in

Table 4. Representative flow rates for modeling and expected corresponding water surface elevation at Site 1.

Discharge Rate		Relative WSE at Site 1 (m) (Bank Top = 0.00 m)	Representing Operational Mode
(m3/s)	(ft3/s)
57	2000	−5.80	Normal, drought
142	5000	−4.35	Normal, spawning
227	8000	−3.16	Normal, spawning
311	11,000	−2.29	Normal
425	15,000	−1.36	Normal
566	20,000	−0.30	Normal, flood
991	35,000	+1.10	Flood

.

4.2. Transient Condition—Peaking

Peaking discharge is a typical release pattern for hydropower dams that generate electricity, which is used to satisfy periods of high-power demand. These peaks in demand typically occur during the hottest and coldest days of the year and daily in the morning and late afternoon (Bob Graham, personal communication, 7 July 2008). During the monitored period between 2005 and 2009, the prolonged peaking event shown in Figure 7a was selected as an example of a regular peaking discharge pattern. This event occurred between 13 and 28 February 2007. Before the peaking period, the discharge rate was maintained at around 70 m³/s for seven days (6–12 February 2007). Then, repetitive peaks ranging from 57 m³/s to 566 m³/s were observed over the next 16 days. The discharge rate generally reached 566 m³/s and dropped instantly, resulting in a typical peaking cycle of about 6 h. However, the discharge often remained at the highest rate for several hours, extending the cycle to up to 12 h. Daily, the monitored peaking events typically occurred twice a day, once starting in the morning around 6 a.m. and again in the early evening around 6 p.m.

4.3. Transient Condition—Drawdown

In general, the external water surface level drawdown in a slope may lead to one of the most critical conditions affecting slope stability. Especially in a slope with cohesive soils, poor drainage causes slow changes in the pore water pressures in the slope. When the WSE drops, the external force that works as a stabilizing factor is lost, while excess pore pressures in the riverbank are slow to dissipate, leaving the slope in a less stable condition. This is called rapid drawdown [1]. Thus, estimating changes in the pore water pressures during transient flow releases is essential for checking the critical conditions of riverbank stability. The drawdown rate selected from the historical records between 19 and 27 January 2009, presented in Figure 7b, represents the fastest drawdown rate between 2005 and 2009.

4.4. Transient Condition—Step-Down

The step-down release is a type of drawdown pattern that includes buffer periods to minimize the environmental impacts induced by drastic flow rate changes. In general, after a prolonged high flow, drawdown is performed with several buffer stages, as shown in Figure 7c, instead of reducing the discharge to base flow instantly. The step-down rates and settings are predetermined and consider the anticipated inflow and outflow of the dam and the ecological and environmental impacts. The step-down scenario presented in this study occurred between 20 June and 9 July 2009. Before the step-down event, the release rate was around 580 m³/s, and the downstream WSE had maintained a full bank condition for 12 consecutive days, providing the initial steady-state condition. The release rate was then dropped to 85 m³/s via four buffer stages. At each stage, the release rate was maintained for about four days on average as a buffer period. The buffer period was longer at the higher rate (4.7 days at 430 m³/s) and shorter at the lower rate (2.8 days at 147 m³/s).

4.5. Fluvial Erosion

Fluvial erosion was estimated by interpreting the USGS erosion pin data. The average annual erosion rate was calculated from the four years between 2005 to 2009 gathered by [9]. As presented in

Table 5. Erosion rate calculated from USGS erosion pin data.

Site	Bank Location	Maximum Erosion		Minimum Erosion
		Rate (cm/year)	Location	Rate (cm/year)	Location
Site 1	Outside	−23.0	Pin 1	−3.0	Pin 6
Site 2	Straight	−5.4	Pin 2	+3.5	Pin 3
Site 3	Inside	−5.9	Pin 2	−0.2	Pin 4
Site 5	Straight	−8.1	Pin 1	−0.05	Pin 6

No erosion pins near Site 4; + deposition, − erosion.

, the maximum erosion was measured at the lowest part of the riverbanks. In contrast, the minimum erosion, or deposition in a few locations, occurred in the middle or upper part of the banks. Using the current bank geometry and the average annual erosion rate, the eroded bank geometry after ten years was predicted, as shown in Figure 8. The slope stability of the riverbank with the new geometry was analyzed for two steady-state, one peaking, and one step-down condition, for a total of four scenarios.

5. Results and Discussion

Transient seepage analysis was applied to calculate the pore water pressures in the riverbanks. The stability of the riverbanks was then analyzed using the pore water pressures for different flow scenarios that represent the field conditions under the various dam operations. As the stability depends directly on the hydraulic conductivity of the soils, the riverbanks were also reviewed for critical cases that used the lowest values determined by either constant head permeability tests or consolidation tests.

5.1. Steady-State

Riverbank stability at all five sites was analyzed for the steady-state flow conditions listed in Table 4, representing typical operational modes in different seasons, including one overbank flow condition. Examples of the slip surface at Site 1 and Site 4 are presented in Figure 9, indicating different depths of failures. Figure 10 illustrates the relationship between discharge rates and riverbank stability under steady-state flow conditions. All slopes were stable and became more stable as the water surface elevation (WSE) rose. The WSE was confirmed to be the dominant factor governing riverbank stability under steady-state conditions.

As shown in Figure 10, the increases in factor of safety as the discharge rate increased are much larger in Sites 1, 3, and 5 than in Sites 2 and 4 This result indicates that the WSE and bank geometry could be factors more critical than the soil properties. The factors of safety, estimated by both saturated and unsaturated shear strengths, are compared in Figure 11, indicating that the actual riverbank stability in the field would be higher than the numerical model due to the unsaturated soil conditions. The influence of unsaturated soil shear strength appeared larger for a steeper slope (Site 4) than for a gentle slope (Site 1).

5.2. Peaking

The riverbank at Site 1 was stable during the peaking events, as shown in Figure 12a. The safety factor values varied, responding directly to changes in the WSE, and this variation was more significant when lower hydraulic conductivity was considered. The FS was lower when the hydraulic conductivity of the soil was more significant, indicating more water seeped into the soils as the WSE rose and the area where matric suction was eliminated became larger, thus decreasing the shear strength of the soil. However, if peaking started when the initial steady-state condition was assumed to be higher, as presented in Figure 12b, the FS dropped more rapidly in the model with lower hydraulic conductivity (k_LAB). This can be attributed to the fact that pore water pressure dissipation in poorly drained soils takes longer. As the rate and magnitude of the WSE changes during peaking events are lower than those in a drawdown event, it was reasonable to assume that the peaking was less critical than the drawdown in terms of slope stability. This transient peaking effectively resulted in a mild case of rapid drawdown, which will be discussed next, but it was not expected to create a rapid rise of downstream WSE.

5.3. Drawdown

For the drawdown event studied here, which occurred in January 2009, the riverbanks at all five study sites experienced a reduction in their factor of safety, as shown in Figure 13. However, the values remained above 1.0 for all areas, indicating stable conditions. This finding seemed to be due to the soil’s relatively high hydraulic conductivity values for the model, which allowed the pore pressure to dissipate. The lower hydraulic conductivities suggested by laboratory tests indicated a much higher sensitivity to WSE changes, but the FS could not be quantified due to several numerical issues. These will be discussed later in more detail.

5.4. Step-Down

Similar to the results of the drawdown analysis, the step-down scenario (Figure 14) resulted in a decreasing factor of safety as the WSE decreased. However, while the factor of safety did reduce during step-down releases, the banks were found to remain stable. The amount of time at each discharge rate seemed sufficient to allow the pore water pressure to dissipate. Larger values of the hydraulic conductivity measured by in situ auger hole tests appear to be responsible for the nearly instantaneous changes in factor of safety. Excess pore pressure dissipates faster in easily drainable soils, and thus the factor of safety was directly related to the magnitude of confining pressure by the river. Additionally, the drawdown rate of the WSE during the step-down scenario was lower than that of the drawdown case, resulting in a more stable riverbank. Therefore, the riverbank appeared to be relatively insensitive to the step-down rate when considering large-scale instability.

5.5. Bank Stability with Fluvial Erosion

Erosion rates were determined for all study sites except Site 4, where no erosion pin data was available. The highest annual average erosion rate was 230 mm/year, as monitored at the lowest pin at Site 1. However, this seemed to be mainly due to mass failures. Erosion of 750 mm out of the 920 mm was measured at the pin between 2006 and 2008 by Schenk et al. [9]. The highest erosion rates at the other sites were 64 mm/year, which were mostly observed at the lowest pins. Some deposits were also measured at upper pins in Site 2. Using the mean erosion rate and erosion profile estimated by the USGS erosion pin data shown in Figure 8 as an example, the bank geometry after ten years was calculated. The stability of the predicted riverbanks was analyzed for 57 m³/s steady-state flow. This flow rate was selected because it produced the lowest factor of safety in the steady-state modeling. As the results in

Table 6. Erosion rate calculated from USGS erosion pin data.

Site	Current		After 10 Years
	Monitored Erosion between 2004–2009 (m)	FS at 57 m3/s	Predicted Cumulative Erosion (m)	FS at 57 m3/s
Site 1	0.92	1.92	2.3	1.58
Site 2	0.22	1.84	0.55	1.83
Site 3	0.2 1	2.07	1.0	1.83
Site 4	N/A	1.74	1.2 2	1.61
Site 5	0.33	1.80	0.83	1.81

¹ 2007–2009; ² Assumed.

show, Sites 1, 3, and 4 are predicted to become more unstable as erosion progresses. The unstable progression of erosion due to the initial bank geometry and predicted erosion mainly originates from the lower sections of the slopes. In contrast, Sites 2 and 5 are expected to maintain a similar FS, primarily due to their mild slopes and very low erosion rates.

5.6. Discussion

In numerical modeling, input parameters are the most critical factors along with the geometry and boundary conditions. Even though soil samples were carefully collected and tested in the field and lab, considering the unsaturated soil conditions, the testing method and interpretation of the results can affect the stability analysis. Direct shear tests usually measure lower friction angles than triaxial tests for anisotropic clays when the failure plane is horizontal, which are common in riverine deposit [1]. In addition, even though practical advantages exist in the multi-stage loading, it may result in reduced shear strength parameters. Therefore, the factors of safety in this study could have been underestimated, which means the actual slope could be more stable than analyzed. In addition, the unsaturated shear strength angle (${\varphi }^b$) showed non-linear behavior by having larger values in the lower suction range [18,46], which might have resulted in underestimation of the factor of safety at lower suction ranges in this study. In this case, the actual slope could be more stable than analyzed as well.

The slope stability analyses were conducted using a linear Mohr-Coulomb failure criterion, whereas the actual failure envelopes for most soils are non-linear, especially at low normal stresses. Jiang et al. [47] concluded that use of the linear approximation of a non-linear strength envelope resulted in significant overestimation of the calculated factor of safety in slope stability analysis, especially where the failure surface was shallow.

Regarding the riverbank stability of the study sites, as the steady state and transient seepage cases are shown in Figure 10, Figure 13 and Figure 14, the factors of safety at Sites 1, 3, and 5 changed significantly as the WSE rose or dropped, whereas those at Sites 2 and 4 did not vary significantly, although they were initially less stable. This can be associated with the sizes of the critical failure envelopes in the models, which were smaller at Sites 2 and 4. This seems to be attributed to the combinations of the bank geometry, soil properties, and WSE. The lower CL layers, in particular, had higher $c^{\prime }$ and ${\varphi }^{\prime }$, and larger thickness within the ranges of the WSE fluctuation at Sites 2 and 4, while the unsaturated soil properties were not much different.

Most slope stability analyses provided good results in terms of the safety factor. However, a few issues such as inconsistent test results, unrealistic body force estimates from equilibrium equations, or unrealistic values of factor of safety were identified with lower values of hydraulic conductivity.

Both transient drawdown and step-down scenarios using the lowest hydraulic conductivity values resulted in less stable conditions. Most of the failure envelopes seemed to be located within the range of normal stresses in the direct shear tests. However, some of the calculated factors of safety values were deemed to be invalid, containing numerical errors that implied a rapid drop in safety factor that may not represent the actual stability of slopes. These errors may occur when the calculated effective normal stress is negative. Among the potential causes for these errors are high pore pressure in the soil, due to the lower hydraulic conductivity, and a steep base angle of the slices in the limit equilibrium method (LEM) [48]. Poor drainage of the soils with lower hydraulic conductivity means that they experience larger excess pore water pressure for some time after a drawdown event, resulting in negative values of normal stress calculated at the bottom of each slice. Consequently, the factor of safety becomes unrealistically low.

Another possible source of error is the development of tensile strength due to cohesion. In cohesive soils, as the tension develops, it induces negative normal forces. The reversal of interslice normal and shear forces due to the cohesion is significant compared to the magnitude of the other forces. These negative interslice forces also produce discontinuities in the line of the thrust. This seems to be especially true for relatively shallow slices in soils, where cohesion dominates the shear strength of soils [49].

Typically, it is recommended that a tension crack be created to eliminate such negative values. However, negative normal forces may be calculated even below the tension zone, and this situation frequently occurs with relatively shallow slip surfaces and high cohesion. In addition, matric suction in unsaturated soil also increases tensile strength values, which ultimately increases the depth of the assumed tension cracks [49]. In low profile slopes like riverbanks, deep tension cracks may not be feasible, as increasing the crack depth decreases the factor of safety significantly, shortening the length of resistible slip lines and increasing the weight of the soil mass. Dry tension cracks may reduce the rate at which FS decreases, but it may still be necessary to eliminate portions of the slip line.

6. Conclusions

The actual reservoir release rates of a hydropower dam were first analyzed to identify representative flow patterns. Then transient seepage and slope stability analyses were performed to highlight the influence of dam operations on downstream riverbank stability using input parameters determined by laboratory tests, in situ tests, and empirical methods. The seven steady-state flow rates modeled represented the range of seasonal and operational modes throughout a typical year and included one overbank flow case. The model provided a valuable picture of the influence of water surface elevation (WSE), unsaturated shear strength parameters, and hydraulic conductivity.

These analyses revealed that the riverbanks at the study sites were stable for all their present flow conditions, including steady-state, peaking, drawdown, and step-down conditions. Considering the unsaturated soil characteristics, the slopes become less stable with wetting events because of a loss of matric suction. With a steeper slope and shallow slip surface, the riverbanks could be more dependent on the changes of the matric suction. However, at the study sites, even if the reduction of shear strength was expected because of wetting, it still seemed to be less critical than the counter effect by rising WSE. For example, a flow rate of 991 m³/s would result in overbank flow throughout much of the lower Roanoke River. Even here, though, the increase in confining pressure would most likely result in the banks remaining stable.

In addition, the banks were predicted to be less stable if they consisted of less permeable soils due to excess pore water pressures. The hydraulic conductivity determined by the laboratory tests tended to be considerably lower than that specified by the in situ tests. Although the lab values were close to the typical values for similar soil types, the in situ hydraulic conductivity represented the flow of water in the riverbanks better, as seen from the GWT modeling results, and can thus be used for modeling. However, as noted earlier, if low permeability soils are present, this may create more unstable conditions, especially when the riverbanks are exposed to transient drawdown or step-down flow conditions. Areas of low permeability soil may exist in the riverbanks, which likely explains the occasional small-scale bank failures observed along the study’s reach.

Although the major analyses conducted for this study provided good results, several numerical issues were identified. These are related to the unrealistic values for the factor of safety and are possibly due to the negative normal stress calculated at the slip line. Cohesion seems to be responsible for the negative values. Thus, the application of unsaturated shear strength may increase the possibility of such errors as the matric suction is regarded as a part of total cohesion. Further investigation of the role and influence of unsaturated shear strength is required for limit-equilibrium slope stability analysis.