Characteristics of vertical velocity in marine stratocumulus: comparison of large eddy simulations with observations

We simulated a marine stratus deck sampled during the Marine Stratus/Stratocumulus Experiment (MASE) with a three-dimensional large eddy simulation (LES) model at different model resolutions. Various characteristics of the vertical velocity from the model simulations were evaluated against those derived from the corresponding aircraft in situ observations, focusing on standard deviation, skewness, kurtosis, probability density function (PDF), power spectrum, and structure function. Our results show that although the LES model captures reasonably well the lower-order moments (e.g., horizontal averages and standard deviations), it fails to simulate many aspects of the higher-order moments, such as kurtosis, especially near cloud base and cloud top. Further investigations of the PDFs, power spectra, and structure functions reveal that compared to the observations, the model generally underestimates relatively strong variations on small scales. The results also suggest that increasing the model resolutions improves the agreements between the model results and the observations in virtually all of the properties that we examined. Furthermore, the results indicate that a vertical grid size <10 m is necessary for accurately simulating even the standard-deviation profile, posing new challenges to computer resources.


Introduction
Aerosol indirect effect (AIE) on climate probably is the most uncertain among known climatic forcings (Intergovernmental Panel on Climate Change 2007); furthermore, its forward and 'inverse' estimates also demonstrate large discrepancies (Anderson et al 2003). An uncertainty analysis pointed out that vertical velocity (W ) was one of the largest uncertainties in the assessment of the first AIE (Chen and Penner 2005). Vertical velocity fields also were highlighted as the most important 3 Address for correspondence: Atmospheric Sciences Division, Building 815E, 75 Rutherford Drive, Brookhaven National Laboratory, Upton, NY, 11973-5000, USA. variable in the parameterization of the AIE in a single-columnmodel inter-comparison study (Menon et al 2003).
Vertical velocity, W , is critical for aerosol activation processes that are the first bridge between aerosols and clouds, and is known to affect both the number concentration (Twomey 1977) and the relative dispersion of cloud droplet size distribution (Liu et al 2006). Updrafts/downdrafts also significantly impact the temporal evolutions of pockets of open cells of various reflective properties (Xue et al 2008), their lifetime, and thus domain average reflectivity and radiative budget consequentially. More importantly, W is known to vary substantially, even in seemingly uniform stratocumulus clouds (Lu et al 2007); such variation could influence the AIE evaluation strongly in both maritime-and continental-clouds (Feingold et al 2003, Leaitch et al 1996. However, most global climate models (GCMs) diagnose W based on largescale convergence/divergence, but their coarse resolution (∼100 km) likely fails to represent subgrid processes that generate small-scale fluctuations, and thereby limiting their ability to characterize W and its variations (Kiehl et al 1998). The high multi-scale variability in W probably is a major source of large uncertainties in the expressions employed to relate cloud properties (e.g., droplet concentrations and relative dispersion) to pre-cloud aerosol properties (Ghan et al 1995, Liu et al 2006, Peng et al 2005. Quantitatively understanding and parameterizing the subgrid variability in W is also the key to improving cloud parameterizations in GCMs in general (Lohmann et al 1999, Intergovernmental Panel on Climate Change 2007, Menon et al 2003.
High-resolution atmospheric models, such as cloudresolving models and large eddy simulation (LES) models, can resolve the spatial-and temporal variations in W much better than GCMs. Considering various observational limitations (such as the instrumental restrictions, sampling bias) and incomplete observational records, researchers suggested using high-resolution simulations to complement such observations (Ackerman et al 2004 in developing the parameterizations of subgrid variability. However, despite the great progress in LES modelling, few studies have been devoted to examining the variability in W in LES against the corresponding observations. This lack is especially true for high-resolution simulations mainly because of limited computer resources. Recognizing the profound significance of W and its variability and the need to evaluate high-resolution models against observations, in this study, we took the advantage of a highly scalable supercomputer (IBM Blue Gene) which is lately available at Brookhaven National Laboratory (BNL) to explore the characteristics of W and its spatial variations within a marine stratus cloud sampled on 19 July, 2005 during the MASE (Marine Stratus/Stratocumulus Experiment, Daum et al 2008). We compared aircraft in situ observations with the LES outputs at different model resolutions. Special focus was on the quantitative measures of variability, and on the agreement between observations and LES outputs of differentorder moments of W (e.g., standard deviation, skewness, and kurtosis), probability density functions (PDF), power spectra, and structure functions of W .

Model set-up and case description
We chose the three-dimensional parallel version of the Goddard cloud ensemble model (GCE) for this study (Juang et al 2007, Tao et al 2003, Tao and Simpson 1993. Turbulence is parameterized via a 1.5-order closure scheme; the microphysical parameterization of clouds is a single-moment bulk scheme; and broadband shortwave-and longwave-radiative transfer modules are interactively coupled with cloud fields. The model domain is 6.4 × 6.4 km 2 horizontally, and 1.25 km vertically. The horizontal ( x) and vertical ( z) grid sizes, respectively, are 25 and 10 m in the base case, denoted as 'X25Z10'. We conducted two sensitivity tests with different vertical grid sizes of 20 and 40 m (denoted, respectively, as 'X25Z20' and 'X25Z40'), mainly because the simulations of this kind of sub-tropical marine boundary layer under strong inversion (during the MASE) are more sensitive to vertical resolution than to horizontal resolution (Guo et al 2008). Our simulations began at 15:30 LST (local standard time) on 18 July, 2005; and the entire simulation period was 32.5 h. (Note: Guo et al (2008) give more detailed descriptions of the model's set-up and the case of interest.) We mainly undertook our analysis from 9:30 to 11:30 LST on 19 July when the aircraft in situ measurements of W were available. During the MASE, the measurements of W were made with a gust probe aboard the Department of Energy G-1 research aircraft. The original sampling frequency was 100 Hz. We utilized the W data at the 10 Hz frequency to reduce measurement noise. The associated observational uncertainty was about 0.07 m s −1 , mainly arising from the interference among different electronic instruments aboard the aircraft. Figure 1 compares the vertical profiles of standard deviation (σ W ), skewness (S W ), and kurtosis (K W ) of W from the aircraft in situ observations with the model outputs at different model resolutions. We note that the horizontal averages of W at different altitudes generally are close to 0 m s −1 , indicating the near cancellation of updrafts and downdrafts over the horizontal domain, and the importance of local fluctuations (such as σ W ) associated with the boundary layer turbulence. According to the airborne observations (the red squares in figure 1(a)), the magnitude of σ W is about a few hundredths of m s −1 , and increases from 0.03 m s −1 at an altitude of 70 m (near the cloud base) to 0.07 m s −1 at an altitude of 360 m (near the cloud top), and then decreases sharply to nearly 0 m s −1 at an altitude of 500 m (above the clouds). These weak vertical motions lead to a small maximum super-saturation rate of about 0.05-0.08% near the cloud base .

Vertical profiles of standard deviation, skewness, and kurtosis
The simulated σ W in 'X25Z10' compares reasonably well with the observed σ W and captures major vertical variations. However, with a coarser resolution, such as in 'X25Z20' and in 'X25Z40', the magnitudes of σ W are underestimated, and this underestimation becomes more significant with decreasing model resolutions. Underestimation is expected because in a lower-resolution simulation, small-scale convection is underresolved, and thus, smeared out by dissipation, so generating weaker updrafts/downdrafts and smaller fluctuations (σ W ). It is noteworthy that this underestimation is more significant near the cloud base than that at the middle of the cloud or near the cloud top. This might be traceable to the following two reasons. First, the cloud base of this stratus is low and close to the sea surface. The simulated cloud properties near the cloud base would, therefore, be more or less influenced by lower boundary conditions, where W is set to be 0 m s −1 . Second, the large-scale forcing is of profound importance for simulation results. Since we used nudging terms to emulate large-scale forcing (Guo et al 2008), we conducted another sensitivity test using a stronger nudging forcing. It turned out that there was a better agreement with observations, and the relative difference in (σ W ) can be reduced to about 30% with the most refined resolution of x = y = 25 m and z = 20 m near the cloud base (not shown here). Hence, large-scale forcing is of importance, and more realistic forcing data are needed to achieve better agreement with observations.
The observed values of skewness (S W ) generally are negative throughout the cloud layer ( figure 1(b)). It was demonstrated that a negative S W often is associated with downdraft cores surrounded by a large area of weak updrafts (Moeng andRotunno 1990, Moyer andYoung 1991). Accordingly, the dominance of the negative S W suggests that cloud top longwave-radiative cooling is an important driver for this marine stratus cloud. Near the cloud base, observations suggest that the S W is slightly positive, signifying positive buoyancy contributions from surface sensible and latent heat fluxes. The S W then drops to −0.4 at the middle of the cloud at an altitude of 0.28 km, and rises to about 0 above the cloud layer (the red squares in figure 1(b)). From the model outputs, the S W is consistently negative and tends to increase from the cloud base to the cloud top (the black curves in figure 1(b)). This monotonic increase of S W with height does not conform to that from the observations, but is compatible with other LES results (Moeng and Rotunno 1990). The discrepancy between the model outputs and the observations may be associated with the upper/lower boundary conditions (used in the model) that change the shapes of eddies as eddies hit the boundary, modify their circulations, and hence, impact the S W (Moeng and Rotunno 1990). On the whole, the S W from 'X25Z10' seems to agree with that deduced in the observations from the lower to upper cloud layers, although both 'X25Z20' and 'X25Z40' underestimate the magnitude of the S W .
The value of kurtosis (K W ) of the observed W is larger than that of the modelled W in 'X25Z10', 'X25Z20', and 'X25Z40', indicating that the observed W has stronger local fluctuations (or a 'fatter' tail) than the modelled one, especially near the cloud base and the cloud top. Within the cloud layer, the difference between the observations and model outputs is relatively small, but not negligible (figure 1(c)). Figure 2 depicts the relative differences of the model outputs in 'X25Z10', 'X25Z20', and 'X25Z40' with respect to the observations. Clearly, for σ W and S W within the cloud layer, the highest-resolution run ('X25Z10') exhibits the smallest difference (∼30%), while the lowest-resolution run ('X25Z40') shows the largest difference (∼100%) (figures 2(a) and (b)). So, with increasing resolution there is better agreement for the second and third order moments of W . For the fourth-order moment (K W ), the magnitudes of the relative difference all are about 100% for different resolution runs; here the superiority of a higher-resolution simulation of 'X25Z10' is not evident (figure 2(c)). It is noteworthy that subgridscale contributions are not explicitly included in the above analysis, because most of the eddy energy is resolved in LES (Moeng et al 1996) and the subgrid contributions are implicitly accounted for via turbulence closure. Furthermore, most of the higher-order moment statistics are not readily obtainable from low-order turbulence closure schemes, for example, 1.5-order closure or the first-order closure scheme, as implemented in most LES models and cloud-resolving models.

Probability density function (PDF)
Preceding analyses suggest that the LES results agree reasonably well with observations only for the lower-order moments of W such as horizontal averages and σ W in high model resolution runs. The agreement between them increasingly worsens with coarser model resolution and for  higher-order moments of W . The poor performance of the model in generating higher-order moments implies the possibility of its failure to simulate strong local fluctuations. To confirm this point, we further examine, in this section, the probability density functions (PDF) of W at different altitudes.  From the cloud base to the cloud top, all the PDFs tend to broaden, suggesting that turbulent vertical motions strengthen with increasing altitudes. The PDFs from the base run of 'X25Z10' compare best with the observations. With decreasing vertical resolution, the PDFs become narrower and the averaged magnitudes of W decline. Then accordingly, the difference between the observations and the model outputs becomes larger, especially over the 'tails' of the PDFs. Near the cloud base, the underestimation of W is more significant than those at the middle of cloud and near the cloud top. Especially in 'X25Z40', the magnitude of W is less than 0.02 m s −1 (figures 1 and 3) while the observed W could reach 0.16 m s −1 . This underestimation is expected to underpredict activated aerosol particles, and would potentially bias the estimate of the first aerosol indirect effect. Nevertheless, the underestimation could be gradually alleviated by refining model resolution to some degree.
Although the model results capture the gross features of the observed PDFs reasonably well, they fail to reproduce the tails of the observed PDFs (figures 3 and 4). Compared to the model outputs, the observations suggest much fatter tails, and these differences increase as the magnitude of W increases (toward the 'tail' of W ). It also is evident that the magnitudes of relative differences between the observations and the model outputs increase up to ∼100% when the magnitude of W reaches or exceeds 0.1 m s −1 (figure 4), suggesting that the model then increasingly under-represents the local but relatively strong vertical motions. This underrepresentation reinforces our findings from the analyses of the standard deviation, skewness, and kurtosis discussed in section 3.1 (figures 1 and 2).

Power spectrum
In addition to the PDFs, another aspect of the variability of a variable is auto-correlation and/or variability at different scales (frequencies). One way to measure this aspect of the variability is via the power spectrum in Fourier space (Mason and Brown 1999). In other words, variables with the same PDF can exhibit very different power spectra. Generally, multi-scale variability can arise from different PDFs and/or auto-correlations. For example, Mandelbrot (1997) refers to the variability with a fat-tailed PDF as the 'Noah effect', and the variability with the long-range dependence (auto-correlation) as the 'Joseph effect', respectively. Some studies have indicated that turbulent clouds may feature both types of variability (Marshak et al 1997, Davis et al 1996. Power spectral analysis provides information on the auto-correlation of the variable (of interest) and locates the dominant frequencies, and hence, complements the PDF analysis. To illustrate this point, figure 5 shows the power spectrum of W at H = 0.28 km where the agreement in PDFs between the airborne observations and the model outputs is the best (figure 3(b)). For the airborne observations, the power spectrum was calculated by performing the periodograms of W at a given sampling altitude. For the three-dimensional model outputs at a given altitude, these one-dimensional power spectra were calculated by performing the periodograms of W in the x-direction and then averaging over the y-direction and over time from 9:30 to 11:30 LST. The power spectrum of the observed W clearly follows the celebrated Kolmogorov's −5/3 power law over a range of wavenumbers from 10 −4 to 3×10 −2 m −1 , so providing some confidence in the observations. However, although the power spectra from the tests of 'X25Z10', 'X25Z20', and 'X25Z40' are encouragingly similar, the power density of the modelled W generally is lower than the observed value, suggesting that the model underestimates the variations in W at various scales (especially at small scales). Nevertheless, with increased model resolutions, the power spectrum tends to show a slope closer to the '−5/3' spectrum before falling off sharply; also the power density is higher (especially at high wavenumbers) and decays more slowly. Hence, higherresolution simulations are desirable to better resolve the smallscale variability and long-range dependence of W , and to reach closer agreement with observations as well as with theoretical predictions. (Note: we conducted spectral analyses at H = 0.07 and 0.37 km as well, and obtained similar spectral behaviours for W from both the observations and the model outputs.) Furthermore, the power density of the modelled W exhibits a sharp spectral fall-off from a certain wavenumber (figure 5), although this sharp fall-off occurs at higher wavenumber in higher-resolution runs. This falling-off is associated with a filter operation that separates resolved scales from subgrid scales. The filter length (λ) often is proportional to the grid's spacing, as is adopted in the turbulence (subgrid) parameterization through which the resolved motions are aware of the subgrid motions (Tao and Simpson 1993). With increased resolution (or smaller grid spacing), λ becomes smaller since λ is positively correlated with the grid's size ( x, y, and/or z), so that less energy then is represented by the subgrid processes that impose a dissipative effect upon the resolved-scale processes. This means that the dissipative effect of the subgrid processes becomes effective at smaller-scale processes with increasing the model resolution. Accordingly, a −5/3 slope behaviour could be better represented in smaller scale (or higher wavenumber) processes, rather than smeared out by the subgrid dissipation. In other words, this rapid dropoff begins at higher wavenumbers in higher-resolution runs.

Structure function
Another way to examine auto-correlation and range-dependence is to examine the so-called qth structure function (Monin and Yaglom 1975), defined as the qth-order statistical moments of the absolute changes in W ( W (r ) q ) across the horizontal distance r : 1, 2, 3, 4, . . .) where · refers to ensemble averages, and here we substituted it with horizontal average; x denotes spatial position. Previous studies have suggested that the statistical moments, such as W (r ) q are likely to be independent of x (see Davis et al 1996, wherein there is a good description of the differences and similarities between power spectral and structure function analysis). We calculated different-order structure functions at H = 0.28 km to uncover any statistical differences of the horizontal variations in W . Figure 6 presents the 2ndorder structure functions of W ( W (r ) 2 ) versus the distance normalized by the grid size (r/ x). For the observed W , the magnitude of W (r ) 2 increases with the distance, r , for 0.025 km < r < 25 km, and then declines for r > 25 km (figure 6), indicating that the observed W is correlated at (spatial) scales from 0.025 to 25 km, beyond which it becomes less well correlated.
The observed W exhibits significant variations on small scales (such as x), but the modelled W likely fails to reveal such variations. In 'X25Z10', the model can capture the variations in W on a scale equal to or larger than 2 5 x = 32 x = 0.8 km (figure 6). In the coarser resolution runs such as 'X25Z20' and 'X25Z40', the model only captures the variations on the scale of 2 6 x = 64 x = 1.6 km or even larger. Moreover, the increase of W (r ) 2 with r/ x is faster in the coarser resolution runs, indicating that the simulated W becomes more and more statistically non-stationary, therefore, to obtain meaningful statistics for W , a longer dataset (and/or over a larger domain) is needed (Davis et al 1996).
The magnitude of W (r ) 2 of the observed W generally is larger than that of the modelled W , that is, the former usually has larger spatial variations, especially at smaller r/ x. The magnitude of W (r ) 2 increases with r/ x, but this increase is much slower for the observed W than for the modelled one (comparing the red squares with black markers in figure 6). Consequently, the smaller the r/ x, the more the W (r ) 2 differs between them. The relative differences of the model outputs with respect to the observations always are positive and larger at smaller r/ x ( figure 7(a)).
The analysis of the 2nd-order structure in physical space basically confirms our findings about the second-order moment from the power spectral analysis in Fourier space. To further examine the long-range behaviour of the higher-order moments, figure 7(b) shows the relative differences of the first-, second-, third-, fourth-order structure functions (i.e., q = 1, 2, 3, and 4) for the base run of 'X25Z10'. As evident from this figure, the relative difference generally decreases with r/ x no matter what the value of q is, suggesting that the model outputs become more comparable (in magnitude) with observations over longer averaging scales (or over larger domains). It is especially noteworthy that the relative differences are bigger for larger q (or higher-order structure functions), implying that even for a run of the highest resolution, the model performs increasingly poorly for higher-order structure functions. These results are understandable since stronger fluctuations at smaller scales carry relatively larger weights and become more distinct in the higher-order moments and/or structure functions. This relationship also underlines the deficiency of the models' capability to capture these small-scale variations and/or the higher-order moments and structure functions of W (such as fourth-order ones). To obtain more robust numerical results, we need more powerful computer resources to perform simulations at higher resolution.
We also calculated the 2nd-order structure functions at H = 0.07 and 0.37 km, and different-order structure functions for the tests of 'X25Z20' and 'X25Z40', and noteworthily, the results were qualitatively similar to the above. Hence, they are not detailed in our discussion.

Concluding remarks
We examined various characteristics of the vertical velocity (W ) within the marine stratus-topped boundary layer during the MASE (Marine Stratus/Stratocumulus Experiment) both from aircraft's in situ observations and from large eddy simulations. Our main purpose was to characterize W and its vertical and horizontal variations, and to examine the fidelity of large eddy simulations at reproducing them.
The horizontal average W at different latitudes are close to 0 m s −1 from both the airborne observations and the model outputs, suggesting that local fluctuations due to boundary layer turbulence are major contributors to local updrafts/downdrafts.
The magnitude of the vertical profiles of the standard deviations is comparable in magnitude between the observations and the model outputs. Generally, the skewness is negative within this cloud layer, indicating that cloud top radiative cooling mainly drives this stratus cloud. The model outputs underestimate the kurtosis compared to the observations. The agreement between observations and model outputs generally deteriorates, with coarser model resolution, for higher-order moments of W , and towards the cloud top and/or cloud base.
Probability density function (PDF) analysis shows that although the model captures the gross features of the observed PDFs reasonably well, it fails to reproduce the 'fatter' tails of the observed PDFs. Power spectral analyses also point to an underestimation of small-scale variations in the model outputs. Further structure function analyses reveal that smallscale variations (here over a length scale of x = 25 m) still are significant, but the model probably under-resolves them and thus smears them out. It also suggests that the scaleindependent statistics might be obtained only over a spatial scale of at least 25 km. Therefore, to achieve meaningful statistics, we need data not only from high-resolution runs and/or samplings down to metres (or even more refined), but also over larger spatial scales up to tens of kilometres or greater.
The combined results indicate that the model does not catch the relatively strong vertical motions at small scales, and its agreement with observations worsens for higher-order moments and/or structure functions of W and with coarser model resolutions. A recent study by Guo et al (2008) found that a minimum vertical resolution of 40 m is required to capture cloud diurnal variations. The results in this study suggest an even higher vertical resolution (i.e., finer than 10 m) is necessary if the variability property of vertical motions is of primary concern. Considering that small-scale variability in real turbulent clouds is ubiquitous, and that both the so-called 'Noah effect' and 'Joseph effect' are likely the norm in any cloud, it is evident that all these aspects of variability in the vertical motion are critical for activating aerosols into cloud droplets and the subsequent droplet-growth processes (Shaw 2003, Kostinski andShaw 2005).
In addition to the local fluctuations, the mean vertical velocity is also of profound significance in the assessment of aerosol indirect effects. For example, (in deep convection) stronger updraft cores in concert with the suppression of precipitation by aerosols can transport more cloud water to a higher altitude, where the cloud water releases more latent heat by freezing, reinforces convection strength, and thus magnifies the aerosol indirect effect (Khain et al 2005).
Current deficiencies of LES models in representing the multi-scale variability in vertical velocity necessitate much finer numerical simulations to resolve small eddies, extending over a larger spatial domain to encompass mesoscale organizations of stratiform and cumulus convection, and runs of longer duration to cover lifecycles of cloud systems.
Furthermore, the statistics of higher-order moments are not readily obtainable via the low-order turbulence closure schemes as commonly implemented in cloud models. These requirements pose new challenges for high-performance computer resources, effective and efficient numerical techniques, and high-resolution simulations under the limited computer resources.
Nevertheless, there exist potential measurement problems, e.g., flow distortion, attack angle, amplification or reduction of pressure fluctuations, and time delay of the circulation generated by the wing Wang 2002a, 2002b). Moreover, for higher-order moments, smaller time (length) scale fluctuations are important. But their smaller power, because of the decrease of the power density with wavenumber to the '−5/3' power (figure 5), requires more demanding resolution and timing, and thus makes the measurements of higher-order moments more challenging, in addition to the required large sampling size (Lenschow et al 1994).