{'fraction_non_alphanumeric': 0.03109827150903088, 'fraction_numerical': 0.01582831617789862, 'mean_word_length': 4.68492961634005, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A new model for mixing by double-diffusive convection (semi-convection). II. The transport of heat and composition through layers

Regions of stellar and planetary interiors that are unstable according to the Schwarzschild criterion, but stable according to the Ledoux criterion, are subject to a form of oscillatory double-diffusive (ODD) convection often called"semi-convection". In this series of papers, we use an extensive suite of three-dimensional (3D) numerical simulations to quantify the transport of heat and composition by ODD convection, and ultimately propose a new 1D prescription that can be used in stellar and planetary structure and evolution models. The first paper in this series demonstrated that under certain conditions ODD convection spontaneously transitions from an initially homogeneously turbulent state into a staircase of convective layers, which results in a substantial increase in the transport of heat and composition. Here, we present simulations of ODD convection in this layered regime, we describe the dynamical behavior of the layers, and we derive empirical scaling laws for the transport through layered convection.


INTRODUCTION
A long-standing problem in modeling stellar and planetary interiors is how to quantify the transport of heat and composition in regions that have a destabilizing (i.e., superadiabatic) thermal stratification and a stabilizing compositional stratification, such that the total density stratification is stable according to Ledoux's criterion (Schwarzschild & Härm 1958). Such regions are present, for instance, in giant planets (Stevenson 1982;Leconte & Chabrier 2012), and in massive stars both main-sequence (Merryfield 1995) and post-main-sequence (Robertson & Faulkner 1972;Langer et al. 1985;Woosley et al. 2002). It has long been recognized that, because the molecular diffusivity of composition is typically only a tiny fraction of the total (radiative plus conductive) diffusivity of temperature, these regions are subject to a form of double-diffusive instability (Stern 1960;Walin 1964;Veronis 1965;Kato 1966). The turbulent convection that arises from this instability, which is often called "diffusive" or "oscillatory" convection, enhances the transport of heat and composition, changing the internal structure of the host object.
Oscillatory double-diffusive convection is closely related to the phenomenon known as "semi-convection," but the latter is often also associated with the effects of convective overshoot and the dependence of thermal diffusivity on composition, neither of which are included in this study. Instead, we consider a localized region away from tsw25@soe.ucsc.edu any radiative-convective boundary, in which the thermal and compositional diffusivities can be approximated as constant, and we study the instability and turbulence arising from the competition between thermal and compositional stratification. Following Paparella et al. (2002) we refer to convection in this situation as "oscillatory double-diffusive (ODD) convection." The linear instability mechanism that leads to ODD convection can be understood from physical arguments. In the absence of thermal and compositional diffusion, vertically displaced parcels of fluid oscillate adiabatically about their equilibrium position with a frequency determined by the total density stratification. But the presence of diffusive processes can amplify the oscillation, producing an instability that resembles an overstable internal gravity wave (Kato 1966). Indeed, if the thermal diffusivity is sufficiently large in comparison with the compositional diffusivity, then a displaced parcel will acquire heat (and therefore entropy) while below its equilibrium position, and thus will be slightly more buoyant than its surroundings when it returns to its equilibrium position. The amplitude of the oscillation grows exponentially until nonlinear effects are sufficient to saturate the instability, leading to a state of weakly turbulent ODD convection.
After saturation, ODD convection initially takes the form of a homogenously turbulent state, in which the timeaveraged stratification remains uniform. Previous work by Rosenblum et al. (2011) and Mirouh et al. (2012) has shown that the turbulent transport of heat and composition in this homogeneous state is fairly weak. Mirouh et al. (2012, hereafter Paper I) derived empirical formulae for the respective fluxes from an extensive suite of numerical simulations performed in local Cartesian domains. However, in certain cases the homogeneous state was found to spontaneously transition into a new state of layered convection, for which the transport of heat and composition is greatly enhanced. Paper I derived a simple semi-analytical criterion that predicts under which conditions spontaneous layering occurs. Under the physical conditions relevant to planetary and stellar interiors, spontaneous layering is expected within a significant fraction of the double-diffusive instability range, specifically in regions that are closer to being Ledoux-unstable (see Figure 1 in Paper I).
In this paper, we present further results from numerical simulations and use them to derive empirical formulae for the fluxes of heat and composition in layered convection. We determine how the fluxes depend on the thickness of the layers and on other governing parameters. We also describe some of the physical properties of layered convection, paying particular attention to the interfaces between the layers, and we compare our results with previous studies of ODD convection. Finally, we discuss our results in the light of recent layered models of giant planets.

THE MATHEMATICAL MODEL AND NUMERICAL SCHEME
The numerical model that we use to study ODD convection was described in detail in Paper I, so here we mention only its most relevant features. The computational domain is a Cartesian box with uniform gravity g = −ge z . We solve the evolution equations for momentum, temperature, and composition in the Boussinesq approximation, which are ∂u ∂t w is the vertical component of the velocity u; all of these quantities are assumed to be periodic in all three Cartesian directions. We use subscript '0's to denote the unperturbed, background fields, which we assume are functions of z alone. The background gradients of temperature, T 0z , and mean molecular weight, µ 0z , are both negative in the situation considered here, as is the adiabatic temperature gradient, We consider a small volume within which the diffusivities of momentum, ν, temperature, κ T , and composition, κ µ , can be treated as constant, and the density perturbations can be expressed as a linear function of the temperature and composition perturbations, Both α and β are positive for any realistic equation of state. We approximate α, β, and the other coefficients in Equations (1)-(3) as constant within our small domain, and we suppose that the background temperature and molecular weight profiles are linear in z, so that the total temperature and composition fields are T tot = T + T 0z z and µ tot = µ + µ 0z z, with T 0z , T ad 0z , and µ 0z all constant. The Schwarzschild and Ledoux stability criteria are then and respectively. We are concerned here with regions which satisfy (9) but not (8), i.e., regions in which R −1 0 > 1, where R −1 0 is the "inverse density ratio," From a linear stability analysis of Equations (1)-(3), it can be shown that such regions are unstable to double-diffusive (Walin 1964;Kato 1966) where R −1 c is the critical value Here, Pr and τ are the diffusivity ratios and In astrophysical objects, both Pr and τ are typically of order 10 −2 or smaller, and so the instability window given by Equation (11) is rather wide. The instability is oscillatory, and has a typical horizontal lengthscale (e.g., Baines & Gill 1969). measure these fluxes in terms of Nusselt numbers, defined as and where the notation · represents a volume average over the entire domain. The total vertical fluxes of heat, F T , and composition, F µ , can be reconstructed from these as follows: and where c p is the specific heat. We note that, in general, the thermal Nusselt number Nu T defined by Equation (16) is not the ratio of the total heat flux to the diffusive (i.e., radiative plus conductive) heat flux, unless the adiabatic temperature gradient T ad 0z is zero. However, as in Paper I, we can introduce the total "potential temperature" which has background gradient ϑ 0z = (T 0z − T ad 0z ) and perturbation ϑ = T , in order to write the thermal Nusselt number as So we may regard Nu T as the ratio of the total and diffusive fluxes of potential temperature. Since ϑ and T are equal, in what follows we will refer to ϑ as simply "the temperature perturbation".
In double-diffusive convection, an important quantity is the "buoyancy flux ratio" γ −1 , which we define as (Radko 2003) 1 . This is closely related to the core erosion "efficiency factor" χ introduced by Guillot et al. (2004), which measures the fraction of a planet's luminosity that is used to erode the stable compositional stratification at the boundary with the core. In the limit where the transport of both heat and composition is dominated by convection (i.e., Nu T,µ ≫ 1) this fraction is precisely γ −1 . The gravitational potential energy released from the thermal stratification must exceed the potential energy used to erode the compositional stratification, which implies that γ −1 < 1. As shown in Paper I, γ −1 plays an important role in the formation of layers, and possibly also controls the evolution of layers after their formation (Radko 2005).
For computational convenience, we work mainly with nondimensional quantities. As in Paper I we nondimensionalize and The governing equations are then 1 Pr where all variables are now dimensionless. The definitions (16)-(17) of the Nusselt numbers then become and We will also make use of alternative definitions of the Nusselt numbers, and which are based on the rates of thermal and compositional dissipation, and are sometimes called "Cox numbers".
We note that nu T and nu µ are guaranteed to be positive at all times, whereas Nu T and Nu µ are not. However, it can be readily verified, from Equations (28)-(30), that the definitions (31)-(32) and (33)-(34) yield the same result after taking a time average in a statistically steady state (Malkus 1954;Shraiman & Siggia 1990). The advantage of retaining both forms is explained in Section 4.1.

TYPICAL PROPERTIES OF LAYERED CONVECTION
In principle, any particular instance of ODD convection is characterized by the three dimensionless parameters R −1 0 , Pr, and τ . However, in Paper I it was shown that, within the instability window given by Equation (11), there is another critical value of R −1 0 , say R −1 L , below which ODD convection spontaneously transitions into layered convection. Paper I presented empirical scaling laws for R −1 L , and for the transport of heat and composition prior to the development of layers. It was shown that, in the parameter regime relevant to astrophysical objects, for which Pr ≃ τ ≪ 1, the critical value R −1 L ∼ Pr −1/2 . This result implies that there is a significant range of density ratios R −1 0 for which layers form spontaneously, and so the transport through layered convection must be considered in models of astrophysical objects. 2 A preliminary study of layered convection by Rosenblum et al. (2011) demonstrated that the turbulent transport is strongly dependent on the layer height, H say, in addition to Pr, τ , and R −1 0 . Their results suggest that Nu T and Nu µ both follow a power law in H, with an exponent between 1 and 4/3. In what follows we present a more comprehensive study of the dependence of the turbulent fluxes on each of the parameters H, Pr, τ , and R −1 0 . In order to provide a clear description of the transport properties of layered convection, and the subtleties involved in measuring the thermal and compositional fluxes in the layered phase, we first revisit one particular simulation that was previously presented in Paper I. This simulation has R −1 0 = 1.5, Pr = 0.03, τ = 0.03, and a domain size of (100d) 3 . As described in Paper I, this simulation spontaneously transitions from homogeneous ODD convection to layered convection, and the three layers that initially form subsequently merge into two layers, then into one layer. This evolution is illustrated in Figure 1, which shows snapshots of the composition perturbation µ during each phase of the simulation, as well as profiles of the horizontally-averaged total "potential density," The profiles show that ρ tot is roughly constant within each layer, but increases slightly with height. The local Ledoux stability criterion can be expressed as ∂ρ tot ∂z < 0 (cf. Equation (9)), so this indicates that the layers are subject to overturning convection. Within each layer, the total composition µ tot and the total potential temperature ϑ tot are both roughly uniform, and there are strong jumps in µ tot and ϑ tot across the interfaces between adjacent layers. with R −1 0 = 1.5, Pr = τ = 0.03, and a domain size of (100d) 3 . The first column shows the homogeneous phase of ODD convection prior to layer formation. The remaining columns show layered convection with 3, 2, and 1 layers respectively. Because the perturbation is plotted in the top row, rather than the total composition, the layers appear as regions with a relatively weak vertical gradient, and the interfaces appear as regions with a strong vertical gradient of the opposite sign. The color bar for each plot is scaled to the maximum of the perturbation µ, which increases as the layers merge. In the profiles, ρtot has been rescaled such that the total change in ρtot across the domain is unity.
Spontaneous layer formation similar to that observed here has previously been reported in simulations of doublediffusive convection in the "fingering" regime, for which the thermal stratification is stabilizing and the compositional stratification is destabilizing (e.g., Stellmach et al. 2011), but only for Prandtl numbers Pr > 1 ).
However, in fingering convection the initial height of the layers is typically much larger than the lengthscale of the linear double-diffusive instability, whereas here we find that layers initially form on a scale comparable to that of the linear instability. In the simulation shown in Figure 1, for example, the fastest growing linear mode has a wavelength of roughly 20d, i.e., 1/5 of the domain size, and the initial layer height is 1/3 of the domain size. Figure 2 shows more detail of the interface structure during the 2-layer phase. We find that the interfaces are highly dynamic and turbulent, exhibiting structure on scales that are much smaller than the wavelength of the fastest growing linear mode. This is in marked contrast to the flat, laminar interfaces found in simulations and laboratory experiments of ODD convection at larger Prandtl number (e.g., Noguchi & Niino 2010a;Carpenter et al. 2012, both of which have Pr ≃ 7). The composition field exhibits structure on finer scales than the temperature field, as expected because of the smaller diffusivity of composition. On larger scales, the temperature and composition fields display similar features, such that the temperature field resembles a coarse-grained version of the composition field. By contrast, the density field has even finer structure, since it depends on the difference between the composition and temperature fields. Regions where the perturbations are small have been made transparent in order to show the interfaces in more detail.
The fluxes through this simulation, measured in terms of Nusselt numbers Nu T,µ and nu T,µ , are plotted in Figure 3 as a function of time, along with the mean-squared velocity, |u| 2 . As mentioned earlier, the mean value of each flux increases significantly when layers first form, and increases still further when layers merge (Rosenblum et al. 2011). The mean kinetic energy also increases, roughly in proportion with the fluxes, indicating that the simulation becomes more turbulent as the layers become thicker. We note that the two definitions for the Nusselt numbers both yield consistent results for the mean flux, although the definition based on the instantaneous fluxes exhibits larger fluctuations in the homogeneous, pre-layer phase of the simulation, during which the dynamics are dominated by wave-like, oscillatory motions (Mirouh et al. 2012).
Between layer merger events the mean fluxes and kinetic energy oscillate quasi-periodically between high and low values, with relative departures from the mean flux of order unity. The oscillations occur as the flow in the simulation cycles between states that are highly turbulent and states that are more quiescent. In the quiescent state, the interfaces are roughly horizontal, whereas in the turbulent state the interfaces have undulations on a horizontal scale comparable to the thickness of the layers. For illustration, Figure 4 shows the composition perturbation, and the horizontallyaveraged total composition µ tot , at two times during the one-layer phase. At the first time, t = 1980, the interface is in a quiescent state, and at the second, t = 2020, the interface is in a more turbulent state. These two times correspond to minima and maxima in the Nusselt numbers, as shown in Figure 3.
The quasi-periodic nature of the oscillations in the fluxes suggests that they are caused by large-scale gravity-wave motions of the interfaces. For an interfacial gravity wave of horizontal scale H, the oscillation frequency is  The initial linear instability grows exponentially until t ≈ 570, then saturates, leading to homogenous ODD convection. Layers first form at t ≈ 1200; the duration of the 3 layer, 2 layer, and 1 layer phases are indicated at the top of the plot. The vertical lines indicate the two times shown in Figure 4. Using nu T,µ rather than Nu T,µ to measure the mean fluxes has the advantage of filtering out fluctuations in the fluxes on short timescales, while retaining the important long timescale information. (Rayleigh 1883), where ∆ρ tot represents the jump in potential density across the interface. For a staircase of layers of thickness H we have ∆ρ tot = Hρ 0 |βµ 0z − αϑ 0z |, and so the characteristic gravity-wave frequency is where N is the buoyancy frequency associated with the mean stratification, which in dimensionless units is N = Pr(R −1 0 − 1). In the simulation shown in Figure 3 this corresponds to a (dimensionless) oscillation period of 2π/ω ≃ 50, which is approximately half of the typical period of the observed flux oscillations. We emphasize, however, that the disturbances to the interfaces seen in Figure 4 are far more complex in structure than linear gravity waves. Moreover, the large-scale motions of the interfaces are not strictly periodic in time -they grow and decay cyclically, but they do not oscillate about zero.
Previous models of layered convection in astrophysical objects (e.g., Spruit 1992) have assumed that the interfaces are horizontal, laminar, and diffusive. If that were the case, the fluxes though the interfaces could be determined simply by measuring the vertical gradients of the horizontally-averaged T tot and µ tot profiles in the center of each interface. Figures 3 and 4 show that this type of model cannot be valid here, because the fluxes are maximal when these gradients are weakest. In fact, since the interfaces are not horizontal, the fluxes depend on the local gradients within the interfaces and on the effective surface area of the interfaces. Our aim is to determine how the mean fluxes of heat and composition, measured in terms of dimensionless Nusselt numbers, depend on the dimensionless parameters R −1 0 , Pr, τ , and the average (dimensionless) layer height, which we denote as H. As described in the last section, during each layered phase the fluxes of heat and composition oscillate quasi-periodically about a mean value, and reliably calculating the mean fluxes requires time averaging over several of these oscillations. We measure the mean fluxes and estimate the errors in our measurements as follows: 1. We first identify the number of layers during each phase of the simulation, by counting the number of "steps" in the vertical profiles of the horizontally averaged fields, and matching each merger event with the corresponding increase in the Nusselt numbers.
2. Once we have identified the time interval corresponding to a particular layered phase, we then divide this interval into four equal subintervals, and calculate the average values of Nu T , Nu µ , nu T , and nu µ during each subinterval using Equations (31)-(34). Thus we obtain eight separate measurements of the thermal and compositional fluxes during each layered phase.
3. We treat each of our eight measurements as independent, and calculate their mean and standard deviation, which yields estimates for the mean flux and the measurement error. (A similar method was used to estimate the measurement error in Paper I.) 4. In the same fashion, we calculate the buoyancy flux ratio γ −1 for each quarter interval using Equation (22), and for both definitions of the Nusselt numbers. As before, we use the resulting eight measurements to estimate the mean value of γ −1 and our measurement error.
We have applied this procedure to each of the simulations from Paper I in which layers formed spontaneously, as well as to many additional simulations that have since been performed. The parameters of each simulation, including those from Paper I that formed into layers, are listed in Table 1.
We present below the results of these measurements and describe their dependence on the parameters Pr, τ , R −1 0 , and H. In principle these four parameters are sufficient to characterize any particular regime of layered convection.
However, in our numerical simulations we also specify the size of the computational domain and the numerical resolution, and so we must check whether these additional parameters have any effect on the results. Appendix A describes some of the tests that we have performed to ensure that our measurements are robust. For reasons explained in the Appendix, only the flux measurements from simulations performed in domains of width 100d are presented below.
We also note that the spontaneous emergence of layers described in Paper I is not the only mechanism by which layers can develop in ODD convection. Laboratory experiments and observations of Earth's oceans indicate that the formation of layers can be "triggered" by external influences (Turner & Stommel 1964;Kelley et al. 2003). Layered convection might even be metastable under conditions for which the linear ODD instability is absent (Veronis 1965;Spiegel 1972). In this paper, however, we only consider the regime in which layers form spontaneously, 1 < R −1 0 < R −1 L , in order to guarantee the reproducibility of our results.

Parameter Dependence of the Mean Fluxes
To facilitate comparison of our results with those of previous studies, we introduce the thermal Rayleigh number for layered convection, which is where Hd is the (dimensional) layer height. Owing to computational limitations, our results cover only a limited range of Pr and τ values. Moreover, by choosing to analyze only those simulations in which layers form spontaneously, we are restricted to a rather narrow range for R −1 0 (see Paper I). On the other hand, the range of Ra T values that can be achieved in each simulation is limited only by the height of the computational domain. Moreover, the fact that layers merge allows us to test the dependence of the Nusselt numbers on Ra T using data from a single simulation. The tallest simulation presented here has a domain of height 400d, and initially formed eight layers. Figure 5 shows the Nusselt numbers in this simulation, plotted as functions of time and Ra T . We find that both Nu T and Nu µ are well approximated by a power law in Ra T with an exponent of 1/3. This result can be explained physically if we suppose that the mean flux depends only on the jumps in potential temperature and composition across the interfaces, ∆ϑ and ∆µ, say, and not on any properties of the layers. More precisely, let us suppose that the (dimensional) buoyancy fluxes gακ T ϑ 0z Nu T and gβκ µ µ 0z Nu µ depend only on the buoyancy jumps gα∆ϑ and gβ∆µ and the diffusivities κ T , κ µ , and ν. On dimensional grounds, the fluxes must then be proportional to (∆ϑ) 4/3 (Turner 1965). In our simulations, the total change in temperature between the top and bottom of the domain is fixed, and so the jump across each interface is inversely proportional to the number of layers, and therefore proportional to the height of the layers. The exponent 1/3 then follows from equation (38).
We now ask whether all of our Nusselt number measurements can be fitted to the same power law. In order to fit all of the data simultaneously, we must make some assumption regarding the dependence not only on Ra T , but also on Pr, τ , and R −1 0 . We consider the thermal and compositional Nusselt numbers in turn, and then the flux ratio γ −1 .  Many analytical and heuristic arguments have been made that the turbulent transport of heat through layered convection follows a power law in both the Rayleigh number, Ra T , and the Prandtl number, Pr (e.g., Spruit 1992;Balmforth et al. 2006). But its dependence on the parameters R −1 0 and τ is less clear. As a first step toward fitting our data to an analytical formula for Nu T , we therefore assume a dependence of the form We then look for a fit of the form The fact that the best fit has exponents a ≈ b suggests that the heat flux depends only on the parameters R −1 0 , τ , and the product Ra T Pr. In that case the heat flux is independent of viscosity, as predicted by Spruit (1992) for layered convection with Ra T ≫ 1 and Pr < 1. However, the model of Spruit (1992) predicts exponents a = b = 1/4, rather than 1/3, which is incompatible with our results.

The Heat Flux (NuT )
Using the result that Nu T − 1 ∝ Pr 1/3 we can directly compare simulations performed with different values of Pr but the same values of τ and R −1 0 . An example is given in Figure 7, which overlays time series of Nu T in three simulations that each have τ = 0.3, R −1 0 = 1.2, and domain size (100d) 3 . For these parameters, the best fit prediction is Nu T − 1 = 0.09H 4/3 Pr 1/3 . Since each of these simulations has the same (dimensionless) domain height, we can directly compare their two-layer (H = 50) and one-layer (H = 100) phases; we find that, as predicted, (Nu T − 1)/Pr 1/3 has the same dependence on H in each of these simulations.  The best-fit values for the coefficients f n , presented in Figure 6, show a weak but well-defined dependence on R −1 0 and τ (for fixed Ra T and Pr). In particular, we find that decreasing either the composition gradient µ 0z or the compositional diffusivity κ µ produces a slight increase in the mean heat flux. We interpret this result as evidence that changing µ 0z and κ µ affects the heat flux only indirectly, by changing the overall level of turbulence.

The Composition Flux (Nuµ)
and fit the data to determine a, b, and f n . This time the best-fit values for the exponents are a = 0.37 ± 0.01 and b = 0.27 ± 0.04. This is roughly consistent with the trend observed in Figure 5b, but indicates that the compositional flux, unlike the heat flux, depends on viscosity even for Pr ≪ 1. This is perhaps not surprising given that the composition field is more sensitive than the temperature field to the small-scale motions of the fluid (e.g., see Figure 2). Figure 8 shows the best fit values for the coefficients f n and a comparison between the measured values of Nu µ and the best fit for Equation (41). As found previously for Nu T , there is a weak increase in Nu µ as R −1 0 is made smaller, i.e., as the composition gradient is reduced. Unlike Nu T , however, Nu µ depends strongly on τ , i.e., on the compositional diffusivity. In fact, we find that Nu µ is roughly inversely proportional to τ , indicating that the mean turbulent flux of composition wµ is approximately independent of compositional diffusivity in the parameter regime considered here.  (41) with exponents a = 0.37 and b = 0.27.

The Buoyancy Flux Ratio (γ −1 )
Taken at face value, the empirical power laws for Nu T and Nu µ derived above imply that the flux ratio γ −1 ∝ Nu µ /Nu T ∝ Ra 0.03 T . However, as discussed in Section 2, γ −1 cannot exceed unity in a statistically steady state, and therefore cannot follow a simple power law in Ra T . As illustrated in Figure 9, we find that γ −1 is generally an increasing function of Ra T , and possibly asymptotes to a constant value at large values of Ra T . Further simulations at high Ra T are required to confirm this trend, and to determine how this asymptotic value depends on the other parameters R −1 0 , Pr, and τ . Figure 9 suggests that γ −1 is a decreasing function of R −1 0 , but again, more data is required to confirm this suggestion.
If γ −1 remains a monotonically increasing function of Ra T for large Ra T then there are important implications for the merging of layers, as we discuss in the next section.

SUMMARY AND DISCUSSION
Ra T 0.6 0.7 0.8 0.9 10 6 10 7 10 8 10 9 10 10 10 11 where the prefactors A and B are weakly dependent on the parameters R −1 0 , τ , Pr, and Ra T , and have typical values A ≃ 0.1 and B ≃ 0.03 respectively. The mean heat and composition flux can be calculated from these expressions using Equations (18) and (19). We emphasize, however, that the parameter values used in our simulations are rather far from true astrophysical values, and so any application of the formulae (42)-(43) to astrophysical conditions requires an extrapolation beyond the range in which these formulae have been derived.
In all our simulations the interfaces between the convective layers exhibit complex dynamical behavior that is very different from the laminar "interface scouring" and entrainment seen in studies of layered convection at higher Prandtl number (e.g., Linden & Shirtcliffe 1978;Fernando 1989). The structure of the interfaces, and their effective surface area, is highly time dependent, which produces order unity variations in the fluxes of heat and composition during each layered phase. Previous models for the transport by ODD convection have assumed that the interfaces remain flat and laminar, implying that the fluxes are limited by the vertical gradients of the horizontally averaged profiles within the interfaces. Our numerical results demonstrate that transport at low Prandtl number is more efficient than predicted by such quasi-one-dimensional models. We argue that the interface dynamics must be considered in future models for transport under astrophysical parameter conditions. Since the transport by layered convection is strongly dependent on the layer thickness, we must consider what determines this thickness in stellar and planetary interiors. Recently, Leconte & Chabrier (2012) have proposed a model for the transport by ODD convection in Jupiter's interior. They use heuristic arguments to obtain upper and lower bounds for the height of layers in Jupiter, and thereby derive bounds for the thermal and compositional transport. In all of our simulations, however, the layers eventually merge into a single layer with the same height as the computational domain, so our results make no prediction for the height of layers in real objects. Whether simulations performed in taller domains would yield a saturation thickness for layers is unknown. For fingering convection, it has been argued by Radko (2005) that there exists such a saturation thickness. The essence of his argument is that the buoyancy flux ratio γ −1 increases every time layers merge, but that there exists a maximum permitted value of γ −1 , which is attained at a finite layer thickness H = H 0 , say. Layers of thickness H H 0 therefore do not undergo further merging. In our simulations of ODD convection, by contrast, γ −1 appears to asymptote toward a limiting value as H → ∞, so this argument does not apply. (This is in contrast to a statement made by Leconte & Chabrier 2012). An alternative picture for layer merging has been proposed by Noguchi & Niino (2010b), but their model also predicts that the merging of layers proceeds indefinitely.
We note that the Boussinesq approximation assumed in our numerical model breaks down once the layer thickness becomes comparable to the local pressure scale height of the planet or star. In practice, this may be what limits the size of layers in real astrophysical objects. At this point, convection within the layers would be almost indistinguishable from "normal" overturning convection, allowing the usual mixing-length approximation for the transport to be applied.
An important issue that has not been explicitly addressed so far is the timescale for layer merging. Our results suggest that this timescale is an increasing function of layer height, which is in accordance with the theory of layer mergers proposed by Radko (2005). If the timescale for merging in an astrophysical object eventually becomes comparable to the global evolution timescale, then the layers will never reach a saturation thickness, and the evolution of the layers must then be modeled as an explicitly time-dependent process in 1D structure and evolution models. Future papers will study the dynamics of layer mergers in more detail, in order to determine what sets the merging timescale, and how it depends on layer height.
We thank B. R. Ruddick, E. A. Spiegel, and M. G. Wells for useful comments and suggestions. P. G. and T. W.
were supported by funding from the NSF (NSF-0933759). Part of the computations were performed on the UCSC Pleiades supercomputer, purchased with an NSF-MRI grant. Others used computer resources at the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the US Department of Energy under contract DE-AC03-76SF00098.

INFLUENCE OF DOMAIN SIZE AND NUMERICAL RESOLUTION
In this appendix we investigate the influence of numerical resolution and domain size on the initial layer height, the timescale for layer formation and merging, the time-averaged fluxes, and the amplitude of oscillations in the fluxes.

Tests of Numerical Resolution
Because of the disparity between the diffusivities of temperature and composition, double-diffusive convection always features a wide range of dynamical length scales, making this a tough problem to model numerically. Furthermore, as mentioned in Section 3, the emergence of layers is associated with an increase in the overall level of turbulence, and increased spectral energy at all scales in a given simulation . Properly resolving the layered phase therefore requires even greater numerical resolution than resolving the homogenous phase. In order to check that the thermal and compositional fluxes in our simulations are not sensitive to the numerical resolution, we have run several simulations that have the same physical parameters and domain size, but different numbers of grid points in the vertical and horizontal directions.
As an example, Figure 10 compares time series of Nu T,µ and nu T,µ from four simulations performed with different numerical resolutions. All four were performed in domains of size (50d) 3 with the same physical parameters as the simulation presented in Section 3, but with 48 3 , 64 3 , 96 3 , and 128 3 Fourier modes respectively. The simulation with the coarsest resolution (48 3 ) only resolves structures on lengthscales d, but even in this simulation the linear double-diffusive instability is well resolved (recall that the fastest growing linear mode has wavelength ≃ 20d at these parameters). Each simulation undergoes a homogeneous phase, then forms three layers, which later merge into two layers, then one layer. The duration of each layered phase varies between the four simulations, suggesting that the merging of layers is, to some extent, a stochastic process. However, the four simulations produce similar fluxes during each phase, and the variations in each flux about its mean are comparable in all four cases. All four cases have Pr = τ = 0.03, R −1 0 = 1.5, and a domain size of (50d) 3 . The two rows show the Nusselt numbers calculated from the instantaneous fluxes (31) and (32), and the instantaneous dissipation rates (33) and (34).
In each simulation, the most stringent test of the numerical resolution occurs during the final merging event, which marks the start of the one-layer phase. Figure 11 compares the density fields at this point in the lowest and highest resolution simulations shown in Figure 10. The lowest resolution simulation is clearly under-resolved, with significant features at the Nyquist scale, whereas the highest resolution simulation displays no numerical artifacts. Nevertheless, the mean fluxes through these two simulations are approximately equal in the one-layer phase, which suggests that the transport of both heat and composition is dominated by scales d in these simulations. Tests of Domain Size and Aspect Ratio Figure 12 compares time series of the Nusselt numbers from four simulations performed in domains of different heights and widths. Each simulation has Pr = τ = 0.3, R −1 0 = 1.15, and the same numerical resolution, i.e., the same (dimensionless) distance between neighboring grid points. These four simulations exhibit the following features, which we find to be typical: • The cases with wider domains (200d) take longer to form layers than the cases with narrower domains (100d).
• The initial layer height is not very sensitive to the size of the domain, provided that the domain is sufficiently wide. In Figure 12, for example, the initial layer height is 50d in all four cases. However, in other simulations (not shown here, but see Rosenblum et al. (2011) for instance) it has been found that in very narrow domains layers form much earlier, and the initial layer height is smaller. This result indicates that the horizontallyperiodic boundary conditions, which encourage horizontal coherence in the simulation, can artificially accelerate the formation of layers when the domain is sufficiently narrow. For this reason, we only use data from simulations that were performed in domains of width 100d when calculating the mean fluxes in Section 4.1.
• The mean fluxes of heat and composition, for layers of a particular height, are consistent between runs with different domain sizes.
• During each layered phase, the cases with wider domains exhibit smaller variations in the value of each flux about its mean value. This suggests that the variations in the fluxes are horizontally localized, and therefore have less effect on the mean flux in wider domains. Similarly, the simulation in the tallest domain (of height 400d) exhibits the smallest variations in the fluxes during each layered phase, because in this simulation the volumetric mean flux represents an average over a greater number of layers and interfaces. Figure 12 also illustrates a practical difficulty that arises when comparing results from simulations performed in domains of different heights, which is that simulations in taller domains usually undergo more layered phases than simulations in shorter domains, and the length of each layered phase is correspondingly shorter. For example, the simulation with a domain height of 400d passes through phases in which the average layer height is 400d/7 and 400d/6; these states cannot be achieved in the simulations with a domain height of 100d. As a result, the time series of the Nusselt numbers in the taller simulation are less obviously step-like, and the mean flux during each layered phase is more difficult to measure accurately. This also makes it difficult to determine whether the domain height affects the rate at which layers merge. The full time series of both Nusselt numbers from this simulation are shown in Figure 5a.
We find that the time between layer merger events increases as the number of layers decreases, which is in accordance with the model of layer merging proposed by Radko (2005).
Another difficulty associated with measuring fluxes in simulations with multiple layers is that the layers are not necessarily all of the same height. In that case, a complete description of the state of the system would require (at least) measuring the height of each layer, and the jumps in temperature and composition across each interface, potentially introducing many new parameters on which the mean flux might depend. Furthermore, in such cases we would not expect the system to be in statistical equilibrium, and so the mean flux might not be well defined. To avoid dealing with this additional complexity, in this study we assume that, in each layered phase, all the layers have approximately the same height, while recognizing that this assumption introduces an additional source of error into our flux measurements.