1 Introduction

Lord Rayleigh [1] and later Sir G.I. Taylor [2] both considered the unstable fluid flow that now bears their names, in which two horizontal immiscible fluids, with a sharp interface separating them, lie horizontally with the heavier fluid above the lighter one. Gravity acts downward, and any small disturbance to the interface results in the unstable growth of fingers of the heavy fluid moving downward and bubbles of the lighter fluid penetrating upward. Rayleigh and Taylor considered planar flow and carried out linearized analyses, based on the assumption that the interface shape never varies greatly from a horizontal plane. Their analysis predicts, however, that perturbations to the interface must grow exponentially in time, whenever the heavier fluid is above the lighter one.

Since linearized Rayleigh–Taylor theory predicts exponential growth of the interface, the linearizing assumption itself becomes invalid within finite time, and nonlinear effects must be included. The difficult question of how nonlinearity affects the shapes of the growing droops and bubbles at the interface has consequently been of great interest and relevance. Earlier numerical solutions for free-flowing inviscid fluids separated by a sharp interface observed that there was a finite time at which the solution abruptly failed, and the review by Sharp [3] mentions this as something of a mystery. However, Moore [4] had shown, in the related Kelvin–Helmholtz instability, that the sharp interface develops a curvature singularity within finite time, and this is the reason for the sudden failure of the inviscid solution at this critical time. Cowley et al. [5] confirmed Moore’s analysis with their own asymptotic theory, and Baker et al. [6] have shown that curvature singularities may also occur within finite time for the classical Rayleigh–Taylor problem.

It might be expected that the inclusion of viscosity in the fluids would avoid a curvature singularity from forming at the interface, but Forbes et al. [7] and Forbes and Bassom [8] have indicated that the situation is a little more nuanced than this. They showed that very large curvature spikes still form at the interface, even when viscosity is included, and to avoid them, a certain amount of fluid mixing is required, at an interface that must be of finite width. This explains why ‘vortex blob’ methods, introduced by Krasny [9] and used later by Baker and Pham [10], are able to continue an otherwise inviscid solution past the critical time at which curvature singularity would be expected, since they effectively replace the sharp interface with a diffuse vortex sheet of finite width.

The Rayleigh–Taylor problem has numerous applications, where it may help explain certain cloud formations in meteorology [11] and finger formation in astrophysics [12], and has thus been the subject of extensive research. A review by Kelley et al. [13] has emphasized the way in which Rayleigh–Taylor instability can occur over vast length scales, from the microscopic to the galactic. Furthermore, flows of this type are not just restricted to two-dimensional planar geometry. They have been studied in cylindrical geometry by Forbes [14], corresponding to outflow from a line source, and recently by Zhao et al. [15]. These flows also occur in spherical-type outflow from a point source, (see Forbes [16, 17]), where they have possible application to the formation of one-sided jet flows in astrophysics, observed by Gómez et al. [18]. The review article by Abarzhi [19] discusses Rayleigh–Taylor flows as mechanisms for fluid mixing, comparing and contrasting them with mixing due to turbulence.

Rayleigh–Taylor-type flows may also occur in porous media, perhaps in situations where leaching fluid is pumped into an ore-bearing rock, in an attempt to dissolve and extract the ore [20]. In such cases, the density difference between the two fluids present in the rock would occur as a result of the presence of dissolved minerals or leaching lixiviants. The presence of such minerals in each fluid of a two-fluid system was taken into account explicitly by Trevelyan et al. [21], in their study of fluid instabilities in porous media. They used a Boussinesq approach, combined with Darcy’s law relating the fluid seepage velocity in the porous medium, the gradient of the fluid pore pressure and the fluid density, and they invoked a type of equation of state in which the density varied linearly with the concentrations of two different chemical species dissolved in the fluid. These concentrations followed prescribed convection–diffusion equations. They carried out detailed quasi-steady stability analyses of various scenarios and also obtained numerical solutions to the nonlinear equations, obtaining long unstable fingers of heavier fluid moving into the region formerly occupied by the lighter fluid. In a recent review article, Hewitt [22] discusses unstable fluid flows in porous media, under conditions where convection strongly dominates diffusion. A state equation is considered, in which the fluid density may depend upon both the concentration of a dissolved chemical species and the temperature.

The present paper considers first the instability of a sharp interface between two horizontal layers of immiscible fluid in porous material, when the medium is perfectly saturated. This is analogous to the classical Rayleigh–Taylor problem, and we briefly discuss a linearized solution before considering the fully nonlinear saturated problem. A Boussinesq model is also discussed, and a numerical spectral method is introduced that allows us to generate accurate solutions for the evolution of very large amplitude disturbances to an interfacial zone between two miscible fluids. Results of both models are presented in Sect. 4, and some concluding remarks are given in Sect. 5.

2 Fully Saturated Flows

Fig. 1
figure 1

Schematic diagram for dimensionless two-layer immiscible fluid system in fully saturated medium

Full size image

In this section, we consider two unbounded fluids in a porous medium such that a sharp interface \(y = \eta \left( x,t \right)\) separates an upper heavier Fluid 1, in \(y > \eta\), from a lighter Fluid 2 in the lower region \(y < \eta\). A Cartesian coordinate system is present, with its y-axis pointing vertically as in Fig. 1, and the acceleration g of gravity is directed downward. Each fluid is presumed to have the same viscosity, but the pore pressure and fluid density in upper Fluid 1 are \(p_1\) and \(\rho _1,\) respectively, while in lower Fluid 2 the pressure is \(p_2\) and density is \(\rho _2\). Each density is constant throughout its entire fluid domain, and for Rayleigh–Taylor instability, we have \(\rho _1 / \rho _2 > 1\). According to Darcy’s Law [23, page 16], the fluid seepage velocities \(\mathbf{q}_1\) and \(\mathbf{q}_2\) in the two fluid regions are found from \(\mathbf{q}_1 = - K \nabla \varPhi _1\) and \(\mathbf{q}_2 = - K \nabla \varPhi _2\) , where \(\varPhi _1 = p_1 + \rho _1 gy\) and \(\varPhi _2 = p_2 + \rho _2 gy\) are the total hydraulic pressure heads in the two fluid domains. The constant K is an hydraulic conductivity and is sometimes expressed in the form \(K = A / \mu\), in which \(\mu\) is the dynamic viscosity of the fluid and A can be thought of as the cross-sectional area of a column of the porous medium through which the fluid percolates. For a fully saturated medium, the fluid incompressibility condition gives \(\nabla \cdot \mathbf{q}_j = 0\), \(j = 1,2\), with the result that the hydraulic heads \(\varPhi _j\) satisfy Laplace’s equation \(\nabla ^2 \varPhi _j = 0\), \(j = 1,2\), in their respective fluid domains. The dynamic condition that the fluid pore pressure must be continuous at the fluid interface gives rise to the simple condition \(\varPhi _1 = \varPhi _2 + \left( \rho _1 - \rho _2 \right) g \eta\) on the interface \(y = \eta \left( x,t \right)\).

Dimensionless variables are now introduced. If a periodic disturbance of fundamental wavelength \(\lambda\) is made to the interface, then \(\lambda / \left( 2 \pi \right)\) is taken to be the length scale from now on. Speeds are referred to the quantity \(\rho _2 g K\), and then time is made dimensionless using the scale \(\lambda / \left( 2\pi \rho _2 g K \right)\). Pressures and the hydraulic heads are scaled with respect to the quantity \(\rho _2 g \lambda / \left( 2\pi \right)\).

In these nondimensional quantities, Darcy’s law in lower fluid 2 becomes

$$\begin{aligned}&\mathbf{q}_2 \equiv u_2 \mathbf{i} + v_2 \mathbf{j} = - \nabla \varPhi _2 \nonumber \\&\varPhi _2 = p_2 + y , \quad \text {in } \quad y < \eta \left( x,t \right) , \end{aligned}$$
(1)

and in upper fluid 1 it takes the form

$$\begin{aligned}&\mathbf{q}_1 \equiv u_1 \mathbf{i} + v_1 \mathbf{j} = - \nabla \varPhi _1 \nonumber \\&\varPhi _1 = p_1 + Dy , \quad \text {in } \quad y > \eta \left( x,t \right) . \end{aligned}$$
(2)

The continuity condition for fully saturated rock gives Laplace’s equations

$$\begin{aligned} \nabla ^2 \varPhi _j = 0 , \quad j = 1,2 \end{aligned}$$
(3)

in each layer. There is only the one dimensionless constant \(D = \rho _1 / \rho _2\) appearing in (2). On each side of the interface, each fluid must obey its own kinematic condition

$$\begin{aligned} v_j = \frac{\partial \eta }{\partial t} + u_j \frac{\partial \eta }{\partial x} \quad j = 1,2 \quad \text {on } y = \eta \left( x,t \right) , \end{aligned}$$
(4)

representing the fact that neither fluid is able to cross this boundary. In addition, there is also a dynamic condition, corresponding to the fact that the pore pressures either side of the interface must be equal. In dimensionless form, this becomes

$$\begin{aligned} \varPhi _1 = \varPhi _2 + \left( D - 1 \right) \eta \quad \text {on } y = \eta , \end{aligned}$$
(5)

where the constant \(D = \rho _1 / \rho _2\) is the ratio of densities, as previously.

2.1 Linearized saturated solution

A linearized solution to the governing equations (1)–(5) can be obtained in a straightforward manner using classical methods. Accordingly, only a brief overview will be given here. If the interface \(\eta \left( x,t \right)\) is initially horizontal, but then subjected to a perturbation of small amplitude \(\epsilon\), it is appropriate to seek a solution in the form

$$\begin{aligned} \varPhi _j \left( x,y,t \right)= & {} \epsilon \varPhi _j^L \left( x,y,t \right) + {{\mathcal {O}}} \left( \epsilon ^2 \right) , \quad j = 1 , 2 \nonumber \\ \eta \left( x,t \right)= & {} \epsilon \eta ^L \left( x,t \right) + {{\mathcal {O}}} \left( \epsilon ^2 \right) . \end{aligned}$$
(6)

The true location \(y = \eta \left( x,t \right)\) of the interface is projected onto the undisturbed plane \(y = 0\), and as a result, conditions (3) are now approximated with

$$\begin{aligned} \nabla ^2 \varPhi _1^L = 0 \quad \text {in } y > 0 \nonumber \\ \nabla ^2 \varPhi _2^L = 0 \quad \text {in } y < 0 . \end{aligned}$$
(7)

The two kinematic conditions (4) at the interface become approximately

$$\begin{aligned} \frac{\partial \eta ^L}{\partial t} = - \frac{\partial \varPhi _1^L}{\partial y} = - \frac{\partial \varPhi _2^L}{\partial y} \quad \text {on } y = 0 \end{aligned}$$
(8)

and the dynamic condition retains its same form (5), but now holds on the approximate plane \(y = 0\) and involves the linearized functions introduced in (6).

For simplicity, we suppose that the initial perturbation to the interface was in the form of a simple cosine function, of period \(2\pi\). Then, the appropriate linearized solutions to the Laplace equations (7) are

$$\begin{aligned} \varPhi _1^L \left( x,y,t \right)= & {} Q_1 (t) \text {e}^{- y} \cos x , \quad \text {in } y > 0 \nonumber \\ \varPhi _2^L \left( x,y,t \right)= & {} Q_2 (t) \text {e}^y \cos x , \quad \text {in } y < 0 . \end{aligned}$$
(9)

Now the two kinematic conditions (8) and the corresponding linearized form of the dynamic condition (5) allow the unknown functions \(Q_1 (t)\) and \(Q_2 (t)\) in (9) to be determined, and a form for the linearized interface perturbation elevation \(\eta ^L \left( x,t \right)\) to be obtained. Combining this information with the expansion (6) therefore yields

$$\begin{aligned} \varPhi _1 \left( x,y,t \right)= & {} \frac{1}{2} \epsilon \left( D - 1 \right) \exp \left[ - y + \frac{1}{2} \left( D - 1 \right) t \right] \cos x \nonumber \\&\quad + {{\mathcal {O}}} \left( \epsilon ^2 \right) \quad \text {in } y > 0 , \nonumber \\ \varPhi _2 \left( x,y,t \right)= & {} - \frac{1}{2} \epsilon \left( D - 1 \right) \exp \left[ y + \frac{1}{2} \left( D - 1 \right) t \right] \cos x\nonumber \\&\quad + {{\mathcal {O}}} \left( \epsilon ^2 \right) \quad \text {in } y < 0 , \nonumber \\ \eta \left( x,t \right)= & {} \epsilon \exp \left[ \frac{1}{2} \left( D - 1 \right) t \right] \cos x + {{\mathcal {O}}} \left( \epsilon ^2 \right) \end{aligned}$$
(10)

for the final form of the linearized solution in the fully saturated material.

2.2 Nonlinear saturated solution

Equations (1)–(5) state the nonlinear mathematical problem to be solved, for the Rayleigh–Taylor instability in a completely saturated medium where a sharp interface is present between two immiscible fluids. Periodicity is again assumed in the x-coordinate, and so from (3), we seek nonlinear solutions

$$\begin{aligned} \varPhi _1 \left( x,y,t \right)= & {} \sum _{n=1}^N A_n (t) \text {e}^{- n y} \cos \left( nx \right) \quad \text {in } y > \eta \left( x,t \right) , \nonumber \\ \varPhi _2 \left( x,y,t \right)= & {} \sum _{n=1}^N B_n (t) \text {e}^{n y} \cos \left( nx \right) \quad \text {in } y < \eta \left( x,t \right) , \nonumber \\ \eta \left( x,t \right)= & {} H_0 (t) + \sum _{n=1}^N H_n (t) \cos \left( nx \right) . \end{aligned}$$
(11)

These expressions (11) become exact as the number N of Fourier modes becomes infinite, and in the numerical work to follow it will be necessary to take N as large as possible. Clearly, (11) satisfy Laplace’s equations (3) identically and represent the generalization of (10) to account for nonlinearity at the interface.

The three sets of coefficients \(A_n (t)\), \(B_n (t)\) and \(H_n (t)\) in the assumed form (11) of the solution are to be found from the kinematic and dynamic conditions at the interface. The spectral method proposed by Forbes et al. [24] is chosen for this purpose, and we use here only the simpler of their two approaches.

To begin, the kinematic condition (4) is taken with \(j = 1\), representing the upper fluid in \(y > \eta\). The fluid velocity components \(u_1\) and \(v_1\) in the x- and y-directions, respectively, are found directly from \(\varPhi _1\) in (11) by differentiation, according to (2). Thus,

$$\begin{aligned} u_1 \left( x,t \right)= & {} \sum _{n=1}^N n A_n (t) \text {e}^{- n y} \sin \left( nx \right) , \nonumber \\ v_1 \left( x,t \right)= & {} \sum _{n=1}^N n A_n (t) \text {e}^{- n y} \cos \left( nx \right) . \end{aligned}$$
(12)

When these two series (12) are substituted into the first kinematic condition (4) with \(j = 1\), the use of the product rule enables this equation to be cast in the form

$$\begin{aligned}&H_0^{\prime } (t) + \sum _{n=1}^N H_n^{\prime } (t) \cos \left( nx \right) = \sum _{n=1}^N A_n (t) \frac{\partial }{\partial x} \nonumber \\&\quad \left[ \text {e}^{- n \eta \left( x,t \right) } \sin \left( nx \right) \right] . \end{aligned}$$
(13)

Following [24], it follows at once from this relation that

Theorem 1

The average interface elevation does not change with time.

Proof

We wish to show that

$$\begin{aligned} \frac{\text {d} \eta _{\text {av}}}{\text {d} t} = 0 \quad \text {where} \quad \eta _{\text {av}} = \frac{1}{2\pi } \int _{- \pi }^{\pi } \eta \left( x,t \right) \, \text {d} x = H_0 (t) . \end{aligned}$$
(14)

To do this, it is only required to integrate (13) over a single period \(- \pi< x < \pi\), giving immediately the result \(H_0^{\prime } (t) = 0\). Thus, the average interface elevation \(\eta _{\text {av}} = H_0\) is constant, and the statement in (14) is proved. \(\square\)

Equation (13) is also Fourier analyzed with respect to the \(\ell\)th mode, \(\ell = 1 , 2 , \dots , N\), by multiplying by the basis function \(\cos \left( \ell x \right)\) and integrating over a single period. Integration by parts gives the expression

$$\begin{aligned} H_{\ell }^{\prime } (t) = \frac{\ell }{\pi } \sum _{n=1}^N S_{\ell ,n}^{(1)} (t) A_n (t) \end{aligned}$$
(15)

for the derivative of the \(\ell\)th Fourier coefficient in the expression for the interface shape.

These calculations may all be repeated for the second kinematic equation, obtained from (4) with \(j = 2\). Velocity components in lower fluid 2 are obtained by differentiation of the series for \(\varPhi _2\) in (11), similarly to (12), and are obtained in the forms

$$\begin{aligned} u_2 \left( x,t \right)= & {} \sum _{n=1}^N n B_n (t) \text {e}^{n y} \sin \left( nx \right) , \nonumber \\ v_2 \left( x,t \right)= & {} - \sum _{n=1}^N n B_n (t) \text {e}^{n y} \cos \left( nx \right) . \end{aligned}$$
(16)

These velocity components (16) are substituted into the second kinematic condition at the interface, to give a result analogous to (13). The \(\ell = 0\) Fourier mode once again confirms Theorem 1, and the higher modes \(\ell = 1 , 2 , \dots , N\) give

$$\begin{aligned} H_{\ell }^{\prime } (t) = - \frac{\ell }{\pi } \sum _{n=1}^N S_{\ell ,n}^{(2)} (t) B_n (t) , \end{aligned}$$
(17)

analogous to the result (15) for the upper fluid.

Similar Fourier analysis is applied to the dynamic condition (5). The series (11) are substituted directly into the condition, which is then evaluated at the interface \(y = \eta \left( x,t \right)\). The expression is multiplied by the basis function \(\cos \left( \ell x \right)\), \(\ell = 1 , 2 , \dots , N\) and integrated over a complete period in x. This gives

$$\begin{aligned} \sum _{n=1}^N C_{\ell ,n}^{(1)} (t) A_n (t) - \sum _{n=1}^N C_{\ell ,n}^{(2)} (t) B_n (t) = \pi \left( D - 1 \right) H_{\ell } (t) . \end{aligned}$$
(18)

In these expressions (15), (17) and (18), it has proven convenient to define

$$\begin{aligned} S_{\ell ,n}^{(1)} (t)= & {} \int _{- \pi }^{\pi } \text {e}^{- n \eta \left( x,t \right) } \sin \left( \ell x \right) \sin \left( nx \right) \, dx \nonumber \\ S_{\ell ,n}^{(2)} (t)= & {} \int _{- \pi }^{\pi } \text {e}^{n \eta \left( x,t \right) } \sin \left( \ell x \right) \sin \left( nx \right) \, dx \nonumber \\ C_{\ell ,n}^{(1)} (t)= & {} \int _{- \pi }^{\pi } \text {e}^{- n \eta \left( x,t \right) } \cos \left( \ell x \right) \cos \left( nx \right) \, dx \nonumber \\ C_{\ell ,n}^{(2)} (t)= & {} \int _{- \pi }^{\pi } \text {e}^{n \eta \left( x,t \right) } \cos \left( \ell x \right) \cos \left( nx \right) \, dx \end{aligned}$$
(19)

to simplify the notation.

The governing equations to be solved for the three sets of Fourier coefficients \(A_n\) , \(B_n\) and \(H_n\) are (15), (17) and (18). However, these form a differential–algebraic inclusion and are not convenient for numerical solution in this present form. We adapt the approach of Forbes et al. [24] and differentiate both the second and third of these equations, to give a simple algebraic system of equations for the derivatives \(A_n^{\prime } (t)\) , \(B_n^{\prime } (t)\) , \(H_n^{\prime } (t)\) of these coefficients. These are then solved, and the entire system of 3N ordinary differential equations is integrated forward in time using standard software.

The second kinematic condition (17) is subtracted from the first condition (15), and the resulting expression is differentiated in time. After a little algebra, this yields

$$\begin{aligned}&\sum _{n=1}^N S_{\ell ,n}^{(1)} (t) A_n^{\prime } (t) + \sum _{n=1}^N S_{\ell ,n}^{(2)} (t) B_n^{\prime } (t) \nonumber \\&\quad = \int _{- \pi }^{\pi } \left[ v_2 u_1 - v_1 u_2 \right] _{y = \eta } \sin \left( \ell x \right) \, dx . \end{aligned}$$
(20)

The dynamic condition (18) is similarly differentiated in time and so gives an expression in the form

$$\begin{aligned}&\sum _{n=1}^N C_{\ell ,n}^{(1)} (t) A_n^{\prime } (t) - \sum _{n=1}^N C_{\ell ,n}^{(2)} (t) B_n^{\prime } (t) \nonumber \\&\quad = \int _{- \pi }^{\pi } \left[ v_2 u_1 - v_1 u_2 \right] _{y = \eta } \frac{\partial \eta }{\partial x} \cos \left( \ell x \right) \, dx + \pi \left( D - 1 \right) H_{\ell }^{\prime } (t) . \end{aligned}$$
(21)

Now the three sets of equations (15), (20) and (21) represent a system of 3N ordinary differential equations for the derivatives of the Fourier coefficients, and we integrate these forward in time using the MATLAB routine for adaptive fourth–fifth-order Runge–Kutta integration. With \(N = 61\) Fourier modes, this runs quickly on a quad-core desktop computer and is straightforward to implement and code.

3 Partially saturated Boussinesq model

If the medium is not fully saturated, and there is in addition a degree of mixing of the fluid in upper layer 1 with that in lower layer 2, then the assumption of a sharp interface between them, as in Sect. 2, may not be appropriate. In that case, the simplest approach is the Boussinesq approximation, in which the fluid density is regarded as a continuously varying function that changes smoothly but rapidly across an interfacial zone of narrow but finite width.

Returning briefly to dimensional variables, the density \(\rho \left( x,y,t \right)\) varies from its value \(\rho _2\) in lower layer 2 continuously to the value \(\rho _1\) in upper layer 1, as y increases from a large negative value to a large positive one. Following Trevelyan et al. [21] and the recent article by De Paoli et al. [25], we regard the larger density \(\rho _1\) of the upper layer as being due to the presence of some additional solute in that upper layer, having concentration \(\left[ C \right] \left( x,y,t \right)\) such that \(\left[ C \right] \rightarrow C_0\) as \(y \rightarrow \infty\). This is therefore relevant to a mineral leaching scenario similar to that envisioned by Forbes [20], for example, in which a mineral of interest is dissolved in situ by a leaching lixiviant. It is assumed here that the fluid density is related to the concentration of the dissolved solute according to the law

$$\begin{aligned} \rho = \rho _2 \left( 1 + \frac{\left( \rho _1 - \rho _2 \right) }{\rho _2 C_0} \left[ C \right] \right) \end{aligned}$$
(22)

and that the solute itself obeys the convection–diffusion equation

$$\begin{aligned} \phi \frac{\partial \left[ C \right] }{\partial t} + \mathbf{q}\cdot \nabla \left[ C \right] = \phi {{\mathcal {D}}} \nabla ^2 \left[ C \right] . \end{aligned}$$
(23)

Here, the constant \({{\mathcal {D}}}\) is a diffusion coefficient for the solute. The quantity \(\phi\) is the porosity of the material and is a dimensionless constant corresponding to a wetted fraction in the medium. The fluid seepage velocity vector \(\mathbf{q}\) is found from Darcy’s law in the form

$$\begin{aligned} \frac{1}{K} \mathbf{q} = - \nabla p - \rho g \mathbf{j} , \end{aligned}$$
(24)

replacing the earlier forms (1), (2) for the fully saturated model. As previously, the continuity equation

$$\begin{aligned} \nabla \cdot \mathbf{q} = 0 \end{aligned}$$
(25)

still holds true, for the fluid in the porous medium.

We return now to dimensionless variables, as previously, and these will be used from now on. As in Sect. 2, all lengths are referred to \(\lambda / \left( 2 \pi \right)\), where \(\lambda\) is the wavelength of a periodic disturbance made to the interface. Speeds are again nondimensionalized by reference to \(\rho _2 g K\), but now it is convenient to make time dimensionless relative to the timescale \(\phi \lambda / \left( 2\pi \rho _2 g K \right)\), so as to absorb the porosity \(\phi\). The density function is referred to the density \(\rho _2\) of the lower fluid, and the pore pressure p is scaled with respect to the quantity \(\rho _2 g \lambda / \left( 2\pi \right)\). The concentration \(\left[ C \right]\) of the solute is scaled relative to its value \(C_0\) in layer 1 far from the interface, and the new dimensionless concentration is written simply as \(C \left( x,y,t \right)\). In these nondimensional variables, the ‘state’ equation (22) becomes simply

$$\begin{aligned} \rho = 1 + \left( D - 1 \right) C \end{aligned}$$
(26)

and the solute concentration satisfies

$$\begin{aligned} \frac{\partial C}{\partial t} + \mathbf{q}\cdot \nabla C = \beta \nabla ^2 C . \end{aligned}$$
(27)

obtained from (23). Darcy’s law (24) becomes simply

$$\begin{aligned} \mathbf{q} = - \nabla p - \rho \mathbf{j} . \end{aligned}$$
(28)

These equations (26)–(28) now involve the two dimensionless parameters

$$\begin{aligned} D = \frac{\rho _1}{\rho _2} \quad \text {and} \quad \beta = \frac{2\pi }{\lambda } \frac{\phi {{\mathcal {D}}}}{\rho _2 g K} . \end{aligned}$$
(29)

The first of these is the density ratio, as before, and the second is a dimensionless diffusion constant for the chemical solute.

From the continuity equation (25), it follows that a streamfunction \(\varPsi \left( x,y,t \right)\) exists, such that

$$\begin{aligned} u = - \frac{\partial \varPsi }{\partial y} \quad \text {;} \quad v = \frac{\partial \varPsi }{\partial x} . \end{aligned}$$
(30)

It then follows from Darcy’s law (28) and the state equation (26) that

$$\begin{aligned} \nabla ^2 \varPsi = - \left( D - 1 \right) \frac{\partial C}{\partial x} . \end{aligned}$$
(31)

The mathematical task is therefore to solve the convection–diffusion equation (27) coupled with (31), for the two functions \(C \left( x,y,t \right)\) and \(\varPsi \left( x,y,t \right)\). The fluid density is then recovered from the state law (26).

3.1 Linearized Boussinesq model

The Boussinesq approximation is predicated on the assumption that the density function \(\rho \left( x,y,t \right)\) does not change much across the different fluid regions. Accordingly, the appropriate small parameter here is the (dimensionless) density difference \(\left( D - 1 \right)\) between the top and bottom levels of the fluid. We therefore postulate

$$\begin{aligned} \varPsi \left( x,y,t \right)= & {} \left( D - 1 \right) \varPsi ^L \left( x,y,t \right) + {{\mathcal {O}}} \left( D - 1 \right) ^2 \nonumber \\ C \left( x,y,t \right)= & {} C^L \left( x,y,t \right) + {{\mathcal {O}}} \left( D - 1 \right) . \end{aligned}$$
(32)

When these approximations are substituted into the governing equations (27), (31) and terms retained only to the first order in \(\left( D - 1 \right)\), the linearized equations

$$\begin{aligned} \frac{\partial C^L}{\partial t}= & {} \beta \left( \frac{\partial ^2 C^L}{\partial x^2} + \frac{\partial ^2 C^L}{\partial y^2} \right) \nonumber \\ \frac{\partial C^L}{\partial x}= & {} - \left( \frac{\partial ^2 \varPsi ^L}{\partial x^2} + \frac{\partial ^2 \varPsi ^L}{\partial y^2} \right) \end{aligned}$$
(33)

are obtained. We seek a solution to these equations, that is \(2\pi\)-periodic in the x-coordinate. For the sake of simplicity, the y-domain will be restricted to the finite-sized interval \(-H< y < H\), and the linearized concentration \(C^L\) in equations (33) solved subject to the conditions

$$\begin{aligned} C^L \left( x, -H , t \right) = 0 \quad \text {;} \quad C^L \left( x, H , t \right) = 1 \end{aligned}$$
(34)

deep within lower fluid 2 and at the top of upper region 1, respectively. Standard methods now give the solution for the linearized concentration in equations (33) to be

$$\begin{aligned} C^L \left( x,y,t \right)= & {} \frac{1}{2} \left( \frac{y}{H} + 1 \right) + \sum _{n=1}^{\infty } A_{0,n} \exp \left( -\beta \varDelta _{0,n}^2 t \right) \nonumber \\&\quad \sin \left( \frac{n\pi }{2 H} \left( y + H \right) \right) \nonumber \\&\quad + \sum _{m=1}^{\infty } \sum _{n=1}^{\infty } A_{m,n} \exp \left( -\beta \varDelta _{m,n}^2 t \right) \nonumber \\&\quad \cos \left( mx \right) \sin \left( \frac{n\pi }{2 H} \left( y + H \right) \right) , \nonumber \\ \end{aligned}$$
(35)

where the additional constants

$$\begin{aligned} \varDelta _{m,n}^2 = m^2 + \left( \frac{n\pi }{2 H} \right) ^2 \end{aligned}$$
(36)

have been defined for convenience. Following Daripa and Hua [26], we observe that an initial condition in which only a single Fourier mode is perturbed continues to involve only that one mode for all later times in the linearized fully saturated model (10). In the Boussinesq model, however, even the linearized solution (35) involves an infinite but discrete spectrum of Fourier modes, when the finite-width computational window \(-H< y < H\) is imposed; this would, in fact, become a continuous spectrum of modes if H were allowed to become infinite (as is the case for the previous linearized solution in Sect. 2.1).

The linearized streamfunction is now obtained from the second of the equations in (33) and becomes

$$\begin{aligned} \varPsi ^L \left( x,y,t \right)= & {} - \sum _{m=1}^{\infty } \sum _{n=1}^{\infty } A_{m,n} \frac{m}{\varDelta _{m,n}^2} \exp \nonumber \\&\quad \left( -\beta \varDelta _{m,n}^2 t \right) \sin \left( mx \right) \sin \left( \frac{n\pi }{2 H} \left( y + H \right) \right) . \end{aligned}$$
(37)

The coefficients \(A_{m,n}\) in these solutions (35), (37) are obtained from initial conditions, and these are now discussed.

Consistent with the linearized saturated solution (10) in Sect. 2.1, we consider the discontinuous initial function

$$\begin{aligned} C \left( x,y,0 \right) = \left\{ \begin{array}{lr} 1 , &{} \text {if } y > \eta _0 \left( x \right) \\ 0 , &{} \text {if } y < \eta _0 \left( x \right) \end{array} \right. \end{aligned}$$
(38)

where we have chosen the initial interface shape to be

$$\begin{aligned} \eta _0 \left( x \right) = K_A \cos x . \end{aligned}$$
(39)

Fourier analysis of (35) gives the coefficients to be

$$\begin{aligned} A_{0,\ell }= & {} \frac{2}{\ell \pi } \cos \left( \ell \pi \right) + \frac{1}{2\pi H} \int _{-\pi }^{\pi } \int _{-H}^H C^L \left( x,y,0 \right) \nonumber \\&\quad \sin \left( \frac{\ell \pi }{2 H} \left( y + H \right) \right) \, dy \, dx \nonumber \\ A_{k,\ell }= & {} \frac{1}{\pi H} \int _{-\pi }^{\pi } \int _{-H}^H C^L \left( x,y,0 \right) \nonumber \\&\quad \cos \left( kx \right) \sin \left( \frac{\ell \pi }{2 H} \left( y + H \right) \right) \, dy \, dx . \end{aligned}$$
(40)

It turns out that the integrals in these expressions can be evaluated in closed form, taking real and imaginary parts of the identity given in Abramowitz and Stegun [27, page 360, formula 9.1.21]. After some algebra, we obtain

$$\begin{aligned} A_{0,\ell }= & {} \frac{2}{\ell \pi } \cos \left( \frac{\ell \pi }{2} \right) J_0 \left( \frac{\ell \pi K_A}{2 H} \right) \nonumber \\ A_{k,\ell }= & {} \frac{4}{\ell \pi } \cos \left( \frac{\pi }{2} \left( k + \ell \right) \right) J_k \left( \frac{\ell \pi K_A}{2 H} \right) , \quad k \ge 1 , \end{aligned}$$
(41)

where \(J_k (z)\) denotes the Bessel function of the first kind of order k.

Fig. 2
figure 2

Initial density profile, reconstructed from the Fourier series (35) with smoothed coefficients (42). The density function is then formed from the state equation (26). Here, \(K_A = -0.2\) and the density ratio is \(D = 1.2\)

Full size image

From the numerical point of view, the use of the discontinuous initial concentration profile (38), with the cosinusoidal profile (39) is undesirable because it gives Fourier coefficients (41) that decay as \(A_{k,\ell } \sim 1 / \ell\) for large \(\ell\). Consequently, Gibbs’ phenomenon [28] means that the Fourier series (35) cannot recreate the initial profile (38). Instead, it replaces the discontinuity with a region containing small-amplitude oscillations. To avoid this difficulty, we use Lanczos smoothing of the Fourier coefficients (41), replacing them with new coefficients

$$\begin{aligned} {\overline{A}}_{m,n} = A_{m,n} \frac{\sin \left( n\sigma \right) }{n\sigma } . \end{aligned}$$
(42)

The Lanczos parameter \(\sigma\) is a small positive constant, and we have taken it to be \(\sigma = 0.05\) here. This process (42) can be shown to be equivalent to replacing the point value of the initial function (38) with a moving average over a window of width \(2\sigma\) centered at that point. Thus, the discontinuous initial condition (38) is replaced by a smoothed, continuous function. This is illustrated in Fig. 2, for an initial profile with \(K_A = -0.2\) and using \(M = N = 101\) modes in each coordinate. The density has been generated from the state law (26) with \(D = 1.2\), and as is evident from the diagram, a smooth but rapid transition from \(\rho = 1\) to \(\rho = D\) is obtained near the interfacial region. A careful examination shows that small wavelets of Gibbs type are still present near the interface, but are too small to be of concern. They can be eliminated entirely by increasing \(\sigma\) slightly, although at the cost of diffusing the interfacial zone still further.

Fig. 3
figure 3

Contours of the density \(\rho\) for the linearized Boussinesq solution, for the four dimensionless times \(t = 4\), 16, 28 and 40. Here, density ratio \(D = 1.2\) , diffusion coefficient \(\beta = 0.001\) and initial wave amplitude \(K_A = -0.2\). The solution (35) was evaluated using \(M = N = 151\) Fourier modes

Full size image

Figure 3 illustrates the linearized solution with concentration \(C^L \left( x,y,t \right)\) calculated from (35) and then density obtained from the state law (26). The solution was evaluated using \(M = N = 151\) Fourier modes on a grid of \(601 \times 601\) points over the computational domain \(-\pi< x < \pi\), \(-H< y < H\) with \(H = 6\). Contours of the density \(\rho\) are shown at the four times \(t = 4\), 16, 28 and 40. In the first two diagrams, there are still some very small wavelets present, as a remnant of Gibbs’ phenomenon present in the choice of initial condition illustrated in Fig. 2; these have no effect upon the solution, although the contouring routine in the graphics package nevertheless detects these small ripples and displays them as the fine corrugations visible at the first two times shown.

It is evident from Fig. 3 that the interfacial zone essentially does not move at all, but remains in place and diffuses slightly due to the effect of the parameter \(\beta\) in equation (33). This is as expected, since the process (32) of linearization results in the (nonlinear) self-convection terms being removed, in the linearized approximation (33) for the concentration function \(C^L\) and hence the density \(\rho\). As a result, the linearized Boussinesq theory in this section is possibly of only limited interest, and the nonlinear theory is now considered.

3.2 Nonlinear Boussinesq model

We now consider extending the mathematical structure revealed by the linearized solution (35), (37) to the full nonlinear equations (26)–(28), thus creating a spectral semi-numerical solution process. An appropriate form for the nonlinear concentration function is

$$\begin{aligned} C \left( x,y,t \right)= & {} \frac{1}{2} \left( \frac{y}{H} + 1 \right) + \sum _{n=1}^N A_{0,n} (t) \sin \left( \frac{n\pi }{2 H} \left( y + H \right) \right) \nonumber \\+ & {} \sum _{m=1}^M \sum _{n=1}^N A_{m,n} (t) \cos \left( mx \right) \sin \left( \frac{n\pi }{2 H} \left( y + H \right) \right) , \end{aligned}$$
(43)

and the streamfunction is represented in the form

$$\begin{aligned} \varPsi \left( x,y,t \right)= & {} - \left( D - 1 \right) \sum _{m=1}^M \sum _{n=1}^N A_{m,n} (t) \frac{m}{\varDelta _{m,n}^2} \nonumber \\&\quad \sin \left( mx \right) \sin \left( \frac{n\pi }{2 H} \left( y + H \right) \right) . \end{aligned}$$
(44)

In these representations (43), (44) of the nonlinear solution, the integers M and N denote the numbers of Fourier modes able to be kept in a numerical implementation of these series; an ‘exact’ solution would require these numbers to become infinite. The constants \(\varDelta _{m,n}^2\) are as defined in (36). Form (44) has been chosen because it automatically satisfies the (linear) Poisson-type equation (31).

It merely remains to determine the coefficient functions \(A_{m,n} (t)\) so as to satisfy the nonlinear convection–diffusion equation (27). The series (43), (44) are substituted in, and the resulting equation is subjected to standard Fourier analysis. After some algebra, this yields the system of nonlinear ordinary differential equations

$$\begin{aligned} A_{0,\ell }^{\prime } (t)= & {} - \beta \varDelta _{0,\ell }^2 A_{0,\ell } (t) - \frac{1}{2\pi H} \int _{-\pi }^{\pi } \int _{-H}^H \nonumber \\&\quad \left( u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} \right) \sin \left( \frac{\ell \pi }{2 H} \left( y + H \right) \right) \, dy \, dx \nonumber \\ A_{k,\ell }^{\prime } (t)= & {} - \beta \varDelta _{k,\ell }^2 A_{k,\ell } (t) - \frac{1}{\pi H} \int _{-\pi }^{\pi } \int _{-H}^H \nonumber \\&\quad \left( u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} \right) \cos \left( kx \right) \sin \left( \frac{\ell \pi }{2 H} \left( y + H \right) \right) \, dy \, dx . \nonumber \\ \end{aligned}$$
(45)

To begin the computation, initial values \(A_{0,\ell } (0)\) and \(A_{k,\ell } (0)\) are chosen, and in fact the values given in (41) are sufficient to generate the initial profile in Fig. 2. This system (45) of nonlinear differential equations is then integrated forward in time using the explicit adaptive Runge–Kutta–Fehlberg routine ode45 in the MATLAB programming environment. The integrals on the right-hand sides of (45) are performed using the Gaussian quadrature package lgwt.m written by von Winckel [29] which, with the number of mesh-points used here, is effectively exact. The two velocity components u and v in the integrands are evaluated directly by differentiation of the streamfunction (44) using the results (30). Although MATLAB is not a particularly fast package, we nevertheless find that we are able to generate highly converged numerical results with \(M = N = 101\) Fourier coefficients, and a numerical grid of \(501 \times 501\) points over the computational domain, in about 21 hours run-time on a modest quad-core desktop computer. We have performed careful convergence tests using 51 and 81 Fourier modes and find that there are no significant differences between results with \(M = N = 81\) and \(M = N = 101\).

4 Results

Fig. 4
figure 4

Density contours for the nonlinear Boussinesq solution, for density ratio \(D = 1.2\), with \(\beta = 0.001\) and initial amplitude \(K_A = -0.2\), at dimensionless time \(t = 1.2\). The dashed line is the interface \(\eta \left( x,t \right)\) calculated with the fully saturated nonlinear model in Sect. 2.2. The scales on the axes are equal

Full size image

When the fully saturated, nonlinear model of Sect. 2.2 is investigated, it is soon found that a disturbance to the interface initially grows rapidly, as predicted by the linearized saturated solution in Sect. 2.1, but then fails abruptly at some finite time. As observed by Sharp [3] for the traditional Rayleigh–Taylor instability, there does not appear to be anything pathological about the nonlinear interface profile at the time at which it fails, and it is only upon an examination of the interfacial curvature that a singularity is revealed, like that first predicted for the Kelvin–Helmholtz instability by Moore [4].

For density ratio \(D = 1.2\) and the initial cosinusoidal disturbance (39) with amplitude parameter \(K_A = -0.2\), the fully saturated nonlinear solution is able to progress up until about \(t = 1.2\), after which a singularity forms abruptly in the curvature, and the saturated solution fails. We show in Fig. 4 the nonlinear interface at that last time \(t = 1.2\), where it is sketched using a heavy dashed line. By way of comparison, we have overlaid this dashed curve on a field of contours of the density function \(\rho \left( x,y,t \right)\) obtained for the nonlinear Boussinesq model and computed from equations (43) and (26), with the same parameter values and the smoothed initial condition in Fig. 2. The agreement between the results of these two different models is very good, and the interface curve for the fully saturated case lies precisely along the mid-density contour \(\rho = \left( 1 + D \right) / 2\) obtained with the Boussinesq model. Nevertheless, the fully saturated results cannot continue significantly past this early time \(t = 1.2\). This could no doubt be addressed using a ‘vortex blob’ approximation [9], which effectively smears the exact interface over a finite-width zone. It may also be possible to avoid singularity formation at the interface by adding extra viscous fluid layers there, as has been suggested by Daripa [30]. Such an approach is possibly equivalent to a ‘vortex blob’ method, to the extent that it effectively makes the interfacial zone of finite width, but might be more difficult to implement in a nonlinear problem. Alternatively, filtering methods and regularization techniques can be used to suppress singularity formation, and an excellent discussion of these is presented by Daripa and Hua [26]; regularization by apodization was in fact used in a viscous interface computation by Forbes et al. [7]. In this paper, however, we continue to later times using the Boussinesq model in Sect. 3.

Fig. 5
figure 5

a Contours of the density \(\rho\) and b contours of the streamfunction \(\varPsi\), for the nonlinear Boussinesq solution, at the four dimensionless times \(t = 4\), 16, 28 and 40. Here, density ratio \(D = 1.2\) , diffusion coefficient \(\beta = 0.001\) and initial wave amplitude \(K_A = -0.2\). The solution was computed using \(M = N = 101\) Fourier modes

Full size image

A nonlinear Boussinesq solution is illustrated in Fig. 5, at the four dimensionless times \(t = 4\), 16, 28 and 40 (as for the linearized solution in Fig. 3). Contours of the density \(\rho\) are shown in part (a) and were computed from (43) and (26) using \(M = N = 101\) Fourier modes with \(501 \times 501\) mesh points over the computational domain. The lower limit of the computational window was set to \(y = - H = - 6\). As for the linearized solution shown in Fig. 3, the density ratio is \(D = 1.2\), the diffusion coefficient \(\beta = 0.001\) and the initial wave amplitude is \(K_A = -0.2\). The initial condition was the smoothed profile illustrated in Fig. 2. Contours of the streamfunction \(\varPsi\) are illustrated in Fig. 5b; since this is an unsteady problem, there is no guarantee that these contours match exactly to the streamlines, but they do serve at least as a guide as to how the fluid velocity patterns develop.

At very early times, the nonlinear solution matches closely to the linearized solution of Sect. 3.1. However, the process (32) of linearization removes the self-convection terms in (27), so that the linearized interface profile then merely remains in about its original position, as seen in Fig. 3. However, the nonlinear interface profiles shown in Fig. 5a exhibit a large-amplitude downward-moving tongue of the heavier fluid from upper layer 1, and to compensate, two bubbles of lighter fluid 2 moving upward either side of it. This is qualitatively similar to nonlinear interfacial behavior in the conventional Rayleigh–Taylor instability [31], except that in the porous medium considered here, there is no pronounced overturning of the tongue near its tip. Some evidence of the very small amplitude ripples from the initial density profile sketched in Fig. 2 is still visible in Fig. 5a, and their appearance has been enhanced somewhat by the contouring routine. They are very small and have no effect on the numerical results.

Contours of \(\varPsi\) for this same solution are drawn for the same four times in part(b). We shall refer to these as ‘streamlines’ subsequently, although this may not be exactly true for these unsteady flows. Clearly, this entire flow comes about as a result of the nonlinear convection terms in (27), and it can be seen from part (b) that the unstably growing tongue and bubbles in part (a) are associated with the formation of a pair of vortices in the flow. The strength of these vortices increases with time, as the unstable interface develops perturbations of larger amplitude.

Fig. 6
figure 6

a Contours of the density \(\rho\) and b contours of the streamfunction \(\varPsi\), for the nonlinear Boussinesq solution, at the four dimensionless times \(t = 4\), 16, 28 and 40. Here, density ratio \(D = 1.5\) , diffusion coefficient \(\beta = 0.001\) and initial wave amplitude \(K_A = -0.2\). The solution was computed using \(M = N = 101\) Fourier modes

Full size image

In Fig. 6, we increase the density ratio to the value \(D = 1.5\). The diffusion coefficient \(\beta = 0.001\) and initial wave amplitude parameter \(K_A = -0.2\) are as previously. Contours of the density \(\rho\) are displayed in part (a) and ‘streamlines’ are drawn in part (b), both for the same four times \(t = 4\), 16, 28 and 40 as in Figs. 3 and 5. Small Gibbs phenomenon wavelets may be visible in the density contours in part (a), but again, they are of sufficiently small amplitude as to have no discernible effect on the numerical results.

At this greater value \(D = 1.5\) of the density ratio, the instability is stronger and so the deformation of the interface is more extreme. By time \(t = 28,\) the downward-moving spike of heavier upper fluid 1 takes up almost all the vertical computational domain \(-H< y < H\) with \(H = 6\), with a corresponding observation for the upward-moving bubbles of lighter fluid 2 either side of it. For interest, we have also included the numerical solution obtained with \(t = 40\) for this case, since it demonstrates that there has now been complete inversion of the two fluid layers, with the lighter fluid 2 now occupying the upper portion and heavier fluid 1 the lower portions of the diagram. Nevertheless, it must be admitted that, while broadly correct, these last profiles at time \(t = 40\) are possibly not entirely what would be observed for porous Rayleigh–Taylor instability in a finite-sized domain with a rigid bottom boundary, since the assumed form (43) still satisfies the original boundary conditions (34), in which \(C = 0\) at the bottom \(y = -H\) of the region, and \(C = 1\) at the top \(y = H\). This explains the presence of the small vertical line from the top of the structure in part (a), at \(t = 40\), connecting it to the top of the computational region. The streamlines, shown in part (b), again demonstrate how buoyancy effects create small vortices in the flow, which grow in size and strength as the instability develops.

So far, the initial condition (38) with the interface shape (39) has involved a single-mode perturbation. For interest, we consider now a multimode disturbance as the initial perturbation, and we choose to replace the simple cosinusoidal shape (39) with the triangular initial perturbation

$$\begin{aligned} \eta _0 \left( x \right) = \left\{ \begin{array}{lr} K_A \left( \frac{| x |}{q\pi } - 1 \right) , &{} \text {if } 0< | x |< q\pi \\ 0 , &{} \text {if } q\pi< | x | < \pi \end{array} \right. \end{aligned}$$
(46)

in which the constant q lies in the interval \(0< q < 1\). The triangular shape (46) is repeated periodically every \(2\pi\) units along the x-axis. Admittedly, this periodic triangular initial perturbation (46) is somewhat unrealistic from a physical point of view, but it is multimodal and does also allow the coefficients of the initial Fourier representation to be calculated in closed form using (40). After some algebra, the coefficients in the representation (43) at \(t = 0\) are determined to be

$$\begin{aligned} A_{0,\ell } \left( 0 \right)= & {} \frac{2}{\ell \pi } \left( 1 - q \right) \cos \left( \frac{\ell \pi }{2} \right) \nonumber \\&\quad + \frac{4 q H}{ \left( \ell \pi \right) ^2 K_A} \left\{ \sin \left[ \frac{\ell \pi }{2} \left( \frac{K_A}{H} - 1 \right) \right] + \sin \left( \frac{\ell \pi }{2} \right) \right\} \nonumber \\ A_{k,\ell } \left( 0 \right)= & {} - \frac{4}{k\ell \pi ^2} \cos \left( \frac{\ell \pi }{2} \right) \sin \left( k q\pi \right) \nonumber \\&\quad + \frac{4Hq}{\ell \pi ^ 2 \left( \ell K_A - 2H k q \right) } \left\{ \sin \left[ \pi \left( \frac{\ell }{2} - kq \right) \right] \right. \nonumber \\&\quad \left. + \sin \left[ \frac{\ell \pi }{2} \left( \frac{K_A}{H} - 1 \right) \right] \right\} \nonumber \\&\quad + \frac{4Hq}{\ell \pi ^ 2 \left( \ell K_A + 2H k q \right) } \left\{ \sin \left[ \pi \left( \frac{\ell }{2} + kq \right) \right] \right. \nonumber \\&\quad \left. + \sin \left[ \frac{\ell \pi }{2} \left( \frac{K_A}{H} - 1 \right) \right] \right\} , \nonumber \\&\quad k \ge 1 . \end{aligned}$$
(47)

Once again, these coefficients (47) are smoothed using (42), giving a smoothed initial profile and minimizing the effects of Gibbs’ phenomenon.

Fig. 7
figure 7

a Contours of the density \(\rho\) and b contours of the streamfunction \(\varPsi\), for the nonlinear Boussinesq solution, at the four dimensionless times \(t = 4\), 16, 28 and 40. Here, density ratio \(D = 1.2\) , diffusion coefficient \(\beta = 0.001\) and initial wave amplitude \(K_A = 0.2\). The solution was computed using \(M = N = 101\) Fourier modes, from the triangular initial perturbation (47) with \(q = 1/3\))

Full size image

The triangular-shaped initial condition (47) with amplitude \(K_A = 0.2\) is used in Fig. 7 to generate porous Rayleigh–Taylor solutions for density ratio \(D = 1.2\), and some results are shown at the four times \(t = 4\), 16, 28 and 40. Figure 7a presents contours of the density \(\rho\) and part (b) illustrates some streamlines at the different times. At the earliest time \(t = 4\) shown, the triangular perturbation shape to the interface is still visible, but as time progresses, it develops into a long finger extending downward. As previously, two upward-moving bubbles of the lighter lower fluid 2 are formed either side of it, so that the average interface location remains approximately constant; this is broadly consistent with Theorem 1, even in spite of the fact that the theorem was derived for the previous fully saturated model.

The streamline patterns in Fig. 7b are similar to those obtained using the cosinusoidal initial condition (39) at earlier times, illustrated in Fig. 5b. The downward-moving spike is enabled by the formation of two vortices in the flow, as can be seen at the first two times \(t = 4\) and \(t = 16\). However, at later times, each vortex further splits into two, giving a system of four vortices, associated with both the spike of heavy fluid moving down as well as the two bubbles moving up. These are visible in the last two frames in Fig. 7b.

Fig. 8
figure 8

a Contours of the density \(\rho\) and b Contours of the streamfunction \(\varPsi\), for the nonlinear Boussinesq solution, at the four dimensionless times \(t = 4\), 16, 28 and 40. Here, density ratio \(D = 1.5\) , diffusion coefficient \(\beta = 0.001\) and initial wave amplitude \(K_A = 0.2\). The solution was computed using \(M = N = 101\) Fourier modes, from the triangular initial perturbation (47) with \(q = 1/3\)

Full size image

For interest, we also include in Fig. 8 the similar solution to Fig. 7, except that now the density ratio has been increased here to the value \(D = 1.5\). Density contours are again shown in part(a) and streamlines are drawn in part (b), for the same four times \(t = 4\), 16, 28 and 40 as previously. As a result of the increased density ratio, the instability at the interface develops more rapidly, to the extent that the downward-moving spike has already reached the bottom of the computational region \(y = -H = -6\) by time \(t = 28\). This is essentially a porous Rayleigh–Taylor flow in a confined region \(-H< y < H\) with an impermeable lower boundary at \(x = - H\), although there is a very narrow boundary layer near both the top and the bottom of the region \(y = \pm H\) where the solution reverses rapidly to its initial boundary conditions \(C = 0\) at \(y = -H\) and \(C = 1\) at \(y = H\), since these are demanded by the form of the solution (43). Thus, the solutions at the last two times \(t = 28\) and \(t = 40\) will not correspond exactly to unstable flow in a confined region because of the presence of these boundary layers at top and bottom, although the other flow features are correct.

Fig. 9
figure 9

Density contours for the nonlinear Boussinesq solution, for density ratio \(D = 1.5\), with \(\beta = 0.001\) and initial amplitude \(K_A = 0.2\), at dimensionless time \(t = 24\). The solution was started from the triangular initial pulse (46) with \(q = 1/3\). The scale is the same on both axes

Full size image

We show in Fig. 9 some contours of the density function \(\rho\) for the same solution as in Fig. 8, but now at the time \(t = 24\). The scale on the horizontal and vertical axes is the same, so that this droplet is as it would actually appear. This is close to the time at which the downward-moving spike of heavier fluid meets the bottom of the enclosure (or else the bottom \(y = -H\) of the computational domain), and it shows how additional, secondary droplets have formed either side of the main spike and are also being convected downward. Two bubbles of lower lighter fluid 2 have also moved deeply up into the above heavier layer 1.

5 Conclusions

This paper is part of a special issue on the general topic ‘Interfaces, Mixing and Non-Equilibrium Dynamics,’ and the authors are grateful for this opportunity to contribute. The Rayleigh–Taylor instability is one of the three most familiar instabilities that can occur at an interface and is known in many applications and over vastly differing length scales [13] (the other two are the Kelvin–Helmholtz and the Richtmyer–Meshkov flows). It is responsible for mixing of the two fluids, where it may eventually exhibit some similarities to turbulent mixing [19].

We have examined two models for Rayleigh–Taylor unstable flow in porous media, both in a linearized approximation in which disturbances to the otherwise horizontal interface are small, and in the fully nonlinear equations pertinent to that model. The first such model was one in which the medium was assumed to be fully saturated with the fluids, and there are strong similarities between that model and the classical Rayleigh–Taylor instability. The linearized theory in that case predicted an interface that is unstable for any density ratio \(D > 1\) for which the upper fluid is heavier than the lower one, so that perturbations grow exponentially with time. The nonlinear results bear this out, but only for early times, since there is a finite time at which the interface develops a curvature singularity and the solution fails. This same behavior also occurs for the classical nonlinear Rayleigh–Taylor flow, as shown by Baker et al. [6].

The second model was one in which the rock is supposed to be only partially saturated, but that a dissolved chemical in the upper fluid was responsible for making it more dense than the lower fluid. This again leads to a Rayleigh–Taylor situation with a heavy fluid above a lighter one, and we used a Boussinesq approximation with an interfacial zone of finite width to continue the numerical solution out for times far later than the critical time at which the corresponding fully saturated model would have failed. Elongated spikes of heavier upper fluid are found to form, alternating with upwardly moving bubbles of the lighter fluid. Unlike the classical Rayleigh–Taylor flow in which the spikes and bubbles form into mushroom-shaped regions with strongly overhanging heads [31], the analogous structures formed here do not possess similar overturning regions near their heads.