Two-fluid formulation of the cloud-top mixing layer for direct numerical simulation

Theoretical and Computational Fluid Dynamics, Dec 2010

A mixture fraction formulation to perform direct numerical simulations of a disperse and dilute two-phase system consisting of water liquid and vapor in air in local thermodynamic equilibrium using a two-fluid model is derived and discussed. The goal is to understand the assumptions intrinsic to this simplified but commonly employed approach for the study of two-layer buoyancy reversing systems like the cloud-top mixing layer. Emphasis is placed on molecular transport phenomena. In particular, a formulation is proposed that recovers the actual nondiffusive liquid-phase continuum as a limiting case of differential diffusion. High-order numerical schemes suitable for direct numerical simulations in the compressible and Boussinesq limits are described, and simulations are presented to validate the incompressible approach. As expected, the Boussinesq approximation provides an accurate and efficient description of the flow on the scales (of the order of meters) that are considered.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

Two-fluid formulation of the cloud-top mixing layer for direct numerical simulation

Juan Pedro Mellado 0 1 2 Bjorn Stevens 0 1 2 Heiko Schmidt 0 1 2 Norbert Peters 0 1 2 0 H. Schmidt Zuse Institute , FU Berlin, Berlin, Germany 1 B. Stevens Department of Atmospheric Sciences , UC Los Angeles, Los Angeles, CA, USA 2 B. Stevens Max Planck Institute for Meteorology , Hamburg, Germany A mixture fraction formulation to perform direct numerical simulations of a disperse and dilute two-phase system consisting of water liquid and vapor in air in local thermodynamic equilibrium using a twofluid model is derived and discussed. The goal is to understand the assumptions intrinsic to this simplified but commonly employed approach for the study of two-layer buoyancy reversing systems like the cloud-top mixing layer. Emphasis is placed on molecular transport phenomena. In particular, a formulation is proposed that recovers the actual nondiffusive liquid-phase continuum as a limiting case of differential diffusion. High-order numerical schemes suitable for direct numerical simulations in the compressible and Boussinesq limits are described, and simulations are presented to validate the incompressible approach. As expected, the Boussinesq approximation provides an accurate and efficient description of the flow on the scales (of the order of meters) that are considered. Phase transition at the cloud boundaries often compounds the difficulty in understanding turbulent entrainment [15, 21]. There are many different aspects of the problem, which can be considered. One of them, the role of buoyancy reversal due to the evaporative cooling that is promoted by the evaporation of the droplets under certain mixing conditions, has been long debated using theory, field and laboratory measurements, and numerical simulations [7, 14, 25, 28, 29, 36, 46, 51, 54, 61]. Although this problem explicitly involves molecular transport processes, in particular, the transport of latent heat, direct numerical simulation (DNS) has not been employed to address questions related to buoyancy reversal as extensively as it has been used in other areas [39]. Continuing and rapid advances in computational capacity is changing this situation and turning DNS into an increasingly attractive tool for this class of problems. This motivates our study, objective of which is to - explore formulations to perform DNS of moist (and potentially saturated) fluids so as to investigate the latent heat effects in the stratocumulus-top turbulent mixing layers that form between the cloud (cool and moist) and the layer on top, which consists of warmer and dry subsiding air. There are several simplified formulations of moist convection in the literature applicable to our questions. Ooyama [40] describes a model for a moist atmosphere based on water phase equilibrium, the thermodynamic state being calculated from transported values of dry air and water bulk mass densities and entropy. This model is extended by Ooyama [41] to account for a precipitating mode of the condensate phase. Satoh [47] rewrites this formulation in terms of the internal energy instead of the entropy. All of these authors neglect molecular processes and introduce some form of closure instead of deriving a formulation that holds at scales where diffusion dominates; in his study, Ooyama employs a low-pass sixth-order filter while Satoh makes use of a fourth-order numerical diffusivity. However, buoyancy reversal due to evaporative cooling is realized as a result of diffusive processes, and thus it is necessary to develop a model that is valid down to the diffusive scales if one wishes to look at the coupling between the fluid dynamics and the evaporative cooling at the cloud boundary. Bannon [5] does essentially just that, as he derives an equation set in which he retains the molecular transport terms in the governing equations using the results from multicomponent gas mixtures. The previous approaches can be further simplified when the cloud-top configuration is investigated. The thermodynamics is then fully reduced to coupling functions between the desired state variable and the mixture fraction , a conserved scalar that represents the amount of matter in a given fluid particle originating from a specific region in the system. This variable is employed when the boundary conditions can be represented as a mixture of two source fluids: in our case, a homogeneous cloud and a cloud-free layer. In addition, the small size over which the transition between cloud and cloud-free air is realized at the top of stratocumulus clouds, of the order of 10 m, and the small velocities, less than 1 m s1 (cf. Caughey et al. [11]), justify the use of an incompressible approach. The small variation of the density, usually less than 5%, further permits the simplifications associated with Boussinesq. Albrecht et al. [1] introduced the mixture fraction into the analysis of cloud-topped mixed layers, and the formulation has been discussed by Bretherton [8] in a very comprehensive manner, with special attention to the thermodynamics. Since then, this methodology has been often employed as a physical model to investigate the role of latent heat effects in the dynamics of cloud interfaces [24,25,28,29,36,4951]. However, a detailed derivation starting from the basic equations governing the conservation of mass, momentum, and energy has not been presented in the literature, and this is one of the two objectives of this study. In particular, the first objective of this article is to demonstrate how to combine the governing equations of gas mixtures [59], approaches from two-phase flows [16,19] and simplifications in thermodynamics [8,40], in a manner similar to Bannon [5] but which yields a final set of equations in terms of the mixture fraction as discussed by Bretherton [8]. In Sect. 2, we clarify and emphasize the required assumptions, not as a criticism to the existing methodology, but as a complement needed for a correct interpretation of the results obtained with this approach. Section 3 concludes the derivation by defining the cloud-top mixing layer as an idealized configuration designed to study the stratocumulus top in terms of that mixture fraction. The second objective is to present a high-order numerical algorithm suitable for solving the system of equations we pose. The need for an accurate description of the small scales has been mentioned above and progress in this direction has been achieved in recent years [4,52]. Section 4 presents a formulation based on compact Pad schemes which is not restricted to periodic boundary conditions along the vertical direction and that has been used previously in the study of compressible free turbulent flows [35,37]. Both compressible and incompressible formulations are explained and numerical simulations using both of them are presented in Sect. 5 and compared one-to-one, the results confirming the expected accuracy of the Boussinesq approximation for the small density variations occurring at the cloud top. 2 Governing equations The species mass fractions of dry air, water vapor, and liquid water will be denoted by qd, qv, and ql, respectively, according to the notation employed in the literature of atmospheric flows [55], and the total density of the mixture will be denoted by . Additional derived quantities often used are qg = qd + qv, mass fraction of the gas phase, and qt = qv + ql, mass fraction of the total water. The quantities qv, ql, and qt will be also referred to as the specific humidities of vapor, liquid, and total water, respectively. Transport equations will be written for the bulk partial densities qd, qv, and ql. The symbols d and v indicate the partial density of the dry air and vapor in the gaseous mixture, respectively, with g = d + v being the density of the gas. Note that, on the gaseous phase, The relations gd = qd and gv = qv will also be used extensively. The common approximation qg g is justified because the relative error is equal to the ratio between the liquid volume fraction and the gas volume fraction, qlg/(qgl) = l/g, of order 106 considering ql = O(103) and g/l = O(103), and its use in the equations derived below will be explicitly indicated when needed. A detailed discussion of these variables and of the corresponding transport equations describing the evolution of the disperse and dilute two-phase flow is presented in the Appendix, introducing general assumptions needed to eventually reduce them to a simplified set of equations appropriate for the mixture fraction formulation. This final set of equations corresponds to that of multicomponent gas mixtures [38,59] with particular terms in the flux vectors depending on the relative velocity of the cloud droplets with respect to the gas, whose mean value is represented by the mean drift velocity VD. Although the resulting transport equations are similar to those presented by Bannon [5], it is of interest to include them in this article for several reasons: (a) to enable a detailed introduction of notation and completeness of the article, (b) to explicitly estimate the order of magnitude of unclosed terms and justify simplifications in the limit of small droplets needed for the mixture fraction, and (c) to derive them from a different perspective using the concept of local volume averages. In this section, those basic assumptions are particularized to the context of stratocumulus tops and, for this purpose, several estimates of the turbulent scales are necessary; if a reference value of the mean turbulent dissipation rate = 103 m2 s3 is considered [11,30], then the Kolmogorov scales are t = (/)1/2 0.13 s, v = ()1/4 11 mm s1, and l = (3/)1/4 1.4 mm, using a kinematic viscosity of air = 1.6 105 m2 s1. Sections 2.1, 2.2 and 2.3 discuss the three main assumptions implicit in the mixture fraction formulation, Sect. 2.4 summarizes the resulting set of governing equations, and Sect. 2.5 explains the treatment of the discontinuities caused by the thermodynamic equilibrium. 2.1 Approximation 1: the two-fluid formulation The two-fluid approximation models the disperse liquid phase as a continuum, in addition to the vapor and the dry air, and solves the corresponding partial differential equations [13,17,19,22]. This formulation can be understood based on either a local volume average applied to the original governing equations of both phases [13,16] or based on an ensemble average [62]. The former approach is particularized in the Appendix to the cloud system using an averaging volume of the order of l3, and it is shown there that the chain of relations The volume fraction of gas, the ratio of the gas volume to the total volume of the fluid particle, is g = qg/g. Equivalently, the volume fraction of liquid is given by l = ql/l, which satisfies l = 1 g and yields must be satisfied for the two-fluid approach to be valid, where d is the droplet diameter. The first inequality of the sequence can be expressed in terms of the droplet number density nl = l/( d3/6) as nll3 1, and it simply states that there is a large number of droplets in a volume of the order of the Kolmogorov scale for the continuum approach to be reasonable. The second inequality implies dilute conditions l 1, meaning that droplets do not interact among themselves. Altogether, this approximation states that the droplets are so small that they can be considered as a small perturbation of the background flow, smooth within a Kolmogorov length scale, and that they are abundant enough for a continuum description. For the cloud-top system, the typical diameter of the cloud droplet is d = 10 m, indeed smaller than the estimated Kolmogorov scale, and l O(106), which also satisfies the dilution constraint. On the other hand, the mean number density is then of the order of 1 mm3, and the condition nll3 1 is not satisfied. This is, therefore, one intrinsic hypothesis of the two-fluid formulation that is not correct, which clearly indicates that a continuum model of the disperse liquid phase at the cloud boundary, like that based on the mixture fraction, may lead to departures between the calculated and the actual evolutions of the physical system. Many more assumptions are to be made before obtaining the simplified formulation. It is shown in the Appendix that the resulting transport equations have unclosed terms, and a major quantity that is needed is the mean drift velocity VD representing the mean motion of the liquid phase with respect to the gas phase, defined by Eq. 56. A transport equation for this velocity can be derived from the equation describing the motion of a single droplet in the limit of low Reynolds numbers [34] and then taking some appropriate mean value [17, 19]. This detailed approach would allow us, for instance, the investigation of preferential accumulation effects, though numerical problems still exist within the two-fluid formulation due to the strong local gradients of l that may be formed [22]. Here, however, these unclosed terms are further simplified, as discussed later. 2.2 Approximation 2: thermodynamic equilibrium Mechanical equilibrium assumes that the Stokes number is small enough [13, 48]. A Reynolds number based on a cloud droplet diameter of 10 m and a slip velocity of the order of v is vd / = 0.007, and the equations for the evolution of the droplet mass, velocity, and temperature are developed in the context of low Reynolds numbers [33, 34]. These equations show that relative changes of the droplet velocity of order 1 occur on a time scale tp = ld 2/(18g) 0.3 ms, where g = 1.9 kg m1s1 is the dynamic viscosity of the gas mixture. Then, the Stokes number tp/t is less than 0.01, which suggests that the cloud droplets reach equilibrium with the local environment faster than this local environment evolves, with a terminal velocity gtp with respect to the gas of the order of 3 mm s1 using a gravity acceleration of 9.8 m s2. This limit of small droplet inertia simplifies considerably the description of the evolution of the mixture because the transport equation for the mean drift velocity VD is completely substituted by the algebraic explicit relation VD = gtp, which equates this mean drift velocity with the terminal velocity of a single droplet. The settling parameter V D/v for the system is then about 0.3. We show in the Appendix how different unclosed terms of the transport equations are proportional to some power of this ratio, and therefore the simplifications in the limit V D 0 seem to be justified. In particular, the leading order errors occur in the mass and heat flux vectors, Eqs. 67 and 81, respectively, and these are of the order of (ql/qv)(V D/v) 0.03. However, although numerical studies [53, 58] and experiments [2] show that the strongest interaction between particles and turbulence occurs at Stokes numbers of order 1, preferential accumulation at much lower values, about 102, are currently being measured in clouds [30], which indicates that inertial effects may still persist at those small values of the Stokes number. The situation is more complicated when phase equilibrium is considered. Phase transition introduces its own time scale (different from tp) into the problem, and this time scale is of the order of a few seconds in typical clouds, i.e., large compared with t [48]. As a consequence, a source term describing phase nonequilibrium appears in the transport equation of certain variables (see Appendix). Moreover, this nonequilibrium may lead to inhomogeneity of the thermodynamic fields inside of the local averaging volume, and then the usual closures for the molecular transport terms within the gas phase and the ideal gas thermodynamic relations applied to the local average quantities are questionable. In sum, complete thermodynamic equilibrium is not achieved in real conditions due to these condensation/vaporization effects even though the Stokes number is small, and this constitutes a second hypothesis of the mixture fraction formulation that is not correct. 2.3 Approximation 3: liquid-phase diffusivity So far, the problem can be understood consistently in the limit d / l 0 because then (a) droplet relaxation times are relatively small and equilibrium holds, and (b) the terminal velocity is relatively small and the settling parameter V D/v vanishes. However, the limit d / l 0 within the context of equilibrium formally implies a droplet size d below which thermal motion overwhelms the gravitational field and the mean drift velocity VD should then be related to the corresponding diffusion flux of the liquid, and not to the terminal velocity. At the same time, it will be shown in Sect. 3 that the mixture fraction formulation requires certain specific conditions on the diffusion terms that cannot be satisfied by the governing equations obtained in the Appendix if we simply substitute VD = 0 in them. This issue is investigated next. The droplet velocity scale vB associated with the thermal motion can be estimated from the relation l(/6)d 3v2B = 3k B T , where k B = 1.381 1023 J K1 is the Boltzmann constant [44, 45]. If we equate this expression with the terminal velocity gtp = gld 2/(18g), the cross-over diameter in air at 300 K is 4 m, with a velocity about 0.5 mm s1: smaller droplets have thermal fluctuations larger than the terminal velocity, and larger droplets have them smaller. Hence, for small enough droplets, the velocity scale vB should be used in the estimates of the Appendix related to second- and third-order moments of the velocity fields inside of the local averaging volume, Eqs. 72 and 79, but those terms involve the liquid volume fraction l as a multiplying factor, and it seems reasonable to still neglect them compared to the other unclosed terms related to VD in view of the small values l = O(106) corresponding to the cloud. (In other cases, they are normally grouped together with the molecular transport terms and then closed defining new transport coefficients [59].) With respect to the terms involving this mean drift velocity VD, it is shown in the Appendix that the leading order corrections appear in the mass and heat flux vectors, Eqs. 67 and 81, respectively, and they describe the transport of latent heat. Ramshaw [44] provides the following expression for the mean relative velocity of the liquid continuum with respect to the total mixture due to Brownian motion where the diffusion coefficients l,T and l,p can be related to l = (1 l)B, and B = kB T /(3 gd) is the Einstein diffusion coefficient [18]. For temperatures of 300 K and droplets of diameter 10 m, this diffusion coefficient is of the order of 1012 m2 s1. This value is very small. First, for the Kolmogorov length scales of the order of 1 mm, which develop in our problem, the corresponding diffusion velocity of the liquid phase is considerably smaller than the settling velocity due to gravity, and it explains that it is always neglected. Second, it is much smaller than the other typical diffusivities entering the equations and consequently, if the terminal velocity is neglected, very thin layers of the field ql can develop, which demands computational resources not yet available. These results highlight how the implicit assumption in the mixture fraction formulation (see Sect. 3) of a liquid-phase diffusivity comparable to that of the gas constituents of the system, though asymptotically consistent (with the limit d/ l 0 and Eq. 3), does not hold in reality at the stratocumulus top, and it constitutes a third simplification (or approximation) to be added to those of liquid-phase continuum and equilibrium. To the extent that one is prepared to accept the aforementioned limitations, the derivation toward a mixture fraction formulation can be continued by adding a diffusion model to the governing equations (as described in the Appendix) such that the liquid phase diffuses with the gas according to Ficks law and with a diffusivity l, This model corresponds to Eq. 4 when only the term proportional to the mass fraction gradient is retained and can be also interpreted as an ad hoc definition of VD in Eq. 56 as VD = (l/qg) ln ql. (The tilde used in the Appendix to represent local average quantities has been dropped in the main text for notational convenience.) The set of relative velocities resulting from this model is completed by VD,d = vd v = v ln qd (l v) ln qg, VD,v = vv v = v ln qv (l v) ln qg. The mass and heat flux vectors defined in Eqs. 67 and 81 can then be expressed using the new mean drift velocity as jt = qvVD,v + qlVD,l = vqt (l v)(qd /qg)ql. = T v(hdqd + hvqv + hlql) (l v)(hl hg)ql, or in terms of the total enthalpy of the mixture Differential diffusion effects appear into the problem when l = v such that the limit V D 0 is recovered asymptotically as l 0. Table 1 Thermodynamic data Rd Rv cp,d cp,v cl el2v73 0.2870 0.4615 1.0070 1.8700 4.2176 2.5016 Gas constants and (constant) heat capacities are in units kJ kg1 K1, and latent heat of vaporization is given at 273.15 K in units MJ kg1 2.4 Nondimensional equations The governing equations presented in the Appendix, with the diffusion terms as formulated in the previous section, are now summarized and nondimensionalized with a reference density 0, velocity U0, length L0, temperature T0, and reference values of the heat capacity c p,0 and dynamic viscosity 0. These equations describe the conservation of mass, Eqs. 65, 66, 68, conservation of momentum, Eq. 73, and conservation of energy, Eq. 79. (The tilde used in the Appendix to represent average quantities has been dropped in the main text for notational convenience.) The pressure p is nondimensionalized with 0U02, the specific internal energy e and enthalpy h with cp,0T0, and ug is a unitary vector in the direction of the volumetric force. The transport coefficients v, l, and /cp are written in terms of and the Prandtl and Schmidt numbers by means of the definitions /cp = / Pr, v = /Sc, and l = /Scl, and Pr, Sc, and Scl are taken as constant. Lewis numbers can be defined as Le = Sc/ Pr and Lel = Scl/ Pr . The transport equations are /t (e) + (ev) = (0 1)M 2 pv + (0 1) with the molecular transport terms and the equilibrium condition ql, Scl Sc (hl hg)ql, Nondimensional thermodynamic relations are written as 0 M 2 p/ = T [1 + qt(Rv/Rd 1) ql Rv/Rd], e = [(qt ql)cv,v + (1 qt)cv,d + qlcl]T ql Q, where the reference ratio of specific heats is defined by 0 1 = 1/(c p,0/Rd 1); Rd/c p,0 and Rv/c p,0 are the nondimensional dry air and vapor gas constants, respectively. Table 1 shows the dimensional thermodynamic constants used in this study. The reference nondimensional numbers appearing in these equations, apart from the ratios already defined before, are Ri = M = Fig. 1 Species mass fractions ql, qv as a function of water content qt for a fixed state (, e). Exact equilibrium (a) and modified equilibrium (b) as given by Eq. 25 with s = qs/16. The reference state corresponds to saturation conditions at p = 940 hPa and T = 10.6C (qs = 8.50 g kg1) A density difference has been assumed to enter into the problem. Re is the Reynolds number, Ri is the Richardson number, M is the Mach number, and Q is a nondimensional latent heat parameter. Other combinations can be used, and additional nondimensional parameters can appear through the boundary or initial conditions, like those related to environmental lapse rates, inversion thicknesses, or the characterization of the initial perturbation. If there is no velocity scale externally imposed on the system, then we can use ( / 0)g L0 in place of U02 in the nondimensionalization and write the problem in terms of Gr = H/L0 = e0 lv , Q = c p,0T0 where Gr = Re2 is the Grashof number, 0 M 2 = (L0/H )( / 0), and Ri = 1; there is one degree of freedom less in the parameter space. The scale-height is defined by H = Rd T0/g; quite often, the Rayleigh number defined as Ra = Gr Pr is used instead of Gr . 2.5 Equilibrium composition We follow here previously described approaches [40,47] to obtain ql,eq . For a given thermodynamic state ( , e, qt), a temperature T can be computed assuming all of the liquid in vapor phase, and the requirement that qt ps(T )/( RvT ) can then be evaluated; if this inequality is satisfied, then the equilibrium state has been already obtained, and if not, then a nonlinear equation needs to be solved by substituting the conditions qv = ps(T )/( RvT ) and ql = qt qv into the energy equation. The saturation pressure ps(T ) is calculated with the polynomial fit provided by Flatau [23]. A ninth-order polynomial is solved with a NewtonRaphson method, and convergence is achieved in three iterations below 107 relative error. The consistency error between the corresponding latent heat obtained from the Clausius Clapeyron equation and the linear relation due to the constant heat capacities is less than 5% within the interval of relevance. Figure 1 shows the variation of the composition with the specific humidity of total water, with constant density and energy. For low qt (unsaturated conditions), as dry air is substituted by water and qt is increased, the temperature decreases because the water vapor has larger heat capacity, and the energy e is fixed. This decrease in temperature implies a decrease in the function ps(T )/( RvT ). At the same time, the vapor content qv increases because qv = qt. This behavior continues until both curves cross with each other, and at this point, the system is saturated, obtaining the value qs( , e). If more water is added to the system, then part of this water is present in the form of liquid instead of vapor and the contribution of the latent heat to the energy implies an increase in the temperature of the mixture, which in turn implies a growth of the vapor content qv = ps(T )/( RvT ). Figure 1a shows that there is always a discontinuity in the derivative ql/qt at constant and e, and the same is true for the function qv, which translates into a discontinuity with respect to the space coordinates of thermodynamic quantities needed in the governing equations. This problem was already identified by Ooyama [40], although here, in our study, it is not only the pressure gradient but also the diffusion terms (where second-order derivatives of the temperature appear) that introduce Gibbs phenomena, whose effect will manifest later in the simulations in the form of a noisy background dilatation field. A smoothing function similar to that used in DNS of the BurkeSchumann problem of infinitely fast chemistry [42] has been considered and is now discussed. Consider the saturation specific humidity qs( , e) where Ts( , e) is the temperature at saturation. This temperature is obtained by solving the nonlinear equation that results after substituting qv = qs and ql = 0 in the energy equation. Second, consider the function ql(qt; , e) and decompose it into a piecewise linear part and a remaining contribution ql in the following manner 0, ql qt qs (qt qs) + qt qs where fixed values of and e are assumed in qs( , e). Then, the derivative of the piecewise linear part can be approximated by a hyperbolic tangent, and the integration provides the following approximation for ql(qt; , e) ql qt qs qt qs The temperature is then recalculated so that the values of density and energy are maintained. This algorithm can be thought of as a replacement of the exact caloric equation of state by a modified one. Figure 1b shows that the resulting modified profile has a continuous derivative at qt = qs( , e), the discontinuity being now only in the second derivative. The smoothing factor s gives the interval in qt space over which the function is smoothed and should be defined as a fraction of min{qs, qt,max qs}. Note that this approach becomes computationally more expensive; the value qs( , e) has to be computed for the unsaturated region of the flow as well, and the calculation of the slope at the saturation point 3 The cloud-top mixing layer ql qt qs = 1 e ed,0 + qs[ed,0 ev,0 + Rv( +1)Ts] 1 = ev,0 el,0 (cl cv,v)Ts 1 The equations obtained in the previous section are now particularized for the cloud-top mixing layer. This system is defined by two infinite (unbounded) horizontal layers, one layer of warm and unsaturated fluid on top of a second layer of cooler and saturated (condensate laden) fluid, the gravity force acting downward, with or without shear. The system is statistically homogeneous along horizontal planes and temporally evolving. This idealized configuration models the top of stratocumulus clouds on a vertical length scale of a few meters with the objective of studying latent heat effects on the local dynamics. The definition of the initial conditions is split into the setup of mean velocity and thermodynamic profiles, and the addition of perturbation fields. The thermodynamic mean profiles must satisfy the hydrostatic equilibrium equation L0 d z + [1 qt + (qt ql)(Rv/Rd)]T p(zc) = pc, where a reference pressure level pc has to be given at the height zc and H = Rd T0/g is the scale height, as introduced before. Two more thermodynamic relations need to be provided, and two options are here considered. The first option is to impose {T (z), qt(z)}, i.e., given profiles of temperature and water content. With this data known, Eq. 28 can be solved analytically if the moist air is unsaturated because then ql 0. On the other hand, if the moist air is saturated, then the composition is not known because qv = qs also depends on p and not only on T , but Eq. 28 leads to L0 d z + (1 qt)T (1 qt)T p(zc) = pc, the solution of which can also be reduced to a quadrature. The second option considers {h(z), qt(z)} given, i.e., the enthalpy is specified instead of the temperature. This case appears naturally when studying the buoyancy reversal instability, as later discussed. The initial hydrostatic equilibrium is then solved iteratively. In any of the two cases, the imposed thermodynamic mean profiles follow an error function f (z) between two levels f0 and f1, f (z) = f0 + f1 f0 z zc Both variables are then conserved scalars with the same diffusivity and obey the same evolution equation. Therefore, with appropriate boundary and initial conditions, the calculation of both qt and h is reduced to the If a mean gradient needs to be set in each of the layers, 0 and 1, respectively, then the following function is added g(z) = x 1 0 2 z zc exp z zc which also satisfies the one-dimensional diffusion equation with diffusivity if d(2)/dt = . Now we turn to discuss the governing equations in Sect. 2.4 within the two-layer system defined above. First, the Mach number M is a handicap of this system because the characteristic velocities of the problem are less than 1 m s1, which implies that M < 1/300 and the time resolution of the acoustic waves dictates a very small time step in an explicit numerical scheme. Some physical problems allow the artificial increase of the reference Mach number up to approximately 0.3 (the exact value depending on the definition and the particular problem) because the results of interest are insensitive to the Mach number below this threshold. However, the cloud-top mixing layer does not allow this approach. In order to understand why, consider an expression for the Mach number in terms of the scale height H and other nondimensional parameters of the problem as The Richardson number Ri and the strength of the inversion / 0 are the dominant terms and need to be reproduced in the simulations; hence, this equation shows that an increase of the Mach number is equivalent to an increase of the length L0 compared to the scale height H , this latter being of the order of 8000 m in the atmosphere. However, the cloud layer (condensate phase) usually varies in a much smaller distance (of the order of 100 m), and therefore L0 cannot be significantly augmented. There is a need for a low Mach number formulation. Second, neglecting differential diffusion effects (Scl / Pr = Sc/ Pr = 1) and accepting the previously stated limitations, the transport equations for the total specific humidity and the enthalpy in the low Mach number limit become e b 1 /b 0.4 Fig. 2 Sketch (a) showing the two-layer structure of the cloud-top mixing layer in terms of the conserved variables and the three-layer structure of the density due to buoyancy reversal. Nondimensional buoyancy mixing function (b) for data in Table 2: dashed, exact thermodynamic equilibrium; solid, approximation Eq. 37 knowledge of one normalized conserved scalar, the mixture fraction , which evolves according to the advectiondiffusion equation as above. In the cloud-top mixing layer without lapse rates, denoting with subscript 0 and 1 the lower and upper layer, respectively, we can choose qt qt,0 h h0 = h1 h0 The mixture fraction so defined indicates the relative amount of matter of the fluid particle that originates from the upper layer. According to the notation here adopted, = 0 in the lower layer and = 1 in the upper layer, as depicted in Fig. 2a. Bretherton [8] refers to variables obeying this relation as linearly mixing variables. In combustion theory, linear relations between thermodynamic variables that are conserved, known as coupling or ShvabZeldovich relations, also lead to the introduction of a mixture fraction as described here [59], and this formulation in terms of has allowed a very significant development in the study and modeling of nonpremixed turbulent reacting flows during the last decades [43]. In order to conclude, note that the assumption that h and qt follow a similar initial condition is not as restrictive as it seems. Both quantities (qt qt,0)/(qt,1 qt,0) and (h h0)/(h1 h0) obey the same advectiondiffusion equation in an infinite two-layer system with the same boundary conditions. Given different initial conditions, there will simply be a (turbulent) transient in which this difference between the normalized enthalpy and water content evolves toward zero. Hence, the fluid particle is advected and its mixture fraction diffuses with the environment, with an approximately constant thermodynamic pressure. The value of at each point and time determines completely the thermodynamic state of the fluid particle through qt( ) and h( ), in particular, the density field (x, t ) = e( (x, t )), the pressure level entering only as a constant external parameter of the problem. Small fluctuations of the density with respect to a mean value, normally related to the fluctuation of two thermodynamic variables (like qt and h here) for a given thermodynamic pressure, can now be written in terms of the conserved scalar [8, 55]. This linear thermodynamic analysis leads to an expression of the buoyancy b, in terms only of as a function be( ) which is piecewise-linear to a very good approximation, as illustrated in Fig. 2b. Partial derivatives in saturated and unsaturated regimes can be calculated with more or less detail, but at the end it comes down to providing three parameters describing the function. A usual choice is [51]: b1/g = (0 1)/0 = / 0, which quantifies the strength of the stable inversion; s, which is the mixture fraction at saturation conditions; and the buoyancy reversal parameter [51] D = Table 2 Buoyancy reversal parameters s, D (defined by Eq. 36) and b1 for different thermodynamic conditions (T0, qt,0) in the lower layer qt,0(g kg1) b1(102g) Pressure level 940 hPa and upper layer at 19.1C and qt,1 = 1.50 g kg1 which characterizes the buoyancy value at saturation, bs = be(s). Figure 2b shows the function be( )/b1 for the cases described in Table 2 as obtained from the equilibrium calculation described in Sect. 2.5. The cases considered in that table are derived from the reference case A1, this last one corresponding to field experimental data of nocturnal marine stratocumulus from the DYCOMS-II campaign [56]. The thermodynamic state of the upper layer is kept fixed at a temperature T1 = 19.1C and with a total-water specific humidity qt,1 = 1.5 g kg1, and the lower state is then varied. The case A0 does not retain evaporation, and cases A2 and A3 augment the buoyancy reversal by means of a higher water content in the lower layer, keeping the strength of the inversion b1 constant. The discontinuity in the derivative observed in those curves at s (saturation conditions) corresponds to that in Fig. 1a and also needs to be smoothed if the high-order schemes normally employed in DNS, and later described, are to be used. Therefore, the buoyancy function is approximated as which corresponds to the profile of the derivative dbe/d following a hyperbolic tangent between two different levels and centered at s (see also Fig. 2). Since is simply a normalized qt, this approximation is the equivalent to Eq. 25, but written here for the density instead of ql and retaining only the piecewise linear contribution. Mixture fractions smaller than the cross-over value c = (s + D)/(1 + D) are negatively buoyant and represent the phenomenon of buoyancy reversal (in our case, due to evaporative cooling). Numerical studies performed to calculate the influence of the smoothing parameter s on growth rates of the thickness of the mixing layer have been presented elsewhere [36] and show that s = s/16 is sufficiently small to regularize the phase transition without noticeably influencing the evolution of the flow for the cases of interest presented in Table 2. In the Boussinesq limit ( / 0 1 and D / 0 1) with constant transport coefficients the nondimensional governing equations are then v = 0 The set of nondimensional parameters is {Re, Pr, Ri, D, s}, apart from those introduced by the boundary and initial conditions, and it reduces to {Gr, Pr, D, s} if there is no velocity scale imposed on the system, as discussed in Sect. 2.4. This set of equations is the final result we were pursuing hitherto. It is similar to that proposed by Bretherton [8], where adiabatic motion was first considered and diffusion was later introduced as a common postulate. The contribution from this article has been the discussion of these diffusion terms and the clarification of the assumptions needed to arrive to this simplified mixture fraction formulation. As a reminder, the major hypotheses are (a) liquid-phase continuum, (b) local thermodynamic equilibrium, and (c) liquid-phase diffusivity equal to the vapor and thermal diffusivities. Estimates in Sect. 2 show that these assumptions are not generally met at the top of stratocumulus clouds; nevertheless, this physical model has proven useful in the study of some aspects of latent heat effects at the boundary between cloud and cloud free-air during the past years, as indicated in the introduction. 4 Numerical algorithms A compressible algorithm used previously for DNS of free turbulent flows [37] has been complemented with the required thermodynamic equilibrium formulae to solve the system of equations of Sect. 2.4 applied to the cloud-top mixing layer described in the previous section. It is based on finite differences, and the spatial firstand second-order derivatives are computed using a sixth-order compact Pad scheme. It is one-sided at the nonperiodic boundaries (top and bottom), having globally fourth-order accuracy [10]. The advancement in time is performed with a low-storage five-stage, fourth-order RungeKutta scheme [9]. The convective term is written in the skew-symmetric form [20] to reduce aliasing errors [27], and the diffusion terms are treated in explicit form. A rectangular domain is considered assuming periodicity in the horizontal planes parallel to x O y and gravity acting downward, i.e., ug = k; the boundary conditions in the nonperiodic directions are implemented using a characteristic nonreflective formulation [32,57]. Compact Pad schemes are thoroughly discussed by Lele [31]. The resulting finite-difference approximations z p and zz p to the first- and second-order derivatives, respectively, of a scalar field p along the O z direction are calculated by solving the linear systems A1z p = B1 p and A2zz p = B2 p. For the sixth-order schemes used here, the matrices Ai are tridiagonal and the matrices Bi are pentadiagonal, and the resolving efficiency measured by the 0.1% error in the corresponding transfer function occurs at about six points per wavelength. An incompressible algorithm for solving the set of equations presented in Sect. 3 has been derived from this compressible formulation following Wilson et al. [60]. In this case, the boundary conditions imposed at the top and the bottom are zero normal velocity (no-penetration) and zero normal derivative of the horizontal velocities (free-slip) and the scalar field (adiabatic). The Neumann boundary conditions for the Poisson equation for the pressure at the top and the bottom are then [26] where the boundary conditions on the velocity have been already applied, and w is the vertical velocity. In addition, one reference value of p has to be given at one point (this value is irrelevant since only p is needed; recall that the effect of the thermodynamic reference pressure is already included in the curve be( )). The discrete Poisson equation for the pressure is written using Fourier decomposition inside the horizontal planes, which leads to the following set of difference equations along the vertical direction O z where pi j is the Nz -dimensional vector that contains the horizontal Fourier modes i = 0, . . . , Nx/2 and j = Ny/2 + 1, . . . , Ny/2 at each vertical position zk , k = 1, . . . , Nz , and In these expressions, Nx , Ny , and Nz represent the number of grid points in the directions O x , O y, and O z, respectively, and x , y, and z indicate the corresponding grid steps. The transfer function f1() of the first-order finite-difference operator needed in the equation above is and f 2() is plotted in Fig. 3a as a function of the scaled wavenumber . The boundary value problem is 1 completely defined by specifying the Neumann boundary conditions {z pi j |1, z pi j |Nz } at the bottom k = 1 and the top k = Nz , values obtained by Fourier transforming Eq. 39. For the case = 0, one of the Neumann boundary conditions has to be substituted by a Dirichlet one, and p00|1 = 0 is used without loss of generality. The three modes (Nx/2, 0), (0, Ny/2), and (Nx/2, Ny/2) are degenerate in the sense that i j = 0, and the problem with Neumann boundary conditions, in general, does not have a solution. This is a consequence of the finite-difference approximation to the first-order derivative, which sets an eigenvalue equal to zero at = , as shown in Fig. 3a, and implies that it is not possible to set the dilatation exactly equal to zero. Nonetheless, the order of magnitude of this error is proportional to the energy of these three modes, which should be small if the simulation is reasonably well resolved, or even zero if some dealiasing technique is used, because they represent part of the 2 x and 2 y plane waves. Fig. 3 Transfer function (a) of several finite-difference operators approximating the second-order derivative: solid, exact; dashed, once sixth-order f2() from Eq. 45; dot-dashed, twice sixth-order f12() from Eq. 42; dot-dash-dashed, twice eighth-order f12() from Eq. 46. Absolute error (b) in the approximation of f2 by f12 Solving the set of discrete equations Eq. 40 poses some difficulties, since each one is a linear system with a full matrix of size Nz Nz . The problem can be simplified by introducing an approximation to the operator z z which leads to a system easier to solve. For instance, introducing the finite-difference operator zz used to calculate the second-order derivative, we have with R p = z z p zz p = ( A11 B1 A1 B1 A1 B2) p and R being a full matrix. Neglecting this term 1 2 [12, 60] leads to a pentadiagonal linear system to obtain an approximation to the solution of Eq. 40 with the same order of accuracy. The error R p introduced by this step in solving the Poisson equation due to the different truncation error between z z and zz is easily analyzed in the case of periodic boundary conditions with help of the corresponding transfer functions and it is presented in Fig. 3b. The transfer function of the finite-difference approximation to the second-order derivative is and the truncation error introduced when solving the Poisson equation is smaller, as shown in Fig. 3. The new scheme is only slightly more expensive, since the matrix B1 is now heptadiagonal instead of pentadiagonal and the maximum CFL number diminishes only by 7%. Further improvement could be achieved developing spectral-like schemes with appropriate constraints [31], but that approach is not pursued here any further. One reason is that, in addition to this truncation error, if nonuniform grids are of interest, the previous formulation does not allow a simple treatment based on the Jacobian of the stretching function because then the problem results again in a full linear system. An improved method to solve Eq. 40 and to cope with both disadvantages can be derived from the solution to the continuous boundary value problem associated with it, a solution that can be written explicitly in terms of quadratures since it is a second-order linear ordinary differential equation with constant coefficients. In particular, the solution can be expressed as b + ezN a e(zN zk ) + a + ezN b ezk + h(1)|k h(2)|k , where is the principal square root of 2 = 0, and it is assumed that z1 = 0, without loss of generality (z N is then equal to the domain size in the O z direction). For notational convenience, the subscript i j has been dropped and the symbol N is chosen to represent Nz . The integration constants {a, b} are to be determined from the specified boundary conditions, and p|k represents the k component of the N -dimensional vector p. The vectors h(1) and h(2) contain the particular solution of the equation and are obtained by solving the two first-order linear problems which then implies that the integration constants of Eq. 47 are given by The derivative along O z of the solution is obtained directly as 1 z p|k = 2 b + ezN a e(zN zk ) a + ezN b ezk + h(1)|k + h(2)|k . For the case = 0, the pressure is obtained solving sequentially, first, for z p with z p|N given, and then for p with p|1 given. The particular form chosen for Eq. 48 instead of, for instance, an explicit quadrature, is motivated by the rounding errors in the floating-point arithmetic calculating the exponential terms at the boundaries for large values of ; the expressions presented here led to the best general behavior. For a similar reason, those specific boundary conditions for h(1) and h(2) where chosen to avoid the formation of a boundary layer in these functions for arbitrary large values of . Note that Eq. 41 implies that O( z) z O(1) if the aspect ratios of the grid are of order 1 and the problem is nondimensionalized with the size of the computational domain. The two systems Eq. 48 are now solved using the same compact schemes as employed before to calculate the first-order derivatives. Let us consider for instance h(1). The linear system to be solved can be written compactly in matrix form as where A indicates the matrix A1 premultiplied by the corresponding Jacobian of the stretched grid. In par1 ticular, the sixth-order schemes used here lead to 0 7/9 1/36 . . . A1 = 1 2 1/4 1 1/4 1/3 1 1/3 . . . where a constant grid step z has been assumed for simplicity. The last row corresponds to the boundary condition h(1)|N = g|N /. The other boundary points are implemented using a forth-order compact scheme at k = 2 and k = N 1 and a third-order biased one at k = 1. The pentadiagonal system can be solved efficiently using an LU decomposition similar to the Thomas algorithm for the tridiagonal case. The operation count for the factorization step is O(10N ), and the forward/backward substitution step requires O(9N ). The eigenvalues k of the (N 1)2 block matrix in B1 A1 are of interest because they characterize the linear system to be solved. They have been calculated using the GNU scientific library and the spectra are shown in Fig. 4a when = 0 for two different grid sizes N . All of them are nonzero, and therefore the matrix is regular and can be inverted. (Note that the boundary condition eliminates the two zero eigenvalues corresponding to = 0 and = shown in Fig. 3a for the periodic case and moves the eigenvalues out of the imaginary axis.) One of them is real 1.6627, independent of the size of the system N (within the accuracy of the numerical calculation). The rest are complex conjugate pairs in the negative half of the complex plane, and, when the real part is scaled with N , they follow closely the curve shown in that figure, where the final point in the real axis is precisely 1.6627N ; i.e., they approach the imaginary axis at a rate N 1. A related quantity that is of interest is the spectral condition number of the matrix, which measures the sensitivity of the solution to perturbations of the independent term. This is shown in Fig. 4b for different grid sizes and two values of z. In this respect, it is noted that z in our problem varies between small values of order 1/Nx or 1/Ny and values of order 1, and therefore only the limiting cases = 0 and z = 1 are included in that figure. The worst case scenario occurs for = 0, corresponding to the matrix B1, where the condition number can be well fitted to the linear behavior 2.80 N . As increases, the condition number decreases. This result implies a loss of about three digits of accuracy for grid sizes of the order of 1024 when solving the Poisson equation due to errors in the forcing term. (For comparison, the condition number of the circulant matrix A1 employed to approximate the first-order derivative with periodic boundary conditions is equal to 5 and independent of N .) In brief, this second formulation based on Eq. 47 improves the first approach Eq. 44 by reducing the dilatation error and, more importantly, allowing stretching in the vertical directionsmaller grids reduce both memory requirements and the overall computational load. The main disadvantage of this new formulation is that two pentadiagonal systems need to be solved, instead of only one. However, for simulations involving large grid sizes and using domain decomposition in parallel computations based on the Message Parallel Interface, this overhead is small compared to the communication time required to transpose the data across all the processors to perform the Fourier decomposition. 5 Simulations The main purpose of this section is the validation of the incompressible algorithm by comparison of results with those obtained from the compressible one. We use as a test case the shear-free configuration under buoyancy reversal conditions discussed by Mellado et al. [36]. Physical details about the flow and the buoyancy reversal instability are also provided therein. A case without an imposed mean shear is chosen because then all of the motion is forced by the initial density distribution and diffusion processes, and it should stress the possible differences between the two formulations. The main characteristic of the system is that the two-layer structure in the mixture fraction (i.e., enthalpy and total-water specific humidity) corresponds to a three-layer system in terms of the density, with the central one heavier than the lower one (cf. Fig. 2a). This configuration has one stable mode of period 2 /b1 (or two dispersive waves with phase velocity b1/(4 )) and Fig. 5 Negative buoyancy field for case A3 showing the evolution (left to right) starting from the initial condition and showing a frame every cycle of linear stable mode: top, compressible formulation; bottom, Boussinesq formulation. The box height shown is only the central 2/3 of the domain employed in the simulation one unstable mode with a characteristic time scale D1/2 times larger than the stable one. The system is then linearly unstable and the corresponding buoyancy reversal instability can lead to a turbulent flow in the lower layer. Figure 5 shows qualitatively the evolution of the instability in terms of the negative buoyancy field every 2 /b1 time units (points x for which b(x, t ) < 0 are visualized using a gray color scale, white for zero, and black for minimum)superimposed on a standing interfacial gravity wave a downdraft develops at the lowest point of the initial perturbation. In particular, we focus on case A3 of Table 2 because it corresponds to the strongest buoyancy reversal and leads to a faster development of the flow. A two-layer mean thermodynamic state is set as described in Sect. 3 and then perturbed by displacing sinusoidally the central isosurface of , in the incompressible case, and of qt and h, in the compressible case, from the hydrostatic equilibrium over a wavelength with an amplitude (a/2)/ = 0.1. The thicknesses of the initial profiles are / = 0.05. Mellado et al. [36] shows that a grid 200 600 allows a reference Grashof number Gr = 3b1/2 = 4 108, with a Prandtl number equal to 1; the reference length used in Sect. 2.4 has been chosen equal to the wavelength, L 0 = . In addition, the incompressible mode needs the corresponding values of D and s of Table 2. For the compressible case we are required to specify b1/g = / 0 = 0.0254, as obtained from the initial thermodynamic state and shown in Table 2, and the scale-height ratio H / = 2.8 103, where the reference thermodynamic state is dry air at the temperature and density of the lower layer. For this ratio, a reference length = 3 m has been chosen, small enough compared to the height-scale H of the density variation in the lower layer. (Note that the viscosity inferred from the previous Gr and the value b1 = 0.25 m s2 measured at the stratocumulus top is then larger than the real one in the atmosphere, but that is irrelevant for the validation purpose of these simulations, and it was more important to retain a higher Mach number to be able to perform a compressible simulation.) The b / /12 )110-3 ( 10-4 s )rm v (10-5 Fig. 6 Temporal evolution (a) of the penetration length hb of the downdrafts for case A3 of Table 2: solid, compressible; dashed, Boussinesq. Effect of smoothing parameter s on the compressible simulations (b) rest of the thermodynamic data needed for the compressible case is given in Table 2, case A3, and the heat release parameter Q introduced in Sect. 2.4 is 8.14. The comparison between the compressible formulation and the Boussinesq limit in Fig. 5 shows the high accuracy of the latter approximation. Relative errors due to compressibility, measured by M 2 = (/ H )( / 0)/0 105, are negligible in our case compared to those introduced by the Boussinesq approximation. The assumption of constant density in the inertial terms leads to errors of the order 3% at the inversion and Fig. 5 shows that these errors are indeed small and the evolution of the Boussinesq fluid follows that of the compressible solution. The compressible case also exhibits a slight background gradient in the buoyancy field due to compressibility, which eventually would modify the evolution if the downdraft penetrates downward enough, but it is not relevant for the sizes of the system considered in this study, being of the order of a few meters. The only observed effect of the Boussinesq approximation is a very small difference in the pattern of the recirculated negatively buoyant fluid inside of the mushroom structure in the last frames. In contrast, the savings in computational cost allowed by the Boussinesq approach are very significant: about 600 iterations are needed in the incompressible case compared to 1.9 106 iterations for the compressible case, and the computational time per iteration required by this latter is about 2.3 times larger than that needed by the Boussinesq formulation. The non-linear evolution of the flow is also quantitatively well represented by the Boussinesq approximation. In order to see this, we define the height of a falling finger hb(t ) based on the mean vertical profile of the mixture fraction. We consider that mean profile and scan from the lowest boundary upward until the mean profile reaches a given threshold; the distance from this point to the inversion plane defines the height hb(t ). A threshold value 0.001 is used, and the results are presented in Fig. 6a. This figure shows the increase with time of the downdraft height, and it confirms the accuracy of the Boussinesq approximation for the usual thermodynamic conditions characterizing the stratocumulus top. The maximum relative difference between both curves is about 3%. Last, the effect of the smoothing parameter s in the compressible algorithm discussed in Sect. 2.5 (cf. Fig. 1) is considered. Recall that this parameter needs to be introduced to avoid the discontinuity in the derivatives of the thermodynamic functions at saturation conditions. This issue is investigated by means of the dilatation v, which is very sensitive to Gibbs phenomena occurring somewhere in the computational domain, and Fig. 6b shows the vertical profile of the fluctuation of the dilatation comparing the cases s = qs/16 and s = 0. The stronger dilatation is found near the inversion, where the stronger mixing events occur (remember that the compressible formulation retains the small density variation due to the heat conduction and the mass diffusion that cause this nonzero dilatation), and that is observed in both cases. However, it is clearly seen that the background noise in the case s = 0 is comparable to those values of dilatation at the inversion, whereas the smoothing parameter reduces the background noise by almost two orders of magnitude. This level of accuracy might be needed for small-scale studies, e.g., those based on enstrophy statistics or gradient trajectories [37]. On the other hand, the effect of s in large-scale quantities like the evolution of hb (not shown) is negligible and the maximum difference observed in this quantity is about 1%. Similar comments apply to the smoothing parameter employed in the Boussinesq formulation (cf. Fig. 2 in Sect. 3), as already discussed by Mellado et al. [36]. 6 Conclusions Two major aspects of the problem of buoyancy reversal by isobaric mixing regularly found at the boundaries of clouds have been addressed in this article. First, the existing mixture fraction formulation has been analyzed in detail to understand the underlying assumptions (or simplifications) necessary for its derivation. So doing facilitates a correct interpretation of the results from simulations based on this framework. For that purpose, a general description of the disperse and dilute two-phase system consisting of liquid and vapor in air has been introduced and then particularized to the cloud-top mixing layer, an idealized configuration representing the top of stratocumulus clouds over a size of the order of meters. It has been argued that the main assumptions intrinsic to this simplified methodology are: (a) the liquid phase can be considered as a continuum, (b) local thermodynamic equilibrium exists, and (c) the liquid-phase diffusivity is equal to that of vapor and dry air. Estimates show that these hypotheses do not hold for the conditions normally measured in stratocumulus clouds. The reasons are: (a) the droplet number density is not large enough to have a smooth liquid-phase field over distances comparable to the Kolmogorov scales, (b) there is no water phase equilibrium, and (c) the diffusivity of the water phase due to the thermal motion is significantly smaller than that of the gaseous constituents. Second, because the buoyancy reversal is ultimately caused by molecular mixing and, therefore, it occurs predominantly at the small scales of a system in a turbulent state, high-order numerical algorithms for the compressible and incompressible cases have been described. They are based on finite differences using sixthorder compact Pad schemes for the spatial derivatives and a low-storage five-stage, fourth-order Runge-Kutta scheme for the time advancement. Smoothing parameters have been introduced to deal with the discontinuity in the partial derivatives of the thermodynamic functions at saturation conditions, following previous study on infinitely fast nonpremixed reacting flows. In the incompressible case, the algorithm to solve the Poisson equation has been written using Fourier decomposition in the horizontal planes with a novel technique to solve the resulting one-dimensional equations. Simulations of both compressible and Boussinesq cases have been presented and compared in the cloud-top mixing layer under real atmospheric conditions. Results confirm the expected benefits of the incompressible formulation, which maintains the details of the small diffusive scales of the problem with a substantial reduction in computational cost. Acknowledgements Partial financial support for this study was provided by the Deutsche Forschungsgemeinschaft within the SPP 1276 Metstrm program. 7 Appendix The mathematical framework based on local volume averages [16] is particularized here to the disperse and dilute two-phase flow formed, on the one hand, by the gas mixture of vapor and dry air, and, on the other hand, by the liquid droplets. The union of the gas and liquid phases will be referred to as total mixture or simply mixture. Details about the general formulation can be consulted in the previous reference, and the only objective of this Appendix is to introduce clearly the physical quantities and the conservation equations used in this article. 7.1 Definitions and notation At a given time t , consider the phase functions Xg(x , t ) and Xl(x , t ) which are equal to 1 if the point x is inside the gas phase or liquid phase, respectively, and zero otherwise. The second tool is a volume average operator, denoted by (x, t ), defined locally at x over a volume V independent of x and t and with only one characteristic length scale V 1/3. Then, g(x, t ) = Xg = Vg(x, t )/ V is the volume fraction of gas at (x, t ), i.e., the volume occupied by the gas divided by the total volume V around x defined above. Similarly l(x, t ) = Xl represents the volume fraction of liquid. The relation g + l = 1 follows from Vg + Vl = V . First, let us consider the material partial density fields, respectively, of dry air and vapor, d(x , t ) and v(x , t ), defined on the gas phase so that the gas density is g = d + v, and of the liquid density, l(x , t ), within the liquid phase (it can be ignored for the time being that the liquid density is constant). We define the phase average partial densities as g = Xgg /g, d = Xgd /g and v = Xgv /g. Then, the relation g = d + v holds. Similarly, that of the liquid is defined as l = Xll /l. Further, the bulk density of the total mixture is defined by = Xgg + Xll and represents the total mass divided by the total volume V at (x, t ). The bulk partial densities of vapor, dry air, gas, and liquid are gv, gd, gg, and ll, respectively, and the following relations exist We can define the bulk mass fractions in terms of the bulk densities by qg = gg/ and ql = ll/, which represent the mass of gas and liquid, respectively, divided by the total mass inside the volume V . Then, the relation qg + ql = 1 holds, and the condition g + l = 1 implies that qg/g + ql/l = 1/. The dry air and vapor bulk mass fractions are introduced as qd = gd/ and qv = gv/. These definitions yield g/ = qg/g = qd/d = qv/v. Besides, the condition g = d + v implies qg = qd + qv and 1 = qg + ql = qd + qv + ql. Second, the linear momentum introduces particular velocities into the problem. The continuum description of multicomponent gas mixtures [38,59] usually considers mass-weighted velocities vd(x , t ) and vv(x , t ), such that the mean velocity of the gas mixture is given by gvg = dvd + vvv. Consistently, we define the phase averages by vd = Xgdvd /(gd) for the dry air inside V and vv = Xgvvv /(gv) for the vapor, and similarly for the gas velocity vg = Xggvg /(gg), retaining then the relation gvg = vvv + dvd among volume average quantities. The governing equations are written in terms of gvg for the gas phase and lvl for the liquid phase, and we still need the volume average of this latter, which is given by vl = Xllvl /(ll). The linear momentum per unit volume of the total mixture is used to define the velocity of the mixture v by Xggvg + Xllvl = v , which then yields the following relations v = qgvg + qlvl = qdvd + qvvv + qlvl. e = qgeg + qlel = qded + qvev + qlel. These paragraphs simply recover mathematically the relations between the mass, momentum and energy of the total mixture and the contributions from each constituent, and could have been equally derived using The relative velocities of the gaseous and liquid phases with respect to the total mixture can be written in terms of the mean drift velocity VD = vl vg of the disperse phase with respect to the gas phase both satisfying qgVD,g + qlVD,l = 0 for any given VD. In principle, VD needs to be obtained from an appropriate governing equation. In equilibrium, it is simply the terminal velocity due to the gravitational field or, if droplets are small enough for thermal motion to become relevant, VD would be related to some diffusion velocity (see Sect. 2.3). If a description of the vapor and dry air content separately within the gas phase is desired, then the corresponding relative velocities = Xgd(vd v ) /(gd) = vd v = qv/qg(vv vd) + VD,g, = Xgv(vv v ) /(gv) = vv v = qd/qg(vv vd) + VD,g, will appear in the formulation. In this case, in addition to VD, a closure for the relative velocity vv vd of the vapor continuum with respect to the dry air needs to be provided. The condition d,v,l qi VD,i = 0 is always satisfied for any given pair {VD, vv vd}. The same rationale used before for the velocity can be applied to other specific quantities, like internal energy, enthalpy, or entropy. For instance, let us consider the specific internal energy. From the fields ed(x , t ), ev(x , t ) for dry air and vapor, and the corresponding gas mixture value eg(x , t ), all of them defined inside the gas phase and related through the equation geg = vev + ded, we can the define phase averages ed = Xgded /(gd), ev = Xgvev /(gv), and eg = Xggeg /(gg), which satisfy the relation geg = ded + vev. The same can be done for the liquid, introducing the quantity el = Xllel /(ll). The internal energy per unit volume of the total mixture is Xggeg + Xllel = e, which then leads to the relations only physical arguments. However, this formalism can be further applied to the governing equations of the gas and the liquid phase to derive the governing equations of the mixture and to obtain a clear expression of the unclosed terms in these transport equations. For instance, let us consider the transport equation of the vapor bulk partial density gv = qv. From the transport equation for v (see e.g., [59]) using the properties of the phase function Xg. The term Mv = Xgvvg gvvg is nonzero because vg is defined in terms of g and constitutes a mass flux to be added to the molecular one jv. Both Mv and jv require modeling. The source term Sv = [v(vg vi )] Xg , where vi is the velocity of the interface and the square brackets here indicate the jump across that interface, is unclosed and describes the production (consumption) of vapor at the interface by vaporization (condensation). A similar approach can be done for each of the governing equations for d, l, vg, vl, eg, and el, and unclosed terms appear due to: (a) interface phenomena, (b) nonlinear terms, and (c) molecular transport. In this article, we consider transport equations for the total mixture quantities v and e, so that the corresponding interface terms cancel each other, but the nonlinear and molecular transport of mass, momentum, and energy need to be closed. In addition, the thermodynamic equations of state, which introduce intensive variables, add further closure problems. We cannot extend the same rationale as used before for the specific internal energy and define average values of intensive quantities like the temperature without additional assumptions. In order to proceed further, we need to introduce certain simplifications. 7.2 Simplifying assumptions The disperse two-phase flow is now simplified in the following way: 1. The volume fraction l(x, t ) is assumed to be a smooth enough function to apply a differential formulation over length scales much larger than the droplet diameter d. This requires that the droplet number density nl(x, t ), related to the liquid volume fraction by l = nl d3/6, is high enough for the number of droplets within the volume V to be large, i.e., nl V 1. 2. The volume V is smaller than or comparable to the Kolmogorov scale l3, so that the nonequilibrium within V is mainly due to the perturbation introduced by the dispersed condensate in the form of small enough droplets with diameter d V 1/3. (Larger volumes would fall in the category of large-eddy simulation.) 3. Diluted conditions nl1/3d 1 are also assumed so that the interactions among droplets are neglected. By definition, this condition is equivalent to l 1. In brief, the following chain of inequalities is hypothesized: and the resulting model is commonly referred to in the literature as the two-fluid formulation [13,17,19,22], in contrast to the EulerianLagrangian approach where the droplets are followed individually. The mean drift velocity VD still needs to be provided, and hypothesis (2) and (3) allow us to use the results from the study of the motion of isolated small particles in a nonuniform flow [34], assuming that 4. the droplets behave as spherical rigid solids. Local thermodynamic equilibrium inside V needs to be discussed next. Mechanical and thermal equilibria are assumed when the time scale ratio between the corresponding droplet inertia and the environment evolution is small enough [48]. Mechanical equilibrium would require uniform pressure within the averaging volume V and, therefore, no relative motion of the droplet with respect to the gas should be allowed. This condition is normally relaxed, and a relative velocity of the condensate is considered along with thermal equilibrium, which could be justified by the very low Mach number associated to the droplet slip velocities. Then, the thermodynamic pressure of the gas equals that of the liquid and decouples from the small hydrodynamic pressure field developed by the droplet motion inside V . On the other hand, phase nonequilibrium is more often retained because the corresponding relaxation time scale is larger than the mechanical one [48]. These condensation/vaporization effects lead to source terms in some of the governing equations that require closure [3,24]. The problem can be even more delicate because, then, thermal equilibrium within the volume V is questionable. In this respect, two-temperature models (one for the liquid, and one for the gas) have been proposed [5], and nonequilibrium theories derived for small perturbations from the single-phase state have been also considered [6]. However, the main purpose of this article is to analyze existing equilibrium formulations and consequently This will allow us to use the standard thermal and caloric equations of state for the gas mixture and liquid water. In addition, 6. surface tension effects are neglected. The assumption of local thermodynamic equilibrium and, therefore, uniform values of the scalar fields in terms of x closes also the molecular diffusion fluxes, as later observed. It also sets Mv = 0 in Eq. 60, and the corresponding term associated to the internal energy e appearing in the energy transport equation, but an unclosed term in the momentum and energy equations still remains because of the nonlinearity in the velocities v(x , t ). For instance, in the momentum equation the following term Mv,g = appears due to the gas phase, and, similarly, Mv,l appears due to the liquid phase. A second closure problem occurs also with the nonlinear terms due to a nonzero mean drift velocity VD. All these issues are discussed next as the transport equations that govern the conservation of mass, momentum, and energy of the total mixture are introduced. As explained at the end of Sect. 7.1, these equations are derived by applying the average operator to the governing equations of each phase separately (see e.g., [59]) and then adding the gas and liquid phase contributions when needed. 7.3 Mass conservation Conservation of mass of each constituent yields = 0, = Sv, /t (ql) + (qlvl) = /t (ql) + qlv + qlVD,l = Sl. The source terms Si are due to the phase transition (both satisfying Sv + Sl = 0), jd = d(vd vg) and jv = v(vv vg) are the molecular mass flux vectors (both satisfying jd + jv = 0), and the relative velocities VD,i have been defined in Eqs. 56 and 57 in terms of the mean drift velocity VD and vv vd. Ficks law [59] is used to model the diffusion between vapor and dry air, jv = jv = gv(qv/qg). Linear combinations of the previous three equations can be solved instead of those equations themselves. In the first place, the sum of all of them is considered, which represents the total mass conservation The second combination is the sum of the second and third equations, the transport equation for the total specific humidity where the total-water mass flux vector jt = qvVD,v + qlVD,l is given by jt = gjv + qdqlVD = qgv(qv/qg) + qdqlVD. The ratio between the second and first term on the right-hand side of the previous equation is of the order of (ql/qv)(V D/v) if v represents an appropriate reference velocity scale associated to the diffusion term, e.g., the Kolmogorov velocity scale in a turbulent state. The third conservation equation could be that describing the evolution of vapor or liquid specific humidity, with an appropriate model for the corresponding source term [3,24]. However, by our equilibrium assumption, we can substitute this last transport equation by 7.4 Momentum conservation The conservation of momentum is written as /t qgvg + qlvl + qgvgvg + qlvlvl + Mv = gg + ll + qgg + qlg. In terms of the mixture velocity v , this equation can be expressed as where, for the second equality, the relations g = pg I + g and l = pl I + l have been used (I is the unit tensor), the pressure inside the liquid pl has been taken equal to the gaseous pressure pg due to equilibrium conditions as discussed before, and the droplet is assumed to have no internal relative motion, so that the viscous stress tensor inside the liquid is zero (l = 0). On top of these two assumptions, the fact that l is negligibly small can also be used. The unclosed nonlinear term Mv is included in the stress tensor D and explained below. The viscous stress tensor inside the gas phase for a Newtonian fluid is given by but vg v = VD,g = qlVD, and, therefore, the viscous stress tensor can be written in terms of the transported variable v within a relative error of order ql(V D/v). Last, the stress tensor D is given by D = D,g + D,l, We can compare the order of magnitude of each of these terms with that of the viscous stress, of the order of v2. From the set of simplifications introduced above, it is easy to show that the first term in the right-hand side of the gas contribution is of order l(V D/v)2, whereas that of the liquid is zero, and the last two terms are of order ql2(V D/v)2 and ql(V D/v)2 for the gas and the liquid contributions, respectively. If these ratios are small enough, then the final transport equation is with gg. As it is customary, the pressure and density fields are decomposed into a background profile and a deviation, p(x, t ) = p (x, t ) + pb(x) and (x, t ) = (x, t ) + b(x), such that pb = bg. 7.5 Energy conservation Similar to the conservation of momentum, conservation of energy is expressed as /t qg(eg + vg2/2) + ql(el + vl2/2) + qgvg(eg + vg2/2) + qlvl(el + vl2/2) + Me = qggvg + qlgvl + gg vg + ll vl gjq,g + ljq,l , where vg = vg and vl = vl and jq,g and jq,l are the molecular heat flux vectors. Introducing the velocity of the mixture v , using the continuity equation and assuming that the liquid and gas pressure are equal, then where the substantial derivative operator is defined by d/dt /t + v and tr( D) indicates the trace of the tensor D given by Eq. 72. The enthalpies hi = ei + pi /i (i = d, v, l) have been introduced; it can be shown that qgeg + g p = qghg and equivalently for the liquid. The viscous stress tensor l of the liquid phase is zero because of uniformity inside the droplet and similarly jq,l = 0, and the molecular heat flux vector within the gaseous phase is [59] The flux vector jD is given by jD = jD,g + jD,l, where +VD,g [ Xgg(vg vg)(vg vg) gg + (qghg tr( D,g)/2)I ], +VD,l [ Xll(vl vl)(vl vl) ll + (qlhl tr( D,l)/2)I ]. This expression for the heat flux vector coincides with that presented by Williams [59] for a multicomponent gas mixture, observing that qghg tr( D,l)/2 = qghg+ Xgg vg vg 2 /2 + qg VD,g 2/2, and similarly for the liquid. Simplifications in the limit V D/v 1 are now discussed in the same way as introduced in the momentum equation. The contribution of the tensor D in the explicit terms in Eq. 76 is, at most, of the order of ql(V D/v)2, as already derived in the previous section. The flux vector jD can be similarly simplified. First, the terms involving second- and third-order central moments of the velocity fields inside the volume V are at most of order l(V D/v)2 and l(V D/v)3, respectively. Second, the kinetic energy associated to the trace tr( D) being transported by the diffusion velocities can be neglected in comparison with the enthalpies transported by the same diffusion velocities because the Mach numbers of those motions are very small, and the same is true for the term gg. We are left with /t (e + v2/2) + (e + v2/2)v = ( pv ) + ( v ) jq + g v , where was previously introduced by Eq. 74 (with an error ql(V D/v), as there discussed). The heat flux vector of the total mixture is defined as jq = gjq,g + qgVD,ghg + qlVD,lhl = T + where gg, the relative velocities are defined by Eqs. 56 and 57, and jq,g is given by Eq. 77 with jv = v(vv vg). Using Eq. 64 yields jq = T qgv(qv/qg)(hv hd) + qlqgVD(hl hg). This expression shows that the ratio between the transport of enthalpy by VD and that by mass diffusion between vapor and dry air is proportional to (ql/qv)(V D/v), as observed in Eq. 67 for the mass flux vector jt. However, here we also have a contribution O(hl hg)/O(hv hd) proportional to the ratio between the latent heat and the variation of the sensible enthalpy inside the system; in some cases, like in the cloud-top mixing layer, this ratio is large, and it can control the evolution of the system. 7.6 Thermal and caloric equations of state These equations follow the standard approach in the literature [47], but are written here explicitly for completeness. The pressure has to be related to the other quantities appearing in the previous conservation equations. Consider first the vapor, which obeys the thermal equation of state of an ideal gas, pv = v Rv T = g(v/g) Rv T = g(qv/qg) Rv T = qv Rv T /g. A similar expression is derived for the partial pressure of dry air pd. Therefore, if the volume fraction of liquid is neglected (relative error l/g O (106)), the thermal equation of state, using Daltons law, is The energy of the mixture is written as e = qi ei (T ) (i = d , v, l, surface tension effects are not considered), the gases are taken as calorically perfect, i.e., ed(T ) = ed(T 0) + cv,d(T T 0) = ed0 + cv,d T and ev(T ) = ev(T 0) + cv,v(T T 0) = ev0 + cv,v T , and for the liquid el(T ) = el(T 0) + cl(T T 0) = el,0 + cl T . Then, the caloric equation of state reads e = [(1 qt)cv,d + qtcv,v + ql(cl cv,v)]T + (1 qt)ed0 + qtev ql elv, 0 0 where a heat capacity of the mixture cv can be defined and el0v = ev0 el0 is approximately equal to the latent heat linearly extrapolated to T = 0 K. A reference state ev0 = ed0 = 0 can be used. as expected, and if Eq. 83 is substituted, then where a heat capacity of the mixture cp can be defined. h = [ (1 qt)cp,d + qtcp,v + ql(cl c p,v )]T + (1 qt)ed0 + qtev ql elv, 0 0 Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

This is a preview of a remote PDF:

Juan Pedro Mellado, Bjorn Stevens, Heiko Schmidt, Norbert Peters. Two-fluid formulation of the cloud-top mixing layer for direct numerical simulation, Theoretical and Computational Fluid Dynamics, 2010, 511-536, DOI: 10.1007/s00162-010-0182-x