A model for Rayleigh-Taylor mixing and interface turn-over

We first develop a new mathematical model for two-fluid interface motion, subjected to the Rayleigh-Taylor (RT) instability in two-dimensional fluid flow, which in its simplest form, is given by $ h_{tt}(\alpha,t) = A g\, \Lambda h - \frac{\sigma}{\rho^++\rho^-} \Lambda^3 h - A \partial_\alpha(H h_t h_t) $, where $\Lambda = H \partial_ \alpha $ and $H$ denotes the Hilbert transform. In this so-called $h$-model, $A$ is the Atwood number, $g$ is the acceleration, $ \sigma $ is surface tension, and $\rho^\pm$ denotes the densities of the two fluids. Under a certain stability condition, we prove that this so-called $h$-model is both locally and globally well-posed. Numerical simulations of the $h$-model show that the interface can quickly grow due to nonlinearity, and then stabilize when the lighter fluid is on top of the heavier fluid and acceleration is directed downward. In the unstable case of a heavier fluid being supported by the lighter fluid, we find good agreement for the growth of the mixing layer with experimental data in the"rocket rig"experiment of Read of Youngs. We then derive an RT interface model with a general parameterization $z(\alpha,t)$ such that $ z_{tt}= \Lambda\bigg{[}\frac{A}{|\partial_\alpha z|^2}H\left(z_t\cdot (\partial_\alpha z)^\perp H(z_t\cdot (\partial_\alpha z)^\perp)\right) + A g z_2 \bigg{]} \frac{(\partial_\alpha z)^\perp}{|\partial_\alpha z|^2} +z_t\cdot (\partial_\alpha z)^\perp\left(\frac{(\partial_\alpha z_t)^\perp}{|\partial_\alpha z|^2}-\frac{(\partial_\alpha z)^\perp 2(\partial_\alpha z\cdot \partial_\alpha z_t)}{|\partial_\alpha z|^4}\right)$. This more general RT $z$-model allows for interface turn-over. Numerical simulations of the $z$-model show an even better agreement with the predicted mixing layer growth for the"rocket rig"experiment.


Introduction
The instability of a heavy fluid layer supported by a light one is generally known as Rayleigh-Taylor (RT) instability (see Rayleigh [7] and Taylor [12]). It can occur under gravity and, equivalently, under an acceleration of the fluid system in the direction toward the denser fluid; in particular, RT is an interface fingering instability, which occurs when a perturbed interface, between two fluids of different density, is subjected to a normal pressure gradient. Whenever the pressure is higher in the lighter fluid, the differential acceleration causes the two fluids to mix. See Sharp [11], Youngs [13,14], and Kull [6] for an overview of the RT instability.
The Euler equations of inviscid hydrodynamics serve as the basic mathematical model for RT instability and mixing between two fluids. This highly unstable system of conservation laws is both difficult to analyze (as it is ill-posed in the absence of surface tension and viscosity) and difficult to computationally simulate at the small spatial scales of RT mixing. As such, our objective is to develop model equations, which can be used to predict the RT mixing layer and growth rate.
In order to derive our RT interface models, we shall assume both incompressible and irrotational flow for the two-fluid Euler equations. Rather than proceeding with a direct approximation of the Euler equations, we shall instead work with the equivalent Birkhoff-Rott singular integral-kernel formulation for the evolution of the material interface. We have found that this formulation possesses a certain robustness in regards to approximations founded upon expansions in various parameters.
In the simplest case, in which the interface is modeled as a graph (α, h(α, t)) of a signed height function h(α, t), Our approach yields a second-order in-time quadratically nonlinear wave equation for the position of the interface, the h-model, which is given by where Λ = H∂ α and H denotes the Hilbert transform. In this h-model equation, A denotes the Atwood number, g is the acceleration, σ ≥ 0 is the surface tension, and ρ ± > 0 denotes the densities of the two fluids.
As we will describe below, our h-model equation for RT instability is both locally and globally well-posed in Sobolev spaces when a stability condition is satisfied, which requires the product of the Atwood number and the initial velocity field to be positive. Under such a stability condition, we also derive the asymptotic behavior of solutions to the h-model as t → ∞. A number of interesting energy laws are also derived, which distinguish between stable and unstable interface dynamics Numerical simulations are performed, which show that the h-model is capable of producing remarkable growth of the interface, followed by (possibly oscillatory) decay to a rest state in the case that the lighter fluid is supported by the heavier fluid. In the highly unstable case, where a heavy fluid is supported by the lighter fluid, we perform the classical "rocket rig" experiment of Read [8] and Youngs [14] for the case of unstable RT mixing-layer growth, and find very good agreement for the growth rates with the predicted quadratic growth rate of Youngs [14], and with both Direct Numerical Simulations and experimental data.
Finally, in order to allow for the interface to turn-over (rather than simply remaining a graph) we develop a more general z-model for the interface parameterization z(α, t) = (z 1 (α, t), z 2 (α, t)) which satisfies Rather than constraining the interface amplitude to grow at the RT instability, the z-model can allow for interface turn-over. We perform numerical simulations to demonstrate the improvement afforded by this more general model, and the even better accuracy in matching the predicted growth of the RT mixing layer for the "rocket rig" experiment. We also perform the so-called "tilted rig" experiment, in which the tank holding the fluid is titled by a small angle away from vertical. Again, our simulations demonstrate good qualitative agreement with DNS. To conclude, we also perform a numerical simulation for a Kelvin-Helmholtz problem, for which the Atwood number is set to zero to prevent an RT instability; starting from a steep mode-1 wave, we show roll-over of the interface with a very localized energy spectrum.

Equations for multi-phase flow
The Euler equations are a system of conservation laws, modeling the dynamics of multiphase inviscid fluid flow. For incompressible two-dimensional motion, the conservation of momentum for an homogeneous, inviscid fluid can be written as where u = (u 1 , u 2 ) denotes the velocity vector field of the fluid, p is the scalar pressure function, ρ is the density, g is the acceleration, ∇ = (∂ x 1 , ∂ x 2 ) is the gradient vector, and e 2 = (0, 1). For incompressible flow, the conservation of mass is given by We shall further assume that the fluid is irrotational so that ω := curl u = 0, where ω is the vorticity of the fluid, and curl u = ∂ 1 u 2 − ∂ 2 u 1 . We assume that there are two fluids with densities ρ + and ρ − , separated by a material interface Γ(t), where t denotes an instant of time in the interval [0, T ]. As shown in Figure  1, the fluid with density ρ + lies above Γ(t), while the fluid with density ρ − lies below Γ(t). The fluid on top of Γ(t) has density ρ + , while the fluid on the bottom has density ρ − .
For each instant of time t ∈ [0, T ], the interface Γ(t) is parameterized by a function z(α, t). We will provide a special choice for the parameterization z(α, t) below.
We denote the jump of a field variable f (x, t) across Γ(t) by We let n and τ denote the unit normal and tangent vectors to Γ(t), respectively; for the RT instability, we have the following jump conditions: Note, that ∂ α z(α, t) is a tangent vector to Γ(t) at the point z(α, t), so that Consequently, the velocity is not continuous at the interface. Due to the fact that the shape of the interface is determined only by the normal component of the fluid velocity, the motion of the interface has a tangential reparameterization symmetry, so that we can add any tangential velocity to the motion of the interface, with the hope that this will simplify the analysis. As such, the evolution equation for the parameterization of Γ(t) is written as where the function c(α, t) will be specified below.

The integral kernel formulation
A number of modal models have been proposed for the evolution of the RT mixing layer; see, for example, Rollin & Andrews [10] and Goncharov [3], and the references therein. These modal models are based on a modal decomposition and approximation of the evolution equations for the parameterization of the interface and the velocity potential; such models consist of a large coupled system of nonlinear ODEs for a finite set of Fourier modes. Rather than developing a modal model for RT mixing and approximating the partial differential equations themselves, we shall take another approach to the development of an RT model, which is founded upon the Birkhoff-Rott integral-kernel formulation for the evolution of an interface, separating two incompressible and irrotational fluids.
In order to introduce the integral-kernel formulation, which consists of singular integrals, we begin by defining the principal value integral of a given function f as The well known Biot-Savart kernel K BS is an integral representation for ∇ ⊥ ∆ −1 where ∇ ⊥ = (−∂ x 2 , ∂ x 1 ), and ∆ = ∂ 2 x 1 + ∂ 2 x 2 denotes the Laplace operator in the plane. In other words, if the two fluids fill the plane, then ∆ −1 is given by the Newtonian potential and the Biot-Savart kernel is given by Similarly, if the fluid flow is periodic in the horizontal variable, then . (6) Due to the characteristics of the irrotational flow, the vorticity is a measure which is supported on the interface Γ(t), written as where δ Γ(t) is the Dirac delta function supported on the interface Γ(t). More precisely, the vorticity ω is a distribution defined as follows: for all smooth test functions ϕ with compact support, The function ̟ is the amplitude of the vorticity along Γ(t). Notice that ̟ is minus the jump of the velocity in the tangential direction: Given the vorticity measure ̟, we can reconstruct the velocity field everywhere; specifically, we have, thanks to the Biot-Savart law, that Since by (2), the velocity is irrotational, there exist a velocity potential function φ : R 2 → R such that u ± = ∇φ ± . The potential φ satisfies Bernoulli's equation and φ is related to ̟ by the singular integral equation As the kernel K BS (x, z(β, t)) is singular for x = z(α, t), it follows that the tangential component of u is discontinuous across the interface Γ(t). By similar reasoning, the potential function φ is also discontinuous across Γ(t). Moreover, the vorticity ̟ along the interface Γ(t) is related to the jump of the velocity potential by Now, using the parameterization z(α, t) for the interface Γ(t), we define the Birkhoff-Rott principal-value integral, denoting R by , as In particular, (11) is explicitly given by The Birkhoff-Rott function K BR (α, t) denotes the average velocity at a given point α ∈ Γ(t), so that

Dynamics of the interface and vorticity amplitude
We now formulate the dynamics of the interface. Since z t = u + c∂ α z, we see that In order to compute ̟ t , we use the identity (10), and shall need to To do so, we shall first compute compute [ [φ] ] across the interface Γ(t). In fact, as we will explain below, it is easier to compute ∂ α [ [φ] ], and to do so, we study the limit of ∂ α φ ± (y, t) as the point y approaches Γ(t) in a direction normal to Γ(t). We recall that ∂ α z(α, t) is a tangent vector to Γ(t) at the point z(α, t) and hence ∂ ⊥ α z(α, t) is a normal vector to Γ(t) at the point z(α, t). For 0 < ǫ ≪ 1, and each time t, we let denote a sequence of points converging to z(α, t) as ǫ → 0 in the normal direction to Γ(t).
Then, the trace of ∂ α φ ± on Γ(t) is defined as the limit as ǫ → 0 of the sequence ∂ α φ ± (y ± ǫ , t). Since u ± = ∇φ ± , the chain-rule shows that In [1], it was shown that Thus, from (14) and (15), it follows that Next, we observe that K BR · ∂ α z can be written as an exact derivative due to the identity which means that due to (17), the identity (16) can be integrated, and We define the Atwood number Using (18) together with Bernoulli's equation (8), we can derive the equation , and using (3), (10) and (13), a lengthy computation shows that

The case that the interface Γ(t) is a graph
We now assume that the interface Γ(t) is the graph of the signed height function h(α, t), so that the parameterization z(α, t) is given by If at time t = 0, Γ(0) is given as the graph (α, h(α, 0)), then we can ensure that Γ(t) stays a graph for future time t > 0 by making an explicit choice of the function c(α, t) in (4). To this end, we define With Γ(t) given by (α, h(α, t), the definition of the Birkhoff-Rott function K BR (α, t) in (11) simplifies to the evolution equation for Γ(t) can be written in terms of the height function h(α, t) as and the evolution equation for vorticity ̟ on Γ(t) takes the form

A simple model equation for RT instability
In this section, we shall derive a model of RT instability in which the interface Γ(t) is a graph of the height function h(α, t). In order to proceed with our derivation, we shall make further approximations to the coupled system of equations (23) and (24).

The equations in the linear regime
In the linear regime, equations (23) and (24) reduce to the following coupled system: where H denotes the Hilbert transform, defined as

The nonlinear regime and the model equation
Having found the linear dynamics, we turn our attention to the the nonlinear regime. Our objective is to derive a model equation which contains the quadratic nonlinearities. To do so, we shall introduce some further notation. We define the operator Λ by Λ̟ = H∂ α ̟ , and the space-integrated vorticity as We shall use the notation f (t) to denote f (β, t)dβ for any integrable function f (β, t).
We make use of the following power series expansions for |ζ| < 1: and ζ 1 + ζ 2 = ζ − ζ 3 + ζ 5 − ζ 7 + · · · , Using these identities together with the approximation the nonlocal terms in equations (23) and (24) can be approximated as If we neglect terms of order O(h 2 ) in equations (23) and (24), we find that Integrating in space, we obtain that the vorticity average ̟ verifies Finally, we can use the Tricomi relation for the Hilbert transform, to obtain The coupled first-order system (29) can be written as one second-order equation for the evolution of the height function: The model equation (30) contains both quadratic and cubic nonlinearities, but we can simply further.
Keeping only the quadratic nonlinearities, we find that the graph of the interface (x, h(x, t)) evolves according to By assuming that our initial vorticity has zero average, we arrive at the nonlinear equa- Equation (32) is supplemented with initial conditions. In particular, we specify the initial interface position and velocity, respectively, by As we explain in the next section, this model equation has a natural stability condition, requiring the product of the Atwood number A and and the initial velocity field h 1 to be positive.
The second-order-in-time nonlinear wave equation (32) is a new model for the motion of an interface under the influence of RT instability. We shall refer to either the system (29) or the wave equation (32) as the h-model.

Well-posedness of the h-model equation (32)
We now prove that our h-model equation (32) is well-posed in Sobolev spaces, when a certain stability condition is satisfied by the data.
The space L 2 (T) consists of the Lebesgue measurable functions on the circle T which are square-integrable with norm h 2 0 = T |h| 2 dx. For s ≥ 0, we define the homogeneous Sobolev spaceḢ For These identities, together with the Poincaré inequality, show that (33) is an equivalent H s (T)-norm. (Recall that we are using the notation f (t) to denote f (β, t)dβ for any integrable function f (β, t).) Theorem 1 (Local well-posedness for the h-model (32)). Let σ ≥ 0, ρ + , ρ − > 0, g = 0 be fixed constants and let (h 0 , h 1 ) denote initial position and velocity pair for equation (32).
then there exists a time 0 < T * (h 0 , h 1 ) ≤ ∞ and a unique classical solution of (32) satisfying Proof.
Step 1. Approximate problem for h ǫ , ǫ > 0. For ǫ > 0, we introduce a sequence of approximations to equation (32). We let P ǫ denote the projection operator in L 2 (T), given by wheref (k) denotes the kth Fourier mode of f . Then, we let h ǫ be a solution to with initial data given by (P ǫ h 0 , P ǫ h 1 ). The projection operator P ǫ commutes with ∂ α and with H (and hence with Λ).
We let the parameter ǫ range in the interval (0, ǫ 0 ] for a constant ǫ 0 ≪ 1 to be chosen later. Then, ODE theory provides a unique short-time solution h ǫ to equation (36) on a time interval [0, T ǫ ]; the solution h ǫ is smooth and can be taken in According to the equation (36), , the fundamental theorem of calculus shows that h ǫ = P ǫ h ǫ and h ǫ t = P ǫ h ǫ t . As such, (36) can be written as Step 2. The higher-order energy norm. We define the higher-order energy norm E(t) by Step 3. Lower bound for h ǫ t . We shall assume that by choosing T ǫ and ǫ 0 sufficiently small, Below, we shall verify that this assumption holds on a time interval [0, T ], where T is independent of ǫ.
Step 4. ǫ-independent energy estimates. We now establish a time of existence and bounds for h ǫ which are independent of ǫ. We multiply equation (32) by Λ 4 h ǫ t , integrate over T, and find that 1 2 where We first write the integral I as The integral I 2 has an exact derivative, so upon integration-by-parts, the last inequality following from the Sobolev embedding theorem. The integral I 3 clearly has the same bound.
With the assumption (38), the integral I 1 provides extra regularity for h ǫ t as follows: The last equality follows from the fact that With our assumed lower bound for h ǫ t in (38), we find that Using the computation 2.5 , and hence Thus, The fundamental theorem of calculus shows that and hence It follows that By Gronwall's inequality, Step 5. Verifying the lower-bound on h ǫ t . As the bound (43) for E(t) is independent of ǫ and σ, there exists a time interval [0, T ], with T independent of ǫ and σ, such that Using the equation (36'), we see that h ǫ tt is bounded in L 2 (0, T ; L 2 (T)). Thus, we can interpolate between the last two inequalities and use the Sobolev embedding theorem to conclude that By choosing T and ǫ 0 sufficiently small, and using (34), we verify (38).
Step 6. Existence of solutions. From (44), for all 1 < p < ∞, By using Rellich's theorem, we can pass to the limit as ǫ → 0 in (36'). Since f L ∞ = lim p→∞ f L p , we find that the limit h is a solution of (32) on [0, T ] and It is easy to prove that t → h(·, t) and t → h t (·, t) are continuous into H 2.5 (T) and H 2 (T), respectively, with respect to the weak topologies on H 2.5 (T) and H 2 (T); furthermore, we can again find a differential inequality that is almost identical to (42) (but this time for the solution h rather than the sequence h ǫ ). It follows that t → h(t) 2.5 and t → h t (t) 2.5 are continuous, hence When σ > 0, by the same argument, we have the better regularity Step 7. Uniqueness of solutions. Uniqueness of solutions follows from a standard L 2 -type energy estimate, and we omit the details.
Remark 1. If the stability condition A h 1 > 0 is not satisfied, the evolution equation (32) may not be well-posed in Sobolev spaces. For analytic initial data, however, the equation does not require a stability condition for a short-time existence theorem. We shall investigate this further in future work.
Having established existence and uniqueness of classical solutions to the RT model (32), we next establish a number of interesting energy laws satisfied by its solutions.
Proposition 2 (Energy laws for solutions to the h-model (32)). Given constants σ ≥ 0, ρ + , ρ − > 0, and g = 0, suppose that h is a solution to (32) with initial data (h 0 , h 1 ). Then h verifies the following energy laws: where the terms D 1 , D 2 , and D 3 are given by Proof. The proof reduces to a careful manipulation of the following integrals: where I j = −D j /2, for j = 1, 2, 3. First, Next, Using the skew-adjointness of the Hilbert transform and Tricomi relation (28), Finally, where we have used the fact that Λ −1 ∂ α = −H, so that H∂ α = Λ.
Corollary 3. If ρ + < ρ − so that the Atwood number A < 0, and if gravity acts downward so that g > 0, then whenever the initial velocity h 1 satisfies the stability condition h 1 < 0, then the energy law (45b) shows that Remark 2. The energy law (45a) provides decay for lower-order norms when A h 1 > 0, while the energy law (45c) may actually cause a growth-in-time t which behaves like t 2 .
There may indeed be other higher-order energy laws to the model equation (32) that have yet to be established.
Our analysis will rely on the fact that since we have assumed the initial velocity satisfies h 1 < 0. It is convenient to introduce a new variable f (α, t) given by It follows that f satisfies the evolution equation with initial data The local existence for h and h t (and consequently, for f and f t ) follows from Theorem 1. Thus, in order to prove that solutions exist for all time, it remains only to establish estimates which are uniform in time, and the desired asymptotic behavior will be established by showing that both f and f t converge to zero. With λ defined by (34), the stability condition (46) for h t reduces to Furthermore, we have that We test equation (48) against Λf t and obtain that We shall make use of the refined Carlson inequality: Using (50) together with the Poincaré inequality shows that w L ∞ < w 1 .
As a consequence, we obtain that Taking a time derivative of equation (48), we find that f t satisfies with initial data given by Testing (52) against Λf tt , and using (50), we obtain that (53) Next, we define the following energy and dissipation functions, respectively, as Then (53) can be written as and we find decay of the energy provided that the initial data satisfies Consequently, due to (47), the initial data f 0 and f 1 satisfy and we have that equations (51) and (53) reduce to Thus, we find the estimates and the asymptotic behavior A model of Rayleigh-Taylor mixing Finally, using (50) and (54), we conclude that and the stability condition (49) is satisfied even in a stricter sense. Finally, we test equation (48) against Λ 4 f t . We find that Using integration by parts, Sobolev embedding and interpolation, we have that where we have used the classical commutator estimate [5] Using Young's inequality and (56), we find that for a certain explicit 0 < δ(h 0 , h 1 ). Using Gronwall's inequality, we conclude that Finally, from the regularity of f and f t , we can conclude the regularity of h and h t .
Remark 3. Let us emphasize that the condition (47) does not require small initial data; rather, we require the data to be sufficiently close to an arbitrarily large homogeneous state. For example, in the case that surface tension σ = 0, we can consider  . Let us set surface tension σ = 0, in which case our initial position h 0 and initial velocity h 1 must be chosen so that Thus, a constraint is placed on the H 1.5 -norm of the initial position function h 0 , even though the nonlinearity acts only on the velocity field h 1 . As can be seen in the proof of Theorem 4, the need to obtain a sign on the derivative of the energy function leads to this constraint, but there is an interesting pointwise argument which also explains the need to constrain the size of Λh 0 . In order to be able to continue our solution for all time, it is necessary to ensure that the initial stability condition h 1 < 0 is propagated, and that we thus maintain h t (·, t) < 0 for all t ≥ 0 in which case we must also have that To give another (and perhaps more transparent) explanation of how this can be achieved, we introduce the time-dependent point x t which satisfies and we define the time-dependent function y(t) by The Hilbert transform, acting on periodic 2π-periodic functions, is defined by We time-differentiate this formula and use integration-by-parts (following the computation in [4]):

Thus, we have that
From (59), it then follows that Now, the condition that h t (·, t) < 0 is equivalent to showing that From the ODE (59), we see that in order for (60) to hold, we must place a size restriction on Λh(x t , t). This size restriction, in turn, requires L ∞ -control of Λh(x t , t); in particular, we must have that 16Λh(x t , t) ≤h 2 1 . It is, therefore, somewhat remarkable that our proof of Theorem 4 requires only control on Λh 0 0.5 rather than Λh 0 L ∞ .
Remark 5. Finally, we note that the asymptotic condition (55) remains true for higherorder Sobolev norms if the condition (56) is replaced with further constraints on the initial data involving higher-order time-derivatives of h evaluated at t = 0. In particular, we do not claim that (58) is sharp.

General interface parameterization with interface turn-over
Our h-model equation (32) for the evolution of the height function h(α, t) uses a special parameterization in which the interface Γ(t) is constrained to be the graph (α, h(α, t)). While this model works well in predicting the mixing layer, in the unstable regime and when the RT instability is initiated, the height function h(α, t) can only grow in amplitude. Our goal is to generalize the h-model equation (32) by using a general parameterization z(α, t) that permits the interface to turn-over. As we will show, the ability for the wave to turn-over, rather than only grow in amplitude, provides an even more accurate prediction of the RT mixing layer, at the expense of a slightly more complicated system of evolution equations.
We now return to the general evolution equations (13) and (19) and set c(α, t) = 0; hence, the interface parameterization z(α, t) = (z 1 (α, t), z 2 (α, t)) evolves according to while the vorticity amplitude satisfies Using the difference quotient approximation for the derivative, we obtain the approximation Substitution of (63) into (61) and (62), and using the Tricomi relation (28), we obtain the general interface RT model evolution equations which allow for wave turn-over: Substituting the time-derivative of (64a) into (64b) yields Notice that (64a) implies that Thus, we obtain the equivalent second-order nonlinear wave system for z(α, t) given by We shall refer to the system (64) or the wave equation (66) as the z-model. The z-model (66) is analogous to the slightly simpler h-model (32).
Remark 6. Note well that due to our choice of setting c(α, t) = 0, the evolving interface solving equation (64) (or equivalently (66)) does not remain a graph, even if the initial data z(α, 0) = (α, h 0 (α)), z t (α, 0) = (0, h 1 (α)) are given as graphs. In particular, it is convenient (especially for the purposes of comparing against the h-model (32)) to prescribe the initial interface position as a graph, and allow the interface to evolve into a non-graph state. As we will show, the z-model (64) captures the turn-over of the RT interface. 9 Numerical study

The algorithm
In order to stabilize numerical oscillations without affecting the amplitude or speed of wave propagation, we shall employ an arbitrary-order artificial viscosity operator for both the h-model (29) (or (32)) and the z-model (64) (or (66)).

Numerical approximation of the h-model
We first consider the h-model, written as a system in (29). We numerically discretize the following approximation: where ǫ > 0 is the artificial viscosity, and s ≥ 2 determines the order of the artificial viscosity employed. Equivalently, we have the approximation for the wave equation given by The term ǫΛ s h ǫ t represents the artificial viscosity operator. The parameter s determines the order of the operator; for example, for s = 2, we recover the classical Laplace operator, while for s > 2, we can study a variety of hyperviscosity operators. We believe that this equation will be an ideal candidate for the C-method artificial viscosity which is localized in both space and time (see [9]), and shall implement this in future work.
We use a Fourier collocation method to solve (67). We note that for the numerical simulations that we consider herein, we choose initial data for which ̟ 0 = 0, in which case the equation that we discretize is equivalent to For N = 0, 1, 2, ..., we define the our Fourier approximation using our projection operator P N (introduced in Section 7), defined as Hence, the mesh size of our algorithm is given by We make use of the following identities in frequency space: It follows that in frequency space, the system (67) can be written as This is a nonlinear system of ordinary differential equations, to which we shall apply the adaptive Runge-Kutta-Fehlberg fourth-order (nominally fifth-order) scheme to advance in time-increments . For the following simulations, we have set the acceleration to a constant value of g, which will either act downward or upward, depending on the type of RT instability that we examine. The initial position of the interface is specified, as well as the initial amplitude of the vorticity, ̟(α, 0).
For our numerical simulations, we shall only consider the case of zero surface tension. We use the same Fourier-collocation method described previously (with N Fourier modes) to approximate the following system of equations: The system (69) employs an artificial viscosity term ǫ∂ 2 α ) to stabilize small-scale noise.
9.2 Simulation 1: h-model, ρ + ρ − = 1 1.5 and Atwood number A < 0 We first study the effect of the density ratio on the interface motion given by the simple RT h-model (67) in the absence of surface tension (σ = 0). For our first simulation, we consider two fluids with density ratio 1/1.5. Specifically, we consider the physical parameters set to g = 9.8 m/s, σ = 0, In order to study convergence of our scheme and the effect of the artificial viscosity operator, we perform a number of different simulations, varying the total number of modes N that are used, as well as the power on the artificial viscosity operator s and the size of the artificial viscosity parameter ǫ. In particular, we consider the following cases: • N = 2 7 , ǫ = 0.01, s = 3 (red solid line in the Figures 2 and 3) • N = 2 7 , ǫ = 0.008, s = 3 (blue solid line with + markers in the Figures 2 and 3) • N = 2 8 , ǫ = 0.008, s = 3 (green dash line in the Figures 2 and 3) • N = 2 7 , ǫ = 0.04, s = 2 (black dash-dot line in the Figures 2 and 3) The results are shown in Figures 2 and 3.  Note that the results are qualitatively similar for various values of parameters, and, in particular, show the convergence of the numerical solutions with mesh size N . A comparison of the blue and the green curves in Figures 2 and 3 demonstrates this nicely. This corresponds to the density of air (ρ + ) and ocean water (ρ − ). We consider the initial data (70) and (71) for (67). The artificial viscosity ǫ = 0.05 and the order of the artificial viscosity operator is s = 3. The number of nodes is N = 2 7 . The results of this simulations are shown in Figure 4, wherein, we once again see the fast growth of the interface h(t) L ∞ over a short time-scale, followed, by decay to equilibrium. The heavier fluid in this simulation, as compared to Simulation 1, reduces the frequency of oscillations, as the interface decays to the rest state. From the intial data, the amplitude grows from 0.7 to 4.3 in a time scale of length 0.2 (between t = 0.2 and t = 0.4). After this remarkable growth, the amplitude of the interface decays and approaches the rest state (see the blue curve in Figure 4). In order to demonstrate the role of the nonlinearity in the growth of the interface, we compare against the linear h-model (25). As expected, the linear model simply decays the interface amplitude; see the red curve in Figure 4.
Finally, in Figure 5 we plot the energy spectrum associated to the energy law (45a) in Proposition 2, as a function of k at t 0 = 0, t 1 = 0.2, t 2 = 0.45 and t 3 = 0.7. We note that outside the Fourier modes k ∈ [−50, 50], the energy spectrum E(k, t) is of order 10 −4 for the time interval considered. Starting from a mode-3 initial interface shape, the energy content is distributed into the smaller scales.    At time t = .45 when the interface is of maximum amplitude, the energy content is welldistributed among all large scales |k| ≤ 20.

Simulation 3: h-model, fingering instability, Atwood number A > 0
We next simulate the highly unstable case of a heavy fluid on top of the lighter fluid, and with the acceleration acting downward. We consider the following physical parameters: g = 9.8 m/s, σ = 0, ρ + = 10 kg/m 3 , ρ − = 1 kg/m 3 , so A = 0.81818 .
We once again use our order one initial data (70) and (71) for (67). The artificial viscosity operator is order s = 3, and the artificial viscosity parameter is set to ǫ = 0.05. The number of nodes is N = 2 7 .
This simulation is intended to demonstrate the ability of the h-model to show fingering phenomenon; indeed, as can be seen in Figure 6, the heavy fluid penetrates the lighter fluid and a strong RT instability is initiated. After a large enough time interval, the absolute value of the interface grows exponentially fast.  To capture the behavior described in Theorem 4, we simulate (68') (instead of (67)). We consider the physical parameters g = 9.8 m/s, σ = 0, ρ + = 0 kg/m 3 , ρ − = 1 kg/m 3 , so A = −1 , and the initial data For this initial data, the homogeneous solution h ∞ (t) is given by Theorem 4 states that h(·, t) − h ∞ (t) converges to zero as t → ∞.
As the initial data satisfies the stability condition, no artificial viscosity is required to stabilize the numerical solution, and we set ǫ = 0. Finally, we fix the number of Fourier modes to N = 2 7 . The results are plotted in Figure 7, where the simulation demonstrates the asymptotic behavior of the theorem. We consider now the situation where two (nearly) incompressible fluids, of densities ρ ± , are subjected to an approximately constant acceleration g normal to the interface separating them, directed from the lighter fluid to the denser fluid. We assume that the initial interface is given by a small and random perturbation of the flat state.
Our goal is to compare the growth rate of the mixing layer using our model equation (32) with that predicted by experiments and numerical simulations of Read [8] and Youngs [14]. Experiments show that if the instability arises in the previous setting, the width of the mixed region grows like t 2 . Actually, as shown by Read [8] and Youngs [14], the mixing region grows as Direct Numerical Simulation in two space dimensions by Youngs [14] has indicated that the parameter δ should range from 0.04 to 0.05, whereas experiments of Read [8] suggest that δ should range from 0.06 to 0.07. In particular, when we have a NaI solution (ρ − = 1.89g/cm 3 ) and Hexane (ρ + = 0.66g/cm 3 ) in a tank, the empirical value is δ = 0.063 (see [8]). Let us emphasize that both the numerical simulations and the physical experiments were run for approximately 70ms.

The h-model
We consider the initial data for the h-model, given by where a j , b j are random numbers following a standard Gaussian distribution, and S denotes a normalization constant such that We consider n = 50, N = 2 7 , and σ = 0. The acceleration g = − 9.8·2·π 0.3 L/s 2 acts upwards, where L is chosen so that 2π/0.3 L = 1m. This gravity force corresponds to the usual gravity force on the surface of the Earth of 9.8m/s 2 . Finally, we use the artificial viscosity parameter ǫ = 0.05 together with a second-order artificial viscosity operator s = 2.
The mixing layer is shown in Figure 8, where the interface position h(x, t j ) is displayed for times t 1 = 0, t 2 = 0.02, t 3 = 0.04, t 4 = 0.06 and t 5 = 0.08.  In Figure 9, we see that up to time t = 150ms, the numerical solution provides a mixinglayer growth rate which agrees well with the predicted growth rate given by (72); the largest difference between the h-model and (72) is given by For large times, due to the strong RT instability present in the equations, the results of numerical simulations are very sensitive to the artificial viscosity parameter (see [2] for the artificial viscosity effects on RT mixing rates); however, for short time (meaning around 70ms, neither viscosity nor nonlinearity plays a critical role in the evolution. To demonstrate this, we perform a numerical simulation of the h-model with zero artificial viscosity ǫ = 0, and with surface tension σ = 0.005. We use the initial data satisfying (73)-(74) with n = 30, and S chosen such that h(0) L 2 = π 1000 .
As can be seen in Figure 10, the growth rate predicted by our model agrees well with the Youngs' growth rate with δ = 0.06 up to around t = 60ms.

The z-model
Now, we repeat the "rocket rig" experiment of Read and Youngs using the z-model, with initial data given by where a j , b j are random numbers following a standard Gaussian distribution, and S denotes a normalization constant such that We consider n = 30 and N = 2 9 . The artificial viscosity for (69) has coefficient ǫ = 0.01. Figure 11 shows the interface evolution for the z-model. The ability of the z-model  parameterization to "fatten" and "finger" produces an even more accurate representation of the mixing region.
We are able to quantitatively validate this statement by comparing the growth of the mixing region of the z-model with the h-model and the quadratic prediction (72) of Read and Youngs. For the comparison, we simulate the h-model with the same physical parameters, and with second-order artificial viscosity with ǫ = 0.05, and with initial data given by We note that the artificial viscosity ǫ is five times larger for the h-model (67) than for the z-model (69). As shown in Figure 12, the z-model mixing region growth rate matches very well with the quadratic prediction (72) of Read and Youngs, and for a relatively long time interval, up to t = 150ms. The width of the mixing region, approximated by  Moreover, the z-model has a longer lifespan than the h-model. Since the h-model is constrained to remain a graph, the amplitude of the interface can only grow rapidly once the RT instability is strongly initiated; on the other hand, as the z-model can turn-over, the instability creates a turn-over and horizontal-fattening of the mixing region. The comparison is shown Figure 13 9.7 Simulation 6: The "tilted rig" experiment Finally, we consider the "tilted" rig experiment of Youngs [14], wherein the tank of the rocket rig experiment is titled by a small angle θ from the vertical.
Again, we have two (nearly) incompressible fluids, of densities ρ ± , with a vertical acceleration g directed from the lighter fluid to the denser fluid. As Youngs notes, the inclination of the initial interface results in a gross overturning motion in addition to the fine scale mixing.
We define θ = 5.7 • and consider the unperturbed (flat) tilted interface given bỹ if π/2 < π ≤ π. We assume that the actual initial interface is given by a small and random perturbation so that z 2 0 (x) =z 2 0 (x) + S n j=1 a j cos(jx) + b j sin(jx), where a j , b j are random numbers following a standard Gaussian distribution, and S denotes a normalization constant such that .
Again, we employ a second-order artificial viscosity operator with ǫ = 0.25, a rather large artificial viscosity, but necessary to stabilize the interface in the highly unstable RT regime. The results are plotted in figure 14.

The z-model
For the z-model, we use the initial conditions δz 1 (α, 0) = 0 , z 2 (α, 0) = z 2 0 (α) , ̟(α, 0) = 0 , with artificial viscosity ǫ = 0.05. As can be seen in Figure 14, as the RT instability is strongly initiated, the amplitude of the h-model starts to grow unboundedly, while the z-model can turn-over and continue to run for a much longer time interval. We show the time evolution of the interface for the h-model and z-model in Figure 15(a). Because the h-model can only grow the amplitude of h, the behavior of the interface is more singular for the h-model than for the z-model. In particular, h(t) L ∞ grows faster than z 2 (t) L ∞ . Moreover, the z-model is more stable, and permits the use of smaller artificial viscosity than the h-model.
Finally, as we have run the titled rig experiment with random initial data, it is interesting to plot the ensemble of runs, and get a clear picture of the mixing region. In particular, as it may be difficult to infer the qualitative description of the mixing region from only one simulation, in Figure 16, we plot the z-model interface location for a number of different random initial conditions at time t = 0.22. As can be seen, the ensemble gives an approximation of the mixing region provided by experiment [14] and DNS [10]. The mixing region grows slightly faster on the left side than it drops on the right side as expected by Youngs [14], and has the qualitative features of the experiment. The z-model can simulate the Kelvin-Helmholtz instability arising in the case of equal densities ρ + = ρ − , equivalently A = 0. We note that when A = 0, due to (64), we have ̟ t = 0, and the problem reduces to describing the interface evolution for z(α, t) via (69a,b). We consider the initial data given by  and use N = 2 8 Fourier modes. We fix the artificial viscosity as ǫ = 0.01. The results are given in Figure 17. The primary effect of the nonlinearity is to increase the length of the interface by starting to roll-over, while keeping the curve smooth. This can be seen by looking at the spectral content of the the L 2 -energy function E(k, t) = |δẑ 1 (k, t)| 2 + |ẑ 2 (k, t)| 2 at different instances of time. In Figure 18 we plot the spectrum of E(k, t) at times t 0 = 0, t 1 = 0.2, t 2 = 0.4 and t 3 = 0.6. For large Fourier modes |k| > 10, the L 2 -energy spectrum E(k, t) is of order 10 −5 for the time interval considered, so we plot E(k, t) versus −10 ≤ k ≤ 10. The energy spectrum remains fairly localized about the k = 1 initial data, with some growth in wavenumbers |k| = 2 and 3, which are responsible for the dilation and rotation of the wave.