## Abstract

We study the Rayleigh–Taylor (RT) mixing layer, presenting simulations in agreement with experimental data. This problem is an idealized subproblem of important scientific and engineering problems, such as gravitationally induced mixing in oceanography and performance assessment for inertial confinement fusion. Engineering codes commonly achieve correct simulations through the calibration of adjustable parameters. In this sense, they are interpolative and not predictive. As computational science moves from the interpolative to the predictive and reduces the reliance on experiment, the quality of decision making improves. The diagnosis of errors in a multi-parameter, multi-physics setting is daunting, so we address this issue in the proposed idealized setting. The validation tests presented are thus a test for engineering codes, when used for complex problems containing RT features. The RT growth rate, characterized by a dimensionless but non-universal parameter *α*, describes the outer edge of the mixing zone. Increasingly accurate front tracking/large eddy simulations reveal the non-universality of the growth rate and agreement with experimental data. Increased mesh resolution allows reduction in the role of key subgrid models. We study the effect of long-wavelength perturbations on the mixing growth rate. A self-similar power law for the initial perturbation amplitudes is here inferred from experimental data. We show a maximum ±5% effect on the growth rate. Large (factors of 2) effects, as predicted in some models and many simulations, are inconsistent with the experimental data of Youngs and co-authors. The inconsistency of the model lies in the treatment of the dynamics of bubbles, which are the shortest-wavelength modes for this problem. An alternative theory for this shortest wavelength, based on the bubble merger model, was previously shown to be consistent with experimental data.

## 1. Introduction

Turbulent mixing is a long-standing problem. We consider the Rayleigh–Taylor (RT) mixing problem [1,2], in which a light fluid accelerates a heavy fluid. Numerical efforts to solve this problem date back more than 60 years, and only recently have systematic numerical solutions based on compressible fluid dynamics been found that address an experimental range of fluid parameters for both immiscible and miscible fluids, the latter over a range of Schmidt numbers [3–6]. The problem is to predict the growth rate *α*, defined in terms of the penetration distance *h* of the light fluid into the heavy fluid, via the equation
1.1where *A*=(*ρ*_{2}−*ρ*_{1})/(*ρ*_{2}+*ρ*_{1}) is the Atwood ratio, a dimensionless measure of the density contrast, and *g* is the acceleration, applied normal to an interface separating fluids of density *ρ*_{1} and *ρ*_{2}.

We identify six dimensionless parameters on which *α* displays sensitivity. These are analysed in [3,4] and summarized in §3. It is argued [7] that the ratio of momentum lost to that gained is statistically less noisy than *α* as an integral characterization of the overall mixing rate.

The scope of this paper is the status of numerical simulations and theoretical models to determine *α* in the regime of laboratory experiments and the comparison of simulation and theory to experiment in this regime.

As documented below, these simulations show agreement with experiment. An additional simulation [8], based on an incompressible code, shows agreement with experiment for moderate Schmidt numbers. We also mention our earlier simulations, for example [9], stretching over a number of years, with increasing agreement to experiment, and we mention the particle-based simulation, with mixing rates within the experimental range [10]. These simulations are parameter-free in their agreement with experiment.

To discuss validation, we first introduce a definition. Validation of a simulation refers to the comparison of a simulation to an experiment, with a small error. Of course, the issues of the size of the error, the range of distinct experiments for which agreement is obtained and the type of measurements used for which approximate agreement is reported are then open issues. For such reasons, validation is seldom complete. We do exclude simulations with adjustable parameters, for which the parameter fitting experiments and the validation checking experiments are not independent.

Other simulations [11,12] show agreement with experiment, but make use of an adjustable parameter (the magnitude of an assumed long-wavelength perturbation in the initial data). Agreement is achieved by the selection of a suitable value for this parameter. As the class of experiments (with a relatively common *α*) used to adjust this parameter coincides with the data being used for comparison to experiment, such simulations do not meet the requirements of validation. Many simulations (see for example [13]) disagree with experiment in their values of *α* by factors of 2 or more, a fact commonly attributed [11,12,14] to unmeasured long-wavelength perturbations in the initial conditions.

Questions have been raised regarding our validation studies, a particularly interesting issue in view of the above noted simulation–experiment discrepancy for many codes.

The following three issues have been raised:

(1) Long-wavelength perturbations could be present in the initial data and could contribute to the self-similar growth constant

*α*by a factor of 2 or more.(2) Simulations that depend on subgrid models leave open the validity of the subgrid models.

(3) Some workers have conjectured the existence of a ‘true’ self-similar regime, universal relative to fluid transport parameters, which may be distinct from regimes of laboratory experiments.

Item 3 will be addressed in future publications.

Item 2 appears to be at most marginally credible, in view of:

— the 20-year history of validated turbulence computations using dynamic subgrid-scale (SGS) models of the type considered here by a large community of computational scientists and engineers and for a wide range of problems;

— our own multiple validated RT simulations with a maximum 6% and usually 0% discrepancy (relative distance outside of reported error bounds) with experiment, including examples with measured initial conditions used in the simulation, in contrast to many non-validated RT simulations using compressible hydro codes with a typical factor of 2 discrepancy to experimental data; and

— the highly resolved simulation shown later in figure 1, with full resolution of the critical Weber length scale for droplet break-up, shows agreement with experimental data and convergence under mesh refinement; additional mesh convergence studies are presented in [3,4].

Item 1 is addressed here, as the main purpose of this paper. Central to the item 1 claims addressed here is the assertion that long-wavelength initial perturbations account for a factor of 2 or more in *α* in the regime of laboratory experiments and that, in the absence of such perturbations, simulation results with *α*≈0.02 or 0.03 are correct.

Our main conclusions regarding item 1 are as follows:

(1a) Multiple simulations achieve

*α*values in highly accurate agreement with experiment. The*α*values are not universal, but depend on specific experimental parameters. The simulations have sufficient accuracy to resolve these differences. A number of parameters showing*α*sensitivity are identified.(1b) Initial data can be reconstructed for the Smeeton–Youngs rocket rig data.

(1c) Even with generous error margins for the reconstruction, the resulting uncertainty in

*α*is ±5%.(1d) The Smeeton–Youngs experiments have significantly weaker initial long-wave perturbation strength than do the Andrews water channel experiments.

(1e) The long-wave perturbations as measured for the water channel experiments are sufficient to explain perhaps a 30% change in the value of

*α*but do not explain a change by a factor of 2.(1f) The theoretical model [14] for bubble dynamics appears to be less suited to describe experimental data than is the bubble merger model [16,17]. We also mention the bubble merger model [18–20], which correctly predicts the mixing growth rate

*α*, but in contrast with [17], incorrectly predicts the experimental values of the bubble height-to-width ratio.

Items (1b–d) and (1f) are newly reported here.

New results presented here include new simulation studies [15] for experiment no. 105, new simulations for the water channel experiments [21], the analysis and inference of initial conditions for experiment no. 105, the error bounds and uncertainty analysis for these initial conditions and the prediction of *α* for this experiment.

## 2. Turbulent mixing

### (a) Front tracking and large eddy simulations

The Navier–Stokes equations solved are the filtered continuity, momentum, energy and concentration equations for two miscible fluid species in an inertial frame, as written in equations (2.1)–(2.4). The filtered quantities are considered to be mesh block averages, and are denoted with an overbar, whereas mass-averaged quantities are denoted with a tilde. Repeated indices are summed:
2.1
2.2
2.3
and
2.4where the SGS variables, *τ*_{ij}, , , and , are expressed as
2.5
2.6
2.7
2.8
and
2.9The viscous stress tensor, *d*_{ij}, in the momentum and energy equations is expressed as
2.10where is the filtered dynamic viscosity. In equation (2.3), and are the partial specific enthalpy of each species defined by
2.11and
2.12where and are the specific internal energy of each species.

The essential features of the algorithms used here are front tracking, to achieve resolution of steep and sharp density gradients, and large eddy simulations (LES) with SGS terms, to model the diffusive transport corrections to the mesh (Reynolds) averaged Navier–Stokes equations. These features are included in the multipurpose simulation code FronTier, as we now describe.

Front tracking assigns computational degrees to a surface, or interface, moving freely through an Eulerian grid. The discretization of the Eulerian grid variables is conventional, using the MUSCL algorithm [22], with the exception that, for stencils cut by the front, extrapolated ghost cell values are used to complete the stencil so that it contains states from a single side of the interface. The front variables acquire a velocity field via a two-sided interpolation of states from the interior, projected onto a single normal velocity field through solution of a Riemann problem [23]. The front is the location of a sharp discontinuity in fluid properties, in the case of immiscible fluids, and it is an isoconcentration contour in the case of miscible fluids with a steep concentration gradient.

It is well known that Eulerian advection schemes suffer from significant levels of numerical mass diffusion, in the case of a steep concentration gradient advected by the fluid relative to a grid. Lagrangian methods can be used to mitigate this problem but they suffer from mesh distortion and remeshing interpolation errors when taken to complex interfaces at late time. Volume of fluids methods [24,25] provide a reconstruction of the interface, in order to overcome these problems. Level sets, in which the level set distance to the interface is advected, have also been used to mitigate Eulerian mass diffusion problems. The ghost cell extrapolation idea, at the heart of the level set approach, is based on the same idea, part of the original front tracking algorithm [26]. Level sets come in several different versions, the simpler of which are of too low quality to be very useful and the more complicated of which offer little advantage in regard to simplicity over front tracking. A comparison study [27] found front tracking to be of higher resolution than either volume of fluids or level sets.

The volume of fluids algorithm was used in one of the *α* group comparison methods [13] and showed no advantage over the untracked codes, with the same factor of 2 disagreement with experiment. Level sets have not been used for turbulent mixing problems. In view of the importance of high-quality resolution for steep concentration gradients, further comparison of these methods is called for, especially for the type of problems considered here, with experimental comparison possible and using the most advanced version of each of these basic methods; for the level sets, the method of [28] should be part of the comparison.

The LES/SGS algorithm is classical [29–31]. The salient point here is that the LES/SGS algorithm assumes a functional form for the turbulent transport terms, basically of a gradient diffusion nature, but the coefficient of this term, a potentially troublesome free parameter, is actually constrained and determined from the simulation itself. The derivation of the Reynolds averaged equations and the SGS terms with the dynamically determined (parameter-free) transport terms is too detailed, with formulae too lengthy to be summarized here, but reference [31] is an excellent self-contained source for this material. The extra equation needed and used to determine the turbulent transport coefficients dynamically comes from looking at two adjacent grid levels, the current grid and a once coarsened grid. On the once coarsened grid, the missing terms and coefficients can be determined in two ways, leading to a dynamic equation for their determination. Assuming an asymptotic relation, the coefficient on the once coarser mesh determines the coefficient on the current mesh. The combination of these two algorithms appears to be very powerful in its ability to overcome difficulties that have impeded progress in a variety of turbulent mixing problems and for a variety of researchers.

In addition to the algorithmic issues discussed above, we have emphasized the careful attention to all details of the experiments being modelled. For example, in comparing with nominally similar fresh–salt water mixing in distinct experiments by different groups, we found that six separate variables changed between the two experiments, and we found [4] sensitivity of *α* to each of these six differences. Our simulations were of sufficient accuracy to reflect the differences between the experiments and to assess each of the six factors that are different between them in terms of their influence on *α*. With the exception of the work of Andrews and associates, experimental modelling issues are often neglected in simulations. In fact, the belief that *α* is universal and insensitive to experimental details (refuted by our simulations) encourages this disregard (see §3).

All simulations reviewed here are in three dimensions. The sides of the computational domain (parallel to the direction of the gravitational field *g*) are periodic; the top and bottom are no-flow boundaries.

### (b) Front tracking Rayleigh–Taylor simulations

Simulations in agreement with 14 experiments were previously reported [3–6] (table 1). We return to three of these experiments with new simulation and data analysis (uncertainty quantification) results.

We report here improved agreement with two experiments [21] on the basis of improved modelling of the initial mass diffusion layer. The initial mass diffusion layer thickness was not reported in [21], while it was reported for the hot–cold water case [8]. Using simple *x*^{2}/*t* scaling based on the diffusion equation, we relate the fresh–salt water initial mass diffusion layer via the square root of the ratio of diffusion constants to the reported initial diffusion layer for the hot–cold water experiment. The initial mass diffusion layer is changed from 0.3 to 0.015 cm, in order to fit the first experimental data point, resulting in a change in *α*_{sym} from 0.075 to 0.084, reducing the discrepancy to experiment from 6% to 0%. This value for the initial mass diffusion layer largely governs the earliest-time experimental data point for the early-time molecular mixing.

We return to experiment no. 105 of [15] and achieve agreement with experiment, and with simulations going through all data points (table 1). Simulations I, II and III presented in figure 1 and related to uncertainty quantification for the initial conditions are new. Simulation III has a much finer grid, of 111 μm, sufficient to fully resolve the Weber length scale associated with droplet break-up, the dominant physics for the immiscible mixing of this experiment. The initial conditions of simulation III are the best estimate reconstruction of the initial conditions from the measured data. The agreement among the three simulations and these with experiment is sufficient to assert the mesh convergence and approximate mesh independence of the results. Further mesh convergence studies were reported in [3,4].

### (c) Numerical mass diffusion and numerical simulations of *α*

It is common practice [12,13] to model RT instability with four or five mesh blocks per elementary initial perturbation mode (bubble). Such resolution, in the absence of front tracking, leads to excessive levels of numerical mass diffusion. We have analysed [34] this effect, comparing a tracked to an untracked simulation [35]. The compared untracked simulation used eight mesh blocks per initial mode wavelength and an artificial compression scheme to reduce numerical mass diffusion. It yielded *α*=0.04, one of the highest values for *α* reported for untracked simulations, probably indicating a lower level of numerical mass diffusion than that which occurs in many untracked simulations. We introduced a time-dependent Atwood number *A*(*t*) defined in terms of the simulation itself, as a measure of the time-dependent density contrast within the simulation bubble region [6,34,36]. We found that the untracked simulation reduced the density contrast, and *A*(*t*), by a factor of ∼2 relative to the nominal *A*=*A*(*t*=0). Defining , we found improved self-similarity and an *α*_{ren} in agreement with typical experimental values and with the tracked simulation. In other words, this low value of *α* and apparently a number of other reported even lower values for *α*≈0.02 or 0.03 are unduly influenced by numerical artefacts. In physical modelling terms, such untracked simulations correspond to low but unquantified numerical Schmidt numbers, perhaps in the range 0.3–0.1. According to the implicit LES (ILES) solution methodology still in common use, physical transport coefficients are set to zero, turbulent transport is ignored, and thus all fluid transport is regulated by the numerical algorithms, with undocumented Schmidt, Prandtl and Reynolds numbers. Conclusions [12] regarding the sensitivity of such simulations to long-wavelength noise may thus pertain to a regime of very low Schmidt number, on the basis of the analogy between numerical and physical mass diffusion.

### (d) The experiments

We consider five immiscible experiments, nos. 56, 63 [37] and 104, 105 and 114 [15]. Immiscible experiments were selected, as we analyse amplitudes of individual bubbles and these are more sharply defined in the experimental record than are the bubbles in the miscible experiments. These five experiments were selected on the basis of sufficient supporting data in the experimental reports. We chose no. 105 for numerical simulations as it had the lowest Atwood number, making it less costly to simulate with a compressible explicit algorithm. Of the five experiments, one (no. 63) was deliberately vibrated to create a low-amplitude standing wave as a part of the initial conditions. This experiment has a clearly defined large central bubble, distinct from the statistical ensemble of bubbles around it. One (no. 104), according to the comments in the experimental report, has mode 2 perturbations initially. One (no. 56) has, according to our own observations, a possible long-wavelength perturbation in its early time steps. Two (nos. 105 and 114) show no visual signs or experimentally reported indications of long-wavelength perturbations. There seems to be no correlation between the observations of long-wave perturbations or their absence and the values observed for *α*. The two experiments nominally free of such perturbations have the largest value of *α* (no. 105: *α*=0.072) and the next to smallest value of *α* (no. 114: *α*=0.060) of the five.

In view of this experimental record, there appears to be no basis for concluding that run-to-run variations in the long-wavelength perturbations explain experimental variability of ±10% in *α*. An assertion that a factor of 2–3 modification in *α* can result from long-wavelength perturbations in the experimental data seems still less plausible. With the two significant digit accuracy (or, more realistically, the agreement within 10%, in the absence of experimental error bars) of our simulations in comparison to experiment, it would seem that the experimental variability in *α* will be explained in some other aspect of the experiments.

### (e) Initial perturbations, *α* and numerical simulations

We have three numerical simulation tests that assess the effects of long-wavelength initial perturbations on *α*. The first two involve the Mueschke–Andrews experiments [21], for which early-time perturbations were measured. The Mueschke–Andrews initial measurement time determination is not an initial time but a small distance downstream, at which point the perturbation amplitudes are visible and are measured. This observation distance (in effect, an observation time) appears not to be very different from the time of the early plates of the rocket rig experiments, in terms of dimensionless perturbation amplitudes. Thus, we believe that the Mueschke–Andrews early-time (‘initial’) data are, in this sense, comparable with the early-time Youngs *et al.* data we consider here.

Mueschke and Andrews have significantly stronger long-wavelength perturbations initially, compared with Youngs *et al.* data at the third experimental plate (the first with data that can be used reliably for bubble amplitudes).

We simulated the Mueschke–Andrews problem [3] with and without the long-wavelength component of the initial data, and found an increase of *α* by 0.023, or 27% of the experimental value 0.085±0.005 due to the long-wavelength perturbations as measured experimentally. This increase does not account for the simulation values of *α*≈0.02 or 0.03 obtained in simulation by others, nor does it allow for a factor of 2–3 variation in *α* due to long-wavelength perturbations.

We report here simulations of experiment no. 105 [15] (figure 2). According to methods explained in §2*g*, we reconstruct the *t*=0 data from early-time experimental plates, and through estimates of the errors in this reconstruction explained in §2*i*, we bound *α* in a range 0.076±0.004. The experimental *α*=0.072 falls within this uncertainty quantification simulation interval. The influence of the long-wavelength initial perturbations on *α* is less than or equal to 5%.

The purpose of the reconstruction and the simulation is validation and uncertainty quantification. We use knowledge of the early-time and inferred initial data to predict the late-time behaviour in agreement with experiment. In a validation study, use of experimental data (which is independent of the observed data) is appropriate. The early-time spectral amplitudes of the perturbations are independent of the long-time growth rate *α*, other than the fact that the two are linked by solutions of the Navier–Stokes equations. It is precisely this linkage (i.e. the Navier–Stokes equation, solved by a front tracking/LES algorithm) that we investigate and validate. For this purpose, the initial or early-time data are an input, not an output, and are logically independent of the dynamical and numerical evolution output being validated. This part of our work is a validation result. We assess uncertainties in the reconstruction and carry forward these errors to late-time results to assess uncertainty in the late time associated with partial knowledge or partial uncertainty of the long-wavelength component of the initial conditions. This is uncertainty quantification. These issues are illustrated in equation (2.13), with the left inference diagram representing a classical validation inference (input plate 1 to output plate 6) and the right, the validation inference diagram we employ. All arrows represent solution algorithms. The first arrow on the right is backward in time and uses linearized equations and analytical formulae from dispersion relations (plate 3 to input plate 1 to final plate 6):
2.13

### (f) Modal analysis of experimental data

In a companion paper to appear in *Physica Scripta*, we show the experimental plates to be analysed below. In the plates, the acceleration is downward, and the inertial force (the effective ‘gravity’) is upward. Thus the light fluid, initially above the heavy, is forced down. The bubbles are elements of light fluid penetrating downward into the lower, heavy fluid and are clearly evident in these experimental plates.

After a brief initial time period, the interface separating the two phases cannot be represented as a (single-valued) function *z*=*z*(*x*,*y*). In the experimental plates, only the upper and lower envelopes of the interface are visible, as viewed along the *x*-axis direction, so that the interface envelope is displayed as a function of *y* alone. Thus, an analysis of the interface based on experimental data must start with some simplifications. Following [38,39], we discard all information other than the sequence of bubble minima at each time. In effect, we are assuming that the bubbles have an equal width and that the height of each bubble minimum has sufficient information for the subsequent analysis.

We start by recording the amplitudes of all bubble minima at the times *t*=*t*_{j} of each of the *j*=0,1,… experimental plates. The decision about which local minima to disregard (small minima in the middle of a spike, representing a bubble already lost from the ensemble or closely adjacent bubbles, partially merged, to be counted as one bubble) is partially subjective. As a check, the data acquisition was repeated independently by two authors, with only minor differences observed. For example in table 3, the independently read data led to average velocity differences of about 10%, with no significant change in the conclusions.

The minima are taken relative to the experimentally recorded extreme minima, so that parallax issues, addressed in the experimental data analysis, are also addressed in our data analysis. The number of such minima is time-dependent, decreasing with time. We regard the amplitudes at the minima as functions of discrete equally spaced locations , or conceptually as a function defined on , the group of integers reduced modulo . A discrete Fourier transform of a function *f*_{j} on is defined by the formula
2.14where [*x*] is the integer part of the real number *x*. Then (2.14) defines a function on the dual group, which is also . Because the *f*_{j} are real, the satisfy an identity , and in the arithmetic of , . Accordingly, such modes, above the Nyquist limit, are omitted, as is indicated in (2.14). The complex Fourier series (2.14) gives rise to independent sine and cosine modes for the same wavenumber, so there are real Fourier modes in the series.

Because we are transforming the bubble minima and not the full bubble envelope, the *n*=0 mode has a special interpretation. As can be seen from equation (2.14), the *n*=0 mode determines the mean bubble penetration height (mean bubble amplitude) at *t*=*t*_{j}. The bubbles themselves are the shortest-wavelength modes present at *t*=*t*_{j}, as is evident from inspection of the experimental plates. Thus, the *n*=0 mode determines the dominant (average) contribution to the amplitude of these short-wavelength modes; in this sense, we can say that the *n*=0 mode describes the dominant short-wavelength behaviour in the bubble ensemble.

According to this analysis, we regard Fourier modes in the range as long-wavelength perturbations.

Bubble height fluctuations, which in the experiments being considered here are reduced by an order of magnitude relative to the mean height, are characterized by the higher mode numbers, . In view of this interpretation, we place this *n*=0 mode at the position , after all the other modes. In this manner, we construct modal amplitudes *A*(*n*,*t*) for for *t*=*t*_{j}, *j*=3,4,5,6, and for five experiments (nos. 56 and 63 [37] and nos. 104, 105, 114 [15]).

In , , so the convention regarding the *n*=0 mode is a mathematical identity. But this distinction between *n*=0 and will be important later, when we interpret the modes for the purpose of numerical simulation initial conditions, as defining sine and cosine modes on the periodic interval [0,1] (§2*i*). In this case, the *n*=0 mode and the mode are distinct.

### (g) A theoretical model for modal dynamics

We introduce a theoretical model for the dynamics of the modes *A*(*n*,*t*), based on three ideas:

(1) Superposition, so that

*A*(*n*,*t*), , moves independently of the other modes.(2) Single mode dynamics, so that

*A*(*n*,*t*), , propagates with the exponential growth law of the linear dynamics of dispersion relations. In principle, nonlinear dynamics for larger amplitudes could be used, but for the experimental data analysed,*A*(*n*,*t*) is small relative to the modal wavelength.(3) Special dynamics for the dominant short-wave modes based on the bubble merger model [16,17,40,41]. The bubble merger model propagates short-wavelength modes from a previous step,

*t*=*t*_{j−1}, into short-wave modes for*t*=*t*_{j}. These*t*_{j−1}modes have a higher mode number, so this propagation is not by superposition within a single mode number. It involves a coupling in the dynamics between distinct short-wavelength modes at different times, thus coupling across distinct mode numbers. According to the bubble merger concept, two or more adjacent bubbles at time*t*_{j−1}combine or are ‘merged’ to produce a larger bubble, with a smaller mode number, but still the short-wave mode for the new time*t*_{j}.

The first, and primary, purpose of the modal dynamics model is to run backward in time and reconstruct the unobserved long-wavelength modal amplitudes *A*(*n*,*t*_{0}) of the initial data. This purpose does not make use of the short-wavelength modes, and so it depends only on ideas 1 and 2 above.

The second purpose of the model is to predict observed data *A*(*n*,*t*_{j}) from adjacent times *t*_{j±1}. This exercise generates errors, which we analyse statistically to ascertain the quality or reliability of the reconstructed *t*=*t*_{0} data. In this way, uncertainty quantification simulation error bars related to initial condition uncertainty can be assessed. As we are interested in errors in the predictions of the initial long-wavelength modes, this purpose is also served by ideas 1 and 2 above only. The third purpose of the theoretical model is to contrast it to a model based on the principle of superposition, in which the motion of each Fourier mode of the bubble interface moves independently of the other modes. Such a model is loosely related to that in [14]. For the purposes of comparison, we also take some liberties with the bubble merger model [17]. The two models differ only in their assumptions related to idea 3, which is thus crucial for this third purpose.

### (h) Long-wavelength modes in initial data

A self-similar initial spectral power law *A*(*n*)^{2}∼*n*^{−4} for the *t*=0 amplitudes has been proposed [11,14,42]. A similar power law was reported in the initial conditions for the surface of a glass inertial confinement fusion (ICF) pellet, with higher amplitudes for the smaller wavelengths [14]. We find a power law *A*(*n*)^{2}∼*n*^{−3.3} (mean exponent averaged over five experiments) in the reconstructed *t*=0 spectra (figure 3*a*). The self-similar power law satisfied by the initial data [15,37] is not sustained in the third and all later experimental plates. The long-wave contribution to the total spectral energy drops from 80% or more initially (when considered over five experiments) to less than 0.5% for *t*=*t*_{3} and all later observational times [38]. See figure 3*b* and table 2 for the data for experiment no. 105.

The *Z*_{n} amplitudes are used to construct Fourier modes; the Fourier modes determine the initial data. The data as analysed have amplitude information only and no phase information. Accordingly, we assign random phases for each amplitude. There are no long-wavelength data in the *y* (depth or into plane) direction in the figures, and we set all long-wavelength amplitudes in this direction to zero. Fourier analysis over [0,1] differs from that over *Z*_{n}, in that there are independent Fourier modes for wavenumbers for any value of *n*, and specifically in the short-wave range . In this range, we assign random amplitudes with a total spectral weight *A*(0,0), for both *x* and *y* directions.

### (i) Uncertainty quantification for Rayleigh–Taylor mixing

Here, we quantify the errors associated with the modal dynamics, and thus with the transfer of data between nearby time steps, as based on the superposition assumption for . The error estimates in the *t*=*t*_{0} reconstructed initial data are used to establish an uncertainty quantification interval for numerical simulations of the RT growth rate *α*.

We compare directly observed data to predicted (transferred by modal dynamics) data for five experiments, for each of approximately five modes, /2. For the comparison of time adjacent experimental plates, we propose a formula
2.15to express the discrepancy between observation and prediction. We attribute *b* to the uncertainty or error in reading bubble height data from scanned images of published experimental plates, with an estimated value of 50 μm, and we estimate *a*=1.0. We attribute *a* to errors in the theoretical model for modal dynamics. The error analysis is illustrated in figure 4*a*. In figure 4*b*, we show a similar scatter plot for the *t*=*t*_{0} data reconstruction. We reconstruct twice, from *t*_{3} and from *t*_{4}, and plot the two predictions relative to their mean. In both cases, the error bounds defined by (2.15) are drawn in the figure. We see that the error bounds encompass over 95% of the data, and thus define 95% confidence intervals for the reconstruction.

We comment that the error bound *a*=1 is large and leads to long-wavelength perturbations nominal ±100% bracketed by (a) 0 and (b) 2× nominal. Thus, choice 0 amounts to pure short-wavelength initial perturbations, and choice 2 allows double the inferred long-wavelength initial perturbation plus short-wavelength perturbations. The short-wavelength spectral amplitude is set to agree with the reconstructed *A*(0,0) value. Even with this large bound, the effect on *α* is ±5%. The proper units for *b* are the same as the errors themselves, and these are in micrometres, and reflect the width (fuzzy edge) of the interface in the experimental plates.

We initialize simulation data with long-wavelength perturbations equal to (a) 0 and (b) 2× the inferred long-wavelength initial data. These simulations will give a prediction of *α* with uncertainty quantification bounds related to the initial long-wavelength initial conditions (table 1 and figure 5*a*). In the initialization, we use the Fourier amplitudes of the modes to define sine and cosine modes on a periodic interval [0,*L*]. The mode and adjacent high-wavenumber modes are replaced by a wave packet (uniform in *k* space), spread over the short-wavelength interval, having the same total spectral energy.

### (j) Comparison of bubble merger and superposition models for the dynamics

For , the superposition model replaces idea 3 (bubble merger) with superposition, with an enhanced nonlinear single mode velocity about double the single mode velocity itself. Such a doubling of the nonlinear single mode velocity was observed in an analysis of experiments [43] and incorporated into the derivation of the bubble merger models. Thus, superposition employs one but not all aspects of the bubble merger dynamics.

The critical difference between the bubble merger and the superposition dynamics defined here lies in the propagation of the shortest-wavelength mode, . Within the dynamics, the key difference is the starting point at each time step for the dynamics.

In the bubble merger model, propagation comes from the shortest wavelength at an earlier time, which has a smaller wavelength and higher mode number. In other words, smaller bubbles merge to form larger ones at the next time step. Thus (for example), modes at length *λ*/2 propagate into modes at length *λ*.

In the superposition model, the modes with wavenumber *n* from an earlier time (when they are not the shortest wavelength) propagate to identical wavenumber mode *n* at a later time when they become the shortest wavelength.

To follow the short-wavelength dynamics, we refer to figure 5*b*, in which data points are extracted from figure 3*b*. We study the dynamics leading to the short-wavelength amplitude for time *t*=*t*_{5}, which in the bubble merger model comes from and in the superposition model comes from . These two dynamics require different changes *ΔA* in the mode amplitudes *A*, namely *ΔA*_{1} and *ΔA*_{2}. Using the evolution relation *ΔA*=*vΔt*, we see that both dynamics share a common *Δt*, and thus must have different velocities *v*_{1} and *v*_{2}. The bubble merger velocity *v*_{1}≈2*αAgt* is of the order of the maximum velocity present in the bubble region. As the superposition change in amplitudes *ΔA* and the resulting required superposition velocity *v*_{2} are approximately two times larger, we conclude that the superposition model does not describe the experimental data.

In table 3, we compare superposition and bubble merger velocities to the actual bubble front velocity, i.e. to 2*αAgt*, for five data points taken from four experiments. Each data point involves predicted dynamics from a time *t*=*t*_{j} to an immediately following time *t*=*t*_{j+1}. All data propagation pairs from the five experiments with sufficient numbers of bubbles to support this analysis at each of the steps are included in the sample. For all cases, the comparison clearly favours the bubble merger model, while the superposition model fails to explain the experimental data. For two of experiments (nos. 56 and 63) and for predictions made between more distant (, etc.) time steps, the same method of analysis occasionally favours the superposition dynamics. For the two experiments (nos. 105 and 114) that are free of long-wavelength perturbations, according to the observations or design of the experimentalists or our own directly visible observations (see the discussion of §2*d*), we find that in all cases considered (including prediction over non-adjacent steps) the bubble merger model is favoured and the superposition model is inconsistent with the observed data.

Predictions of self-similar mixing and of *α* [11,14,42] have been based on a model for modal dynamics and power-law initial modal amplitudes. This modal dynamics model assumes superposition (independence of propagation of distinct Fourier modes). It is thus loosely related to the superposition model discussed here. As the superposition model discussed here fails to describe the modal dynamics in a manner consistent with experimental data, especially in cases not deliberately agitated or which do not show visible signs of long-wave perturbations in the early time steps, we regard predictions based on superposition as problematic. In place of such predictions, we prefer the bubble merger model predictions and the front tracking numerical simulations with detailed fidelity to experimental detail, both of which do agree with experimental data.

## 3. Non-universal but approximately universal growth rates

We identified six dimensionless parameters associated with RT mixing that influence *α* as defined within the experimental regime [3,4]. Three of these describe fluid parameters, the Grashof, Schmidt and Prandtl numbers, and three describe initial conditions, the ratios of long- to short-wavelength initial amplitudes, the ratios of short-wavelength determined from the most rapidly growing mode (dispersion theory) to the observed (or inferred) initial short-wavelength perturbation, and to the width of an initial mass diffusion layer. As the level of agreement between simulation and experiment continues to improve, new dimensionless parameters will no doubt be seen to play a role, such as the Reynolds and Atwood numbers. Thus, the observed (experimental) approximate universality of *α* is only that: approximate universality. How does this observed approximate universality arise? A number of these controlling dimensionless parameters take on preferred values in an experimental setting. Usually the Reynolds number and the Grashof number are large, and within the range of present experiments appear to be convergent. Usually, the ratio of the lengths of the short-wavelength initial perturbations to that of the maximally growing wavelength is unity. Although the long-wavelength perturbations may be large, usually (but not always) they have amplitudes very significantly within the linear regime and are quickly suppressed in importance relative to the rapidly growing short-wavelength modes. Usually, the initial mass diffusion layer is as thin as the experimentalist is able to achieve. This leaves the dependence on the Schmidt and Prandtl numbers, which we have shown to be strictly non-zero but not large [3]. Conventional (ILES) RT simulations do not control most of these dimensionless parameters, resulting in significant variation from one ILES simulation to another, not to mention comparison to experiment. Even supposedly exact direct numerical simulations, when not tied to an experimental choice of parameters, face possible run-to-run variation within the experimental time scale of the present analysis or the need to set these parameters in some uniform manner across multiple simulations.

## 4. Conclusion

We have achieved the first systematic validated simulations of the RT mixing rate *α* across a wide range of explicitly defined experimental situations. We use a compressible code applicable to a variety of compressible turbulent mixing problems. We have addressed a concern regarding these simulations, namely the role of long-wavelength perturbations in the initial conditions. Our simulation results are in disagreement with the claim that long-wavelength initial perturbations can affect the *α* values of presently performed experiments by factors of 2 or 3.

The success of these simulations is based on two principles: careful modelling of all physical detail in the experiments, and the use of an effective, high-resolution numerical algorithm. Because *α* is not universal, it is sensitive to both the modelling issues and the algorithmic issues. Regarding the latter, the essential features of our algorithmic strategy are twofold: tracking to control numerical mass diffusion (sharp interfaces or steep gradients), and LES with dynamic subgrid models to account for the effects of the unresolved scales on the resolved ones. A future issue is to determine the degree to which the present results are partially achieved as the quality of these innovations is partially adopted, as well as to find alternative numerical strategies that duplicate or improve on these results.

## Funding statements

This research was supported in part by the US Department of Energy grants nos. DE-FC02-06-ER25770, DE-FG07-07ID14889, DE-FC52-08NA28614 and DE-AC07-05ID14517 and the Army Research Organization grant no. W911NF0910306. The work of D.H.S. was supported by the Department of Energy under contract no. DE-AC52-06NA25396. This research used resources of the Argonne Leadership Computing Facility at Argonne National Laboratory, which is supported by the Office of Science of the US Department of Energy under contract DE-AC02-06CH11357.

## Acknowledgements

This paper has been co-authored by Brookhaven Science Associates, LLC, under Contract no. DE-AC02-98CH10886 with the US Department of Energy. The United States Government retains, and the publisher, by accepting this paper for publication, acknowledges, a worldwide licence to publish or reproduce the published form of this paper, or allow others to do so, for the United States Government purposes. Use of the computing facilities from the Stony Brook Galaxy, the Stony Brook/BNL NewYorkBlue and the ANL Intrepid is gratefully acknowledged. Los Alamos National Laboratory Preprint LA-UR 11-00423. Stony Brook University Preprint SUNYSB-AMS-11-01. It is a pleasure to thank David Youngs for helpful comments and for providing unpublished information related to his experiments.

## Footnotes

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond: non-equilibrium processes from atomistic to astrophysical scales II’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.