Controlling entrainment in the smoke cloud using level set-based front tracking

BMeteorologische Zeitschrift,Vol.23, No. 6, 661–674 (published online November 14, 2014) Open Access Article © 2014 The authors Controlling entrainment in the smoke cloud using level set-based front tracking Eckhard Dietze1∗,HeikoSchmidt 1, Bjorn Stevens2and Juan Pedro Mellado2 1Brandenburg University of Technology Cottbus-Senftenberg, Germany 2Max-Planck-Institute for Meteorology, Hamburg, Germany (Manuscript received February 28, 2014; in revised form July 23, 2014; accepted July 24, 2014) Abstract Although large-eddy simulation (LES) has been shown to produce a reasonable representation of the turbulent circulations within the stratocumulus-topped boundary layer, it has difficulties to accurately predict cloud-top entrainment rates. In this paper, we present a front-tracking algorithm for LES to untangle the numerical and physical contributions to entrainment. Instead of resolving the cloud-top inversion, we treat it as a discontinuity separating the boundary layer from the free atmosphere and use the level set method to track its location. We apply our method to the smoke cloud test case as presented by Bretherton et al. (1999) which is simpler than stratocumulus in that it is only driven by radiative cooling avoiding evaporative feedbacks on entrainment. We present three-dimensional LES results with and without use of the level set method varying the grid resolution and the flux limiter. With the level set method, we prescribe zero entrainment and use this case to evaluate our method’s ability to maintain a non-entraining smoke-cloud layer. We use an empiricallybased entrainment law to estimate numerical errors. With the level set method, the prescribed entrainment rate was maintained with errors about one order of magnitude smaller than the entrainment errors found in the standard LES. At the same time, the dependence of the entrainment errors on the choice of the limiter was reduced by more than a factor of 10. Keywords: Entrainment, Stratocumulus, Large-eddy simulation, Level set method 1 Introduction Marine stratocumulus clouds extend over large regions of the subtropical eastern oceans where their annual mean coverage exceeds 40 % (Wood, 2012). Due to their high albedo and high frequency of occurrence, stratocumulus clouds significantly affect Earth’s radiative balance making them one of the climatologically most important cloud systems. However, simulation of such clouds and their associated feedbacks remain key uncertainties in current climate models (Bony, 2005). The wide range of temporal and spatial scales characteristic of these clouds make it still impossible to simulate all important details despite increasing computational resources. The stratocumulus-topped atmospheric boundary layer (STBL) consists of a layer of cool moist air which is capped by relatively warmer and drier air. The boundary layer is well mixed due to turbulent convection and is topped by a stratocumulus cloud. The convection in the boundary layer is driven from the cloud top, mainly by radiative cooling and evaporative cooling. Entrainment at the cloud top is tied to the inversion layer which separates the mixed cloudy layer and the free atmosphere above the cloud. The thickness of this inversion layer ∗Corresponding author: Eckhard Dietze, Brandenburg University of Technology Cottbus-Senftenberg, Siemens-Halske-Ring 14, 03046 Cottbus, Germany, e-mail: [email protected] is measured in metres giving rise to small-scale mixing processes. Large-eddy simulation (LES) has been used extensively to improve understanding of the dynamics of the STBL. However, important quantities such as entrainment remain dependent on the grid resolution and grid spacing aspect ratio, even in recent high-resolution LES with vertical grid spacings of 2.5 m (Yamaguchi and Randall, 2012). The reasons for the difficulty of simulating the STBL with high accuracy using LES are both physical and numerical. The physical aspect is that small-scale mixing processes are not explicitly simulated but typically parameterized using standard closures, the underlying assumptions of which are not satisfied at the cloud top. The numerical aspect is that the steep gradients of total water content and temperature in the inversion layer are insufficiently resolved. As a result, important quantities such as cloud-top entrainment can develop leading-order errors (Dietze et al., 2013). The fact that many physical and numerical uncertainties interact makes it hard to separate their individual contributions and make the solution dependent on details in the numerics as well as in the microphysical and turbulence model (c.f. Moeng et al., 1996). Lilly (1968) introduced a simplified model of the STBL, the smoke cloud. It is similar to the STBL in that the convective boundary layer (CBL) is driven by radiative cooling. It is simpler than the STBL in that it is dry and, thus, avoids evaporative feedback on entrain- © 2014 The authors DOI 10.1127/metz/2014/0595 Gebrüder Borntraeger Science Publishers, Stuttgart, www.borntraeger-cramer.com 662 E. Dietze etal.: Controlling entrainment in the smoke cloud using level set-based front tracking Meteorol. Z., 23, 2015 ment. Instead of water, in the form of vapour or droplets, it contains radiatively active smoke. This simplification has two main advantages. First, it makes the boundary layer problem accessible to experiments. McEwan and Paltridge (1976) as well as Sayler and Breidenthal (1998) conducted tank experiments the latter of which was also numerically reproduced by Schmidt etal. (2012) using a 1D stochastic model. Secondly, it allows more direct investigation of the uncertainty associated with differences in the numerical methods used in LES. Bretherton etal. (1999) carried out LES of the smoke cloud comparing various LES codes. Because they specified the radiation using the identical dependence on the smoke concentration, they could attribute differences among the various codes to differences in the numerical algorithms and the choice of subgrid-scale models. Small-scale studies are also possible. de Lozar and Mellado (2013) carried out direct numerical simulations of the smoke cloud-top interface, but the link to the large-scale boundary-layer motions remains difficult. In this paper, we present and evaluate an algorithm for tracking the cloud-top boundary with the main goal being to untangle the numerical and physical contributions to entrainment. Using the level set method, we represent the cloud-top boundary as a discontinuity on the LES grid and by supplying internal boundary conditions on both sides we avoid discretization over the discontinuity. A similar front-tracking algorithm has been used before by one of the authors in the context of combustion modelling for tracking flame fronts (Schmidt and Klein, 2003). The entrainment process can be included in the form of an additional velocity component in the level set transport. In a similar effort to improve the representation of cloud boundaries and to reduce numerical errors in their vicinity, e.g. Margolin et al. (1997) and Kao et al. (1999) have used the volume-of-fluid method. The main advantage of the level set method is, however, that there is no substantial logic required in order to reconstruct the topology of the interface. Rather, it is given by an isosurface of the associated level set scalar. In this paper, we apply our model to the case of the radiatively driven smoke cloud as presented by Bretherton et al. (1999). We consider the case of vanishing entrainment. The present paper is structured as follows. In section 2, we describe the smoke cloud case, the governing equations, and initial and boundary conditions. We describe our numerical method in the third section,and discuss the simulation results thereafter in section 4.We summarize and present our main conclusions in section 5. 2 Formulation 2.1 Governing equations We consider the smoke cloud case as described by Bretherton et al. (1999). It is a dry CBL filled with radiatively active smoke. The CBL is topped by clear Table 1: Physical parameters Smoke absorptivity Ka0.02 m2kg−1 Specific gas constant of dry air R287 J (kg K)−1 Isobaric heat capacity of dry air cp1004 J (kg K)−1 Gravitational acceleration g9.8 m s−2 Reference pressure p00 1000 hPa Reference temperature Θ0291.5 K Reference density ρ01.1436 kg m−3 and relatively warmer air forming a temperature inversion between the two layers. Convection is solely driven by radiative cooling from the top of the smoke layer. There are no surface fluxes of heat or smoke. We formulate the problem in terms of the potential temperature θand the non-dimensional smoke concentration S.The latter is bounded by 0 and 1 and can also be interpreted as a mixture fraction relating mass of air from the CBL relative to the total mass of a fluid parcel (c.f. Mellado etal., 2010). We use the same 1D, column-wise radiation model as presented by Bretherton et al. (1999) where radiation is incorporated in terms of the vertical radiative flux Frad(z) whose vertical gradient contributes to the temperature tendency according to ∂θ ∂trad =−1 cpρ0 ∂Frad ∂z.(2.1) Assuming the average temperature of the smoke cloud stays close to its initial value over the time considered, the cloud cools at a constant net rate limited by the radiative flux F0at the top of the domain. With this assumption, the radiative flux at any given height zis given by Frad(z)=F0exp −KaH z=z ρ0Sdz .(2.2) Here, His the height of the domain and Kais the smoke absorptivity (c.f. Table 1). The remaining prognostic variables are the Cartesian velocity vector v=(u,v,w) and the dynamic Exner pressure π=(p/p00)R/cp(pis the physical dynamic pressure). The governing equations are the LES-filtered, Ogura-Phillips-type (Ogura and Phillips, 1962) anelastic equations ∂¯v ∂t+1 ρ0∇·(ρ0¯v⊗¯v)=cpΘ0∇¯π+g¯ θ Θ0k +1 ρ0∇·(ρ0τ) (2.3) ∂¯ θ ∂t+1 ρ0∇·(ρ0¯ θ¯v)=1 ρ0∇·(ρ0γθ) −1 ρ0cp ∂Frad ∂z(2.4) Meteorol. Z., 23, 2015 E. Dietze etal.: Controlling entrainment in the smoke cloud using level set-based front tracking 663 Figure 1: Vertical structure and initial conditions of the smoke cloud case. The left side shows a sketch of a vertical cut through the domain. The boundary layer filled with smoke (grey) has an initial depth hand is driven by a constant radiative heat flux F0. The right hand side shows the initial profiles of the non-dimensional smoke concentration S, the potential temperature θ, and the level set function φ. ∂¯ S ∂t+1 ρ0∇·(ρ0¯ S¯v)=1 ρ0∇·(ρ0γS) (2.5) ∇·(ρ0¯v)=0,(2.6) with the bar representing the filter operator. In these equations, ρ0is the hydrostatic density profile of the isentropic atmosphere at temperature Θ0,¯ θis the temperature anomaly (¯ θ−Θ0), and kis the vertical unit normal vector. Table 1shows the values of the physical parameters used here. These are the isobaric heat capacity and the gas constant of dry air, cpand R, respectively, the gravitational acceleration g, the basic state temperature Θ0, and reference pressure p00. From the filtering of the equations, the additional turbulent fluxes of momentum, τ, and scalars, γθ,S, arise. These additional unknowns are approximated using a Smagorinsky-Lilly model (see Stevens et al., 2000). The anelastic continuity equation (2.3) is obeyed by solving the Poisson equation for ¯π ∇·(ρ0∇¯π)=1 cpΘ0∇·−∇·(ρ0¯v⊗¯v) +ρ0¯ θg Θ0k+∇·(ρ0τ).(2.7) In the present case of the radiatively driven smoke cloud, two scalar transport equations are solved, one for S, the non-dimensional smoke concentration, and another one for θ, the potential temperature. Both equations are coupled via the radiative cooling term in the temperature equation. 2.2 Initial and boundary conditions Following Bretherton etal. (1999),weconsidera3D domain extending over 3.2 km in the two horizontal directions and from 0 to 1.25 km in the vertical. The bottom of the domain is filled with smoke over a depth of 700 m (c.f. Fig.1). The initial profiles of the potential temperature θand the smoke Sare: S(z)=⎧⎪⎪⎪⎨⎪⎪⎪⎩ 1ifz∈[0,687.5] m 1−0.04( z m−687.5) K if z∈(687.5,712.5) m 0ifz∈[712.5,1250] m and θ(z)=⎧⎪⎪⎪⎨⎪⎪⎪⎩ 288 K if z∈[0,687.5] m 288 K +0.28 ( z m−687.5) K if z∈(687.5,712.5) m 295 K +10−4(z m−712.5) K if z∈[712.5,1250] m . When the level set method is used, we omit the linear transitional layer and set up the profiles such that they are piecewise constant with only one vertical cell having an intermediate value. In order to accelerate the spin-up of the boundary layer (BL) convection, all temperature values below 650m where perturbed by a spatially uncorrelated uniform random noise. The amplitude of the noise was set to ±0.1K. The boundary layer is initially at rest. The domain is periodic in the two lateral directions and free-slip boundary conditions are imposed at the top and bottom. For scalars, zero-gradient boundary conditions are set at the top and bottom. A sponge layer is used occupying the top ten grid levels. We simulated the smoke cloud for a period of four hours. 3 Numerical methods 3.1 The UCLA-LES As a basis for the numerical algorithm, we use the UCLA-LES. The code has been widely used for simulation of various problems in the realm of atmospheric convection, such as shallow cumulus clouds (Matheou et al., 2011;Stevens and Seifert, 2008;Van- Zanten et al., 2011) and stratocumulus clouds (Ackerman et al., 2009;Stevens et al., 2005)aswellas transitions between cloud types under changing largescale conditions (Bellon and Stevens, 2012;Sandu and Stevens, 2011). We augmented it by a level setbased front-tracking algorithm the details of which we describe in the following section. The UCLA-LES solves the anelastic system of Eqs. (2.3) to (2.5) using finite volumes on a staggered Cartesian grid where velocities are located half a point up-grid in the direction of the respective velocity component from the rest of the variables (θ, S,π). An example grid cell is shown in Fig. 2. For the present case, the grid spacing is equidistant in all three directions. Momentum advection uses a fourth-order directionally-split central method, and scalar advection uses a second-order 664 E. Dietze etal.: Controlling entrainment in the smoke cloud using level set-based front tracking Meteorol. Z., 23, 2015 Figure 2: Sketch of a cell of the staggered grid. The prognostic variables are located at the cell centre (◦), velocities are located at the cell faces ( ) and the level set scalar is located at the cell corners (•). directionally-split upwind method. The code offers various algorithms to monotonize scalar advection. The ones we use in this study are the Minmod, Superbee, and the Monotonized Central (MC) limiter. The flow is advanced in time using a third-order Runge-Kutta method where the new time step is computed as a weighted average of tendencies at intermediate steps. The intermediate steps of any variable, here for the smoke variable S,are calculated according to S∗=Sn+α1Δt∂Sn ∂t(3.1) S∗∗ =S∗+Δtα2∂Sn ∂t+β2∂S∗ ∂t (3.2) Sn+1=S∗∗ +Δtα3∂S∗∗ ∂t+β3∂S∗ ∂t (3.3) where ∂S ∂tare the numerical approximations of the smoke tendencies at the steps n,∗,and∗∗, respectively. The weights are (α1,α 2,α 3)=(8/15,−17/60,3/4) and (β1,β 2,β 3)=(0,5/12,−5/12). Fourier decomposition in the two periodic directions is used to reduce the Poisson Eq. (2.7) to a second-order ordinary differential equation (ODE) in the vertical. The resulting ODE is solved directly to machine accuracy using a tridiagonal solver. 3.2 The level set method Theory Level set methods have been successfully used in a variety of fields to describe the evolution of interfaces. The concept goes back to Osher and Sethian’s (1988) seminal paper and since then many methods have been developed for a variety of problems. These include multiphase and compressible flows, combustion modelling, as well as image processing and computer vision. An overview of the spectrum of applications can be found in Osher and Fedkiw’s (2003) book. Level set methods are based on the idea of implicit surfaces. Rather than explicitly keeping track of a number of points connecting to a surface, level set methods describe the surface geometry in terms of the zero isosurface, or level set, of a space filling scalar function φ(x,t). Hence, the surface φ0is given as the set of implicit points where φvanishes: φ0={x:φ(x)=0}. The function φ(x,t) is referred to as level set function. Once defined, the location of the implicit surface is readily obtained by linear interpolation. Dynamics of the surface are included by solving a transport equation for the level set function ∂φ ∂t+(v+wen)·∇φ=0,(3.4) where wenaccounts for entrainment and nis the normal direction of the interface. Note that in the present paper we consider we≡0. Also note that vand φabove are the LES-filtered variables and the bars have been dropped for convenience. Equation (3.4) is equivalent to the G equation concept used in the combustion literature (e.g. Markstein, 1964,Schmidt and Klein, 2003). Here, we use the level set concept to track the location of the cloud-top interface to the free atmosphere. Specifically, we track the jumps of the smoke and temperature at the cloud top. Both are discontinuous by definition when the level set method is used. The major advantage of the level set method over other interface tracking methods is that the level set function carries specific information about the topology of the interface. 1. The location of the interface is given by the φ=0 isosurface. 2. Both sides of the interface are directly identified by the sign of φ, which divides the domain of interest Ω into the three regions Ω+={x:φ(x,t)>0}(3.5) φ0={x:φ(x,t)=0}(3.6) Ω−={x:φ(x,t)<0}.(3.7) 3. Local information, such as the interface normals n and curvature κare given by the local derivatives n=∇φ |∇φ|,κ=∇·n=∇2φ |∇φ|. In addition, if a relatively smooth level set function is chosen, the interface can be evolved with high accuracy as opposed to the case when scalar fields with steep gradients are advected directly. Discretization The level set Eq. (3.4) is discretized using finite differences on the vertices of the finitevolume grid cells where φis located (see Fig. 2). This location was chosen in order to allow for direct evaluation of the interface intersection points on the cell edges. The interface locations are later used to evaluate face Meteorol. Z., 23, 2015 E. Dietze etal.: Controlling entrainment in the smoke cloud using level set-based front tracking 665 fractions and volume fractions which we explain in section 3.3. First, the velocities at the cell corners – uφ ijk,vφ ijk, and wφ ijk – are linearly interpolated from the LES velocity field. Then, we use a directionally-split upwind method based on second-order upwind polynomials. For example, the first spatial derivative in the xdirection is approximated as ∂φ ∂xi =⎧⎪⎨⎪⎩ 1 2Δx(−1φi−2+4φi−1−3φi)ifuφ i≥0 1 2Δx(3φi−4φi−1+1φi+2)ifuφ i<0(3.8) where indices jand khave been dropped for convenience. The approximations in the other two dimensions are the same with indices and velocity components changed accordingly. The level set Eq. (3.4) is evolved in time using the same third-order Runge-Kutta scheme as used for the other variables (Eqs. (3.1) to (3.3)). Level set reinitialization The level set function, φ, may have any sufficiently smooth shape as long as it satisfies Eqs. (3.5) to (3.7). It is numerically convenient to initialize φinto a signed-distance function of the interface which satisfies the eikonal equation |∇φ|=1. Generally, the level set function will not maintain its signed-distance property as the flow evolves. Especially, in the present case of interfacial convection, the level set function will quickly steepen in the vicinity of the stably stratified interface. This will generate the very difficulties in accurately simulating the evolution of zero level set which we desire to circumvent in the scalar transport. This problem can be avoided by frequently reinitializing φinto a signed-distance function. While the most intuitive way is to directly set φvalues to the shortest distance to the interface, it is more efficient to use a partial differential equation (PDE). Sussman et al. (1994) proposed to iterate the reinitialization equation ∂φ ∂τ =sign(˜ φ)(1 −|∇φ|) (3.9) in virtual time τto steady state. Here, ˜ φis the solution of the current LES time step and held constant during the reinitialization process. They propose discretizing ∇φusing upwind differences, upwind here meaning with a stencil biased towards the interface. Russo and Smereka (2000) note that, applying this method may considerably displace the location of the interface, which they attribute to the choice of the discretization stencil. On points adjacent to the interface, Sussman et al.’s (1994) method would discretize |∇φ| across the interface, locally – only on these points – violating the upwind principle. Russo and Smereka (2000) propose a modification making the method strictly upwind by approximating |∇φ|based on geometrical considerations in points adjacent to the interface. They call their modification the subcell fix and showed (i) that the subcell fix greatly minimizes the spurious displacement of the interface, and (ii) that the maximum displacement error is independent of the number of iterations. In the present work, this is what we use to reinitialize φ. 3.3 Level set/FV coupling So far, the evolution of the prognostic variables and the level set function are coupled only one-way, i.e. the level set is passively advected with the flow. In order to couple the evolution of the prognostic variables with the level set, we follow the concept by Smiljanovski et al. (1997). In their compressible framework, they used a level set method to track the position of a flame front, and they used Rankine-Hugoniot jump conditions to supply internal boundary conditions at the interface. The reconstruction allows for computing fluxes and source terms associated with the two states individually. The individual fluxes and source terms can then be superimposed to obtain the net effect in the particular cell. In our present anelastic system, we retain only the latter part of the approach and use Fedkiw et al.’s (1999) ghost fluid method to supply internal boundary conditions. In the following we describe the process of coupling the FV method with the level set. It can be divided into three steps: 1. Reconstruct the two coexisting states (S,θ)0,(S,θ)1. 2. Evaluate fluxes and source terms associated with the two states. 3. Compute total fluxes. In the present formulation, we apply the method only to scalar transport where we, so far, only couple advective fluxes. For the purpose of illustrating the method in detail, let us consider only the smoke S; the same method is applied to the potential temperature θ. In step one, we use the ghost fluid method by Fedkiw et al. (1999). The idea is to generate ghost fluids on both sides of the interface in order to supply internal boundary conditions at the interface. Fedkiw et al. (1999) generate the ghost fluids by iterating a PDE in virtual time extrapolating information across the interface. After the process, there are two fluids present on both sides of the interface: one physical fluid and one ghost. Let us look at the generation of the ghost fluid in the top region. For this, we carry over the notation of Eqs. (3.5) to (3.7) into the discrete sense: We denote the set of cut cells as Ω0, the set of cells for which all φ>0asΩ+, and the set of cells for which all φ<0asΩ−. First, the physical field Sis copied to a new one, S0, serving as initial condition. Then, in the top region (Ω0Ω+), i.e. in all cut cells and the uncut cells above the interface, we iterate the extrapolation equation ∂S0 ∂τ +n·∇S0=0 (3.10) 666 E. Dietze etal.: Controlling entrainment in the smoke cloud using level set-based front tracking Meteorol. Z., 23, 2015 Figure 3: Sketch of a cut cell to illustrate the ghost fluid method (top) and flux superposition (bottom). The dashed line indicates the location of the interface. Shaded regions are immersed in smoke. The primed indices abbreviate i=i+1/2andj=j+1/2. in virtual time τwith the values in Ω−serving as Dirichlet boundary conditions. Having defined φso it is positive above the interface, this carries information at speed 1 from the bottom upwards. After the process, the field S0contains physical values of Sin Ω−and ghost values in (Ω0Ω+). The process is then repeated in (Ω0Ω−) to generate S1from another copy of S. The only difference is that the plus sign in Eq. (3.10) is inverted to a minus sign so information travels from the top downwards. The same is done for the potential temperature θ, so that we end up with (S,θ)0and (S,θ)1 both of which consist of physical and ghost fluid. The index signifies the origin of the ghost information: Index 0 stands for “originating from the bottom”, and index 1 stands for “originating from the top” (see the top Fig. 3). In the second step, following Smiljanovski et al. (1997), we compute advective fluxes based on the extrapolated fields. Consider the smoke transport Eq. (2.3) in 2D. The integration over one grid cell Vij and the application of Gauss’s theorem lets us write the integral form ∂S ∂t=Fi−1,j−Fi,j Δx+Gi,j−1−Gi,j Δy(3.11) where Sis the cell average of the smoke and Fi,jand Gi,jare the total horizontal and vertical fluxes of Sat the right and top cell faces, respectively. The unprimed indices signify cell centred locations, and the primed indices signify cell face locations (abbreviating i= i+1/2andj=j+1/2). Based on the extrapolated fields S0and S1, we compute fluxes associated with the cloudy air, F(S0), and the free atmosphere, F(S1), individually. We obtain the total flux as the weighted sum Fij=(1 −˜ βij)Fij(S0)+˜ βijFij(S1),(3.12) where the ˜ βijare the face fractions at the respective cell face (see Fig. 3) time-averaged over one LES time step. This approach conveniently enables us to use the same flux scheme both at the interface as well as in the interior of the flow. 3.4 Level set synchronisation An interesting detail of the method so far is its overdetermination in terms of how the interface is described. By introducing the level set Eq. (3.4) we added one more equation describing the evolution of the interface which, in terms of this information, is redundant with the scalar transport equations. If the location of the interface φ0 and the resulting volume fractions and face fractions are not synchronised with the interface location represented by the scalar fields, errors in the form of overshoots and undershoots may accumulate to leading order. One approach to resolve this redundancy is to combine the level set method with volume-of-fluid techniques. Schneider (2001) formulates a correction method in terms of an elliptical equation connecting all cells along the interface. The solution of the elliptical equation are corrections of the face fractions to be used in Eq. (3.12) which distribute the corrections along the cut cells. Similar techniques are subject of current research (Waidmann, 2013, private communication). A second way is to let the level set govern the location of the interface and reassign cell averages in cells cut during an LES time step according to the current volume fraction given by φ0: S θcut =αcut S θ1 +(1 −αcut)S θ0 Here, (S,θ)0,1are the extrapolated fields as discussed in Sec. 3.3. In the present case where there is no mixing across the interface (we=0), we know the exact reconstruction in cut cells which are S0≡1andS1≡0which we use for the synchronisation instead of the extrapolated fields. This is the method we use in the present paper. The method is local and, thus, computationally simpler than the first approach. However, note that this correction method is not globally conservative. The smoke mass change in the simulations presented below was +1.4%and+0.8% in the low and high resolution case, respectively. 4 Simulations and results 4.1 Goals and setup We present results of our simulations of the smoke cloud using the Level Set LES (LS-LES) and the standard LES as a reference. The two main questions we address in this paper are: 1. How accurately does the LS-LES maintain the prescribed zero entrainment and decouple the entrainment process from the BL convection? This is a requirement for (super-)parameterizing entrainment. Meteorol. Z., 23, 2015 E. Dietze etal.: Controlling entrainment in the smoke cloud using level set-based front tracking 667 Figure 4: The smoke cloud after 3 hours as simulated by the LS-LES at the double resolution. The grey sheet indicates the location of φ0. Colours on the left slice show the smoke concentration, colours on the right slice show potential temperatures. Both slices are overlaid with streamlines in the respective planes colour-coded by velocity magnitude. In the background is an arbitrary potential temperature isosurface in red (θ≈287.3 K). 2. Does the method minimize the dependency of flow statistics on grid resolution and details of numerical methods used? We hypothesize, that by avoiding discretization over the interface, part of the problem can be removed. For this, we ran a series of simulations where we modify two parameters. The first one is grid resolution. We use equidistant grids with Bretherton et al.’s (1999) standard resolution of 64 ×64 ×50 grid cells, denoted ‘S’, as well as double that resolution which we denote by ‘D’. The standard resolution corresponds to a grid spacing of 50 m in the horizontal directions and 25 m in the vertical direction. The double resolution reduces these numbers by a factor of 2. The second parameter is the choice of the flux limiter. We use the Minmod, Superbee, and Monotonized Central (MC) limiter. These two parameters span a total of 12 simulations, 6 for both the standard LES and the LS-LES. In addition, we ran two simulations with the standard LES with the limiter switched off. Figure 4shows a snapshot of the ‘D’ LS-LES after 3 hours with the grey sheet indicating the φ=0 isosurface. The instantaneous streamlines and the θisosurface (red) indicate a complex turbulent flow in the BL. At the same time, there is a large-scale motion with strong vertical flow in the BL interior and horizontal redirection near the cloud top and bottom. In the following sections, we discuss horizontally averaged statistics where we specifically focus on entrainment and how it is affected by the choice of the numerical parameters. In addition, we look at their effect on the turbulent kinetic energy (TKE), its evolution, and distribution in the vertical and horizontal components. For a more extensive discussion of the details of the flow we refer the reader to Bretherton et al.’s (1999) original paper. The entrainment rate, in the context of interfacial convection, is typically defined as the time derivative of the height of an interface, zi, chosen to separate the turbulent rotational flow from the irrotational flow: we=d dtzi. Here, we define zias the height of the S=0.5iso- surface and  denotes the horizontal averaging operator. In order to evaluate the accuracy at which the LSLES maintains the prescribed entrainment, we use the simplest possible case: We prescribe zero entrainment (we≡0inEq.(3.4)). Thus, we separate numerical errors from uncertainties in the modelling of entrainment and we can attribute any net entrainment to numerical errors. In both cases, the standard LES and the LS-LES, an estimate for the exact entrainment is needed in order to define the error and a common reference entrainment. Here, we use Bretherton et al.’s (1999) formula wth e(A)=1.25A 1+1.25A gF0(1 −(2/Kaρ0zi)) ρ0cpΘ0Δb(4.1) as an estimate (‘th’ stands for theoretical). It is derived from Sayler and Breidenthal’s (1998) Richardson number scaling using a mixed layer model. Ais an 668 E. Dietze etal.: Controlling entrainment in the smoke cloud using level set-based front tracking Meteorol. Z., 23, 2015 Figure 5: Evolution of the horizontally averaged cloud top height zicomparing various limiter choices for the standard LES. ziis defined as the height of the S=0.5 isosurface. Left: standard resolution, right: double resolution. empirical model parameter with values of 0.2 to 0.4 as suggested by Sayler and Breidenthal’s (1998) experiments. For the physical parameters, we use the values presented in Table 1and as the buoyancy jump Δbwe use the value 0.25 m s−2which is the average buoyancy difference across the inversion during the third hour of the D simulations. For the given range of the A parameter, Eq. (4.1) predicts entrainment rates between wth e,min =wth e(0.2) =1.24mm s−1and wth e,max =wth e(0.4) =2.06mm s−1. With the LES generally overestimating entrainment, the two values define the worst case and best case errors, respectively. We will use these two values as a reference to define two relative errors in the form Δwe wth e =we−wth e wth e .(4.2) 4.2 Standard LES In Fig. 5, we show the evolution of the inversion height for the standard LES for the S and D grid and the above mentioned limiter choices. All simulations show an initial transient, after which the BL depth grows quasi linearly with superimposed oscillations with a time-scale of one to two hours. The initial transient as well as the amplitude of the oscillations are reduced in the doubleresolution case. Relative to the results obtained with the MC limiter, inversion heights are consistently higher with the Minmod limiter and lower with the Superbee limiter. Due to the oscillations in the standard-resolution case, this order is not always reflected in the hourly averaged entrainment rates. In order to ensure comparability with Bretherton et al.’s (1999) results, we look at statistics averaged over the third hour (i.e. 2h-to-3 h averages). However, as noted by those authors, the oscillations in the standard resolution runs render statistics averaged over just one hour less reliable. For this reason we also include 2 h-to-4 h averages in our considerations. Average entrainment rates, as presented in Table 2, are in agreement with Bretherton etal.’s (1999) intercomparison. Across our simulations, we measure entrainment rates between 3.1 and 4.1mm s−1which is within the measured range of the 3D runs of the intercomparison. As mentioned before, there is a consistent correlation between higher and lower entrainment and the limiter used. This is best seen in the 2h-to- 4h averages in the second column. Minmod consistently produces greater entrainment than MC (+7.37% and +13.41% in the S and D run, respectively), and Superbee produces lower entrainment (−6.12 % and −10.51%). We observe a trend of increasing entrainment rates over time, consistent with Bretherton etal.’s (1999) results, which results from the gradual weakening of the inversion and a slight increase of cloud-top mixing due to increasing TKE over time. The increase of the entrainment rate is stronger in the standard-resolution simulations, which in part may be attributed to the phase of the oscillations of the inversion height and in part to greater weakening of the inversion on the coarser grid. The last two columns show how entrainment rates change with grid resolution. Entrainment rates reduce as the resolution is increased which is also consistent with Bretherton et al.’s (1999) observation. The reductions of more than 10 % in the 2-to-3-hours averages should be considered less reliable than the 2-to-4- hours averages due to the above mentioned oscillations. The trend, however, clearly remains when considering 2h-to-4 h averages. Figure 6shows the evolution of the horizontal average of the vertically integrated TKE density (VTKE) for the standard LES. The subgrid-scale part of the TKE is diagnosed from the filtered velocity field consistent with the Smagorinsky model (see e.g. Stevens et al., 2000). As with the evolution of the inversion height, the VTKE Meteorol. Z., 23, 2015 E. Dietze etal.: Controlling entrainment in the smoke cloud using level set-based front tracking 669 Figure 6: Evolution of the horizontal average of the vertically integrated TKE density (resolved + subgrid-scale) comparing various limiter choices for the standard LES. Left: standard resolution, right: double resolution. Table 2: Time-averaged entrainment rates we(left) and their relative changes (right) with grid resolution and flux limiter choice for the standard LES. The percentages in brackets are the relative changes with respect to the MC limiter. we[mm s−1] Relative changes (I) 2 h–3 h average (II) 2 h–4 h average (I) →(II) S →D (2 h–3 h avg.) S →D (2 h–4 h avg.) (S) Standard resolution Minmod 3.4817 (−0.74 %) 4.1339 (+7.37%) +18.73 % MC 3.5076  3.8503  +9.77 % Superbee 3.3815 (−3.60 %) 3.6148 (−6.12 %) +6.90 % (D) Double resolution Minmod 4.0396 (+15.84 %) 4.0727 (+13.41 %) +0.82 % −14.22 % −1.48% MC 3.4841  3.5565  +2.08 % −13.45 % −7.63% Superbee 3.0984 (−11.00 %) 3.1519 (−10.51 %) +1.73 % −16.70 % −12.80 % exhibits an initial transient and a decaying oscillation. It eventually settles at a relatively stable magnitude after about 1.5 to 2 hours. In Table 3, we list 2 h-to-4 h averaged values. We observe a consistent ordering for the three limiters with Superbee exceeding the MC limiter by +10.12% (S) and +9.47% (D), respectively, and Minmod going below it by −23.6% (S) and −12,88% (D), respectively. This trend continues to exist in the individual profiles of the horizontal and vertical velocity variances which are shown in Fig. 7. The non-vanishing velocity variances above the cloud layer in the no-limiter case are due to scalar overshoots caused by the nonmonotone advection scheme. 4.3 LS-LES In Fig. 8, we show the evolution of the inversion height for the LS-LES. It exhibits a similar initial transient over a period of roughly 30 minutes after which it evolves in a quasi linear way. While there is a dependence on the limiter visible, it is minimal compared to the standard LES runs. Overall, entrainment is drastically reduced. In Table 4we compare the 2 h-to-4 h averaged entrainment rates for both the standard LES and the LSLES with the rates predicted by Eq. (4.1). On the left Table 3: 2 h-to-4 h averages of the vertically integrated TKE VTKE [kgs−2] Limiter Standard LES LS-LES (S) Standard resolution Minmod 244.80 580.08 MC 320.40 533.94 Superbee 352.83 572.76 (D) Double resolution Minmod 302.11 579.63 MC 346.78 541.30 Superbee 379.62 535.79 side are the entrainment rates compared to the lower entrainment prediction; on the right side we compare to the greater predicted value. We focus on the maximum and minimum values from the high resolution ‘D’ runs. This range of entrainment rates is representative of the overall range observed with the standard LES. It also contains the lowest one observed which is closest to theoretical predictions. The combination of lowest and highest observed entrainment and lowest and highest predicted entrainment defines the four cases for which the relative errors are presented.