Semiconvection: theory

# Semiconvection: theory

H.C. Spruit Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany
July 14, 2019
###### Key Words.:
stars: semiconvection – stars: mixing – convection: double diffusive, convection: thermohaline

A model is developed for the transport of heat and solute in a system of double-diffusive layers under astrophysical conditions (viscosity and solute diffusivity low compared with the thermal diffusivity). The process of formation of the layers is not part of the model but, as observed in geophysical and laboratory settings, is assumed to be fast compared to the life time of the semiconvective zone. The thickness of the layers is a free parameter of the model. When the energy flux of the star is specified, the effective semiconvective diffusivities are only weakly dependent on this parameter. An estimate is given of the evolution of layer thickness with time in a semiconvective zone. The model predicts that the density ratio has a maximum for which a stationary layered state can exist, . Comparison of the model predictions with a grid of numerical simulations is presented in a companion paper.

## 1 Introduction

In stellar evolution, ‘semiconvection’ denotes the situation where a thermally unstable stratification is stabilized against (adiabatic) overturning by a gradient in composition (called solute in the following; typically the Helium concentration, increasing with depth). It was first recognized as a complication in the calculation of stellar structure by R. J. Tayler (1953). Uncertainty whether the stabilizing gradient should be ignored (‘according to Schwarzschild’), included in the condition for overturning (‘according to Ledoux’), or something in between has led to a number of different recipes for mixing of the Helium gradient. The evolution of the star subsequent to the semiconvective phase is sensitive to these differences (e.g. Weiss 1989, Langer et al. 1989, Langer 1991, Stothers & Chin 1994, Langer and Maeder 1995). The fluid mechanics encountered in geophysics in the same case of a thermally unstable stratification stabilized by a stabilizing solute is called double diffusive or thermohaline convection.

The observations show that such a stratification forms a stack of many thin layers, called a ‘staircase’, each consisting of overturning convection sandwiched between stable steps in composition and temperature (Turner & Stommel 1964, Padman & Dillon 1987, Schmid et al. 2010 and references therein). In effect, this is a ‘compromise between Schwarzschild and Ledoux’. Correspondingly, the net transport coefficients (of heat and solute) are intermediate between those of convection and diffusion.

The reason for this layer formation are understood (for references, cf. Spruit 1992, hereafter S92, and Zaussinger & Spruit 2013, hereafter ZS13). The main theoretical contribution in this context has been the work of Proctor (1981), who showed analytically that a finite amplitude layered state exists for conditions when the stratification is still linearly stable: layering is the result of a subcritical instability. His analysis applies to the case of vanishing diffusivity of the solute, and Prandtl number not exceeding , which is the astrophysically relevant case. In this limit, layered states exist whenever the Rayleigh number exceeds the critical value for normal convection, independent of the strength of the stabilizing solute gradient. An energy argument (S92) shows that the layered state in this case can be reached with an initial perturbation of vanishing amplitude as the layer thickness decreases.

Attempts have been made to address the astrophysical problem with direct numerical simulations (Merryfield 1995, Biello 2001, Rosenblum et al. 2011). This encounters two problems: the very high thermal diffusivity in a stellar interior (very small Prandtl number) cannot be matched without some form of approximation. More importantly the quantity of main interest, the effective mixing rate, depends on the thickness of the double diffusive layers formed, a quantity that is not a stable outcome of the simulations.

It makes sense to disentangle the ‘semiconvection problem’ into two parts: on the one hand the physics that determines the thickness of the layers, on the other hand the effective mixing rate for a given layering state. This separation is especially meaningful because observations in geophysics and laboratory experiments show layer thickness to be a slowly changing function of time compared with the overturning times within the layers. In the following we concentrate on the second question, that is, the mixing rate is studied as a function layer thickness.

It turns out that the simple observation of a layer structure consisting of an overturning zone between stable zones is sufficient to derive a predictive model for the effective transport coefficients (Sect. 3). It is sufficiently quantitative to be tested against the results from numerical simulations. The translation of the model to the astrophysical case of a compressible stratification can be done exactly in the limit where the layers are thin compared with the pressure scale height. It is given in sec. 6, and yields an easily implementable prescription for the mixing rate.

In the present treatment, the transition in composition and temperature between neighboring layers has a finite width, as opposed to S92, ZS10, where it was approximated as a step. This generalization turns out to make no significant difference for practical astrophysical application, but is essential for meaningful comparison of the theory with numerical results. Section 7 finally gives a (less quantitative) estimate of the layer thickness and its evolution in time.

The astrophysical term ‘semiconvection’ will be used in the text interchangeably with the terms ‘double diffusive’ or ‘thermohaline’ convection used in laboratory and geophysical literature111Existing nomenclature is somewhat ambiguous. The opposite case of a destabilizing composition gradient in a stable thermal gradient is called ‘saltfingering’, but in geophysics also ‘thermohaline convection’, where our semiconvective case is often called ‘diffusive convection’ (e.g. Schmitt 1994). ‘Double diffusive’ is usually meant to cover both cases..

## 2 Physics of a semiconvective layer

### 2.1 Notation and definitions

As in the above, we will use the term ‘layer’ for an individual double-diffusive step in the semiconvective staircase, and the term ‘solute’ for the stabilizing component (e.g. Helium in Hydrogen). Such a layer consists of a zone of overturning convection between adjacent stagnant zones. In the stagnant zones overturning convection is suppressed by a stable () density gradient, and transport of heat and solute takes place by diffusion.

As illustrated in Fig. 1, depth through the layer is counted from the middle of a stagnant zone. The thickness of the layer as a whole is , that of the stagnant zone . The temperature difference across the layer is , the solute difference . The stabilizing influence of the solute is expressed conveniently in terms of the density ratio, the ratio of density differences caused by temperature and solute differences and across the layer:

 Rρ=βΔS/αΔT, (1)

where , the thermal expansion coefficient, is the relative density decrease per unit temperature and , called the haline contraction coefficient, is the relative density increase per unit of solute concentration. Since density increases with and decreases with , the density differences are of the opposite sign when and both increase with depth (in the direction of gravity) as in our semiconvective case.

The flow in the overturning zone of the layer is driven by a temperature difference () and opposed by a solute difference (). Associated with these is an ‘internal density ratio’, related to the overall density ratio by

 ~Rρ≡β~ΔS/α~ΔT=Rρ~ΔSΔSΔT~ΔT. (2)

Obviously, the values of and depend on where we define the boundaries between stagnant zone and overturning zone. Since the internal layer is actually convecting, the boundary has to be set at a point in the profile where (see also Sect. 3).

The overturning zone is characterized by a Rayleigh number; for convection in a layer of thickness , with temperatures and at top and bottom, the Rayleigh number is

 Ra=gα(Tb−Tt)D3κTν, (3)

where is the acceleration of gravity, the thermal diffusivity, and the (kinematic) viscosity. In our case (Fig. 1), has the value . The critical value for onset of convection is of order Ra (for no-slip boundary conditions). If viscosity is low, (Prandtl number Pr), the heat flux at large Rayleigh number becomes independent of viscosity, (a fact that is used implicitly in the ‘mixing length’ formalism for convection in a stellar interior). At high Ra, a more relevant quantity to characterize the heat flux in this case is the modified Rayleigh number Ra:

 Ra∗=PrRa=gα(Tb−Tt)D3κ2T. (4)

Its square root can be read as the ratio of the thermal diffusion time to the free fall time of the density contrast over the distance .

Let be the (time averaged) heat flux across the layer, and the flux in the absence of convection, i.e. when the temperature profile between and is determined by diffusion only. The Nusselt number is then defined as

 NuT=F/Fd. (5)

Similarly, there is a Nusselt number for the solute flux :

 NuS=FS/FSd. (6)

In the absence of a solute, i.e. for normal (unstratified) laboratory convection, Nu is a function of Pr and Ra only (apart from boundary conditions and geometry). In the case of semiconvection, the Nusselt numbers are functions of two additional parameters characterizing the solute: the density ratio and the Lewis number Le, the ratio of solute diffusivity and thermal diffusivity :

 Le=κS/κT. (7)

### 2.2 Transport coefficients: classical results

The heat flux is determined by the boundary layers at top and bottom. As the overturning flow passes along the boundary, the temperature contrast with the boundary diffuses into the flow, and it is this temperature difference that carries the heat flux. If is the time during which the flow is in contact with the boundary before descending/ascending into the bulk, the depth over which the temperature difference penetrates is . At high Rayleigh numbers, this depth is small compared with the layer thickness . In a simple two-dimensional view of the flow the amount of fluid carrying this difference also scales with . The same flow along the boundary determines how much solute contrast diffuses into the flow. The solute diffusion depth is thus , and the amount of solute flowing into the overturning zone is proportional to and inversely proportional to the flow speed along the boundary.

The ratio of the convective fluxes of solute and heat is thus expected to scale as . On the other hand, the (diffusive) fluxes over the thickness of the overturning layer, in the absence of convection, scale as and . In terms of the Nusselt numbers, this implies, in the limit such that the boundary layers are thin compared with the layer thickness:

 NuS=qLe−1/2NuT, (8)

where is a numerical factor expected to be of order unity. If the heat flux is kept fixed, the effective diffusivity of the solute then scales as

 κs,eff=κSNuS=q(κSκT)1/2NuT. (9)

These are classical scalings, as found in previous analyses of double diffusive convection, e.g. in Turner (1980, 1985). In this form they do not address the dependence of the fluxes on the density ratio R. The argument above implicitly assumes that all the density contrast in the solute boundary layer is carried across with the flow. In reality, only a fraction of the boundary layer can flow across, namely the fraction that is not too buoyant to be carried down by the flow, respectively too heavy to be carried up. Taking this into account (Spruit 1992, hereafter S92), yields a correction to (8):

 NuS=qRρLe−1/2NuT.(NuT≫1) (10)

This scaling applies to the overturning zone; it would be valid for the layer as a whole only if the stagnant zone were absent. More precisely, it assumes that the thickness of the stagnant zone is smaller than the thickness of the boundary layers . This is not always the case (cf. Sect. 5 below). In the following, the presence of the stagnant zone is taken into account explicitly, by consistently taking the distinction between between Nusselt numbers for the overturning zone and those for the layer as a whole into account. This results in a generalization of (10).

#### 2.2.1 ‘Erosion’

To determine the numerical factor in (10) more quantitatively than order unity, the hydrodynamic interaction of the flow with the stagnant zone has to be considered in more detail. This is somewhat beyond the scope of the model to be developed here, but we can identify a process involved, and use this to show that is probably somewhat larger than unity. I will call the ‘erosion factor’.

The unsteady flow in the overturning zone induces perturbations in the stagnant zone. Since it is stably stratified, these take the form of internal gravity waves. These can contain internal structure on length scales (perpendicular to the interface) less than the thickness of the stagnant zone. Diffusion of solute on these length scales increases the amount that has sufficiently low buoyancy to be carried with the overturning flow, hence we may expect . In the absence of a more detailed theory, its value can in principle be used as a fitting parameter, as long as it is not taken larger than order unity. In the following, however, I will ignore this option, and simply set

 q=1. (11)

For astrophysical applications, the effective solute transport will turn out to be so low that tuning the erosion factor by order unity would have little effect anyway. At low Le, this erosion process takes place well inside the thermal boundary layer. It consequently affects only the transport of the solute; its effect on the transport of heat can be neglected.

## 3 Model

The layer of thickness now consists (Fig. 1) of a stagnant zone of thickness , and a overturning zone occupying the remainder of . An incompressible (Boussinesq) fluid is assumed, to be generalized to the compressible astrophysical case in Sect. 6.

In the stagnant zone, we assume that the transport of heat and solute is by diffusion only (i.e. ignoring the possible ‘erosion’ effect discussed above). In the overturning zone of the layer, the flow is approximated as convection as it would take place in the absence of a solute. The presence of the solute affects the flow somewhat, but the amount of solute carried, and hence its influence on the flow, vanish in the limit of low solute diffusivity (S92, Schmitt 1994). Apart from the boundary conditions it sets, the stagnant zone has little effect on the flow inside the overturning zone. In the limit considered here, the effect of the solute on the flow of heat in the overturning zone can thus be neglected.

To describe transport of heat in the overturning zone, we use a fit from laboratory measurements for the Nusselt number as a function of the temperature difference. Since the thickness and temperature difference are different from those of the layer as a whole, the Rayleigh and Nusselt numbers of the overturning zone are distinguished here with a  . With the notation

 ϵ=~ΔT/ΔT,δ=ds/d, (12)

they are related to Ra and Nu by (cf. eq. 4 and the geometry sketched in Fig. 1):

 ~Ra∗=Ra∗ϵ(1−δ)3, (13)
 ~NuT=NuT (1−δ)/ϵ. (14)

In the measurements of Niemela et al. (2000), at Ra up to , the Nusselt number is well fit by a power law of slope 0.309 and amplitude 0.124. The Prandtl number in these experiments is order unity, so Ra and Ra are approximately the same. Combining the above, I approximate the heat flux in the overturning zone as given by

 ~NuT=1+a(~Ra∗−Ra∗c)b,(a=0.124, b=0.309). (15)

Since, as argued in 2.1 the dependence on Prandtl number is expected to be weak below Pr, (15) will be assumed approximately valid for all Pr (this can be checked with numerical simulations, see Zaussinger and Spruit 2013, hereafter ZS13). The constant 1 and the critical Rayleigh number have been added to better approximate the behavior at low Ra.

The buoyancy of the solute is opposite to that of the driving temperature difference, and the flow will only be driven by density differences of the unstable sign. This implies that the solute concentration contrast in the overturning flow, say, is limited by the temperature contrast through

 Rρ~S/ΔS≤~T/ΔT. (16)

At the boundary between the overturning and the diffusive zones, these contrasts are , . Define then the location of this boundary as the point where the equals sign holds in (16), that is, the point where the buoyancy of the solute is just low enough for it to be carried with the thermally driven flow in the overturning zone (cf. discussion above). This yields the relation

 ~ΔSΔS=1Rρ~ΔTΔT. (17)

In the overturning zone, the original scaling (8) applies, i.e. (assuming ), since is the density difference that can just be carried with the flow. At low Ra, should approach unity at the same time as . To accomplish this I adopt a slight modification:

 ~NuS−1=Le−1/2(~NuT−1), (18)

which is now assumed to hold for all .

This completes the definition of the model. It defines the fluxes of heat and solute uniquely as functions of the external parameters.

## 4 Effective diffusivities

In a stationary state, the fluxes of heat and solute are constant with depth through the layer. The fluxes by diffusion in the stagnant zone are equal to the fluxes in the overturning zone. Expressing this in terms of the Nusselt numbers, we first need to write and in terms of the Nusselt numbers for the layer as a whole. For temperature this is given by (14):

 NuT=~NuTϵ(1−δ). (19)

The equivalent relation for the solute is

 NuS=~NuSϵS(1−δ), (20)

where (using 17):

 ϵS≡~ΔS/ΔS=ϵ/Rρ. (21)

The fluxes in the stagnant zone are given by:

 NuT=(1−ϵ)/δ,NuS=(1−ϵS)/δ. (22)

(since carried by diffusion alone). With (18), eqs. (19)-(22) yield two equations for the unknown and :

 (1/δ−1)(1/ϵ−1)=~NuT, (23)
 (1/δ−1)(Rρ/ϵ−1)=1+Le−1/2(~NuT−1), (24)

with given by (15). Eliminating between these two yields

 1−δ−ϵ/Rρ=(1−δ−ϵ)/Q, (25)

where

 Q=RρLe1/2. (26)

Solving this for :

 ϵ=(1−δ)1−Q1−Q/Rρ. (27)

Since must be a positive number, and is larger than 1 (otherwise we would not be ‘below Ledoux’), it follows that in a stationary state as envisaged must be less than 1, or . This is a necessary condition (independent of the Rayleigh number), but not sufficient. The actual, somewhat lower, value of the critical density ratio has to be determined from (19)-(20) by solving for as well as . To see how this solution comes about, consider the left-hand and right-hand sides of (23) separately. They represent the heat flux in the stagnant and the overturning zones, respectively, and in a stationary state they are equal. Fig. 2 shows the two as functions of , with taken from (27) and from (15).

The two intersection points of the curves are potential solutions for a steady state. To see which of the two is the relevant one, consider the slopes of the curves. At the equilibrium point with the smaller value of the flux in the stagnant zone decreases more rapidly with than the flux in the overturning zone. A small decrease of away from the equilibrium would increase the flux in the stagnant zone relative to that in the overturning zone. This would result in a temperature deficit at the boundary between the two, causing the temperature to decrease there. This would reduce the heat flux in the stagnant zone again, the opposite of the assumed perturbation. This equilibrium point is thus the stable one; the intersect at the larger is an unstable point.

The necessary condition for existence of the layered state, also figures prominently in Proctor (1981). In his analysis it shows up as a necessary condition for its validity.

### 4.1 Maximum density ratio

The maximum density ratio can be determined as a function of Ra as the value for which the two curves in Fig. 2 just touch. The result is shown in Fig. 3. At large Ra, (slowly) approaches the value . The behavior of Nu and near is shown in Fig. 4 for an illustrative case.

At the critical density ratio the value of is about 0.2 (cf. Fig. 2), while the Nusselt number reaches a minimum value of the order of a few. At this density ratio, the model predicts a jump from an overturning state () to a purely diffusive state . An example is shown in Fig. 4. The presence of this jump suggests that the behavior of the system near the maximum density ratio needs a closer look than is possible with the present model. This question is explored with numerical simulations in ZS13.

The maximum on can also be read as a minimum on the Rayleigh number. For an observed density ratio of 4 at Le for example, one should not expect to see double-diffusive layering with layer thickness less than corresponds to a Ra of (from Fig. 3). In terms of the reference length defined below (Eq. 42), the layers must have a thickness , for this combination of and Le.

### 4.2 Solute flux

The net solute Nusselt number can be expressed in terms of Nu . Using (18)–(22) and (26)–(27), this yields the simple expression

 NuS=1+1RρLe1/2(NuT−1). (28)

The stagnant zone of finite thickness included here thus leads to the same relation between Nu and Nu as the simpler model in S92 (cf. eq. 10 above). There is a difference, however, in the relation between Nu and Ra, and in the presence of a maximum density ratio (cf. Figs. 3,4).

### 4.3 Asymptotic behavior

The Nusselt number is, with (19), (27):

 NuT=~NuTP (29)

where

 P≡1−Q1−Q/Rρ=1−Le1/2Rρ1−Le1/2, (30)

and is given in terms of by (15). For , and with not close to , the stagnant zone is thin, . as a function of and Le is then, with (14), (27):

 ~Ra∗≈Ra∗P.(Ra∗≫1, Rρ

The solute Nusselt number follows from (28),

 NuS≈NuTRρLe1/2.(Ra∗≫1, Rρ

It is instructive to compare to the thickness of the boundary layers , of the overturning zone. These can be written in terms of the Nusselt numbers Nu and Nu since at the boundaries with the stagnant zone the fluxes are carried by diffusion. For the thermal boundary layer for example this implies

 DT/(d−dS)=1/~NuT. (33)

Using (23), (27):

 ~NuT=(Rρ/P−1)/δ, (34)

hence in the limit :

 δ≈DTd(RρP−1)=DTdRρ−11−Q. (35)

Unless is close to unity, the width of the stagnant zone is thus of the same order as the boundary layers of the overturning zone, and the distinction between the two becomes somewhat academic. At large density ratio, or under the fixed heat flux conditions discussed in the next Section, however, the distinction becomes more significant.

## 5 Fixed heat flux conditions

In laboratory situations and theoretical analyses it is usual to consider the temperature difference or temperature gradient as given and the heat flux or Nusselt number as the object to be determined. In some natural systems, however, the conditions under which double diffusive convection occurs are closer to ones where the heat flux is the given quantity. In the east African volcanic lakes, for example, the heat flux imposed by the influx at the bottom of the lake is probably more of a given than the temperature at the bottom of the lake. The same is the case in semiconvective zones of stars.

The temperature difference across the layers adjusts to an imposed heat flux. Since both the density ratio and the Rayleigh number depend on , neither of these can be used as control parameter any more. This requires a change of perspective on the problem.

The stratifications of both temperature and solute change with time, under the effective transport properties of the semiconvective process. In the limit of low solute diffusivity, the time scale for changes in the solute profile is long compared with the time scale on which the temperature profile adjusts. In this thermally quasisteady state, the solute profile can be taken as fixed. Assume therefore that the mean solute gradient is given. As new control parameters, use the heat heat flux and layer thickness .

To make the imposed heat flux practical as control parameter, measure it with respect to a reference flux that can be expressed in terms of the solute gradient. Take for this the diffusive heat flux that would be present in the linear temperature profile that is just marginally stable against adiabatic overturning (‘Ledoux’). If is the thermal conductivity, the heat flux of this stratification is

 F0=KdT0dz, (36)

and its density ratio is

 Rρ0=βdSdz/αdT0dz. (37)

Since we have assumed that the stratification is marginally stable, . The reference heat flux is thus

 F0=KβαdSdz. (38)

The actual layered state has a Nusselt number Nu and a density ratio ; by definition of the Nusselt number, its heat flux is

 F=NuTKΔTd. (39)

With :

 F/F0=NuTαΔTβΔS=NuT/Rρ. (40)

To replace the Rayleigh number as control parameter, note that it can be written in terms of solute gradient and density ratio as (using 1)

 Ra∗=gβκ2TRρd4dSdz. (41)

The quantity

 d0=(κ2Tgβ/dSdz)1/4=(κT/NS)1/2, (42)

is a characteristic length scale of the problem: the distance over which temperature diffuses on the buoyancy time scale of the stable solute gradient, . The Rayleigh number can then be written as

 Ra∗=1Rρ(d/d0)4. (43)

The Nusselt number in (40) is a function of and Ra, which can be evaluated from Eqs. (15,22,23,24) in Sect. 4 above. Together, eqs. (40) and (43) thus determine Ra and R, and other quantities of interest, as functions of the new control parameters and . An example of the dependence on of Ra, , the Nu’s, and is shown in Fig. 5 for .

### 5.1 Asymptotic dependences for fixed heat flux

The limiting case (corresponding to Ra), has a pleasingly simple form. As expected from the discussion in Sect. 4 (cf. Fig. 3), the density ratio approaches its maximum value Le in this limit. With (40) and (28) the Nusselt numbers approach the value

 NuS≈NuT≈Le−1/2F/F0(d/d0≫1). (44)

For , this is in fact a good approximation also at lower values of , as Fig. 5 (lower left) shows. The corresponding effective solute diffusivity becomes

 κs,eff=NuSκS≈(κSκT)1/2F/F0, (45)

in agreement with the classical ‘geometric mean of diffusivities’ scaling. The thickness of the stagnant zone is related to by eq. (22). In the present limit, , hence

 δ≈NuS−1≈Le1/2(F/F0)−1(d/d0≫1). (46)

The relative thickness of the stagnant zone is thus small for low Lewis numbers or at high heat flux. In the opposite case of modest Le and conditions closer to marginal, it can be a significant fraction of the layer thickness, however. This may be the explanation for the relatively large thickness of the stagnant zones observed in lake Kivu (Schmid et al. 2010). The heat flux measured there corresponds to Nusselt numbers of order 2. A stagnant zone thickness of order 30% is observed, larger than in other natural cases like the double-diffusive steps under the arctic ice sheet (e.g. Timmermans et al. 2008 and references therein). Though the present analysis does not apply directly to this case because the assumption Pr does not hold (Pr for water), the thickness of the stagnant zones in lake Kivu may be an indication that it is actually close to the marginal state to be expected at imposed low heat flux conditions.

#### 5.1.1 Thickness of the stagnant zone

As Fig. 5 shows, the density ratio approaches its maximum in the limit (cf Sect. 4.1):

 Rρ≈Le−1/2→1−Q≪1. (47)

The actual value of as a function of the parameters and is a somewhat complicated expression involving the coefficients in (15) (it will not be needed in the following).

Comparing with the thermal boundary layer thickness of the overturning zone then yields, from (35) :

 δ≈DTdLe−1/2−11−Q. (48)

In this asymptotic case with fixed heat flux, the stagnant zone is thus much wider than the boundary layers of the overturning zone (but still thin compared with the layer thickness ). Its presence manifests itself in the approximate equality of Nu and Nu (eq. 44) as opposed to the original estimate (8).

## 6 Astrophysical conditions

For the compressible gas in a stellar interior, an incompressible approximation (Boussinesq) can still be used for flows on length scales small compared with the pressure scale height and time scales short compared with the sound travel time over this length, provided two factors are taken into account (e.g. Massaguer & Zahn 1980). First, the adiabatic lapse rate has to be subtracted from the temperature difference driving the flows. The modified Rayleigh number for a layer of thickness is now

 (49)

where for the ideal gas with constant ratio of specific heats assumed here. In terms of the logarithmic gradient :

 Ra∗=gd4κ2TH(∇−∇a), (50)

where is the scale height of the gas pressure ,

 H=dz/dlnp. (51)

This length scale is not present in the Boussinesq case. An equivalent of the quantity used in Sect. 5 still plays a role as well, but in the astrophysical context is a more conventional choice as reference length. can then be written as

 Ra∗=f2(dH)4(∇−∇a), (52)

where the dimensionless number ,

 f=(gH3/κ2T)1/2, (53)

typically a large number, is the ratio of the thermal diffusion and free fall time scales over a pressure scale height. The dimensionless layer thickness will be the control parameter replacing above.

A second difference concerns the radiative heat flux, since it is driven by the temperature gradient itself rather than the potential temperature222Note that in the presence of a solute the superadiabaticity is not equivalent to the entropy gradient any more. Instead, it measures the gradient of potential temperature, the temperature on adiabatic displacement, in pressure equilibrium, relative to a given reference level.. This affects the definition of Nusselt number, as well as the choice of control parameter specifying the heat flux. Define a thermal conduction constant (a function of the opacity and the thermodynamic state of the gas), in terms of the the radiative heat flux :

 Fr=k∇. (54)

Following standard notation, represent the total heat flux by the radiative gradient :

 F=k∇r. (55)

will be the astrophysical control parameter specifying heat flux. The ratio of density changes due to composition and temperature changes upon adiabatic displacement in the stratification, i.e. the density ratio, is

 Rρ=∇μ∇−∇a, (56)

where

 ∇μ=dlnμ/dlnp, (57)

and the the mean atomic weight per particle. Define the Nusselt number Nu through

 F=k∇a+NuT k(∇−∇a), (58)

so that it measures the part of the heat flux that is due to (only) the superadiabatic part of the temperature gradient, instead of the total heat flux. This has to be clearly distinguished from the ratio of heat flux to radiative heat flux in the stratification, which one might consider the more logical definition of a Nusselt number. The latter includes, however, a flux which is irrelevant for the convective flow (the radiative heat flux due to the adiabatic part of the stratification); it would contribute even if there is no flow at all333The significance of this distinction is not apparent in some of the numerical work on semiconvection, e.g. Biello (2001).. From (55), (56), (58) we obtain a relation between Nusselt number and density ratio for our astrophysical, imposed heat flux conditions:

 ∇r−∇a∇μRρ=NuT. (59)

To complete the model, a prescription for the Nusselt number as a function of the control parameters is needed. For this, we assume that the layer thickness is small enough that an incompressible approximation is valid, so the results of Sect. 4 can be used. Eqs. (18)–(26) then determine as a function of and . Together with from (56) and from (50), this forms a set of equations for the superadiabaticity as a function of the control parameters and .

### 6.1 Asymptotic results

This is a somewhat implicit algebraic problem, but the limiting case equivalent to in Sect. 5 again has a simple form. This limit corresponds to

 f2(d/H)4=gd4κ2TH≫1, (60)

or

 d≫l0≡(κ2TH/g)1/4, (61)

i.e. large compared to the length on which the thermal diffusion time scale equals the free fall time over a scale height. As before, tends to its maximum, in this limit, so that

 NuT≈∇r−∇a∇μLe−1/2, (62)

while the superadiabaticity follows from (56):

 ∇−∇a≈Le1/2∇μ. (63)

The solute Nusselt number follows as in Sect. 5,

 NuS≈NuT, (64)

so the effective solute diffusivity is

 κSeff=(κSκT)1/2(∇r−∇a)/∇μ, (65)

the same as in the simpler model of S92 and ZS10 (apart from a factor involving radiation pressure). It is independent of the layering thickness , within the range of validity of the approximations,

 l0≪d≪H. (66)

The relative thickness of the stagnant zones is

 δ≈1/NuS=Le1/2∇μ/(∇r−∇a). (67)

For Le the stagnant zone is thus thin compared with the layer thickness (but still significantly thicker than the boundary layers of the overturning zone, cf. Sect. 5.1.1).

## 7 Evolution of the layer thickness

The main uncertainty of any model for semiconvection is the thickness of the double diffusive layers. In a stellar interior the layer thickness has little influence on the quantities of interest (eqs. 63, 65), but it can become important for the structure of the star when layer thickness becomes macroscopic, i.e. comparable with the pressure scale height.

In geophysical and laboratory cases, it is observed that evolves secularly, on a long time scale compared with the thermal diffusion time. Layer thickness is therefore a quantity that cannot be discussed independently of the history of the system. The process has not been studied very extensively (but see Wirtz & Reddy 1979, McDougall 1981, Young & Rosner 2000, Ross & Lavery 2009, and the coffee table experiment in ZS13). The layer thickness increases by a process of merging of neighboring layers. Two mechanisms are observed: vanishing contrast, and drift of interface position. This is illustrated schematically in Fig. 6.

Which of these two dominates, why, and at what rate the merging takes place appears to be not well understood. From the model presented here, however, an estimate for the merging rate which is independent of this uncertainty can be derived, as well as the resulting layer thickness as a function of time.

Both merging mechanisms involve the redistribution of solute between neighboring layers. For example, if an interface between two layers transmits solute somewhat faster than average, the contrast between them decreases. The time scale on which such changes take place is limited by the effective solute diffusivity, hence is of the order

 τ≈d2/κSeff. (68)

The rate of merging is then

 dlnd/dt≈1/τ≈κSeff/d2, (69)

so that

 d≈(2κSefft)1/2, (70)

if is constant in time. The layer thickness is thus predicted to increase as the square root of time. This agrees qualitatively with the dependence inferred empirically in some laboratory experiments. In most of these, however, the layering is induced by lateral heating of the solute gradient rather than vertical heating, hence is not directly comparable with the astrophysical case. Wirtz & Reddy (1979) for example find an initial square-root-like dependence, which saturates when the thickness reaches the separation between the lateral boundaries.

In terms of a merging time scale , (70) becomes

 tm≈d2/2κSeff. (71)

For a quantitative check, compare this to the measurements from lake Kivu, where the observed layer thicknesses are 0.3-0.5 m, the relative interface thickness of order 0.4 (Schmid et al. 2010). With the diffusivity of CO in water, ms, the effective solute diffusivity is

 κSeff=NuSκS≈κS/δ≈710−9 m2/s. (72)

For a layer of thickness 0.4 m this predicts a merging time of order 8 months. This agrees with the observed time scale for changes in the layering, of the order of several months. This comparison has to be taken with a grain of salt, of course, since the present analysis is not strictly valid for the Prandtl number of water.

## 8 Summary

The theory presented expands on the previous analyses in S92 and ZS10. It improves on these by including the effect of a stagnant zone of finite thickness. That is, the region over which heat and solute are transported by diffusion is not limited to just the boundary layers of the overturning interior of a double-diffusive step, but can in principle be an arbitrary fraction of the layer thickness. The need for this extension arose from observations of geophysical examples of thermohaline layering and results from numerical experiments.

A simple 2-zone model consisting of a stagnant and an overturning zone, and using an experimental fitting formula for convective heat transport in the overturning zone produces a clear physical picture for the dependence of the effective transport properties of double diffusive layered convection (Sect. 4). The results predict the existence of a maximum to the density ratio for which such a layered state can exist. It is of the order and approaches this value from below with increasing Rayleigh number (or layer thickness).

Of special interest is the astrophysical case where the heat flux rather than a temperature difference is given. In this case the dependence of effective solute diffusivity on solute stratification and heat flux has the simple form (65). This is the same as before in ZS10 and S92 (except for the effect of radiation pressure on the equation of state not included here). Due to the presence of the stagnant zone the value of the superadiabaticity (63) differs from that in S92, ZS10. For practical conditions in a stellar interior, however, is so small that its exact functional form makes little difference for the temperature stratification. The thickness of the stagnant zone is small relative to the layer thickness, as a consequence of the low value of the Lewis number.

The main conceptual difference with respect to ZS10 and S92 is the existence of a maximum to the density ratio. In the astrophysical case of imposed heat flux it plays only a rather implicit role, however. The differences are more significant when conditions are not in the astrophysically relevant limiting case. In numerical simulations in particular, which are necessarily much closer to marginal conditions for double diffusive layering, the stagnant zone and the maximum density ratio have a substantial effect on the results. A comparison with such simulations is presented in the companion paper ZS13.

Semiconvection is only one of the potential mixing processes in stars. Rotation induced mixing and magnetic processes are likely to be relevant as well, and could actually be more effective.

###### Acknowledgements.
The author thanks F. Zaussinger for extensive discussions and suggestions, and the referee for many corrections and detailed comments which have led to significant improvement.

## References

• (1) Biello, J. A. 2001, Ph.D. Thesis, U. of Chicago
• (2) Langer, N., El Eid, M. F., & Baraffe, I. 1989, A&A, 224, L17
• (3) Langer, N. 1991, A&A, 252, 669
• (4) Langer, N., & Maeder, A. 1995, A&A, 295, 685
• (5) Massaguer, J.M. & Zahn, J.-P. 1980 A&A 87, 315
• (6) McDougall, T. 1981, Progress in Oceanography, 10, 91
• (7) Merryfield, W. J. 1995, ApJ, 444, 318
• (8) Niemela, J., Skrbek, L., Sreenivasan, K.R. & Donnelly, R. 2000, Nature 404, 837.
• (9) Padman, L. & Dillon, T.M., 1987, J. Geophys. Res., 92, 10799
• (10) Proctor, M.R.E. 1981, Journal of Fluid Mechanics, 105, 507
• (11) Rosenblum, E., Garaud, P., Traxler, A., & Stellmach, S. 2011, ApJ, 731, 66
• (12) Ross, T., & Lavery, A. 2009, Experiments in Fluids, 46, 355
• (13) Schmid, M., Busbridge, M., Wüst, A. 2010, Limnol. Oceanogr. 55, 225
• (14) Schmitt, R. W. 1994, Annual Review of Fluid Mechanics, 26, 255
• (15) Spruit, H.C. 1992, A&A 253, 131 (S92)
• (16) Stothers, R. B., & Chin, C.-W. 1994, ApJ, 431, 797
• (17) Tayler, R.J. 1953, Doctoral Dissertation, Cambridge University
• (18) Turner, J. S., & Stommel, H. 1964, Proceedings of the National Academy of Science, 52, 49
• (19) Timmermans, M.-L., Toole, J., Krishfield, R. and Winsor, P. 2008, Journal of Geophysical Research, 113, C00A02
• (20) Turner, J.S. 1980, Buoyancy Effects in Fluids, Cambridge University Press
• (21) Turner, J.S. 1985, Ann. Rev. Fluid Mech. 17, 11
• (22) Young, Y., & Rosner, R. 2000, Phys. Rev. E, 61, 2676
• (23) Weiss, A., 1989, A&A, 209, 135
• (24) Wirtz, R.A. & Reddy, C.S., 1979, Journal of Fluid Mechanics, 91, 451
• (25) Zaussinger, F., & Spruit, H.C. 2010, arXiv:1012.5851 (ZS10)
• (26) Zaussinger, F., & Spruit, H.C. 2013, submitted to A&A (ZS13)
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters