The term “breaching” refers to the slow, retrogressive failure of a steep subaqueous slope, so forming a nearly vertical turbidity current directed down the face. This mechanism, first identified by the dredging industry, has remained largely unexplored, and yet evidence exists to link breaching to the formation of sustained turbidity currents in the deep sea. In this paper we model a breach-generated turbidity current using a three-equation, layer-averaged model that has as its basis the governing equations for the conservation of momentum, water, and suspended sediment of the turbidity current. In the model, the turbidity current is divided into two regions joined at a migrating boundary: the breach face, treated as vertical, and a quasi-horizontal region sloping downdip. In this downstream region, the bed slope is much lower (but still nonzero), and is constructed by deposition from a quasi-horizontal turbidity current. The model is applied to establish the feasibility of a breach-generated turbidity current in a field setting, using a generic example based on the Monterey Submarine Canyon, offshore California, USA.

Turbidity currents are sediment-laden underflows that occur in oceans and lakes. These currents are driven downslope by the density excess arising from the presence of suspended sediment in the flow. Not only are these currents responsible for most submarine canyons on continental shelves and slopes (Pantin, 1979; Parker, 1982; Imran et al., 1998), they are also important conduits delivering sediment to deeper waters. While active turbidity currents have been monitored in freshwater lakes since the time of the work of Normark and Dickson (1976a, 1976b) and Lambert et al. (1976) (see also e.g., Chikita, 1989; Best et al., 2005), they have proved much more difficult to measure in the field (Inman et al., 1976; Hay, 1987a, 1987b; Xu et al., 2004). This notwithstanding, current properties can often be inferred from evidence left on the ocean and lake floors. For example, it has been inferred that these currents are capable of traveling at surprisingly swift velocities, e.g., as high as 8–14 m/s (Krause et al., 1970). The size of the submarine canyons (the Zaire Canyon is ∼15 km wide and as deep as 1300 m; Babonneau, 2002) underlines the erosive power of these currents.

Turbidity currents are initiated through a number of triggering mechanisms. Hyperpycnal flows are perhaps the most commonly identified source of sustained turbidity currents in lakes and reservoirs. In oceans, however, with the notable exception of the Yellow River in China, these hyperpycnal flows are for the most part restricted to active margins where there is sufficient sediment supply for the river water entering the sea to be heavier than the sediment-free seawater (e.g., Mulder and Syvitski, 1995). The identification of self-accelerating turbidity currents (Pantin, 1979; Parker, 1982) supports the possibility of multiple triggering mechanisms appropriate for active and passive margins. Sustained turbidity currents need not be produced by sufficient sediment supply; neither need they be produced by continuous, sustained events. Under the right conditions, even a small turbidity current plume can accelerate and entrain sediment from its bed in a self-reinforcing cycle, so strengthening and elongating (e.g., Sequeiros et al., 2009). A notable case is Scripps Submarine Canyon (offshore California, USA), where the field research of Inman et al. (1976) suggested that nearshore sediment plumes suspended near the head of the canyon by storms are responsible for sustained turbidity currents in the canyon. The numerical model by Fukushima et al. (1985) indicates that such conditions are capable of producing upstream conditions sufficient to ignite a self-accelerating current. Other triggering mechanisms also discovered or suggested include suspension of material from submarine debris flows, breaking internal waves, slope failure, and breaching in fine-grained sands. This last mechanism was presented in van den Berg et al. (2002); it is considered in more detail in the following.

Breaching as a Failure Mechanism

Breaching is a generic term that is often used to describe the process of retrogressive erosion and eventual failure subsequent to the overtopping of embankments, dams, dikes, and even sand barriers in closed estuaries (Coleman et al., 2002; Visser, 1998; Stretch and Parkinson, 2006). Here, however, breaching is used in a much more restrictive sense, i.e., that of van den Berg et al. (2002) and Mastbergen and van den Berg (2003), who referred specifically to a kind of relatively slow, subaqueous retrogressive failure that was apparently first identified by the Dutch dredging industry. The process is illustrated in Figure 1A: a dredge removes sediment from the base of a slope face that is retreating as it fails. The failure does not occur as a sudden collapse; instead, grains spall off one by one to form a turbidity current traversing a steep (in some cases nearly vertical) slope. Once this mode of failure is instigated, it can continue for many hours. This breaching mode of failure can in many cases be sustained even if the suction from dredging is halted, resulting in a continuous turbidity current propagating down the base of the slope.

The phenomenon of breaching may have escaped the notice of many partly because it is may often be mistaken for the more common liquefied slope failure. Both failures produce similar morphologies, occur under the same environmental conditions, and are often identified with fine sand. Van den Berg et al. (2002) sought to distinguish between the two failure patterns. While liquefaction tends to occur in loosely packed contractant sand (with porosity >43%), breaching is common with densely packed sands that tend to become dilatant. Dilatant sand undergoing shear deformation undergoes an increase in pore volume, leading to a negative pore pressure buildup (with respect to hydrostatic pressure) and an increase in the shear resistance (van Rhee and Bezuijen, 1998; van den Berg et al., 2002). The erosion process is thus severely retarded, leading to a gradual (in some cases, grain by grain) rather than sudden slope failure. These sand grains typically mix with the ambient fluid, forming a turbidity current that runs along the steep slope surface and down the gentler slope that forms at the base of the breach.

Liquefied slope failures, on the one hand, typically result in hyperconcentrated or slurry-like density flows, which may dilute farther downstream into a turbidity current. Breaching, on the other hand, results in a well-defined, dilute turbidity current from the upstream end. These observations may provide a basis for distinguishing between deposits from the two types of flow.

Mastbergen and van den Berg (2003) indicated that breaching failures of levees or deltaic deposits could occur naturally; for example, they quoted the work of Torrey (1995) on natural bank failures in the well-packed sand deposits of the Mississippi River in this regard. Mastbergen and van den Berg (2003) also provided evidence that the breaching of deltaic deposits likely constituted a triggering mechanism for sustained turbidity currents in the deep sea.

Breaching as a Mechanism for Generating Sustained Turbidity Currents

As noted here, the breaching process can be accompanied by the formation of a turbidity current. As the sand grains slowly peel off from the slope surface, they tend to form a flowing sand-water mixture with the ambient water. With the exception of a major research effort in the Netherlands, which was largely performed in secret for proprietary reasons (D. Mastbergen, 2002, personal commun.), turbidity currents generated by breaching have remained largely unexplored. Although the Dutch experiments were focused on the retrogressive failure associated with breaching rather than the turbidity current generation, van Rhee and Bezuijen (1998) observed turbidity currents in the laboratory with velocities as high as 1 m/s.

Mastbergen and van den Berg (2003) used this research as a tentative link between breach-generated turbidity currents in Scripps Submarine Canyon and canyon formation. Using a one-dimensional, depth-averaged formulation based on momentum, water and sediment continuity equations, Mastbergen and van den Berg (2003) developed a numerical model of a breach-generated turbidity current. The model incorporated a relation for erosion modified from the work of Winterwerp et al. (1992) to consider the effect of breaching. The model showed that with reasonable estimates of an initial sediment transport rate, wall velocity, sediment concentration, and Froude number, the breaching mechanism is capable of producing turbidity currents of sufficient strength to excavate a submarine canyon.

The three-equation, layer-averaged model proposed in this work derives from the work of Fukushima et al. (1985) and Parker et al. (1986). Like the model of Mastbergen and van den Berg (2003), it has a basis in the governing equations for the conservation of momentum, water, and suspended sediment of the turbidity current, but the three-equation model provides a more formal and complete treatment.

The present implementation is designed to capture at the same time both the turbidity current as it runs down the steep face of the breach and the same current on a much gentler slope farther downstream. To this end, the turbidity current is divided into two regions joined at a migrating boundary corresponding to the base of the breach face, i.e., the breach face and a quasi-horizontal region sloping downdip (see Fig. 1B). (The precise definition of quasi-horizontal is provided herein; however, here it suffices to note that a slope of 0.01 [5.71°] or less is within the definition.) In the present formulation, the breach face is treated as vertical for convenience. This assumption can be relaxed to include varying steep inclinations. Over the downstream region, the bed slope is much lower (but still nonzero), and is assumed to have been sculpted by either deposition from or erosion into a quasi-horizontal turbidity current. On the breach face, the rate of erosion of sediment into the vertical turbidity current is explicitly linked to the rate of retrogression of the breach face. In the turbidity current downstream, sediment is allowed to freely exchange with the bed via erosion and deposition. This condition may not in general be satisfied in the field, because canyons likely erode into shelf-slope material that has developed cohesive strength. The assumption of free exchange with the bed, however, allows capture of the relevant tendencies.

A similarity solution of a set of equations based on the three-equation model of turbidity current dynamics and breaching yields a formulation for the vertical current running down the breach face. The flow computed from this solution at the breach face serves to specify the upstream boundary condition for computation of the more nearly horizontal, down-canyon turbidity current beyond using either the full three-equation model or the four-equation model of Parker et al. (1986).

The three-equation model is subject to limitations. It is capable of predicting the onset of self-acceleration, but typically overaccelerates farther downstream because it does not have within it a mechanism for preventing excessive sediment entrainment. Parker et al. (1986) partially resolved this problem with a four-equation model that includes a description of damping of turbulence due to density self-stratification. This mechanism does not operate in the case of vertical or very steep flows, so the use of the three-equation model on the breach face is justified. In the down-canyon region below the breach face, we have used the three-equation model solely for the purpose of predicting the onset of self-acceleration. In modeling actual field cases, we have used the four-equation model in down-canyon regions. The fourth equation, i.e., for conservation of turbulent kinetic energy, provides the necessary mechanism for limiting sediment entrainment (and hence self-acceleration), yielding more realistic results.

The present three-equation model has been tested and validated against experimental data (Eke et al., 2009; Eke, 2008). These experiments, unlike those of van Rhee and Bezuijen (1998), were designed to model breaching as a deep-water phenomenon, and not the predominantly shallow-water breaching observed during dredging for which the present numerical may not apply. The relevant aspects of the experiments and of the validation of the numerical model are briefly described and discussed.

Consider a steep subaqueous breach wall of height Zb(t) (where t = time), as shown in Figure 2. For simplification, we assume a constant, nearly vertical slope angle θv along the wall. As the breach wall fails by spalling, a turbidity current is generated that runs downstream, interacting with the wall and the channel bed surface below the base of the wall. This results in erosion of the breach wall and, depending upon conditions, net deposition or erosion on the bed downstream of the wall. Therefore, we consider the failure process to be divided into two separate but connected regions: (1) the nearly vertical region where the current of grains runs down the steep breach face; and (2) the down-canyon region where the current runs over a much milder bed slope. In the first region, the flow is oriented in the vertically downward z direction. The origin of the z coordinate is at the top of the breach face. At the base of the breach face [i.e., at z = Zb(t)], the flow makes a sharp turn, becoming quasi-horizontal in the sense that the bed slope angle θh satisfies the conditions
where S is the bed slope. The streamwise coordinate in the quasi-horizontal region of flow is denoted as x and the breach face is located at x = su(t), where su defines a boundary that moves with the retreating breach face and t is time.

For the case of sediment consisting of a mix of grain sizes, the active layer approximation (Hirano, 1972) is useful in describing the interaction between the current and the bed in the down-canyon region as deposition progresses. In such a formulation, sediment in transport exchanges directly with a surface or active layer of the bed of specified thickness. Sediment in the active layer then exchanges with sediment in the substrate below through bed aggradation and/or degradation (e.g., Parker, 1991). As grain sizes associated with breaching are usually within the fine sand to coarse silt range, it is assumed that essentially all the sediment moves in suspension, so that bedload transport can be neglected. The various regions and layers described above are depicted in Figure 2.

Each layer is assumed to be a mix of N grain size classes with characteristic diameter Di. The fraction of sediment in each grain size class is denoted as Fi and fi in the active layer and substrate layer, respectively. The active layer thickness La is assumed to be a linear function of the characteristic grain size D90 of the surface layer (i.e., La = 2D90; Parker, 1991), which may change in time as the flow interacts with the bed.

To describe the dynamics of the turbidity current flow, the layer-averaged three-equation model of Fukushima et al. (1985) and Parker et al. (1986) is extended to encompass nonuniform sediments. This model captures the interaction of the developing turbidity current with the ambient water and sediment bed. To capture bed and breach wall morphodynamics, the three-equation model is coupled with the Exner equation for bed sediment conservation, as shown in the following. In the coupling, the quasi-steady flow approximation is adopted, so that temporal variation is considered only in the bed elevation. Strictly speaking, this approximation applies only to Froude subcritical flow (De Vries, 1965), but it is a reasonable approximation in the absence of propagating shock waves or moving zones of transcritical flow.

Governing Equations for Turbidity Currents

The approximate layer-averaged three-equation model of Fukushima et al. (1985) and Parker et al. (1986) is derived from the basic conservation laws of momentum balance, fluid, and sediment mass balance. The dependent variables in the down-canyon region are the current thickness H, the layer-averaged velocity in the stream-wise direction U, and the layer-averaged volume suspended sediment concentration C; these same parameters are denoted as Hv, W, and Cv, respectively, for the nearly vertical flow on the breach face. The top-hat approximation is adopted, where concentration and velocity distributions are replaced by layer-averaged values (Turner, 1973). The volume concentration C is assumed to satisfy the condition C << 1 everywhere within the flow, so that the suspension is dilute (e.g., volumetric concentration measurements of turbidity currents in Lake Superior created by the disposal of mine tailings were measured to be ∼1.8 10−5; Normark, 1989). This justifies the assumption of a kinematic viscosity ν that is equal to the value for clear water (Garcia, 1990), as well as the use of the Boussinesq approximation (Fukushima and Parker, 1990). The density of the suspension can be taken to be equal to ρ(1 + RC), where ρ is the density of the ambient sediment-free water and R is the submerged specific gravity of the sediment in suspension (∼1.65 for quartz; Garcia, 1990). The slender flow approximations are also applied and the current is assumed to be fully turbulent, with all viscous terms negligible (Garcia, 1990). The general flow equations resulting from the above approximations are presented in the following.

Three-Equation Model for the Down-Canyon Current

For the quasi-horizontal turbidity current in the down-canyon region, the derivation of the layer-averaged flow equations from the conservation equations of momentum, fluid mass, and sediment mass were detailed in Parker et al. (1986).

Upon modification for multiple grain sizes, these equations are more conveniently expressed as a set of backwater flow equations. These equations take the form:
(e.g., Salaheldin et al., 2000), where S denotes the bed slope, cf denotes a dimensionless bed friction coefficient, ew denotes a dimensionless coefficient describing the entrainment of ambient clear water into the current, and Ci denote the volume concentration of sediment in the ith grain size range, which has characteristic size Di. In addition, vsi denotes the fall velocity of size Di (here computed using the relationship of Dietrich, 1982), Esi denotes a coefficient of entrainment of bed sediment into suspension for size Di, and roi denotes the ratio between the near-bed concentration of suspended sediment in the ith range and the layer-averaged value Ci. Ri is the dimensionless bulk Richardson number. For simplicity, the concentration ratio roi is assumed to be constant and equal to 2 (e.g., Garcia, 1990). The concentrations Ci are related to the total concentration Ct as
In addition, qi = UCiH denotes the volumetric suspended sediment transport rate per unit width in the ith grain size range,
denotes the volumetric suspended sediment transport rate per unit width summed over all grain size ranges, and
denotes the equilibrium value of qi at which neither erosion nor deposition would occur. Finally, Ri, the dimensionless bulk Richardson number, is given as

The sediment entrainment function, Esi, developed by Wright and Parker (2005a, 2005b) for sediment mixtures at field scale, is used here (see also Kostic and Parker, 2006; Fildani et al., 2006). The relation for the dimensionless water entrainment parameter ew is adopted from Parker et al. (1987).

Three-Equation Model for the Nearly Vertical Current on the Breach Face

For the nearly vertical turbidity current on the breach face, the governing equations differ somewhat from that of the down-canyon region. Here, the pressure gradient term is dropped from the momentum equation as a consequence of the assumption of a vertical face (an assumption that remains valid even for the case of a very steep but nonvertical face). In addition, the flow is considered to be purely erosive, as is appropriate for a vertical or nearly vertical face. The three-equation model thus reduces to the forms:
where Hv and W are the layer thickness and layer-averaged downward vertical velocity, respectively. In addition, Cvi and Cvt denote the layer-averaged suspended sediment concentration in the ith grain size class and the total concentration summed over all grain sizes in the layer, respectively, Ebi denotes the volume rate of spalling of sand in the ith size class from the breach face per unit area normal to the face, and cfb is the wall friction coefficient for the nearly vertical flow on the breach face.
Assuming that the vertical current thickness, velocity, and concentration can be represented as power functions of the vertical breach height z, we arrive at the following similarity solutions:
Hence flow parameters for the nearly vertical current can be computed directly without resorting to a numerical model.

The parameter Ebt denotes the total volume rate of sediment loss per unit area due to spalling from the breach face. It must be constant in z if the nearly vertical face is to be maintained as it retreats, and is related to the breach wall velocity vwall as Ebt = vwall (1 – λpr), where λpr is the porosity of undisturbed sediment reservoir material of the breach wall. Here Ebi is calculated as Ebi = EbtFvi, where Fvi denotes the fraction of material of grain size Di in the breach wall sediment. For simplification purposes, the water entrainment ewv is assumed to be constant and equal to the maximum value of 0.075, since the entrainment of clear water into the vertical turbidity current is not impeded by turbulence suppression due to self-stratification (in that the mixing occurs orthogonally to the vertical).

Equations 9, 10, and 11 evaluated at the point z = Zb corresponding to the base of the breach face provide the upstream boundary conditions on equations 1, 2, and 3, which describe the evolution of the turbidity current on the quasi-horizontal region (see Fig. 3B).

Governing Equations for Bed Evolution

Change in Bed Elevation

In a complete model of turbidity current morphodynamics, it is also necessary to track the evolving bed elevation in terms of the Exner equation of conservation of bed sediment. This relation in its one-dimensional form is given as
where η denotes the bed elevation, t is time, and λp denotes the porosity of the sediment deposit.

Change in Grain Size Fractions

Solution of the following equations for each of the grain size classes provides the time rate of change in the grain size fractions Fi in the active layer (Parker, 1991):
In equation 13, fIi denotes the grain size fractions at the interface between the active and substrate layers, here given as
where α is a value between 0 and 1. The value of 0.5 was used in our analysis for lack of better information. More details about this formulation can be found in Hoey and Ferguson (1994) and Toro-Escobar et al. (1996). The grain size fractions in the substrate fi of the down-canyon region are assumed to be varying in space but constant in time.

Change in Breach Height

Consider the diagram of Figure 3A. The sediment surface above the top of the breach face has (fixed) elevation ηt at x = 0 and a fixed slope St. The bed elevation η at the toe of the breach is thus given as:
Noting that the upstream end of the down-canyon reach must migrate upstream at the same speed as the breach wall, the following condition must be satisfied:
Taking the derivative of both sides of equation 15 with respect to time and employing equation 16,
Note that Su denotes the bed slope of the quasi-horizontal region at the base of the breach face.

Numerical Modeling Procedure

The numerical model presented here involves a grid system that changes in time. The nodes are fixed so the grid spacing (∆x) does not change. However, as the breach migrates updip and the domain length increases, a node is added upstream. The position of this node is allowed to shift in time until the grid spacing for the new node is ≥∆x. At this point a new node is introduced (Fig. 4). This system was adopted to avoid the interpolation problems inherent in a moving boundary coordinate system.

The input parameters required include the top slope St above the breach face, the initial breach height Zbi, the initial horizontal position of the breach face sui, the initial bottom slope Si below the base of the breach face, the initial reach length Li, the grain size fractions of the vertical wall, the bed fraction redistribution coefficient α, the friction coefficient of the bed and breach surface cf and cfb, respectively, the porosity of the bed and breach wall λp and λpr, respectively, and the wall velocity vwall.

The numerical solution procedure is as follows.

  1. The similarity solution for the vertical current is first computed.

  2. The flow values of Hv, Cvi, Cvt, and W at the toe of the vertical face are then input as the upstream boundary conditions for the quasi-horizontal current. The assumption here is one of complete transfer of flow properties from the vertical region to the horizontal region at the toe of the breach wall.

  3. A simple Euler time-stepping scheme is adopted for integrating the ordinary differential equations 1–7, 12, 13, and 17 describing flow in the down-canyon region. The new bed slope, breach height, and fraction of each individual size class in the active layer Fi are then updated in the flow equations.

The single-layer nature of this model prevents its use in shallow water. That is, it can be used only for breaching associated with systems that are sufficiently deep for turbidity current thickness to be a small fraction of water depth. The present model is capable of treating as many as 10 different grain size ranges. It has been coded in Visual Basic, and placed in the repository of the Community Surface Dynamics Modeling System as “1DBreachingTurbidityCurrent” (

Experimental Validation

Small-scale experiments were conducted at the Hydrosystems Laboratory, University of Illinois at Urbana-Champaign, to investigate the possibility of breach-generated turbidity currents. In previous experiments carried out by van Rhee and Bezuijen (1998), the failure was instigated by a sustained suction with a tube, a configuration designed to mimic the effect of dredging. In our work, however, the turbidity current is allowed to develop freely from a breach face that fails naturally under the effect of gravity, and flows downslope under its own strength. This configuration was intended to resemble more closely what might occur at the head of a submarine canyon. In the experiments, breaching from a nearly vertical face was achieved by removing a vertical wall between the sediment and the ambient water. Figure 5A shows a schematic of the experimental set-up and Figure 5B shows an experiment in progress. (For details on the experimental procedure and results, see Eke, 2008; Eke et al., 2009.)

A validation of the three-equation model with experimental results was performed (Eke et al., 2009) for the case of uniform sediment. For the case of sediment mixtures, the major calibration parameters are the concentration ratio roi for each grain size range and the sediment redistribution coefficient α. The parameter roi was calibrated by comparing the bed profile predicted by the model with experimental results, and α was calibrated by comparing measured grain size distribution profiles for the bed surface with those predicted by the model. It was found that a good match for roi was obtained using the Garcia (1990) relation, but with a maximum limit set at 2 (Garcia, 1990). The choice of α = 0.8–1.0 appeared to be the most appropriate to match model predictions with experimental data. This is in line with Wright and Parker (2005a), wherein it was observed that the weighting of the interface fractions is more biased toward the active layer fractions in fine sands than in gravel beds because the bed material grain size distribution is generally narrower (Wright and Parker, 2005b).

In general, the value of parameter α is not well known. For pea gravel in laboratory experiments, it has been found that α = 0.2–0.3 (Toro-Escobar et al., 1996; Viparelli et al., 2010). At field scale, a value α = 0.3 has been used to design gravel augmentations on the Trinity River, California (Viparelli et al., 2011), and in Wright and Parker (2005b) α = 1 was used to predict downstream fining in a generic sand-bed stream based on the Mississippi and Atchafalaya Rivers. It has been found that the sediment in the active layer becomes coarser as α decreases, and although it affects the downstream fining pattern, it does not significantly affect the overall equilibrium bed slope (Wright and Parker, 2005b; Viparelli et al., 2010). For the field case in this study, we use α = 0.5 instead of α = 1 for a sand-bed stream because, unlike in Wright and Parker (2005a, 2005b), where active layer thickness was scaled with dune height, the active layer in our model is scaled with the characteristic grain size D90 of the bed surface. A more significant sediment exchange between the flow and the substrate is expected as a result of this much thinner active layer.

The ∼176-km-long Monterey Submarine Canyon extends from Moss Landing, California, into the Pacific Ocean, and is one of the larger submarine canyons along the Pacific margin of North America (Greene et al., 2002). Normark and Carlson (2003, p. 176, table 1) noted the following concerning this canyon: “cross-section at shelf edge is greater than Grand Canyon width and relief.” The Monterey Canyon is a part of the greater Monterey Bay Canyon System, which consists of the Monterey, Soquel, and Carmel Canyons. Soquel Canyon joins Monterey Canyon 18 km seaward of Monterey Canyon's head at a water depth of 970 m. Farther offshore and ∼30 km seaward from its head, Monterey Canyon joins Carmel Canyon at a depth of 2000 m of water (Fig. 6).

Approximately 200,000–250,000 m3 of sediment is delivered to the head of the Monterey Canyon each year via littoral drift (Best and Griggs, 1991). Extensive research on the canyon suggests that it is actively transporting this sediment part way down its axis from the coast toward the Monterey Bay Fan via sediment transport events that occur several times a year (Xu et al., 2004; Paull et al., 2005; Smith et al., 2007). Triggers for these flushing events have been identified as earthquakes and high surf conditions; however, several of these events have been recorded in the absence of either of these conditions (Smith et al., 2007).

One possible mechanism for the formation of sustained turbidity currents under such conditions is breaching (van den Berg et al., 2002; Mastbergen and van den Berg, 2003). That is, as fine sand delivered by littoral drift to the head of the canyon accumulates, it may eventually undergo failure, resulting in sustained breaching. While the breaching shown in Figure 1A was initiated by dredging, van den Berg et al. (2002) noted that breaching events can occur naturally without dramatic forcing, e.g., as depicted in Figure 1B. Here turbidity currents initiated by breaching are studied parametrically.

The model simplifies the sediment pocket emplaced by littoral drift at the head of the canyon so as to take the form of a prismatic sediment reservoir with a vertical breach face of uniform thickness in the transverse direction. Figure 7 shows the long profile of the channel bottom (Paull et al., 2005; Smith et al., 2007). The bed has a 3% bed slope from a depth of 200 to 1200 m, and a 4.5% slope farther down canyon. The points R1, R2, and R3 indicated in the figure are mooring sites where velocity measurements were recorded (Xu et al., 2004). The median grain size D50 of the sediment in the sediment pocket as well as the canyon axis is taken to be 159 μm (Paull et al., 2005). The sediment is assumed to be quartz with a specific gravity of 2.65.

Numerical Experiments

Several scenarios were simulated numerically in order to analyze at field scale the effect of change in breach height, downslope friction coefficient, and grain size distribution on the turbidity current generated and its corresponding depositional signature. All simulations were conducted for a 1 km reach length, with porosities λp = λpr = 0.4, and assuming a horizontal top slope of the breach face (St = 0). The breach wall velocity of 0.0014 m/s used here was estimated based on field observations reported in Smith et al. (2007) and information in Mastbergen and van den Berg (2002).

Although the median grain size remains the same (i.e., D50 = 0.159 mm), two sediment mixtures are considered here: mix 1, consisting of sand only, and mix 2, which has a silt tail. The grain size distributions, which are shown in Figure 8, were selected in part based on information in Combellick and Osborne (1977). Conservative estimates of initial breach heights ranging from 2 to 5 m were considered. Calculations were performed using two friction coefficients, cf = 0.001 and cf = 0.008 (e.g., Normark, 1989), in order to bracket the range of behavior.

The results of four numerical runs are reported here. As documented in Table 1, in runs 1, 2, 3, and 4 the parameters (Zbi, cf, Fsilt) corresponding to initial breach height, friction coefficient, and fraction silt in breach material take the respective values (5 m, 0.001, 0.2), (5 m, 0.008, 0.2), (2 m, 0.008, 0.2), and (2 m, 0.001, 0). Also listed are the temporal and spatial intervals, Δt and Δx respectively. The results are shown in Figures 9–12 in terms of layer averaged downslope velocity U, current thickness H, volumetric concentration of suspended sediment C, geometric mean diameter of the sediment in the active layer of the bed Dg, and geometric mean diameter of the sediment in suspension, Dg susp.

Effect of Bed Friction Coefficient

The only difference between runs 1 and 2 is the bed friction coefficient cf. A comparison of the numerical results shows that as the bed friction coefficient increases from 0.001 in run 1 to 0.008 in run 2, the turbidity current changes from one that decelerates and thickens more rapidly downstream to one that accelerates and thickens less rapidly (Figs. 9 and 10), as seen in the plots of the downslope variation of velocity and current thickness. In particular, run 2 depicts an increasingly erosional down-canyon turbidity current that shows incipient ignition (Pantin, 1979; Parker, 1982). The entrainment rate of bed sediment into the current increases downstream, thus increasing the driving gravity force of the flow. This is in contrast to the decelerating current in run 1 that gradually dies as it deposits its suspended sediment load along the bed. The corresponding current is considerably thinner in run 2; after 1 h of breaching, it shows signs of decreasing in thickness downstream (Fig. 10).

The increasing erosive power of the current in run 2 leads to a coarser active layer, as shown in Figure 11. The same figure shows that (1) in both runs the active layer tend to become coarser in the downslope direction, and (2) in run 1 the active layer tends to become finer in time. An interpretation of this result can be made in terms of the downslope variation of volumetric concentration and geometric mean diameter of suspended sediment represented in Figure 12. In particular, the volumetric concentration of suspended sediment during run 1 decreases in the downslope direction and in time, but the geometric mean diameter of the sediment in suspension remains nearly constant in space and time. This implies that as the current travels in the down-canyon direction, its erosive power decreases and only the finer sediment is entrained in suspension. As time passes, the coarser sediment remains in the upslope region of the canyon and only the finer component is transported downslope, leading to an active layer that becomes finer in time.

In the case of the ignitive turbidity current of run 2, however, both the volumetric concentration of suspended sediment and the geometric mean diameter of the sediment in suspension increase in the downslope direction. This implies that as the current accelerates, its erosive power increases and more coarse sediment is entrained in suspension, with a consequent downslope coarsening of the active layer. The temporal step length that guarantees the stability of the model in the numerical run with the higher friction coefficient that results in an ignitive turbidity current (i.e., run 2) is one order of magnitude smaller than the temporal step length used in the case of a smaller friction coefficient that results in a net-depositional, decelerating turbidity current (i.e., run 1).

Effect of Grain Size Distribution

The only difference between runs 1 and 4 is the silt fraction in the breach material, which is 0.2 in run 1 and 0 in run 4. Both runs are net depositional turbidity currents that decelerate downstream. As seen in Figures 9 and 10, however, the model predicts a slower and thicker current as Fsilt is reduced from 0.2 (run 1) to 0 (run 4). This result, which indicates the tendency for finer sediment to augment the transport of coarser sediment, is in accordance with the results of Salaheldin et al. (2000).

The current of run 4 in particular is observed to decay downstream much more rapidly and die in time much more quickly than run 1. The plots of active layer median grain size for runs 1 and 4 (Fig. 11) show the same pattern of downstream fining in the upstream part of the reach, followed by coarsening as the current deposits and dies along the channel axis. Run 4, however, shows a much shorter region near the breach face showing downstream fining, as the current deposits its sediment more rapidly than run 1.

Effect of Breach Height

Changing the breach height effectively changes the upstream boundary conditions for the flow. As can be seen from equations 9–11, an increase in the breach height corresponds to an increase in upstream sediment feed and hence a stronger turbidity current. The effect of change in breach height is seen by comparing runs 2 and 3, where the breach height is reduced from 5 m to 2 m while keeping all other parameters constant. As shown in Figures 9–11, while the current in run 2 is clearly accelerating, the current with the reduced breach height, i.e., that of run 3, is slowly decelerating and is completely dead within an hour. From the plots of geometric mean diameter of the active layer (Fig. 11), we see that after 24 min, the upslope portion of the canyon bed is significantly finer for the case of the current that completely dies after 1 h (i.e., run 3). In the downslope part of the canyon (i.e., x ≥ 500 m), the difference between the geometric mean diameters between the two runs is negligible.

Active Layer Mean Size Versus Substrate Mean Size

Figure 11 shows that in the case of run 3, the median size in the active layer is coarser than the substrate in the downstream part of the modeled reach. In run 1, however, the sediment in the active layer is everywhere finer than the substrate. The relatively coarse material in the active layer observed for the case of run 3 is likely due to the higher shear velocity (u* = Ucf 0.5) compared to run 1 and thus the higher capacity of the current to entrain sediment into suspension. The result is a tendency for the decelerating current of run 3 to emplace coarser material in preference to finer material.

Model Simulation of a Flushing Event in Monterey Canyon

The model simulation presented here is loosely based on a flushing event in Monterey Canyon between September 2003 and September 2004, as reported by Smith et al. (2005, 2007) and inferred by Xu et al. (2004). Smith et al. (2007) reported the loss of ∼230,000 m3 of sediment from a zone near the canyon head and estimated a corresponding shoreward retreat of the canyon rim of 25 m over the 12 month period. Based on information in Smith et al. (2007), the width of this rim at the channel head was estimated to be 400 m. These numbers yield a rough estimate of the breach height Zb of 20 m, a value adopted for the present simulation.

The three-equation model, while adequately capturing the mechanism of self-acceleration, often shows a tendency toward overacceleration. This results from the failure of the model to include damping of turbulence due to self-stratification (Parker et al., 1986). Other limitations of the three-equation model when applied at field scale were outlined in Kostic and Parker (2006). With this in mind, the four-equation model has been used to model the region down canyon of the breach, in a field-scale simulation loosely based on the flushing event of November 2003 in Monterey Canyon. In the four-equation model, the friction coefficient cf is not constant, but varies according to the relation
where K is the layer-averaged turbulent kinetic energy per unit mass and αk is a parameter that may vary from 0.05 to 0.5 (Parker et al., 1986; Kostic and Parker, 2006). As shown in Parker et al. (1986), the system of equations is relatively insensitive to the value of αk, which in this analysis was assumed equal to 0.1 (e.g., Kostic and Parker, 2006). Details of the four-equation model can be found in Parker et al. (1986). It suffices to mention that the first three equations of the model account for momentum balance, fluid mass balance, and suspended sediment mass balance in forms analogous to that of the three-equation model. The fourth equation describes the layer-averaged balance of turbulent kinetic energy and is expressed in backwater form analogous to equations 1, 2, and 3 as:
where β is the function of Ri and a pseudo-friction coefficient cf* specified as:
In the simulation presented here, the value of cf* has been set equal to 0.004, i.e., a value commonly used at field scale (e.g., Normark, 1989; Kostic and Parker, 2006). In determining an upstream boundary condition for K, we have assumed that the upstream friction coefficient cf0 (equation 19) to be given by
which is set equal to cf*.

The results in terms of layer-averaged depth of the current interface and downslope velocities are plotted in Figure 13. Figure 13A shows the profile of the bed, as well as the predicted profiles of the interface at 2 h and 5 h after initiation of the event. Figure 13B shows predicted profiles of layer-averaged velocity at 2 h and 5 h after the initiation of the event. Also shown in the figure are estimates of layer-averaged downslope velocities based on the velocity measurements taken at mooring sites R2 and R3 by Xu et al. (2004) during the event of November 2003. At mooring R2, layer-averaged velocities were estimated from the measurements at times 2, 3, 4, and 5 h, resulting in values ranging from 0.7 m/s to 0.5 m/s. At this mooring site, the flow showed a clear deceleration over time. At mooring R3, however, the layer-averaged velocity increased from 0.4 m/s at 3 h to 1.1 m/s at 4 h, and then decreased to 0.9 m/s at 5 h.

The numerical run was performed with a downslope step length of 5 m and a time step of 0.1 s. The estimated location of the interface of the turbidity current was always found to be sufficiently deep below the water surface to ensure that the formulation is appropriate to the flow being modeled. Figure 13B illustrates that the computed downslope layer-averaged velocities at the measuring locations are similar to those measured in the field. The predicted lines pass about as well through the measured data as the scatter allows. Figure 13B also shows that the predicted streamwise velocity profiles show rapid drops over the first 2 km of the canyon, and then become roughly constant in space until the break in bed slope evident in Figure 13A is reached. The flow then shows slow spatial acceleration farther downstream.

In the model, the current accelerates slowly in time. The field data, however, show a current at R2 that decelerates in time; at R3 it accelerates between the 3rd and the 4th hour and then decelerates. This discrepancy may in part be due to the assumption of a fully erodible canyon bed in the model. In order to fully capture field conditions, in addition to a limitation induced by the reduction of turbulent kinetic energy due to self-stratification, there should also be a limitation on sediment entrainment associated with insufficient cover of loose sediment over a bed that is otherwise resistant to erosion (e.g., bedrock or stiff consolidated mud).

The analysis presented here supports the hypothesis that breaching may be an important source of turbidity currents in subaqueous canyons or any submarine setting where sufficiently fine sand can accumulate and densify (e.g., some river deltas). In the model presented, the turbidity current is divided into two reaches joined at a migrating boundary: the breach face, which is treated as vertical in the model, and a quasi-horizontal region sloping downdip (down-canyon reach). Three additional sets of equations were developed to describe the flow of the vertical current on the breach face in addition to the three previously derived by Fukushima et al. (1985). The evolving bed was tracked using the Exner equation for conservation of sediment, and the bed grain size fractions were tracked using the formulation detailed in Toro-Escobar et al. (1996). Short numerical experiments have been used to demonstrate that the breaching mechanism can produce the right conditions for a self-accelerating current to develop in a submarine canyon, even where the breach height is limited.

In a field-scale application to a 50 km reach of the Monterey Canyon system, the four-equation formulation of Parker et al. (1986) was used in the down-canyon reach to model turbidity current evolution. Although it is likely that the net retreat of the sand pocket mentioned in Smith et al. (2007) may have been the result of a series of turbidity current events and not just the one reported by Xu et al. (2004), the results predicted by the model based on reasonable assumptions are within the scatter of the data.

C: Layer-averaged volumetric suspended sediment concentration for uniform sediment

Ci, Cvi: Layer-averaged volumetric suspended sediment concentration of sediment in the ith size range for the down-canyon region and breach face, respectively

cf: Down-canyon friction coefficient

cfb: Breach-face friction coefficient

Ct, Cvt: Layer-averaged volumetric concentration of summed over all grain sizes for the down-canyon region and breach face, respectively

Di: Diameter of sediment in the ith size range

Esi: Dimensionless sediment entrainment coefficient of material in the ith size range

ew: Water entrainment coefficient

Fbi: Near bed volumetric upward normal Reynolds flux of sediment in the ith size range

Fi: Grain size fractions of the ith size range in the active layer

Fvi: Grain size fractions of the ith size range in the sediment reservoir

fIi: Grain size fractions of the ith size range in the active layer-substrate interface

g: Acceleration due to gravity, 9.81 m/s2

H, Hv: Turbidity current thickness in the down-canyon region and on breach face, respectively

La: Active layer thickness

qei: Equilibrium volume suspended sediment transport rate per unit width in the ith size range

qi: Volume suspended sediment transport rate per unit width in the ith size range

qtotal: Total volume suspended sediment transport rate per unit width summed over all grain size ranges

R: Submerged specific gravity of sediment, Gs–1

Ri: Bulk Richardson number

roi: Ratio of between the near-bed concentration of suspended sediment in the ith range and the layer-averaged value Ci

S: Down-canyon bed slope

Si: Initial down-canyon slope

St: Breach top slope

t: Time

U, W: Layer- averaged flow velocity components in the x and z directions, respectively

vsi: Fall velocity of sediment in the ith size range

vwall: Breach wall retreat velocity

x: Longitudinal down-canyon coordinate

z: Upward-normal (quasi-vertical) coordinate

Zb: Breach height

α: Mixing parameter (0 ≤ α ≤ 1)

λp: Canyon bed sediment porosity

λpr: Porosity of undisturbed sediment in reservoir and/or breach face

ν: Kinematic viscosity of water

θ: Canyon bed angle

θv: Slope angle of breach face

ρ: Density of water

ρs: Density of sediment

This research represents a contribution of the National Center for Earth-surface Dynamics, a Science and Technology Center funded by the U.S. National Science Foundation (agreement EAR-0120914). We acknowledge valuable discussions with Dick R. Mastbergen, Jan H. van den Berg, and Tetsuji Muto. We also acknowledge the contribution of the reviewers Svetlana Kostic and Alessandro Cantelli in improving this paper, and Davide Motta for discussion on the implementation of the four-equation model.