Stability Analysis of Gravity Currents of a Power-Law Fluid in a Porous Medium

We analyse the linear stability of self-similar shallow, two-dimensional and axisymmetric gravity currents of a viscous power-law non-Newtonian fluid in a porous medium. The flow domain is initially saturated by a fluid lighter than the intruding fluid, whose volume varies with time as t. The transition between decelerated and accelerated currents occurs at α = 2 for two-dimensional and at α = 3 for axisymmetric geometry. Stability is investigated analytically for special values of α and numerically in the remaining cases; axisymmetric currents are analysed only for radially varying perturbations.The two-dimensional currents are linearly stable for α < 2 (decelerated currents) with a continuum spectrum of eigenvalues and unstable for α = 2, with a growth rate proportional to the square of the fluid behavior index. The axisymmetric currents are linearly stable for any α < 3 (decelerated currents) with a continuum spectrum of eigenvalues, while for α = 3 no firm conclusion can be drawn. For α > 2 (two-dimensional accelerated currents) and α > 3 (axisymmetric accelerated currents) the linear stability analysis is of limited value since the hypotheses of the model will be violated.


Introduction
A thorough comprehension of numerous class of flows in porous media is firmly based on approximations and simplifying assumptions leading to differential problems amenable to analytical or quasi-analytical solutions. A notable example is the shallow water assumption, with negligible vertical velocity and a horizontal length scale much larger than the vertical one [1][2][3][4]; the approximation leads in turn to modelling a wide class of propagation problems as onedimensional ones. An important class of solutions arising in several one-dimensional transient flows driven by gravity is constituted by self-similar solutions, representable in terms of functions of one variable and describing the "intermediate asymptotic" behavior of the system in a region where the solution is no longer dependent on the specific initial and/or boundary conditions adopted [3]. Whenever self-similar solutions can be obtained in analytical or in approximate analytical form, they represent an excellent test to verify model correctness by means of comparison with experiments; they also constitute a reliable benchmark for numerical integration.
Noteworthy examples of self-similar solutions to singlephase gravity driven flows in porous media include flow over a horizontal surface in plane [2] and radial geometry, Lyle et al. [5]; these solutions were extended to incorporate twolayer flow [6], the effect of a sloping bottom [7], the action of impermeable confining boundaries [8], drainage effects [9], and non-Newtonian power-law fluids [10], also in a vertically graded porous medium [11]. The dipole solution, originally derived by Barenblatt [3], was extended by King and Woods [12] to include drainage, and by Mathunjwa and Hogg [13] to incorporate a vertical variation in permeability.
An important step following the derivation of a self-similar solution is the check of its stability. In fluid mechanics, such an analysis can essentially be reduced to the study of the spatiotemporal evolution of disturbances applied to the system. In a broad sense, the evolution of instabilities leads to new solutions of the differential problem describing the flow and provides hints on possible stable solutions.

Mathematical Problems in Engineering
Stability analyses were performed in the context of several filtration problems amenable to self-similar solutions. The linear stability of the Barenblatt-Pattle (B-P) self-similar solutions of the porous medium equation, describing a number of flows including viscous and porous media gravity currents, was investigated by Grundy and McLaughlin [14] for plane and axisymmetric geometry. Their analysis was later extended to nonradially symmetric perturbations by Mathunjwa and Hogg [15]. The linear stability of a class of self-similar solutions to the filtration-absorption equation was demonstrated by Barenblatt et al. [16] and Chertok [17]. Mathunjwa and Hogg [13] showed the linear stability of the dipole self-similar solution of the first kind [3]. With respect to all these analyses the present study is new since it considers non-Newtonian fluids.
The study of gravity currents in porous media was recently extended to non-Newtonian fluids by Di Federico et al. [18,19], respectively, for two-dimensional and axisymmetric geometry; self-similar solutions in analytical and numerical form were obtained for a power-law fluid spreading in a uniform medium subject to the injection of a volume of intruding fluid increasing with time with the exponent . The motivation for their study was the nonlinear rheological nature displayed by fluids flowing in porous media in a variety of industrial and environmental applications ( [20] and references therein).
In this paper, we analyse the stability of the aforementioned solutions. Section 2 presents a linear stability analysis for two-dimensional geometry. First, the special cases = 0 (correspondent to an instantaneous release of a finite volume of fluid) and = 2 (correspondent to a linear increase in time of injected volume) are discussed, deriving a closedform expression for the eigenvalues of the associated Sturm-Liouville (SL) problem. Second, the stability of the general case < 2 with ̸ = 0 is demonstrated numerically. In Section 3, a similar analysis is conducted for axisymmetric geometry, considering radially varying perturbations. Stability is investigated analytically for = 0 and numerically for 0 < < 3. Concluding remarks are outlined in Section 4. Further mathematical details are included in the Appendices.

Stability Analysis in Two-Dimensional Geometry
Di Federico et al. [18] (Figure 1) consider the motion of shallow two-dimensional gravity currents of a purely viscous and relatively heavy non-Newtonian fluid of uniform density in a homogeneous porous layer saturated with a lighter fluid of uniform density −Δ . The flow is driven by the release from a source at the boundary = 0 of a volume ∝ , above a horizontal impermeable bottom. The non-Newtonian fluid in simple shear is described by the rheological Ostwald-DeWaele power-law model: where is the shear stress,̇the shear rate, the fluid consistency index, and the flow behavior index. When 0 < < 1 the model describes a shear-thinning behaviour, when  The flow law for the fluid is a modified Darcy's law taking into account the nonlinearity of the rheological equation (1) proposed by Bird et al. [21]: where = + is the generalized pressure, the pressure, the vertical coordinate, the fluid density, the specific gravity, u the Darcy velocity, the intrinsic permeability coefficient, and eff the effective viscosity. The mobility ratio / ef is [22] where denotes the porosity. It is convenient to write (2)-(3) as [23] highlighting the dependence upon permeability . The Darcy velocity in the direction ( , ) is given by the one-dimensional counterpart of (4), which reads, for horizontal flow In the shallow water approximation, the pressure gradient driving the flow is thus related to the gradient of the unknown free surface of height ℎ( , ) by Mathematical Problems in Engineering 3 where any capillary forces (surface tension) between the two liquids have been neglected. Introducing (6) into (5), we find that For one-dimensional transient flow, the local continuity condition takes the form Substituting (7) into (8), we obtain the nonlinear partial differential equation for height ℎ( , ) as: Consider the flux to be such that the total volume per unit width of the dense fluid in the porous medium at any time is given by , where both and are constant. The global mass balance equation then reads The mathematical statement of the problem is completed by the boundary condition: Equations (9)-(11) may be nondimensionalized by setting where the time, space, and velocity scales are * = ( 2 ) thereby generalizing to non-Newtonian case the choice of Vella and Huppert [7], valid for ̸ = 2 (stability of the special case = 2 is discussed later in Section 2.2). Equations (9)-(11) may be written in nondimensional form as where and are dimensionless space and time, is the dimensionless current depth, and ( ) represents the length of the gravity current. By introducing the similarity variable = −( + )/(2+ ) , and denoting the value of for = ( ) by , the similarity solution takes the form with Φ( ) being the solution of subject to Φ(1) = 0, while is given by Another necessary boundary condition is derived for the first derivative at the front end of the current. In the limit → 1 a Frobenius series expansion gives [18,19] and the solution is determined only by the front end conditions.
Once is determined, the gravity current length is given by ( ) = ( + )/(2+ ) ; the case = 2 acts as a transition between decelerated and accelerated currents. In the cases = 0 and = 2, the problem is amenable to an analytical solution. For a generic value of , approximate solutions are available in Frobenius series or in quasi-analytical form; else, a numerical integration of (16) is performed to derive the shape function Φ( ) [18].
To analyze the linear stability of the self-similar solution expressed in functional form by (15), we note that (14a) can be written as The first step is a variable transformation with = / ≡ /( ( + )/(2+ ) ) in order to avoid the difficulties in managing a movable boundary and to reduce the domain to −1 ≤ ≤ 1. In the new independent variable , (19) becomes 4

Mathematical Problems in Engineering
The perturbed solution can be expressed as where 0 is the basic state, is a small parameter, and ℎ is a perturbation. Due to perturbation, the boundary is displaced and we deal with bounded perturbed support −1 − 2 ( ) ≤ ≤ 1 + 1 ( ), where 1 ( ) and 2 ( ) are the perturbations of the boundaries at the edge of the current. To evaluate the structure of 1 ( ) and 2 ( ), we can use the analytical expression of the basic state available for = 0; that is, . Substituting in one of the boundary conditions (1 + 1 ( ), ) = 0 and linearizing results in , it follows that the displacement is of ( ) and can be expressed as assuming a generic nonsymmetric perturbation; in (22a) and (22b), + and − are the perturbed edges of the domain in the positive and negative directions, and 1 + and 1 − are the first coefficients of their series expansion. Substituting the perturbed solution (21) in (20) and comparing powers of , at (1) the equation for the basic state is recovered. At ( ), the following expression in the perturbed variable ℎ is derived: The basic state 0 is expressed by (15) and we seek a solution for ℎ in (23) by separating the variables according to the following expression: In (24a), (24b), and (24c), the function ( ) is the eigenfunction corresponding to the th eigenvalue , while , + , and − are the coefficients of the series. If the eigenvalues exist and are nonnegative, and the set of eigenfunctions is complete, the perturbation decays in time and the basic state is stable.
The boundary conditions at ( ) become in terms of ℎ After substitution of (24a), (24b), and (24c) into (23), the perturbation is the solution of the following differential problem, where the prime symbol indicates the derivative with respect to : with and the following boundary conditions: Equation (26) with boundary conditions (28) can be converted in the classical self-adjoint SL form [24]: with ( ) = ( ) ( ) , Mathematical Problems in Engineering 5 The SL structure of the differential problem guarantees that (i) the set of eigenfunctions is complete and (ii) different and orthogonal eigenfunctions correspond to different eigenvalues. In the following, we first evaluate the stability of the solution for the two special cases = 0 and = 2, for which a basic solution is available in closed form; we then discuss the stability of the general solution valid for ̸ = 0 and < 2. We do not explore the case > 2, which is of limited practical interest and subject to the limitations implicit in the shallow water approximation, namely, the condition / ≪ 1. The spatial gradient of the current height is readily obtained by (15) as showing that the former condition is time dependent and that the spatial gradient within the domain tends to decay with time if < 2, to remain constant if = 2, and to increase with time if > 2. In the latter case, the basic solution obtained by solving (16)- (17) is valid only for a limited time interval in a portion of the domain, but sooner or later it will not satisfy the shallow water approximation in the entire domain. Equation (31) can be reformulated introducing the average gradient of the shape function at a given time, equal to Imposing that ℎ/ = ≪ 1 (a small quantity) at time requires that < 2 + 3 2 ln ( /2) ln .
(33) Figure 2 shows that at short times the limiting condition is even more stringent than the asymptotic condition for → ∞.
2.1. The Case = 0. For = 0, the self-similar solution has the closed-form expression and the SL problem describing the perturbation can be handled with analytical functions. The differential problem for the perturbed variable takes the form (29) with and it is possible to reduce the differential problem to a standard hypergeometric one by using the transformation → 1/( +1) . The solution in terms of hypergeometric functions is then We further note that replacing by − the differential problem admits the solution Consequently, where the second coefficient 2 has been dropped to ensure the regularity in the origin. The spectrum is continuum with real and positive values [25] and is given by Hence, the eigenvalues are always positive, and the basic solution is stable. The spectrum is linear and the minimum eigenvalue is equal to 0 = 2 /( + 2). Perturbations in shear-thickening fluids gravity currents decay faster than perturbations in shear-thinning fluids gravity currents.
2.2. The Case = 2. If = 2 a different scaling is required since the characteristic time scale defined by (13) breaks down and an additional velocity scale arises, namely ( / ) 1/2 , in addition to the velocity scale [18]. Defining an arbitrary time scalẽ * and spatial length scales̃ * = ( / ) 1/2̃ * and dimensionless variables as̃= /̃ * ;̃= /̃ * ,̃= /̃ * , and̃= ℎ/̃ * , the original dimensional problem can be expressed by the following differential equation in dimensionless form: where is the ratio between the two velocity scales in the problem. By introducing the similarity variable = −1 , the similarity solution is of the form and the dimensionless length of the gravity current at a given time is given by ( ) = once is determined. Equation (40) has the following analytical solution [18]: In the linear stability analysis for = 2, the approximations are similar to those already applied for ̸ = 2. The details are reported in Appendix B. We look for a solution to (B.1) in the separable forms (24a), (24b), and (24c), with boundary conditions given by (25a), (25b), and (25c).
Substituting (43) into (B.1) gives with With the substitution = 1 + / , (45) can be reduced to a Laguerre form admitting a solution in terms of Laguerre functions with eigenvalues equal to = − 2 . The point = 1 is a regular singularity. For integer eigenvalues, the Laguerre functions become Laguerre polynomials, whose orthogonality requires nonnegative eigenvalues; hence, < 0, and the system is linearly unstable. The less unstable eigenvalue corresponds to shear-thinning fluids. Note that, ( ) = 0, in this case indicating that the value = 2 can represent the limiting value between two different behaviours of the sign of the eigenvalues: negative for < 2 and positive for > 2.
then there exist an infinite number of eigenvalues, which are all real and positive. The theorem can be extended to a singular SL problem in a finite domain with the singularity arising from the conditions (−1) = (1) = 0, provided that the existence of the eigenvalues and the completeness of the set of eigenfunctions are demonstrated. The singularity at the domain edges renders the boundary conditions redundant but requires the boundedness of ( ) at the edges. A sufficient, but not necessary, condition to guarantee that the set of eigenfunctions is complete is that the integral is finite (see, e.g., [27]), with being the associated Green function. For ̸ = 0 and < 2, a numerical approach is required to define the functions involved in the Sturm-Liouville problem. Starting from the numerical solution of the basic state, the functions , , and Υ are computed. and Υ are always positive and < 0 for < 2. The results are shown in Figure 3 for = 0.5, 1.0, 1.5, and = 1.5.
The boundedness condition of ( ) is satisfied, the integral in (47) is finite, and the system is linearly stable. This result is also supported by the numerical integration described in Appendix C. The results for = 1 (Newtonian fluid) are coincident with a similar analysis due to Grundy and McLaughlin [14].

The Axisymmetric Case
The axisymmetric case can be analysed by using the same approximations adopted in the two-dimensional case. The geometry is the one sketched in Figure 1 problem in nondimensional form is expressed for ̸ = 3 by the following equations [19]: with (± , ) = 0, representing the front position of the gravity current. Equations (48a) and (48b) do not include terms containing the derivatives with respect to azimuthal coordinate , since the present stability analysis is restricted to axisymmetric perturbations, modulated only along the radial coordinate. By means of the similarity variable = −( + )/(3+ ) , the solution to (48a) and (48b) was derived in the form The current is accelerated for > 3 and decelerated for < 3.
The analytical solution of the basic state is available only for = 0: For arbitrary (for the special case = 3, see [19]) (49)-(51a) can be solved numerically with (51b) and a second boundary condition, which can be computed developing the asymptotic solution in terms of a Frobenius series expansion near the current tip ( = 1) obtaining Since results = / ≡ / ≡ / ( + )/(3+ ) , −1 ≤ ≤ 1, in analogy to the plane case, (48a) and (48b) become Expressing the perturbed solution as ( , ) = 0 ( , ) + ℎ( , ) and substituting in (54), the following equation in the perturbed variable ℎ is obtained at ( ): Assuming that the basic state 0 is expressed by (49), we seek its solution by separating the variables with the following boundary conditions:

Mathematical Problems in Engineering
After substitution, the perturbation is the solution of the following differential problem: with which can be converted in the classical SL form [24]. For the special case = 0, the differential equation (57) for the perturbations becomes in SL form: with It is then possible to reduce the differential problem to a standard hypergeometric one by using the transformation → 1/( +1) . The solution is in terms of hypergeometric functions: while the spectrum is continuum with real and positive values and is given by Hence, the eigenvalues are always positive and the perturbations decay in time for any value of the flow behaviour index . As for the two-dimensional case, the spectrum is linear and the minimum eigenvalue is equal to 0 = 2 / ( + 3). Perturbations in shear-thickening fluids gravity currents decay faster than perturbations in shear-thinning fluids gravity currents. For = 3, a different scaling is requested since the time scale breaks down and the differential problem is similar to (48a) and (48b) except for a coefficient on the left hand side (see [19]), in analogy with (40) in the plane case.
Since an analytical solution is not available for ̸ = 0 and < 3 (decelerated currents) no analytical solution for the perturbation can be found, but on the same arguments used in Section 2.3 we found that the conditions for stability (46) are met. In fact the bounded nature of the term ( ) at the edges of the domain can be demonstrated and the completeness of the set of eigenfunctions is verifiable. Also a numerical integration conducted with the same methodology adopted for plane geometry (not shown) confirms that the system is stable for < 3. The results for = 1 (Newtonian fluid) are coincident with a similar analysis due to Grundy and McLaughlin [14] and to Mathunjwa and Hogg [15]. For = 3 the basic state solution is not available and no analytical solution for the perturbation can be found. In this particular case (since it marks the transition from decelerated to accelerated currents) the term ( ) is bounded at the edges of the domain but the integral in (47) is not finite. Since the finiteness of the integral is a sufficient, but not necessary, condition to guarantee that the set of eigenfunctions is complete, no firm conclusion can be drawn, even though in analogy with the two-dimensional case a linear instability could be inferred.
For > 3 the spatial gradient, given by grows in time; hence, the shallow water approximation will be eventually violated. With the same reasoning adopted for the two-dimensional currents, < 3 is the limit of validity of the present linear stability analysis. In addition, recasting (63) in terms of similarity variables indicates that the limits at short times are more stringent than the asymptotic limit for → ∞ (see Section 2 for a similar analysis in the twodimensional case).

Conclusion
A novel linear stability analysis of self-similar solutions to non-Newtonian power-law gravity currents intruding in a porous layer saturated by a lighter fluid has been performed for two-dimensional and axisymmetric flows; in the latter case, the analysis was limited to radial perturbations. Our investigation leads to the following key results.
(1) For two-dimensional flows, closed-form expressions of the stability equation can be obtained in the special cases = 0 (instantaneous injection of a finite volume of fluid) and = 2. For = 0 and any , the corresponding differential problem is solved in terms of hypergeometric functions with a continuous spectrum of positive eigenvalues; hence, the system is linearly stable for all values of . The spectrum of the eigenvalues is linear and the decay rate for perturbation in shear-thickening fluids gravity currents is larger than the decay rate for perturbations in shear-thinning fluids gravity currents. For = 2 the differential problem is solved in terms of Laguerre polynomials with a discrete spectrum of negative eigenvalues; hence, the system is linearly unstable and the eigenvalues are proportional to − 2 with a higher growth rate in gravity currents of shearthickening than in those of shear-thinning fluids. For ̸ = 0 and < 2 (decelerated currents), the existence of a complete set of eigenfunctions is demonstrated and the flow is linearly stable for both shearthinning and shear-thickening fluids, including the Newtonian case. For > 2 (accelerated currents) the present linear stability is of limited value since accelerated currents are characterised by a spatial gradient growing in time; hence, the shallow water approximation will be eventually violated.
(2) For axisymmetric flows, an analytical solution is available only for = 0. In this case the differential problem shows a continuum spectrum of real and positive eigenvalues. Hence, the flow is linearly stable. The spectrum of the eigenvalues is linear and decay rate for perturbation in shear-thickening fluids gravity currents is larger than the decay rate for perturbations in shear-thinning fluids gravity currents. For ̸ = 0 and < 3 (decelerated currents), the flow is always linearly stable. For = 3 no firm conclusion can be drawn, even though in analogy with the plane case we would expect the arising of instability. For axisymmetric accelerated currents ( > 3) as for the plane case, the present linear stability is of limited values since the limiting approximations of the shallow water approximation will be eventually violated.

A. Stability Condition for a Regular Sturm Liouville Problem
Given the Sturm-Liouville problem

C. Numerical Integration for the Plane Case
In order to check the validity of our conjecture, we perform the numerical integration of (20) at a given time starting from a self-similar solution plus a perturbation. The integration is based on a fourth-order, four-stage explicit Runge-Kutta scheme [28], with centred-space finite difference scheme for the spatial derivatives. The numerical calculations are carried out for = 1.9, 1.0 and = 0.5, 1.0, and 1.5, simulating the perturbation with periodic disturbances satisfying the integral boundary condition. Though the numerical solutions conserve volumes to a high degree of accuracy, due to numerical dissipation a volume loss of less than 0.01% is computed after a time interval Δ = 0.1, starting at 0 = 10.0. Figure 4 shows the results of the integration for = 0.5 and = 1.9 and 1.0, respectively, in (a) and (b). The perturbation is a sine function, given by ℎ ( , 0 ) = 0.00001 sin ( 2 0.2 ) . (C.1) Inspection of (a) and (b) reveals that the initial perturbation tends to decay with time. Similar results are obtained for = 1 and = 1.5 (not shown). Stability is thus confirmed.