Introduction

How microscopic interactions affect the macroscopic flow behavior of complex fluids is at the core of soft matter physics. Recently, it has been shown that shear-thickening in dense particulate suspensions corresponds to a frictional transition at the microscopic scale; when the imposed shear stress exceeds the inter-particle short-range repulsive force, the grain contact interaction transits from frictionless to frictional1,2,3,4,5,6,7,8. During this transition, the proliferation of frictional contacts can be so massive that it triggers a remarkable macroscopic rheological response: the rate of shear of the suspension decreases when the imposed shear stress is increased. As a result, highly concentrated shear-thickening suspensions have peculiar S-shape rheological laws9,10, which have been rationalized by a frictional transition model3,11,12. So far, the consequences of the frictional transition and its associated S-shape rheology have been essentially investigated in rheometers, where instabilities, shear bands and spatiotemporal patterns have been documented11,13,14,15. By contrast, very little is known about the behavior of shear-thickening suspensions in real hydrodynamic flow configurations beyond rheometry, in spite of the numerous applications16,17,18.

An archetypical case, which is widely encountered in industrial and geophysical applications, is the incline plane flow configuration. As previously reported19 and illustrated in Fig. 1a (see also Supplementary movie 2), when a thin layer of shear-thickening suspension flows down an inclined plane, surface waves of wavelengths much larger than the thickness can develop spontaneously and grow as they propagate downstream. This longwave free-surface instability may seem reminiscent of the Kapitza instability observed when a thin liquid film flows down a slope20,21, or more generally of the so called “roll waves” instability observed from turbulent flows in open channels22,23,24, to avalanches of complex fluids like mud25,26 or granular media27. These latter two instabilities rely on the same primary mechanism: the amplification of kinematic surface waves at high velocity owing to inertial effects28. For a Newtonian liquid in the laminar regime, the destabilization occurs only when the Reynolds number of the flow, Re = ρu0h0/η, where ρ is the fluid density, u0 its mean velocity, h0 the flow thickness and η the fluid viscosity, exceeds the Kapitza threshold, \(R{e}_{K}=5/(6\tan \theta )\), which is typically much larger than 1 for a small tilting angle θ of the incline29,30,31. By contrast, the growth of surface waves observed in Fig. 1a for a dense shear-thickening suspension occurs at a Reynolds number of only  ≈1, i.e., far below the Kapitza threshold ReK ≈ 5 predicted for θ = 10 (see also ref. 19). This suggests that a different instability mechanism is at play for dense shear-thickening suspensions, yet its origin remains an open question.

Fig. 1: Experimental characterization of the instability onset.
figure 1

a Non-inertial surface waves emerging spontaneously when a concentrated suspension of cornstarch particles flows down an incline (volume fraction ϕ = 0.45, inclination angle θ = 10 and normalized flow Reynold number Re/ReK ≈ 0.2). b Sketch of the experimental setup. We use the progressive drainage of the reservoir to quasi-steadily vary the flow rate. The instability onset is determined by measuring the wave amplitudes both at the top and at the bottom of the incline with two laser sheets and cameras. For ϕ ≤ 0.4, an oscillation of the gate is added to impose a controlled perturbation. c Spatiotemporal plots of the laser sheet transverse-position versus time, indicating the vertical oscillations (blue and red arrows) of the free surface, at the top and at the bottom of the incline (ϕ = 0.33, θ = 2, Re ≈ 37). d Reynolds number of the flow, Re, and e amplitude of the perturbation at the top, Δh1, and at the bottom, Δh2, during the drainage of the suspension reservoir (ϕ = 0.36, θ = 3). The instability onset (Δh1 = Δh2) is given by Rec ≈ 28 (black-dashed-line).

Here, we investigate the origin of this instability by studying the flow of a shear-thickening suspension down an inclined plane over a wide range of volume fractions and flow rates. We confirm that this instability is not inertial and fundamentally different than the classical Kapitza or Roll waves instabilities. We provide experimental evidence together with a theoretical explanation, which show that this destabilization arises from the coupling between the flow free surface and the non-monotonic (S-shape) rheological laws of shear-thickening suspensions.

Results and discussion

Evidence of an instability distinct from the classical Kapitza or roll waves instabilities

We perform experiments with shear-thickening aqueous suspensions of commercial native cornstarch (Maisita®, http://www.agrana.com). We vary the particle volume fraction over a wide range (0.30 < ϕ < 0.48) and characterize the onset of stability (the value of ϕ refers to the dry volume of cornstarch computed from its dry weight and density, 1550 kg m−3). We use a 1 m long and 10 cm wide inclined plane covered with a diamond lapping film (663-3M with roughnesses of  ~45 μm) to insure rough boundary conditions. The suspension is released from a reservoir at the top of the plane through a gate with an adjustable aperture (Fig. 1b). A scale, placed at the bottom end of the incline (not shown in the schematic), provides the instantaneous flow rate q of the suspension. To probe the stability of this free-surface flow for moderate volume fractions ϕ ≤ 0.4, the gate is mounted on a translating stage imposing a small sinusoidal modulation of its aperture (3 Hz,  ±100 μm). At large volume fractions however (ϕ 0.4), no forcing is required because the flow is so unstable that it is dominated by noise amplification of its most unstable mode. Two low-incident laser sheets and two cameras are used to measure the mean film thickness h0~2−10 mm, and the crest-to-crest amplitude of the waves upstream and downstream of the incline (Fig. 1c). The calibration of the laser sheet projections on the film surface yields a local measurement of h0 with a precision of 10 μm.

We use a protocol designed to characterize the instability onset with a single experiment for each volume fraction ϕ and inclination θ. The flow rate is varied quasi-steadily, either using the progressive drainage of the reservoir (decreasing flow rate) or by slowly increasing the gate aperture (increasing flow rate). At each instant, an effective Reynolds number of the flow is computed from the instantaneous flow rate q and from the mean film thickness h0 using the relation \(Re=3{q}^{2}/(g{h}_{0}^{3}\sin \theta )\), where g is the gravitational acceleration. Note that the flow rate is varied sufficiently slowly that at each instant, the flow rate is constant along the incline. With such a definition based on mean quantities only, the Reynolds number can be applied to any rheology and is directly related to the Froude number \(F={u}_{0}^{2}/(g{h}_{0}\cos \theta )\) commonly used to described roll waves24 by \(Re={3F}/{\tan} \theta\), where u0 = q/h0 is the mean flow velocity. Re also reduces to the standard expression Re = ρu0h0/η for a Newtonian fluid of viscosity η and density ρ in the laminar regime (Nusselt velocity profile). Figure 1d shows the concomitant evolution of Re and of the wave amplitude upstream (Δh1) and downstream (Δh2) the incline during the drainage of the reservoir, starting from an unstable situation where Δh2 > Δh1. The stability onset is precisely reached when Δh1 = Δh2, which provides us with the critical Reynolds number Rec (dashed line in Fig. 1c), the critical flow thickness hc, and the critical basal shear stress, \({\tau }_{c}=\rho g{h}_{c}\,\sin \theta\). We have verified that the same instability onset is obtained by carrying successive steady-state measurements at various constant flow rates. A new freshly prepared suspension of cornstarch is used for each measurement. Experiments are repeated four times for each volume fraction.

  Figure 2a shows the critical Reynolds number Rec, normalized by the Kapitza threshold for a Newtonian fluid, as a function of ϕ. For the lowest volume fraction investigated, ϕ = 0.30, the stability threshold is close to ReK (it is typically 50% above, owing to the finite forcing frequency and finite width of the plane32,33). As the volume fraction is increased, Rec becomes increasingly larger relative to ReK, reaching  ≈5ReK at ϕ = 0.41. This behavior is actually expected from Kapitza’s inertial mechanism for a medium that is continuously shear-thickening34, like our cornstarch suspension over that range of volume fractions (0.35 ϕ 0.41). More strikingly, for ϕ 0.41 the relative critical Reynolds number drops drastically, down to two orders of magnitude below the Kapitza threshold at the largest volume fraction investigated (ϕ = 0.48, for which Re ≈ 0.15 and F ≈ 10−2). Clearly, in this domain, the flow destabilization can no longer be explained within the Kapitza framework since inertial effects are negligible. Note that inertia is also negligible at the particle scale, since the Stokes number, \(St \sim {(d/{h}_{0})}^{2}Re\), where d ~ 10 μm is the particle size, is  ~105 times smaller than Re. As shown in Fig. 2b, this qualitative change in the onset of instability around ϕ = 0.41 is also observed in the evolution of the critical shear stress τc with ϕ. Similarly, the surface wave speed at the instability threshold changes significantly and abruptly. Its value cc, normalized by the mean flow velocity u0 drops from 3, which is Kapitza’s prediction for a Newtonian layer in the long wavelength limit, to ~2 for volume fractions exceeding  ≈0.41 (see Fig. 2c). These results confirm that above ϕ ≈ 0.41, a longwave free-surface instability, which is fundamentally distinct from the Kapitza instability, emerges. In the following, we call this instability, which to our knowledge as no equivalent in classical fluids, “Oobleck waves”.

Fig. 2: Onset of instability.
figure 2

a Critical Reynolds number of the instability normalized by the Kapitza threshold Rec/Rek versus volume fraction ϕ. Inset: Rec/Rek versus inclination angle θ for ϕ = 0.45. b Critical shear stress τc versus ϕ. c Normalized critical wave speed cc/u0 versus versus ϕ. Inset: cc/u0 versus speed of the kinematic waves ckin/u0, the black line shows that the waves propagate at the speed of the surface kinematic waves. Dashed-blue-line: Kapitza prediction (inertia+Newtonian fluid). Solid-red-line: prediction of the linear stability analysis (\(A={\rm{d}}\tilde{\dot{\gamma }}/{\rm{d}}{\tilde{\tau }}_{b}{| }_{{\tilde{\tau }}_{b} = 1}=0\)) describing the coupling between the flow free-surface and the suspension shear-thickening rheology. Different symbols indicate different inclination angles θ: 2, 3, 6, 9, 10, 22. The error bars indicate the standard deviation between experimental measurements. Different background colors highlight which instability emerges: Kapitza (gray) or Oobleck waves (red).

Oobleck waves arise from the S-shape rheology and kinematic wave propagation

To understand the origin of Oobleck waves, we characterize the rheology of the cornstarch suspension in a cylindrical-Couette rheometer (Fig. 3). We find that the volume fraction at which Oobleck waves appear (ϕ ≈ 0.41) corresponds precisely to ϕDST, the volume fraction at which the shear-thickening transition becomes discontinuous. Indeed, we analyze our rheological data along Wyart & Cates’s model, which assumes that the effective viscosity of the suspension, \(\eta (\phi ,\tau )={\eta }_{s}{({\phi }_{J}(\tau )-\phi )}^{-2}\), diverges at a critical volume fraction, ϕJ, that depends on the applied shear stress, τ, according to \({\phi }_{J}(\tau )={\phi }_{0}(1-{e}^{-{\tau }^{* }/\tau })+{\phi }_{1}{e}^{-{\tau }^{* }/\tau }\), where (ηs, τ*, ϕ0, ϕ1) are material constant. Here, ηs is a prefactor proportional to the solvent viscosity, τ* is the short-range repulsive stress scale above which the frictional transition occurs, which may be tuned by changing the particle roughness or surface chemistry, ϕ0 (resp. ϕ1) is the jamming volume fraction at which the suspension viscosity diverge at low (resp. large) stress, with ϕ1 being dependent on the inter-particle friction coefficient3. By fitting our measurements with this model, we find that the rheological curve (shear stress τ versus shear rate \(\dot{\gamma }\)) becomes S-shaped when ϕ ≥ 0.41 ± 0.005 = ϕDST (see Fig. 3 and Methods for the fitting procedure). This suggests that negatively sloped portion in the rheological curve (\({\rm{d}}\dot{\gamma }/{\rm{d}}\tau \, < \, 0\)) is a key ingredient of the instability. S-shaped flow curves, and more generally rheograms with a negatively sloped region, are known to produce unstable flow conditions35,36. Previous studies on shear-thickening suspensions have based their analysis on this feature to explain for instance the emergence of random fluctuations37 reported initially by Boersma et al.38, and the oscillations observed when an object moves in a shear-thickening fluid39,40 or in rheometric configurations15,41. However, all these models require inertia to predict an instability. By contrast, here, the instability seems to be of a fundamentally different nature. First, it can occur at very low Reynolds and Froude numbers, for which inertial effects are negligible. Second, at the instability onset, the unstable mode propagates at the speed of the surface kinematic waves defined by ckin ≡ dq/dh028 (Fig. 2c inset). This indicates that the coupling between the flow and the free-surface deformation, which was not considered in previous studies, is essential to explain the emergence of Oobleck waves.

Fig. 3: Rheograms of the aqueous cornstarch suspension.
figure 3

Shear-stress τ versus shear rate \(\dot{\gamma }\) for various volume fractions ϕ. Solid lines: fit by Wyart & Cates rheological laws setting the jamming volume fraction for frictionless and frictional particles to ϕ0 = 0.52 ± 0.005 and ϕ1 = 0.43 ± 0.005, respectively, the short-range repulsive stress scale above which the frictional transition occurs to τ* = 12 ± 2 Pa and the prefactor to ηs = 0.91 ± 0.01 mPa s. The rheograms are negatively sloped (\({\rm{d}}\dot{\gamma }/{\rm{d}}\tau <0\)) in the region highlighted in blue.

Oobleck waves instability mechanism

We now show how the negative slope in the rheology coupled with a gravity-driven free-surface flow can yield to an instability without invoking inertial effects. In the zero-Reynolds number limit, the force balance on a slice of suspension, as depicted in Fig. 4a, imposes that the basal stress, τb, is equal to the sum of the projected weight of the slice, \(\rho gh\sin \theta\), where h(xt) is the local flow thickness, and of the longitudinal pressure gradient induced by the free-surface deflection. For wavelengths much larger than the flow thickness and the capillary length, the pressure profile perpendicular to the plane can be assumed to be hydrostatic \(P(x,z,t)=\rho g\cos \theta (h(x,t)-z)\), where x is the flow direction and z the height within the flowing layer in the perpendicular direction21. The depth-averaged force balance is then given by

$${\tau }_{b}=\rho gh\sin \theta -\rho gh\cos \theta \frac{\partial h}{\partial x}.$$
(1)

Let us now consider a perturbation of a base flow of constant thickness, as illustrated in Fig. 4b. A local increase of the flow thickness implies that ∂h/∂x becomes positive upstream of the perturbation and negative downstream. To satisfy the force balance (1), the basal shear stress τb upstream must therefore decrease, whereas it must increase downstream. However, owing to the S-shape rheology of the suspension, when \({\rm{d}}\dot{\gamma }/{\rm{d}}\tau \, < \, 0\), a decrease (resp. increase) in τb implies a local increase (resp. decrease) of the flow rate \(\dot{\gamma }\). Therefore, the shear rate increases upstream and decreases downstream, inducing a net inward mass flux underneath the bump and the amplification of the initial perturbation (red arrows in Fig. 4b).

Fig. 4: Instability mechanism.
figure 4

a Depth-averaged forces acting along the flow x direction on a slice of suspension of width dx (shaded in blue): basal force τbdx, projected weight of the slice \(\rho ghdx\sin \theta\) and hydrostatic pressure \(\rho g\cos \theta {h}^{2}(x)/2\), where τb is the basal shear stress, g is gravity, h(x) is the flow thickness, θ is the plane inclination angle. b Positive feedback for a shear-thickening suspension with a S-shaped rheological curve: a local increase of the flow thickness h implies that ∂h/∂x is positive (resp. negative) upstream (resp. downstream) of the perturbation. Force balance (see Eq. (1)) then implies that the basal shear stress τb upstream (resp. downstream) must decrease (resp. increase). When the suspension rheogram is negatively sloped (\({\rm{d}}\dot{\gamma }/{\rm{d}}\tau <0\)), this yields a local increase (resp. decrease) of the flow rate \(\dot{\gamma }\). The combination of these two feedback cycles (gray arrows) induces a net inward mass flux towards the bump (red arrows) that amplifies the initial perturbation (red vertical arrow).

Quantitative depth-averaged model without inertia

To go beyond this qualitative picture, we perform a linear stability analysis of a steady uniform flow of thickness h0 and depth-averaged velocity u0 using the Saint-Venant approximations (long wavelength limit27,42) and neglecting inertia (lubrication approximation21). The equations are written using the dimensionless variables \(\tilde{h}=h/{h}_{0}\), \(\tilde{x}=x/{h}_{0}\), \(\tilde{u}=u/{u}_{0}\), \(\tilde{t}=t{u}_{0}/{h}_{0},\) \({\tilde{\tau }}_{b}={\tau }_{b}/\rho g{h}_{0}\sin \theta\) and linearized by writing \(\tilde{h}=1+{h}_{1}\), \(\tilde{u}=1+{u}_{1}\) and \({\tilde{\tau }}_{b}=1+{\tau }_{1}\) with (h1u1τ1 1). Under these conditions, the mass conservation, \({\partial}_{\tilde{t}} h+{\partial}_{\tilde{x}}(hu)=0\), and the force balance (1) become \({\partial}_{\tilde{t}} {h_{1}}+{\partial}_{\tilde{x}}{h_{1}}+{\partial}_{\tilde{x}}{u_{1}}=0\) and \({\tau }_{1}={h}_{1}-\tan {\theta }^{-1}{\partial }_{\tilde{x}}{h}_{1}\), respectively. The linearization of the normalized shear-rate \(\tilde{\dot{\gamma }}(\tilde{{\tau }_{b}})\equiv \tilde{u}/\tilde{h}\), gives Aτ1 = u1 − h1, where \(A={\rm{d}}\tilde{\dot{\gamma }}/{\rm{d}}{\tilde{\tau }}_{b}{| }_{{\tilde{\tau }}_{b} = 1}\) is the slope of the rheological curve for the base state basal stress \({\tau }_{b}=\rho g{h}_{0}\sin \theta\), obtained from the integration of the flow velocity profile (see Methods). Taking the spatial derivative of the force balance and substituting τ1 and \({\partial}_{\tilde{x}}{u_{1}}\) using the rheology and mass balance, respectively, lead to a single partial differential equation for the free-surface perturbation h1:

$$\frac{\partial {h}_{1}}{\partial \tilde{t}}+\tilde{c}\frac{\partial {h}_{1}}{\partial \tilde{x}}=\frac{A}{\tan \theta }\frac{{\partial }^{2}{h}_{1}}{\partial {\tilde{x}}^{2}},$$
(2)

where \(\tilde{c}={c}_{{\rm{kin}}}/{u}_{0}=2+A\) is the dimensionless speed of the kinematic waves28. Interestingly, the perturbation amplitude is found to follow a diffusion equation in the reference frame of the kinematic waves, with an effective “diffusion coefficient” \(A/\tan \theta\). When A < 0, i.e., when the slope of the rheological law \({u}_{0}/{h}_{0}=\dot{\gamma }({\tau }_{b})\) is negative, anti-diffusion occurs, which leads to an amplification of all perturbations, whereas for A > 0 the flow is stable. The onset of instability is thus given by A = 0. This criteria can be expressed in terms of a critical Reynolds number Rec and a critical basal shear stress τc using the Wyart & Cates rheological laws (see Methods). Figures 2a, b show that, for ϕ ≥ ϕDST, these predictions (red-solid lines) capture very well the value of the critical Reynolds Rec and its dramatic drop over two decades when increasing ϕ, as well as the order of magnitude and the drop of the critical shear stress τc with ϕ. The decrease of the instability onset with increasing ϕ is a direct consequence of Wyart and Cates’ rheological laws3, where the DST onset stress also decreases with ϕ. Physically, it comes from the fact that, when approaching the maximal possible packing fraction ϕ0, less and less frictional contacts are required to reach the DST region. The model also predicts a weak dependence of Rec on the plane inclination angle, as observed experimentally (Fig. 2a inset). This overall agreement is all the more conclusive that it is rooted on physically based constitutive laws, of which the rheological parameters are measured independently, without further fitting. Another strong prediction of the model is that, at the instability onset (A = 0), the speed of the unstable mode is equal to the speed of the kinematic waves cc/u0 = 2. This prediction is fully consistent with the drop and value of the normalized wave speed observed experimentally for ϕ > ϕDST (see Fig. 2c). These results conclusively show that surface waves can emerge from the coupling between a negatively sloped rheology and the flow free surface, without the need for inertial effects.

The Oobleck waves instability mechanism highlighted in this study is not limited to shear-thickening suspensions and could be extended to any other complex fluids having a rheology with a negatively sloped region (e.g., granular materials and geomaterials exhibiting velocity-weakening rheology43,44, concentrated polymers or surfactant solutions35,36, liquid crystals45, active self-propelled suspensions46). More generally, our analysis shows that gravity forces, which are usually stabilizing for gravity-driven free-surface flows, can become destabilizing in the presence of a non-monotonic rheology. Our result could thus be extended to other stabilizing forces such as capillary forces arising from the free-surface deformation. We thus anticipate that other interesting instabilities may be explained directly, or in the light of our study.

Finally, in a broader context, our study reveals that kinematic waves can be unstable in an overdamped medium, where inertia is negligible. These waves, which are observed in a wide range of situations (e.g., traffic and pedestrian flows47, sediment transport48, fluidized bed49, surge, and floods23), result from mass conservation and a general relationship between a local flow rate and a local concentration50,51 (e.g., number of vehicles or pedestrians on a road, solid packing fraction in a suspension, depth of the flow). However, in all these systems, their emergence from a uniform state requires an inertial lag between the flow rate and the concentration28. In our system, the situation is very different as these waves are unstable not from inertia, but from the intrinsic constitutive flow rule of the material. Whether this description of kinematic waves could be extended to more complex systems, such as crowds47 or active self-propelled particles46, are interesting topics to address in future studies.

Methods

Rheological data fitting procedure

Rheograms of the aqueous suspension of cornstarch are measured for various volume fractions in a narrow-gap cylindrical-Couette cell (Fig. 5a) using a rheometer (Anton Paar MCR 501). The height of the shear-cell (40 mm) is sufficiently large to neglect sedimentation effects of the particles during the measurement. Similarly to the procedure followed by Guy et al.4, the viscosity below (ττ*) and above (ττ*) the shear-thickening transition are extracted and plotted versus ϕ (Fig. 5b). The low viscosity branch (frictionless branch) is first fitted with \(\eta (\phi )={\eta }_{s}{({\phi }_{0}-\phi )}^{-2}\), with ηs and ϕ0 as fitting parameters. This yields ηs = 0.91 ± 0.01 mPa.s and ϕ0 = 0.52 ± 0.005. The large viscosity branch (frictional branch) is then fitted with \(\eta (\phi )={\eta }_{s}{({\phi }_{1}-\phi )}^{-2}\), using the previous estimation of ηs, and letting ϕ1 as the only fitting parameter. Note that the rheograms are obtained using both very rough (square symbols) and rough walls (circle symbols), by covering the cell-walls with sand papers of different grades (roughnesses of  ≈80 μm and  ≈15 μm, respectively). The two measurements overlap, except in the frictional branch at high volume fraction (shaded symbols in Fig. 5b). For instance, the data from ϕ = 0.4 and 0.41 are included in the fitting procedure, while 0.42 and 0.43 are not, because the first two points overlap, independently of the roughness of the boundaries, whereas for 0.42 and 0.43 systematic deviations are observed indicating slippage or other artefacts. All data points which are interpreted as biased measurements (transparent symbols) are discarded from the fitting procedure; this yields ϕ1 = 0.43 ± 0.005. Once the values of ηs, ϕ0 and ϕ1 are set, we determine the value of τ* by fitting the full rheograms \(\tau (\dot{\gamma })\) with Wyart & Cates laws: \(\tau ={\eta }_{s}{({\phi }_{J}(\tau )-\phi )}^{-2}\dot{\gamma }\), with \({\phi }_{J}(\tau )={\phi }_{0}(1-{e}^{-{\tau }^{* }/\tau })+{\phi }_{1}{e}^{-{\tau }^{* }/\tau }\). The best fit, shown in Fig. 5c, is obtained for τ* = 12 ± 2 Pa. The value of τ* represents the critical shear stress required to overcome the inter-particle repulsive force and activate frictional contacts between particles. For an inter-particulate force f and a particle size d, the critical shear stress is expected to be of order f/d2. The value we obtain (≈12 Pa) is consistent with the values already reported in the literature for cornstarch in water.

Fig. 5: Rheological data fitting procedure.
figure 5

a Cylindrical-Couette geometry used to characterize the rheology of the aqueous suspension of cornstarch (yellow). The gray area represents the cylinder rotating in the direction indicated by the black arrow. b Viscosity (Pa.s) below (ττ*) and above (ττ*) the shear-thickening transition versus ϕ. Blue-solid-line: frictionless branch, red-solid-line: frictional branch, black-dashed-lines: jamming volume fractions for the frictionless and frictional branches defining ϕ0 and ϕ1, respectively. Square and circle symbols correspond to measurements performed with different wall roughness. The data points made transparent are not used in the fitting procedure as for these points, the measurement depend on the wall roughness. c Shear stress τ versus shear rate \(\dot{\gamma }\) measured at various volume fraction ϕ. Solid lines: Wyart & Cates rheological laws. The blue shaded area highlights the region where the rheograms are negatively sloped (\({\rm{d}}\dot{\gamma }/{\rm{d}}\tau <0\)).

Computation of τ c and R e c

To compute the critical shear stress τc and the critical Reynold number Rec from the instability criteria resulting from the linear stability analysis \(A\equiv {\rm{d}}\tilde{\dot{\gamma }}/{\rm{d}}\tilde{{\tau }_{b}}{| }_{\tilde{{\tau }_{b}} = 1}=0\), we need to relate the shear rate \(\dot{\gamma }\equiv {u}_{0}/{h}_{0}\), defined as the ratio of the depth-averaged flow velocity to the flow thickness, to the basal shear stress τb and the basal suspension viscosity η(τb).

For a steady uniform flow down an inclined plane of slope θ, the momentum equation applied to a surface layer of thickness h0 − z gives

$$\tau (z)=\rho g\sin \theta ({h}_{0}-z)=\eta (z)\frac{{\rm{d}}\hat{u}(z)}{{\rm{d}}z},$$
(3)

where the second equality uses the definition of viscosity, \(\eta =\tau /({\rm{d}}\hat{u}/{\rm{d}}z)\), and \(\hat{u}(z)\) is the local velocity parallel to x. From the proportionality between τ and h0 − z, the local velocity can be expressed as

$$\hat{u}(\tau )=\frac{1}{\rho g\sin \theta }\mathop{\int}\nolimits_{\tau }^{{\tau }_{b}}\frac{\tau ^{\prime\prime} }{\eta (\tau ^{\prime\prime} )}{\rm{d}}\tau ^{\prime\prime} .$$
(4)

Using the definition of the depth-averaged flow velocity, \({u}_{0}=\mathop{\int}\nolimits_{0}^{{\tau }_{b}}\hat{u}(\tau ^{\prime} )\ {\rm{d}}\tau ^{\prime} /{\tau }_{b}\), we obtain the expression of the depth-averaged shear rate

$$\dot{\gamma }=\frac{{\tau }_{b}}{3\eta ({\tau }_{b})}{\mathcal{G}}({\tau }_{b}),$$
(5)

where

$${\mathcal{G}}(\tau )=\frac{3\eta (\tau )}{{\tau }^{3}}\mathop{\int}\nolimits_{0}^{\tau }\mathop{\int}\nolimits_{\tau ^{\prime} }^{\tau }\frac{\tau ^{\prime\prime} }{\eta (\tau ^{\prime\prime} )}{\rm{d}}\tau ^{\prime\prime} {\rm{d}}\tau ^{\prime} ,$$
(6)

embeds the shear-thickening of the suspension. By definition, \({\mathcal{G}}=1\) for a Newtonian fluid.

From (5) and (6) we obtain

$$A\equiv \frac{{\rm{d}}\tilde{\dot{\gamma }}}{{\rm{d}}\tilde{{\tau }_{b}}}{| }_{\tilde{{\tau }_{b}} = 1}=\frac{{\tau }_{b}}{\dot{\gamma }}\frac{{\rm{d}}\dot{\gamma }}{{\rm{d}}{\tau }_{b}}{| }_{{\tau }_{b} = \rho g\sin \theta {h}_{0}}=\frac{3}{{\mathcal{G}}({\tau }_{b})}-2,$$
(7)

where we have used \(\frac{{\rm{d}}}{{\rm{d}}\tau }(\mathop{\int}\nolimits_{0}^{\tau }\mathop{\int}\nolimits_{{\tau }^{\prime}}^{\tau }\frac{{\tau }^{^{\prime\prime} }}{\eta ({\tau }^{^{\prime\prime} })}{\rm{d}}{\tau }^{^{\prime\prime} }{\rm{d}}{\tau }^{\prime})={\tau }^{2}/\eta (\tau )\). Note that in (7), the basal stress τb, the shear rate \(\dot{\gamma }\) and the derivative are taken at the base state, i.e., for \(\dot{\gamma }={u}_{0}/{h}_{0}\) and \({\tau }_{b}=\rho g\sin \theta {h}_{0}\).

Finally, the critical shear stress τc, at which the flow destabilizes (black-solid-line plotted in Fig. 2b of the main text), is obtained numerically by finding the value of τb for which A(τc) = 0, i.e., \({\mathcal{G}}({\tau }_{b}={\tau }_{c})=3/2\). From the value of τc, we obtain the critical Reynolds number (black-solid-line plotted in Fig. 2(a) of the main text)

$$R{e}_{c}\equiv \frac{3{u}_{0}^{2}}{g{h}_{0}\sin \theta }=\frac{3{{\tau }_{c}}^{3}{\left({\phi }_{J}({\tau }_{c})-\phi \right)}^{4}}{9{\eta }_{s}^{2}\rho {g}^{2}{\sin }^{2}\theta }{[{\mathcal{G}}({\tau }_{c})]}^{2}=\frac{3{{\tau }_{c}}^{3}{\left({\phi }_{J}({\tau }_{c})-\phi \right)}^{4}}{4{\eta }_{s}^{2}\rho {g}^{2}{\sin }^{2}\theta }.$$
(8)