Study of low-order numerical effects in the two-dimensional cloud-top mixing layer
Juan Pedro Mellado
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. 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 . 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 K km1) giving rise to the terminology inversion layer, as the temperature gradient is inverted relative to
the tropospheric norm. Already the mixed-layer theory by Lilly  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  that relates minimum and maximum buoyancy bs and b1,
respectively. This parameter is typically small. For instance, in the first research flight (RF01) of the DYCOMS
II field study , 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  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  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  and the one-dimensional turbulence model  but also direct numerical simulations
(DNS), based on the Boussinesq approximation and a continuous two-fluid formulation . Mellado et al. 
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
 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 = /we forms and travels upwards at a constant mean speed we (|bs c2)1/3, where is the
molecular diffusivity and c is 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 . 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.
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.  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)  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
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, hb and 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.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
qt and specific enthalpy h transition 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 . Assuming equal thermal and mass diffusivities, the thermodynamic coordinates h and
qt obey 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 scalarthe 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 h and qt as
h(x, t ) h0
qt (x, t ) qt,0
Herein, x = (x , z)T denotes the location and t is 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,
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. 
= 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
g denoting 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 = 00s g. Shy and Breidenthal  introduce the buoyancy reversal
b(x, t ) =
to describe the strength of the reversal, relating buoyancy of the heaviest mixture to the buoyancy jump
b1 = 001 g across 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
introduced by Mellado et al.  is adopted. Figure 1 shows 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 s was introduced to smooth the discontinuity between the piecewise
linear sections. Obviously, the smoothing damps the maximum negative buoyancy bs and, thus, the effective
reversal parameter D. Previous work  shows that a smoothing parameter of s = s /16 leads to a
sufficient independence of the solution from changes in s . This value of s is adopted here for consistency in the
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 p and the mixture fraction .
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 = 0.
The LES code, whose native formulation solves Ogura-Phillips  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), k is the vertical unit vector, and p is a modified pressure, normalized
by the density.
The Poisson equation for the dynamic pressure
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) = 2
z z0(x )
Herein, z0 denotes 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
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 ) = 2 cos
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, a measures 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 . 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.  introduce the reference Grashof number
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 . For the reference simulation, the kinematic viscosity = 1.5 105 m2/s and Gr = 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
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.1 The numerical methods
We use the cloud-top problem defined in Sect. 2 as 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 . 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. ).
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.
The flow is advanced in time using a third-order Runge-Kutta scheme. Throughout all simulations, Courant
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 first- and 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. .
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. 3 through 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.  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 hb of the falling
finger and the perturbation thickness ht that measures the upwards growth of the moist layer. Both are defined
with the help of the horizontally averaged mixture fraction profile (z) as shown in Fig. 4. Taking the mean
initial height of the = 0.5 isoline, z0, as a reference, both lengths are defined as distances to the vertical
Fig. 4 Definition of the penetration length hb and the upper perturbation thickness ht using the horizontally averaged mixture
fraction profile . The dashed lines indicate the vertical position where departs from a given threshold (103) of the upper and
lower state, respectively
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 ht the 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
103 is chosen.
The temporal evolution of these lengths is presented in Figs. 5 and 6, respectively, where the solid lines
indicate the results from the DNS by Mellado et al. . The evolution of the penetration length hb observed
in Fig. 5 exhibits 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 ht also 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).
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 p of the UCLA-LES according to the L2 and maximum norm L
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 p in terms of the L 2 norm 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, and 4h,
respectively, the norms
L 2(h) =
are computed, where xi = (xi , zi ) is the position of the i th 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
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 L 2 norm.
In the following, we want to quantify the absolute errors of the solution. We do this by comparing the
evolution of hb and ht in the DNS and LES code. Besides the reference DNS solution, Figs. 5 and 6 show the results
obtained from the LES on the six meshes. For both hb and 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 hb are up to 30 % and decrease to 5 % at reference resolution
and to 2.4 % on the finest grid. Relative errors of ht tend 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.
Table 2 Relative errors of the ht growth rate ELES in the low-order simulations with respect to the DNS solution
Fig. 7 Amount of fluid entrained in the section of the domain below the lowest point of the initial wave of the = 0.5 isoline.
It is measured by the integral I of the mean mixture fraction profile normalized by as in Eq. (8)
The mean vertical growth of the upper perturbation thickness ht is an effect that resembles the cloud-top
entrainment in real stratocumulus. We quantify this mean growth E with help of a linear regression based on
the data points after the first period until the ninth period of each individual data set. The latter corresponds to
the last ht maximum in the DNS reference data. The regression lines for the DNS solution and the finest and
coarsest LES solution are shown in Fig. 6. Relative errors with respect to the DNS solution are presented in
Table 2. While at four times the reference resolution the error in the mean vertical growth of the upper
perturbation thickness reduces to 1.4 %, it increases to 89 % on the coarsest grid with 1/8th the reference resolution.
The latter corresponds to a mesh spacing of about 1 cm, still much finer than in current high-resolution LES
of stratocumulus, where typically the vertical resolution is between 1 and 5 m at the cloud top [19, 24]. This
result shows clearly that state-of-the-art LES does not resolve the sharp gradients at the cloud top and that this
under-resolution has important effects on the thermodynamics of the system and helps explain the difficulties
encountered in understanding local processes of the cloud top .
In Fig. 7, we show the total amount of fluid entrained from the upper layer into the section of the domain
below the lowest point of the initial wave. This portion of fluid is simply measured by the (vertical) integral
I of the mean mixture fraction profile in that subdomain and that portion of fluid monotonously increases
when the grid is coarsened. This indicates that the more rapid increase in hb on the coarser grids is, in fact,
associated with a larger amount of fluid from the upper layer carried downwards by the falling finger. It is not
solely attributed to the stronger spurious (numerical) spread of the border of the falling finger, but also to a
modification of the global dynamics of the system. Here, numerical errors act as to increase the entrainment
of drier fluid into the cloudy layer.
3.4 Computational efficiency
Both the DNS and UCLA-LES algorithms were originally implemented to solve three-dimensional problems.
For the work discussed so far, the codes have been run on a two-dimensional grid in order to reduce
computational cost and focus on some details of the cloud top in a very controlled manner. However, that causes
computational overheads that depend on the respective methods and their implementations. To enable for a
fair comparison in terms of computational efficiency, three-dimensional runs were carried out. We extended
DNS results are based on a low-storage fourth-order Runge-Kutta scheme and UCLA-LES results on a third-order Runge-Kutta
the 128 256 grid (1/4th of the reference resolution) in the lateral to 128 128 256. Both codes were
compiled with the same compiler and run on one core of an AMD Opteron model 2384, 2.7 GHz of a Sun
X2200 M2 machine. In Table 3, we compare the computational time per iteration, T1, and the computational
time T2 needed per characteristic time of the system /b1 (averages over 100 iterations were considered in
order to eliminate I/O times and other overheads).
The DNS is more expensive per iteration, as expected because of the fourth-order Runge-Kutta scheme
(which uses 5 stages because of the two-level storage property) and an implicit spatial discretization. When the
DNS code is run with a third-order Runge-Kutta scheme, similarly to the LES code, the time per iteration is
22 seconds. This is about 40 % more expensive than the low-order code; this 40 % difference is due purely to
spatial discretization and accounts for the overhead due to the solution of the linear systems being associated
with the implicit compact schemes, in contrast to the explicit character of the second-order central schemes.
However, it is arguably more interesting to consider the computational cost per characteristic physical time
of the system, that is, how much it costs to advance the transport equations during a given time interval. This
is the second parameter shown in Table 3, the time T2. The difference with the measurement per iteration
T1 discussed in the previous paragraph is due to the fact that the low-storage fourth-order Runge-Kutta has a
stability region significantly larger than that of the third-order Runge-Kutta and allows a larger CFL number.
Comparisons of the values of T2 between the DNS and the UCLA-LES codes show that the DNS is still slower,
but only by 10 %.
This last result is very important. As we have shown in Sect. 3.3, using the low-order numerical method
of the UCLA-LES code, higher resolution is needed to converge to the reference solution. For the particular
problem we are studying, the buoyancy reversal instability at the cloud top, using the low-order algorithm leads
to errors of the order of 20 % and is only about 10 % faster. Let us imagine we desire to reduce errors using
the UCLA-LES code to about 10 %. Then, according to Table 2, we need to increase the resolution by about
a factor of two. This step already implies a higher computational cost; a reasonable rough estimate would be
a factor of 8 in three-dimensions, or a factor of 4 in two dimensions, assuming a linear relation between the
operation count and the number of grid points and neglecting deviations due to parts of the algorithm which
do not scale linearly with the number of points, like fast Fourier transforms or cache-size management. On top
of that, the time increment t imposed by the explicit scheme, so stability is satisfied, multiplies further this
factor by 2 or 4, according to Eqs. (9) or (10), respectively. In our particular problem, it has been checked that
the latter dominates because the diffusion constraint Eq. (10) becomes the most restrictive and t scales as
( x )2. Altogether, this analysis means about a factor of 10 or more increase in computational cost using a
loworder algorithm, compared to a factor 1.1 increase using the implicit spectral-like compact finite differences,
to obtain a solution within a 10 % error.
It is appropriate to note that the discussion of the last paragraph refers to a serial implementation of the
algorithm. The advantage of low-order spatial schemes is that they are explicit and parallelizations thereof
are known to scale almost linearly, so that the computational cost remains approximately independent of the
number of processors, whereas the implicit schemes suffer a strong overhead from the communication among a
relatively large number of processors. Hence, there is a trade-off and a crossover number of processors beyond
which a low-order scheme is more efficient than a high-order implicit scheme. However, this crossover value
is a relatively high number. For instance, the DNS algorithm used here as a reference is run typically using
1,024 or 2,048 MPI tasks (e.g. simulations presented in ) in the German Climate Computer Centre. In
those architectures, the time spent in communication among processors varies between 40 and 60 %, meaning
an overhead of the order of a factor of 2 with respect to the reference computational cost. The same behaviour
has been observed in scaling studies of the DNS code up to 8,192 tasks in the Jlich Supercomputer Centre.
This overhead is still significantly smaller than the factor of 10 estimated above for a second-order scheme
used in the UCLA-LES code to obtain solutions within 10 %. Last but not least, we need to consider as well
the concomitant additional memory requirements using low-order schemes because of the larger number of
points (this can actually be the most critical aspect in large-scale three-dimensional simulations).
Hence, as a summary, high-order schemes, in spite of the higher complexity in implementation, might be
attractive for LES of stratocumulus, and other high Reynolds number flows if small-scale processes might
play a relevant role or simply to better represent the smooth portions of the flow, such as waves, that might
exist at some multiple of the grid scale. This has implications for the current trend to increase resolution in
LES models and the application of adaptive solvers to such problems. These conclusions hold also for parallel
implementation, at least up to the order of 10,000 processors, which current architectures provide.
We presented in detail the two-dimensional evaporatively driven cloud-top mixing layer as a test case for the
study of low-order numerical effects. We did that comparing the UCLA-LES code, which uses an explicit
second-order scheme, against a DNS solution, based on a spectral-like compact sixth-order accurate method. The
subgrid-scale turbulence model of the UCLA-LES code was switched off, and a constant molecular viscosity
was used instead so that only numerical effects were addressed in this work. We presented a grid convergence
study where we investigate the sensitivity of the solution to numerical errors and the resolution requirements
for the UCLA-LES algorithm to achieve convergence.
Analysis of the solutions on the various grids showed that the LES numerical algorithm is self-convergent
meeting the anticipated effective second-order rate of convergence. In order to quantify numerical errors, we
compared the evolution of the penetration length hb of the downdraught and the upper perturbation thickness
ht between the LES and DNS solutions. The coarsest grid had a vertical spacing of approximately 1 cm which
is at least two orders of magnitude finer than in typical LES. For these quantities, 30 and 40 % relative errors,
respectively, were measured on the coarsest mesh; almost 100 % relative errors in the entrainment rate. Hence,
we have shown that very high-resolution or high-order schemes are needed to capture the small-scale molecular
mixing. If these small-scale effects are not accurately represented, important large-scale characteristics, such
as the vertical growth of the inversion height, resembling cloud-top entrainment, may contain order one defects
caused by numerical artefacts.
The analysis of the computational efficiency showed that the cloud-top problem is solved more efficiently
(in terms of total CPU time) by the high-order scheme used in the reference DNS than by the low-order scheme
of the UCLA-LES code. Based on this finding, we conclude that increasing the order of numerical methods
might be better suited than further increasing the resolution of low-order methods to improve the LES of flows
that are not smooth on the grid scale (such as LES of stratocumulus). The present study suggests that the larger
computational cost per integration step is by far compensated by a smaller number of grid points, a larger
possible temporal stride, and reduced magnitude of numerical errors. The difference is so remarkable that the
same conclusion holds in a parallel implementation of the algorithms, at least up to the order of 1-10 thousand
processors, despite the communication overhead associated with the implicit schemes.
Acknowledgments Financial support for this study was provided by the Deutsche Forschungsgemeinschaft within the SPP 1276
MetStrm and the Max Planck Society for the Advancement of Science. Computational resources have been provided by the
German Climate Computing Centre, Hamburg.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use,
distribution, and reproduction in any medium, provided the original author(s) and the source are credited.