The paper theoretically investigates the stability properties of the dam-break wave of a fluid with power-law rheology. Assuming the long-wave approximation, a depth-averaged flow model is considered. The linear stability analysis of the wave is carried out to individuate the marginal stability conditions. To this aim, the multiple-scale technique is applied with reference to the kinematic wave solution, which formally limits the validity of the theoretical achievements to relatively long time scales. Both shear-thinning and shear-thickening fluids are considered. Similarly to the case with uniform conditions, the analysis indicates that stable conditions can be associated with a marginal value of the Froude number. However, differently from the uniform conditions, the marginal Froude number is shown to be a function not only of the power-law index but also of the streamwise gradient of the base flow velocity and of the disturbance wavelength. The critical Froude number is found to be larger than the corresponding one in uniform conditions. Numerical solutions of the full model confirmed the outcomes of the linear stability analysis for both shear-thinning and shear-thickening fluids.

## I. INTRODUCTION

A fluid layer flowing on a sloping bed may exhibit free surface instabilities, which are often referred to as roll-waves. Their presence may cause overflows and intermittent forces when impacting against obstacles, increasing the risk of infrastructure damages.^{1} The presence of roll-waves has been observed in water but also in non-Newtonian fluids, such as mud, debris, and lava flows.^{2} The majority of literature works studied surface instabilities under parallel-flow conditions, even if they may develop also over non-parallel and transient base flows, such as in dam-break phenomenon, represented by a rapid release of a mass of fluid, which can be water or exhibit a non-Newtonian behavior.^{3} The analysis of the kinematics of dam-break still has a great interest.^{4}

This study is focused on the definition of the unstable conditions of a non-parallel time-varying base flow of a non-Newtonian fluid, represented by a dam-break wave generated by the rapid release of a mass of power-law fluid. The power-law model is suitable in describing shear-thickening (with power-law index *n* > 1) and shear-thinning (*n* < 1) behavior of different polymeric liquids, such as colloids and suspensions,^{5} but also of geophysical flows involving mud or lava.^{6,7}

The formation of roll-waves in non-Newtonian fluids has been widely studied within the framework of the long-wave approximation and using depth-integrated models. This approach offers the advantage that the dynamics of the flowing layer can be analyzed without requiring the detailed knowledge of the internal structure of the flow.^{8,9} In this regard and considering the power-law rheology, the instability development has been early studied performing the linear stability analysis of a parallel (or uniform) flow over an inclined plane.^{10} Berezin *et al.*^{11} investigated the effect of the inertial terms on the instability of a thin film of a fluid obeying to a power-law constitutive relation of the Ostwald–de Waele type through a linear stability analysis and a non-linear numerical solution. Lin and Hwang^{12} approximated the governing equations in the long-wave regime and applied the method of perturbation with multiple scales to derive the non-linear stability conditions of a power-law film. The results showed that the stability conditions are strongly affected by the power-law exponent *n* with the system becoming more unstable as *n* decreases. Dandapat and Mukhopadhyay,^{13,14} using the method of integral relations, derived an evolution equation revealing the presence of both kinematic and dynamic wave processes, which may either act together or singularly dominate, depending on the value of the problem parameters. Miladinova *et al.*,^{15} using the Benney long-wave model, studied the non-linear evolution of falling uniform films of a power-law fluid flowing over an inclined plane. The results showed that a saturation of non-linear interactions occurs, resulting in a finite-amplitude permanent wave. Amaouch *et al.*^{16} developed a method that combines the lubrication theory and the weighted residual approach to describe the non-linear behavior a thin film of a power-law fluid both far from and close to the instability threshold in the presence of a small-to-moderate Reynolds number. Ng and Mei^{10} and Longo,^{17} using the simplified von Kármán's momentum integral approach, studied finite-amplitude permanent roll-waves in a shallow layer of shear-thinning and shear-thickening fluids, respectively. The conditions for the existence of a periodic discontinuous solution, constituted by smooth profiles with a monotonically increasing flow depth between periodic shocks, have been individuated for both fluids. Zayko and Eglit^{18} investigated the effect of oblique perturbations that propagate with an arbitrary angle to the velocity of the undisturbed fluid flow moving down an incline plane, showing that, for a power-law fluid and under certain conditions, oblique perturbations can lead to flow instability. Chesnokov^{19} performed a non-linear stability analysis on a film of a power-law fluid down an inclined plane using both two-dimensional governing equations and a depth-averaged hyperbolic simplification. The comparison of the results of the two models indicates that the amplitude of roll-waves obtained by the 2D equations is slightly larger than that for the 1D model and for certain flow parameters the small perturbations of the basic solution grow for the 2D equations and decay for the depth-averaged model. Other studies have addressed the stability conditions of power-law films including the effects of wind stress,^{20} bed porosity,^{21–25} electric and magnetic fields,^{26,27} thermocapillarity,^{28,29} and oscillation of the inclined plane.^{30}

All of the above cited works address instabilities developing over an unperturbed state represented by the uniform flow of a power-law fluid. The stability problem of the steady solutions of the flow equations different from the uniform one has been investigated for several gradually varying profiles^{31} and in the presence of a wavy impermeable^{32} and porous^{33} bottom wall.

Focusing the attention on transient flow conditions, roll-waves have been observed in clear-water flows such as in the dam-break experiments carried out by Logan and Iverson.^{34} It has been found that instabilities develop after some time, i.e., at some distance from the dam where the bulk of the flow reaches a quasi-uniform and quasi-steady state. Motivated also by the above-mentioned experimental evidence, Bohorquez^{35} carried out a linear stability analysis to predict the occurrence of surface instabilities superposed to a dam-break water wave. A multiple-scale technique^{36,37} was applied to the 1D shallow-water model considering a space and time varying base flow. The study assumes the kinematic approximation to be valid. The validity of such an assumption, although restricted to relatively long time scales, has been shown by the author himself and confirmed by several further studies.^{38–41} Bohorquez^{35} concludes that the non-parallel and time-dependent characteristics of the wave play a fundamental role in causing the growth/decay of surface disturbances, as the result of the competition between kinematic and dynamic waves. Numerical simulations carried out using the fully dynamic 1D model confirmed the theoretical achievements.

Similarly to the clear-water case, surface instabilities have been experimentally observed even in association with dam-break waves of non-Newtonian fluids, e.g., Cochard and Ancey^{42} and Iverson *et al.*^{43} However, to the authors' knowledge, their description in terms of an instability mechanism has never been attempted.

The present paper aims to investigate the occurrence of unstable conditions in dam-break waves of power-law fluid according to the wave competition dynamics. To this aim, the depth-averaged 1D model of Ng and Mei^{10} is considered. Similarly to Bohorquez,^{35} the multiple-scale technique has been applied assuming as base flow the wave predicted by the kinematic model, widely used also with non-Newtonian fluids.^{3,44,45} The applicability of the kinematic model is first investigated. Then, the linear stability analysis is performed, assuming as base flow the space/time varying solution of the kinematic model. To verify the theoretical achievements, the predictions of the linear stability analysis are finally compared with the outcomes of numerical simulations of the full model. The study furnishes different original outcomes. First, the multi-scale technique is applied for the first time to a non-Newtonian fluid to investigate the development of hydrodynamic instabilities in dam-break waves, which constitutes a non-modal stability problem of great interest since the non-parallel time-varying characteristics of the base flow may affect the stability conditions. Moreover, the results permit to pinpoint the intrinsic differences between water and non-Newtonian power-law rheology behavior of free surface small disturbances. Finally, it furnishes indications on the development of roll-waves in geophysical flows such as those involving mud or lava, which is an aspect of practical and technical interest for the assessment of the correlated risk and for the design of appropriate countermeasures.^{46–48}

The paper is organized as it follows. Section II presents the formulation of the problem, while in Sec. III the applicability of the kinematic approximation is discussed. In Sec. IV, the linear stability analysis is performed and the results are illustrated. In Sec. V, theoretical predictions of the linear stability analysis are compared with the numerical results of the fully non-linear model. Finally, conclusions are given in Sec. VI. The full expression of the key parameters required for the evaluation of the stability conditions of the flow is given in the Appendix.

## II. PROBLEM FORMULATION

*θ*, with $ tan \u2009 \theta \u2264 O ( 1 )$] respect to the horizontal plane, without lateral inflows or outflows. Neglecting the surface tension and considering a laminar gradually varied flow where spatial variations occur over length scales larger than the flow depth, the dimensional depth-averaged mass and momentum conservation equations are as follows:

^{10}

^{,}

*g*and

*ρ*are the gravity and the fluid density, respectively. The expressions of the momentum correction coefficient

*β*and of the bottom shear stress $ \tau \u0303 b$ in laminar conditions read

^{10}

*μ*and

_{n}*n*represent the consistency and the rheological index of the power-law fluid, respectively. For shear-thinning fluids the rheological index is smaller than one, while values larger than one represent shear-thickening fluids.

*B*being the basal drag coefficient:

_{d}## III. KINEMATIC APPROXIMATION

### A. Kinematic model equation

^{8}Let us denote with $ l \u0303 b$ the length scale of the flood wave (base flow) in the

*x*direction and with $ \epsilon b = \eta \u0303 0 / l \u0303 b$ the corresponding shallowness parameter. To describe the flow at sufficiently far distance from the dam, i.e., for $ \epsilon b \u226a 1$, it is convenient to introduce the following slower space and time variables:

^{49}

^{3,44,50}can be viewed as the solution at the zeroth order of $ \epsilon b$ of the problem (13) and (14)

^{35}To be consistent with volume conservation, the shallowness parameter $ \epsilon b$ is assumed to be

^{51}

^{,}

^{52}

^{,}

^{52}

^{,}

*et al.*

^{52}

The conditions for the applicability of the kinematic approximation are addressed in Sec. III B.

### B. Kinematic model applicability

^{10,31}

Independently of the *n* value, Eq. (27) reveals that the kinematic approximation loses its validity for small values of the Froude number, i.e., $ \varphi * \u2192 0$ when $ F \u2192 0$. Conversely, the limiting value $ \varphi *$ diverges when the Froude number equals the $ F *$ value and reduces when *F* increases. Figure 2 shows, for three different values of *n*, namely, $ n = 0.5 , 1.0 , 1.5$, the limiting value $ \varphi *$ as a function of the Froude number. Considering that the kinematic wave model may be applied whenever $ \varphi \u226a \varphi *$, i.e., in the region below the curves, Fig. 2 indicates that the increase in the power-law index favors the applicability of the kinematic approximation, consistently with the theoretical findings of Hogg and Pritchard.^{53} Indeed, Hogg and Pritchard^{53} have shown that, for sufficiently high values of *n*, i.e., *n* > 0.26, a reduction of the power-law exponent implies a bottom drag increase, which induces a reduction of the dam-break wave front celerity and therefore of the length of the flood wave $ l \u0303 b$.

*B*, Eq. (10), the following non-linear equation in $ \eta \u0303 0$ is obtained:

_{d}*θ*and the released fluid volume $ A \u0303$, the value of the reference length scale $ \eta \u0303 0$ is evaluated by solving Eq. (33), and therefore, the values of

*B*, $ \epsilon b$, and $ t \u0302 c$ [see Eqs. (10), (17), and (32), respectively] are known. Finally, the value of the minimum time for the kinematic model applicability in fast coordinate follows from $ t c = t \u0302 c / \epsilon b$.

_{d}## IV. LINEAR STABILITY ANALYSIS

This section presents the linear stability analysis of the dam-break wave of a power-law fluid, by applying the multi-scale technique This technique has been fruitfully employed to analyze the stability conditions of two-dimensional laminar boundary layer^{54–56} and recently of water dam-break wave in turbulent conditions.^{35}

We want to analyze the stability properties of a base flow slowly varying in the space and in the time, i.e., $ h b ( x \u0302 , t \u0302 )$ and $ u b ( x \u0302 , t \u0302 )$. In addition to the characteristic length scales $ \eta \u0303 0$ and $ l \u0303 b$ (with $ \eta \u0303 0 \u226a l \u0303 b$), the perturbation length $ l \u0303 p$ has to be introduced as a third length scale. We further assume that not only the characteristic length scale of the base flow, $ l \u0303 b$, but also the perturbation length scale, $ l \u0303 p$, is much larger than $ \eta \u0303 0$.^{57}

^{35}we expand the flow depth $ h ( x , t )$ and the velocity $ u ( x , t )$ with respect to the perturbation parameter $ \epsilon p = l \u0303 p / l \u0303 b \u226a 1$ in the form:

^{54}and Citro and Luchini

^{57,58}in the derivation of the stability bounds for the steady and unsteady boundary layer, respectively, and to Bohorquez

^{35}with reference to the water dam-break problem.

*t*larger than

*t*. Consequently, we assume

_{c}**A**,

**B**, and

**C**matrices and of the source term

**f**, all of which depend on $ x \u0302$ and $ t \u0302$ only, in (38) are as follows:

**A**,

**B**,

**C**, and

**f**on the (

*x*,

*t*) pair, the Laplace transform is applied to Eq. (38) for

*x*> 0, so that the following o.d.e. is deduced:

**u**in $ x \u2208 [ 0 , + \u221e )$ and

**D**matrix can be diagonalized and let us denote with $\Lambda $ and

**V**the corresponding eigenvalues and eigenvectors matrices, respectively, i.e., $ \u2212 D = V \Lambda V \u2212 1$. Under these assumptions, it is easy to verify that Eq. (43) admits the following solution:

**u**at

*t*= 0 and

**I**the identity matrix.

^{35,59}Next, with the aim of performing the temporal stability analysis of the flow model (38), the complex number $ s = s r + i \u2009 s i$ is assumed to have null real part, i.e.,

*s*= 0, and an imaginary part $ s i = 2 \pi / \Lambda 0$, with $ \Lambda 0$ the disturbance wavelength. Moreover, it is convenient to introduce the local wavenumber

_{r}*α*defined as follows:

**D**admits two distinct eigenvalues $ \lambda \u0302 \xb1 = 2 n ( 3 n + 2 ) ( n + 1 ) u kin \u2009 sin \u2009 \theta \lambda \xb1$, with

Owing to the negative definiteness of $ \u211c ( \lambda \u2212 ) , \u2009 \lambda \u0302 \u2212$ does not play any role in defining the stability conditions; i.e., only a temporal decay of the perturbation is associated with the $ \lambda \u0302 \u2212$ eigenvalue.

Conversely, $ \u211c ( \lambda + )$ may be positive or negative. Therefore, either a growth or a decay of the perturbation may take place. Setting $ \u211c ( \lambda + ) = 0$ in (50) and solving for *F* provides the expression of the threshold stability Froude number, *F _{m}* (marginal Froude number), as a function of the rheological exponent

*n*, the parameter $\varphi $, and the disturbance wavenumber

*α*. The lengthy expression of

*F*is given in the Appendix, along with a more concise approximation, accurate for moderate values of the $\varphi $ parameter.

_{m}Figure 3 reports in the ( $ F , \alpha $) plane the contour lines of $ \u211c ( \lambda + )$ for *n* = 0.5 (a), *n* = 1.0 (b), *n* = 1.5 and (c), assuming $ \varphi = 0$. Independently of the *n* values, Fig. 3 shows that the neutral stability curves, i.e., the ( $ F , \alpha $) pairs such that $ \u211c ( \lambda + ) = 0$, correspond to a vertical line in the ( $ F , \alpha $) plane. Therefore, the marginal stability condition of a parallel flow is independent of the disturbance wavenumber and the critical Froude number *F _{c}*, defined as the minimum value of

*F*, is a function of the power-law index only. Moreover, for

_{m}*F*close to

*F*the growth rate is so small that a remarkably long channel would be needed for permitting the development of unstable waves up to a perceivable threshold. Finally, Fig. 3 suggests that

_{c}*F*coincides with $ F *$ (red lines), in agreement with the results of normal mode

_{c}^{10}and near-front expansion

^{31}analyses.

Figure 4, which is the counterpart of Fig. 3, refers to the $ \varphi = 10 \u2212 3$ case. Figure 4 puts in evidence that the presence of a velocity gradient strongly modifies the linear stability conditions. Independently of the *n* value, the neutral stability condition in the ( $ F , \alpha $) plane does not show up as a straight vertical line. Indeed, *F _{m}* is found to be a function, not only of the

*n*value as in the case of parallel flow, but also of the disturbance wave number

*α*. At very high values of

*α*,

*F*tends to a critical value

_{m}*F*, which represents the lower bound in terms of the Froude number for the instability. A reduction of

_{c}*α*leads to an increase in

*F*. Moreover, there exists a minimum value of

_{m}*α*below which the base flow is always stable, for all Froude number values (horizontal asymptote, $ \alpha c , \u221e$). A similar result has been found by Bohorquez

^{35}in analyzing the stability conditions of a dam-break water wave in turbulent regime. Finally, Fig. 4 shows that the increase in

*n*induces an increase in the $ \alpha c , \u221e$ value and then an enlargement of the stability region. Similarly to the uniform flow conditions, the stabilizing effect due to the increase in the power-law exponent may be mainly ascribed to the increase in the effective viscosity of the film flows with

*n.*

^{12}

With the aim of deeply analyzing the marginal stability conditions in the presence of a spatial velocity variation, Fig. 5 reports the marginal Froude number as function of *α*, for different values of $\varphi $ and *n*.

Independently of *n*, Fig. 5 shows that two branches exist. The former refers to very low values of the marginal Froude number and it is particularly evident at high $\varphi $ values, i.e., $ \varphi > 10 \u2212 2$. However, this branch does not belong to the region in which the kinematic model is applicable [see Eq. (27) and Fig. 2], and therefore, its interest is marginal. Conversely, the second branch develops always for $ F > F *$, and therefore, it refers to conditions in which the kinematic model can confidently approximate the full model.

Independently of the *n* values, Fig. 5 confirms the dependence of the marginal condition on $\varphi $ and the existence of the vertical (*F _{c}*) and horizontal ( $ \alpha c , \u221e$) asymptotes, both depending on the ( $ \varphi , n$) pair. The analytical expressions of both bounds

*F*and $ \alpha c , \u221e$, reported in the Appendix, have been found starting from Eq. (50), accounting for (51) and (52), setting $ \u211c ( \lambda + ) = 0$ and taking the limit for $ \alpha \u2192 + \u221e$ and $ F \u2192 + \u221e$, respectively.

_{c}Figure 6 highlights, for three different *n* values, the dependence of both *F _{c}* [Fig. 6(a)] and $ \alpha c , \u221e$ [Fig. 6(b)] on $\varphi $. Figure 6(a) shows that

*F*is always larger than $ F *$ and it tends to $ F *$ as $ \varphi \u2192 0$. Therefore, the unsteady non-parallel nature of the base flow leads to increase the stability region expressed in terms of the Froude number.

_{c}Moreover, also the minimum value of the wave number below which the flow is stable for all Froude number values increases with $\varphi $ starting from 0 [Fig. 6(b)]. It is also evident that the $ F c / F *$ ratio decreases as the rheological index *n* increases, while the opposite is observed for $ \alpha c , \u221e$.

Present results indicate that similarly to the turbulent clear-water flow,^{35} also for the power-law model, the temporal evolution of small disturbances superimposed on a kinematic wave strongly differs from the corresponding one in the presence of a plane-parallel flow. Indeed, there exists a minimum value of the disturbance wavenumber below which the disturbance does not grow; i.e., a cutoff wavenumber for the spectrum of the unstable perturbations occurs.

^{35}In fact, by substituting Eq. (19) in Eq. (48), it is easily verified that the Froude number monotonically increases along the flow direction, attaining its maximum value,

*F*, at the wave front, i.e.,: for $ x \u0302 = x \u0302 s$, which reads

_{s}Moreover, as stated before also $\varphi $ increases in the flow direction until the maximum value $ \varphi s$ expressed by Eq. (29) is reached at the downstream front. Since the increase in *F* and $\varphi $ has opposite effect in terms of stability, in the present case, the detection of unstable conditions requires, in general, the analysis of the local fulfillment of the stability criterion along the dam-break wave.

Under unstable conditions, it is expected that the instability predicted by the above analysis is observed even in the solution of the fully non-linear governing equations, provided that a suitable disturbance is applied to the unperturbed flow. Section V is devoted to verify this issue. The spatial/temporal evolution of small perturbations superimposed on a flood wave is analyzed by numerically solving Eqs. (8) and (9).

## V. NUMERICAL SIMULATIONS

As follows, the results of the stability analysis developed in the previous section are compared with the outcomes of numerical solution of the fully non-linear model. Indeed, the triggering of instability as the consequence of the suitable choice of the disturbance under linearly stable and unstable conditions is verified. To this aim, the propagation on an inclined plane with a constant slope *θ* of the wave arising from the sudden release of a mass of power-law fluid initially confined in a triangular wedge (Fig. 1) is simulated. The initial conditions in dimensionless variables are given by Eq. (11). The non-linear hyperbolic system (8) and (9) has been numerically solved with a finite volume scheme, which is second order in both time and space.^{31} The solution of the associated Riemann problem, i.e., the evaluation of the convective fluxes at the interface between adjacent finite volumes, is performed with the HLL approximation, whereas the temporal integration is performed with a Runge–Kutta scheme.^{60,61} Local absorbent boundary conditions have been imposed at the upstream and downstream limits of the channel.^{62} Additional details about the numerical method can be found in Campomaggiore *et al.*^{31}

Two power-law fluid rheologies are considered. The former describes a lava (herein denoted as *LF*) with a shear-thinning behavior characterized by the following parameters: *n* = 0.763, $ \rho = 2780 \u2009 kg / m 3 , \u2009 \mu n = 51.88 Pa s n$.^{63} The latter refers to a shear-thickening mud $ ( denoted as \u2009 M F )$ with the following rheology: *n* = 1.53, $ \rho = 1200 \u2009 kg / m 3 , \u2009 \mu n = 0.13 Pa s n$.^{64} In both fluids, the initial dimensional volume (per unit width) $ A \u0303$ is assumed equal to $ 10 \u2009 m 2$. The bottom inclination angles are as follows: $ \theta = 10 \xb0$ (*LF*) and $ \theta = 3 \xb0$ (*MF*). Therefore, the initial depths are $ h \u0303 0 = 1.88 \u2009 m$ and $ h \u0303 0 = 1.02 \u2009 m$, for the *LF* and *MF*, respectively.

The stability analysis developed in Sec. IV assumes that the base profile may be approximated (with a tolerance $\Psi $) by the solution of the kinematic model when the critical time *t _{c}* is overwhelmed. Indeed, preliminary tests have been carried out in order to numerically verify this conclusion, along with the theoretical estimate of

*t*given in Sec. III A.

_{c}As first step, for both fluids, the tolerance has been preliminary fixed using the following $\Psi $ values: 0.1 and 0.01, for the *LF* and *MF* tests, respectively. The solution of Eq. (33) therefore gives the value of the reference length scale $ \eta \u0303 0$ as $ \eta \u0303 0 = 0.15 \u2009 m$ (*LF*) and $ \eta \u0303 0 = 0.06 \u2009 m$ (*MF*). For both tests, Table I reports the computed values of the Basal drag coefficient *B _{d}*, of the shallowness parameter $ \epsilon b$ and of the minimum applicability time

*t*, deduced through Eqs. (10), (17), and (32), respectively.

_{c}Fluid . | B
. _{d} | $ \epsilon b$ . | t
. _{c} |
---|---|---|---|

LF | $ 1.56 \xd7 10 \u2212 1$ | $ 2.24 \xd7 10 \u2212 3$ | 160 |

MF | $ 4.23 \xd7 10 \u2212 2$ | $ 3.42 \xd7 10 \u2212 4$ | 1600 |

Fluid . | B
. _{d} | $ \epsilon b$ . | t
. _{c} |
---|---|---|---|

LF | $ 1.56 \xd7 10 \u2212 1$ | $ 2.24 \xd7 10 \u2212 3$ | 160 |

MF | $ 4.23 \xd7 10 \u2212 2$ | $ 3.42 \xd7 10 \u2212 4$ | 1600 |

All the simulations have been carried out keeping fixed the following values of both mesh spacing and time step: $ \Delta x = 0.25 , \Delta t = 5 \xd7 10 \u2212 4$. The results have been verified to be mesh-independent upon repetitions of the simulations with smaller values of both mesh spacing and time step.

Figure 7 compares the velocity profiles predicted by the fully dynamic (8)–(11) and the kinematic (18)–(20) models at different instants for both fluids. As theoretically expected, the agreement between the two solutions strongly improves as the time proceeds. Indeed, at *t* = *t _{c,}* the deviation between the slopes of the two velocity profiles is lower than 5%, in agreement with the results of Bohorquez.

^{35}Moreover, it has been verified that the difference on the peak value of the velocity is lower than 3% and 4% and the difference on the front position is lower than 2% and 3% for

*LF*and

*MF*, respectively. Therefore, present results confirm the validity the theoretical estimate of Sec. III A. Indeed, it may be concluded that in both tests the kinematic model can be confidently applied for $ t \u2265 t c$.

*t*=

*t*between the two models, simulations have been carried out assuming as initial condition the solution of the kinematic model at

_{c}*t*=

*t*,

_{c}*t*=

*t*), respectively. The values of $ x s ( t c )$, given by Eq. (20), are as follows: $ x s ( t c ) = 630$ (

_{c}*LF*) and $ x s ( t c ) = 4700$ (

*MF*).

For a given fluid, the linear stability analysis suggests that, at given Froude number *F* and the non-parallelism parameter $\varphi $ values, a small disturbance superimposed to the base flow may trigger or not the instability depending on the local wavenumber *α* and therefore on the disturbance wavelength $ \Lambda 0$.

*x*and

_{min}*x*, with

_{max}*x*> 0 and $ x max < x s ( t c )$ has been superposed to the longitudinal velocity profile at

_{min}*t*=

*t*. The imposed velocity disturbance is as follows:

_{c}*U*

_{0}the disturbance amplitude, assumed as $ U 0 = 5 \xd7 10 \u2212 3$ for both fluids.

At *t* = 0 and for all abscissas in the range $ x min \u2264 x \u2264 x max$, the values of *F* and $\varphi $ are known from Eqs. (48) and (22), respectively. Moreover, for a given value of $ \Lambda 0$, also the *α* values are known [Eq. (47)]. Indeed, a proper choice of the wavelength $ \Lambda 0$ allows to define perturbed initial conditions, i.e., the $ ( F , \varphi , \alpha )$ spatial distribution at *t* = 0, for which the linear theory predicts stable or unstable conditions.

For each fluid, two values of the perturbation wavelength have been considered, namely, $ \Lambda 0 , s$ and $ \Lambda 0 , u$. The former has been chosen to simulate linearly stable flow conditions, while the latter refers to unstable ones. The extents of the perturbed region, i.e., the values of *x _{min}* and

*x*, along with the wavelengths of disturbance $ \Lambda 0 , s$ and $ \Lambda 0 , u$ assumed in the numerical simulations are reported in Table II.

_{max}Fluid . | x
. _{min} | x
. _{max} | $ \Lambda 0 , s$ . | $ \Lambda 0 , u$ . |
---|---|---|---|---|

LF | 100 | 500 | 400 | 10 |

MF | 2000 | 4000 | 1000 | 100 |

Fluid . | x
. _{min} | x
. _{max} | $ \Lambda 0 , s$ . | $ \Lambda 0 , u$ . |
---|---|---|---|---|

LF | 100 | 500 | 400 | 10 |

MF | 2000 | 4000 | 1000 | 100 |

Figure 8 represents the spatial distribution of $ \u211c ( \lambda + )$ for the perturbed initial conditions obtained assuming the $ \Lambda 0$ values given in Table II. As shown in Fig. 8, $ \u211c ( \lambda + )$ corresponding to the perturbation $ \Lambda 0 , s$ is everywhere negative for both fluids, and therefore, the linear theory predicts stable conditions. Conversely, the spatial distribution of $ \u211c ( \lambda + )$, evaluated with reference to the $ \Lambda 0 , u$ values, shows the presence of $ \u211c ( \lambda + )$ positive values and therefore of unstable flow conditions. Thus, a perturbation growth (respectively decay) is expected in response to the disturbance characterized by wavelength $ \Lambda 0 , u$ (respectively $ \Lambda 0 , s$).

Figure 9 reports the flow depth [Fig. 9(a)] and the velocity profiles [Fig. 9(b)] for the test *LF* at three different times for $ \Lambda 0 , s$, confirming that the imposed perturbation does not grow, nor in time neither in space. Based on the numerical results for $ \Lambda 0 , u$, Fig. 10 depicts the flow depth [Fig. 10(a)] and the velocity profiles [Fig. 10(b)] at three different times. The results indicate that the instabilities become apparent at time *t* = 50, when small amplitude waves are observed on the free surface. Present results, therefore, confirm the predictions of the stability analysis for the lava fluid. Consistently with the linear analysis, the amplitude of the perturbation increases from *t* = 50 to *t* = 200. In the same time interval, the disturbances are observed to descend the dam-break wave up to reach the downstream front *x* = *x _{s}*. Moreover, a non-linear wave dynamics is clearly observed between

*t*= 50 and

*t*= 200, with a progressive coarsening of the perturbations and their evolution toward the form of a shock-wave train. This dynamics is typical of roll-waves superposed to uniform conditions.

^{65}

With reference to the numerical simulations of the mud fluid with $ \Lambda 0 , s$, Fig. 11, which reports the flow depth [Fig. 11(a)] and the velocity profiles [Fig. 11(b)] at three different instants, indicates the absence of disturbance growth. Similarly, Fig. 12 shows the flow depth [Fig. 12(a)] and the velocity profiles [Fig. 12(b)] at three different instants with reference to $ \Lambda 0 , u$. The growth of the disturbance on the free surface is confirmed, as also indicated from the velocity evolution [Fig. 12(b)]. The amplitude of the perturbation grows in time while descending the dam-break wave, as it is evident from the profiles at *t* = 1000, where the perturbation has reached the downstream front. Also, in this case, the fully non-linear analysis shows the progressive coarsening and steepening of the initially sinusoidal disturbance, peculiar of the roll-waves.^{65}

In conclusion for both shear-thinning and shear-thickening fluids, the predictions of the linear stability analysis are confirmed by the numerical simulations performed with the fully non-linear model.

## VI. CONCLUSIONS

The present paper has investigated the stability of dam-break waves of a non-Newtonian fluid propagating on an inclined plane. A depth-averaged shallow-water model with a power-law rheology has been adopted and both shear-thinning and shear-thickening fluids have been analyzed. The linear stability of non-uniform time-dependent conditions has been studied by applying the multiple-scale technique, widely used to investigate the stability conditions of two-dimensional laminar steady and unsteady boundary layer. The analysis has been performed assuming as initial condition the solution of the kinematic model, i.e., the outer solution of matched asymptotic expansion of the governing equations. Therefore, present results are formally restricted to relatively long time-scales, i.e., $ t > t c$ with *t _{c}* the critical time, for which the kinematic approximation is valid within a prescribed tolerance. It has been found that the critical time

*t*strongly depends on the power-law index, increasing as

_{c}*n*decreases. For both shear-thinning and shear-thickening fluids, the linear stability analysis has shown that the non-uniform time-dependent character stabilizes the kinematic wave and increases the marginally stable Froude number

*F*required for instability. Moreover, the marginal Froude number

_{m}*F*is not only a function of the power index

_{m}*n*, as in uniform conditions. In fact, in the kinematic wave case,

*F*depends also on the local wavenumber

_{m}*α*, which involves the disturbance wavelength, and on the measure of the spatial velocity variation $\varphi $. The analysis concludes that there exists a lower bound of the local wavenumber $ \alpha c , \u221e$, which depends on both

*n*and $\varphi $, below which the kinematic wave is stable for all Froude numbers. The $ \alpha c , \u221e$ value increases as the spatial velocity variation and the rheological index increase. Furthermore, the critical Froude number

*F*is larger than the value in uniform condition. The closed expressions of both

_{c}*F*and $ \alpha c , \u221e$ have been provided. Numerical simulations have been carried out solving the full model using a finite volume scheme, with second-order accuracy both in time and space. Both shear-thinning and shear-thickening fluids have been considered for analyzing the spatial/temporal evolution of small perturbations superposed to a kinematic flood wave in its applicability conditions. Numerical results have confirmed the theoretical achievements of the linear stability analysis.

_{c}Future research will be devoted to extend the range of validity of the assumed base flow, including higher-order terms in the solution.^{35} Furthermore, the definition of the convective/absolute character of the stability, along with the upset of finite-amplitude disturbances, will be worthy of exploration. As far as the former issue is concerned, the Briggs criterion could be applied,^{36} while for the latter the weakly non-linear stability analysis of the problem^{66,67} could be carried out. Moreover, a possible transient growth caused by the non-modal character of the spatial differential operators would deserve investigation, as well.^{68,69} Finally, additional efforts could be devoted to extend the presented study to a two-dimensional steep-slope flow.^{70}

## ACKNOWLEDGMENTS

The authors are indebted to Professor Patricio Bohorquez of Universidad de Jaén for fruitful discussions and suggestions. The role of the anonymous reviewers, which allowed to improve the quality of the paper, is also gratefully acknowledged.

The paper was supported by the SEND intra-university project, financed through the V:ALERE 2019 program of the University of Campania “Luigi Vanvitelli” (Grant No. B68D19001880005).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Cristana Di Cristo:** Data curation (equal); Validation (equal); Writing – original draft (equal). **Michele Iervolino:** Data curation (equal); Funding acquisition (equal); Validation (equal); Writing – original draft (equal). **Andrea Vacca:** Conceptualization (equal); Methodology (equal); Writing – original draft (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: MARGNAL BOUNDS EXPRESSIONS

*n*, of the non-parallelism parameter, $\varphi $, and of the spatial wavenumber of the perturbation

*α*. A relatively concise expression may be worked out as follows:

*F*involves up to the eighth power of the $\varphi $ parameter. However, it is easy to verify that in both the presented test-cases, the values of $\varphi $ are far below the unity [e.g., $ \varphi s ( t c ) = 1.2 \xd7 10 \u2212 2$ and $ \varphi s ( t c ) = 4.1 \xd7 10 \u2212 3$ for

_{m}*LF*and

*MF*, respectively]. Under such conditions, assuming the non-parallelism parameter $ \varphi \u226a 1$ and developing in series Eqs. (A2) and (A3) around $ \varphi = 0$, the following more handy form is obtained:

*F*,

_{c}Both Eqs. (A8) and (A12) may be employed to estimate the asymptotic boundaries for sufficiently small values of the $\varphi $ parameter.

Figure 13 compares the exact (lines) and the approximated (lines plus circles) values of $ F c / F *$ [Fig. 13(a)] and $ \alpha c , \u221e$ [Fig. 13(b)] for the same *n* values reported in Fig. 6. In the investigated range of *n* values, Fig. 13 shows that the provided approximations perform with reasonable accuracy (5% of the exact value) up to about $ \varphi = 0.05$ [Eq. (A8)] and up to approximately $ \varphi = 0.0125$ [Eq. (A12)].

## REFERENCES

*Stability and Transition in Shear Flows*

*Mineral Matter and Ash in Coal*