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

Theoretical and Computational Fluid Dynamics, Mar 2013

Large-eddy simulation (LES) has been extensively used as a tool to understand how various processes contribute to the dynamics of the stratocumulus layer. These studies are complicated by the fact that many processes are tied to the dynamics of the stably stratified interface that caps the stratocumulus layer, and which is inadequately resolved by LES. Recent direct numerical simulations (DNS) of isobaric mixing due to buoyancy reversal in a cloud-top mixing layer show that molecular effects are in some instances important in setting the cloud-top entrainment rate, which in turn influences the global development of the layer. This suggests that traditional LES are fundamentally incapable of representing cloud-top processes that depend on buoyancy reversal and that numerical artefacts can affect significantly the results. In this study, we investigate a central aspect of this issue by developing a test case that embodies important features of the buoyancy-reversing cloud-top layer. So doing facilitates a one-to-one comparison of the numerical algorithms typical of LES and DNS codes in a well-established case. We focus on the numerical effects only by switching off the subgrid-scale model in the LES code and using instead a molecular viscosity. We systematically refine the numerical grid and quantify numerical errors, validate convergence and assess computational efficiency of the low-order LES code compared to the high-order DNS. We show that the high-order scheme solves the cloud-top problem computationally more efficiently. On that basis, we suggest that the use of higher-order schemes might be more attractive than further increasing resolution to improve the representation of stratocumulus in LES.

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

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

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

Eckhard Dietze Juan Pedro Mellado Bjorn Stevens Heiko Schmidt Large-eddy simulation (LES) has been extensively used as a tool to understand how various processes contribute to the dynamics of the stratocumulus layer. These studies are complicated by the fact that many processes are tied to the dynamics of the stably stratified interface that caps the stratocumulus layer, and which is inadequately resolved by LES. Recent direct numerical simulations (DNS) of isobaric mixing due to buoyancy reversal in a cloud-top mixing layer show that molecular effects are in some instances important in setting the cloud-top entrainment rate, which in turn influences the global development of the layer. This suggests that traditional LES are fundamentally incapable of representing cloud-top processes that depend on buoyancy reversal and that numerical artefacts can affect significantly the results. In this study, we investigate a central aspect of this issue by developing a test case that embodies important features of the buoyancy-reversing cloud-top layer. So doing facilitates a one-to-one comparison of the numerical algorithms typical of LES and DNS codes in a well-established case. We focus on the numerical effects only by switching off the subgridscale model in the LES code and using instead a molecular viscosity. We systematically refine the numerical grid and quantify numerical errors, validate convergence and assess computational efficiency of the low-order LES code compared to the high-order DNS. We show that the high-order scheme solves the cloud-top problem computationally more efficiently. On that basis, we suggest that the use of higher-order schemes might be more attractive than further increasing resolution to improve the representation of stratocumulus in LES. Marine boundary layer clouds, of which marine stratocumulus is an important archetype, have been identified as one of the primary sources for the uncertainty in model-based estimates of the equilibrium climate sensitivity. Thus, an improved representation of these clouds is thought to substantially increase the reliability of climate change predictions [1]. One characteristic feature of marine stratocumulus that challenges our ability to model them is the strongly stratified layer at the cloud top. This layer has a scale of metres over which the temperature gradient can be positive (keeping in mind the temperature gradient for neutral stratification is 10 K km1) giving rise to the terminology inversion layer, as the temperature gradient is inverted relative to - the tropospheric norm. Already the mixed-layer theory by Lilly [5] identifies cloud-top entrainment, associated with the mixing of air across this region of strong stratification, as an important determinant for the evolution of the cloud. The main processes that drive entrainment are (1) radiative cooling, (2) shear due to large-scale boundary layer vortices and (3) evaporative cooling. Above the cloud top of subtropical marine stratocumulus, the environment is relatively warm and dry, as it is characteristic of the large-scale environment of subsiding air in which the stratocumulus-topped boundary layer (STBL) prevails. The STBL itself is more characteristic of the underlying ocean surface in that it is much cooler and moister, and as the name indicates, cloud-topped. On the one hand, warmer air that is mixed into a cloud-top parcel raises its temperature. On the other hand, because the warmer air is drier, evaporation takes place which cools the mixture. Both effects compete and if the latter dominates, mixtures can develop that are heavier than either one of the two unmixed states. This effect is known as buoyancy reversal [16,23] which in this case is caused by evaporative cooling. In past studies of this problem, the influence of effects at the molecular scale and below has largely been ignored; these include droplet dynamics, phase changes and molecular mixing. The question whether buoyancy reversal alone can cause a runaway instability by positive feedback of evaporative cooling, mixing, and entrainment leading to cloud break-up [2,14] has long been debated. The strength of the reversal is typically expressed by the buoyancy reversal parameter D = bs /b1 [16] that relates minimum and maximum buoyancy bs and b1, respectively. This parameter is typically small. For instance, in the first research flight (RF01) of the DYCOMS II field study [18], a buoyancy reversal parameter of D = 0.031 was measured. Large-eddy simulation (LES) has been extensively used as a tool to improve understanding of the processes at the cloud top. The fidelity of LES has increased tremendously over the last decades. Spatial resolution increased from 50 m in the early studies [2] to 5 m and less in recent LES [11,19,22,24]. Unfortunately, general circulation models and even LES currently have, and in the foreseeable future will have, insufficient grid resolution to resolve the cloud top and the small-scale mixing processes that give rise to buoyancy reversal. Even at the relatively high-resolution characteristics of present studies, the evolution of stratocumulus remains sensitive to the parameterization of subgrid-scale (SGS) processes, or the numerical representation thereof. In some cases, there are order one errors in important properties of simulated stratocumulus such as the entrainment rate. By making an effort to limit mixing, for instance by switching the SGS terms off and using numerical limiters as a closure model, the results improve in some cases [19] but for unphysical reasons. This is unsatisfactory, since the good comparison depends among other on the grid size and the order of the numerical method. This suggests that (1) SGS processes play a significant role in the evolution of the cloud and that (2) numerical errors are of leading order, thus making it hard to distinguish between numerical and physical effects. This paper studies this second aspect of the problem, that is, low-order numerical effects. High-resolution methods, which allow for a representation of molecular effects, have been used to improve the understanding of the entrainment mechanism. Examples include one-dimensional stochastic models such as the linear eddy model [3] and the one-dimensional turbulence model [23] but also direct numerical simulations (DNS), based on the Boussinesq approximation and a continuous two-fluid formulation [10]. Mellado et al. [8] show by means of two-dimensional DNS that complex flow patterns can be generated by weak buoyancy reversal typical for stratocumulus. The DNS also shows that the turbulent mixing dominates the downward development of the cloudy layer and that there is no significant enhancement of turbulent entrainment of upper fluid in the shear-free case. Furthermore, linear stability analysis identifies two time scales: the time scale associated with the downdraughts of heavy fluid due to buoyancy reversal and the time scale of the restoring force of the inversion. The ratio between both is D, and the conjecture is that for typical conditions (0 < D 1) buoyancy reversal alone cannot dissipate the cloud because the inversion returns to equilibrium much faster than the time required by downdraughts to descend deep enough into the cloud layer. Three-dimensional DNS [7] confirm this conjecture. In the simulation, buoyancy reversal promotes convection leading to turbulence with an intensity that is too weak to overturn the inversion. The system evolves into a self-preserving state that retains the vertically layered structure. Above the convection layer, a thin diffusion layer with a constant thickness h = /we forms and travels upwards at a constant mean speed we (|bs c2)1/3, where is the molecular diffusivity and c is the crossover mixture fraction, equal to the interval of negatively buoyant mixtures. Analysis of the budgets of turbulent kinetic energy shows that a strong conversion of vertical to horizontal momentum occurs in a thin region between the inversion layer and the convection layer [7,10]. Probability density functions of velocity in the turbulent region show a strong vertical variability and often non-Gaussian behaviour, even in the horizontal component [9]. This indicates strong anisotropy and inhomogeneity of the turbulence which, in addition to the relevance of molecular transport, further complicates turbulence modelling in LES for this kind of flows. Despite the obvious advantages of using DNS, bounded computational resources mean that LES is still required for global, complex geophysical systems like the STBL. As a step towards better LES of stratocumulus, we study the behaviour of a typical LES numerical algorithm for a test problem that we have defined to capture the processes thought to be relevant to buoyancy reversal and its contribution to cloud-top entrainment. We choose the two-dimensional cloud-top mixing problem by Mellado et al. [8] and use their DNS solution as a reference. This particular case is explored because it embodies a number of important features such as a stably oscillating inversion, molecular diffusion, baroclinic production of vorticity, and evaporative cooling leading to buoyancy reversal, which merit the introduction of a new test problem for small-scale atmospheric flow solvers. We use the University of California, Los Angeles LES (UCLA-LES) [1921] where we switched off the SGS model and introduced a constant molecular viscosity. It is emphasized that, in this paper, the term LES will be used solely to refer to the UCLA-LES numerical algorithm. Where, in contrast, LES modelling is meant it will be clear from the context. The reference DNS uses a sixth-order compact scheme, whereas in LES low-order methods are commonly used, typically second-order as in this study. The emphasis on low-order effects on LES scales arises because the flow is not smooth on the grid scale, thereby obviating the advantages of higher-order methods. Low-order methods also use a smaller stencil and thus are more easily to implement efficiently in a highly parallel environment. What we are interested in here is how sensitive solutions are to numerical errors and what the resolution requirements are as compared to a high-order DNS solver. For that, we carry out a grid convergence study where we systematically refine the computational grid until convergence is achieved. By doing that, we also verify the UCLA-LES code for the test case. Because the choice of highversus low-order methods often involves trade-offs, we also investigate how the computational costs of both codes compare. The remainder of the paper is organized as follows. In Sect. 2, the formulation of the test case is described. In Sect. 3, we briefly describe the numerical algorithms of both the reference DNS and the LES code and present results of the convergence study showing self-convergence of the LES code and convergence of the LES code to the DNS solution. Two lengths, hb and ht , are introduced, which quantitatively characterize the evolution of the flow, and are used to quantify convergence. Also, we address computational efficiency by comparing the computational cost of the DNS and the LES code packages as they are normally used. Section 4 concludes the present contribution by summarizing the main findings and discussing the implications of this study for LES of stratocumulus in general. 2 Formulation 2.1 Reference problem We study a two-layer system in which a cold moist layer, which is saturated and has liquid water, is topped by a relatively warmer, drier layer. In the mixing region between the two layers, the total water mixing ratio qt and specific enthalpy h transition between the two values characteristic of the upper unsaturated and lower saturated layer, respectively. Initially, the system is at rest, gravity acting downwards. The system so defined is a simplified surrogate of the stratocumulus top in order to investigate certain physical phenomena associated with relatively small-scale dynamics. We describe the multiphase flow using a two-fluid formulation. The main underlying hypotheses are the following: (1) The liquid water can be considered a continuous phase, (2) the mixture always stays in local thermodynamic equilibrium and (3) the diffusivity of the liquid phase equals that of dry air and water vapour. These assumptions are not identically satisfied under real conditions in stratocumulus but provide a good approximation to study latent heat effects. A more thorough discussion of these hypotheses and their applicability can be found in [10]. Assuming equal thermal and mass diffusivities, the thermodynamic coordinates h and qt obey the same advection-diffusion equation. This implies that, with appropriate boundary conditions, both quantities have the same normalized solution after an initial transient. Thus, their evolution can be described by that of a single conserved scalarthe mixture fraction . With the indices 0 and 1 indicating the states in the lower and upper layer, respectively, the mixture fraction can be expressed in terms of h and qt as h(x, t ) h0 h1 h0 qt (x, t ) qt,0 qt,1 qt,0 Herein, x = (x , z)T denotes the location and t is the time coordinate. Hence, the mixture fraction simply measures the fraction of mass of fluid from the top layer in an arbitrary fluid parcel in the flow; by definition, /b1 0.4 b Fig. 1 Non-dimensional buoyancy mixing function for the present case (solid), and two fictional cases: a non-reversing (dashed) and a strongly reversing case (dotted). The latter one corresponds to the case A3 by Mellado et al. [8] = 0 in the bottom layer and = 1 in the top layer. The thermodynamic state of the system, in particular the density field, is then completely defined by the field of . We use the Boussinesq approximation and describe gravitational forces in terms of the buoyancy g denoting the magnitude of the gravitational acceleration and is the density. This is justified because differences in density across the cloud-top inversion are small. Depending on the initial data, evaporative cooling can cause saturated mixtures to be heavier than either one of the unmixed reference states. The term buoyancy reversal is introduced to identify this eventuality [16,23]. The heaviest mixture 0 s 1 occurs at saturation conditions and has buoyancy bs = 00s g. Shy and Breidenthal [15] introduce the buoyancy reversal parameter b(x, t ) = to describe the strength of the reversal, relating buoyancy of the heaviest mixture to the buoyancy jump b1 = 001 g across the inversion. With the above assumptions and for small density differences 0 1, buoyancy can be accurately approximated by a piecewise linear function of the mixture fraction. Here, the formulation introduced by Mellado et al. [8] is adopted. Figure 1 shows three examples of the graph of b( ): namely a non-reversing and a strongly reversing case for illustrative purposes, and the buoyancy mixing function for the present case. The parameter s was introduced to smooth the discontinuity between the piecewise linear sections. Obviously, the smoothing damps the maximum negative buoyancy bs and, thus, the effective reversal parameter D. Previous work [8] shows that a smoothing parameter of s = s /16 leads to a sufficient independence of the solution from changes in s . This value of s is adopted here for consistency in the comparison. 2.2 Model equations With the assumptions made in Sect. 2.1, the set of prognostic variables in two dimensions consists of the two Cartesian velocity components u = (u, w)T, the kinematic pressure p and the mixture fraction . Fig. 2 Vertical profiles of the mixture fraction and buoyancy as following from Eqs. (3) and (8) The reference DNS solves the Boussinesq limit of the Navier-Stokes equations, given by u = 0. The LES code, whose native formulation solves Ogura-Phillips [12] type anelastic equations, was modified to solve exactly the same system. In Eq. (4), is the constant kinematic viscosity of the mixture. The buoyancy mixing function b( ) is given by Eq. (3), k is the vertical unit vector, and p is a modified pressure, normalized by the density. The Poisson equation for the dynamic pressure couples conservation of mass, see Eq. (6), to the momentum Eq. (4) and completes the system of equations. 2.3 Initial and boundary conditions We consider a density-stratified system that is at rest in the beginning. Profiles of specific enthalpy and total water mixing ratio change in the vertical from values characteristic of the saturated and unsaturated layer over a height of the order of 4 following an error function. Employing the mixture fraction formulation, Eq. (1) allows for a specification of the initial thermodynamic conditions by 1 (t = 0, x , z) = 2 z z0(x ) Herein, z0 denotes the vertical position of maximum stratification that corresponds to the height where = 0.5. The combination of Eqs. (3) and (8) gives the vertical profile of the buoyancy. Both profiles are plotted in Fig. 2. In a non-reversing scenario (cf. dashed graph of Fig. 1), the normalized buoyancy would match the mixture fraction profile. In contrast, under conditions of buoyancy reversal, a small peak to the negative buoyancy axis develops just below the = 0.5 isoline. Thus, the system is unstable with respect to the lower fluid in a thin layer. The instability process is triggered by imposing one mode of a sinusoidal wave on the inversion height z0 with the wave length and amplitude a/2 according to a z0(x ) = 2 cos Fig. 3 Time series of the buoyancy field at times t / T = (0.0, 1.05, 2.18, 3.32, 4.45) with T = 2 /b1, colours range from zero (white) to b1 (black). UCLA-LES, 512 1,024 grid cells Here, a measures the thickness of the perturbed mixing layer. The ratio of amplitude to wave length is (a/2)/ = 0.1, and the initial thickness of the error function profile is / = 0.025. We consider a twodimensional rectangular domain with the extent ( x , z ) = (, 2). The mean height of the inversion, given by the mean height of the = 0.5 isoline, is located at 60 % of the vertical domain. The leftmost plot in Fig. 3 shows the initial situation. We consider conditions that correspond to the STBL as observed in the first research flight of the second DYCOMS field programme [18]. The data show weak buoyancy reversal with a buoyancy reversal parameter of D = 0.031. The mixture fraction at saturation conditions is s = 0.09, which completes the description of the buoyancy function (3) (cf. Fig. 1). Mellado et al. [8] introduce the reference Grashof number Gr = which can be interpreted as the square of the ratio of two characteristic time scales. The first one is the viscous time that scales as a2/, the second one is the period of the oscillation of the interface and is proportional to /b1 [8]. For the reference simulation, the kinematic viscosity = 1.5 105 m2/s and Gr = 6.4 105 were used. That corresponds to, keeping the above ratio of (a/2)/ in mind, a wave length 0.71 m of the sinusoidal perturbation. We impose cyclic boundary conditions at the lateral and free-slip conditions at the top and bottom boundaries. For scalars, zero gradients in the vertical are set at the top and bottom. The simulations are stopped before finite-size effects associated with the computational domain become important. 3 Simulations 3.1 The numerical methods We use the cloud-top problem defined in Sect. 2 as a test case to investigate the numerical algorithm used in the UCLA-LES code. As a reference, we use the accurate numerical solution obtained with a DNS code and already employed by some of the authors in previous work [8]. We emphasize again that the focus in this work is the analysis of low-order numerical effects, and not the subgrid-scale models intrinsic to LES (cf. [13]). For this purpose, the UCLA-LES code is modified so that it solves Eqs. (3)(6), exactly the same equations solved with the reference DNS code. Practically, this is achieved by substituting the subgrid-viscosity in the UCLA-LES code with a constant molecular viscosity . The UCLA-LES model discretizes the system of Eqs. (4), (5), and (7) in a finite-volume manner on a static Arakawa-C grid (regular Cartesian), where velocities are staggered half a grid point from the rest of the (cell-centred) variables ( , p). The grid spacing is uniform in both the horizontal and the vertical direction. Second-order spatial discretization is used, scalar fluxes are limited using the monotonized central flux limiter. The flow is advanced in time using a third-order Runge-Kutta scheme. Throughout all simulations, Courant numbers of CFLadv = 0.5 were maintained. The UCLA-LES requires periodic boundary conditions in the lateral as it uses fast Fourier transform in the horizontal direction to transform the Poisson Eq. (7) to a second-order ordinary differential equation (ODE). The ODE is solved using a tridiagonal solver. The reference DNS uses finite differences and a regular Cartesian mesh. Spectral-like sixth-order compact Pad schemes are used to discretize the first- and second-order spatial derivatives (closed at the top and bottom with third-order biased formulae, which leads to a global fourth-order accuracy in space). A low-storage 5-stage fourth-order Runge-Kutta scheme is employed for the time advancement. This scheme has a stability region larger than that of the third-order counterpart, and the Courant numbers used in the DNS are CFLadv < 1.2 and CFLdiff < 1.2/4. Similarly to the UCLA-LES, the solenoidal constraint is imposed using Fourier decomposition along the periodic horizontal planes in order to reduce it to a set of one-dimensional second-order equations which are solved using those same compact schemes. The boundary conditions imposed at the top and the bottom are the same as those in the UCLA-LES code. Further details of the numerical algorithm can be found in Mellado et al. [10]. The major difference between the DNS and the UCLA-LES code is the resolving efficiency. This resolving efficiency is normally quantified in terms of the number of points per wavelength (PPW) required to maintain errors in the corresponding transfer function below a specified level (or, equivalently, a specified error in the dispersion velocity of the linear advection problem). In these terms, an error of 1 % requires 4 PPW in the case of the implicit compact scheme used in the DNS, whereas the second-order explicit scheme requires 25 PPW [4,6] (this difference is even higher if the common reference error of 0.1 % is retained, for which the previous finite difference schemes require 6 PPW and 100 PPW, respectively). However, these are theoretical values based on linear analysis of the algorithm: in this work, we ascertain the relative performance of these two numerical schemes in an applied case very relevant to the physics of the STBL, namely the stratocumulus top. 3.2 Analysis of the flow The qualitative evolution of the flow is illustrated in Fig. 3 through successive snapshots that visualize the negative buoyancy field. The data shown were obtained from the UCLA-LES run on the reference grid with 512 1, 024 cells. The leftmost panel corresponds to the sinusoidal initial condition, as previously discussed in Sect. 2.3. The stratified system is initially at rest. In the vicinity of the inversion, net vertical accelerations arise from horizontal buoyancy differences. A persistent oscillation of the interface sets in that, at high Gr , is only slowly damped by molecular viscosity. At the same time, a spike or finger starts to form of the thin negatively buoyant sheet at the lowest point of the oscillation. The spike grows downwards and later rolls up to a vortex pair. Molecular diffusion mixes fluid of the falling finger with the ambient moist fluid and, thus, spreads buoyancy reversal. According to the buoyancy mixing function (cf. Fig. 1), moist fluid ( = 0) that is mixed with negatively buoyant fluid is negatively buoyant itself. This process continues as the pulsating motion of the inversion pumps more negatively buoyant fluid downwards. The instability process is driven by the transformation of potential energy (here subsumed within the enthalpy of vaporization) to kinetic energy. As the unsaturated fluid in the upper layer mixes with the saturated fluid in the lower layer, more negatively buoyant fluid is produced due to evaporative cooling, and part of it is ejected downwards. This gradually leads to a vertical growth of both, the size of the saturated bottom layer and the depth of the mixing region. Mellado et al. [8] quantify the instability process by means of two characteristic lengths that measure the thickness of the mixing layer to the top and to the bottom. These are the penetration length hb of the falling finger and the perturbation thickness ht that measures the upwards growth of the moist layer. Both are defined with the help of the horizontally averaged mixture fraction profile (z) as shown in Fig. 4. Taking the mean initial height of the = 0.5 isoline, z0, as a reference, both lengths are defined as distances to the vertical Fig. 4 Definition of the penetration length hb and the upper perturbation thickness ht using the horizontally averaged mixture fraction profile . The dashed lines indicate the vertical position where departs from a given threshold (103) of the upper and lower state, respectively S ,DN0.3 b h /hb0.2 Fig. 5 History of the penetration length (left) and the relative error with respect to the reference DNS (right). Grid resolutions are given with respect to the reference grid with 512 1,024 grid points position where the mean mixture fraction departs from a given threshold of the value of the respective layer. For hb, the first position looking from the bottom, for ht the first position looking from the top are detected, respectively. Because mixing in the lower layer is limited to small fractions of below 0.1, a threshold of 103 is chosen. The temporal evolution of these lengths is presented in Figs. 5 and 6, respectively, where the solid lines indicate the results from the DNS by Mellado et al. [8]. The evolution of the penetration length hb observed in Fig. 5 exhibits a superposition of the stable oscillatory mode of the inversion and the unstable mode of the falling finger. If buoyancy reversal was not present (cf. dashed graph in Fig. 1), only the stable mode would be visible. The evolution of the upper perturbation thickness ht also has the stable mode; it is superimposed by a quasi-linear growth that results from molecular mixing of at the top of the inversion and the baroclinic production of vorticity. The transport of mixed fluid away from the inversion due to buoyancy reversal steepens local scalar gradients which in return enhances molecular mixing. In the following, we use these two lengths to quantify numerical errors and the convergence of the UCLA-LES code. 3.3 Convergence behaviour In order to study the convergence behaviour of the UCLA-LES code, a grid convergence study was carried out. The reference DNS employs a grid with 512 1,024 grid points; this will be referred to as reference resolution or reference mesh in the following. Here, simulations were carried out on six different meshes, ranging from 1/8th to 4 times the reference resolution (see Table 1). Fig. 6 History of the upper perturbation thickness (left) and the relative error with respect to the reference DNS (right). Straight lines (left) are least squares regressions of the growth rate from the 2nd to the 9th local maximum. Grid resolutions are given with respect to the reference grid with 512 1,024 grid points Table 1 Order of convergence p of the UCLA-LES according to the L2 and maximum norm L The factors in the mesh row refer to the number of mesh cells, for example, 1/8 refers to an eight times coarser mesh The UCLA-LES code, used as explained in Sect. 3.1, has an anticipated second-order convergence in space. To confirm this, simulations on all meshes have been performed with a common time increment, and the order of convergence p in terms of the L 2 norm and the maximum norm L of the mixture fraction field have been determined. Based on the solution of on three meshes resolved by x = z = h, 2h, and 4h, respectively, the norms S h,tND0.3 /ht 0.2 L 2(h) = i=1 are computed, where xi = (xi , zi ) is the position of the i th cell-centred point of the mesh. Linear interpolation on the staggered mesh is used to interpolate values from the finer mesh to the coarser mesh. The order of convergence is computed according to p2, = ln L 2,(2h) L 2,(h) This procedure is repeated for all subsets of three neighbouring meshes. The results are presented in Table 1 and confirm the anticipated second-order convergence in both, the maximum and L 2 norm. In the following, we want to quantify the absolute errors of the solution. We do this by comparing the evolution of hb and ht in the DNS and LES code. Besides the reference DNS solution, Figs. 5 and 6 show the results obtained from the LES on the six meshes. For both hb and ht , a clear trend is obvious: with coarsening the mesh spacing the growth of both quantities is increased. This can mainly be attributed to increased numerical mixing, yielding more generation of negatively buoyant fluid that enhances the convective mixing process in the cloud layer. The temporal evolution of the relative errors is explicitly showed in the second panel of Figs. 5 and 6. On the coarsest grid, relative errors of hb are up to 30 % and decrease to 5 % at reference resolution and to 2.4 % on the finest grid. Relative errors of ht tend to be greater on the coarsest grid but vanish more rapidly as the grid is refined. On the finest grid, that is, four times the number of grid points in each direction than in the reference DNS, the relative error remains below 1 % throughout the simulation. 1/8 1/4 1/2 [512 1,024] 2 4 Table 2 Relative errors of the ht growth rate ELES in the low-order simulations with respect to the DNS solution Fig. 7 Amount of fluid entrained in the section of the domain below the lowest point of the initial wave of the = 0.5 isoline. It is measured by the integral I of the mean mixture fraction profile normalized by as in Eq. (8) The mean vertical growth of the upper perturbation thickness ht is an effect that resembles the cloud-top entrainment in real stratocumulus. We quantify this mean growth E with help of a linear regression based on the data points after the first period until the ninth period of each individual data set. The latter corresponds to the last ht maximum in the DNS reference data. The regression lines for the DNS solution and the finest and coarsest LES solution are shown in Fig. 6. Relative errors with respect to the DNS solution are presented in Table 2. While at four times the reference resolution the error in the mean vertical growth of the upper perturbation thickness reduces to 1.4 %, it increases to 89 % on the coarsest grid with 1/8th the reference resolution. The latter corresponds to a mesh spacing of about 1 cm, still much finer than in current high-resolution LES of stratocumulus, where typically the vertical resolution is between 1 and 5 m at the cloud top [19, 24]. This result shows clearly that state-of-the-art LES does not resolve the sharp gradients at the cloud top and that this under-resolution has important effects on the thermodynamics of the system and helps explain the difficulties encountered in understanding local processes of the cloud top [17]. In Fig. 7, we show the total amount of fluid entrained from the upper layer into the section of the domain below the lowest point of the initial wave. This portion of fluid is simply measured by the (vertical) integral I of the mean mixture fraction profile in that subdomain and that portion of fluid monotonously increases when the grid is coarsened. This indicates that the more rapid increase in hb on the coarser grids is, in fact, associated with a larger amount of fluid from the upper layer carried downwards by the falling finger. It is not solely attributed to the stronger spurious (numerical) spread of the border of the falling finger, but also to a modification of the global dynamics of the system. Here, numerical errors act as to increase the entrainment of drier fluid into the cloudy layer. 3.4 Computational efficiency Both the DNS and UCLA-LES algorithms were originally implemented to solve three-dimensional problems. For the work discussed so far, the codes have been run on a two-dimensional grid in order to reduce computational cost and focus on some details of the cloud top in a very controlled manner. However, that causes computational overheads that depend on the respective methods and their implementations. To enable for a fair comparison in terms of computational efficiency, three-dimensional runs were carried out. We extended DNS results are based on a low-storage fourth-order Runge-Kutta scheme and UCLA-LES results on a third-order Runge-Kutta scheme the 128 256 grid (1/4th of the reference resolution) in the lateral to 128 128 256. Both codes were compiled with the same compiler and run on one core of an AMD Opteron model 2384, 2.7 GHz of a Sun X2200 M2 machine. In Table 3, we compare the computational time per iteration, T1, and the computational time T2 needed per characteristic time of the system /b1 (averages over 100 iterations were considered in order to eliminate I/O times and other overheads). The DNS is more expensive per iteration, as expected because of the fourth-order Runge-Kutta scheme (which uses 5 stages because of the two-level storage property) and an implicit spatial discretization. When the DNS code is run with a third-order Runge-Kutta scheme, similarly to the LES code, the time per iteration is 22 seconds. This is about 40 % more expensive than the low-order code; this 40 % difference is due purely to spatial discretization and accounts for the overhead due to the solution of the linear systems being associated with the implicit compact schemes, in contrast to the explicit character of the second-order central schemes. However, it is arguably more interesting to consider the computational cost per characteristic physical time of the system, that is, how much it costs to advance the transport equations during a given time interval. This is the second parameter shown in Table 3, the time T2. The difference with the measurement per iteration T1 discussed in the previous paragraph is due to the fact that the low-storage fourth-order Runge-Kutta has a stability region significantly larger than that of the third-order Runge-Kutta and allows a larger CFL number. Comparisons of the values of T2 between the DNS and the UCLA-LES codes show that the DNS is still slower, but only by 10 %. This last result is very important. As we have shown in Sect. 3.3, using the low-order numerical method of the UCLA-LES code, higher resolution is needed to converge to the reference solution. For the particular problem we are studying, the buoyancy reversal instability at the cloud top, using the low-order algorithm leads to errors of the order of 20 % and is only about 10 % faster. Let us imagine we desire to reduce errors using the UCLA-LES code to about 10 %. Then, according to Table 2, we need to increase the resolution by about a factor of two. This step already implies a higher computational cost; a reasonable rough estimate would be a factor of 8 in three-dimensions, or a factor of 4 in two dimensions, assuming a linear relation between the operation count and the number of grid points and neglecting deviations due to parts of the algorithm which do not scale linearly with the number of points, like fast Fourier transforms or cache-size management. On top of that, the time increment t imposed by the explicit scheme, so stability is satisfied, multiplies further this factor by 2 or 4, according to Eqs. (9) or (10), respectively. In our particular problem, it has been checked that the latter dominates because the diffusion constraint Eq. (10) becomes the most restrictive and t scales as ( x )2. Altogether, this analysis means about a factor of 10 or more increase in computational cost using a loworder algorithm, compared to a factor 1.1 increase using the implicit spectral-like compact finite differences, to obtain a solution within a 10 % error. It is appropriate to note that the discussion of the last paragraph refers to a serial implementation of the algorithm. The advantage of low-order spatial schemes is that they are explicit and parallelizations thereof are known to scale almost linearly, so that the computational cost remains approximately independent of the number of processors, whereas the implicit schemes suffer a strong overhead from the communication among a relatively large number of processors. Hence, there is a trade-off and a crossover number of processors beyond which a low-order scheme is more efficient than a high-order implicit scheme. However, this crossover value is a relatively high number. For instance, the DNS algorithm used here as a reference is run typically using 1,024 or 2,048 MPI tasks (e.g. simulations presented in [10]) in the German Climate Computer Centre. In those architectures, the time spent in communication among processors varies between 40 and 60 %, meaning an overhead of the order of a factor of 2 with respect to the reference computational cost. The same behaviour has been observed in scaling studies of the DNS code up to 8,192 tasks in the Jlich Supercomputer Centre. This overhead is still significantly smaller than the factor of 10 estimated above for a second-order scheme used in the UCLA-LES code to obtain solutions within 10 %. Last but not least, we need to consider as well the concomitant additional memory requirements using low-order schemes because of the larger number of points (this can actually be the most critical aspect in large-scale three-dimensional simulations). Hence, as a summary, high-order schemes, in spite of the higher complexity in implementation, might be attractive for LES of stratocumulus, and other high Reynolds number flows if small-scale processes might play a relevant role or simply to better represent the smooth portions of the flow, such as waves, that might exist at some multiple of the grid scale. This has implications for the current trend to increase resolution in LES models and the application of adaptive solvers to such problems. These conclusions hold also for parallel implementation, at least up to the order of 10,000 processors, which current architectures provide. 4 Conclusions We presented in detail the two-dimensional evaporatively driven cloud-top mixing layer as a test case for the study of low-order numerical effects. We did that comparing the UCLA-LES code, which uses an explicit second-order scheme, against a DNS solution, based on a spectral-like compact sixth-order accurate method. The subgrid-scale turbulence model of the UCLA-LES code was switched off, and a constant molecular viscosity was used instead so that only numerical effects were addressed in this work. We presented a grid convergence study where we investigate the sensitivity of the solution to numerical errors and the resolution requirements for the UCLA-LES algorithm to achieve convergence. Analysis of the solutions on the various grids showed that the LES numerical algorithm is self-convergent meeting the anticipated effective second-order rate of convergence. In order to quantify numerical errors, we compared the evolution of the penetration length hb of the downdraught and the upper perturbation thickness ht between the LES and DNS solutions. The coarsest grid had a vertical spacing of approximately 1 cm which is at least two orders of magnitude finer than in typical LES. For these quantities, 30 and 40 % relative errors, respectively, were measured on the coarsest mesh; almost 100 % relative errors in the entrainment rate. Hence, we have shown that very high-resolution or high-order schemes are needed to capture the small-scale molecular mixing. If these small-scale effects are not accurately represented, important large-scale characteristics, such as the vertical growth of the inversion height, resembling cloud-top entrainment, may contain order one defects caused by numerical artefacts. The analysis of the computational efficiency showed that the cloud-top problem is solved more efficiently (in terms of total CPU time) by the high-order scheme used in the reference DNS than by the low-order scheme of the UCLA-LES code. Based on this finding, we conclude that increasing the order of numerical methods might be better suited than further increasing the resolution of low-order methods to improve the LES of flows that are not smooth on the grid scale (such as LES of stratocumulus). The present study suggests that the larger computational cost per integration step is by far compensated by a smaller number of grid points, a larger possible temporal stride, and reduced magnitude of numerical errors. The difference is so remarkable that the same conclusion holds in a parallel implementation of the algorithms, at least up to the order of 1-10 thousand processors, despite the communication overhead associated with the implicit schemes. Acknowledgments Financial support for this study was provided by the Deutsche Forschungsgemeinschaft within the SPP 1276 MetStrm and the Max Planck Society for the Advancement of Science. Computational resources have been provided by the German Climate Computing Centre, Hamburg. Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.

This is a preview of a remote PDF:

Eckhard Dietze, Juan Pedro Mellado, Bjorn Stevens, Heiko Schmidt. Study of low-order numerical effects in the two-dimensional cloud-top mixing layer, Theoretical and Computational Fluid Dynamics, 2013, 239-251, DOI: 10.1007/s00162-012-0263-0