# Mass exchange between geophysical flows and beds: Idealised computational modelling using a Herschel-Bulkley rheology

Saoirse Robin Goodwin, Guillaume Piton, Guillaume Chambon

► To cite this version:

Saoirse Robin Goodwin, Guillaume Piton, Guillaume Chambon. Mass exchange between geophysical flows and beds: Idealised computational modelling using a Herschel-Bulkley rheology. 8th International Conference on Debris Flow Hazard Mitigation (DFHM8), 2023, Torino, Italy. pp.01008, 10.1051/e3sconf/202341501008. hal-04687601

HAL Id: hal-04687601

https://hal.inrae.fr/hal-04687601v1

Submitted on 5 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mass exchange between geophysical flows and beds: idealised computational modelling using a Herschel-Bulkley rheology

Saoirse Robin Goodwin<fnref n="1" />,<fnref n="2" />*, Guillaume Piton<fnref n="1" /> and Guillaume Chambon<fnref n="1" />

**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 some consequences 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 the state of the art for modelling entrainment is essentially still empirical. In this study, we implement a Herschel-Bulkley (non-Newtonian) rheology into an existing open-source Smoothed Particle Hydrodynamics solver (DualSPPhysics). This rheology can reasonably represent clay-rich flows, typical of those observed in the French Prealps. We hence undertake a highly-idealised, quantitative investigation of entrainment mechanisms for flows overriding non-fixed beds. For the beds, we vary the yield stress and the depth. 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 (but not 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 climate change, as well as a continued propensity for human developments in mountainous areas worldwide [1]. Ultimately, this means that engineers and scientists increasingly need to identify which areas are likely most at risk, and to develop countermeasures accordingly.

One key tool for doing so is explicit numerical modelling of potential flows, using approaches based on fluid mechanics. In general, numerical solvers implement expressions describing conservation of both mass and momentum, either using the full Navier-Stokes equations, or the Saint-Venant (shallow water/depth-averaged) equations. Various constitutive models can then be implemented, for instance by modifying 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 of flows, as well as mass exchange mechanisms (entrainment and deposition) that are well-known to cause 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 full Navier-Stokes equations – they are of high interest to practising 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 of mass that includes a term for mass exchange *b*. The sign of *b* distinguishes entrainment and deposition.

The crux of the matter rests on how we should be computing 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 was a fine starting point, solvers were essentially blind to the local 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 earlier works, empirical relationships tend to be site-specific, and in any case rely on the availability of statistically-significant 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 and the bed, thus circumventing the need for ad-hoc relations. Nonetheless, such expressions for *b* would require evaluation against experimental or computational data for mass-exchange processes (the latter being explored in this Extended Abstract).

It is also worth noting that up until relatively recently, depth-averaged solvers have not accounted for the effects of mass exchange on expressions for conservation of momentum. This, as well as the state of the art for mass-exchange problems continuing to rely on empirical relationships, belies a fundamental lack of understanding of mass-exchange mechanisms.

<footnote n="1">¹Univ. Grenoble-Alpes, INRAE, UR ETNA, 38402 Saint-Martin-d'Hères, France</footnote>

<footnote n="2">²Karlsruher Institut für Technologie, 76131 Karlsruhe, Germany</footnote>

<footnote>* Corresponding author: srgoodwin@protonmail.com</footnote>
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 (non-Newtonian) rheology, to model both flows and non-fixed beds. Non-Newtonian rheologies constitute a reasonable simplification of clay-rich flows and advantageously comprise relatively few parameters. Bates & Ancey studied a range of bed depths and lengths, 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 we used 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 work of Bates & Ancey. The objectives of the parametric study are (i) to delineate regimes for different behaviour, and hence (ii) to provide a basis for developing/evaluating physically-based expressions for b.

# 2 Numerical model

## 2.1 Rheological model

We used a Herschel-Bulkley model [19], to match the rheology of the Carbopol used by Bates & Ancey [18]. In tensor form, the shear stresses $\underline{\tau}$ can be written as:

$$ \begin{aligned} \dot{\underline{\tau}} &= 0 & \text{if } |\underline{\tau}| < \tau_0 \\ \underline{\tau} &= \left( \frac{\tau_0}{|\dot{\underline{\tau}}|} + k|\dot{\underline{\tau}}|^{n-1} \right) \dot{\underline{\tau}} & \text{if } |\underline{\tau}| \ge \tau_0 \end{aligned} \quad (1) $$

where $\tau_0$ is the yield stress [Pa]; k is the consistency [Pa·sⁿ], $\dot{\underline{\tau}}$ is the shear rate [s⁻¹], and n represents the flow index [-]. Eqn. (1) represents a viscoplastic material – below the $\tau_0$, the material behaves as though it is a rigid solid (i.e. “infinitely” viscous), and above $\tau_0$, viscous flow can occur. Note that below $\tau_0$, the equation is indeterminate, and specific techniques are needed to produce a solution. Furthermore, note that values of n < 1 correspond to shear-thinning liquids.

For the flow, we set $\tau_0 = 58$ Pa, $k = 35$ Pa·sⁿ and $n = 0.33$ to match the parameters used in [18]. For the bed, we used the same values of k and n, but varied $\tau_0$ in the range $10 < \tau_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 Particle Hydrodynamics (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 the Navier-Stokes equations. For our case, conservation of mass and momentum can be expressed as:

$$ \frac{\partial \rho}{\partial t} = -\rho \underline{\nabla} \cdot \underline{U} \quad (2) $$

$$ \frac{\partial \underline{U}}{\partial t} = \frac{1}{\rho} \left( -\underline{\nabla} P + \underline{\nabla} \cdot \underline{\tau} \right) + \underline{g} \quad (3) $$

where $\rho$ is density, $P$ is pressure, $\underline{\tau}$ is the extra-stress tensor and $\underline{g}$ is acceleration due to gravity.

To implement the Herschel-Bulkley rheology, we used the well-established technique of regularisation [22], wherein we compute an equivalent viscosity $\mu_{eff}$ which can be capped arbitrarily at a threshold $\mu_{max}$. (We set $\mu_{max} = 10$ kPa, following a sensitivity study.) The expression for the rheology can thus be re-cast as:

$$ \left. \begin{aligned} \mu_{eff} &= \mu_{max} & \text{if } \mu_{eff} \ge \mu_{max} \\ \mu_{eff} &= \frac{\tau_0}{|\dot{\underline{\tau}}|} + k|\dot{\underline{\tau}}|^{n-1} & \text{if } \mu_{eff} < \mu_{max} \end{aligned} \right\} \quad (4) $$

The advantages of this approach are twofold: it means that $\mu_{eff}$ is determinate everywhere, and the minimum timestep required to solve the simulation can be limited (since the timestep is inversely related to viscosity terms). However, a limitation of the approach is that true unyielded zones cannot be resolved. In other words, even material where $\mu_{eff} = \mu_{max}$ can continuously yield. The rate of this yielding is determined by $\mu_{max}$, and should ideally be as low as possible, to approximate “true” unyielded material.

## 2.3 Numerical setup

The numerical setup for our study is modelled on that 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 a certain distance downstream from the gate. For all of the simulations presented in this Extended Abstract, the length 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: it greatly reduced the amount of time required for each simulation, and it also reduced the chances of numerical instabilities manifesting.

Fig. 1a shows the starting condition for the dam-break simulation. We stopped this simulation before the flow material reached the non-fixed bed. Each of the mass-exchange simulations started by regenerating the flow material according to the flow bounds, as well as its bulk velocity (see Fig. 1b).

The density of the flow material was set at 1000 kg/m³, but this value could vary slightly given that we were using a slightly-compressible formulation. The channel inclination was fixed at 20°.

## 2.4 Boundary conditions

The conditions for the fixed boundaries were tricky to implement, as non-Newtonian flows feature an approximate no-slip condition at the base, particularly as the flows undergo creeping motion. To approximate a no-slip condition, approximating the ‘stickiness’ between the flowing boundary layer and the rigid boundary, we
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Numerical setup: (a) shows the initial setup for the very 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 be undertaken with varying rheological parameters for the bed.</figcaption></figure>

generated the boundaries using several layers of SPH material points, and ascribed a very high yield stress to them (with $\tau_0$ of around 1000 Pa, determined from a sensitivity study). This method is not ideal, because the timestep for the simulations then becomes governed by the boundaries. However, at the time of undertaking the simulations, 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 of beds. Each screenshot shows the flow highlighted by (i) the different phases and (ii) the local velocity. Snapshots at $t=2$ and $t=8$ s are shown for each case.

Fig. 2 shows a flow on a ‘fixed’ bed, i.e. where the yield stress of the bed was set to be around 1 kPa. The velocity of the flow does not change substantially from when 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 conditions to the one in Fig. 2. This time, the flow is passing over a bed with a yield stress $\tau_0$ equal to that of the flow. At $t=2$ s, just after the flow has arrived on the bed, it can already 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$ s seems 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 a bed where the yield stress has been set low enough that it is on the point of yielding under its own weight. In other words:

$$ \frac{\tau_0}{\rho g h \sin \theta} \approx 1 \qquad (5) $$

Here we see that the initial velocity remains at or near 4 cm/s at $t=2$ s. At $t=8$ s, we see that the flow has

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Snapshots of a flow on a fixed bed.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Snapshots of a flow on a bed with depth 6 mm and yield stress of 58 Pa (i.e. the same as the flow).</figcaption></figure>

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, it means 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 lose less energy, as compared to interfacing directly with the no-slip base (see Fig. 2).

Fig. 5 summarises the link between the bed yield stress and the “flow reach”. The “flow reach” is defined as the frontmost point of the material that originally constituted the flow (i.e. material belonging originally to the bed is always ignored). Data for a single timestep is shown, i.e. $t=7.5$ s. The grey region corresponds to the mobile bed. The non-linearity of the mass-exchange processes is revealed even by the very simple metric of
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Snapshots of a flow on bed with depth 6 mm and yield stress of 18 Pa (i.e. almost yielding under its own weight).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Position of the flow (as defined by phase) for a variety of bed depths and bed yield stresses.</figcaption></figure>

"flow reach": at the low end of $\tau_0$ for the bed, the flow is almost at the end of the bed, due to the aforementioned lubrication effect; the flow reach is minimum for around $\tau_0$ = 80 Pa, due to the strong transfer of momentum into the bed for this value; and the flow reach starts to increase again for higher values of $\tau_0$, since progressively 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 deeper beds ($D$ = 18 or 30 mm) is less than the shallower beds ($D$ = 6 or 12 mm). This is related to the distance that the flow material is able to fall on the deeper beds, which is a function of both the yield stress and depth of the bed.

## 4 Conclusions

In this study, we implemented a regularised Herschel-Bulkley 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,

even for highly idealised flows and beds. Future work is needed to convert the phenomena identified in this study 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. M. J. Froude, D. N. Petley. NHESS **18** (2018)

2. 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. Flow **132** (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. non-Newt. Fluid Mech **163**, 1-3 (2009)

23. R. M. Iverson. J. Geophys. Res.: Earth Surf. **110**, F2 (2005)