Geoscientific Model Development Formulation of the Dutch Atmospheric Large-Eddy Simulation ( DALES ) and overview of its applications

The current version of the Dutch Atmospheric Large-Eddy Simulation (DALES) is presented. DALES is a large-eddy simulation code designed for studies of the physics of the atmospheric boundary layer, including convective and stable boundary layers as well as cloudy boundary layers. In addition, DALES can be used for studies of more specific cases, such as flow over sloping or heterogeneous terrain, and dispersion of inert and chemically active species. This paper contains an extensive description of the physical and numerical formulation of the code, and gives an overview of its applications and accomplishments in recent years.


Introduction
Modern atmospheric research relies on a spectrum of observational and modeling tools.Among the numerical tools that are commonly used for the most detailed studies of atmospheric flow, an important spot is taken by Large-Eddy Simulations (LES).This type of modeling is widely used for atmospheric boundary layer (ABL) studies and provides, in combination with observations, the basis for many cloud and boundary layer parameterizations in models on the other side of the spectrum, such as General Circulation Models.
Correspondence to: T. Heus (thijs.heus@zmaw.de) The principle of LES is to resolve the turbulent scales larger than a certain filter width, and to parameterize the smaller, less energetic scales.This filter width, in practical applications often a function of the grid size of the LES, and ranges typically between 1 m for stably stratified boundary layers, to 50 m for simulations of the cloud-topped ABL.In such a typical LES setup, up to 90% of the turbulence energy resides in the resolved scales.For applications of LES like the ones presented in this paper, LES has the advantage over larger-scale models that it relies less on parameterizations.In comparison with observational studies, LES has the advantage of providing a complete data set, in terms of time and space, and in terms of variables that can be diagnosed.Especially the combined use of LES and observations is a popular methodology in process studies of the ABL.In comparison with the yet finer Direct Numerical Simulations (DNS) that aim to resolve all turbulence scales, LES has the advantage of being able to cover larger domains than a few meters.
LES modeling of the ABL started in the late sixties (e.g., Lilly, 1967;Deardorff, 1972); cloudy boundary layers were first simulated by Sommeria (1976).From Nieuwstadt and Brost (1986) onward, several cycles of intercomparison studies compare state-of-the-art LES codes with observational studies and with each other.The aim of these studies was not so much to determine which LES code performs best in which situation, but more to determine the strengths and weaknesses of LES.Two particularly active cycles are organized under the umbrella of the Global Energy Published by Copernicus Publications on behalf of the European Geosciences Union.
The Dutch Atmospheric Large-Eddy Simulation (DALES) has joined virtually all of the intercomparisons mentioned in the previous paragraph.Beyond these intercomparison studies, that discuss convective, stable and cloud-topped boundary layers, DALES has been used on a wide range of topics, such as for studies of shear driven flow, of heterogeneous surfaces, of dispersion and of turbulent reacting flows in the ABL, and of flow over sloped terrain.Whenever appropriate, results from DALES have been compared to observational data to provide additional validation for the less standard use cases.In a recent effort, DALES is being used in the KNMI Parameterization Testbed (Neggers et al., 2010), which allows for a day-by-day comparison between observational data, LES, and large-scale model results.As such, DALES is one of the most all-round tested available LES codes for studies of the ABL.In this paper, we aim to describe and validate DALES 3.2, the current version of DALES.
In the remainder of this paper, we first give a thorough description of the code in Sect. 2. In Sect.3, an overview of studies conducted with DALES is given, both as a validation of the code as well as an overview of the capabilities of an LES like DALES.In Sect.4, an outlook is given on future studies that are planned to be done with DALES, as well as an outlook on future improvements.

Generalities
DALES is rooted in the LES code of Nieuwstadt and Brost (1986).Cuijpers and Duynkerke (1993) first used DALES for moist convection, and provided a general description of an older version of DALES.Large parts of the code have been changed ever since and contributions of many people over a number of years have resulted in the current version 3.2 of DALES.Currently, DALES is maintained by researchers from Delft University, the Royal Netherlands Meteorological Institute (KNMI), Wageningen University and the Max Planck Institute for Meteorology.
Notable changes in comparison with the version that has been described by Cuijpers and Duynkerke (1993) include: different time integration and advection schemes, revised subfilter-scale, surface and radiation schemes, addition of a cloud-microphysical scheme, capabilities for chemical reactive scalar transport and for Lagrangian particle dispersion, for flow over heterogeneous and for flow over sloping terrain.These revisions in DALES result in faster simulations and higher stability, and in an easier and more extendable user interface.Due to the modular setup of the code, newly written code for specific applications of DALES can easily improve the code as a whole.This makes DALES suitable as a community model; besides the actively developing core users, the code is currently used in several other institutes across the world.DALES 3.2 is released under the GPLv3 license.It is available at dales.ablresearch.org.Documentation is also available there.Although the code is completely free to use, to modify and to redistribute, it is regarded courtesy to share bug fixes and extensions that can be of general interest, and to keep in contact with the core developers.Given the experimental character of the code, it is also appreciated to discuss co-authorship in case of publications coming out of research conducted with DALES.
To improve compatibility and portability of the code, we make an effort to stay as close as possible to standard Fortran 95.To create makefiles and compile the code, Kitware's Cmake (www.cmake.org) is being used.This system is installed on most modern systems (and is usually installable with user permissions if that is not the case).Cmake facilitates flexible handling of compiler options, various platforms and library locations.It also automatically keeps track which source code needs to be (re)build.For the communication between multiple processes, DALES relies on the Message Passing Interface (MPI).For purposes of input and output, the Network Common Data Form (NetCDF) version 3 or higher is an optional dependency.Code for Fourier transformations is part of the DALES package, leaving the code as portable as possible.To the best knowledge of the authors, DALES runs on all common combinations of platform architecture, compiler, and MPI implementation.Currently, an effort is being made to port DALES to Compute Unified Device Architecture (CUDA), to be able to run simulations on graphical processors as a fast and cost efficient solution.
The prognostic variables of DALES are the three velocity components u i (i = 1,2,3), the liquid water potential temperature θ l , the total water specific humidity q t , the rain water specific humidity q r , the rain droplet number concentration N r , and up to 100 passive or reactive scalars.The subfilterscale turbulence kinetic energy (SFS-TKE, e) is an additional prognostic variable, and is being used in the parameterization of the sub-filter scale dynamics.To decrease simulation time, only calculations of u i , e, and θ l are obligatory; all the additional scalars need not to be calculated when these variables are not used.
Given that ice is not currently implemented in the model, the total water specific humidity is defined as the sum of the water vapor specific humidity q v and the cloud liquid water specific humidity q c : q t = q v + q c . (1) Note that this definition of q t does not include the rain water specific humidity q r .Any conversion between rain water on the one hand, and cloud water or water vapor on the other hand, will therefore enter the equations for q t and for θ l as an addition source term.We use the close approximation explained by Emanuel (1994): with θ being the potential temperature (related to the absolute temperature T following T = θ ), L = 2.5 × 10 6 J kg −1 the latent heat of vaporization, c pd = 1004 J kg −1 K −1 the heat capacity of dry air, and being the exner function: in which R d = 287.0J kg −1 K −1 is the gas constant for dry air and p 0 =10 5 Pa is a reference pressure.
In the absence of precipitation and other explicit sources, θ l and q t are conserved variables.The liquid water virtual potential temperature θ v is in good approximation defined as (Emanuel, 1994): with R v = 461.5 J kg −1 K −1 being the gas constant for water vapor.The most important thermodynamical constants that are used throughout this paper are summarized in Table 1.DALES is run on an Arakawa C-grid (see Fig. 2).The pressure, the SFS-TKE, and the scalars are defined at grid cell center, the three velocity components are defined at the West side, the South side, and the bottom side of the grid cell, respectively.
Hereafter, quantities that are averaged over the LES filter width are denoted with a tilde •, time averages with a overbar • , and averages over the two horizontal directions of the domain with angular brackets • (slab average).The prognosed scalars can often be treated in an identical manner as the generic scalar field ϕ∈{θ l ,q t ,q r ,N r ,s n }.Primes denote the subfilter-scale fluctuations with respect to the filtered value.Double primes indicate local deviations from the horizontal slab average.To remain consistent with notational conventions used in literature and also in the source code of DALES, some symbols can have different meaning between different subsections.In such cases, the immediate context should always make it clear what each symbol stands for in a particular section.Vertical velocities and fluxes are in general positive when directed upward; only the radiative and

R v
Gas constant for water vapor 461.5 J kg −1 K −1 R d Gas constant for dry air 287.0 J kg −1 K −1 L Latent heat release for vaporization 2.5 × 10 6 J kg −1 c pd Heat capacity for dry air 1004 J kg −1 K −1 sedimentation fluxes are positive when pointing downward, following conventions.
In the following sections, different components of the code are described one by one.Sections 2.2-2.7 describe the physical and numerical components that are necessary to conduct a minimal experiment with DALES.After that, Sects.2.8-2.12describe various forcings and source terms that extend the core of DALES for use in more specific applications.Finally, Sect.2.13 describes the most relevant statistical routines in DALES.

The governing equations
DALES assumes the Boussinesq approximation, with the reference state θ 0 ,ρ 0 ,p 0 equal to the surface values of liquid water potential temperature, density and pressure, respectively.For an extended treatment see for example Wyngaard (2004).
Within the Boussinesq approximation the equations of motion, after application of the LES filter, are given by where the tildes denote the filtered mean variables.Molecular transport terms have been neglected.The z-direction (x 3 ) is taken to be normal to the surface, π = p ρ 0 + 2 3 e is the modified pressure, δ ij the Kronecker delta, and F i represents other forcings, including large scale forcings and the Coriolis acceleration where is the Earth's angular velocity.Source terms for scalar ϕ are denoted by S ϕ , and may include of microphysical (S mcr ), radiative (S rad ), chemical (S chem ), large-scale (S ls ), and relaxation (S rel ) terms.The subfilter-scale (SFS), or residual, scalar fluxes are denoted by R u j ,ϕ ≡ u j ϕ− u j ϕ, i.e., the contribution to the resolved motion from all scales below the LES filter width.The deviatoric part of the subgrid momentum flux: A schematic overview of how the different processes affect the different variables is given in Fig. 1.

Subfilter-scale model
In DALES, the SFS fluxes are modeled through an eddy diffusivity as: and where e = 1 2 ( u i u i − u i u i ) is the subfilter-scale turbulence kinetic energy (SFS-TKE) and K m and K h are the eddy viscosity/diffusivity coefficients.In DALES, these eddy diffusivity Fig. 2. The Arakawa C-grid as used in DALES.Pressure, SFS-TKE and the scalars are defined at cell-center, the 3 velocity components at the face of the cell.The level of cell center is called the full level (denoted with an "f"); the level where w is located is called the half level (an "h").The (variable) vertical grid spacing z is defined centered around the belonging level.The grid spacing in the horizontal directions ( x and y) are constant over the entire domain.
coefficients are modeled in two ways: either as a function of the SFS-TKE e (Deardorff, 1973) (which is the default), or using Smagorinsky closure (Smagorinsky, 1963).

SFS-TKE model
Following Deardorff (1980), the prognostic equation for e is adopted in the form: with ε the SFS-TKE dissipation rate.The first right-handside term is solved, and the second term (the production of SFS-TKE by shear) can be calculated with Eq. ( 11).The other right-hand-side terms need to be parameterized to close the equation.Following Deardorff (1980), we express the third term, the SFS-TKE production due to buoyancy, as: with coefficients A and B depending on the local thermodynamic state (dry or moist): Geosci.Model Dev., 3, 415-444, 2010 www.geosci-model-dev.net/3/415/2010/ where q s is the saturation specific humidity at the given temperature.At a cloud interface, it is a matter of choice whether to use the dry or the moist coefficients in calculation of the buoyancy production.Especially in situations where the properties of the cloud deck are around the buoyancy reversal criterion, this choice proves to be critical (Randall, 1980;Bretherton et al., 2004;de Roode, 2007).To determine whether a parcel that contains a mixture of saturated and unsaturated air is saturated itself, we calculate the amount of unsaturated air that is needed to evaporate all the liquid water in a mixed air parcel.In particular the saturation mixing ratio χ sat defines the ratio of cloudy to total air mass for a just saturated mixed air parcel (Stevens, 2002): where θ l = θ l (z + z)− θ l (z − z) and q t = q t (z + z)− q t (z − z) are the differences over the cloud interface.If turbulent mixing occurs, it is assumed that at level z k the mass mixing fraction is If χ < χ sat , the mixed parcel will be saturated and consequently the coefficients for saturated air (Eq.15) will be used.
The fourth and fifth term in Eq. ( 12) are together parameterized as To model the dissipation rate ε, we again follow (Deardorff, 1980): with α = 1.5 the Kolmogorov constant and c f λ the filter width.
The eddy diffusivity for heat and scalars is modeled similarly as K h = c h λe 1/2 , and for the dissipation ε we write: Still following Deardorff (1980), the SFS parameters are depending on the stability of the flow: 1.5 2.5 0.19 0.51 0.12 1 2 0.76 with denoting the Brunt-Väisälä frequency, and c N = 0.76.Now all parameters of the subfilter-scale parameterization of DALES are defined; they are summarized in Table 2.
Substituting the closure relations and parameters into Eq.( 12) gives the following prognostic equation for e 1/2 , which is implemented in DALES: which closes the system.

Smagorinsky SFS modeling
The Smagorinsky model that is implemented in DALES is as follows (Wyngaard, 2004): the Smagorinsky constant and the strain tensor, respectively.Pr = K m /K h = 0.33 is the Prandtl number, equivalent to c m T. Heus et al.: The Dutch Atmospheric LES but does not resolve the flow up to the surface-roughness scale.The surface fluxes enter the domain at subfilter-scale, since by definition the resolved fluctuations in the vertical velocity at the surface are equal to zero.In the remainder of this section we define an arbitrary surface flux of variable φ as F s,φ = wφ − w φ.
We followed the common way of parameterizing turbulent fluxes in atmospheric models by applying the transfer laws of Louis (1979).In DALES we assume that the first model level is in the atmospheric surface layer.We apply Monin-Obukhov similarity theory for the computation of the spatially averaged fluxes F s,φ and gradients at the bottom boundary of the model.
The procedure for determining the bottom boundary conditions starts with the evaluation of the Obukhov length.This value is approximated using a Newton-Rhapson method for solving the implicit equation that relates the bulk Richardson number to the Obukhov length (see Eq. 28). with and where Ri B is the averaged bulk Richardson number of the layer between the surface and the first full level z 1 , L is the Obukhov length, z 0m and z 0h are the roughness lengths for momentum and heat, H and M are the integrated stability functions as provided by Beljaars (1991) for the stable atmosphere and Wilson (2001) for the unstable atmosphere, θ v0 is the spatially averaged filtered surface virtual potential temperature, θ v1 is the spatially averaged filtered virtual potential temperature at the first model level, U 1 is the magnitude of the horizontal wind vector at the first model level, defined as , κ is the Von Karman constant and w θ v0 is the horizontally averaged surface virtual temperature flux.Subsequently, the calculated Obukhov length is used in the computation of the slab averaged friction velocity u * 0 and scalar scales ϕ * 0 = − F s,φ u * 0 , based on the scaling arguments of Businger et al. (1971); Yaglom (1977).Now, we can calculate the drag coefficients C M and C ϕ : Although all locations in the horizontal use the same drag coefficient, we calculate local fluxes and gradients that average to the values computed in our evaluation of the Obukhov length.The subfilter-scale momentum fluxes are calculated by decomposing u 2 * 0 along the two components of the horizontal wind vector (Eqs.33, 34), whereas Eq. ( 35) gives the scalar flux.This results in For land surfaces where moisture is not freely available, such as a vegetated land surface or a bare soil, an additional step has to be made before the similarity relation as in Eq. ( 35) can be applied to the specific humidity.Here, we define the aerodynamic resistance r a as C ϕ U 1 −1 and introduce the surface resistance r s that takes into account the limited water supply at the land surface.The value for r s is either prescribed or calculated using the Jarvis-Stewart model (Jarvis, 1976), where the correction functions for radiation, soil moisture and vapor pressure deficit are taken from the ECMWF Integrated Forcast System (cycle 31R1), and the correction function for temperature from Noilhan and Planton (1989).A urface value can be computed: Note that the drag coefficients and resistances are based on slab averaged values, to assure that the spatially averaged fluxes and gradients are consistent with Monin-Obukhov similarity theory.In DALES there is also the option available to work with locally computed values.We are aware that this method overpredicts gradients at the first model level (Bou-Zeid, 2005).We, however, use this method for exploratory experiments over heterogeneous land surfaces, because here a universal surface model formulation is still lacking (Bou-Zeid, 2005).

Overview of surface boundary options in DALES
DALES has four options to provide the surface momentum and scalar fluxes and surface scalar values to the model, with different degrees of complexity.
1. Parameterized surface scalar and momentum fluxes, parameterized surfaces values.Here, a Land Surface Model (LSM, see Sect.2.4.3)calculates the surface temperature and the stomatal resistance which enters in the evaporation equation based on the vegetation type that is assigned to the grid cell.The variables u * 0 , L and ϕ * 0 are determined iteratively to get the drag coefficients.This is the method that represents a fully interactive land surface.Combined with the radiation model, this options allows for the simulation of full diurnal cycles, in which both the surface fluxes and the surface temperature are free variables.
2. Parameterized surface scalar and momentum fluxes, prescribed scalar values at the surface.In this option u * 0 , L and ϕ * 0 are solved iteratively to get the drag coefficients.The surface momentum and scalar fluxes are computed using the prescribed scalar values at the surface and the acquired drag coefficients.This option is commonly used as the surface boundary condition for simulations of marine boundary-layers.It is also applied in the simulation of stable boundary layers.For simulations over land, a fixed surface resistance r s can be prescribed.
3. Prescribed surface scalar fluxes, prescribed u * 0 .In this option no iterations are necessary and the scalar surface values ϕ 0 are calculated diagnostically.This is an option that is used in idealized simulations in which the surface drag is preferred to be controlled, thereby neglecting that u * 0 is an internal parameter of the flow.
Here u * 0 and L are calculated iteratively, whereas ϕ * is diagnostically calculated as a function of the prescribed scalar fluxes and the calculated u * 0 .This is the most commonly used option for simulation of daytime convection over land.
Prescribed fluxes or surface values may depend on time; linear interpolation is then performed between the given "anchor" points.
In addition to the previous description which assumed homogeneous surfaces, DALES is also able to simulate heterogeneously forced ABLs.Under such conditions, only the prescribed scalar fluxes boundary conditions are available.Scalar fluxes are defined per grid cell, whereas the momentum flux is dynamically computed.In these conditions, local values of u * 0 , L and ϕ * 0 are used.

Land surface model
DALES has the option to use a land surface model (LSM).The LSM has two components, namely a solver for the surface energy balance and a four layer soil scheme which calculates the soil temperature profile for each grid cell.The following surface energy balance equation is solved: in which C sk is the heat capacity per unit of area of the skin layer (see Duynkerke, 1999), T 0 is the surface temperature, Q * is the net radiation and G is the ground heat flux.If the LSM is used, the surface resistance r s in Eq. ( 36) is calculated using the Jarvis-Stewart parameterization (Jarvis, 1976).
The ground heat flux is parameterized as: in which is a bulk conductivity for the stagnant air in the skin layer (Duynkerke, 1999) depending on the type of surface, and T soil1 is the temperature of the top soil layer.
The soil consists of four layers in which the heat transport is solved using a simple diffusion equation in which both the conductivity and the heat capacity are functions of the properties of the soil material and of the soil moisture content.The temperature at the bottom of the lowest soil layer is prescribed.

Boundary conditions: the sides and top
In comparison with the boundary conditions at the bottom boundary, the boundary conditions at the top and the sides of the domain are relatively straightforward.In the horizontal directions, periodic boundary conditions are applied for all fields.At the top of the domain, we take: Fluctuations of velocity and scalars at the top of the domain (for instance due to gravity waves) are damped out by a sponge layer through additional forcing/source terms (added to the right-hand-side of the transport equations): with t sp a relaxation time scale that goes from t sp 0 =1/(2.75×10−3 ) s≈6 min at the top of the domain to infinity at the bottom of the sponge layer, which is by default a quarter of the number of levels, with a minimum of 15 levels.

Pressure solver
To solve for the modified pressure π, the divergence ∂ ∂x i of Eq. ( 6) is taken.Subsequently, the continuity equation (Eq.5) is applied (both divergence and continuity equation are applied in their discrete form).As a result, the left hand side of the equation is equal to zero.Rearranging the terms leads to a Poisson equation for the modified pressure: Since computations are performed in a domain that is periodic in both horizontal directions, the Poisson equation is solved by applying a Fast Fourier Transform in the lateral directions followed by solving a tri-diagonal linear system in the z-direction using Gaussian elimination.In the latter, the pressure gradients at the upper and lower boundary are set to zero.An inverse Fast Fourier Transform in both lateral directions is applied to the result of the Gaussian elimination to obtain the modified pressure.
T. Heus et al.: The Dutch Atmospheric LES

Numerical scheme
A Cartesian grid is used, with optional grid stretching in the z-direction.For clarity, an equidistant grid is assumed in the discussion of the advection scheme.The grid is staggered in space as an Arakawa C-grid; the pressure, the SFS-TKE and the scalars are defined at x+ 1 2 ( x, y, z), the u is defined at x+ 1 2 (0, y, z), and similar for v and w.The level of cell center is called the full level (denoted with an "f"); the level where w is located is called the half level (an "h").The (variable) vertical grid spacing z is defined centered around the belonging level (see Fig. 2).The grid spacing in the horizontal directions ( x and y) is constant over the entire domain.
To be able to use multiple computational processes, thus decreasing the wall clock time of experiments, DALES 3.2 has been parallelized by dividing the domain in separate stripes in the y-direction.Tests show that this method is computationally efficient as long as the amount of processes is smaller than a quarter of the number of grid points in the y-direction.In the near future, we plan to also divide the domain in the x-direction, leaving narrow columns to be calculated by each process, and ensuring that the maximum number of processes would scale with the total number of grid points in each slab, thus allowing for much larger experiments.
Time integration is performed by a third order Runge-Kutta scheme following Wicker and Skamarock (2002).With f (φ n ) the right-hand side of the appropriate equation of Eqs.(6-7) for variable φ = { u, v, w,e 1/2 , ϕ}, φ n+1 at t + t is calculated in three steps: with the asterisks denoting intermediate time steps.The size of the time step t is determined adaptively, and is limited by both the Courant-Friedrichs-Lewy criterion ( CFL) and the diffusion number d (see Wesseling, 1996).
The numerical stability and accuracy depends on the spatial scheme that is used.Furthermore, additional terms, such as chemical or microphysical source, may require more stringent time stepping.Therefore, the limiting CFL and d numbers can be manually adjusted to further optimize the timestep.By default CFL and d are set well below the stability levels known from the literature of the respective combinations of spatial and temporal integration scheme (see Wicker and Skamarock, 2002).Depending on the desired properties (like high accuracy or monotonicity), several advection schemes are available.With advection in the x-direction discretized as with the convective flux of variable φ through the i − 1 2 plane; the i − 1 2 plane is the plane through the location of velocity u i (i), perpendicular on the direction of velocity u i (i).Since we are using a staggered grid, the velocity is available at i − 1 2 without interpolation (see Fig. 2).Second order central differencing can be used for variables where neither very high accuracy nor strict monotonicity is necessary: and similar for F 2nd i− 1 2 .A higher-order accuracy in the calculation of the advection is reached with a sixth order central differencing scheme (see Wicker and Skamarock, 2002): Starting with this sixth order scheme, a nearly monotonous fifth order scheme can be constructed by adding a dissipative term to For advection of scalars that need to be strictly monotone (for example chemically reacting species) the κ scheme (Hundsdorfer et al., 1995) has been implemented: in case u > 0. Following Hundsdorfer et al. (1995), κ i−1/2 serves as a switch between third order upwind advection in case of small upwind gradients of φ, and a first order upwind scheme in case of stronger gradients.This makes the scheme monotone, but also more dissipative, effectively taking over the role of the SFS-scheme in regions of strong gradients.

Cloud microphysics
The cloud-microphysical scheme implemented in DALES is a bulk scheme for precipitating liquid-phase clouds.The hydrometeor spectrum is divided in a cloud droplet and rain drop category.The cloud liquid water specific humidity q c is Fig. 3. Representation of the prognostic thermodynamical variables θ l , q t , the microphysical parameter and variables N c , q c , N r , q r , and the microphysical processes relating these variables.
diagnosed using a classic saturation adjustment (Sommeria and Deardorff, 1977).The cloud droplet number concentration N c (with dimensions m −3 ) is a fixed parameter that can be adjusted according to the degree of pollution of the cloud.Two precipitation schemes have been implemented, both 2moment bulk schemes: rain drop spectra are characterized by the rain drop number concentration N r and the rain water specific humidity q r .The first one is based on Seifert and Beheng (2001, hereafter SB01) (and will be referred to SB01 scheme) two-moment bulk scheme developed for heavy precipitating warm clouds and the second one on Khairoutdinov and Kogan (2000, hereafter KK00), valid only for stratocumulus clouds.
For each prognostic variable modified in microphysics, the source term due to microphysical processes S mcr can be described as effects of autoconversion (au), accretion (ac), rain drop selfcollection (sc), break-up (bu), rain sedimentation (ser), cloud droplet sedimentation (sec), and of rain evaporation (evr): Microphysical tendencies in θ l can be expressed directly as function of q t tendencies: with the exner function based, from here, on the slab averaged pressure p .The prognostic thermodynamical variables, microphysical variables, processes and parametrizations are summarized in Fig. 3 and are described in the next sections.The conversion rates that impact rain formation and evolution are parametrized according to KK00 or according to SB01, Seifert andBeheng (2006, hereafter SB06), andSeifert (2008) (hereafter S08).The cloud water specific humidity is diagnosed from the cloud condensation and evaporation scheme.

Cloud droplet condensation and evaporation
The cloud water specific humidity q c is diagnosed from the resolved pressure, the temperature and the total specific humidity using an "all or nothing" cloud adjustment scheme i.e. it is assumed that there is no cloud water present in an unsaturated grid box, while all moisture exceeding the saturation value q s is cloud water: To calculate q s ≡ q s ( T , p ), an implicit equation needs to be solved, because T is not directly available and has to be diagnosed from the prognostic variables θ l and q t .The temperature T is approximated with help of the liquid water temperature T l , which is equal to: Following Sommeria and Deardorff (1977), q s ( T , p ) is found through a Taylor expansion around q sl ≡ q s ( T l , p ): and the higher order terms are neglected.For ideal gases, the saturation specific humidity is expressed through the saturation vapor pressure as: with R v = 461.5 J kg −1 K −1 denoting the gas constant for water vapor.It can be solved in very good approximation as: with constants e s0 = 610.78Pa, T trip = 273.16K, a = 17.27 and b = 35.86K.After having substituted in Eqs.(56-58) into the truncated Taylor expansion Eq. ( 55) we obtain for the saturated specific humidity: and finally the cloud water specific humidity can be calculated with Eq. ( 53).When higher accuracy is necessary, the procedure can be applied iteratively.

Cloud droplet sedimentation
The cloud droplet sedimentation process has an impact on cloud evolution by reducing entrainment at stratocumulus cloud top (Ackerman et al., 2004;Bretherton et al., 2007).Its cloud water specific humidity source term can be expressed as the derivative of a sedimentation flux.The latter is parametrized by assuming a Stokes law to calculate the cloud droplets terminal velocity and a log-normal distribution to represent the cloud droplet spectrum (Ackerman et al., 2009), which lead to the following expression: with ρ w = 1000 kg m −3 , k St = 1.2×10 8 m −1 s −1 and the lognormal geometric standard deviation parameter σ gc is set to 1.3 (Geoffroy et al., 2010).

Rain drop processes
The precipitation processes are parameterized as functions of the local thermodynamical state.Thus they are valid only for simulations where microphysical fields are explicitly resolved, as is the case in LES.To be able to neglect subfilterscale fluctuations, resolution must not be more than 200 m horizontally and a few ten of meters vertically.
In slightly precipitating clouds, most of the falling mass is contained in particles smaller than 50 µm in radius, also referred to as drizzle.In Khairoutdinov and Kogan (2000) scheme, the limit between the cloud category and the rain category is set at the radius value of 25 µm which permits consideration of drizzle in the precipitating category, which can have significant impact on the evolution of the boundary layer.This scheme is empirically based: it has been tuned with spectra derived from 3-D simulations of stratocumulus clouds using a coupled LES-bin microphysics model.Thus it is valid only for stratocumulus clouds.Because a description of that scheme is fully given in Khairoutdinov and Kogan (2000), it is not described here.The SB01 scheme assumes the limit at the separating mass value x 0 of 2.6×10 −10 kg which corresponds to a separating radius r 0 of the order of 40 µm.Thus the SB01 scheme is more suitable for heavily precipitating clouds, in which most of the falling mass is contained in millimeter size particles.The parametrized collection rates are expressed in function of the microphysical variables by analytical integration of the stochastic collection equation (SCE) and assuming analytical distributions to represent the hydrometeor spectra.A correction function is added to the autoconversion and accretion rate, that take in account the evolution of the cloud droplet spectra due to conversion of cloud water in rain water.
The rain drop size distribution (RDSD) is assumed to be a Gamma distribution: where r is the rain drop radius.N 0 and the slope parameter λ r can be expressed as a function of the prognostic variables and µ r .In autoconversion and accretion parametrizations, µ r has been set to the Marshall and Palmer (1948) value (i.e.0) and is fixed because the parametrizations have been tuned with such a value using spectra derived from 1-D simulations using a coupled LES-bin microphysics model.The value of the shape parameter µ r is parametrized in function of the rain water content (Geoffroy et al., 2010):

Autoconversion from cloud droplets to rain drops
Autoconversion is the process that initializes the rain drop spectra.
After analytical integration of the SCE and adding of correction function, SB06 obtained the following parametrized expressions: with: and k au = 10.58×10 9 m 3 kg −2 s −1 (Pinsky and Khain, 2002, SB06) and x c is the mean mass of the cloud droplet distribution.ν c is parametrized according to Geoffroy et al. (2010): New drizzle drops are assumed to have a radius equal to the separating radius r 0 .Thus the rain number concentration source term due to autoconversion is:

Accretion of cloud droplets
The growth rate of rain drops by collecting cloud droplets is taken to be a function of the cloud and rain water contents (SB06): with: and k acc =5.25 m 3 kg −1 s −1 .

Rain drop selfcollection
The rain number concentration decreases because of the selfcollection process, i.e. interaction between rain drops together to form larger rain drops.Its parametrization is expressed as the following (SB06): with k sc = 7.12 m 3 kg −1 s −1 and κ sc = 60.7 kg −1/3 .

Break-up of rain drops
The break-up of rain drops into smaller rain drops is applied for spectra with a mean volume radius r vr larger than 150 µm following (SB06): with k br = 2000 m −1 and r eq = 550 µm.When r vr becomes larger than r eq the break-up process becomes predominant over the selfcollection process.The strong increase of the break-up process for large mean volume radius is not taken in account.

Rain drop sedimentation
Assuming the Rogers et al. (1993) dependence of rain drop terminal velocity in function of the drop radius, the flux of the rain number concentration and the flux of the rain water specific humidity are (Stevens and Seifert, 2008): with a = 9.65 m s −1 , b=9.8 m s −1 , c = 1200 m −1 (see S08).

Rain drop evaporation
The tendency of the rain water specific humidity due to evaporation is expressed by integration of the drop growth rate by vapor diffusion formulation (S08): where (x) is the gamma distribution, Sc, the Schmidt number, a v and b v are ventilation factor coefficients with the following values: Sc = 0.71, a v = 0.78 and b v = 0.308 (Pruppacher and Klett, 1997).
The tendency of the rain drop number concentration due to evaporation is assumed to be (S08): with γ = 0.7 (A.Seifert, personal communication, 2008).Note that a 0 value of γ means that no rain drop disappear during evaporation.A value larger than 1 would be possible if a large number of little rain drops totally evaporate in presence of large drops.

Radiation schemes
The net radiative heating is equal to the (downward pointing) radiative flux divergence (per unit wave length) integrated over all wavelengths ν: DALES includes a multi-waveband radiative transfer model.It needs information of the vertical profiles of temperature, humidity and ozone up to the top of the atmosphere.To reduce the computational cost of the radiative transfer calculation, DALES has implemented the Monte Carlo Spectral Integration (Pincus and Stevens, 2008), where at each grid point and at each time step the radiative flux is approximated by the radiative flux of one randomly chosen waveband, or a randomly chosen part of that waveband where all absorption coefficients are similar.A complete discussion of the radiative transfer model can be found in Fu and Liou (1992); Fu et al. (1997).DALES also includes a simple parameterization for the vertical component of the longwave radiation and of the www.geosci-model-dev.net/3/415/2010/Geosci.Model Dev., 3, 415-444, 2010 T. Heus et al.: The Dutch Atmospheric LES shortwave radiation through computationally cheap analytic approximations of the Mie theory, that maintain sufficient accuracy for most purposes.In the parameterized radiation scheme, radiative transfer is computed at every single column of the LES, neglecting horizontal radiative transfer.

The GCSS parameterization for longwave radiation
The absorptivity of longwave radiation is controled by the liquid water path (LWP), The net longwave radiative flux F rad L is linked to the liquid water path through an analytic formula, where k is the absorption coefficient, and F L (z top ) and F L (0) represent the total net longwave radiative flux at the top of the cloud and the cloud base, respectively.Larson et al. (2007) discuss the validity of this parameterization in detail.They conclude that when the parameterization constants are optimized for individual stratocumulus cases like the ones set up by Duynkerke et al. (1999Duynkerke et al. ( , 2004)), and Stevens et al. (2005), the formula yields remarkably accurate fluxes and heating rates for low clouds.
To study the effect of longwave radiative on the generation of turbulence, but in the absence of latent heat release effects that occur in a real liquid water cloud, one can substitute the liquid water path by the vertical integral of a passive scalar.This so-called "smoke" cloud has an initial concentration set to unity in the boundary layer and zero above (Bretherton et al., 1999b).The liquid water path in the longwave radiation Eq. ( 81) is then replaced by the smoke path, which can be computed by substituting q c by the smoke concentration s in Eq. ( 80).For a smoke absorptivity k = 0.02 m 2 kg −1 one obtains similar cooling rates as in stratocumulus (Bretherton et al., 1999b).It should be noted that unlike liquid water, smoke is a conserved quantity.This means that if smoke is transported by turbulence into the inversion layer, it will cause a local cooling tendency in this layer.

The delta-Eddington model for shortwave radiative transfer
In the shortwave band the cloud optical depth τ is the most important parameter defining the radiative properties of clouds, Here r e defines the cloud droplet effective radius, i.e. the ratio of the third moment to the second moment of the droplet size distribution (Stephens, 1984).A constant value for r e is used, and for marine boundary layer clouds is r e = 10 µm, which was observed for stratocumulus over the Pacific Ocean off the coast of California during FIRE I (Duda et al., 1991).Cloud droplets scatter most of the incident radiation into the forward direction.This asymmetry in the distribution of the scattering angle is measured by the first moment of the phase function, and is commonly refered to as the asymmetry factor g which is taken g = 0.85.The radiative transfer for shortwave radiation in clouds is modeled by the delta-Eddington approximation, in which the highly asymmetric phase function is approximated by a Dirac delta function and a two term expansion of the phase function (Joseph et al., 1976).The physical interpretation of this approach is that forward scattered radiation is treated as direct solar radiation.
The ratio of the scattering coefficient Q s to the extinction coefficient Q e is called the single scattering albedo ω 0 = Q s /Q e , and is unity for a non-absorbing medium.Following Fouquart (1985), with τ t the total optical depth in a subcloud column.Due to multiple scattering this small deviation of ω 0 from unity leads to a non-negligible absorption of shortwave radiation.The delta-Eddington equations are exactly the same as the Eddington equations (Joseph et al., 1976) with transformed asymmetry factor g, single-scattering albedo ω 0 and optical depth τ substituted by primed quantities: ) For constant ω 0 and g the delta-Eddington equation can be solved analytically (Shettle and Weinman, 1970;Joseph et al., 1976): with: and µ 0 = cosα 0 for a solar zenith angle α 0 .The values of the constants C 1 and C 2 in Eq. ( 87) are calculated from the Geosci.Model Dev., 3, 415-444, 2010 www.geosci-model-dev.net/3/415/2010/boundary conditions.A prescribed value for the total downward solar radiation (parallel to the beam) determines the upper boundary condition at the top of the cloud F 0 .In addition, it is assumed that at the ground surface a fraction of the downward radiation reaching is reflected back by a Lambertian ground surface with albedo A g .See for further details Shettle and Weinman (1970) and Joseph et al. (1976).
The delta-Eddington solution is applied in every column using the local cloud optical depth.A study by de Roode and Los (2008) of the cloud albedo bias effect showed a good agreement between results obtained with the delta-Eddington approach and from the I3RC Monte-Carlo model (Cahalan et al., 2005) that utilizes the full three-dimensional structure of the cloud field.

Other forcings and sources
Large-scale forcings and sources, such as the mean geostrophic wind u g , the large-scale subsidence w s , and the horizontal advective scalar transport may depend on height and time.The effects of large-scale subsidence are calculated using the slab-averaged scalar profiles and a prescribed subsidence velocity w s (z,t) (see, e.g.Siebesma et al., 2003): Optionally, the slab-averaged prognostic variables can be nudged with a relaxation time scale t rel to a prescribed (time depending) value ϕ rel : analogous to large-scale forcings in single column models (see Neggers et al., 2010).The application of S rel ϕ to the horizontal mean ϕ , instead of to the individual values of ϕ, ensures that room for variability within the LES domain remains, and the small-scale turbulence will not be disturbed by the nudging.

Flow over tilted surfaces
To simulate flow over a sloped surface under an angle α (> 0), a coordinate transformation is performed; computations are then done in a system (s,y,n) , with s and n are the coordinates in the down slope direction and perpendicular to the slope, respectively.Under the assumption that the flow can be considered homogeneous along the slope (see Sect. 3.5), only the buoyancy force is directly dependent on s.The original gravitational forcing g θ 0 θ v δ i3 needs to corrected.The reformulated gravitational forcing is equal to: As of yet, the SFS model is not adjusted, which limits the accuracy of the simulations, especially for bigger slope angles and very stable conditions.To accommodate the periodic horizontal boundary conditions for slope flow, we follow Schumann (1990) in splitting each scalar field ϕ in an ambient component ϕ a that incorporates the z dependency of the mean state, and a deviation ϕ d with respect to ϕ a .
The ambient profile ϕ a (s,y,n) is equal to slab averaged value of the scalar ϕ (z), where the brackets still stand for an average over the x-and y-directions.Defining z = 0 at s,n = {0,0}, the ambient value of the scalar is equal to: The deviation ϕ d is now homogeneous along the slope surface direction, and periodic boundary conditions can be applied on it.Currently, this splitting procedure is not implemented for the total specific humidity, focusing slope flow studies exclusively on the dry boundary layer for now.

Chemically reactive scalars
DALES is equipped with the necessary tools to study the dispersion of atmospheric compounds using the Eulerian and Lagrangian framework and their chemical transformation.The Lagrangian framework is explained in Sect.2.13.2.In the Eulerian approach, a line or surface source of a scalar or a reactant is included to mimic the emission of an atmospheric constituent in the ABL flow allowing the calculation and analysis of the diagnostic scalar fields (Nieuwstadt and de Valk, 1987).If the atmospheric compounds react, the source or sink term in Eq. ( 7) needs to be included in the numerical calculation.For a generic compound ϕ l , this reaction term reads: The respective terms P(t,ϕ m ) and L(t,ϕ m ) are nonnegative and represent production and loss terms for atmospheric constituent ϕ l reacting on time t with the n number of species ϕ m it is reacting with.
In DALES, we compute the chemical source term S ϕ l using the chemical solver TWOSTEP extensively described and tested by Verwer (1994) and Verwer and Simpson (1995).In short, this chemical solver uses an implicit second-order accurate method based on the two-step backward differentiation formula.Since in atmospheric chemistry we are dealing with chemical system characterized by a wide range of chemical time scales, i.e., stiff system of ordinary differential equations, the time step is adjusted depending on the chemical reaction rate.
A simple chemical mechanism serves us as an introduction of the specific form of P(t,ϕ k ) and L(t,ϕ k ).Atmospheric www.geosci-model-dev.net/3/415/2010/Geosci.Model Dev., 3, 415-444, 2010 chemistry mechanism are composed of first-and secondorder reactions.Third-order reactions normally involve water vapor or an air molecule, for instance nitrogen or oxygen.Due to the much larger concentration of these compounds than the reactant concentration, third-order reaction rates are normally expressed as a pseudo second-order reaction, i.e., where [M] is a molecule of H 2 O or air.Therefore, the generic atmospheric chemical system consisting of a first-and a second-order reaction reads: where a, b and c are atmospheric compound concentrations, j and k are the first-and second-order reaction rate.For reactant a the L and P are, respectively: The photodissociation rate j depends on the ultraviolet actinic flux and specific photodissociation properties of the atmospheric compound.Therefore, in DALES j is a function on the diurnal variability (latitude, day of the year) and the presence of clouds.The j -values are updated every time step.The cloud influence on the actinic flux is implemented using a function that depends on the cloud optical depth (Eq.82) (Vilà-Guerau de Arellano et al., 2005).The reaction rate k depends on the absolute temperature, on the water vapor content and the pressure.Depending on the reaction, several reaction rate expressions can be specified at DALES.Moreover, the generally very low concentrations of chemical species in the atmosphere allows us to neglect the heating contribution of the reactions on the liquid water potential temperature θ l , or on the water content q t and q r .For the chemical solver, it is essential that the concentration of the species is non-negative.Therefore, the entire numerical discretization for the reactants, spatial and temporal integration of advection and diffusion and temporal integration of the chemistry, has to satisfy the following three numerical properties: it has to be conservative, monotone and positive definite.Of the advection schemes that are implemented in DALES, the kappa scheme (see Sect. 2.7) fulfills these properties.
The chemistry module is designed to be flexible in order to allow study of different chemical mechanisms.Required input parameters include the number of inert scalars, and of chemical species, their initial vertical profiles and surface fluxes, and a list of chemical reactions, together with the reaction rate functions.More information on the chemistry module can be found at Vilà-Guerau de Arellano et al. (2005) and Vilà-Guerau de Arellano et al. (2009).

Statistics
In DALES, standard output includes time series and slabaveraged profiles of the main variables, the (co-) variances, and of the resolved and SFS-modeled fluxes.The modular set-up of the code facilitates inclusion of many other statistical routines, specifically aimed at the purposes of a particular research question.Sharing such code with the community leaves the code base with a rich palette of statistics, including specific routines that focus on the details of, for example, radiation, cloud microphysics, or the surface layer.A few examples of the statistical capabilities of DALES are given below.

Conditional sampling
Conditionally averaged profiles can be found by defining a mask M, which is equal to 1 or 0, depending on whether a set condition is true or false, respectively.Frequently used sampling conditions are, for instance, clouds (q l > 0), areas of updrafts ( w > 0), areas of positive buoyancy ( θ v > θ v ), and any combination of these conditions.New definitions of the mask M are possible with small adjustments of the code.

Lagrangian statistics
While the Eulerian formulation of the LES favors a Eulerian frame of reference for statistics, many problems can greatly benefit from a Lagrangian approach.This holds in particular for studies of entrainment and detrainment, since these problems can often be stated as a study on the evolution in time of a parcel of air.To this end a Lagrangian Particle Dispersion Model (LPDM) has been implemented into DALES.Within this model, massless particles move along with the flow.Since each of the particles is uniquely identifiable, the origins and headings of the particles (and of the air) can be captured.
The position of a particle x p is determined using: where u is the LES-resolved velocity linearly interpolated to the particle position, and u is an additional random term that represents the SFS-velocity contribution.This term is especially important in regions where the SFS-TKE is relatively large, such as near the surface or in the inversion zone.The calculation of u follows Weil et al. (2004), and was tailored for use in LES with TKE-closure.It is implemented in DALES as follows: where C 0 is the Langevin-model constant (Thomson, 1987) that has been set to 6; f s is the slab-averaged ratio between SFS-TKE and total TKE.dξ is a Gaussian noise to mimic the velocity field associated with the subfilter turbulence.Boundary conditions are periodic in the horizontal directions, and emulate the LES boundary conditions at the top and bottom of the domain.Particles are reflected (w p changes sign) should they hit the top or bottom.For time integration, the third order Runge Kutta scheme is again used.The LPDM was validated by Heus et al. (2008) for a cumulus topped boundary layer and additionally by Verzijlbergh et al. (2009) for a scalar point source emission in different clear and cloud-topped boundary layer flows.

Transport, tendencies and turbulence
To study the mechanisms driving the development of the ABL, tendency statistics are included that diagnose slab average profiles of every forcing and source term in Eqs. ( 6) and ( 7).Where necessary, the individual terms of the underlying equations can also be diagnosed, such as for the SFS-TKE, radiation or microphysical components.Fluxes and co-variances of the main variables are also calculated.
To understand the turbulence in the boundary layer it is interesting to analyze the resolved turbulence kinetic energy budget E, which, under horizontal homogeneous conditions and neglecting subsidence, can be based on the turbulence kinetic energy budget as given by e.g.(Stull, 1988): where the double prime indicates a deviation from the slabaverage, θ 0 is a reference virtual potential temperature.The left-hand side term represents the total tendency of turbulent kinetic energy.The right-hand side terms are, respectively, the shear production, the buoyancy production, the turbulent transport, the pressure correlation, and the viscous dissipation term: where K m is the SFS eddy diffusivity.Due to the staggered grid used in DALES each variable entering in the budget terms is evaluated at a different position.In order to correctly build up the different terms, several interpolations have to be performed, which have to be consistent with the spatial discretization of the model.Due to these numerical issues, the budget is not fully closed, although the residual is small compared to the physical terms (see Fig.  for the budget of E in a sheared CBL).In order to further reduce this residual term, a method based on Gao et al. (1994) is currently in development.

Convective boundary layer
One of the most common test cases for an atmospheric LES is the dry convective boundary layer (CBL).In a CBL a positive heat flux at the surface destabilizes the air resulting in a vigorous turbulence which mixes (thermo)dynamic quantities like heat and momentum over the entire depth of the boundary layer, and which comprises eddies that vary over a wide range of scales, i.e. from the depth of the boundary layer (∼km) down to the Kolmogorov-scale (∼mm).But because the largest scales of motion control most of the vertical transport (e.g. the vertical fluxes of heat and momentum), it is reasonable to fully resolve the large scales on a resolution of ∼ 10−100 m, and account for the scales of motion smaller than the grid scale using the subgrid model (such as Eq.25).
Probably the most defining feature of a CBL is the fact that the mixed layer is not confined by a rigid lid (such as Rayleigh-Bénard convection), but that it is capped by an inversion, a marked increase of the potential temperature with height.As such the mixed-layer depth z i (the top of which is often defined as the height of the maximum gradient in temperature, or as the height of the minimal buoyancy flux), is not fixed, but grows in time: thermals impinging on the inversion cause overlying free tropospheric air to be entrained into the mixed-layer, the depth of which therefore increases.The rate of growth is called the entrainment rate w e , a key unknown in weather, climate and air quality  Sullivan et al. (1998).For extra information see Table 3.All results are averages over hour 3-4 and spatially averaged in the horizontal.(a, d): average temperature profile (thin line in a shows the initial temperature profile).(b, e): normalized heat-flux profiles, resolved (thin line), subgrid (dashed) and total (solid line).cf.turbulence statistics of resolved velocities σ 2 u = ũ ũ (solid line), σ 2 v = ṽ ṽ (thin line), σ 2 w = w w (dashed line), and subgrid contribution 2e/3 (dotted line).
Fig. 6.Vertical profile of the total w θ v + w θ v (solid line) and subfilter-scale contribution w θ v (dashed) of the virtual potential temperature flux obtained after four hours of simulation by DALES2.0 (black) and DALES 3.2 (red) with the same physical conditions and advection scheme (2nd order central differences).
models.Large-Eddy Simulation provides a powerful tool to make a comprehensive study of entrainment (see e.g., Sullivan et al., 1998;Fedorovich et al., 2004a) and investigate Table 3. Simulation details of the two simulated CBLs: weak inversion case (W06) and strong inversion case (S24).z i (0) and θ(0) denote the initial mixed-layer depth and initial temperature jump, respectively.the dependencies on for example the inversion jump θ v , the surface heat flux F s,θ v and the actual mixed-layer depth z i (t).Rather than studying the entrainment rate directly, one can also focus on the entrainment flux of heat, in particular the value of the heat flux at the inversion.This approach is followed below for DALES.
To test the performance of DALES for dry convective boundary layers, we simulated two of the cases studied by Sullivan et al. (1998), one with a weak inversion θ v ∼ 0.5 K (their case number W06), and one with a strong inversion θ v ∼ 5 K (case S24).The corresponding surface heat flux, initial mixed-layer depth z i (0) and stratification d θ v /dz of Geosci.Model Dev., 3,2010 www.geosci-model-dev.net/3/415/2010/ the overlying layer, are given in Table 3.In both cases there is no mean wind and hence no (mean) shear.Note that W06 was initiated without an inversion jump.For S24 the initial inversion thickness amounted to 120 m (linear interpolation between 300 K and 308 K over 120 m).Both simulations were conducted on a grid of N x = N y = 64,N z = 96, using the same resolution as in the original simulations, x = y = 100 m, z = 20 m.Time-step was variable, and for the advection of all variables the fifth-order scheme (see Sect. 2.7) was chosen.
In Fig. 5 we present the results averaged from hour 3 to 4. Turbulence statistics are normalized using the convective velocity scale where Q 0 = F s,θ v is the surface kinematic heat flux in K m/s, and z i the mixed layer depth (see Table 3).The figures are formatted such that they can be directly compared with the original study by Sullivan et al. (1998).Although W06 was initiated without an inversion, the CBL dynamics is such that it creates its own inversion, as can be seen in Fig. 5a showing the characteristic "steepening" of the temperature profile in the entrainment zone.The strength of the resulting inversion is the same as observed by Sullivan et al. (1998).The same holds for case S24 (Fig. 5d).For both cases also the normalized heat flux profiles displays a value of roughly −0.15 in the entrainment zone, indicative of the entrainment process (Fig. 5b, e).The SFS contribution to the heat flux is rather small in the mixed-layer and near the inversion.The SFS contribution to TKE, on the other hand, extends over the entire layer (Fig. 5c, f); again the magnitude and shape of the SFS-TKE are in very good agreement with the results reported by Sullivan et al. (1998).

Generation of mesoscale fluctuations
where the terms on the rhs indicate the production, transport and dissipation of scalar variance, respectively.If the vertical flux of a scalar changes sign at some level in the mixed layer, and if the mean vertical gradient of the scalar changes sign at a different height, the vertical flux will be counter the mean vertical gradient.According to the variance equation this implies that the production of variance will become negative.The large-eddy simulations show a clear correlation between the depth of the countergradient flux layer and a decrease in the production of variance.
In addition, it was found that the growth rate of mesoscale fluctuations of passive scalars is tightly connected to the shape of the vertical profile of the buoyancy flux and its magnitude in particular.This was demonstrated with aid of simulations of a smoke cloud that was radiatively cooled from its top.This cooling drives top-down convection and produces positive buoyancy fluxes in the bulk of mixed layer, where a larger cooling rate causes larger buoyancy fluxes.The largest growth rates of mesoscale fluctuations of passive scalars were found for the case with the largest values for the buoyancy flux.The main message that these studies carry is that turbulence alone is capable to generate mesoscale fluctuations of passive scalars, indicating that for these quantities a spectral gap does not exist.

Sheared convective boundary layer
To analyze the influence of wind-shear characteristics on the evolution of the CBL, long simulations and large domains are necessary to fulfill a quasy-stationarity flow pattern that matches with the prescribed surface fluxes, and to resolve the expected pattern for forced convection (Khanna and Brasseur, 1998).With DALES, resolutions up to 25 m and 6 m in the horizontal and vertical directions, respectively were considered.
The studies of the sheared CBL focus on the influence of the wind shear on the boundary layer growth due to the modification of the entrainment fluxes (Pino et al., 2003); on identification and parameterization of the main physical mechanisms that control the entrainment heat flux (Kim et al., 2006;Pino et al., 2006b); on the role of shear and the inversion strength in the decay of convective turbulence during sunset (Pino et al., 2006a); and most recently on how to parameterize the different terms of the TKE budget by means a first order jump mixed layer model (Pino and Vilà-Guerau de Arellano, 2008).In an intercomparison study of the sheared CBL in different wind regimes by Fedorovich et al. (2004b), a previous version of the model showed larger entrainment fluxes than the other codes.Consequently, a drier and warmer boundary layer was obtained.In comparison with this older version of DALES, DALES 3.2 shows smaller entrainment fluxes (see Fig. 6).Presumably, this is due to the improved numerical scheme; as the κ advection scheme performs better in combination with the Runge Kutta integration scheme than with the previously used leap frog scheme Hundsdorfer et al. (1995).
Among the results mentioned above we would like to emphasize first the influence of the shear in the boundary layer www.geosci-model-dev.net/3/415/2010/Geosci.Model Dev., 3, 415-444, 2010 growth by using LES and observations (Pino et al., 2003), and second the influence of the wind shear in the characteristics length scales during afternoon decaying convective turbulence (Pino et al., 2006a).It was shown there that the enhancement of the entrainment heat flux caused by the wind shear at the inversion zone is responsible for an increased boundary layer growth rate.Neglecting this wind shear in the parameterizations of the entrainment heat fluxes would result in a significant underestimation of the boundary layer depth (see Fig. 7).

Stable boundary layers
In the context of LES, one of the characteristics of stable boundary layers (SBLs) is the mere absence of large eddies (as compared to the height above the surface, or the depth of the boundary layer; see e.g. the spectra presented in Kaimal and Finnigan, 1994).The stable stratification suppresses vertical motion and transfers turbulence kinetic energy into turbulence potential energy (defined as 1 2 ∂θ ∂z −1 θ 2 , see Zilitinkevich et al., 2007) through the buoyancy destruction term in the TKE equation.Part of that potential energy is released back as turbulence kinetic energy but part is dissipated through the dissipation of temperature variance.Due to these two aspects, the role of the subfilter-scale model tends to be much larger in LES of SBLs than it is for convective or neutral (but sheared) boundary layers.This implies that for the SBL generally much higher resolutions are used than for other simulations.
The first application of (a previous version of) DALES to stable boundary layers was reported by Galmarini et al. (1998) where a slightly different version of the subfilter-scale model was used.
In the context of the GEWEX Atmospheric Boundary Layer Study (GABLS, Holtslag, 2006), a series of model intercomparisons has been organized for SBL cases.In all intercomparisons a single-column model intercomparison case was defined, whereas an LES case was defined in the first and third intercomparison.The first case (Beare et al., 2006) was inspired by the setup of the simulations of Kosović and Curry (2000): a moderately stable boundary layer (with z i /L ≈ 2, where z i is the depth of the SBL (here defined as the height of vanishing shear stress) and L the Obukhov length).The domain size was set to 400 m in all three directions.The roughness length for momentum z 0 was set to 0.1 m and for heat the same roughness length was applied.The lower boundary condition for heat was imposed as a constant cooling rate of 0.25 K/h for the surface temperature.The flow was forced with pressure gradient representative of a geostrophic wind of 8 ms −1 at a latitude of 73 • N.
In total 11 models participated in the intercomparison, being run at resolutions from 12.5 m down to 1 m for some models.DALES participated in the intercomparison at resolutions of 12.5 and 6.25 m.For this paper, the case was re-run at a resolution of 3.125 m.The results are shown in Fig. 8.The results of DALES are clearly within the range of the other models, although the mean shear is stronger than in most models close to the surface and weaker at higher levels in the SBL.Furthermore, the strength of the low-level jet (or, more precisely, the super geostrophic jet, i.e., the wind maximum at the top of the SBL) seems to be slightly less than in the other models.

Cloud topped boundary layer
If there is sufficient moisture in the convective boundary so that the total specific humidity q t exceeds its saturation value q s , condensation processes will initiate and clouds will start to form.Since q s increases exponentially with temperature and as temperature decreases with 10 K/km in the convective boundary layer, clouds typically start to form at the top of the convective boundary layer.They are often referred to as boundary layer clouds, as long as the capping inversion at the top of the boundary layer is strong enough to encapsulate them.As a result they have a limited vertical extent of around 3 km which makes the use of LES well suitable to study the dynamics of boundary layer clouds.
Stratocumulus and shallow cumulus are the two main types of boundary layer clouds that have been simulated extensively in the past with DALES and the schematics of these different types of boundary clouds are depicted in Fig. 9.
Stratocumulus (see Fig. 9b) clouds are low-lying, stratiform clouds often covering the sky completely, with a thickness of only several hundreds of meters, capped by a strong inversion.The turbulence that maintains the well-mixed profiles of the conserved variables q t and θ l is mainly driven from the top of the stratocumulus deck due the longwave radiative cooling in addition to local cooling and heating due to condensation and evaporation of cloud droplets.In contrast, shallow cumulus clouds (see Fig. 9c) occur as a population of separated small cauliflower shaped clouds with a cloud base height at around 1 km and a maximum vertical extend of around 2 km.These clouds generally only cover 10 to 30% of the sky.Shallow cumulus clouds usually form on top of the dry rising thermals in the subcloud layer and are dynamically characterized by strong vertical motions due to the condensational heating resulting in inner cloud cores that are positively buoyant with respect to the (dry) environment.As a result the stratification in terms of the mean profile of θ v is stable with respect to vertical displacements of unsaturated test parcels and unstable with respect to saturated test parcels.This effect, often referred to as conditional instability is unique for moist convection and has no counterpart in dry convection and is responsible for the strong intermittant behaviour of cumulus updrafts.The dynamics of stratocumulus might appear to be simpler than shallow cumulus due to the fact that it is horizontally more homogeneous than shallow cumulus and that it is well mixed in the vertical so that it can be conceptually well described by a simple mixed layer model.However, it is actually harder to to simulate stratocumulus clouds in an LES model due to the strong inversion at the top of the stratocumulus deck where temperature jumps of 10 K over 100 m are not uncommon.Such strong inversions result from the radiative cooling and are difficult to resolve with LES techniques, resulting in unwanted numerical diffusion over this Fig. 10.Hourly-averaged entrainment rates from LESs for four different GCSS cases.The "DALES old" and "other LES" indicate entrainment results as obtained from previous versions of the code.The observed entrainment rates with their uncertainties are also plotted.Because the Eurocs FIRE case is based on a monthly mean climatology, no observed entrainment rates are available for this case.The "DALES old" results were all obtained with the kappa scheme.The results with the current DALES version were obtained with the kappa, second-order and fifth-order advection schemes as indicated in the legend.
interface which can dominate the transport over the inversion interface.On the other hand, in the case of shallow cumulus clouds the interaction with the radiation is not so strong due to the short life span of the clouds, and due to the low cloud fraction.As a result, shallow cumulus clouds are not topped by such strong inversions which simplifies the numerical simulations.Another related simplifying factor is that because of the low cloud fraction the interaction between the clouds and the radiation is not so critical that an interactive treatment of both processes would be essential.
DALES has participated in numerous LES intercomparison studies organized over the last 15 years by the GEWEX Cloud System Studies (GCSS).These intercomparison studies have been set up to serve several purposes.They provide critical evaluations of the participating LES codes and unique data sets to obtain further insights in the dynamics of the cloud topped boundary layer.More specifically these LES data sets have helped in improving the parameterized formulation of these processes in large scale Numerical Weather Prediction (NWP) and climate models.In the coming 2 subsections examples are presented how research with DALES have contributed to the improved knowledge of the physics and dynamics of shallow cumulus and stratocumulus.

Stratocumulus
One of the most critical phenomena in the dynamics of stratocumulus is the entrainment of dry air at the top of the cloud layer.Following the flux-jump relation (Lilly, 1968), the entrainment rate (w e ) determines the turbulent flux at the top of the boundary layer ( w ϕ e ), w ϕ e = −w e ϕ , (105) with ϕ the jump across the inversion.This equation is valid for an infinitesimally thin inversion layer and shows the importance of the entrainment rate on the turbulent fluxes at the top of the boundary layer.The representation of turbulent transport by an LES code therefore critically depends on its capability to produce realistic entrainment rates.
Figure 10 shows the modeled entrainment rates by previous versions of DALES and by other models as reported in intercomparison studies by Duynkerke et al. (1999), Duynkerke et al. (2004), andStevens et al. (2005).Results obtained with the current version of DALES using three different advection schemes are also shown.To facilitate a direct comparison to the previous findings, the new simulations were all run nearly identically to the original case descriptions.The entrainment rates obtained from the LES models represent hourly-average values from the third (ASTEX) or fourth hour of simulations (Eurocs FIRE and DYCOMS RF01).
The entrainment rates from previous versions of DALES were all rather large in comparison to results from other participating LES codes used in the intercomparison studies.An initially large entrainment rate generated by the previous version of DALES led to a rapid thinning and subsequent breakup of the DYCOMS RF01 stratocumulus deck.As a consequence the longwave radiative forcing at the boundary-layer top decreased, explaining the small value for the entrainment rate as shown in the figure for older versions of DALES.It is clear from the figure that entrainment rates are reduced in the current version of DALES.Entrainment rates in simulations with the monotone kappa advection scheme tend to be slightly larger in comparison to the second or fifth order advection schemes.Only a simulation with the latter scheme is capable to maintain a solid stratocumulus cloud deck for the DYCOMS RF01 case.Less sensitivity can be seen for the Eurocs FIRE case, since in this case the stratocumulus deck is generally much more stable as in the Dycoms RF01 case.

Shallow cumulus
A number of interesting and well-documented shallow cumulus cases based on observational studies have been simulated by DALES over the last 10 years.These studies include: non-precipitating steady-state marine shallow cumulus based on the Barbados Oceanographic and Meteorological Experiment (BOMEX) (Siebesma et al., 2003) and on the Atlantic Trade Wind Experiment (ATEX) (Stevens et al., 2001), diurnal cycles of shallow cumulus over land observed on 21 June 1997 at the Southern Great Plains (SGP) site (Brown et al., 2002) and during the Small Cumulus Microphysics Study (SCMS) (Neggers et al., 2003a) and more recently precipitating marine shallow cumulus (van Zanten et al., 2010) such as observed during the Rain in Cumulus over the Ocean (RICO) field study Rauber et al. (2007).All these cases have been used to critically evaluate the DALES results against observations and to help developing and testing theories, conceptual models and parameterizations of shallow cumulus convection.In this section we will give a short overview of the results of these studies.
The first category of these studies is related to cloud geometrical issues.In Siebesma and Jonker (2000) it has been shown that the simulated cumulus cloud boundaries have self-similar or fractal properties that can be characterized by a fractal dimension D f = 7/3.These results are in excellent agreement with observational studies and therefore provide a critical test of the capability of DALES to simulate realistic cumulus clouds.Moreover, these results helped in constructing theoretical scaling to explain why cloud boundaries appear to be self-similar with a dimension of 7/3.Another intriguing cloud geometrical topic is related to the question: what is the shape of the cumulus cloud size distribution?It is well known that shallow cumulus cloud ensembles consist of many small clouds and lesser large clouds but the precise shape of the cloud size distribution is still an open issue.Extensive numerical studies with DALES show that the cloud size density of the simulated cloud populations is described well by a power-law from scales smaller than the standard grid-spacing (50 m) up to scales of typically 1000 m with a power-law exponent of −1.7 (Neggers et al., 2003b).This exponent is comparable to values found in observational studies (Cahalan and Joseph, 1989;Rodts et al., 2003).No convincing theory for the power-law behaviour nor for the scale break has yet been put forward.Finally, more recently analyses with DALES of up-and downdrafts in and around individual cumulus clouds have shown that strong updrafts in individual cumulus clouds are typically surrounded by so-called subsiding shells with persistent downdrafts (Heus and Jonker, 2008).These downdrafts are driven by negative buoyant forces that result from the evaporative cooling of the cloud water.As they surround the clouds along their entire perimeter, the subsiding shells cover a significant area and are therefore found to be responsible for a large part of the downward mass transport (Jonker et al., 2008).
The second category studies is related to transport due to cumulus convection which is one of the important processes that needs parameterization in large scale Numerical Weather Prediction (NWP) and climate models.The time evolution of a moist conserved variable ϕ due to moist convection can be written as where F ϕ is the (upward pointing) turbulent flux.A popular method to parameterize this turbulent flux is through the use of a so called mass flux approach where ρ is the density and the subscript c refers to cloud averaged values of ϕ, ϕ the average of ϕ over the entire horizontal slab, and the mass flux is defined as M ≡ ρa c w c Betts (1975), i.e. essentially the product of the cloud averaged vertical velocity times w c and the fractional cloud area a c .Usually a cloud model is derived to obtain equations for M and ϕ c Within this cloud model the key variables are the fractional entrainment ε and fractional detrainment δ rate.These inverse length scales are measures of the rate of dilution of the cloud ensemble (entrainment) and the rate of air leaving the cloud ensemble (detrainment) and LES results from DALES have been used extensively to diagnose ε and δ on the basis of Eqs. ( 108) and ( 109) (Siebesma and Cuijpers, 1995).This approach has initiated considerable research in developing theories and models of these exchange mechanisms between clouds and environment.From these studies it has become clear that the fractional entrainment rate can be well estimated by the inverse cloud depth (Siebesma et al., 2003).
The fractional detrainment rate δ is typically larger than ε as a result of the fact the cloud fraction a c is in general decreasing with height.
Another useful additional equation often used in cloud models is the vertical velocity equation for the cloud ensemble (Simpson and Wiggert, 1969) 1 2 which describes how buoyancy forces and entrainment processes influence the vertical velocity in the clouds.Adjustable prefactors a and b are introduced in this equation to incorporate pressure perturbation effects and incloud turbulence effects in an implicit way.By using Eqs.( 110) and ( 108) we can derive alternative expressions for the entrainment that are more linked to the dynamics than a reciprocal dependency on the height is: In Fig. 11 we compare the entrainment and detrainment rates based on Eqs. ( 110), (111), and (112) for which estimates of a = 0.6 and b = 1 are used for a large variety of different LES experiments.
T. Comparison of LES derived fractional entrainment and detrainment rates ε q and δ q using φ = q t based on Eq. ( 110) (horizontal axis) versus LES estimates of these rates ε w and δ w based on the vertical velocity equation (Eqs. 111 and 112).
The fact that the results fall reasonably well on the diagonal shows that Eqs. ( 110) and ( 108) are consistent, so that the subscripts of ε w and δ w can be removed and Eqs. ( 111) and ( 112) can be used as well to interpret the exchange rates.It can also be observed that ε can vary considerably between values of 1 ∼ 4×10 −3 m −3 indicating that parameterizations that use a constant value for ε is not a good option.Furthermore it should be noted that the range of variability for δ is much larger 1 ∼ 20 × 10 −3 m −3 .More detailed analysis shows that this large variability is mainly due to the gradient of the cloud fraction with height in Eq. ( 112).This indicates that, in order to have a good estimate of the mass flux M, it is more relevant to have a good parameterisation of δ rather than for ε, a statement already emphasized in de Rooy and Siebesma (2008).In that respect it is surprising to see that most of the research efforts have been concentrated on entrainment rather than on detrainment.

Heterogeneous surfaces
DALES has contributed to the understanding of flow over thermally heterogeneous terrain.The study of van Heerwaarden and Vilà-Guerau de Arellano (2008) addressed the question whether convective cloud formation is more likely to form over a land surface that has a heterogeneous surface flux than over a land surface that is homogeneously heated.
Heterogeneous land surfaces were simulated by creating two stripes of 3.2 km wide at the land surface in the model, as this is the spatial scale at which heterogeneity is considered to modify the turbulent structure of the overlying CBL the most.All fluxes at the land surface were prescribed, with a constant u * of 0 m s −1 .Note, that this study did not take into account the impact of the induced circulation on surface friction, which was not possible in the previous version of DALES that did not allow for local fluxes.Both stripes had the same available energy (sum of sensible and latent heat fluxes), but a different Bowen ratio.The left stripe was characterized by a small Bowen ratio, whereas the right stripe had a large ratio.In the different runs in this study, the Bowen ratio was varied.The LES model was run for four hours; statistics were calculated over the last hour.
The main findings of the study are summarized in Fig. 12 that shows the relative humidity in the CBL and the wind vectors in a case where the free atmosphere is moist (left panel) and in a case where the free atmosphere is dry.
In both cases a secondary circulation (see wind vectors) distributes heat and moisture towards the area that has a relatively large sensible and a small latent heat flux.At these hot spots, strong but moist thermals rise, resulting in a large relative humidity over the area that has the smallest latent heat flux.In case of a dry free troposphere (right panel), the secondary circulation can transport very dry free tropospheric air downwards to the land surface.Therefore, a very low relative humidity is found over the area that has the largest latent heat flux.
To conclude, the study showed that heterogeneity results in a situation that is more favorable for cloud formation, regardless of the specific humidity of the free troposphere.Using a similar set up, Górska et al. (2008) use DALES to determine the role of heterogeneity on the carbon dioxide distribution.The study points out of the need to redefine aircraft measurements strategies above non-uniform surfaces.

Atmospheric flow over sloping surfaces
Compared with the many successful large-eddy simulations of the boundary layer over flat terrain, as of yet only a few simulations of the ABL over sloping surfaces have been carried out.One of the problems concerning the simulation of slope flow, is that the potential temperature as well as the Axelsen and van Dop (2010) performed a model validation by comparing simulation results to observations from two glaciers.They found that the simulated profiles of temperature and downslope velocity were quantitatively in agreement with the observations.An example is given in Fig. 13.Near the surface the downslope velocity increases with height and reaches a maximum at a height of 4 m.Above the wind maximum height, the downslope velocity decreases with height.The figure shows that near the surface the simulated and observed velocity profiles agree, but above the wind maximum the model underestimates the velocity.The profile of the simulated potential temperature is also seen to agree rather well with the mast measurements, but that there is a systematic offset between the balloon measurements and the simulated potential temperature.

Dispersion and chemically reacting flows
We summarize here the main research results achieved in the field of turbulent dispersion and chemical transformations using DALES.The plume dispersion main characteristics T. Heus et al.: The Dutch Atmospheric LES and statistics under different ABL flow conditions have been thoroughly investigated using DALES.Dosio et al. (2003Dosio et al. ( , 2005) ) and Dosio and Vilà-Guerau de Arellano (2006) investigated the plume dispersion in the dry CBL from Eulerian and Lagrangian perspectives.Based on DALES results, they derived a parameterization to include the effect of shear on the plume spreading, studied the validity of Taylor's diffusion theory for horizontal and vertical dispersion, and separated the contributions of small-and large-scales on the plume evolution, both from an absolute coordinate system as well as relative to the plume's center of mass.Verzijlbergh et al. (2009) extended this study to determine the influence of stratocumulus and shallow cumulus on the turbulent dispersion properties and related to turbulent structures like skewness of the vertical velocity.As an example, Fig. 14 shows the vertical concentration characteristics and the location of the maximum concentration under different ABL conditions: dry convective boundary layer, stratocumulus and shallow cumulus Verzijlbergh et al. (2009).Similarly to turbulent dispersion, the chemical transformations in the ABL are influenced by the characteristics of the turbulent flow.This turbulence control is particularly important when the turbulent time scale (τ t ) and the chemistry time scale (τ c ) have similar values, i.e., the order of magnitude of the Damköhler number (τ t /τ c ) is O(1).Under this regime, the species are chemically transformed at a different reaction rate depending on the way species are introduced in the ABL, premixed or non-premixed, and the turbulence intensity to mix chemical species.Key tropospheric chemical reactions involving species such as nitric oxide and certain biogenic hydrocarbons like isoprene are therefore controlled by turbulence.Following Schumann (1989), Petersen et al. (1999) and Petersen and Holtslag (1999) studied by means of LES how the transport and mixing of reactants in the CBL is influenced by the presence of vigorous thermals and subsidence motions.Based on the DALES results, they suggested a parameterization to represent the fluxes and covariance of reactants in large scale chemistry transport models.The research was extended to study more complex mechanism under non-uniform emissions of the reactants (Krol et al., 2000;Vilà-Guerau de Arellano et al., 2004).To further study the influence of the reactivity on high-order moments, a spectral analysis showed that the reactant variability (variance) depends strongly on the reaction rate (Jonker et al., 2004).The analysis was done using the DALES simulation of a turbulent flow reacting according to the scheme (R1-R2).These results showed large variations in the characteristic length scale as a function of the Damköhler number and the state of the chemical equilibrium.
To improve parameterizations in large-scale atmospheric chemistry models, Vinuesa andVilà-Guerau de Arellano (2003, 2005) proposed an expression of an effective reaction rate (k eff ) that takes into account explicitly the influence of turbulent mixing on the reaction rate.
The moist and optically thick boundary layer clouds can also influence atmospheric chemistry.DALES was used to study the combined effect of turbulence and radiation on simple chemical mechanism in a dry smoke cloud (Vilà-Guerau de Arellano and Cuijpers, 2000) and shallow cumulus (Vilà-Guerau de Arellano et al., 2005).Figure 15 shows the cloud water content and the photostationary state in a CBL developed over land characterized by the presence of shallow cumulus.This state quantifies the effect of the physical processes (turbulence and radiation) on the atmospheric chemistry.For such reactants as nitric oxide (NO), ozone (O 3 ) and nitrogen dioxide (NO 2 ), it is defined as = (k[NO][O 3 ])/(j [NO 2 ]).Departure from the value = 1 indicate perturbations of the chemical equilibrium either by radiation or turbulent processes.

Resolution dependencies and convergency
Given the pretention of LES to resolve the significant scales of turbulence, one would expect that the outcomes could be independent from the exact resolution.Indeed, in the bulk of a typical convective boundary layers at common resolutions (∼ 10 m), the energy contained in the unresolved turbulence is an order of magnitude smaller than the energy of the resolved turbulence.However, given the limited amount of computational time available and the numerous complex processess we often want to simulate, full convergency is rarely reached or even desired.For the idealized atmosphere that DALES attemps to model, qualitative agreement with higher resolutions usually suffices.Nevertheless, it is important to at least realize where DALES results converge and where they do not.Especially in the vicinity of strong (local) gradients, such as close to the surface and around cloud interfaces convergence is not easily reached.
Although, due to a lack of strong gradients in the crucial regions, cumulus topped boundary layers only need moderate resolutions to get the lower moments of the domain-averaged statistics converge, one needs to be more carefull when investigating more detailed around the cloud interface.For example, Heus et al. (2009) showed that to obtain correct values for the in-cloud vertical mass flux in cumulus clouds, much finer resolutions are necessary than are commonly used in for instance intercomparison studies.When the correct spatial mass flux distribution is desired, DALES shows convergence at a horizontal resolution of 25 m (see Fig. 16).
As shown by Stevens et al. (2005), the strong temperature and moisture gradients at the top of a stratocumulus layer are difficult to mimick in LES.A further complication of the matter is that since stratocumulus convection is for a significant driven by cloud-top cooling, a misrepresentation of these cloud-top gradients can severly affect the state of the entire cloud layer.As an illustration, in Fig. 17  Departure from the standard resolution as prescribed by the intercomparison also yields improvements.The original aspect ratio of a grid cell was much smaller than 1 ( z/ x = 5 m/50 m = 0.1, and the current length-scale formulation of the SFS scheme gives more mixing in such situations.A coarser grid with better aspect ratio (dark gray line; x = 12.5, z = 10 m) can give results that are close to the observed values for similar computational cost.At a resolution of x = 25 m, z = 20 m, the length scale λ = ( x y z) 1/3 equals to λ of the intercomparison setup, and the liquid water path is comparable or slightly better, but at much less computational cost.

Outlook
As was shown in this paper, DALES can provide reliable results for a multitude of atmospheric conditions, and there are many alleys of study that can be pursued with DALES 3.2.In the field of cloudy boundary layers, very fine grid spacing can be used to reliably resolve most of the dynamics within and around the cloud.Simulations on relatively large horizontal domains (∼ 25 km) can mimic the physics in an area similar to a single column of a regional or global model.On that scale, LES is well capable of variability studies that are necessary to improve the GCMs, and to study the impact of GCM grid refinement.For other studies, LES can provide spatial and temporal turbulence characteristics that cannot be easily retrieved from measurements alone.This is always a role that LES can play, but it can be especially important in spatially anisotropic or inhomogeneous situtations, such as in the fields of flow over sloping or heterogeneous surfaces.
While there are many plans to use DALES in its current state, ongoing improvement of the code is also planned.In the near future, we aim to be able to run DALES in more diverse and more realistic scenarios than what was shown in this paper.Furthermore, we aim to focus on studies that makes integrated use of several of the features of DALES.
As was shown throughout many parts of the applications section, LES could still benefit from a better representation of anisotropic turbulence around steep gradients and inversion layers, specifically in stable boundary layers, dry convective boundary layers, and stratocumulus layers.Increasing computer power and resolution could end up simply resolving these gradients in the future, but more intelligent subfilter-scale modeling could also give a significant contribution in solving this problem.This is especially important in critical stratocumulus cases, where entrainment of relatively dry and warm air leads to buoyancy reversal.
To study the interactions between the various components of the model, we strive to have the modules as interactive as possible.This could for instance lead to better under-standing of coupling mechanisms between radiative forcings and the surface conditions, coupling between radiation and chemistry, or between chemistry and cloud and aerosol formation.
Jonker et al. (1999)  performed a simulation of a dry convective boundary layer on a large horizontal domain and demonstrated that the fields of passive scalars can be dominated by fluctuations at length scales much larger than the boundary layer depth.However, it was found that scalars for which the entrainment to surface flux ratio is close to ∼ −0.25, as is a typical value for the buoyancy flux in a clear convective boundary layer, do not exhibit mesoscale fluctuations.In subsequent studies byde Roode et al. (2004a) andde Roode et al. (2004b) the tendency equation for the variance of scalars was analyzed.For an arbitrary passive scalar ϕ it reads

Fig. 7 .
Fig. 7.Boundary layer height z i observed by radiosondes launched at different facilities of ARM campaign (symbols) and obtained by means of LES: without shear (black), including a constant geostrophic wind of 10 ms −1 in the east-west direction (green), and prescribing the observed mean wind (red).Adapted fromPino et al.  (2003).

GeosciFig. 8 .
Fig. 8. Profiles of mean wind speed (top) and potential temperature (bottom) for the first GABLS1 LES intercomparison (average over 9th hour of simulation).Solid black line: DALES result at 3.125 m resolution and c f = 2.0; grey lines: results of other participants at 3.125 resolution.

Fig. 9 .
Fig. 9. Schematic overview of the different types of boundary layer clouds.

Fig. 12 .
Fig. 12.Cross section of the 1-h-averaged relative humidity RH for a case with a moist free troposphere (left) and a case with a dry free troposphere (right).The horizontal coordinates are scaled by the patch size λ and the vertical coordinates are scaled by the CBL height z i .Vectors indicate the wind direction and magnitude.From van Heerwaarden and Vilà-Guerau de Arellano (2008).

Fig. 13 .
Fig. 13.Mast profiles (squares), balloon data (dots) and LES profiles of downslope velocity (a) and potential temperature (b) in flow over a sloped surface.

Fig. 14 .
Fig. 14.Evolution on time of the vertical concentration (crosswind integrated) of a plume released as a function of the releasing time (dimensional and non-dimensional) for (top) dry convective conditions, (middle) stratocumulus topped boundary layer, and (bottom) shallow cumulus topped boundary layer.Concentration has been multiplied by a factor 1000 to obtain a convenient scale.The crosses indicates the position of the maximum concentration.

Fig. 15 .Fig. 16 .
Fig. 15.Instantaneous vertical cross section of the cloud water q c (g/kg) content and the photostationary state ( ) calculated using the NO, NO 2 and O 3 mixing ratios.At chemical equilibrium = 1.

Fig. 17 .
Fig. 17.The liquid water path for DYCOMS2 using the 5th order advection scheme for various resolutions.Runs are performed without χ * -correction, which in it self significantly enhances the LWP.
we show the time-dependent liquid water path (LWP in the DYCOMS2(Ackerman et al., 2009, DYnamics and  Chemistry Of Marine Stratocumulus experiment) RF02 intercomparison of a stratocumulus-topped boundary layer with a strong inversion ([10]K across a grid cell).It should be noted that the runs www.geosci-model-dev.net/3/415/2010/Geosci.Model Dev., 3, 415-444, 2010 are performed without χ * -correction, which in it self significantly enhances the LWP to values of 110, but are not shown here.

Table 1 .
The main thermodynamical constants used throughout this paper.

Table 2 .
An overview of the parameters used in the SFS scheme of DALES.Not all parameters are independent.
Heus et al.:The Dutch Atmospheric LES