Stability of boulders on cobble streambeds

The placement of boulders in streams enhances aquatic habitat by increasing the heterogeneity of flow conditions. Practical design must ensure the stability of individual boulders, requiring calculation of their incipient movement conditions. The stability of a boulder on a cobble bed is shown experimentally to depend on its size, the bed material size, the degree of embedment of the boulder, the slope of the bed and the flow velocity and depth. An equation is derived through a pivoting analysis for predicting the relationship between the critical ambient depth-averaged velocity and the critical flow depth; this, together with a resistance equation, can be used to predict the flow conditions for boulder stability. The equation is used to develop a simpler form for unsubmerged spherical boulders and cobbles. The stability equation is tested against the experimental data, using an experimentally determined drag coefficient relationship.


INTRODUCTION
The placement of rocks or boulders is common practice in stream stabilization and rehabilitation projects.Boulders increase resistance (James, 2021), thereby reducing scouring velocities and inducing sediment deposition (Dermisis and Papanicolaou, 2008).Sparsely placed boulders can have a stabilizing effect on the channel as a whole and are usefully employed in gradient-controlling rock ramps (Pagliara and Chiavaccini, 2006).
Boulders have been found to be an effective and low-cost alternative solution for river rehabilitation (Rutherfurd et al., 2000;Van Zyll De Jong et al., 1997;Dolinsek et al., 2007).The heterogeneous flow conditions around boulders provide suitable habitats for various life stages and activities of biota (Harvey and Clifford, 2009).The low-velocity wake zones provide important refuge areas for invertebrates and fish for resting and feeding, as well as enhancing the resilience of biotic communities to sudden floods (Engström et al., 2009;Huusko and Yrjänä, 1997;Merz et al. 2004).The turbulence and vortices in these zones also create cover for fish from visual predators by diffusing sunlight (Dermisis and Papanicolaou, 2008) and inducing local scour (Fischenich and Seal, 1999).Riverine Salmonid species use flow obstructions (including boulders) as velocity shelters to minimize exertion and save energy while migrating upstream; boulders are therefore valuable features in natural fishway design.
Various recommendations have been published for the use and placement of boulders for habitat enhancement (Fischenich and Seal, 1999;Rutherfurd et al., 2000;Schueler and Brown, 2004).These give guidance for the size, location, spacing and distribution pattern of the boulders.As an aid for habitat definition and rehabilitation design, Hamuy-Blanco and James (2014) and Heyneke (2019) described the 'zone of influence' around an individual boulder or boulder cluster with quantitative descriptions of the variations of velocity upstream, downstream, and laterally away from a boulder.
Apart from providing the required habitat conditions, the boulder size must be selected so as to also ensure that it remains in place and effective under anticipated flow conditions.Failure may occur by the boulder moving out of position, or by becoming deeply embedded if the cobbles surrounding it are excessively scoured.Various approaches are available for predicting incipient movement conditions for sediments.Probably the most common is the Shields criterion (Shields, 1936), which is incorporated in the moment-equilibrium based factor of safety calculation recommended by Fischenich and Seal (1999) for boulders.Such criteria are intended for application to relatively uniform-sized and well-submerged bed material rather than for the considered case of a single large particle on a bed of very much smaller ones.They are also based on comparison of the shear stress on the bed with a critical value, which is unrealistic for situations where the flow depth is comparable with the size of the potentially moving particle and the force actually moving the particle is not well described by the bed shear.A stability criterion in terms of velocity is therefore preferable.Erosion criteria in terms of velocity have previously been proposed by, for example, Novak (1948, cited in Vanoni, 1975), Neill (1968) and Prakash (2010).A method is presented here for predicting the critical ambient depth-averaged velocity and flow depth for incipient movement of a boulder on a cobble bed.Failure by self-embedding is not considered.

Analytical model
It is assumed that at the condition of incipient motion the boulder will begin to move by pivoting over an axis formed by contact points with supporting particles (as by White, 1940;Bagnold, 1941;Slingerland, 1977;Komar and Wang, 1984;Komar and Li, 1986;James, 1990).The development here follows the analysis of James ( 1990), but with the prediction made in terms of the depth-averaged velocity rather than the boundary shear stress.
The situation analysed is of a boulder with nominal diameter D resting on bed particles with nominal diameter K (Fig. 1).The centroid of the boulder is a distance μD above the centroids of the supporting particles, and the distance between the centroids of the two supporting particles defining the pivot axis is λK.The boulder may be embedded with its bottom a distance E below the tops of the supporting particles.
At incipient motion the forces on the boulder will include its submerged weight (W), friction (F r ) and reaction (R) forces only at points of contact with the supporting particles defining the pivot axis, and the hydrodynamic lift (F L ) and drag (F D ) forces (Fig. 1).The moments of the disturbing and stabilizing forces are in equilibrium: in which ϕ is the pivot angle and a, b and c are the distances from the pivot axis to the boulder centroid, line of action of the drag force, and line of action of the lift force, respectively.
The submerged weight of the boulder is its unsubmerged weight less the weight of water it displaces, which depends on whether it is completely or partially submerged:.
in which g is gravitational acceleration, V disp is the volume displaced, ρ s is the boulder material density, ρ w is the density of water, and k 1 is a proportionality constant for the relationship between the boulder volume and the cube of its diameter (= π/6 for a sphere).
The pivot angle (ϕ) depends on the sizes of the boulder and the supporting particles, their spatial arrangement, the slope of the bed (α), and the embedment (E) of the boulder.
The distance between the boulder centroid and the pivot axis, a, is assumed to be proportional to the size of the boulder: with k 2 depending also on the size and arrangement of the supporting particles.
The drag force is given by: in which C D is the drag coefficient, A n is the boulder's exposed cross-sectional area projected in the flow direction, and V cr is the ambient depth-averaged velocity at incipient movement.The drag force is assumed to act through the centroid of the exposed area of the boulder projected in the flow direction.The lever arm, b, is then given by: in which k 3 depends on the sizes of the boulder and supporting particles, the spacing of the supporting particles, the embedment of the boulder and whether the water surface is above or below the top of the boulder.
The lift force, assumed to occur only when the boulder is completely submerged, is given by: in which C L is the lift coefficient and A p is the boulder's crosssectional area projected normal to the flow direction.The lift force is assumed to act through the centroid of the area of the boulder projected normal to the flow direction.The lever arm, c, is then given by: Substituting the boulder weight (Eq.2), centroid-pivot distance (Eq. 3) and the drag and lift force expressions (Eqs 4 and 6) into Eq. 1 yields an equation for the critical velocity at incipient motion: Further manipulation and introduction of the submerged specific weight (g(ρ s − ρ w )) leads to an expression for the dimensionless critical velocity: with S s = ρ s /ρ w .Equation 9 provides a relationship between the critical velocity and the critical flow depth (y cr ) at incipient motion, the flow depth being implicit in V disp , C D, A n and b.A second relationship between depth and velocity is therefore required for its solution, which is provided by a conventional resistance equation, such as the Manning equation.

Experimental investigation
Two laboratory investigations were carried out.The first was to measure the incipient movement conditions for a single sphere resting on a bed of spheres in order to verify Eq. 9.The second was to determine the drag coefficient for the exposed sphere to provide the necessary value for applying Eq. 9.

Incipient motion experiments
The incipient motion experiments were carried out in a 12.0 m long, 2.00 m wide, tilting flume with a working width of 0.90 m.The discharge (Q) was controlled by a computer-operated valve in a line leading from a constant head tank, and measured by an inline flow meter as well as a V-notch downstream from the flume.Flow depth (y) was controlled by a weir at the end of the flume, and measured on several scales along its length.The flume bed was lined with alluvial pebbles with a mean diameter of 25 mm, and a patch of 25 mm spheres was placed in the immediate vicinity of the test sphere.
Test spheres were placed in a recess halfway along the flume.A surface within the recess was made with 'Lego' toy building blocks to enable the sphere mounting level to be adjusted.Test spheres were placed initially on three 25 mm supporting spheres in closest packing, with the front two aligned perpendicular to the flow direction.Further tests were done with the test sphere embedded at a lower level and with the front two supporting spheres spaced 31.3 mm apart, centre-to-centre.
Test spheres with different diameters and densities were used, with a number of different flume slopes, as listed in Table 1.The different densities were obtained by filling the hollow steel 80 mm sphere with different substances.
The experimental procedure was to place a boulder sphere in position and to increase the discharge in small increments for 5-min durations until the sphere was dislodged.The sphere was then removed from the flume and the uniformity of the flow checked by measuring the flow depth at 3 locations along the length of the flume.If the flow was not uniform, the weir was adjusted to make it so.The flume was then drained, the sphere placed back in position, and the process repeated with the correct weir setting so that the dislodging occurred with uniform flow.The ambient depth-averaged critical velocity was then calculated by dividing the discharge by the cross-sectional flow area.
Results are presented in Table 1, with the critical flow depth (y cr ) measured from the theoretical bed, assumed to be at the level of the centroids of the supporting particles (James, 1990).
The results show that the stability of the boulder depends on its size (D), density (ρ s ), embedment (E) and the longitudinal bed gradient (S o = tan α).(The supporting particle size (K) was constant in the experiments.)The incipient motion condition is defined by the critical velocity (V cr ) and the critical flow depth (y cr ) pair.Because each measured critical velocity occurs with a different critical depth, the results are interpreted in terms of the critical unit-width discharge (q cr = V cr y cr ) rather than simply V cr .
Data for assessing the effect of boulder size on stability are provided by the results of the experiments with D = 80 mm, 100 mm, 120 mm and 150 mm, with S o = 0.00645 and E = 0.The densities of the different-sized boulders were not identical, ranging from 1 558 kg/m 3 to 1 682 kg/m 3 , but this range was shown to be associated with a variation in q cr of less than about 6%.The boulder stability depends on the size of the bed particles as well as of the boulder itself.As the ratio D/K increases the pivot angle reduces, tending to decrease stability, but for D/K > ~1 this tendency is more than offset by the greater boulder weight as D increases.
For the fixed value of K in the experiments, the critical discharge increases considerably with D (Fig. 2).
The effect of embedding is shown by the experiments with D = 80 mm and ρ s =1 645 kg/m 3 , for which 5 embedded conditions were used.Two embedded conditions were also used for experiments with D = 80 mm and ρ s = 1 877 kg/m 3 , D = 100 mm and ρ s = 1 591 kg/m 3 , and D = 120 mm and ρ s = 1 682 kg/m 3 .S o was the same (0.00645) for all these experiments.As shown in Fig. 3, the effect of embedding is considerable: the critical discharge is doubled by embedding the boulder by just a tenth of its diameter.
The effect of channel slope on stability is shown by the experiments with D = 80 mm and ρ s = 1 645 kg/m 3 , 1 877 kg/m 3 and 2 335 kg/m 3 , and with D = 100 mm and ρ s = 1 591 kg/m 3 .Slopes tested ranged from 0.00134 to 0.0071.As shown in Fig. 4, the critical discharge is quite insensitive to channel slope over the range tested.
The effect of boulder density on stability is shown by the experiments with D = 80 mm with S o = 0.00391, 0.00519 and 0.00647, and with ρ s = 1 645 kg/m 3 , 1 877 kg/m 3 and 2 335 kg/m 3 .
Figure 5 shows that density increases stability (similarly for all channel slopes), but the effect would be relatively minor for the range of densities experienced in practice.

Drag coefficient experiments
The experiments to determine the drag coefficient on the exposed spheres were carried out in a 0.90 m wide horizontal flume.The bed was prepared in the same way as for the incipient motion experiments, except that the test sphere did not rest on supporting spheres but was suspended just clear of them.The drag force was measured by means of moment equilibrium.The test sphere was attached to a rod, which pivoted about a pin, inducing a force measured on a tension sensor placed above the pin.Details of the apparatus and its calibration are described by Stols (2018) and Jackson (2016).Discharge (Q) was measured using the inline flow meter as well as volumetrically in the sump at the end of the flume; the water level was controlled by a weir at the end of the flume and measured at the test location using a digital point gauge.Twelve boulder conditions were tested, 10 for an 80 mm diameter sphere, including two embedded cases, and two for a 120 mm diameter sphere.Relative depths ranging from 0.65 to 1.83 were tested.
The conditions are listed in Table 2.The sphere Reynolds number (Re = VD/ν, where V is the ambient depth-averaged velocity and ν is the kinematic viscosity of water) was in the range 9.2 x 10 4 to 1.65 x 10 5 , within which C D would be expected to be practically independent of Re.
For each test, a target water depth was set on the point gauge.
The discharge in the flume was then slowly increased to a desired value and once this was stable the weir was adjusted to produce the target water depth at the test location.This condition was maintained for 5 min for the force sensor data collection.The point gauge was then lowered to a new target level and the weir readjusted accordingly.
Values of C D (as listed in Table 2) were calculated from the measured force by Eq. 4 with A n calculated by Eqs 21 to 25.The flow depth (y) was defined as the distance from the water level to the theoretical bed (K/2 below the tops of the bed particles).The relationship between C D and relative depth is shown in Fig. 6 with R 2 = 0.65.The data are limited, and the embedded and non-embedded cases are combined for defining this relationship although they appear to be slightly distinct.It is also assumed that for high submergences the value of C D will level off (Coleman (1972) found C D for well-submerged spheres to be approximately equal to the value in free fall); a minimum of 0.77 is suggested by the highest non-embedded y/D value, which was not included in the regression, and is assumed to apply for y/D > 1.5.

Model confirmation
The experimental results are used to assess the performance of Eq. 9 for predicting the depth-averaged velocity at incipient motion.
The geometry of the experimental conditions with a spherical boulder and supporting particles is well defined, enabling the variables to be described precisely (Fig. 7).
For a spherical boulder k 1 = π/6.The volume of water displaced by the boulder depends on whether it is completely or only partially submerged.For the submerged condition: and for the unsubmerged condition: in which h is the distance of the top of the boulder above the water surface, given by: with E being the embedment and y the flow depth measured from the theoretical bed.
The factor k 2 , defining the distance (a) between the boulder centroid and the pivot axis (Fig. 1), depends on the spacing of the two supporting particles forming the pivot axis, λK.This is given by: λ = 1 for the non-embedded condition and 1.252 for the experimental embedment.
The pivot angle (ϕ) (Fig. 1) (here in radians) depends on the supporting particle arrangement as well as the longitudinal bed slope, as: The distance μD (the height of the boulder centroid above the supporting particle centroids) is given by: For a boulder resting on 3 particles in contact (λ = 1) the embedment is a minimum, given by: For λ ≠ 1, E is specified as required.
The drag force is calculated using Eq. 4 with the drag coefficient given by Eq. 10.The projected area for drag is the total projected area of the spherical boulder less the area below the theoretical bed (A E ) and less the area above the water surface (A s ), if applicable (Fig. 7): The subtracted area A E is given by: θ E is the angle subtended at the boulder centroid by the width at the The area A s is given by: θ s is the angle subtended at the boulder centroid by the width at the water surface, and is given by: The drag force lever arm is given by Eq. 5.The distance k 3 depends on whether the boulder is partially or completely submerged and the degree of embeddedness.For a submerged boulder it is given by: where V 1 is the volume of the boulder below the theoretical bed (at K/2 below the tops of the supporting particles), given by: with h / being the distance of the bottom of the boulder below the theoretical bed according to: z 1 is the distance of the centroid of V 1 below the centroid of the boulder, given by: For an emergent boulder: where V 2 is the volume of the boulder above the water surface, given by: z 2 is the vertical distance between the centroids of the whole boulder and its portion above the water surface, given by Eq. 26 without the minus sign and with h / substituted by h.
The lift force is calculated using Eq. 6. Lift coefficient values have not been well established (James, 1990;Marsh et al., 2004).It is assumed that lift is zero if the boulder is not fully submerged (as by Carling et al., 2002).For large, fully submerged particles C L has been found to be fairly constant; Marsh et al. (2004) recommend a value of 0.2.The lift coefficient has also been related to the drag coefficient (see James, 1990) and for large exposed particles different values of C L /C D have been proposed; for example, 0.1 by Aksoy (1973), 0.5 by James (1990) and 0.2 by Carling et al. (2002).These must depend on the definition of the velocity used in Eq. 6, which is not always clear.For the experimental conditions, good agreement was obtained using C L /C D = 0.5 in terms of the ambient depth-averaged velocity.The projected area for lift is πD 2 /4, and the lever arm is: Equations 11 to 29 are used to represent the experimental conditions for calculating the dimensionless critical velocity by Eq. 9 for each measured critical flow depth.
Predicted values of V cr * are listed in Table 1 and compared with the measured values in Fig. 8.The average absolute prediction error is 5.3%.

Practical stability assessment
Practical application of Eq. 9 would be laborious, requiring approximation of the geometric parameters.A single equation involving the key parameters was therefore developed by regression of synthetic data generated by Eq. 9 for a range of conditions.A sphere with a diameter of 0.40 m and a density of 2 650 kg/m 3 was considered.The size of supporting spheres (K) was selected to give a range of relative sizes (D/K) from 4 to 20.Flow depths (y) were selected to give an upper limit of y/D equal to 1.0 (the change in V cr * was found to be small once the boulder was submerged).Very shallow flow depths (y/D − 0.5K/D < 0.15) were for y ≤ D (for y > D, using y = D -E + K/2 would be a reasonable approximation).Predictions by Eqs 9 and 30 are compared in Fig. 9. Equation 30 reproduces Eq. 9 with R 2 = 0.98, and an average absolute difference of 8.6%.
Both Eqs 9 and 30 assume spherical boulders and supporting particles.Departures from sphericity would tend to increase stability by interlocking and increasing frictional resistance.The drag coefficient equation (Eq.10) also applies to spherical shapes, and values would be higher for angular boulders, tending to decrease stability.For emerging cylinders, Jackson (2016) found C D to increase by about 50% for square compared to circular shapes, and more than 100% for diamond shapes.
Assessment of stability through Eqs 9 or 30 requires simultaneous solution with a resistance equation, such as Manning's, because of the dependence of both on flow depth.The stability of the boulder also depends on the stability of the supporting particles.If the supporting particles can be eroded at conditions lower than able to move the boulder, then the boulder is likely to become further embedded and therefore more stable, although less effective for modifying flow conditions.The scour potential of the supporting particles can be assessed approximately by the Shields criterion, but may be enhanced by vortex formation around the boulder (as for scour around bridge piers).
As an example, the stability of a 0.55 m boulder embedded by 0.010 m in a bed of 0.050 m cobbles with a channel slope of 0.0050 is considered.The relationship between V cr and flow depth (y) is described by Eq. 9.The relationship between the actual velocity (V) and y is described by the Manning equation, with n = K 1/6 /26, and assuming the channel to be wide.The two relationships are plotted in Fig. 10, showing that the boulder will be stable for all flow depths less than 0.43 m, when V < V cr .It also shows that V cr decreases with depth towards a constant value at full submergence.It is possible that self-embedding could occur if the cobbles are able to be scoured at lower flow depths.There appear to be no methods for predicting local scour around boulders, but bridge pier scour predictions can give a rough indication by assuming the boulder to be a cylinder with the same diameter.Chiew (1995) proposed a procedure for estimating the approach velocity at which scour around piers would occur (V o ).The relationship between V o and y according to this procedure is also plotted in Fig. 6, suggesting that scour, and hence self-embedding, would occur for depths greater than 0.37 m, when V > V o .This is a conservative estimate, as scouring vortices would be undeveloped with flow occurring below part of the boulder.

CONCLUSIONS
The stability of an individual boulder on a cobble bed depends on the size and density of the boulder, the size of the cobbles, the embedment of the boulder in the cobbles, the channel slope, and the flow velocity and depth.The critical velocity for incipient movement reduces with increasing flow depth, reaching an effective minimum at the top of the boulder.
Stability can be predicted through a classical pivoting analysis.
Equation 9 provides the relationship between the critical dimensionless velocity and the critical flow depth (implicit in the boulder submergence condition).Equation 30 provides a simpler, direct relationship for unsubmerged spherical boulders on spherical cobbles.Practical predictions require a second relationship between velocity and flow depth, which is provided by a resistance equation such as Manning's.
The drag coefficient on a boulder depends strongly on its degree of submergence.The relationship presented by Eq. 10 is provisional, and further investigation for the transition to full submergence and the influence of boulder shape is necessary.
The stability of a boulder depends strongly on its size, and is increased considerably by even small degrees of embedding.
Prediction capability has been confirmed by experiments with spherical boulders and cobbles.Non-sphericity would tend to increase stability through particle interlocking and increasing friction, but reduce it through an increase in the drag coefficient.
Self-embedding is likely if the bed material around the boulder erodes at flow conditions lower than would move the boulder.This occurrence needs further investigation.

AUTHOR CONTRIBUTIONS
K Stols: Experimental methodology and design; all data collection; initial analysis and interpretation of results.
CS James: Conceptualization of the study; substantial re-analysis and interpretation of data; extension of results application; manuscript writing.

Figure 1 .
Figure 1.Schematic diagram of boulder on cobble bed

Figure 2 .Figure 3 .
Figure 2. Effect of boulder size on stability

Figure 4 .Figure 5 .
Figure 4. Effect of channel slope on boulder stability

Figure 7 .
Figure 7. Geometry of spherical boulder on spherical supporting particles

Figure 8 .Figure 9 .
Figure 8.Comparison of predicted and measured dimensionless critical velocity

c
A E Cross-sectional area below theoretical bed of boulder normal to flow direction A n Exposed cross-sectional area of boulder in flow direction A p Cross-sectional area of boulder normal to flow direction A s Cross-sectional area above water surface of boulder normal to flow direction a Distance from pivot axis to boulder centroid b Distance from pivot axis to drag force line of action Distance from pivot axis to lift force line of action