## Abstract

A theoretical and experimental investigation of the propagation of free-surface, channelized viscous gravity currents is conducted to examine the combined effects of fluid rheology, boundary geometry and channel inclination. The fluid is characterized by a power-law constitutive equation with behaviour index *n*. The channel cross section is limited by a rigid boundary of height parametrized by *k* and has a longitudinal variation described by the constant *b*≥0. The cases *b*>0, the cases *b*=0 implies a constant cross-sectional channel. For a volume of released fluid varying with time like *t*^{α}, scalings for current length and thickness are obtained in self-similar forms for horizontal and inclined channels/fractures. The speed, thickness and aspect ratio of the current jointly depend on the total current volume (*α*), the fluid rheological behaviour (*n*), and the transversal (*k*) and longitudinal (*b*) geometry of the channel. The asymptotic validity of the solutions is limited to certain ranges of parameters. Experimental validation was performed with different fluids and channel cross sections; experimental results for current radius and profile were found to be in close agreement with the self-similar solutions at intermediate and late times.

## 1. Introduction

Viscous buoyancy-driven gravity currents are observed in a wide range of natural phenomena and industrial processes [1–6]. The dynamics of these flows were first analysed considering two-dimensional or radial currents of Newtonian fluids flowing on a horizontal or sloping bottom [7,8]. More recently, a variety of applications have suggested exploring the influence of fluid rheology and geometry of the flow domain on the propagation of the current.

Consideration of non-Newtonian effects is essential for studying several environmental, industrial and biological phenomena, such as the spreading of muds in submarine settings [9]; the motion of natural suspensions containing sand, silt, clay and organic fractions [10,11]; debris flows in mountainous areas [12]; the propagation of pyroclastic flows produced by volcanic eruptions [13,14]; magma flow at sub-liquidus temperatures, owing to gas bubbles and to the presence of crystals [15]; the disposal of mine tailings in the minerals industry [16]; the flow of biomass slurries [17]; blood motion in arterioles [18]. In a number of instances, the fluid rheological behaviour is appropriately described by a Cross or Carreau–Yasuda model, but can be approximated over a certain range of shear rates by the simpler Ostwald–de Waele power-law model [19] if the yield stress is limited. Viscous gravity currents of power-law fluids in unidirectional planar or radial geometries have been investigated theoretically by several authors [20–23]. These studies were mainly aimed at obtaining the rate of propagation over time of a current generated by the release of a constant or time variable volume of fluid. Their results were substantiated by the experiments conducted by various researchers [24–26]. The accumulated body of knowledge allows, in principle, using controlled experiments as rheometers to determine the fluid parameters [26,27].

Gravity currents frequently advance within confining boundaries instead of freely flowing on a flat surface. This is the case for lava [28] or mud flows [29], the first acting by melting the substrate or by creating levees, the second by eroding the substrate or by flowing in previously generated channels. Extensive studies of channelized gravity currents [30,31] have demonstrated that the shape of the cross section influences the velocity of the front and modifies the overall behaviour of the current. Significant qualitative differences were shown between wide sections, termed channels, that mainly extend in the horizontal direction and narrow sections, termed fractures, that mainly extend in the vertical direction. A further generalization of the scheme [32] included longitudinal cross-sectional variations, with the channel transversal length scale increasing with the distance from the source and producing widening channels or squeezing fractures. This model mimics natural scenarios for lava flows [33], which are strongly dependent on the volume of magma released during eruption, the effusion rate and the cross-sectional variation of the channel [34]. Another important application is flow in rock fractures. Real fractures exhibit a complex geometry with multiple intersecting channels and preferential flow paths [35], yet if there is a decreasing trend in the aperture field, a single fracture with a squeezing cross section may encapsulate the overall system behaviour.

To the best of our knowledge, we are not aware of studies combining non-Newtonian effects with longitudinal variations of the cross section for free-surface gravity currents. Similar studies investigated the effect of non-Newtonian rheology on porous gravity currents in media with spatially varying permeability [36]. Hence, the motivations of this study are to develop closed-form solutions for free-surface gravity currents of power-law fluids flowing in horizontal or inclined widening channels or squeezing fractures with gently varying section and to test the models experimentally. This approach constitutes a generalization of the results derived by Takagi & Huppert [32] on one hand, and by Longo *et al.* [37] on the other. With respect to Takagi & Huppert [32], this study (i) deals with non-Newtonian power-law fluids, an extension which covers several cases of practical interest; (ii) presents an in-depth discussion of the model limits of applicability and of the interrelations among its parameters; (iii) reports on novel experiments with gently widening channels (which also test the Newtonian fluid case treated by Takagi & Huppert [32]). With respect to Longo *et al.* [37], the analysis is extended to longitudinally varying cross sections, a quite common feature of natural channels in realistic conditions. Different combinations of channel shape, longitudinal variation and fluid rheology mimic an array of situations of environmental interest.

For horizontal channels and fractures, we derive self-similar solutions of the first-kind, with the exponents of the self-similar variables computed by dimensional analysis (scaling). We need to mention that for converging horizontal channels and squeezing horizontal fractures a self-similar solution of the second-kind is suggested [38], with the exponents computed as a solution of a nonlinear eigenvalue problem in the frame of a phase-plane formalism [39]. However, the correct solution (self-similar of the first-kind or of the second-kind) should be tested with experiments which are not trivial and are left for future work. Our experiments confirm the proposed solutions for horizontal and inclined diverging channels.

In §2, governing equations are derived for gravity currents of time-varying volume in rectilinear channels or fractures with longitudinally varying cross sections. Scalings for current length and thickness are obtained in self-similar forms for horizontal and inclined channels/fractures. Section 3 discusses the dependence of the main features of the current on the problem parameters and derives the bounds on parameters necessary to adhere to the model assumptions. In §4, theoretical results are compared with laboratory experiments conducted at constant volume or influx rate with shear-thinning, Newtonian and shear-thickening fluids in a linearly divergent circular section. A set of conclusions (§5) closes the paper.

## 2. Formulation

Consider a rectilinear channel inclined at an angle *β* with the horizontal plane (figure 1) and having a cross section whose rigid boundary is described by a power function *d*(*x*,*y*)=*a*|*y*/*a*|^{k}, where *a* is the length scale controlling the width and the positive constant *k* characterizes the shape [32] (note that their *n* is our *k*). The cases *k*>1 and *k*<1 correspond to wide channels and narrow fractures, respectively, the former mainly extending in the horizontal direction and the latter mainly extending in the vertical direction. Relevant special cases are (i) the symmetric triangular section (*k*=1) of vertex angle 2*θ*, which is wide for *m*≪1, and belongs to neither category for generic values of *m*; (ii) the circular section of assigned radius *a*/2, which locally approximates the case *k*=2 if the width of the current is limited compared to the radius and (iii) the rectangular section of width 2*a*, corresponding to *x* along the channel axis by assuming *a*(*x*)=*rx*^{b}, with *r* a positive constant and *b* a non-negative constant. The width variation per unit length of channel is moderate, so that the normal stress at the lateral wall has a negligible component along the flow direction. The cases *k*>1 and *k*<1 with *b*>0 (figure 2) describe widening channels and squeezing fractures, respectively; the sub-case *b*=0 identifies a channel or fracture of constant cross section. A free-surface current of a non-Newtonian fluid with uniform density *ρ* and having a total volume varying with time *t* as *t*^{α} is released at *t*=0 at the upper end of the channel. After some time from the initiation of the motion, the width 2*W*(*x*,*t*) and height *h*(*x*,*t*) of the current become small when compared with its length *x*_{N}(*t*), allowing us to consider its one-dimensional motion along the channel axis *x*, described by the single velocity component *u*(*x*,*y*,*z*,*t*). At the same time, any existing height differences in the transverse direction-*y* are smoothed out, and the current height becomes uniform across the channel. Further, viscous forces are assumed to prevail over inertia. Neglecting capillary effects, the pressure distribution in the intruding current is hydrostatic and given by *z* is the vertical coordinate, *p*_{0} the atmospheric pressure at the free surface and *g* gravity. Under the previous assumptions, the momentum equation reduces to
*τ*_{yx} and *τ*_{zx} are the shear stresses. Inside wide channels, it is assumed that the current height is much smaller than its width (*h*≪*W*∼*a*), so that ∂*τ*_{zx}/∂*z*≫∂*τ*_{yx}/∂*y*. Inside narrow fractures, the opposite is true as *h*≫*W*. It is thus seen that in either case one term in equation (2.1) is negligible, thus simplifying its formulation. We note in passing that either assumption remains valid along the current length under the hypothesis *b*≥0; considering the case *b*<0 would lead to initially wide channels narrowing with *x*, or to initially narrow fractures widening with *x*.

The rheological equation for a power-law fluid with simple shear is *n* the fluid behaviour index (*n*<1 for shear-thinning, *n*=1 for Newtonian, *n*>1 for shear-thickening fluids), and *u*/∂*z*>0 or ∂*u*/∂*y*>0 and taking the hydrostatic hypothesis into account, yields for wide and narrow sections, respectively

In deriving equation (2.2), the additional assumption was made that motion is induced by the free surface slope for a horizontal channel (*β*=0) and solely by the down-slope component of gravity for an inclined channel (*β*≠0). The latter requirement implies *β* is small enough to make the two actions comparable [30]. Hence the source term *S* takes the two alternative forms
*A* occupied by the fluid [30,37], yields
*Q* is the flow rate. Finally, the global continuity equation requires
*q*>0 and *α*≥0 are constants; *α*=0 indicates the instantaneous release of a given fluid volume, and *α*=1 indicates the continuous injection of a fluid at a constant flow rate. For *β*=0, the boundary condition at the current front is
*β*≠0, as the order of the differential equation of motion is lower; in this case, the effect of surface tension, neglected in the present model, brings the current height to zero at the front.

### (a) Widening channels

Wide channels (*k*>1) that retain their shape (*b*=0) or gently widen with distance from the source (*b*>0) are described by the first equation in (2.2) with boundary conditions expressing the no-slip condition at the channel boundary and the null value of the shear stress at the free surface, i.e.

#### (i) Horizontal channels

In the horizontal case (*β*=0), equation (2.4) yields
*H*=*h*/*x**, *X*=*x*/*x**, and *T*=*t*/*t** are dimensionless variables obtained upon introducing the following length and time scales:
*H*[*X*_{N}(*T*),*T*]=0 and
*η*_{N} is the value of *η* at *X*_{N}, *ζ*=*η*/*η*_{N} the reduced similarity variable, and the shape function *ψ* is the solution of
*ζ*. The global continuity equation (2.16) then becomes

For *α*=0, equations (2.19) and (2.20) are amenable to a closed-form solution given by
*n*=1 and *b*=0, are equivalent to the expressions given by Takagi & Huppert [30]. For *α*≠0, equations (2.19) and (2.20) are integrated numerically; see [37] for details. Equation (2.17) indicates that for the horizontal case (*h*) and widening cross sections (*w*), the current front propagates as *X*_{N}∝*T*^{chw} with *c*_{hw}≡*F*_{1} and a velocity *U*_{XN}∝*T*^{chw−1}. For a Newtonian fluid (*n*=1), *c*_{hw} reduces to that derived by Takagi & Huppert [32], their eqn (15). For a channel of constant cross section (*b*=0), the exponent *c*_{hw} reduces to that predicted by Longo *et al.* [37], their eqn (63). Figure 3 shows the shape function *ψ*(*ζ*) and prefactor *η*_{N} for a diverging channel with *k*=2 and *b*=1, for constant volume (*α*=0) and flux (*α*=1). The values of the shape function become larger as *α* increases and *n* decreases, while the prefactor decreases with *α* and increases with *n*. Similar plots for different values of *k* and *b* (not shown here) indicate that (i) decreasing *k* brings about modest changes in *ψ*(*ζ*) and *η*_{N}; (ii) for a nearly triangular channel (*ψ*(*ζ*) and *η*_{N} are insensitive to the exponent *b*; and (iii) changes in the exponent *b* bring about a decrease in *η*_{N}, whereas *ψ*(*ζ*) is only slightly affected.

#### (ii) Inclined channels

In the inclined case (*β*≠0), equation (2.4) becomes in dimensionless form
*X*_{0} indicating the coordinate where *H*=0. This is clearly unphysical since several model hypotheses are violated at *X*=*X*_{0}; in fact, equation (2.25) describes the profile of the current for *X*_{0}≪*X*<*X*_{N}. At *X*=*X*_{N}, we assume that the profile is ending abruptly and can be smoothed by including the surface tension [7]. Note that no further boundary condition is needed.

The constraint represented by equation (2.10) gives the length of the current for *X*_{N}≫*X*_{0} in inclined (*i*) widening (*w*) channels as
*U*_{XN}∝*T*^{ciw−1}. For *n*=1, the results of Takagi & Huppert [32] are recovered. For a channel of constant cross section (*b*=0), the exponent *c*_{iw} reduces to that predicted by Longo *et al.* [37], their eqn (64).

### (b) Squeezing fractures

Narrow fractures (*k*<1), which retain their shape (*b*=0) or gently squeeze with distance from the source (*b*>0), allow neglecting the term ∂*τ*_{zx}/∂*z* in equation (2.2), which is then integrated, imposing
*W*(*z*)=*z*^{1/k}*r*^{(1−1/k)}*x*^{b(1−1/k)}, 0<*k*<1, obtaining the velocity as

Substituting the above equation in equation (2.4) and using the constraint represented by the mass conservation (2.5), dimensional arguments show that the front end advances as *X*_{N}∝*T*^{c}, with the exponent *c* given for horizontal (*h*) or inclined (*i*) squeezing (*s*) fractures by
*n*=1 and to those predicted by the eqns (65) and (66) of [37] for *b*=0.

### (c) Channels with symmetric triangular cross section

The symmetric triangular section with an angle 2*θ* at the vertex is the dividing case between wide and narrow sections, yet it cannot be analysed as a limiting case for *τ*_{zx}/∂*z* and ∂*τ*_{yx}/∂*y* in (2.1) are of the same order. For a Newtonian fluid (*n*=1), a general closed-form solution for channels with a spatially constant cross section (*b*=0) is available, providing the velocity and the flow rate in the form of an infinite series; the flow rate is approximated by the analytic expression *Q*≈(0.137*m*^{3})*h*^{4}*S*/(1+*m*^{2}) [30], where *S* is the source term defined in (2.3). Using this expression, a self-similar solution was shown to exist [30] for *b*=0; on the contrary, for widening channels (*b*>0), equation (2.4) does not in general allow a self-similar solution, which exists for the two asymptotic cases of nearly flat (*m*≫1) or very narrow sections (*m*≪1) [32].

Nonlinearity prevents a closed-form solution of this type for *n*≠1. For channels of constant cross section (*b*=0), Longo *et al.* [37] adopted an empirical approach and showed that the flow rate can be expressed as
*K*_{t}(*n*,*m*) includes a dimensionless coefficient *C*=*C*(*m*) that is a function of the vertex angle 2*θ* and is determined experimentally. For a right triangular section with 2*θ*=90° and a power-law fluid, *C*=14.6 on the basis of the experiments by Burger *et al.* [40]. Expression (2.31) evaluated for *n*=1 gives a coefficient *K*_{t} approximately equal to the value 0.137 predicted by Takagi & Huppert [30]. For a different value of 2*θ*, the coefficient *C* in equation (2.31) will take a different numerical value, but the dependence of *K*_{t} on *m* and *n* will not be altered. However, as in the case of Newtonian fluids, a self-similar solution to the problem for channels with varying cross sections (*b*≠0) exists only for the two asymptotic cases of nearly flat sections (*m*≫1) or very narrow sections (*m*≪1). In these two cases, the coefficient *m*^{(2n+1)/n}/(1+*m*^{2})^{(n+1)/(2n)} in equation (2.31) reduces to ≈*m* and ≈*m*^{2+1/n}, respectively. Using these simplified expressions, substituting the discharge for a triangular section given by (2.31) in (2.4) and using the mass conservation constraint (2.5), the time rate of advancement *c* of the current front is given for horizontal (*h*) or inclined (*i*) channels of triangular (*t*) nearly flat ( *f*) cross section by
*n*) sections, it is equal to
*n*=1. Further, they coincide with equations (2.17), (2.26) and (2.30) when written for *k*=1 and *b*=0; for *b*>0, a discontinuity arises in the values of the time exponent. In nearly flat triangular channels, the propagation speed decreases with *b*; the reverse is true for narrow triangular fractures. The behaviour of these exponents can be evaluated in a similar manner as the general case (§3); for the narrow case, a limit value of *b* must not be exceeded for the exponent *c* to remain positive.

## 3. Discussion of results

An analysis of the theoretical results requires examining the dependence on model parameters of: (i) the rate of advancement *F*_{1}=*c*; (ii) the time exponent of the thickness *F*_{2} (equation (2.18)) and (iii) the difference between the two aforementioned factors, *F*_{2}−*F*_{1}, which describes the evolution over time of the current aspect ratio (and of the average free-surface gradient, driving the motion for horizontal currents). The factors *c*, *F*_{2}, and *F*_{2}−*F*_{1} are discussed in the following subsections, together with the respective limits placed on problem parameters by model assumptions. For porous gravity currents, this concept was earlier elucidated for flow in confining channels [41] and in variable-conductivity media [42].

### (a) Propagation rate

The current propagation rate *F*_{1}=*c* is summarized in tables 1 and 2. The *c* expressions obtained for widening channels and narrow fractures coincide for *k*=1. For each case, the current front accelerates or decelerates depending on whether the rate of increase of the volume, *α*, is smaller or larger than a critical value *α*_{a} relative to acceleration. Expressions for *α*_{a} are reported in tables 1 and 2 for widening channels and squeezing fractures, respectively. For horizontal channels or fractures, *α*_{a} depends on both *k* and *n*; for the inclined case, *α*_{a} is independent of rheology. For *b*=0, the values of *α*_{a} reduce to those predicted by Longo *et al.* [37]; their counterparts for *n*=1 are easily deduced from Takagi & Huppert [32]. Physically, for widening channels, *α*_{a} increases for increasing *b* because propagation is controlled by mass balance, and the storage capacity per unit length of the channel increases with *b*. For squeezing fractures, *α*_{a} also increases for increasing *b*, but this increase is due to the dominance of the viscous resistance. As the cross-section becomes narrower along the flow direction, viscous resistances increase, and therefore a larger inflow rate (*α*) is needed to sustain the current acceleration.

Next, we study the dependence of *c* on the current volume, described by *α*, on the geometry of the channel/fracture, parametrized by *k* and *b*, and on the fluid rheology, encapsulated in the behaviour index *n*. By considering the partial derivatives of equations (2.17), (2.26) and (2.30), it is first observed that *c* increases with *α* for any value of *n*, *k* and *b*. This result follows directly from mass balance considerations. Moreover, two further critical values emerge for each case, *α*_{k} and *α*_{n}, which render the speed of the current front independent of *k* (the channel/fracture shape) and of *n* (the fluid rheology), respectively. These values are also reported in tables 1 and 2.

We consider *α*_{k} first and note that it assumes different values for the horizontal and inclined cases but has the same value for widening channels and squeezing fractures with the same inclination. The critical value *α*_{k} depends on *n* and *b* and reduces to that predicted by Longo *et al.* [37] for a constant cross section (*b*=0); further considering a Newtonian fluid (*n*=1) yields *α*_{k} is defined only below a critical value *b*_{k} of *b*, listed in tables 1 and 2 as well; *b*_{k} depends only on the channel/fracture inclination but again assumes the same value for widening channels and squeezing fractures. The tendency of the current front to increase or decrease its velocity with *k* is governed by the interplay of the critical values *α*_{k} and *b*_{k} as reported in tables 1 and 2. Physically, this tendency reflects the relative importance of the competing effects of mass balance and viscous resistances. In channels that widen rapidly (*b*≥*b*_{k}), mass balance prevails and *c* decreases with *k* for any value of *α*. If the widening is moderate (*b*<*b*_{k}), then mass balance prevails for a moderate influx rate (*α*<*α*_{k}), while viscous resistances dominate above this threshold (*α*≥*α*_{k}). These considerations are reversed for squeezing fractures.

The second critical value *α*_{n} depends on *k* and *b* irrespective of fluid rheology; it reduces to 3 and 1, for the horizontal and inclined case, respectively, in wide rectangular channels of constant cross section. The current advances less rapidly for an increasing flow behaviour index *n* when *α*>*α*_{n}; the reverse is true for *α*<*α*_{n}, as reported in tables 1 and 2.

Finally, we observe that the current speed always decreases with increasing *b* for any *α*, *k*, *n*. This behaviour reflects the dominance of the mass balance effects in widening channels and of the viscous resistances in squeezing fractures.

To better grasp the behaviour of the propagation rate, we depict in figure 4*a* *c*_{h} as a function of *k* for different values of *b* and a Newtonian fluid (*n*=1) flowing in horizontal channels (*h*). Near *k*=1, the values used are qualitative due to the approximation involved in the expression of *c*_{h} for the triangular section. For *k*=1, a maximum value of the exponent *c*_{h} is attained when *b*<*b*_{k} for 0<*α*<*α*_{k}, or when *b*>*b*_{k} for any *α*≥0, and a minimum value is attained when *b*<*b*_{k} for *α*>*α*_{k}. In the special case *α*=*α*_{k}, *c*_{h} is independent of the shape of the channel. Moreover, *c*_{h} is a decreasing function of *b*, more so for lower values of *α*. A dependence of *c*_{h} on *α* and *k*, which is very similar to that illustrated in figure 4, is obtained for different values of *n* and is not shown here. For *n*<1, values of the exponent *c*_{h} are higher than those obtained for *n*=1 if *α*<*α*_{k} and lower if *α*>*α*_{k}; the reverse is true for *n*>1. The differences in *c*_{h} due to the variation of *n* are modest. Figure 4*b* shows the time exponent *c*_{i} as a function of *k* and different values of *b* for a Newtonian fluid (*n*=1) flowing in inclined channels/fractures (*i*); again, the diagrams are approximate near *k*=1. All of the qualitative observations made for the horizontal case hold true, albeit with higher values of the exponent due to the larger current speed. Figure 5 shows the possible combinations of parameters *n* and *b* giving rise to the fastest propagation rate for the triangular cross section when *α*<*α*_{k} and *b*<*b*_{k}.

Only strictly positive values of the propagation rate *c* are physically meaningful. Upon examining the values in tables 1 and 2, we observe that the propagation rate is positive for all parameter values; therefore, no restriction is placed on the problem parameters on this account.

### (b) Current thickness

The dependence of the rate of temporal variation of the current thickness *F*_{2} on model parameters is illustrated in table 3 for horizontal widening channels and squeezing fractures. The first row of table 3 summarizes the expressions for *F*_{2}, while the second row lists the conditions under which the thickness decreases or remains constant with increasing time. Only this case is asymptotically meaningful, as for positive values of *F*_{2}, the current eventually reaches the top of a given channel or fracture. For widening channels, *F*_{2}≤0 for *α*≤*α*_{t}, as the larger channel volume available due to widening allows the current thickness at a given point to diminish if the volume growth rate is below a certain threshold. For squeezing fractures, the analysis is more complex as a second threshold *b*_{t} becomes important, depending on the channel shape *k*. For *b*≥*b*_{t}, the current thickness increases with time irrespective of the value of *α*; for *b*<*b*_{t}, the current thickness does not increase with time only when *α*≤*α*_{t}. Physically, this behaviour occurs because in squeezing fractures the available volume shrinks along the flow direction; therefore, both the rate of growth of the current volume and the squeezing parameter must be limited for the current to become thinner over time.

Table 3 then shows that *F*_{2} increases with *α* irrespective of the other model parameters by virtue of mass balance. Furthermore, the dependency of *F*_{2} on *k* and *n* is governed by the critical values *α*_{k} and *α*_{n}; the first makes *F*_{2} independent on *k* and the second on *n*. The critical value *α*_{k} is defined only below a critical value of *b*, *b*_{k}, that is identical for widening channels and squeezing fractures. Analogously, *α*_{n} is defined only below a critical value *b*_{n} but only for squeezing fractures. Interestingly, the critical values *α*_{k}, *b*_{k} and *α*_{n} coincide with those derived earlier for the propagation rate *c*.

For a large widening or squeezing rate (*b*≥*b*_{k}), *F*_{2} decreases with increasing *k* for any *α*, reflecting the effect of the larger channel volume available for a larger *k*. For a moderate longitudinal variation rate (*b*<*b*_{k}), *F*_{2} decreases with *k* for a moderate influx (*α*≤*α*_{k}); the reverse is true for *α*≥*α*_{k}.

When the dependency on *n* is considered, *F*_{2} is observed to increase/decrease with *n* in widening channels for *α* above/below *α*_{n}. For squeezing fractures, a second critical value *b*_{n} becomes important, and *F*_{2} becomes independent of *n* when either *α* or *b* is equal to its critical value. *F*_{2} increases or decreases with *n* for different combinations of the ranges of parameters *α* and *b*. Finally, the exponent *F*_{2} decreases with *b* for widening channels and increases for squeezing fractures, as larger values of *b* imply an increase or decrease of available channel volume for these two cases, respectively.

The dependency of *F*_{2} on model parameters for the inclined case follows a similar pattern: in this case, *H*∼*T*^{−n/(n+1)} and *H*∼*T*^{−kn/(n+1)}, respectively, for widening channels and squeezing fractures. Therefore, the exponent *F*_{2} is always negative, and the thickness decreases with time for any combination of model parameters.

### (c) Current aspect ratio

A similar analysis can be readily developed for the time exponent of the current aspect ratio, *F*_{2}−*F*_{1}, based on the results of the previous subsections. The main result is that *F*_{2}−*F*_{1} is always negative except for the case of horizontal squeezing fractures. For this case, *F*_{2}−*F*_{1}<0 if (∀ *α*∧*b*<*b*_{g}) and (*α*<*α*_{g}∧*b*_{g}<*b*<(2*k*+1)/(1−*k*)), with *b*_{g}=[*n*(1−*k*)+1]/[(1−*k*)(*n*+1)]) and *α*_{g}=[*bn*(1−*k*)−*n*(1+2*k*)]/[*b*(*n*+1)(*k*−1)+*n*(1−*k*)+1]. Physically, this means that when the motion is driven by gravity, the current becomes more elongated with time irrespective of channel geometry and rheology. The same is true in horizontal widening channels, as the viscous resistances decrease along the flow direction, and a smaller free-surface gradient is needed to drive the flow. In horizontal squeezing fractures, the aspect ratio decreases with time only if the squeezing is modest, or for intermediate rates of squeezing and moderate influx rates. For pronounced squeezing, or for high influx rates with intermediate squeezing, *F*_{2}−*F*_{1} is positive, which renders the thin-current approximation asymptotically invalid. The validity of the thin-current approximation can also be checked at finite times according to, e.g. Ciriello *et al.* [43].

The analysis of the dependence of *F*_{2}−*F*_{1} on model parameters is similar to that developed for *F*_{1} and *F*_{2}; e.g. the time exponent of the aspect ratio always increases with *b*, which indicates that horizontal or inclined widening channels or squeezing fractures (*b*>0) experience a time-decreasing aspect ratio but at a smaller rate with respect to the constant cross-sectional case (*b*=0).

## 4. Experiments

The theoretical formulation was tested by conducting a series of experiments in horizontal (*β*=0) or inclined (*β*>0) channels for constant volume (*α*=0) or volume flux (*α*=1), and with Newtonian (*n*=1), shear-thinning (*n*<1) and shear-thickening (*n*>1) fluids. Experiments were conducted for circular channels (*k*=2) with a radius linearly increasing downstream (*b*=1), as widening channels are common in environmental applications. The case of non-Newtonian flow in channels of constant cross section (*b*=0) is covered by the experiments of Longo *et al*. [37].

The layout of the experimental apparatus is depicted in figure 6. A 200 cm long channel was built with a 0.75 mm thick Lexan sheet in the shape of a truncated half-cone, supported by a series of evenly spaced ribs. The inlet and outlet sections had radii *R*=5 cm and *R*=25 cm, respectively, hence *r*=0.2. For tests with *β*>0, the channel was lifted up at the upstream end and then fixed to obtain the desired inclination angle, which was measured with an electronic spirit level with an accuracy of 0.1°. A lock gate located 20 cm from the inlet section allowed for performing tests with *α*=0; the fluid volume was calculated by weighing the mass filling the upstream reservoir and dividing by the mass density, with a relative accuracy ±1.5%. A syringe pump, delivering the fluid through a plastic tube, or a Mariotte bottle, with a solenoid valve controlled by a personal computer, were employed to feed the channel for tests with a constant volume flux (*α*=1). The two devices were used for tests with Newtonian or shear-thinning fluids, respectively, for flow rates below and above 3.7 ml s^{−1}, with a relative accuracy of ±1%. Shear-thickening fluids were directly poured in the channel through funnels of different sizes, refilling the fluid manually to ensure a constant delivery rate. The mean flow rate was estimated by weighing the released fluid, with a relative accuracy ±1.5%.

The position of the current front was detected by employing the acquisition system and software described in Longo *et al.* [37], adapted to the experiments performed here. The corresponding accuracy was ±1.5 mm. In a few cases, the images were processed by hand, detecting the interface manually on the pictures when the fluid was not clearly distinguishable. In addition, for some tests, the current thickness in a single section was measured with an ultrasonic device, sampling at 100 Hz and with an overall accuracy within 0.3 mm. The accuracy in time measurements was taken to be 1/50 s or negligible depending on whether a video camera or a set of two photo cameras were used.

The Newtonian fluid used in the experiments was a mixture of water, glycerol and ink. Xanthan gum was added to obtain a shear-thinning fluid. The shear-thickening fluid was a mixture of water (45 wt%) and cornstarch (55 wt%). The rheological behaviour of all fluids employed was independently assessed using strain-controlled rheometers (coaxial cylinders or parallel plate). The flow behaviour index *n* and the consistency coefficient *m* were obtained by fitting an Ostwald–de Waele power law to the measured data in the expected shear rate range (less than 5 s^{−1}). The relative uncertainties associated with *n* and *m* were ≈2.2% and ≈3.8%, respectively. The mass density was determined using a hydrometer or by direct weighing, with a relative uncertainty of ≈1%. The temperature was measured at regular intervals during each experiment by a mercury-in-glass thermometer (accuracy 0.02°C) and proved to be significantly constant during each test and uniform across tests.

Table 4 lists the main parameters of the tests performed in diverging channels. Figure 7 shows experimental results compared to theoretical predictions, respectively, for linearly diverging horizontal and inclined channels with a semicircular cross section. The scaled dimensionless distance is depicted as a function of dimensionless time, expressing all factors in the various expressions of *X*_{N} as *f*(*X*_{N}). Deviations from the theoretical predictions occur at early times; at intermediate and late times, the agreement with theory is very good for both constant volume (*α*=0) and constant volume flux (*α*=1). In the former case, there are several reasons that justify the discrepancy between theory and experiments. First of all, the fluid is initially in a lock of finite length, whereas the theory assumes that the volume release is described by a Dirac function at *t*=0. In addition, the lock is not located in the virtual origin but has its initial section at *x*=50 cm; moreover, the gate of the lock is removed in a finite time and possibly with an imperfect movement, which perturbs the initial phase of the motion. Similar problems were faced by other researchers dealing with experiments modelling an instantaneous fluid release. However, they found that the imperfect boundary and initial conditions were soon forgotten by the current [8,25], that after a finite time evolves into a self-similar behaviour representing an intermediate asymptotic in the sense of Barenblatt & Zel’dovich [44]. A preliminary indication of the time shift is given by the theoretical model, but a more correct evaluation is here provided using a best-fitting procedure with data corresponding to the intermediate-time propagation of the current. No instability was observed during the experiments.

The model used to interpret the experimental data is represented by the general equation *h* at a given location. The problem parameters are affected by measurement or estimate errors; therefore, the propagation of uncertainty to the output variables through the variance formula allows deriving the associated confidence limits under the assumption that the input variables are independent. Figure 8 depicts the current length against time for test 13 compared with the theoretical results and the 95% confidence limits, evaluated by assuming a standard deviation for each parameter on the basis of its accuracy (see figure caption). The individual contributions of the input parameters to the overall standard deviation of the current length are shown in figure 8 against time. For the present set of parameters, the uncertainty of the consistency coefficient *q*. At early times, the uncertainty of the fluid index *n* also accounts for a significant fraction of the total uncertainty. As noted earlier, for two tests with constant flow rates (11 and 13), the thickness of the current was measured in a single cross section; results are shown in figure 9. The left panel refers to a test in a horizontal channel, which shows a satisfactory agreement with the theory even though the two curves exhibit a time lag and a slightly different shape. The right panel refers to the case of an inclined channel. The arrival time of the front end of the current and the maximum height are correctly predicted, the thickness of the current at late times is within the confidence interval. On a final note, some preliminary experiments were also conducted in inclined converging channels with circular cross-sectional shapes, obtained by inverting the orientation of the experimental channel; the geometric parameters thus were *k*=2 and *b*=−1. These tests yielded a current speed significantly slower, by a factor of 2, than that predicted by the theoretical model in this case. We tentatively interpret these discrepancies as a demonstration that the effects of the normal stress exerted on the fluid by the walls, and acting against the flow, are significant and cannot be neglected. Further experiments are needed to support this conjecture. For horizontal converging channel, a self-similarity solution of second kind could arise [38,39]. This analysis is left for future work.

## 5. Conclusion

We have theoretically and experimentally investigated the flow of laminar gravity currents of power-law fluids in horizontal or inclined channels or fractures, with different cross-sectional shapes and longitudinally variable width. The theoretical scalings allow for the evaluation of the position of the current front and the thickness of its profile.

Laboratory experiments to validate the theory were conducted with Newtonian, shear-thinning or shear-thickening fluids in a linearly divergent channel of circular cross section. The results for the current front show a fair agreement with theoretical predictions, well within the 95% confidence limits; the measurements of the current height variation over time in a given cross section are in satisfactory agreement with the model. However, since a single geometry was tested, no general conclusion can be drawn on the performance of the theory, especially for the converging channels and squeezing fractures. An uncertainty analysis on the propagation rate of the current reveals a stronger dependence on the variability associated with *α* and *k* than with *n* and *b*.

The behaviour of the current during its spreading is encapsulated in three key factors, *c*=*F*_{1}, *F*_{2}, *F*_{2}−*F*_{1}, which govern the rate of change over time of the current speed, thickness and aspect ratio. These factors depend on the model parameters *k*, the transversal shape of the channel (*k*>1) or fracture (*k*<1); *b*, the rate at which a channel widens and a fracture narrows along the flow direction; *n*, the flow behaviour index describing rheologically non-Newtonian power-law fluids; and *α*, the rate of growth of the current volume over time.

The dependency of these factors on model parameters is complex and governed by the influence of mass balance on one hand and of viscous resistances on the other. Several critical values of *α* and *b* are identified and determine the tendency of the current to: (i) accelerate/decelerate, (ii) become thicker/thinner and (iii) increase or decrease its speed, thickness or aspect ratio with increasing *k*, *b* and *n*. Some of these critical factors act as thresholds to limit the asymptotic validity of the resultant solutions to within certain ranges of parameters.

In particular, the dependence on influx rate (*α*) and longitudinal geometric variations (*b*) is stronger than the dependence on cross-sectional shape (*k*) and fluid rheology (*n*). The triangular shape is associated with the fastest propagation rate for (i) strongly divergent channels or squeezing fractures (*b*>*b*_{k}) irrespective of influx rate or fluid rheology and (ii) moderately divergent channels or squeezing fractures (*b*<*b*_{k}) coupled with a moderate influx rate (*α*<*α*_{k}). For (*b*>*b*_{k}) and *α*>*α*_{k}, the triangular section has the lowest propagation rate. Inclined channels exhibit faster currents.

Our model can be helpful for interpreting a variety of environmental or industrial channelized viscous flows of fluids with negligible or modest yield stress. Many of these, e.g. lava flows, involve pseudoplastic or Newtonian fluids (*n*≤1) flowing in uniform or mildly diverging channels (*k*>1, *b* of the order of a few per cent), with an influx rate constant or decreasing with time (*α*≃1 or lower), except possibly in the initial stages of flows. In these cases, *α*_{k} is larger than unity, while *b*_{k} is of the order of *α*≃1) in fractures which typically squeeze away from the source (*k*<1, *b*>0). In this case, *b*_{k} is again of order unity; hence, depending on the rate of squeezing, the current speed in nearly triangular fractures is either the fastest or the slowest.

## Data accessibility

The datasets supporting this article have been uploaded as part of the electronic supplementary material.

## Funding statement

Support from Università di Bologna RFO (Ricerca Fondamentale Orientata) 2011 and 2012 and from the Università di Parma FIL 2011 is gratefully acknowledged.

## Author’s contributions

L.C. and S.L. developed the experimental apparatus and techniques. V.D.-F. and S.L. gave an equal contribution to all other aspects of the work. All authors gave final approval for publication.

## Conflict of interests

We have no competing interests.

## Acknowledgements

The authors thank the anonymous referees for their comments on this manuscript.

- Received February 2, 2015.
- Accepted April 13, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.