Study of low-order numerical effects in the two-dimensional cloud-top mixing layer

Theor. Comput. Fluid Dyn. (2013) 27:239–251 DOI 10.1007/s00162-012-0263-0 ORIGINAL ARTICLE Eckhard Dietze ·Juan Pedro Mellado ·Bjorn Stevens · Heiko Schmidt Study of low-order numerical effects in the two-dimensional cloud-top mixing layer Received: 30 July 2011 / Accepted: 7 March 2012 / Published online: 11 April 2012 © The Author(s) 2012. This article is published with open access at Springerlink.com Abstract Large-eddy simulation (LES) has been extensively used as a tool to understand how various processes contribute to the dynamics of the stratocumulus layer. These studies are complicated by the fact that many processes are tied to the dynamics of the stably stratified interface that caps the stratocumulus layer, and which is inadequately resolved by LES. Recent direct numerical simulations (DNS) of isobaric mixing due to buoyancy reversal in a cloud-top mixing layer show that molecular effects are in some instances important in setting the cloud-top entrainment rate, which in turn influences the global development of the layer. This suggests that traditional LES are fundamentally incapable of representing cloud-top processes that depend on buoyancy reversal and that numerical artefacts can affect significantly the results. In this study, we investigate a central aspect of this issue by developing a test case that embodies important features of the buoyancy-reversing cloud-top layer. So doing facilitates a one-to-one comparison of the numerical algorithms typical of LES and DNS codes in a well-established case. We focus on the numerical effects only by switching off the subgridscale model in the LES code and using instead a molecular viscosity. We systematically refine the numerical grid and quantify numerical errors, validate convergence and assess computational efficiency of the low-order LES code compared to the high-order DNS. We show that the high-order scheme solves the cloud-top problem computationally more efficiently. On that basis, we suggest that the use of higher-order schemes might be more attractive than further increasing resolution to improve the representation of stratocumulus in LES. Keywords Cloud-top mixing layer ·Buoyancy reversal instability ·Large-eddy simulation 1 Introduction Marine boundary layer clouds, of which marine stratocumulus is an important archetype, have been identified as one of the primary sources for the uncertainty in model-based estimates of the equilibrium climate sensitivity. Thus, an improved representation of these clouds is thought to substantially increase the reliability of climate change predictions [1]. One characteristic feature of marine stratocumulus that challenges our ability to model them is the strongly stratified layer at the cloud top. This layer has a scale of metres over which the temperature gradient can be positive (keeping in mind the temperature gradient for neutral stratification is −10 Kkm−1) giving rise to the terminology inversion layer, as the temperature gradient is inverted relative to Communicated by R. Klein. E. Dietze (B )·H. Schmidt Faculty of Mechanical, Electrical and Industrial Engineering, Cottbus University of Technology, Siemens-Halske-Ring 14, 03046 Cottbus, Germany E-mail: [email protected]us.de J. P. Mellado ·B. Stevens Max Planck Institute for Meteorology, Hamburg, Germany 240 E. Dietze et al. the tropospheric norm. Already the mixed-layer theory by Lilly [5] identifies cloud-top entrainment, associated with the mixing of air across this region of strong stratification, as an important determinant for the evolution of the cloud. The main processes that drive entrainment are (1) radiative cooling, (2) shear due to large-scale boundary layer vortices and (3) evaporative cooling. Above the cloud top of subtropical marine stratocumulus, the environment is relatively warm and dry, as it is characteristic of the large-scale environment of subsiding air in which the stratocumulus-topped boundary layer (STBL) prevails. The STBL itself is more characteristic of the underlying ocean surface in that it is much cooler and moister, and as the name indicates, cloud-topped. On the one hand, warmer air that is mixed into a cloud-top parcel raises its temperature. On the other hand, because the warmer air is drier, evaporation takes place which cools the mixture. Both effects compete and if the latter dominates, mixtures can develop that are heavier than either one of the two unmixed states. This effect is known as buoyancy reversal [16,23] which in this case is caused by evaporative cooling. In past studies of this problem, the influence of effects at the molecular scale and below has largely been ignored; these include droplet dynamics, phase changes and molecular mixing. The question whether buoyancy reversal alone can cause a runaway instability by positive feedback of evaporative cooling, mixing, and entrainment leading to cloud break-up [2,14] has long been debated. The strength of the reversal is typically expressed by the buoyancy reversal parameter D=−bs/b1[16] that relates minimum and maximum buoyancy bsand b1, respectively. This parameter is typically small. For instance, in the first research flight (RF01) of the DYCOMS II field study [18], a buoyancy reversal parameter of D=0.031 was measured. Large-eddy simulation (LES) has been extensively used as a tool to improve understanding of the processes at the cloud top. The fidelity of LES has increased tremendously over the last decades. Spatial resolution increased from 50 m in the early studies [2] to 5 m and less in recent LES [11,19,22,24]. Unfortunately, general circulation models and even LES currently have, and in the foreseeable future will have, insufficient grid resolution to resolve the cloud top and the small-scale mixing processes that give rise to buoyancy reversal. Even at the relatively high-resolution characteristics of present studies, the evolution of stratocumulus remains sensitive to the parameterization of subgrid-scale (SGS) processes, or the numerical representation thereof. In some cases, there are order one errors in important properties of simulated stratocumulus such as the entrainment rate. By making an effort to limit mixing, for instance by switching the SGS terms off and using numerical limiters as a closure model, the results improve in some cases [19] but for unphysical reasons. This is unsatisfactory, since the ‘good’ comparison depends among other on the grid size and the order of the numerical method. This suggests that (1) SGS processes play a significant role in the evolution of the cloud and that (2) numerical errors are of leading order, thus making it hard to distinguish between numerical and physical effects. This paper studies this second aspect of the problem, that is, low-order numerical effects. High-resolution methods, which allow for a representation of molecular effects, have been used to improve the understanding of the entrainment mechanism. Examples include one-dimensional stochastic models such as the linear eddy model [3] and the one-dimensional turbulence model [23] but also direct numerical simulations (DNS), based on the Boussinesq approximation and a continuous two-fluid formulation [10]. Mellado et al. [8] show by means of two-dimensional DNS that complex flow patterns can be generated by weak buoyancy reversal typical for stratocumulus. The DNS also shows that the turbulent mixing dominates the downward development of the cloudy layer and that there is no significant enhancement of turbulent entrainment of upper fluid in the shear-free case. Furthermore, linear stability analysis identifies two time scales: the time scale associated with the downdraughts of heavy fluid due to buoyancy reversal and the time scale of the restoring force of the inversion. The ratio between both is √D, and the conjecture is that for typical conditions (0 <D1) buoyancy reversal alone cannot dissipate the cloud because the inversion returns to equilibrium much faster than the time required by downdraughts to descend deep enough into the cloud layer. Three-dimensional DNS [7] confirm this conjecture. In the simulation, buoyancy reversal promotes convection leading to turbulence with an intensity that is too weak to overturn the inversion. The system evolves into a self-preserving state that retains the vertically layered structure. Above the convection layer, a thin diffusion layer with a constant thickness h=κ/weforms and travels upwards at a constant mean speed we∼(κ|bsχ2 c)1/3,whereκis the molecular diffusivity and χcis the crossover mixture fraction, equal to the interval of negatively buoyant mixtures. Analysis of the budgets of turbulent kinetic energy shows that a strong conversion of vertical to horizontal momentum occurs in a thin region between the inversion layer and the convection layer [7,10]. Probability density functions of velocity in the turbulent region show a strong vertical variability and often non-Gaussian behaviour, even in the horizontal component [9]. This indicates strong anisotropy and inhomogeneity of the turbulence which, in addition to the relevance of molecular transport, further complicates turbulence modelling in LES for this kind of flows. Study of low-order numerical effects in the two-dimensional 241 Despite the obvious advantages of using DNS, bounded computational resources mean that LES is still required for global, complex geophysical systems like the STBL. As a step towards better LES of stratocumulus, we study the behaviour of a typical LES numerical algorithm for a test problem that we have defined to capture the processes thought to be relevant to buoyancy reversal and its contribution to cloud-top entrainment. We choose the two-dimensional cloud-top mixing problem by Mellado et al. [8] and use their DNS solution as a reference. This particular case is explored because it embodies a number of important features such as a stably oscillating inversion, molecular diffusion, baroclinic production of vorticity, and evaporative cooling leading to buoyancy reversal, which merit the introduction of a new test problem for small-scale atmospheric flow solvers. We use the University of California, Los Angeles LES (UCLA-LES) [19–21] where we switched off the SGS model and introduced a constant molecular viscosity. It is emphasized that, in this paper, the term ‘LES’ will be used solely to refer to the UCLA-LES numerical algorithm. Where, in contrast, LES modelling is meant it will be clear from the context. The reference DNS uses a sixth-order compact scheme, whereas in LES low-order methods are commonly used, typically second-order as in this study. The emphasis on low-order effects on LES scales arises because the flow is not smooth on the grid scale, thereby obviating the advantages of higher-order methods. Low-order methods also use a smaller stencil and thus are more easily to implement efficiently in a highly parallel environment. What we are interested in here is how sensitive solutions are to numerical errors and what the resolution requirements are as compared to a high-order DNS solver. For that, we carry out a grid convergence study where we systematically refine the computational grid until convergence is achieved. By doing that, we also verify the UCLA-LES code for the test case. Because the choice of highversus low-order methods often involves trade-offs, we also investigate how the computational costs of both codes compare. The remainder of the paper is organized as follows. In Sect. 2, the formulation of the test case is described. In Sect. 3, we briefly describe the numerical algorithms of both the reference DNS and the LES code and present results of the convergence study showing self-convergence of the LES code and convergence of the LES code to the DNS solution. Two lengths, hband ht, are introduced, which quantitatively characterize the evolution of the flow, and are used to quantify convergence. Also, we address computational efficiency by comparing the computational cost of the DNS and the LES code packages as they are normally used. Section 4 concludes the present contribution by summarizing the main findings and discussing the implications of this study for LES of stratocumulus in general. 2 Formulation 2.1 Reference problem We study a two-layer system in which a cold moist layer, which is saturated and has liquid water, is topped by a relatively warmer, drier layer. In the mixing region between the two layers, the total water mixing ratio qtand specific enthalpy htransition between the two values characteristic of the upper unsaturated and lower saturated layer, respectively. Initially, the system is at rest, gravity acting downwards. The system so defined is a simplified surrogate of the stratocumulus top in order to investigate certain physical phenomena associated with relatively small-scale dynamics. We describe the multiphase flow using a two-fluid formulation. The main underlying hypotheses are the following: (1) The liquid water can be considered a continuous phase, (2) the mixture always stays in local thermodynamic equilibrium and (3) the diffusivity of the liquid phase equals that of dry air and water vapour. These assumptions are not identically satisfied under real conditions in stratocumulus but provide a good approximation to study latent heat effects. A more thorough discussion of these hypotheses and their applicability can be found in [10]. Assuming equal thermal and mass diffusivities, the thermodynamic coordinates hand qtobey the same advection-diffusion equation. This implies that, with appropriate boundary conditions, both quantities have the same normalized solution after an initial transient. Thus, their evolution can be described by that of a single conserved scalar—the mixture fraction χ. With the indices 0 and 1 indicating the states in the lower and upper layer, respectively, the mixture fraction can be expressed in terms of hand qtas χ(x,t)=h(x,t)−h0 h1−h0=qt(x,t)−qt,0 qt,1−qt,0 ,(1) Herein, x=(x,z)Tdenotes the location and tis the time coordinate. Hence, the mixture fraction simply measures the fraction of mass of fluid from the top layer in an arbitrary fluid parcel in the flow; by definition, 242 E. Dietze et al. -0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b/b1 χ -D χs RF01-DYCOMS II strongly reversing case non-reversing case Fig. 1 Non-dimensional buoyancy mixing function for the present case (solid), and two fictional cases: a non-reversing (dashed) and a strongly reversing case (dotted). The latter one corresponds to the case A3 by Mellado et al. [8] χ=0 in the bottom layer and χ=1 in the top layer. The thermodynamic state of the system, in particular the density field, is then completely defined by the field of χ. We use the Boussinesq approximation and describe gravitational forces in terms of the buoyancy b(x,t)=ρ0−ρ(x,t) ρ0 g,(2) gdenoting the magnitude of the gravitational acceleration and ρis the density. This is justified because differences in density across the cloud-top inversion are small. Depending on the initial data, evaporative cooling can cause saturated mixtures to be heavier than either one of the unmixed reference states. The term ‘buoyancy reversal’ is introduced to identify this eventuality [16,23]. The heaviest mixture 0 ≤χs≤1 occurs at saturation conditions and has buoyancy bs=ρ0−ρs ρ0g. Shy and Breidenthal [15] introduce the buoyancy reversal parameter D=ρs−ρ0 ρ1−ρ0=−bs b1 to describe the strength of the reversal, relating buoyancy of the heaviest mixture to the buoyancy jump b1=ρ0−ρ1 ρ0gacross the inversion. With the above assumptions and for small density differences ρ0−ρ1, buoyancy can be accurately approximated by a piecewise linear function of the mixture fraction. Here, the formulation b(χ)/b1=−D χsχ+1+D 1−χs+D χsδsln exp χ−χs δs+1(3) introduced by Mellado et al. [8] is adopted. Figure 1shows three examples of the graph of b(χ): namely a non-reversing and a strongly reversing case for illustrative purposes, and the buoyancy mixing function for the present case. The parameter δswas introduced to smooth the discontinuity between the piecewise linear sections. Obviously, the smoothing damps the maximum negative buoyancy −bsand, thus, the effective reversal parameter D. Previous work [8] shows that a smoothing parameter of δs=χs/16 leads to a sufficient independence of the solution from changes in δs. This value of δsis adopted here for consistency in the comparison. 2.2 Model equations With the assumptions made in Sect. 2.1, the set of prognostic variables in two dimensions consists of the two Cartesian velocity components u=(u,w) T, the kinematic pressure pand the mixture fraction χ. Study of low-order numerical effects in the two-dimensional 243 -5 -2.5 0 2.5 5 0 0.5 1 (z - z0)/δ χ, b/b1 mixture fraction χ buoyancy b/b1 Fig. 2 Vertical profiles of the mixture fraction and buoyancy as following from Eqs. (3)and(8) The reference DNS solves the Boussinesq limit of the Navier-Stokes equations, given by ∂u ∂t+u·∇u=−∇p+ν∇2u+b(χ)k(4) ∂χ ∂t+u·∇χ=ν∇2χ(5) ∇·u=0.(6) The LES code, whose native formulation solves Ogura-Phillips [12] type anelastic equations, was modified to solve exactly the same system. In Eq. (4), νis the constant kinematic viscosity of the mixture. The buoyancy mixing function b(χ) is given by Eq. (3), kis the vertical unit vector, and pis a modified pressure, normalized by the density. The Poisson equation for the dynamic pressure ∇2p=∇·−u·∇u+ν∇2u+b(χ)k(7) couples conservation of mass, see Eq. (6), to the momentum Eq. (4) and completes the system of equations. 2.3 Initial and boundary conditions We consider a density-stratified system that is at rest in the beginning. Profiles of specific enthalpy and total water mixing ratio change in the vertical from values characteristic of the saturated and unsaturated layer over a height of the order of 4δfollowing an error function. Employing the mixture fraction formulation, Eq. (1) allows for a specification of the initial thermodynamic conditions by χ(t=0,x,z)=1 21+erf z−z0(x) 2δ.(8) Herein, z0denotes the vertical position of maximum stratification that corresponds to the height where χ=0.5. The combination of Eqs. (3)and(8) gives the vertical profile of the buoyancy. Both profiles are plotted in Fig. 2. In a non-reversing scenario (cf. dashed graph of Fig. 1), the normalized buoyancy would match the mixture fraction profile. In contrast, under conditions of buoyancy reversal, a small peak to the negative buoyancy axis develops just below the χ=0.5 isoline. Thus, the system is unstable with respect to the lower fluid in a thin layer. The instability process is triggered by imposing one mode of a sinusoidal wave on the inversion height z0 with the wave length λand amplitude a/2 according to z0(x)=a 2cos 2πx λ. 244 E. Dietze et al. Fig. 3 Time series of the buoyancy field at times t/T=(0.0,1.05,2.18,3.32,4.45)with T=2√π√λ/b1, colours range from zero (white)to−b1(black). UCLA-LES, 512 ×1,024 grid cells Here, ameasures the thickness of the perturbed mixing layer. The ratio of amplitude to wave length is (a/2)/λ =0.1, and the initial thickness of the error function profile is δ/λ =0.025. We consider a twodimensional rectangular domain with the extent (x, z)=(λ, 2λ). The mean height of the inversion, given by the mean height of the χ=0.5 isoline, is located at 60 % of the vertical domain. The leftmost plot in Fig. 3 shows the initial situation. We consider conditions that correspond to the STBL as observed in the first research flight of the second DYCOMS field programme [18]. The data show weak buoyancy reversal with a buoyancy reversal parameter of D=0.031. The mixture fraction at saturation conditions is χs=0.09, which completes the description of the buoyancy function (3)(cf.Fig.1). Mellado et al. [8] introduce the reference Grashof number Gr =a4b1 ν2λ, which can be interpreted as the square of the ratio of two characteristic time scales. The first one is the viscous time that scales as a2/ν, the second one is the period of the oscillation of the interface and is proportional to √λ/b1[8]. For the reference simulation, the kinematic viscosity ν=1.5×10−5m2/sandGr =6.4×105 were used. That corresponds to, keeping the above ratio of (a/2)/λ in mind, a wave length λ≈0.71 m of the sinusoidal perturbation. We impose cyclic boundary conditions at the lateral and free-slip conditions at the top and bottom boundaries. For scalars, zero gradients in the vertical are set at the top and bottom. The simulations are stopped before finite-size effects associated with the computational domain become important. 3 Simulations 3.1 The numerical methods We use the cloud-top problem defined in Sect. 2as a test case to investigate the numerical algorithm used in the UCLA-LES code. As a reference, we use the accurate numerical solution obtained with a DNS code and already employed by some of the authors in previous work [8]. We emphasize again that the focus in this work is the analysis of low-order numerical effects, and not the subgrid-scale models intrinsic to LES (cf. [13]). For this purpose, the UCLA-LES code is modified so that it solves Eqs. (3)–(6), exactly the same equations solved with the reference DNS code. Practically, this is achieved by substituting the subgrid-viscosity in the UCLA-LES code with a constant molecular viscosity ν. The UCLA-LES model discretizes the system of Eqs. (4), (5), and (7) in a finite-volume manner on a static Arakawa-C grid (regular Cartesian), where velocities are staggered half a grid point from the rest of the (cell-centred) variables (χ, p). The grid spacing is uniform in both the horizontal and the vertical direction. Second-order spatial discretization is used, scalar fluxes are limited using the monotonized central flux limiter. Study of low-order numerical effects in the two-dimensional 245 The flow is advanced in time using a third-order Runge-Kutta scheme. Throughout all simulations, Courant numbers of CFLadv =max[u(x), w(x)]t x≤0.5(9) and CFLdiff =tν (x)2≤0.5/4(10) were maintained. The UCLA-LES requires periodic boundary conditions in the lateral as it uses fast Fourier transform in the horizontal direction to transform the Poisson Eq. (7) to a second-order ordinary differential equation (ODE). The ODE is solved using a tridiagonal solver. The reference DNS uses finite differences and a regular Cartesian mesh. Spectral-like sixth-order compact Padé schemes are used to discretize the firstand second-order spatial derivatives (closed at the top and bottom with third-order biased formulae, which leads to a global fourth-order accuracy in space). A low-storage 5-stage fourth-order Runge-Kutta scheme is employed for the time advancement. This scheme has a stability region larger than that of the third-order counterpart, and the Courant numbers used in the DNS are CFLadv <1.2 and CFLdiff <1.2/4. Similarly to the UCLA-LES, the solenoidal constraint is imposed using Fourier decomposition along the periodic horizontal planes in order to reduce it to a set of one-dimensional second-order equations which are solved using those same compact schemes. The boundary conditions imposed at the top and the bottom are the same as those in the UCLA-LES code. Further details of the numerical algorithm can be found in Mellado et al. [10]. The major difference between the DNS and the UCLA-LES code is the resolving efficiency. This resolving efficiency is normally quantified in terms of the number of points per wavelength (PPW) required to maintain errors in the corresponding transfer function below a specified level (or, equivalently, a specified error in the dispersion velocity of the linear advection problem). In these terms, an error of 1 % requires 4 PPW in the case of the implicit compact scheme used in the DNS, whereas the second-order explicit scheme requires 25 PPW [4,6] (this difference is even higher if the common reference error of 0.1 % is retained, for which the previous finite difference schemes require 6 PPW and 100 PPW, respectively). However, these are theoretical values based on linear analysis of the algorithm: in this work, we ascertain the relative performance of these two numerical schemes in an applied case very relevant to the physics of the STBL, namely the stratocumulus top. 3.2 Analysis of the flow The qualitative evolution of the flow is illustrated in Fig. 3through successive snapshots that visualize the negative buoyancy field. The data shown were obtained from the UCLA-LES run on the reference grid with 512 ×1,024 cells. The leftmost panel corresponds to the sinusoidal initial condition, as previously discussed in Sect. 2.3. The stratified system is initially at rest. In the vicinity of the inversion, net vertical accelerations arise from horizontal buoyancy differences. A persistent oscillation of the interface sets in that, at high Gr, is only slowly damped by molecular viscosity. At the same time, a spike or finger starts to form of the thin negatively buoyant sheet at the lowest point of the oscillation. The spike grows downwards and later rolls up to a vortex pair. Molecular diffusion mixes fluid of the falling finger with the ambient moist fluid and, thus, spreads buoyancy reversal. According to the buoyancy mixing function (cf. Fig. 1), moist fluid (χ=0) that is mixed with negatively buoyant fluid is negatively buoyant itself. This process continues as the pulsating motion of the inversion pumps more negatively buoyant fluid downwards. The instability process is driven by the transformation of potential energy (here subsumed within the enthalpy of vaporization) to kinetic energy. As the unsaturated fluid in the upper layer mixes with the saturated fluid in the lower layer, more negatively buoyant fluid is produced due to evaporative cooling, and part of it is ejected downwards. This gradually leads to a vertical growth of both, the size of the saturated bottom layer and the depth of the mixing region. Mellado et al. [8] quantify the instability process by means of two characteristic lengths that measure the thickness of the mixing layer to the top and to the bottom. These are the penetration length hbof the falling finger and the perturbation thickness htthat measures the upwards growth of the moist layer. Both are defined with the help of the horizontally averaged mixture fraction profile ¯χ(z)asshowninFig.4. Taking the mean initial height of the χ=0.5 isoline, ¯z0, as a reference, both lengths are defined as distances to the vertical 246 E. Dietze et al. -2 -1 0 1 2 0 0.5 1 (z - – z0)/a – χ ht/a hb/a Fig. 4 Definition of the penetration length hband the upper perturbation thickness htusing the horizontally averaged mixture fraction profile ¯χ.Thedashed lines indicate the vertical position where ¯χdeparts from a given threshold (10−3) of the upper and lower state, respectively 0.0 1.0 2.0 3.0 4.0 5.0 6.0 0 5 10 15 20 hb/a t(b1/λ)1/2 DNS, Mellado LES (1/8 res.) LES (1/4 res.) LES (1/2 res.) LES (1x res.) LES (2x res.) LES (4x res.) 0.0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 Δhb/hb,DNS t(b1/λ)1/2 1/8 res. 1/4 res. 1/2 res. 1x res. 2x res. 4x res. Fig. 5 History of the penetration length (left) and the relative error with respect to the reference DNS (right). Grid resolutions are given with respect to the reference grid with 512 ×1,024 grid points position where the mean mixture fraction departs from a given threshold of the value of the respective layer. For hb, the first position looking from the bottom, for htthe first position looking from the top are detected, respectively. Because mixing in the lower layer is limited to small fractions of ¯χbelow 0.1, a threshold of 10−3is chosen. The temporal evolution of these lengths is presented in Figs. 5and 6, respectively, where the solid lines indicate the results from the DNS by Mellado et al. [8]. The evolution of the penetration length hbobserved in Fig. 5exhibits a superposition of the stable oscillatory mode of the inversion and the unstable mode of the falling finger. If buoyancy reversal was not present (cf. dashed graph in Fig. 1), only the stable mode would be visible. The evolution of the upper perturbation thickness htalso has the stable mode; it is superimposed by a quasi-linear growth that results from molecular mixing of χat the top of the inversion and the baroclinic production of vorticity. The transport of mixed fluid away from the inversion due to buoyancy reversal steepens local scalar gradients which in return enhances molecular mixing. In the following, we use these two lengths to quantify numerical errors and the convergence of the UCLA-LES code. 3.3 Convergence behaviour In order to study the convergence behaviour of the UCLA-LES code, a grid convergence study was carried out. The reference DNS employs a grid with 512×1,024 grid points; this will be referred to as reference resolution or reference mesh in the following. Here, simulations were carried out on six different meshes, ranging from 1/8th to 4 times the reference resolution (see Table 1). Study of low-order numerical effects in the two-dimensional 247 0.4 0.6 0.8 1.0 1.2 0 5 10 15 20 ht/a t(b1/λ)1/2 DNS, Mellado LES (1/8 res.) LES (1/4 res.) LES (1/2 res.) LES (1x res.) LES (2x res.) LES (4x res.) 0.0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 Δht/ht,DNS t(b1/λ)1/2 1/8 res. 1/4 res. 1/2 res. 1x res. 2x res. 4x res. Fig. 6 History of the upper perturbation thickness (left) and the relative error with respect to the reference DNS (right). Straight lines (left) are least squares regressions of the growth rate from the 2nd to the 9th local maximum. Grid resolutions are given with respect to the reference grid with 512 ×1,024 grid points Table 1 Order of convergence pof the UCLA-LES according to the L2and maximum norm L∞ Mesh. 1/8 1/4 1/2 [512 ×1,024]2×4× p[L2]–1.9991 1.9975 1.9987 1.9976 – p[L∞]–1.9948 1.9951 1.9874 1.9508 – The factors in the mesh row refer to the number of mesh cells, for example, 1/8 refers to an eight times coarser mesh The UCLA-LES code, used as explained in Sect. 3.1, has an anticipated second-order convergence in space. To confirm this, simulations on all meshes have been performed with a common time increment, and the order of convergence pin terms of the L2norm and the maximum norm L∞of the mixture fraction field have been determined. Based on the solution of χon three meshes resolved by x=z=h,2h,and4h, respectively, the norms L2(h)=1 NN i=1[χ2h(xi)−χh(xi)]2 and L∞(h)=max(|χ2h(xi)−χh(xi)|) are computed, where xi=(xi,zi)is the position of the ith cell-centred point of the mesh. Linear interpolation on the staggered mesh is used to interpolate values from the finer mesh to the coarser mesh. The order of convergence is computed according to p2,∞=ln L2,∞(2h) L2,∞(h)/ln(2). This procedure is repeated for all subsets of three neighbouring meshes. The results are presented in Table 1 and confirm the anticipated second-order convergence in both, the maximum and L2norm. In the following, we want to quantify the absolute errors of the solution. We do this by comparing the evolution of hband htin the DNS and LES code. Besides the reference DNS solution, Figs. 5and 6show the results obtained from the LES on the six meshes. For both hband ht, a clear trend is obvious: with coarsening the mesh spacing the growth of both quantities is increased. This can mainly be attributed to increased numerical mixing, yielding more generation of negatively buoyant fluid that enhances the convective mixing process in the cloud layer. The temporal evolution of the relative errors is explicitly showed in the second panel of Figs. 5 and 6. On the coarsest grid, relative errors of hbare up to 30 % and decrease to 5 % at reference resolution andto2.4%onthefinestgrid.Relativeerrorsofhttend to be greater on the coarsest grid but vanish more rapidly as the grid is refined. On the finest grid, that is, four times the number of grid points in each direction than in the reference DNS, the relative error remains below 1 % throughout the simulation.