E3S Web of Conferences 415, 01008 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341501008

Mass exchange between geophysical flows and beds: idealised
computational modelling using a Herschel-Bulkley rheology
Saoirse Robin Goodwin1,2*, Guillaume Piton1 and Guillaume Chambon1
1Univ. Grenoble-Alpes, INRAE, UR ETNA, 38402 Saint-Martin-d’Hères, France
2Karlsruher Institut für Technologie, 76131 Karlsruhe, Germany

Abstract. A key mechanism by which geophysical flows evolve is mass exchange with the underlying bed,
either by entraining material from the bed, or by depositing material. Although it is known that someconsequences of these mass exchange processes include changes in the volume, momentum and local
rheology of the flow, the circumstances under which specific changes occur are not well-established. Given
the enormous number of competing mechanisms present in geophysical flows, it is not surprising that thestate of the art for modelling entrainment is essentially still empirical. In this study, we implement aHerschel-Bulkley (non-Newtonian) rheology into an existing open-source Smoothed Particle
Hydrodynamics solver (DualSPHysics). This rheology can reasonably represent clay-rich flows, typical of
those observed in the French Prealps. We hence undertake a highly-idealised, quantitative investigation ofentrainment mechanisms for flows overriding non-fixed beds. For the beds, we vary the yield stress and thedepth. Preliminary results reveal a rich variety of behaviours that can be obtained for different bed properties,
including both acceleration and deceleration of the flow material. These mechanisms are reminiscent (butnot identical) of observations from other studies where geo-materials were used.

1 Introduction
There is an increasing risk due to geophysical flows
across the globe, as a direct consequence of climatechange, as well as a continued propensity for humandevelopments in mountainous areas worldwide [1].
Ultimately, this means that engineers and scientists
increasingly need to identify which areas are likely mostat risk, and to develop countermeasures accordingly.
One key tool for doing so is explicit numericalmodelling of potential flows, using approaches based on
fluid mechanics. In general, numerical solvers
implement expressions describing conservation of bothmass and momentum, either using the full NavierStokes equations, or the Saint-Venant (shallowwater/depth-averaged) equations. Various constitutive
models can then be implemented, for instance bymodifying the stresses contributing to the resisting terms(e.g. the Voellmy model), or by including earth pressure
coefficients for the pressure gradient terms [2, 3]. In
principle, these models can predict both the trajectory offlows, as well as mass exchange mechanisms(entrainment and deposition) that are well-known tocause enormous changes in the flow volume [4-9].
Given that models based on the Saint-Venant
equations run fastest – since they are depth-averaged,
they are solving for one fewer dimension than the fullNavier-Stokes equations – they are of high interest topractising engineers. Saint-Venant solvers can consider
mass exchange between the flow and the underlying bed
through the expression for the conservation of mass [2,3, 10, 11]. Indeed, for a depth-averaged incompressible
flow, we can write an expression for the conservation ofmass that includes a term for mass exchange b. The signof b distinguishes entrainment and deposition.
The crux of the matter rests on how we should becomputing the term b. For earlier pioneering work using
Saint-Venant solvers, b was set to have a constant value
for the entire duration of flows [2, 10]. Whilst this wasa fine starting point, solvers were essentially blind to thelocal underlying topography. More recently, the term b
has been adapted to account for empirical entrainment
relationships (see [12-14] for examples of such
relationships, and [15] for an implementation).
Although this represents a logical extension of earlierworks, empirical relationships tend to be site-specific,
and in any case rely on the availability of statisticallysignificant and well-characterised datasets. It is worth
noting that in principle, it should be possible to derive
an expression for b from the rheologies of the flow andthe bed, thus circumventing the need for ad-hocrelations. Nonetheless, such expressions for b would
requireevaluation
againstexperimentalor
computational data for mass-exchange processes (thelatter being explored in this Extended Abstract).
It is also worth noting that up until relatively
recently, depth-averaged solvers have not accounted forthe effects of mass exchange on expressions forconservation of momentum. This, as well as the state of
the art for mass-exchange problems continuing to relyon empirical relationships, belies a fundamental lack ofunderstanding of mass-exchange mechanisms.

* Corresponding author: srgoodwin@protonmail.com© The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0(https://creativecommons.org/licenses/by/4.0/).

E3S Web of Conferences 415, 01008 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341501008

Authors such as Hu et al. [16], Song & Choi [17]and Bates & Ancey [18] have done sterling work to
study mechanisms of entrainment using small-scale
flume apparatus. In particular, [18] used Carbopol,
which can be modelled using a Herschel-Bulkley (nonNewtonian) rheology, to model both flows and nonfixed beds. Non-Newtonian rheologies constitute a
reasonable simplification of clay-rich flows andadvantageously comprise relatively few parameters.
Bates & Ancey studied a range of bed depths andlengths, and hence evaluated newly-proposed equations
describing the evolution of both the flow and the bed.
In this study, we use results from Bates & Ancey[18] to evaluate a Smoothed Particle Hydrodynamics(SPH) model, into which we implemented a regularised
Herschel-Bulkley rheology. The SPH model that weused solves the full 2D Navier-Stokes equations. We
then perform a broad parametric study on the rheological parameters of the bed (including the yield stress),
as well as the depth of the bed, thus extending the workof Bates & Ancey. The objectives of the parametricstudy are (i) to delineate regimes for different behaviour,and hence (ii) to provide a basis for developing/evaluating physically-based expressions for b.

𝜕𝜌= −𝜌∇ ∙ 𝑈(2)𝜕𝑡𝜕𝑈 1= (−∇𝑃 + ∇ ∙ τ) + 𝑔(3)𝜕𝑡 𝜌
where ρ is density, P is pressure, τ is the extra-stresstensor and g is acceleration due to gravity.
To implement the Herschel-Bulkley rheology, weused the well-established technique of regularisation[22], wherein we compute an equivalent viscosity µeffwhich can be capped arbitrarily at a threshold µmax. (Weset µmax = 10 kPa, following a sensitivity study.) The
expression for the rheology can thus be re-cast as:
𝜇eff = 𝜇max
if 𝜇eff ≥ 𝜇max
𝜏0}(4)𝑛−1𝜇eff =+ 𝑘|𝛾̇ |if 𝜇eff < 𝜇max|𝛾̇ |The advantages of this approach are twofold: itmeans that µeff is determinate everywhere, and theminimum timestep required to solve the simulation canbe limited (since the timestep is inversely related toviscosity terms). However, a limitation of the approachis that true unyielded zones cannot be resolved. In otherwords, even material where µeff = µmax can continuously
yield. The rate of this yielding is determined by µmax,
and should ideally be as low as possible, to approximate“true” unyielded material.

2 Numerical model
2.3 Numerical setup
2.1 Rheological model

The numerical setup for our study is modelled onthat of [18]. The dam-break technique was used,
wherein material is placed initially in the storage area at
the top of the flume. A non-fixed bed was included at acertain distance downstream from the gate. For all of thesimulations presented in this Extended Abstract, thelength of the bed was set at a constant 300 mm.
The simulations were split up into two distinct
simulations: (i) dam-break; and (ii) mass exchange
between the flow and the bed. The single simulation for
dam-break was recycled for all of the mass-exchange
simulations. This approach had two advantages: itgreatly reduced the amount of time required for eachsimulation, and it also reduced the chances of numericalinstabilities manifesting.
Fig. 1a shows the starting condition for the dambreak simulation. We stopped this simulation before theflow material reached the non-fixed bed. Each of themass-exchange simulations started by regenerating theflow material according to the flow bounds, as well as
its bulk velocity (see Fig. 1b).
The density of the flow material was set at 1000kg/m3, but this value could vary slightly given that wewere using a slightly-compressible formulation. The
channel inclination was fixed at 20°.

We used a Herschel-Bulkley model [19], to match therheology of the Carbopol used by Bates & Ancey [18].
In tensor form, the shear stresses τ can be written as:𝛾̇ = 0if |𝜏| < 𝜏0}(1)
𝜏0𝜏 = ( + 𝑘|𝛾̇ |𝑛−1 ) 𝛾̇ if |𝜏| ≥ 𝜏0|𝛾̇ |where τ0 is the yield stress [Pa]; k is the consistency[Pa⋅sn], 𝛾̇ is the shear rate [s-1], and n represents the flowindex [-]. Eqn. (1) represents a viscoplastic material –below the τ0, the material behaves as though it is a rigid
solid (i.e. “infinitely” viscous), and above τ0, viscous
flow can occur. Note that below τ0, the equation isindeterminate, and specific techniques are needed toproduce a solution. Furthermore, note that values ofn < 1 correspond to shear-thinning liquids.
For the flow, we set τ0 = 58 Pa, k = 35 Pa⋅sn andn = 0.33 to match the parameters used in [18]. For thebed, we used the same values of k and n, but varied τ0 in
the range 10 < τ0 < 1000 Pa, covering beds initially near
the solid-fluid transition, to those which are rigid.
2.2 Numerical solver
We use a weakly-compressible Smoothed ParticleHydrodynamics (SPH) solver [20], given its ability for
dealing efficiently and intuitively with (i) large
deformations and (ii) interactions between multiple
phases. We specifically adopted the GPU-parallelised
open-source SPH solver DualSPHysics [21].
The main underlying equations for the SPH are theNavier-Stokes equations. For our case, conservation ofmass and momentum can be expressed as:

2.4 Boundary conditions
The conditions for the fixed boundaries were tricky toimplement, as non-Newtonian flows feature an approximate no-slip condition at the base, particularly as theflows undergo creeping motion. To approximate a noslip condition, approximating the ‘stickiness’ between
the flowing boundary layer and the rigid boundary, we2

E3S Web of Conferences 415, 01008 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341501008

Fig. 1. Numerical setup: (a) shows the initial setup for thevery first (dam-break) simulation; (b) shows qualitatively
how the flow was frozen just before reaching the non-fixed
bed, so as to enable mass-exchange simulations to beundertaken with varying rheological parameters for the bed.

Fig. 2. Snapshots of a flow on a fixed bed.

generated the boundaries using several layers of SPH
material points, and ascribed a very high yield stress tothem (with τ0 of around 1000 Pa, determined from asensitivity study). This method is not ideal, because thetimestep for the simulations then becomes governed by
the boundaries. However, at the time of undertaking thesimulations, we were not aware of alternative/more
efficient numerical methods for simulating a no-slip
condition within DualSPHysics.

3 Preliminary results
Figs. 2 to 4 show the kinematics of flows on a variety ofbeds. Each screenshot shows the flow highlighted by (i)the different phases and (ii) the local velocity. Snapshotsat t = 2 and t = 8 s are shown for each case.
Fig. 2 shows a flow on a ‘fixed’ bed, i.e. where theyield stress of the bed was set to be around 1 kPa. The
velocity of the flow does not change substantially fromwhen the flow front arrives at t = 2 s, compared to further
on at t = 8 s. The velocity stays around 4 cm/s.
Fig. 3 shows a flow with identical starting conditionsto the one in Fig. 2. This time, the flow is passing overa bed with a yield stress τ0 equal to that of the flow. Att = 2 s, just after the flow has arrived on the bed, it canalready be seen that some deceleration is taking place,with the frontal moving at under 3 cm/s. Further
deceleration is apparent at t = 8 s, with the plug moving
at less than 2 cm/s. Notably, the screenshot at t = 8 sseems to show the material originally from the bed being
pushed and behaving as a mobile dam.
Fig. 4 shows another flow, this time passing over abed where the yield stress has been set low enough that
it is on the point of yielding under its own weight. In
other words:
𝜏0(5)≈1𝜌𝑔ℎ sin 𝜃
Here we see that the initial velocity remains at ornear 4 cm/s at t = 2 s. At t = 8 s, we see that the flow has

Fig. 3. Snapshots of a flow on a bed with depth 6 mm andyield stress of 58 Pa (i.e. the same as the flow).

already reached the far end of the numerical domain(with a velocity that is in excess of 4 cm/s). There is a“lubrication” effect caused by the underlying bed having
a low yield stress (see also [23]). The arrival of the flow
kick-starts it; the rheology being shear-thinning, itmeans that the bed can then start to flow very rapidly,
even though it was initially at rest. This layer of material
with a low yield stress enables the flow material to loseless energy, as compared to interfacing directly with theno-slip base (see Fig. 2).
Fig. 5 summarises the link between the bed yieldstress and the “flow reach”. The “flow reach” is definedas the frontmost point of the material that originallyconstituted the flow (i.e. material belonging originallyto the bed is always ignored). Data for a single timestep
is shown, i.e. t = 7.5 s. The grey region corresponds tothe mobile bed. The non-linearity of the mass-exchange
processes is revealed even by the very simple metric of3

E3S Web of Conferences 415, 01008 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341501008

even for highly idealised flows and beds. Future work isneeded to to convert the phenomena identified in thisstudy into quantitative relationships that can be
implemented into depth-averaged solvers.
This study used funding from the MOPGA 2021/22 scheme(MOPGA-976501H), from the ANR. The authors thank Dr
Georgios Fourtakas of the University of Manchester well as
Suzanne Lapillonne of INRAE for their input.

References
1.
2.

M. J. Froude, D. N. Petley. NHESS 18 (2018)
S. McDougall, O. Hungr. Can. Geotech. J. 41(2004)
3. S. R. Goodwin, C. E. Choi. Comp. & Geotech.
141, 104503 (2022)
4. C. E. Choi, M. Nikooei. Comp. Geotech 154(2022)
5. M. Nikooei, C. E. Choi. Comp. Geotech 148(2022)
6. M. Nikooei, M. T. Manzari. Comp. Geosci. 134(2020)
7. S. P. Pudasaini, J.-T. Fischer. Int. J. Mult. Flow132 (2020)
8. S. P. Pudasaini, Krautblatter. Nat. Comm. 12(1)(2021)
9. S. Cuomo, M. Pastor, V. Capobianco, L. Cascini.
Eng. Geol. 212 (2016)10. S. McDougall, O. Hungr. Can. Geotech. J. 42(2005)
11. S. Cuomo, M. Pastor, L. Cascini, G. C. Castorino.
Can. Geotech. J. 51, 11 (2014)12. H. Chen, G. B. Crosta, C. F. Lee. Géotechnique,56, 5 (2006)
13. F. Frank, B. W. McArdell, C. Huggel, A. Viel.
NHESS 15 (2015)14. H. Y. Hussin, B. Q. Luna, C. J. van Westen, M.
Christen, J.-P. Malet, T. W. J. van Asch NHESS
12 (2012)
15. G. B. Crosta, S. Imposimato, D. Roddeman. Eng.
Geol. 109 (2009).
16. K. Hu, P. Li, X.-k. Wang. J. Mount. Sci. 13, 6(2016)17. P. Song, C. E. Choi. JGR: Earth Surf. 126, 5(2021)
18. B. M. Bates, C. J. Ancey. Non-Newt. Fluid Mech.
243 (2017)19. W. H. Herschel, R. Bulkley. Kolloid Zeitschrift 39(1926)
20. J. J. Monaghan. Comp. Phys. Comm. 48 (1988)
21. A. J. C. Crespo, J. M. Dominguez, B. D. Rogers,
M. Gomez-Gesteira, S. Longshaw, R. Canelas, R.
Vacondio, A. Barreiro, O. Garcia-Feal. Comp.
Phys. Comm. 187 (2015)
22. A. Putz, I. A. Frigaard, D. M. Martinez. J. nonNewt. Fluid Mech 163, 1–3 (2009)
23. R. M. Iverson. J. Geophys. Res.: Earth Surf. 110,F2 (2005)

Fig. 4. Snapshots of a flow on bed with depth 6 mm and yieldstress of 18 Pa (i.e. almost yielding under its own weight).

Fig. 5. Position of the flow (as defined by phase) for a variety
of bed depths and bed yield stresses.

“flow reach”: at the low end of τ0 for the bed, the flowis almost at the end of the bed, due to the aforementioned
lubrication effect; the flow reach is minimum for aroundτ0 = 80 Pa, due to the strong transfer of momentum into
the bed for this value; and the flow reach starts toincrease again for higher values of τ0, sinceprogressively less momentum is being lost into the bed.
A similar trend is shown for all of the bed depths plotted;
this is reminiscent of (but not identical to) findings from[17], albeit where unsaturated geo-materials were used.
Other types of non-linearity are present as well.
Specifically, the minimum flow reach on the deeperbeds (D = 18 or 30 mm) is less than the shallower beds
(D = 6 or 12 mm). This is related to the distance that theflow material is able to fall on the deeper beds, which isa function of both the yield stress and depth of the bed.

4 Conclusions
In this study, we implemented a regularised HerschelBulkley model into the open-source SPH solver
DualSPHysics. We showed results for flows overriding
mobile beds, wherein various types of mass-exchange
processes were observed. Several regimes are possible,

4