An extended Boussinesq theory for interfacial fluid mechanics

Boussinesq theory can model quite accurately viscous flows that involve multiple fluids with interfaces between them, so long as there is not much difference between the densities of the various fluids. However, the Boussinesq approximation is generally poor when the density ratio between the fluids is large. Here, we propose an Extended Boussinesq approximate equation, that allows for large density ratios, while still remaining straightforward to implement. Examples are given for planar Rayleigh–Taylor instability, where the Boussinesq and the novel Extended Boussinesq models are compared with the predictions of an SPH fluid dynamics code, to confirm this approach.

Much has been learned about free-surface and interfacial flows by ignoring the complicating effects of viscosity. Then the fluid flows irrotationally, and so its velocity vector q can be expressed as the gradient of a scalar potential . Furthermore, the incompressibility of the fluid results in Laplace's equation ∇ 2 = 0 for the potential, and this is a linear partial differential equation. Nevertheless, the overall problem is still nonlinear, since the location of the free surface remains unknown. However, for small-amplitude disturbances of the surface, the problem can often be linearized about some base flow, and solved in closed form. The classical solutions for free-surface flow over topography in the book by Lamb [3], and for wave patterns behind moving ships in Wehausen and Laitone [1] have been obtained in this manner. Furthermore, the fact that Laplace's equation is obeyed throughout the fluid allows elegant analytical techniques to be applied in designing numerical methods even for nonlinear problems involving large surface deformations, and some of these are surveyed by Forbes et al. [4] in a recent review article on free-surface flow over bottom-mounted obstacles.
When there are two inviscid fluids present, with an interface between them, there are situations under which the interface is unstable to small perturbations, and therefore develops large-amplitude plumes or overhanging mushroom-shaped profiles as time progresses. Of these, the three most well known are the Kelvin-Helmholtz, the Rayleigh-Taylor and the Richtmyer-Meshkov instabilities. The first of these occurs when two inviscid fluids flow past one another at different speeds. This is thus a shearing instability, under which the interface develops rolled-up "cat's eye" spirals. The text by Chandrasekhar [5,Chap. 11] discusses Kelvin-Helmholtz instabilities under various circumstances; linearized theory shows that the amplitude of the interface grows exponentially, but for inviscid fluids, Moore [6] has shown that nonlinear effects cause a curvature singularity to be formed at the interface within finite time. The second of these interfacial instabilities is the Rayleigh-Taylor situation, in which a heavy fluid overlies a light one. Small disturbances to the interface grow exponentially with time until, again, a Moore curvature singularity forms at some finite critical time in the inviscid theory. Forbes [7], among many others, demonstrated numerically that when viscosity and finite interface thickness are taken into account, the curvature singularity is replaced by a small patch of high vorticity in the fluid, and this is responsible for the interface rolling over to form a mushroom-shaped segment. The Rayleigh-Taylor instability can be found over vast differences in length and time scales, from the microscopic to the galactic [8]. The third well-known interfacial instability is the Richtmyer-Meshkov type, which arises when an interface between two fluids is perturbed by a shock wave. A review of this topic is presented by Brouillette [9] and Matsuoka and Nishihara [10] have computed some nonlinear interface profiles.
Curvature singularities of the type discovered by Moore [6] have been avoided in the "vortex-blob" method of Krasny [11]. This heuristic technique effectively "smears" the interface over a finite thickness, and allows mushroomshaped structures to develop in an otherwise inviscid formulation. Matsuoka and Nishihara [10] followed a similar approach, adding a "regularizing" parameter δ 2 in the denominator of their integral equation, effectively smearing the interface over a width 2δ. Nevertheless, Baker and Pham [12] demonstrated that different types of vortexblob methods can give different outcomes as time progresses. In order to overcome the heuristic nature of these techniques, Forbes et al. [13] investigated a two-fluid system with a sharp interface, and included the effects of viscosity. They found, however, that viscosity alone is insufficient to avoid a type of curvature singularity at the interface, and that some sort of diffusion of the interface was also required. This explains, in part, the success of the viscous Boussinesq theory used by Forbes [7] to study the Rayleigh-Taylor instability, since it allows a small degree of mixing of the fluids at the interfacial zone.
Boussinesq theory is derived on the assumption that the two fluids in the system have densities that are very nearly equal. Lee and Kim [14] undertook an investigation of the results produced by Boussinesq theory and compared them with numerical solutions obtained using a more elaborate variable-density model in which the fluid density was a function of a "phase variable" φ that satisfied a nonlinear fourth-order advection-diffusion equation based on the Cahn-Hilliard equation from reaction chemistry. They found that the two theories agreed well when the density ratio D = ρ 2 /ρ 1 of the two fluids was nearly one, but that the results differed strongly as D increased. In particular, for the Rayleigh-Taylor problem studied by Forbes [7], Lee and Kim [14] showed that, while the interface profiles predicted by Boussinesq theory remain symmetric when viewed either from above or below for all density ratios D, their profiles generated with the Cahn-Hilliard model become less and less symmetric as D increases. Furthermore, their findings are at least in qualitative agreement with the experiments conducted by Waddell et al. [15]. Nevertheless, it is perhaps not obvious that the Cahn-Hilliard equation is necessarily appropriate for all pairs of fluids exhibiting a Rayleigh-Taylor instability event, and Yakovenko [16] has instead carried out similar computations using a system of equations in which the density and the fluid viscosity are functions of a volume fraction f , that is itself convected through the fluid. They similarly find that the top-bottom symmetry of the interfacial profiles that occurs for D near one is progressively lost as D increases. Talat et al. [17] simulate the Rayleigh-Taylor instability using a "phase-field" method based on the Cahn-Hilliard equation, in a formulation similar to that of Lee and Kim [14], and they likewise obtain interfacial profiles that become less symmetric as D increases.
In the present paper, our aim is similarly to derive a model equation that allows us to compute interfacial flows in multi-fluid systems. The correct approach is, of course, to apply the Navier-Stokes equations to each fluid and the viscous interfacial boundary conditions on each interface, given in Batchelor [18, p. 150]. However, such an approach is prohibitively difficult, and particularly in overturning interfaces or cases in which the interface forms a separate, detached structure. We seek to retain the simpler, approximate, approach afforded by the Boussinesq system, but allow larger density ratios D and to remain as faithful as possible to the original Navier-Stokes system. We briefly review the Boussinesq approach in Sect. 2, and use this as a basis for our Extended Boussinesq system in Sect. 3. The system is then applied to two-dimensional Rayleigh-Taylor flow in Sect. 4 and the results are shown in Sect. 5, where they are also compared with the output from a Smoothed-Particle Hydrodynamics (SPH) program, described in Sect. 5.2. The paper concludes in Sect. 6 with a discussion.

Review of Boussinesq theory
We consider a system of two fluids, as illustrated in Fig. 1. We suppose they have the same dynamic viscosity μ, but differing densities ρ 1 and ρ 2 and there is an interface separating them. In Boussinesq theory, the two densities are regarded as sufficiently similar that their density ratio D = ρ 2 /ρ 1 is close to one. In addition, Boussinesq theory approximates the two-fluid system as a single fluid in which the density varies slightly throughout the medium; this means that an interface between the two fluids is represented approximately as a narrow zone across which the density varies smoothly, although rapidly, from ρ 1 on one side to ρ 2 on the other. The density of this approximating medium is thus expressed as a continuous function ρ (x, y, z, t), which is then further represented in the form with ρ taking values between zero and ρ 2 − ρ 1 . The perturbed density ρ is thus considered a "small" quantity. If q (x, y, z, t) denotes the fluid velocity vector, the continuity equation of mass conservation is approximated, in Boussinesq theory for an incompressible fluid, by "splitting" it into a component div q = 0 for the quantities of order one in magnitude, and a remaining component for the "small" quantities. The zero on the right-hand side of (2) has been replaced with the term K ∇ 2 ρ to model diffusion, and the constant K is related to a Prandtl number, as indicated by Farrow and Hocking [19]. An alternative way of deriving these relations is to assume that the fluid is genuinely incompressible, so that (3) holds without approximation, and then to assume that density differences in the fluid are due to the presence of some dissolved mineral that has concentration C. If it does not react with the fluid, then the concentration C of this mineral obeys a convection-diffusion equation identical to (4) for the quantity ρ. The Boussinesq approximation is equivalent to supposing that the concentration and the density are related by a simple linear "equation of state" which is consistent with the models of Trevelyan et al. [20] and Forbes et al. [21] for Rayleigh-Taylor flows in porous media. Finally, the Navier-Stokes equation of viscous fluid flow is approximated by the simpler equation Here, the body force per mass acting on the fluid is f. The only difference between this and the full Navier-Stokes equation is that the first term (in square parentheses) in (5) is multiplied only by the constant ρ 1 rather than by the full density (ρ 1 + ρ). In many studies, the pressure p is eliminated from (5) by taking the vector curl of the equation. The incompressibility condition (3) is satisfied identically using a vector velocity potential A. In addition, the usual vorticity vector ζ is introduced, so that q = curl A and ζ = curl q.
When the flow is two dimensional in the (x, y) plane, it is sufficient simply to take A = k with unit vector k pointing in the positive z-direction and the single component function (x, y, t) as the streamfunction. This approach was adopted by Forbes [7] for the viscous Rayleigh-Taylor problem. When, however, the flow is fully three dimensional, the vector velocity potential A needs at least two components, so that the velocity is expressed in terms of two streamfunctions and χ and this technique was used by Walters and Forbes [22], for example, to study three-dimensional Rayleigh-Taylor flow. Now the vector curl of (5) gives The last term on the right-hand side of this equation will vanish if the body force is conservative, since then curl f = 0.
A solution to this Boussinesq model consists of finding the streamfunction(s) that make up the components of the vector potential A and the density perturbation ρ that satisfy (4) and (7).

Derivation of extended Boussinesq model
Boussinesq theory might be expected to work well for density ratios D = ρ 2 /ρ 1 near one. However, since the approximation (ρ 1 + ρ) ≈ ρ 1 was invoked to obtain the Boussinesq momentum equation (5), this theory cannot retain accuracy if D becomes large. We therefore wish to derive a more accurate theory that can allow larger variations in density, while still retaining the simplifying feature of Boussinesq theory, that the two-fluid system is represented as a single medium in which the density changes smoothly and rapidly at an interface. It is, in fact, possible to retain the full Navier-Stokes equation, without the approximation described above, and again eliminate the pressure p by taking the vector curl. This gives a difficult type of vorticity equation, in which the time-derivative term ∂ζ /∂t on the left-hand side of (7) no longer appears by itself, but is multiplied by other unknown functions. Furthermore, an additional term containing the time derivative ∂q/∂t is also present, and is similarly multiplied by a complicated unknown function. We have designed and implemented an apparently successful numerical scheme based around the use of this equation; however, it is prohibitively complicated and slow, since integration forward in time requires the creation and inversion of a very large matrix, at each new time step. Such an approach evidently offers no advantages over the use of package computational fluid dynamics (CFD) programs. Instead, the aim here is to retain the basic philosophy and advantages of the Boussinesq approach, while also extending the approximation to some higher order of accuracy, and all within moderate computational requirements.
With the Boussinesq view (1) of the fluid density, the exact Navier-Stokes equation is As in Sect. 2, we take the vector curl of (8) and obtain, after some algebra Once again, the last term vanishes if the body force per mass is conservative. So far, this equation (9) is exact, and does not involve the severe approximation made to the full density term that multiplies the time derivative in (5). However, the pressure p is still present, and some approximation is therefore required.
To proceed, we have supposed that since only the gradient of the pressure is involved in (9), we can eliminate this quantity using the hydrostatic approximation

Consequently, this gives the approximate vorticity equation
The incompressibility condition (3) and transport equation (4) still apply.

The Rayleigh-Taylor instability
We now investigate the application of the standard Boussinesq theory of Sect. 2 and the Extended Boussinesq theory of Sect. 3 to the classical Rayleigh-Taylor instability in two-dimensional flow. We consider the model problem in Forbes [7] in which there are two horizontal layers of fluid; lower fluid 1 has density ρ 1 and dynamic viscosity μ and is located over −H 1 < y < 0, and above it is upper layer 2 with density ρ 2 and the same viscosity μ, located in 0 < y < H 2 . There is an interface along the x-axis (y = 0) and the flow is confined between the lower and upper walls at y = −H 1 and y = H 2 . We now suppose that the interface is perturbed by some small periodic disturbance, and takes up the shape y = A I cos (π x/λ). This disturbance has initial amplitude A I and wavelength λ and is symmetric about the y-axis, and so we consider only the segment −λ/2 < x < λ/2. Dimensionless variables are introduced now, using λ/(2π) as the length scale and √ λ/ (2πg) as the scale for time. All speeds are then referenced to the quantity √ gλ/(2π). The density ρ 1 of the lower layer is used to make the density perturbation variable ρ dimensionless, and pressure is referred to (ρ 1 gλ) /(2π). Here, g is the (constant) acceleration of gravity, directed along the negative y-axis, so that the body-force vector in Sects. 2 and 3 is f = −gj. A sketch of this system, in these new dimensionless coordinates, is given in Fig. 1. Solutions to this problem, in either of the two models above, are thus dependent on the six dimensionless parameters The first two constants h 1 and h 2 give the locations of the two horizontal walls, and the amplitude and density ratio D appear in the initial condition The constant σ is an effective Prandtl number [19] and R e is a Reynolds number for the flow. For this two-dimensional flow, illustrated in Fig. 1, the velocity vector can be represented q = ui + vj and the vorticity vector has only a single component ω, and so can be expressed as ζ = ωk. It then follows from (6) that in which (x, y, t) is the streamfunction. Both models then satisfy the transport equation (4) in the form In the Boussinesq model, the vorticity equation (7) for this Rayleigh-Taylor flow takes the form Alternatively, the extended Boussinesq model of Sect. 3 produces a vorticity equation (10), which now yields with for this planar Rayleigh-Taylor flow. As in Forbes [7] we will ignore the boundary layers on the upper and lower walls, for simplicity, and thus assume that the fluid, although viscous, is nevertheless effectively "slipping" along these two surfaces. We seek bilaterally symmetric solutions in the semi-analytic (spectral) representation for the streamfunction, in which it is convenient to define constants After some experimentation, an appropriate representation for the density perturbation function has been found to be using the same constants β n defined in (19). The two integers M and N in the upper limits of the summations should be chosen as large as practicable. With this representation, the single component of the vorticity takes the form with extra constants 2 m,n = m 2 + β 2 n .
The two components u and v of the velocity are calculated from (13) by direct differentiation of (18). Differential equations are now derived for the Fourier coefficients C m,n (t) in (20), by multiplying the transport equation (14) by basis functions cos(kx) cos (β (y + h 1 )), for integers k = 0, 1, 2, . . . , M and = 0, 1, 2, . . . , N , and integrating over the rectangular region −π < x < π, −h 1 < y < h 2 . After a little algebra this yields Here, the symbol δ k,0 is the Kronecker delta, and is zero unless k = 0 when it takes the value one. Similarly, δ 0, is zero unless = 0 when it becomes one.
For the Boussinesq model in Sect. 2, the vorticity equation (15) is similarly multiplied by basis functions sin(kx) sin (β (y + h 1 )) and integrated over the domain. This gives This is a system of differential equations for the coefficients A m,n (t). For the more complicated Extended Boussinesq model in Sect. 3, the appropriate vorticity equation is (16), and when subject to Fourier analysis it leads to the similar system of differential equations. Here, the function H E is as defined in (17). The numerical solution of the classical Boussinesq model of Sect. 2 therefore consists of integrating the two sets of ordinary differential equations (23) and (24) forward in time, to create the required coefficients C m,n and A m,n . This is a total of (M + 1)(N + 1) + M N nonlinear equations, and we use the Runge-Kutta-Fehlberg differential equation package ode45 from Matlab for this purpose. The integrals appearing in these equations are evaluated using the Gaussian quadrature routine lgwt written by von Winckel [23]. With more than about 200 numerical points, this process calculates these integrals essentially exactly. If, on the other hand, the Extended Boussinesq model of Sect. 3 is used, the systems of ordinary differential equations (23), (25) are integrated forward numerically. Once coefficients C m,n and A m,n have thus been obtained, in either model, the solution variables ρ and are then recreated from (20) and (18).
We impose initial conditions (12) directly on the series representation (20), and then obtain initial values for the coefficients C m,n using Fourier analysis, as described above in the derivation of the system (23). The result is similar to that described in Forbes [7]. After some algebra, and using formulae from Abramowitz and Stegun [24, p. 360], we obtain and In addition, In these expressions, the symbol J ν denotes the first-kind Bessel function of order ν. These coefficients (26)-(29) are obtained by Fourier analysis of the discontinuous initial condition (12), and so the function reconstructed with these coefficients will replace the exact discontinuity with oscillations either side of it, according to the Gibbs phenomenon [25]. Accordingly, we smooth the initial profile using Lanczos smoothing, in which the coefficient C m,n is multiplied by the factor sin (nσ L ) / (nσ L ), n = 1, 2, . . . , N . This has the effect of replacing the original discontinuous function (12) with a continuous function that possesses smooth but rapidly-varying sections at the sites of the original discontinuities. Furthermore, the Lanczos parameter σ L can be chosen to eliminate nearly all of the undesirable Gibbs oscillations in the profile, and we have typically chosen σ L = 0.05. Figure 2 illustrates the profile (12) for the density perturbation ρ at the initial time t = 0, after reconstruction from the Fourier representation (20) with coefficients (26)-(29), followed by Lanczos smoothing. Here, the density ratio is D = 10. The channel walls are placed at y = −h 1 = −2 and y = h 2 = 2 and the initial amplitude of the perturbation to the interface is = 0.03. It is evident from Fig. 2 that the discontinuity at the original interface has been approximated by a smooth, continuous, narrow segment across which ρ changes rapidly from zero to (D − 1) = 9. Some very small oscillations related to Gibbs' phenomenon might be visible near the interfacial region; they could be reduced still further with slightly more vigorous damping (increasing the Lanczos parameter σ L ), although at the cost of widening the interfacial zone further. This diagram therefore illustrates the effects of smoothing on the exact discontinuity (12).

Results for symmetric profiles
We begin this presentation of numerical results by considering first the case of symmetric solutions, based around the use of the spectral solution representations (18)- (22) in Sect. 4. The key aim here is to compare the behaviour of solutions to the unstable Rayleigh-Taylor outflow, as predicted by the classical Boussinesq solution of Sect. 2 and the Extended Boussinesq model of Sect. 3.
We consider two sets of solutions in Fig. 3, both obtained with the same parameter values. These figures display contours of the density perturbation function ρ; however, the results in part (a) were generated using the classical The scales on the several axes in Fig. 3 are equal, so that the solution profiles are as they would actually appear. To remove spurious effects due to the horizontal bounding walls at y = −h 1 and y = h 2 , we have chosen large values h 1 = h 2 = 10, so as to position these walls away from the region of interest near the evolving plume. The density ratio has the small value D = 1.05, for which the pure Boussinesq model would be expected to be accurate. The Reynolds number is R e = 10 4 , the diffusion coefficient in (11) has the value σ = 10 −4 and the initial amplitude of the disturbance to the interface was = 0.03. The two sets of solutions in parts (a) and (b) were obtained at the same four (dimensionless) times t = 40, 52, 60 and 68.
There is clearly very close agreement between the predictions of these two models, for this small density ratio D = 1.05, as is to be expected. The results are almost identical for the first two time frames, at t = 40 and t = 52, and they both show a bubble rising and turning over to form a mushroom-shaped structure. In addition, the plumes are top-bottom symmetric, in the sense that the rising bubble has the same shape as the downward-moving spike. Even for the later two times t = 60 and t = 68 shown, the two sets of solutions obtained with the two different models are extremely similar, with only minor differences between them. This close agreement shows the robustness of predictions made at low values of the density ratio D.
Some comment on the factors affecting the numerical accuracy of the solutions in Fig. 3 is required. Within the physical assumptions of each of the models described, there are three potential sources of numerical error. The first is that, although we are using spectral representations (18) and (20) that are capable of arbitrarily high precision, we are in practice restricted to some maximum numbers M, N of Fourier modes that can be employed. This represents the most severe limitation on the accuracy of our scheme. A second source of error involves the  (23) and (24) or (25) can be performed. Here, however, these integrations are carried out essentially exactly, due to the large numbers of mesh points and the highly accurate Gauss-Legendre quadrature algorithm [23] that are used. The third potential source of error arises when the initial density profile is allowed to be discontinuous across an effective interface, as in the condition (12) illustrated here. In that case, it is well known that a Fourier series of smooth basis functions is not able to represent the discontinuity exactly, but instead, converges at the location of the discontinuity itself to the average value (D − 1)/2, with undulations either side, according to Gibbs' phenomenon [25]. This effect can be mitigated by replacing the exact discontinuous condition (12) with one in which the density varies smoothly, and as discussed in relation to Fig. 2, this has been done here using Lanczos smoothing. Nevertheless, if the initial density gradient is large, some remnant of the small oscillations due to the Gibbs phenomenon persist, although their amplitude decreases with increasing numbers M, N of Fourier modes.
These features of the numerical solution are illustrated in Fig. 4. These are both solutions computed using the Extended Boussinesq Theory at the particular time t = 60 and with density ratio D = 1.05, as in Fig. 3b. The diagram on the left has been obtained using (M, N ) = (51, 101) Fourier modes, whereas (M, N ) = (71, 141) modes were used in generating the figure on the right. This second diagram required a very substantial allocation of computing resources to create, although it clearly agrees very closely with the less accurate diagram on the left, even for a solution such as this, in which the interfacial zone involves the formation of elaborate and complex structures. Most apparent, however, is the fact that the more accurate picture on the right has its density contours focussed far more closely around the interfacial zone, with almost no contour lines further away. For the less accurate case (M, N ) = (51, 101) on the left, the very small-amplitude oscillations in the initial profile, associated with the smoothed discontinuity and illustrated in Fig. 2, diffuse throughout the solution domain, and the many such contour lines are evident still at time t = 60 in the picture on the left of Fig. 4. When the number of modes is increased to (M, N ) = (71, 141), the corresponding oscillations in the initial density profile are even smaller, and are damped more rapidly by viscosity by time t = 60, as can be seen in the right-hand picture. Contouring software can also appear to exaggerate the importance of very small-amplitude undulations, as is perhaps evident in the initial profile illustrated in Fig. 2.
In strong contrast to Fig. 3, the results of applying the two different models to a situation in which the density ratio D = 10 is large are shown in Fig. 5. The results in part (a) are again obtained using the classical Boussinesq theory of Sect. 2, and they show florid plume growth with extensive overturning regions, to such an extent that the plume cap has evidently detached at the later times shown, forming a separate rising bubble of lighter fluid. However, the predictions of the Extended Boussinesq model of Sect. 3 in part (b), at this high density ratio, indicate that the heavier upper fluid has suppressed the growth of the upwardly-moving bubble, which only begins to overturn at about dimensionless time t = 5. Of course the classical Boussinesq theory of Sect. 2 has been derived on the assumption that the two fluids have comparable densities, and so that theory cannot be expected to be valid for this density ratio D = 10. It was found that, for this large density ratio, the Extended Boussinesq model in Sect. 3  produces overturning regions at the interfacial zone on extremely fine length scales, that appear to correlate roughly with O (1/R e ). Consequently, when the same parameter choice 1/R e = 10 −4 , σ = 10 −4 was made (as in Fig. 3), the interface developed a small region of extremely high curvature, reminiscent of the Moore curvature singularity [6] of inviscid free-surface flow, and the numerical time integration of the differential equations (23), (25) fails at some finite time, associated with the formation of this region. For this reason, we have decreased the Reynolds number in Fig. 5 to the value R e = 10 2 , so that the numerical solution was capable of resolving the overturning sections of the interface. Contours of the streamfunction in (18) are displayed in Fig. 6, for the same solutions as shown in Fig. 5, and at the same four times. This set of diagrams highlights the strong differences in flow patterns produced by the pure Boussinesq theory in part (a), and the Extended Boussinesq theory shown in part (b). The Extended theory produces two bilaterally symmetric vortices near the interface in part (b), which, as time progresses, become more intense and rotate from a somewhat horizontal orientation to a more vertical alignment as the interface begins to overturn. By contrast, the Boussinesq results in part (a) indicate that more and more vortices are generated as time progresses, associated with the prediction by that theory of the elaborate plume structure displayed in Fig. 5a.
In Fig. 7, another set of solutions is displayed, again for the large density ratio D = 10, and with h 1 = h 2 = 10. Here, however, we have increased the Reynolds number to R e = 10 4 , but have also altered the density diffusion parameter to the value σ = 10 −2 . As discussed above, the use of larger Reynolds numbers in the Extended Boussinesq theory can produce very small regions of very large curvature in the interfacial zone when D is large, mimicking the Moore curvature singularity [6] of inviscid theory, even although viscosity is present [13], and this is sufficient to overwhelm the accuracy requirements of a numerical solution technique. In order, therefore, to allow the  numerical scheme to continue to the longer times shown in Fig. 7, it was necessary to use the larger value σ = 10 −2 in these computations. From the mathematical point of view, the two parameters 1/R e and σ act to regularize the solution, by smoothing out features at very small length scales; however, from the physical viewpoint, the use of the larger Reynolds number in Fig. 7 (as opposed to the smaller value in Figs. 5 and 6) is possibly more appropriate in many circumstances, although the larger value of σ may be more difficult to justify. Figure 7 again shows how, at larger density ratio D, the classical Boussinesq theory of Sect. 2 predicts the rapid development of an overturning plume. The results for that case, shown in part (a) are not significantly affected by the choice of Reynolds number, and the contours of density perturbation function ρ shown at time t = 4 in Fig. 7a are almost identical to those at the same time in Fig. 5a. Mixing occurs rapidly in the classical Boussinesq model and there is almost complete interchange of the two fluids by time t = 8 in Fig. 7a; mixing, too, occurs on progressively smaller length scales, and we found it was no longer practical to continue the computations much beyond t = 8, since smaller and smaller time intervals were needed in the calculations, resulting in excessively long computational times.
By contrast, the Extended Boussinesq theory in Fig. 7b illustrates the way in which the heavier upper fluid restrains the growth of the upwardly-moving bubble of lighter fluid. This is perhaps most easily seen at time t = 6. However, the Boussinesq theory in part (a) predicts an elaborately-shaped plume that has almost risen to the top wall at y = h 2 = 10, with spikes that have likewise descended almost to the bottom wall at y = −h 1 = −10, the Extended Boussinesq theory in part (b) suggests that the interface has only just begun to overturn at this same time t = 6.

Comparison with smoothed-particle hydrodynamics
For the case of Rayleigh-Taylor instability, it has been observed by a number of authors that the classical Boussinesq model tends to exaggerate the rapidity of interface overturning and the rate of mixing [14][15][16][17] when the density ratio D is large. This is discussed in Sect. 1, and is not surprising since the classical approximation is based on the assumption D ≈ 1. Our results have indeed confirmed that the Boussinesq and the Extended Boussinesq theories agree very closely for small density ratios, such as the results presented in Fig. 3 with D = 1.05, but disagree quite markedly at larger density ratios, as illustrated in Fig. 5 for D = 10. It remains, however, to confirm that the Extended Boussinesq model is indeed producing results that are consistent with a full hydrodynamical theory.
In this section, we attempt to verify the Extended Boussinesq model by comparing its results with those generated by a smoothed-particle hydrodynamical (SPH) simulation. In this approach, the smoothed representation of some scalar attribute of the fluid, A, at the vector position r is given by the interpolation integral (see for example [26]) where d 2 r is the differential for an area integral, and W (..., ε) is a Gaussian kernel function of smoothing length ε (see [27]). The smoothed representation of some vector attribute, v, is similarly defined. The integral is approximated using a Riemann sum over the discrete locations of the particles (that is, r i for i = 1, . . . n), and analytically differentiated to yield expressions for the gradient of scalar attributes, ∇ A(r), and divergence of vector attributes, ∇ · v(r). The value of a scalar or vector attribute at the location of the i-th fluid particle is found by evaluating the relevant equation at r = r i for each particle. The motion of each fluid particle is described by equations for continuity, conservation of momentum and the equation of state for water-like fluids [28], expressed in terms of the weighted Riemann sum over the other particle properties [29,30].
We use the PyPSH framework (see [31]) to manage efficiently the particles and numerical integration for the smoothed-particle hydrodynamical simulation in this work. Their default second-order predictor-corrector integrators are used to update the density, velocity and position based on the time-derivatives in the continuity equation and equation of motion (the position is updated using the velocity). The pressure is then calculated using the updated value for the density and the equation of state.
The simulations in this approach involve three types of particles: (1) a dense, slightly compressible, low-Reynolds number fluid (ρ = D, dimensionless sound speed c 0 = 100, R e = 100); and (2) a lower density, slightly compressible, low-Reynolds number fluid (ρ = 1, c 0 = 100, R e = 100); and (3) boundary particles. The initial conditions for the fluid particles are set to match the semi-analytic methods described above. Meanwhile, the boundary particles numerically enforce the nonpenetration condition q ·n = 0 on the periodic boundaries at x = ±π and walls at y = ±h (see [32]). We first ran a simulation with 640, 000 fluid particles initially spaced uniformly in the Cartesian grid; this number of particles is required to resolve adequately the initial perturbation in the two fluid densities. However, the curved interface is represented as a "grainy" boundary, containing many small discrete steps caused by the presence of the particles themselves. Unfortunately, these small discontinuities serve as seeding points for the growth of Rayleigh-Taylor instabilities, and a similar situation was encountered by Walters and Forbes [22], who found that even the errors produced by single-precision arithmetic could grow into instabilities, and required double-precision accuracy to prevent. We overcame these numerical shortcomings here with the following two modifications: Firstly, we increased the smoothing area of the kernel by a factor of 2.8 from the default values assumed in the technique; that is, we set ε = 2δx, in which δx represents the initial particle spacing. Secondly, we transformed the coordinate system (for the initial particle placement) so that y = 0 lies along the interface at t = 0. This ensures that particles are cleanly allocated between the two fluids. This modified simulation comprises a reduced system of only 80, 100 fluid particles and 9, 300 boundary particles, successfully producing results with very reduced numerical artefacts in a much shorter computational time. Figure 8 shows the results produced by the three different fluid models discussed here, at the five different times t = 1, 2, 3, 4 and 5, for density ratio D = 10 and with initial amplitude = 0.03, as for Fig. 5. Here, however, we have set h 1 = h 2 = 1 so as to reduce the number of SPH fluid particles needed in the simulation, and thus to minimize the computer run-time required. The predictions of the classical Boussinesq model are shown in part (a), and the results obtained using the Extended Boussinesq model are shown in part (b). The third set of results in part (c) were generated using the SPH code, modified as described above.
The Boussinesq theory in part (a) predicts an interfacial zone that grows rapidly and forms florid overturning structures quite early. In that respect, it is the outlier among the three sets of results shown in Fig. 8, although of course, this model assumes a density ratio near one, and so cannot be expected to be accurate with D = 10. Thus, by dimensionless time t = 3, Boussinesq theory has created a large plume of lower density fluid that has nearly reached the upper wall at y = h 2 = 1, and which possesses well-developed rolled-up overhanging portions. The Extended Boussinesq model used in Fig. 8b, by contrast, generates an interfacial zone which, at time t = 3, still retains most of the features of the initial co-sinusoidal profile. This is very similar to the interface shape generated by the SPH computation in part (c), although this has a slightly lower amplitude. As discussed above, the particle interactions in the SPH simulation are modelled using a smoothing kernel, resulting in an effective damping parameter in the method; this will certainly contribute to the lower wave amplitudes in Fig. 8c. Overall, however, there is significant qualitative agreement between the results of the Extended Boussinesq theory in part (b) and the SPH results in part (c).
The predictions of the three different models are shown in slightly greater detail in Fig. 9, for the same parameter values as Fig. 8, and focussing on the last two times t = 4 and t = 5, in parts (a) and (b), respectively. At time t = 4 in part (a), the Boussinesq theory has produced an elaborate plume that has hit the upper wall, and by time t = 5 in part (b) it has mixed substantially with the ambient fluid, resulting in a rather indistinct structure. The predictions of the Extended Boussinesq model and the SPH calculations are of more interest, however, and while they broadly agree, there are nevertheless some interesting differences between them. These can be seen most clearly at the later time t = 5 shown in Fig. 9b. The Extended Boussinesq theory has produced a plume with some degree of topbottom asymmetry, whereas the SPH results predict a plume with a broad flat (bubble) top but much more sharply Fig. 8 Calculation of a Rayleigh-Taylor instability at five different times t = 1, 2, 3, 4 and 5 using the three fluid models, for large density ratio D = 10 and initial perturbation amplitude = 0.03. The Reynolds number is R e = 10 2 , and h 1 = h 2 = 1. Maps of the density perturbation ρ are shown for each model. a Boussinesq model; b Extended Boussinesq model; and c SPH model spiked downward-facing spikes of the heavier fluid. This is consistent with the results for high density obtained by Yakovenko [16], which perhaps suggests that the hydrostatic approximation used in one term of the Extended Boussinesq theory might still need further refinement in some applications. The Extended Boussinesq theory also predicts a very narrow overturning portion at the interface for t = 5; this is not present in the SPH solution, although the effects of damping may possibly affect this prediction on such a small length scale. Nevertheless, Fig. 9 gives Fig. 9 Detailed comparison of the three models at two times, with density ratio D = 10, channel heights h 1 = h 2 = 1 and starting from initial perturbation amplitude = 0.03. The Reynolds number is R e = 10 2 . Maps of the density perturbation ρ are shown for each model. Results are presented at a time t = 4 and b time t = 5 some confidence that the Extended Boussinesq theory is capable of giving acceptable results for large density ratio D, whereas the classical Boussinesq theory clearly cannot.

Asymmetric profiles
So far, it has been assumed that the Rayleigh-Taylor plumes are bilaterally symmetric, but of course this may not necessarily be the case. Any sort of disturbance to the background fluid medium is capable of distorting the plume and giving rise to an asymmetric RT plume.
We investigate this briefly here, by generalizing the spectral forms of our solutions (19) and (20) for the RT instability, to become and ρ (x, y, t) = M m=0 N n=0 C m,n (t) cos (α m (x + π )) cos (β n (y + h 1 )) . Here, it is convenient to define constants α m = m/2 and β n are the same as defined in (19). These new forms (30), (31) are substituted into the governing equations and Fourier analysed, as before. For the initial conditions (12), the spectral representation (31) is again used to generate initial forms C m,n (0) for the coefficients, to start the calculation. In Fig. 10, we show a sample result for an asymmetric RT instability computed using the Extended Boussinesq theory, with a large density ratio D = 10. Here, the upper and lower walls are placed at h 1 = h 2 = 10 and the Reynolds number and diffusion coefficients are chosen to be R e = 10 4 and σ = 10 −2 , as in Fig. 7. The solution was run forward in time until it failed at dimensionless time a little larger than about t = 9.5, due to interfacial structures that develop on extremely small length scales. However, it is accurate and well converged for earlier times, and four such profiles are displayed in Fig. 10 at times t = 2, 4, 6 and 8. The disturbed background flow in this example was created simply by setting A 1,1 (0) = V 0 in (30) to create an initial velocity that is periodic in x, and the results in Fig. 10 used V 0 = 0.1.
Density profiles are displayed in part (a) at these four times. Under the effect of the initial background flow, the initially symmetric interfacial profile begins to distort towards the left of the diagram, and this is already evident by time t = 4. Eventually a plume with overturning segments is formed, although it is now strongly asymmetric, as can be seen at the last time t = 8 shown. In addition, small wavelets begin to appear on the upper surface of the stem of the plume, and these appear to be related to Kelvin-Helmholtz instabilities, since they continue to grow and eventually make the calculations fail at about t = 9.5.
Contours of the streamfunction are shown in part (b) at these same four times. While these do not necessarily correspond exactly to streamlines, for a time-dependent flow such as this, they nevertheless give a reasonable overview of the flow, and the first set of streamlines at time t = 2 indicate how the initial choice A 1,1 (0) = V 0 has produced a type of vortex in each wavelength. As time progresses and the interface distorts to the left, a secondary vortex begins to appear in the upper right-hand side of each diagram, associated with the overturning interface in part (a).

Summary and conclusion
This paper has assessed the classical Boussinesq theory as an approximate method for computing viscous freeboundary flows. The model has the enormous advantage that it eliminates the need to cope with an exact interface, upon which the correct interfacial boundary conditions given in Batchelor [18, p. 150] are imposed. Instead, a multi-fluid system is approximated as a single fluid, within which the density changes smoothly but rapidly in interfacial zones. This is equivalent to treating the different fluid densities as having been caused by the presence of a heavy contaminant dissolved in the fluid, as discussed in Sect. 2.
Boussinesq theory is derived on the assumption that the density differences between different fluids in the system remain small. Unsurprisingly, we have found that the theory performs very well when precisely this is the case, but that it becomes increasingly inappropriate as the density ratio D between different fluids in the system becomes larger. To overcome this, while still maintaining the Boussinesq philosophy of treating the system as a single fluid of varying density, so as to simplify calculations enormously, we have proposed an Extended Boussinesq theory that avoids the severity of assuming the density ratio D is approximately one; however, it does approximate the pressure gradient using a fairly simplistic hydrostatic law (involving the variable density of the approximate single fluid), and further refinements to that assumption may perhaps be possible in the future. The Extended Boussinesq theory, when used at high density ratios, gives results more consistent with those obtained by other authors [14][15][16][17] using a variety of different methods, and avoids the exaggerated interfacial profiles that Boussinesq theory would predict. We have been able to give broad confirmation of the results produced by the Extended Boussinesq model by comparison with the results of an SPH simulation in Sect. 5.2, although this comparison has nevertheless revealed that the Extended Boussinesq theory may still not produce the degree of top-bottom asymmetry expected in the interface profiles.
We have also demonstrated the use of Extended Boussinesq theory for highly nonsymmetric unstable interfaces in Fig. 10. In addition, different initial conditions have also been investigated, although these are not shown here in the interest of space. The Extended Boussinesq approach developed here should be capable of being developed further for use in magnetohydrodynamic situations in astrophysics, where a background magnetic field can drive asymmetries in the plume, similar to those shown in Fig. 10, and examples of this effect have been presented by Chambers and Forbes [33].