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

. 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 (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 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.


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][5][6][7][8][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, * Corresponding author: srgoodwin@protonmail.com 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][13][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 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 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.
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 nonfixed 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.

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 τ can be written as: where τ0 is the yield stress [Pa]; k is the consistency [Pa⋅s n ], ̇ is the shear rate [s -1 ], and n represents the flow index [-]. Eqn. (1) represents a viscoplastic materialbelow 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 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 τ0 = 58 Pa, k = 35 Pa⋅s n 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 τ0 in the range 10 < τ0 < 1000 Pa, covering beds initially near the solid-fluid transition, to those which are rigid.

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: (3) where ρ is density, P is pressure, τ is the extra-stress tensor and 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 µeff which can be capped arbitrarily at a threshold µmax. (We set µmax = 10 kPa, following a sensitivity study.) The expression for the rheology can thus be re-cast as: The advantages of this approach are twofold: it means that µ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 µ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.

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 dambreak 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 3 , but this value could vary slightly given that we were using a slightly-compressible formulation. The channel inclination was fixed at 20°.

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 noslip condition, approximating the 'stickiness' between the flowing boundary layer and the rigid boundary, we https://doi.org/10.1051/e3sconf/202341501008 , 01008 (2023) E3S Web of Conferences 415 DFHM8 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. generated the boundaries using several layers of SPH material points, and ascribed a very high yield stress to them (with τ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.

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 τ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: 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  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  "flow reach": at the low end of τ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 τ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 τ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.