This paper discusses a strategy to initialize a two-dimensional (2D) Reynolds-averaged Navier–Stokes model [LANL's Besnard–Harlow–Rauenzahn (BHR) model] in order to describe an unsteady transitional Richtmyer–Meshkov (RM)-induced flow observed in on-going high-energy-density ensemble experiments performed on the OMEGA-EP facility. The experiments consist of a nominal single-mode perturbation (initial amplitude $ a 0 \u2248 10$ and wavelength $ \lambda = 100 \u2009 \mu $m) with target-to-target variations in the surface roughness subjected to the RM instability with delayed Rayleigh–Taylor in a heavy-to-light configuration. Our strategy leverages high-resolution three-dimensional (3D) implicit large eddy simulations (ILES) simulations to initialize BHR-relevant parameters and subsequently validate the 2D BHR results against the 3D ILES simulations. A suite of five 3D ILES simulations corresponding to five experimental target profiles is undertaken to generate an ensemble dataset. Using ensemble averages from the 3D simulations to initialize the turbulent kinetic energy in the BHR model (*K*_{0}) demonstrates the ability of the model to predict the time evolution of the interface as well as the density-specific-volume covariance, *b*. To quantify the sensitivity of the BHR results to the choice of *K*_{0} and the initial turbulent length scale, *S*_{0}, we execute a parameter sweep spanning four orders of magnitude for both *S*_{0} and *K*_{0}, generating a parameter space consisting of 26 simulations. The Pearson's correlation coefficient is used as a measure of discrepancy between the 2D BHR and 3D ILES simulations and reveals that the ranges $8$ $ \u2272 S 0 \u2272 20 \u2009 \mu $m and 10^{9} $ \u2272 K 0 \u2272 10 10$ cm^{2}/s^{2} produce predictions that agree best with the 3D ILES results.

## I. INTRODUCTION

In the context of inertial confinement fusion (ICF), the mixing between initially separated different materials plays an important role in the performance of capsule implosion.^{1} For example, the thermonuclear yield may be reduced when the outer ablator material mixes with the central hotspot of the capsule, thus limiting the efficiency of fusion reactions. Material mixing may originate from the growth of hydrodynamic instabilities, such as the Rayleigh–Taylor (RT),^{2,3} Richtmyer–Meshkov (RM),^{4,5} and Kelvin–Helmholtz (KH)^{6,7} instabilities. These instabilities develop due to small imperfections (perturbations) present at the surface of the materials (e.g., due to manufacturing defects), which may grow under the influence of shear (KH), material accelerations (RT), and/or shock waves (RM). The turbulence (which drives mixing) ensuing from these instabilities and its effect on capsule performance are still an active area of research.

Laboratory experiments, especially under high-energy-density (HED) conditions, face many challenges in studying turbulence in part due to the low repetition rate of the experiments—resulting in only partially converged statistics—and the wide range of temporal and spatial scales the experiments are performed under. Furthermore, due to limitations in the diagnostics, experiments typically only provide integrated [i.e., two-dimensional (2D)] measures of three-dimensional (3D) complex flow dynamics, potentially ignoring important physics and mechanisms. Relying on numerical simulations is, therefore, critical and plays a complementary role to the experiments.

Simulating KH-, RT-, or RM-driven mixing using direct numerical simulations (DNS) can be prohibitively expensive because all the spatial and temporal scales have to be resolved, making approaches such as Reynolds-averaged Navier–Stokes (RANS) often the more favorable engineering option for applications like ICF. To this day, no DNS has been performed for ICF even on the largest supercomputers, as the resolution requirements are too stringent.^{8} However, RANS approaches also face their own challenges: RANS is based on the Reynolds decomposition of a flow where mean quantities are intended to represent an average over an ensemble of realizations, which is often replaced by a spatial average due to the scarcity of ensemble datasets. Replacing ensemble averages by space averages may be appropriate for flows that are in homogenous-, isotropic-, and fully developed turbulent states in which spatial, temporal, and ensemble averaging are often equivalent. However, most HED hydrodynamic experiments involve transitional periods in which the flow is neither homogeneous nor isotropic nor fully developed but may contain large-scale unsteady dynamics; thus, the equivalency of averaging can no longer be assumed. Yet, RANS models often still require to be initialized in such states of turbulence, and knowing how and when to initialize them in a transitional state is, therefore, challenging and is still poorly understood.

The goal of this paper is to develop a strategy allowing the initialization of a RANS model to describe an unsteady transitional RM-induced flow. We seek to examine how ensemble-averaged quantities evolve during the transition to turbulence based on some of the first ensemble experiments repeated under HED conditions. Our strategy involves using 3D high-resolution implicit large eddy simulations (ILES) to supplement the experiments and both initialize and validate the RANS model. We use the Besnard–Harlow–Rauenzahn (BHR) model,^{9–12} specifically designed to predict variable-density turbulent physics involved in flows like RM. Previous studies have considered different ways of initializing the BHR model. In early generations of the model (original model,^{9} BHR2,^{11,13} and BHR3^{10}), initialization is done in a heuristic manner, in which the turbulent kinetic energy is initially set to be a fraction of the laminar kinetic energy of the flow—although the meaning of the chosen fraction remains unclear. Later, Rollin and Andrews^{14} developed a modal model for the RT instability based on the analytic work of Goncharov,^{15} in which non-linear evolution of RT modes are tracked and used to initialize the BHR model at a time the flow is assumed to have transitioned to turbulence. To the author's knowledge, this approach has only been reported by Braun and Gore^{12,16} for the classical RM and RT instability^{16} and reversed-gravity RT,^{12} the latter introduced in the latest version of the model (BHR4). Other approaches that do not involve heuristic guesses or a modal model have also proven useful. For example, Haines *et al.*^{17} have leveraged the use of ILES simulations to validate the BHR model on an inverse chevron shock tube experiment^{18} and showed that reasonable agreement is obtained when the initial value of the turbulent length scale, *S*_{0}, is of the same order of magnitude as that of the standard deviation of the surface roughness of perturbations. Haines *et al.*^{19,20} have applied this methodology to laser-driven shear experiments and showed reasonable agreement with the experiments. Flippo *et al.*^{21} have also demonstrated experimentally that *S*_{0} is related to the initial surface roughness. Although the above-mentioned studies have laid down a path for BHR initialization strategies, most of these studies focus primarily on RT or KH instabilities and are conducted in conditions relevant to traditional fluid mechanics. Thus, it is not apparent whether these strategies apply to the RM instability and under HED conditions.

In this work, we perform ensemble 3D high-resolution ILES simulations of on-going RM HED experiments performed on the OMEGA-EP laser facility^{22} (the ModCons experiments). One the unique features of these experiments is that they are repeated multiple times to generate an ensemble dataset. Although the repetition rate is low—this campaign has taken ∼20 shots over 2 years—these experiments are one of the first ensemble datasets generated under HED conditions.

This paper is organized as follows: first, we describe the ModCons experiments and show the differences between each target in initial conditions and late-time experimental radiographs in Sec. II. Next, we show results from the 3D simulations by first describing the flow evolution from one individual simulation and then considering the ensemble suite of simulations to investigate the order by which the 3D fields are reduced to 2D fields in Sec. III. Using the 3D simulations, we initialize 2D simulations using the BHR model, which are compared to the 3D results in Sec. IV. To evaluate the results from the 2D simulations, a sensitivity analysis with respect to the model parameters *S*_{0} and *K*_{0} is performed in Sec. IV C, and conclusions are drawn in Sec. V.

## II. ENSEMBLE HED HYDRODYNAMIC INSTABILITY EXPERIMENTS

The ModCons experiments are one of the experimental campaigns designed by Los Alamos National Laboratory in an effort to expand our current understanding of turbulence and mixing in HED conditions. Originally designed to study the RM instability in single-mode^{23} and multi-mode^{24} configurations, it has also served as a modeling platform for combined instabilities^{25–27} and laser-induced time-dependent accelerations.^{23} The supporting experiments in this paper focus on single-mode RM.

A typical ModCons target is shown in Fig. 1(a) and consists of a cylindrical layer of polyamide-imide (PAI of density 1.45 g/cm^{3}) adjacent to a layer of foam (density 0.1 g/cm^{3}), both of which having a 1100 *μ*m diameter. For visualization purpose, a 100 *μ*m thick strip of iodine-doped plastic is inserted in the PAI adjacent to the interface between the PAI and the foam. The iodinated plastic is preferentially opaque to the x ray used to radiograph the system. Ripples at the interface are then machined into the PAI and serve as perturbation profile described by a single-mode sinusoidal wave with wavelength $ \lambda = 100 \u2009 \mu $m (wavenumber $ k = 2 \pi / \lambda $) and initial amplitude $ a 0 = 10 \u2009 \mu $m such that $ y ( x ) = a 0 \u2009 sin ( k x )$, as shown in Fig. 1(b) from the top view. The system is then irradiated from the brim of the PAI (called the ablator) by three OMEGA-EP laser beams depositing 14.3 kJ of energy, launching a shock wave with a post-shock pressure $ \u2273 10$ Mbar. The shock arrives at the interface at $ \u223c 13$ ns after the lasers turn on, which causes the RM-instability to grow. At $ \u223c 10$ ns, the lasers turn off causing an expansion wave to travel into the system, which reaches the interface at $ \u223c 19$ ns. Therefore, the system is RM-unstable alone for $ \u223c 6$ ns. When this expansion wave reaches the interface, it decompresses and decelerates the interface, making the system RT-unstable,^{23,26} such that both RM and RT contribute to perturbation growth past $ \u223c 19$ ns.

A unique feature of the on-going ModCons experimental campaign is the ability to repeat the experiments multiple times (of the order of 20) generating some of the first HED RM ensemble datasets used to calculate average quantities. Although each realization of the experiments is characterized by the same nominal perturbation profile (single-mode sinusoidal wave described in the previous paragraph), they differ between each other in the random character of the surface roughness imposed by the machining process. Figure 2 shows an example of a rendered perturbation profile over one wavelength [Fig. 2(a)] obtained from point coordinates of one experimental target (target 1) along with its associated spectral content [Fig. 2(b)]. For comparison purpose, Fig. 2(b) also includes the spectral content of four other targets. The spectral content of each target is obtained by taking a Fourier transform of the perturbation profile, $ y \u0302$, averaged over two different directions: the first is in the direction into the page [z direction, solid lines in Fig. 2(b)] and the second is in the spanwise direction [x direction, dashed lines in Fig. 2(b)]. As expected, the Fourier transform for the z-averaged profile shows the presence of a dominant mode (low-wavenumber high-amplitude) corresponding to the nominal perturbation, along with higher harmonics (high-wavenumber small-amplitude) corresponding to the surface roughness. By design of the experiments, all targets have the same dominant mode but slightly different higher harmonics due to the randomness of the surface roughness. The Fourier transform of the x-averaged profile also has a dominant mode corresponding to the mean location of the perturbation but its amplitude is different for each target due to slight variations in the location of the mean interface. The amplitudes of the higher harmonics are similar to those in the z-averaged direction.

The primary diagnostic for the experiment is line-integrated absorption radiography, and only one high quality radiograph can be taken per shot.^{23} The images are timed the same on every shot, and an ensemble dataset of density maps is obtained at that time, where each image is a slightly different instantiation of the system evolution due to physical variations in target fabrication, like surface roughness. Figure 3 demonstrates this time-matching procedure over six different realizations of the same experiment. The resolution of the density map is limited by the line-integration through the target (CHI strip thickness), the resolution, and the motion blur. Density fluctuation statistics can still be obtained out the images above a certain scale,^{28} but these diagnostics limitations motivate us to use ILES simulations to examine other averaged variables or ways of averaging.

## III. 3D ILES SIMULATIONS

The purpose of using RANS approaches in applications like ICF is to be able to predict results in an engineering-relevant timeframe that would otherwise be impractical to obtain with approaches like DNS. RANS models typically perform averages over a fully turbulent state, where spatial, temporal, and ensemble averages are often equivalent. During the transition to turbulence, however, this equivalency is not straightforward and we seek to examine how averaged quantities evolve with different types of averaging during the transition. To this end, ensemble 3D ILES simulations are performed, which allows us to access both ensemble and spatial averages (we do not consider the temporal average since the flow is in a transitional regime). Moreover, the 3D ILES simulations allow for BHR-relevant parameters (such as the initial turbulent kinetic energy) to be directly computed, which are then used to initialize 2D BHR simulations (Sec. IV). The ensemble set consists of five 3D ILES simulations which simulate the RM-only portion of perturbation growth, i.e., for 6 ns after the arrival of the shock at the interface.

### A. Problem setup and numerical methods

*μ*m

^{2}cross section is selected, as depicted by the red square in Fig. 1(b). The computational domain is shown in Fig. 4, where an initial perturbation separates a layer of PAI from a layer of foam. The initial perturbation is generated from the point coordinates of the selected experimental cross section obtained from a stereolithography (STL) file and then reconstructed into a closed surface, as shown in Fig. 2(a). Table I indicates the number of point coordinates used to reconstruct the closed surface for each of the five targets. The amplitude of the initial perturbation, defined as the halfway distance between peak and valley, is $ a 0 \u2248 10 \u2009 \mu $m and its wavelength $ \lambda = 100 \u2009 \mu $m. All simulations are performed with the radiation-hydrodynamics code xRAGE,

^{29}which solves the multi-material compressible Euler equations of fluid motion as follows:

*ρ*is the material density,

*U*(

_{i}*i*= 1, 2, 3) is the three-dimensional material velocity,

*p*is the pressure, $ E \u2261 1 2 U i 2 + e$ is the total energy (with

*e*being the specific internal energy),

*Y*is the mass fraction of species

_{n}*n*, and

*δ*is the delta function. These equations assume that body forces, viscous stresses, and heat flux are negligible (further justification is given in the next paragraph). Although the materials are initially in the solid phase, they vaporize upon shock passage and are characterized using a modeled equation of state corresponding to polystyrene (at density 1.45 g/cm

_{ij}^{3}for the PAI and 0.1 g/cm

^{3}for the foam) obtained from the SESAME database.

^{30}Although the shock in the experiments is driven by laser irradiation onto the PAI, simulating such physics (e.g., via a laser model) in three dimensions is computationally prohibitive. Instead, the shock is artificially driven by setting up inflow conditions obtained from a one-dimensional precursor simulation with a laser model, giving the post-shock density $ \rho s = 4.5$ g/cm

^{3}, post-shock pressure $ p s = 13.5$ Mbar, and post-shock material velocity $ u p = 24 \u2009 \mu $m/ns. This procedure has been described in detail previously by Di Stefano

*et al.*

^{24}and is, therefore, not repeated here; however, Appendix A shows profiles of the density, pressure, and velocity obtained from this precursor simulation. In all our simulations (except for the 1D simulation of Appendix A), time is taken relative to the arrival time of the shock on the interface and not relative to lasers turn-on.

Target number . | Target 1 . | Target 2 . | Target 3 . | Target 7 . | Target 9 . |
---|---|---|---|---|---|

Number of point coordinates | 146 675 | 141 746 | 136 639 | 130 320 | 133 631 |

Target number . | Target 1 . | Target 2 . | Target 3 . | Target 7 . | Target 9 . |
---|---|---|---|---|---|

Number of point coordinates | 146 675 | 141 746 | 136 639 | 130 320 | 133 631 |

The assumptions made to arrive at Eq. (1)—viscous effects, heat flux, mass diffusion, and laser-related physics (such as radiation and three-temperature effects) are negligible—are justified by the fact that the Reynolds number (*Re*) and Péclet number (*Pe*) are high and that laser-related physics are important mainly when simulating the laser-energy deposition on the target. To further justify these assumptions, we performed a supporting 2D simulation of the entire target with a laser model^{31–33} (data not shown), and concluded that effects of the lasers on the evolution of the bulk hydrodynamics are minimal; including such effects would not change the results of this paper. Concerning viscous effects, we estimate the Reynolds number $ R e \u223c 10 5$ based on the post-shock velocity (24 *μ*m/ns), initial perturbation amplitude (10 *μ*m), and a plasma kinematic viscosity (0.014 cm^{2}/s). The latter was previously used in HED KH experiments^{34,35} with similar materials than the ones of the present study. Concerning mass diffusion, the Péclet number—the ratio of advective mass transfer to diffusive mass transfer—is estimated as $ P e \u223c 10 5$ based on the study of Bender *et al.*^{36} who simulated an RM flow also with similar materials than the ones of the present study.

The system of Eq. (1) is discretized using a second-order Godunov finite-volume scheme on a cartesian grid generated by an adaptive-mesh-refinement procedure. The latter is chosen to have four levels of refinement, where level 0 corresponds to the cell size away from the interface (set to $ 1 \u2009 \mu $m) and level 3 (which is the minimum cell size in the domain) corresponds to the cell size at the interface (set to 1/8 *μ*m). Inflow and outflow boundary conditions are imposed at the entrance and exit of the shock tube, respectively, and the four boundaries on the sides are reflecting. Each 3D simulation employs approximately 10^{8} cells and runs approximately for 10^{4} CPU-hours over 180 cores. A fourth-order Runge–Kutta scheme is used to march the equations in time with a varying time step of the order of 0.1 ps.

### B. Individual simulation

Figure 5 shows the time evolution of the perturbation for an individual target at three different times. Because the shock travels from a heavy to a light material, the initial interface undergoes a phase inversion, i.e., the troughs of the perturbation become peaks and vice versa. The type of phase inversion (direct,^{37} indirect,^{38} or marginal^{25,26}) is determined by the relative speed of the post-shock interface with respect to the speed of the incident shock. In our case, since the latter is approximately the same as the speed of the post-shock interface, the initial interface undergoes a marginal phase inversion. As a result, the interface is relatively flat immediately after the passage of the shock, as shown in Fig. 5(a). Growth of the small scale roughness is evident at that time. As the nominal portion of the perturbation grows into a large bubble-and-spike structure—where the spike corresponds to the heavy material penetrating the light one and the bubble to the light material penetrating the heavy one—the perturbations associated with the small scale roughness also grow into their own (but smaller) bubble and spike, as depicted by the arrows in Fig. 5(b). Eventually, large KH roll-ups form along the part of the interface between the nominal bubble and spike, as seen in Fig. 5(c).

Figure 6 compares the time evolution of the amplitude of the nominal perturbation between the 3D simulations and the experiments of Di Stefano *et al.*^{23} and shows good agreement between the two. Error bars of $ \xb1 3 \u2009 \mu $m in amplitude and ±0.5 ns in time are added to reflect the uncertainty made during the measurements. The agreement between the simulation and the experiments gives us confidence that the bulk hydrodynamic parameters (such as the shock strength and bubble-and-spike evolution) can be reproduced even with our artificial initialization of the shock.

*E*, which is reduced to a 1D function of

*κ*by averaging the values of

*E*contained in the cylindrical strip defined by the radius

*κ*, length

*N*(where

_{y}*N*is the number of points in the

_{y}*y*direction), and thickness $ r i \u2264 \kappa \u2264 r i + 1$, where

*r*is a sample radial coordinate with $ i = 1 : N x$ (

*N*being the number of points in the

_{x}*x*direction). Figure 7 shows the TKE spectrum as a function of the radial wavenumber at different times. At all times, the lowest wavenumber (i.e., largest wavelength) mode dominates; this mode corresponds to the large nominal perturbation present at all times and growing in a bubble-and-spike structure. The high-wavenumber modes ( $ \kappa \u2273 10 \u2009 \mu $m

^{−1}) do not change the TKE levels significantly over time, as evidenced by the merging of all the curves into one, indicating that lower wavenumber modes are the ones responsible for TKE increase in the flow. Indeed, there is an amplification of the TKE levels by all modes contained in the range $ 1 \u2272 \kappa \u2272 10 \u2009 \mu $m

^{−1}, although this amplification is significantly reduced past $ t \u2273 2$ ns for modes in the range $ 5 \u2272 \kappa \u2272 10 \u2009 \mu $m

^{−1}, indicating an increase in TKE levels due to lower and lower wavenumber modes. These modes are representative of the formation of the large KH roll-ups in the flow. Although the flow is neither homogeneous, isotropic, nor fully turbulent, the theoretical framework of homogeneous isotropic turbulence is often used as a reference (in particular the Kolmogorov power law $ \kappa \u2212 5 / 3$). This reference is added to Fig. 7 (solid black line) and shows that part of the spectrum adopts the $ \kappa \u2212 5 / 3$ trend over time, indicating the presence of an inertial subrange in which the large vortical structures associated with the KH roll-ups break down and transfer their kinetic energy to smaller eddies. It is interesting to note that the spectrum suggests that a dissipation range may be also present even though no physical viscosity is present (since only the Euler equations are solved); we attribute this dissipation to numerical viscosity introduced by the ILES numerical framework. In Appendix B, we show that although numerical viscosity is present, the ILES framework provides a reasonable description of the flow. The trends shown in Fig. 7 reveal that the flow is indeed in a transitional state but that there is evidence of an evolution toward the Kolmogorov scales of developed turbulence, thus further motivating the need to understand the evolution of averaged quantities (which are used in BHR) in a transitional state.

### C. Ensemble simulations

*b*, relevant to BHR. We consider two cases: one case consists of space averaging first then ensemble averaging (referred to as space-ensemble and denoted $ \u27e8 \u27e8 \xb7 \u27e9 s \u27e9 e$) while the other case consists of ensemble averaging first then space averaging (referred to as ensemble-space and denoted $ \u27e8 \u27e8 \xb7 \u27e9 e \u27e9 s$). The quantity

*b*is a measure of density fluctuations in the flow and is a pertinent quantity to choose because the density field is the most accessible data from the experiments. It is defined as

*f*), and the brackets represent either a space average (denoted by $ \u27e8 \xb7 \u27e9 z$) or an ensemble average (denoted by $ \u27e8 \xb7 \u27e9 N$), where

*N*is the total number of realizations. By definition,

*b*is non-negative and its values range between 0 and 1, where

*b*= 0 signifies that the flow is fully mixed. The procedure for calculating

*b*constitutes a series of steps that are dependent upon the type of averaging used.

#### 1. The space-ensemble case: $ \u27e8 \u27e8 \xb7 \u27e9 s \u27e9 e$

The space-ensemble case corresponds to the method of averaging closest to the averaging method performed in the experiments. Figure 8 shows the steps followed to obtain a 2D field of *b* from the 3D fields of the density. The first row shows the 3D density field for each target, which is first averaged in space in the z direction, leading to a 2D field in the (*x*, *y*) plane for each individual instantiation (step 0, second row), similar to line-absorption radiography. The ensemble average is then taken leading to a 2D ensemble-averaged density field, $ \rho e ( x , y )$ (step 1). The large bubble and spike of the nominal profile are similar for each target but the small scale structures are different. Ensemble averaging conserves the shape of the nominal perturbation but results in smearing the small-scale structures. In the limit $ N \u2192 \u221e$, these small-scale structures would disappear and only the nominal perturbation would remain.

The third row shows the 2D field of the density fluctuations, $ \rho \u2032$, obtained by subtracting $ \rho e ( x , y )$ from the density field for each instantiation: $ \rho \u2032 ( x , y ) = \rho ( x , y ) \u2212 \rho e ( x , y )$ (step 2). The same is performed on the specific volume to obtain $ v \u2032$.

The fourth row shows the 2D field of the product $ \u2212 \rho \u2032 v \u2032$ (step 3), which is only positive because regions of positive $ \rho \u2032$ correspond to regions of negative $ v \u2032$ and vice versa.

The resulting 2D field of *b*(*x*, *y*) is obtained by ensemble averaging the fourth row and is shown in the last column of the fourth row (step 4).

#### 2. The ensemble-space case: $ \u27e8 \u27e8 \xb7 \u27e9 e \u27e9 s$

In the ensemble-space case, all operations are performed on the 3D fields until the last step to obtain a 2D field of *b*. As in the space-ensemble case, Fig. 9 shows the steps followed to obtain a 2D field of *b* from the 3D fields of the density. The first row shows again the 3D density field for each target, which is now ensemble averaged in 3D first, leading to a 3D ensemble-averaged density field, $ \rho e ( x , y , z )$ (step 1).

The second row shows the 3D field of the density fluctuations obtained by subtracting $ \rho e ( x , y , z )$ from the density field for each instantiation: $ \rho \u2032 ( x , y , z ) = \rho ( x , y , z ) \u2212 \rho e ( x , y , z )$ (step 2). The same is performed on the specific volume to obtain $ v \u2032$.

The third row shows the 3D field of the product $ \u2212 \rho \u2032 v \u2032$ (step 3). Ensemble averaging the third row leads to a 3D field of *b* (step 4), which is then space averaged in the z direction to obtain a 2D field of *b* (step 5, fourth row, last column). There is also another way to obtain this 2D field of *b*, which is to space average the 3D product $ \u2212 \rho \u2032 v \u2032$ for each target, which leads to 2D fields of *b* for each target, which are then ensemble averaged to obtain the final 2D field of *b*. This case is encircled in a dashed box in Fig. 9.

Each step is summarized in Table II for both methods of averaging. Only the steps for *ρ* are outlined but the same procedure follows for *v*. The steps are included in red in Figs. 8 and 9. The 2D fields of *b* obtained from both methods of averaging are compared in Fig. 10 and shows that both methods lead to similar results. To measure their differences, we calculate the root mean square (RMS) of the difference between the *b* obtained from the space-ensemble case and the one obtained from the ensemble-space case, which gives $ RMS \u2248 0.01$. This low RMS value indicates that the method of averaging performed on the 3D ILES simulations is not important, i.e., it does not affect subsequent comparison with 2D-BHR simulations. This result shows two things: first, when compared to 2D-BHR simulations, one method is not better than the other but both are valid comparison, and second, that both methods connect the 3D ILES simulations to a 2D representation valid for the experiment. Note that this low RMS value is obtained only with five instantiations; it is interesting to conjecture that the RMS would tend to zero as the number of instantiations increases.

. | Quantity calculated . | Space–ensemble $ \u27e8 \u27e8 \xb7 \u27e9 s \u27e9 e$ . | Field dimension . | Ensemble–space $ \u27e8 \u27e8 \xb7 \u27e9 e \u27e9 s$ . | Field dimension . |
---|---|---|---|---|---|

Step 0 | ρ | $ \u27e8 \rho [ x , y , z , N ] \u27e9 z \u2192 \rho [ x , y , N ]$ | 2D | – | – |

Step 1 | ρ _{e} | $ \u27e8 \rho [ x , y , N ] \u27e9 N \u2192 \rho e [ x , y ]$ | 2D | $ \u27e8 \rho [ x , y , z , N ] \u27e9 N \u2192 \rho e [ x , y , z ]$ | 3D |

Step 2 | $ \rho \u2032 = \rho \u2212 \rho e$ | $ \rho \u2032 [ x , y , N ] = \rho [ x , y , N ] \u2212 \rho e [ x , y ]$ | 2D | $ \rho \u2032 [ x , y , z , N ] = \rho [ x , y , z , N ] \u2212 \rho e [ x , y , z ]$ | 3D |

Step 3 | $ \u2212 \rho \u2032 v \u2032$ | $ \u2212 \rho \u2032 [ x , y , N ] v \u2032 [ x , y , N ]$ | 2D | $ \u2212 \rho \u2032 [ x , y , z , N ] v \u2032 [ x , y , z , N ]$ | 3D |

Step 4 | $ \u2212 \u27e8 \rho \u2032 v \u2032 \u27e9$ | $ \u27e8 \u2212 \rho \u2032 [ x , y , N ] v \u2032 [ x , y , N ] \u27e9 N \u2192 b [ x , y ]$ | 2D | $ \u27e8 \u2212 \rho \u2032 [ x , y , z , N ] v \u2032 [ x , y , z , N ] \u27e9 N \u2192 b [ x , y , z ]$ | 3D |

Step 5 | – | – | – | $ \u27e8 b [ x , y , z ] \u27e9 z \u2192 b [ x , y ]$ | 2D |

. | Quantity calculated . | Space–ensemble $ \u27e8 \u27e8 \xb7 \u27e9 s \u27e9 e$ . | Field dimension . | Ensemble–space $ \u27e8 \u27e8 \xb7 \u27e9 e \u27e9 s$ . | Field dimension . |
---|---|---|---|---|---|

Step 0 | ρ | $ \u27e8 \rho [ x , y , z , N ] \u27e9 z \u2192 \rho [ x , y , N ]$ | 2D | – | – |

Step 1 | ρ _{e} | $ \u27e8 \rho [ x , y , N ] \u27e9 N \u2192 \rho e [ x , y ]$ | 2D | $ \u27e8 \rho [ x , y , z , N ] \u27e9 N \u2192 \rho e [ x , y , z ]$ | 3D |

Step 2 | $ \rho \u2032 = \rho \u2212 \rho e$ | $ \rho \u2032 [ x , y , N ] = \rho [ x , y , N ] \u2212 \rho e [ x , y ]$ | 2D | $ \rho \u2032 [ x , y , z , N ] = \rho [ x , y , z , N ] \u2212 \rho e [ x , y , z ]$ | 3D |

Step 3 | $ \u2212 \rho \u2032 v \u2032$ | $ \u2212 \rho \u2032 [ x , y , N ] v \u2032 [ x , y , N ]$ | 2D | $ \u2212 \rho \u2032 [ x , y , z , N ] v \u2032 [ x , y , z , N ]$ | 3D |

Step 4 | $ \u2212 \u27e8 \rho \u2032 v \u2032 \u27e9$ | $ \u27e8 \u2212 \rho \u2032 [ x , y , N ] v \u2032 [ x , y , N ] \u27e9 N \u2192 b [ x , y ]$ | 2D | $ \u27e8 \u2212 \rho \u2032 [ x , y , z , N ] v \u2032 [ x , y , z , N ] \u27e9 N \u2192 b [ x , y , z ]$ | 3D |

Step 5 | – | – | – | $ \u27e8 b [ x , y , z ] \u27e9 z \u2192 b [ x , y ]$ | 2D |

## IV. 2D BHR SIMULATIONS

Having highlighted the differences between space-averaging and ensemble-averaging of the 3D ILES simulations, we now validate 2D BHR simulations against the 3D ILES simulations. The initial perturbation profile in the 2D simulations corresponds to the nominal perturbation of the 3D profile, i.e., $ \lambda = 100$ and $ a 0 = 10 \u2009 \mu $m without the surface roughness. Such initial profile corresponds to the profile obtained if the ensemble average of the full 3D profile (i.e., with the surface roughness) was taken over an infinite amount of realizations.

### A. BHR initialization

*f*, as in a Reynolds decomposition, is decomposed in a mean and fluctuating parts, $ f = f \u0303 + f \u2032$, but where $ f \u0303$ is a density-weighted average, i.e., $ f \u0303 = \u27e8 \rho f \u27e9 N / \u27e8 \rho \u27e9 N$, where $ \u27e8 \xb7 \u27e9 N$ represents the ensemble average. Decomposing all variables in Eq. (1) this way and averaging these equations over an ensemble of realizations leads to a set of averaged transport equations that the BHR model evolves.

^{10}In practice, three parameters are required to initialize the BHR model: the initial TKE,

*K*

_{0}, the initial turbulent length scale,

*S*

_{0}, and the time at which the model is “activated,”

*t*. These parameters are evaluated by directly computing their values from the 3D ILES simulations. The TKE is defined as

_{BHR}First, the ensemble average of the 3D density field is computed, giving $ \u27e8 \rho \u27e9 N$. Then, $ U \u0303 i$ is computed by taking the product of the 3D density field with the 3D velocity field for each target, $ \rho U i$, taking the ensemble average of the product ( $ \u27e8 \rho U i \u27e9 N$), and dividing by $ \u27e8 \rho \u27e9 N$ to give $ U \u0303 i$. The 3D velocity fluctuations are then computed for each target via Eq. (5); then, the triple product $ \rho u i \u2033 u i \u2033$ is also computed for each target and ensemble averaged to finally give a 3D field of the TKE.

Since *K*_{0} is a global quantity (i.e., it only takes one value), the 3D TKE field must be reduced to a single value, which is done by computing its volumetric average in the entire field, $ K \xaf = ( 1 / \u222c d x dydz ) \u222c K \u2009 dxdydz$, and its time evolution is shown in Fig. 11. The TKE sharply increases between the initial time and $ t \u2248 1$ ns due to the passage of the shock over the interface (increasing velocity fluctuations) and then decreases with some slight oscillations. The fact that the TKE has a peak value indicates that most of the TKE is produced during the interaction of the shock with the interface, suggesting that turbulent motion becomes significant only after $ t \u2248 1$ ns. Therefore, the model “activation” time is chosen to be *t _{BHR}* = 1 ns and the peak value $ K \xaf = 2.8 \xd7 10 9$ cm

^{2}/s

^{2}is chosen as initial value for

*K*

_{0}. Concerning

*S*

_{0}, Haines

*et al.*

^{17}suggested that

*S*

_{0}is related to the standard deviation of the initial material interface. Using the point coordinates of the initial 3D surface from an STL file, the standard deviation is about 0.2

*μ*m, which is the value chosen for

*S*

_{0}. Note that if an ILES simulation is not available, strategies using potential flow models have been previously adopted,

^{12,14,39}in which the evolution (amplitude and velocity) of initial modes present along the interface is tracked using linear perturbation analysis or vorticity-based equations. These potential flow models are run until an appropriate time typically determined by a critical Reynolds number, above which the flow may be turbulent.

The values $ S 0 = 0.2 \u2009 \mu $m, $ K 0 = 2.8 \xd7 10 9$ cm^{2}/s^{2}, and *t _{BHR}* = 1 ns constitute our baseline case and is discussed below. The model is “activated” in the entire domain with constant values for

*S*

_{0}and

*K*

_{0}, and the other model parameters (such as the turbulent mass flux and the density-specific-volume covariance) are initialized at zero.

### B. The baseline BHR simulation

Our starting reference in the validation process is the case when $ K 0 = 2.8 \xd7 10 9$ cm^{2}/s^{2}, $ S 0 = 0.2 \u2009 \mu $m, and *t _{BHR}* = 1 ns. Figure 12 compares the density field, the

*b*field, and the TKE field between the 2D BHR simulations and the 3D ILES simulations at different times. As in Fig. 9, the 3D ILES simulations are first ensemble averaged over all realizations and then spatially averaged in the z direction to obtain a 2D field in the (

*x*,

*y*) plane. Initially, the interface is sharp in the 2D simulations but has some thickness in the 3D simulations, as seen on the density field, which is the result of ensemble averaging only over five realizations. The

*b*value in the 2D simulation is initialized at zero but is close to one over the thickness of the interface in the 3D simulation (i.e., the materials are segregated), reflecting the fact that the location of the mean interface is slightly different for each realization. The TKE has the initial value of

*K*

_{0}in the 2D case and zero in the 3D case (as seen on the logarithmic colorbar). Slightly after shock passage over the interface (

*t*= 1.1 ns), the interface has flattened (with growth of the small scale roughness in the 3D case) and an expansion wave is reflected into the PAI, as seen by the high-density values ahead of the wave. Both

*b*and the TKE increase in the 3D simulation across the interface and the shock due to an increase in density and velocity fluctuations associated with the surface roughness. In the 2D simulation,

*b*is still close to zero everywhere and the TKE increases but is approximately two orders of magnitude lower than in the 3D simulation. At

*t*= 3.7 ns, the bubble and spike are clearly visible in both the 2D and 3D simulations with a non-uniform distribution of

*b*along the interface: higher values near the spike than along the rest of the interface, as seen in the 3D simulation. These higher values are again the result of the low number of realizations in the ensemble average and are not visible in the 2D simulation. The TKE in the 3D simulation also shows spatial variation along the interface but with relatively similar magnitude on each half of the domain. Evidence of shock refraction patterns near the interface is also visible but not in the 2D simulation. At

*t*= 6.0 ns, the roll-ups have now formed, causing an increase in both

*b*and the TKE in the roll-up regions, as seen in the 3D simulation. Although this increase is slightly visible in the 2D simulation, the magnitude is lower by approximately a factor of two for

*b*and an order of magnitude for the TKE.

From this comparison between the 2D and the 3D simulations, the 2D-BHR simulation overall predicts well the behavior of the nominal perturbation at all times, both in terms of the shape of the interface and the growth of the bubble and spike. The interface in the 2D simulation stays smooth over time whereas growth of the small scale structures is shown in the 3D simulation (due to the relatively low number of realizations in the ensemble average). Concerning the parameter *b* and the TKE, the agreement between the 2D and 3D simulations is not as good: the 2D simulation predicts a thinner area over which *b* and the TKE are distributed, lower values in magnitude, and does not capture the high-*b* regions near the spike seen in the 3D simulation. To understand the sensitivity of these results with the initial values of *K*_{0} and *S*_{0}, we report next on a parameter sweep with respect to *K*_{0} and *S*_{0}.

### C. Parameter sweep with respect to *K*_{0} and *S*_{0}

Figure 13 shows the density, the *b* parameter, and the TKE at *t* = 6 ns in the parameter space (*S*_{0}, *K*_{0}) from the 2D simulations. For comparison purpose, the fifth column shows the solution from the 3D simulation (ensemble averaged first then space averaged in the *z* direction) and is repeated in each row. Increasing *S*_{0} or *K*_{0} has the general trend of diffusing the interface. Increasing *K*_{0} too much tends to diffuse the interface so much that the roll-ups are no longer visible (e.g., when $ S 0 = 10 \u2009 \mu $m and $ K 0 = 10 12$ cm^{2}/s^{2}), while increasing *S*_{0} too much tends to underpredict the magnitude of *b* and the TKE (e.g., when $ S 0 = 100 \u2009 \mu $m and $ K 0 = 10 9$ cm^{2}/s^{2}). The high-*b* regions near the spikes are best predicted when $ S 0 = 1 \u2009 \mu $m and $ K 0 \u2265 10 9$ cm^{2}/s^{2} but the interface is too thin compared to the 3D case. Increasing *S*_{0} to $ S 0 = 10 \u2009 \mu $m thickens the interface to thicknesses closer to the 3D case but the peak values near the spike are then reduced. Concerning the TKE, its distribution along the interface as well as its order of magnitude for $ S 0 = 10 \u2009 \mu $m and $ K 0 = 10 10$ cm^{2}/s^{2} relatively corresponds to the ones for the 3D case. These qualitative observations suggest that the ranges $ 1$ $ \u2264 S 0 \u2264 10 \u2009 \mu $m and 10^{9} $ \u2264 K 0 \u2264 10 12$ cm^{2}/s^{2} may be the most appropriate for comparison with the 3D results. To further refine these ranges, we conducted an additional six simulations in a refined region of the parameter space corresponding to $ S 0 = 3 , \u2009 5 , and 7 \u2009 \mu $m and $ K 0 = 10 11 and 10 12$ cm^{2}/s^{2}, and are shown in Appendix C. For all three values of *S*_{0}, the case $ K 0 = 10 12$ cm^{2}/s^{2} results in an interface that is too diffuse and the peak value of *b* near the spikes is not reproduced. For $ K 0 = 10 11$ cm^{2}/s^{2}, all three values of *S*_{0} give similar results and would be appropriate values to use.

*b*variable:

*M*is the total number of pixels in the

*b*image, and $ b \xaf$ is the mean value of

*b*taken over the entire image. Since the

*b*-fields are 2D but the summation 1D, the 2D images are flattened into a 1D array. The Pearson's coefficient represents the ratio of the covariance between the 2D and 3D datasets, $ cov ( b 2 D , b 3 D )$, to their respective standard deviations, $ \sigma b 2 D$ and $ \sigma b 3 D$. Its values range from –1 to 1 with 0 meaning the two datasets are uncorrelated. When the coefficient is 1 or –1, the two datasets are linearly correlated. Figure 14 shows the correlation coefficient in the parameter space (

*S*

_{0},

*K*

_{0}) at

*t*= 6.0 ns and reveals that values closest to 1 lie in the region $ 8 \u2009 \mu $m $ \u2272 S 0 \u2272 20 \u2009 \mu $m and 10

^{9}cm

^{2}/s

^{2}$ \u2272 K 0 \u2272 10 10$ cm

^{2}/s

^{2}. Note that this range of values for

*S*

_{0}slightly differs from the values found from the qualitative analysis above, but the order of magnitude is the same.

*S*

_{0}and

*K*

_{0}separately. We consider two quantities as follows:

*S*

_{0}(with fixed

*K*

_{0}) and

*K*

_{0}(with fixed

*S*

_{0}), respectively.

Figure 15 shows that $ S 0 = 1 \u2009 \mu $m and $ K 0 = 10 9$ cm^{2}/s^{2} give the closest match to the 3D ILES simulations since the behavior of the mean TKE widely deviates from the 3D simulations for the other values of *S*_{0} and *K*_{0}. Figure 16 shows that the evolution of *b* disagrees by similar factors for all values of *S*_{0} but the trend is relatively following the one from the 3D simulations. This disagreement in the magnitude of *b* comes from initializing *b* = 0 in the 2D-BHR simulations, which is not the case as demonstrated by the non-zero value of *b* at all times. When varying *K*_{0}, the same observation of wrong initial value but right trend applies, with the best agreement again for $ K 0 = 10 9$ cm^{2}/s^{2}.

This analysis helps us conclude that there exists a case of a 2D-BHR simulation ( $8$ $ \u2272 S 0 \u2272 20 \u2009 \mu $m and 10^{9} $ \u2272 K 0 \u2272 10 10$ cm^{2}/s^{2}) that is a good comparison to the 3D ILES simulations and ultimately to the experiments.

## V. CONCLUSIONS

In this work, we develop a strategy to initialize a 2D RANS model (LANL's BHR model) in the context of on-going HED ensemble experiments of laser-driven RM with delayed RT instabilities. Our strategy leverages 3D ILES simulations to initialize BHR-relevant parameters used in 2D BHR simulations, which are subsequently validated against the 3D ILES simulations.

By performing a suite of five high-resolution 3D ILES simulations corresponding to five experimental target profiles, we are able to provide an ensemble dataset used to evaluate the averaging method performed in the experiments, the latter taking an ensemble average of measured 2D mix quantities. We show that the order by which our 3D results are reduced to 2D fields (i.e., either by taking a space average of the 3D fields first followed by an ensemble average, or by taking an ensemble average of the 3D fields first followed by a space average) is not important when calculating mix quantities. This result is key because it shows that a direct comparison between experimental data and subsequent 2D BHR simulations can be achieved without additional constraints for the BHR model.

The availability of a simulated 3D ensemble dataset further allows us to initialize 2D BHR simulations by directly computing the initial parameters needed for BHR (initial TKE, *K*_{0}, initial turbulent length scale, *S*_{0}) from the 3D ILES simulations. By initializing *K*_{0} with the peak value of the TKE at the time immediately after shock passage and *S*_{0} with a value of the same order of magnitude than the standard deviation of the initial surface roughness, good agreement is obtained with the 3D ILES simulations for the shape of the interface but not for the spatial distribution of *b* nor the *TKE*. We show that the latter can be reproduced by BHR by increasing both *S*_{0} and *K*_{0} but within a reasonable range: increasing *S*_{0} too much underpredicts the magnitude of *b* while increasing *K*_{0} too much overdiffuses the interface to the point where roll-up structures are no longer visible. By investigating the distribution of the Pearson's correlation coefficient in the parameter space (*S*_{0}, *K*_{0}) generated by 26 2D BHR simulations, we show that the ranges $8$ $ \u2272 S 0 \u2272 20 \u2009 \mu $m and 10^{9} $ \u2272 K 0 \u2272 10 10$ cm^{2}/s^{2} produce predictions that agree best with the 3D ILES results. This range for *S*_{0} is an order of magnitude higher than the standard deviation of the initial surface roughness, which is in agreement with the previous findings of Haines *et al.*^{17} This parameter sweep with respect to *K*_{0} and *S*_{0} allows us to show that a case of a 2D BHR simulation corresponding to the experiments exist, which gives us confidence in future direct comparison with experimental data.

This work sets a precedent for practical BHR predictions of future HED experiments in which only a few realizations can usually be afforded. Further validation against the ModCons experimental dataset is warranted and will require the implementation of a laser model in 3D to capture the late-time delayed RT physics from laser turn-off. On-going research is tackling this task and will be used in a future study.

## ACKNOWLEDGMENTS

We wish to acknowledge Noah Braun and Brian Haines for helpful conversations about the BHR model and details on numerical methods, as well as Fernando Zigunov for fruitful conversations. Los Alamos National Laboratory, an affirmative action/equal opportunity employer, is operated by Triad National Security, LLC for the National Nuclear Security Administration of U.S. Department of Energy under Contract No. 89233218CNA000001.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Sam Pellone:** Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). **Alexander Martin Rasmus:** Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Supervision (equal); Writing – review & editing (supporting). **Carlos Alex Di Stefano:** Conceptualization (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Supervision (supporting); Writing – review & editing (supporting). **Elizabeth C. Merritt:** Conceptualization (supporting); Funding acquisition (lead); Methodology (supporting); Project administration (lead); Supervision (supporting); Writing – review & editing (supporting). **Forrest W. Doss:** Funding acquisition (supporting); Supervision (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: 1D PRECURSOR SIMULATION WITH A LASER MODEL

To give more details about the inflow conditions used in the 3D simulations, we perform a 1D precursor simulation with a laser model to track the evolution of the density, pressure, and velocity induced by the laser-generated shock. The energy deposition by the laser onto the target is modeled using a laser package,^{31–33} which aims to reproducing the 10 ns, square-shaped drive pulse present in the experiments. Figure 17 shows pressure values in the (*t*, *y*) plane from the precursor simulation. Time is taken relative to when the lasers turn on. Initially, material ablation from lasers turn-on launches a shock in the PAI material, as seen by the high-pressure region behind the shock (10 Mbar $\u2272$ pressure $ \u2272 15$ Mbar). After further examination, the pressure immediately behind the shock is $ \u223c 13.5$ Mbar. This shock reaches the PAI-foam interface at ∼13 ns, which causes the interface to be displaced and a shock to be transmitted in the foam material, as seen by the pressure jump in the latter. The interface location is obtained by tracking the position of a passive tracer added at the location of the initial interface. At 10 ns, the lasers turn off, as seen by the sharp decrease in the pressure, which launches an expansion wave in the domain. To obtain further information on the evolution of the 1D system, Fig. 18 shows the evolution of the density [Fig. 18(a)] and the velocity [Fig. 18(b)] profiles at four different times: the initial time (*t* = 0 ns), a time after lasers turn-on (*t* = 5 ns), a time after lasers turn-off (*t* = 12 ns), and a time corresponding to the end of the RM-only phase (which lasts $ \u223c 6$ ns after lasers turn-off, *t* = 16 ns). The location of the PAI-foam interface (referred to as $ y int$) is shown by black dashed lines. Since the shock does not arrive at the interface until *t* = 13 ns, the interface location is the same for all times and corresponds to its initial location expect for *t* = 16 ns. Initially, the density profile corresponds to the three-material system used to initialize the simulation: the PAI, the foam, and a background vacuum modeled as a near-zero density material. This background region is $ \u223c 4300 \u2009 \mu $m thick behind the PAI and serves as a region in which the laser drive is initialized. At this time, the velocity is zero everywhere since the lasers have not turned on yet. After the lasers turn on, a shock travels into the PAI with a tail of decreasing density corresponding to the ablation region of the PAI material, as seen at *t* = 5 ns. This tail is also visible on the velocity profile with large negative values corresponding to the ablated material traveling in the direction opposite of the incoming laser beams. Note the magnitude of the density and velocity of the shocked PAI material: $ \rho s \u2248 4.5$ g/cm^{3} and $ u s \u2248 24 \u2009 \mu $m/ns, respectively. At *t* = 12 ns, there is a change of slope in the density and velocity profiles in the tail region between $ 300$ $ \u2272 y \u2272 380 \u2009 \mu $m, corresponding to the expansion wave generated by the lasers turn-off at *t* = 10 ns. This expansion wave travels in the same direction as the shock, as shown by the positive values of the velocity. At *t* = 13 ns, the shock arrives at the interface, causing a shock to be transmitted into the foam and an expansion wave to be reflected into the PAI. At *t* = 16 ns, this traveling shock and expansion are evident.

### APPENDIX B: ILES NUMERICAL VISCOSITY

*ν*is the kinematic molecular viscosity and

*ν*is the numerical viscosity. From

_{n}*ν*, an effective Reynolds number can then be calculated $ R e eff = U L / \nu eff$, where

_{eff}*U*and

*L*are characteristic velocity and length scales, respectively. An estimation of the numerical viscosity can be performed

^{40–42}based on the dissipation rate of the TKE,

*ε*, and the strain-rate tensor, $ S i j \u2261 1 2 ( \u2202 U i \u2202 x j + \u2202 U j \u2202 x i )$,

*D*is a dimensionless ratio approaching a constant value when the Reynolds number is large enough ( $ D \u2248 0.5$). Following the methodology of Zhou

*et al.*,

^{42}

*L*is taken to be the thickness of the mixing zone (defined in terms of mixedness), and

*U*is of the order of the root mean square turbulence intensity, $ U = ( 2 \u27e8 TKE \u27e9 3 ) 1 / 2$. To compare

*ν*to the dissipation associated with the BHR model, we consider the turbulent viscosity based on the values of the initial TKE (

_{eff}*K*

_{0}) and the initial turbulent length scale (

*S*

_{0}), $ \nu t = 0.28 \u2009 S 0 K 0$. Figure 19 shows the time evolution of the effective viscosity and Reynolds number [Fig. 19(a)], along with a map of the turbulent viscosity in the parameter space (

*S*

_{0},

*K*

_{0}) [Fig. 19(b)]. The effective Reynolds number toward the end of the simulation is $ R e \u223c 10 5$, indicating that the effective value of the viscosity is not as significant and the ILES approximation is reasonable. Around 1 ns (when the shock impacts the interface), there is a rapid decrease in the Reynolds number to $ R e \u223c 10 2$, suggesting that numerical viscosity may have more of an impact on the solution around that time, and reflecting the dissipative nature of shock-capturing schemes in the ILES framework. When comparing the numerical viscosity to the turbulent viscosity [Fig. 19(b)], the latter is multiple orders of magnitude higher than the former, reflecting the fact that RANS models are more dissipative than an ILES description.

### APPENDIX C: REFINED PARAMETER SWEEP: $ S 0 = 3 , \u2009 5 , \u2009 7 \u2009 \mu $m, $ K 0 = 10 11 , \u2009 10 12$ Cm^{2}/s^{2}

Figure 20 compares the 2D density field, *b* field, and TKE field to the 3D simulations for the refined values of *S*_{0} and *K*_{0}.

## REFERENCES

*Proceedings of the Ninth International Workshop on the Physics of Compressible Turbulent Mixing*

*Implicit Large Eddy Simulation*