A grid refinement study of trade wind cumuli simulated by a Lagrangian cloud microphysical model: the super‐droplet method

The impact of spatial resolution on the simulation of trade wind cumuli was investigated. The super‐droplet method, an efficient stochastic Lagrangian cloud microphysical model, was used to reduce uncertainties due to the empirical parameterisation of cloud microphysics and numerical diffusion for advection, which is inevitable in an Eulerian cloud microphysical model. We showed for the first time that cloud cover numerically converged with a grid resolution of 12.5 m. Our grid refinement analysis elucidated a significant contribution of small cumulus clouds to total cloud cover, as such cumuli are generated by small‐scale updrafts that can be resolved only at a fine resolution.


Introduction
Shallow clouds are one of the main sources of uncertainty in climate prediction (IPCC AR5; Stocker et al., 2013). General circulation models (GCMs), which have been used for climate prediction usually express clouds by parameterisations (e.g. Tiedtke, 1993;Considine et al., 1997;Kain, 2004), but they cannot simulate shallow clouds explicitly. Cloud-resolving models based on the large-eddy simulation (LES) technique are powerful tools for reproducing small-scale phenomena associated with clouds and have often been used to improve cloud parameterisation (e.g. Golaz et al., 2007;Bretherton and Park, 2009;Kogan, 2013).
Among these issues, the sensitivity to resolution is a key issue of LES. It is still not clear how small a grid must be to obtain an accurate numerical solution of shallow cumuli. For example, in terms of cloud microphysics, Grabowski and Jarecka (2015) determined that a vertical grid resolution of less than 10 m is required for correctly simulating supersaturation and condensation growth at the cloud base. As well as the effects of the grid resolution on cloud microphysics, Matheou and Chung (2014) investigated the effects on the turbulent structure within shallow cumuli. They indicated that 90% of the total turbulent kinetic energy (TKE) should be resolved to obtain a numerically converged mean profile of trade wind cumulus and 20-m horizontal grid spacing is required. Brown (1999) indicated that the ensemble-mean statistics in shallow cumulus are sensitive not to grid spacing, but to the size of individual clouds, which is directly related to cloud cover and has a strong dependency on grid resolution. Matheou et al. (2011) showed that the grid convergence of the mean liquid water profile was achieved with their finest resolution, namely 20-m grid spacing, for a non-precipitating case. However, they could not confirm the convergence for cloud cover. Because cloud cover is a crucial quantity for radiative transfer calculations in GCMs, it is an important quantity, as a reference solution, to construct parameterisation for GCMs.
The cloud microphysical scheme is also crucial for simulated clouds. Bulk schemes and spectral bin schemes have been widely used in LES models. These models are Eulerian cloud models (ECMs), in which prognostic variables associated with cloud are defined on grid points. This treatment has uncertainty caused by numerical diffusion, which is artificially added or implicitly included for model stability and is inevitable in ECMs. (Henceforth, we refer to this as numerical diffusion).
Recently, several Lagrangian cloud microphysical models (LCMs) have been developed (e.g. Andrejczuk et al., 2008;Shima et al., 2009;Sölch and Kärcher, 2010;Riechelmann et al., 2012). Because LCMs explicitly solve the microscopic time evolution equations of particles, they can reduce the uncertainties derived from empirical parameterisations and assumptions related to several cloud microphysical processes (e.g. autoconversion, accretion, saturation adjustment and the size-distribution function of cloud particles). Furthermore, the cloud microphysical variables in LCMs are defined not as Eulerian grid variables but as Lagrangian particles over the calculation space. Due to this particle-based representation, LCMs can also reduce the error derived from the numerical diffusion of cloud microphysics.
To understand the influence of the grid spacing of LES, the grid refinement approach using LCMs is useful because it can reduce uncertainties in cloud microphysics. Arabas and Shima (2013) investigated the influence of LES resolution on shallow trade wind cumuli using the super-droplet method (SDM) (Shima et al., 2009), which is a computationally efficient stochastic LCM. They used a 25-m grid spacing, but the grid convergence on the cloud microphysical properties could not be fully confirmed.
Considering the background information presented above, the aim of this study was to clarify the grid spacing required to produce a well-converged numerical solution of shallow cumuli using an LES model coupled with the SDM. We systematically investigated the influence of resolution on shallow non-precipitating trade wind cumuli, especially their cloud cover, through the grid refinement approach.

Model and experimental setup
The LES model used in this study was constructed from the Scalable Computing for Advanced Library and Environment (SCALE) library Sato et al., 2015). The SDM (Shima et al., 2009), an LCM using a Monte Carlo scheme for the stochastic coalescence process, was implemented into SCALE (an Eulerian dynamical core). In the SDM, aerosol/cloud/precipitation particles were represented by Lagrangian particles called super-droplets (SDs). In this study, gravitational settling, activation/deactivation, condensation/evaporation and coalescence processes were considered. Nucleation scavenging and rain out of aerosol particles were considered, i.e. the aerosol number concentration changed over time. The influence of sub-grid-scale (SGS) turbulence on the SDs was ignored.
The experimental setup was based on an intercomparison study called the Barbados Oceanographic and Meteorological Experiment (BOMEX) (Siebesma et al., 2003). The domain size of the model was 6.4 2.1 × 10 10 2.6 × 10 9 3.3 × 10 8 4.0 × 10 7 5.1 × 10 6 (SD number per one grid) a The grid aspect ratio of horizontal to vertical grid spacing was kept at 0.8 because a large aspect ratio results in artificial error in the turbulence structure . b The time step of the dynamics was determined by the Courant-Friedrichs-Lewy (CFL) condition for an acoustic wave. c Microphysics includes aerosol activation/deactivation, condensation, evaporation and coalescence processes.
× 7.2 × 3.5 km 3 , which was almost the same as that used by Siebesma et al., (2003). To investigate the effects of resolution on cloud microphysical properties and turbulence structure, we conducted sensitivity experiments for spatial resolution. Both horizontal and vertical grid resolution were varied without changing the grid aspect ratio. A summary of the resolution and time step of each experiment is shown in Table 1. To maintain model stability, numerical diffusion with a fourth-order Laplacian operator was added to the Eulerian prognostic variables (density, three-dimensional momentum, density-weighted potential temperature and vapour mixing ratio). The diffusion coefficient was given non-dimensionally as = 1 × 10 −3 , based on a sensitivity experiment of shallow cumuli . For all cases, the number concentration of SDs was fixed as 30 per grid box at the initial time, which was sufficient to simulate cumuli using the SDM (Arabas and Shima, 2013). The total number of SDs used is also summarised in Table 1. The aerosol size distribution at the initial time was set by another intercomparison study, which targeted shallow cumuli (van Zanten et al., 2011). Figure 1 shows the temporal evolution of cloud cover and the liquid water path (LWP) averaged over the whole calculation domain, together with the results of a previous intercomparison study (Siebesma et al., 2003). Figure 1 also shows the LWP averaged over whole cloudy grid. Cloud cover was defined as the proportion of cloudy columns in the total column that had more than one grid with a cloud water mixing ratio greater than 0.01 g kg −1 . Cloud cover with a coarse grid resolution (Δx = 50 and 100 m), which was close to that of the previous study, was mostly within the range of the previous study. However, the cloud cover increased as the resolution increased. This dependency of cloud cover on resolution was similar to that reported in a previous study targeting shallow cumulus (Matheou et al., 2011) using an ECM. The authors reported that, for the non-precipitating case, the convergence of cloud cover was not achieved with their finest resolution although LWP was converged. The results of the present study, with the finest resolution of 6.25 m, indicate that cloud cover was also converged for the 12.5 m resolution. In contrast to cloud cover, resolution had a small effect on LWP (Figure 1(b)), even though LWP averaged over whole cloudy grids has dependency on the resolution (Figure 1(c)). This indicates that the total mass of cloud water generated by the convection triggered by the surface flux was the same regardless of resolution, whereas the three-dimensional spatial distribution of liquid water was sensitive to resolution. Figures 2(a)-(e) shows the spatial distribution of the LWP simulated at each resolution. The cloud patterns gradually changed as the resolution increased. As Brown (1999) indicated, the size of an individual cloud is sensitive to resolution. Figure 3(a) shows the frequency histogram N(S) of the horizontal area S of a cloud using uniform bins on the log (S) axis. N(S) can be fitted by a power law S -1/2 , which agrees fairly well with the previously reported fractal properties of fair weather cumuli analysed from Landsat satellite data (Cahalan and Joseph, 1989). From this power law, the contribution of clouds in the size range log (S) to log (S) + Δlog (S) to cloud cover can be evaluated as SN(S) ∼ S 1/2 . Because S 1/2 is a slowly increasing function of S, our analysis suggests that the largest clouds contribute the most to total cloud cover, but smaller clouds also make a notable contribution. Cloud cover values originated from clouds with a small (large) area averaged over the last 1 h, whose S was smaller (larger) than 0.1 km 2 , were 0.061, 0.063, 0.047, 0.029 and 0.016 (0.220, 0.220, 0.150, 0.109 and 0.0726) for Δx = 6.25, 12.5, 25, 50 and 100 m, respectively. These results support our argument that there is a notable contribution (about 20%) of small clouds to cloud cover. As well as small clouds, the cloud cover from both small and large clouds numerically converged, with a grid resolution of Δx = 12.5 m.

Results and discussion
The dependency of the spatial distribution of clouds on resolution originated from differences in the structure of turbulence, which was visualised by the vertical velocity field (w) near the surface (z = 50 m; Figures 2(f)-(i)). The w field had a striped pattern, with a narrow and strong upward velocity in fine-resolution runs (Δx ≤ 12.5 m), which is a characteristic of shear-driven turbulence (Moeng and Sullivan, 1994). On the other hand, a broader hexagonal shape with weaker upward velocity, which is typical in buoyancy-driven turbulence (Moeng and Sullivan, 1994), appeared in coarse-resolution runs (Δx ≥ 50 m). The dependency of the spatial pattern of LWP and w is seen during last 4 h of the simulation (see Videos S1-S10, Supporting information).
These differences in the turbulence structure can be explained through an analysis of the vertical profile of the liquid water mixing ratio (q l ) (Figure 3(b)), vertical flux of total water (Figure 3(c)), and skewness of the vertical velocity of grid-resolved turbulence (w ′ ; Figure 3(d)). With a coarse resolution, a positive value of q l was reached at the surface, corresponding to surface precipitation. The surface precipitation amounts averaged over the last 4 h were 1.12, 0.35, 0.03, 0.002 and 0.016 mm day −1 with Δx = 100, 50, 25, 12.5 and 6.25 m, respectively. Although the precipitation with Δx = 12.5 m was smaller than that with Δx = 6.25 m, the difference (0.014 mm day −1 ) was small enough that the precipitation amounts could be regarded as  equal. The surface precipitation with Δx = 12.5 and 6.25 m was regarded as the non-precipitating. Thus, the surface precipitation gradually decreased as the resolution became finer, and the non-precipitating cumuli of the BOMEX case could be reproduced with Δx = 12.5 and 6.25 m. The hexagonal structure of w shown in Figures 2(i) and (j) could be produced due to the lack of grid resolution (Piotrowski et al., 2009). The strong upward flux of total water below the cloud layer, shown in Figure 3(c), produced strong precipitation. The strong surface precipitation triggered the mass divergence in the precipitating area near the surface. Consequently, the circulation enhanced the strong convection and precipitation via positive feedback, as suggested by Feingold et al. (2010). This circulation enhanced the strong convection and precipitation via positive feedback. The circulation triggered by the surface precipitation resulted in a strong upward flux of total water below the cloud layer, as shown in Figure 3(c). The strong surface precipitation with coarse grid resolution would reduce after implementing the influence of SGS turbulence on the SDs. On the other hand, the striped pattern of w was constructed by a narrow and strong upward velocity that appeared in the high-resolution runs, which generated small clouds. The skewness of w ′ clearly shows this distinctive feature of the w field. Although the skewness was included among the results of the previous study, it was found to increase as the resolution became finer. Its numerical convergence was also achieved with Δx = 12.5 m (Figure 3(d)). These convergences of the vertical velocity field contributed to the convergence of cloud cover. From these results, we can conclude that the fine-resolution simulation (Δx ≤ 12.5 m) accurately captured the detailed structure of the vertical velocity field induced by turbulence. This led to changes in cloud structure and precipitation properties. Figure 4 shows the turbulence kinetic energy spectra E(k) at three altitudes, z = 250, 700 and 1200 m. These spectra were calculated from the velocity profile (u(x,y), the Fourier image of the velocity profile. E(k) was obtained by integrating E(k x ,k y ) along a circle of radius k = √ k 2 x + k 2 y and then averaged over the last half-hour of the simulation (i.e. t = 5.5 to 6 h) using snapshots taken every minute. The slope of E(k) in the cloud layer is flatter than Kolmogorov's − 5/3 power law, which is characterised by a constant energy cascade in wave number space (Kolmogorov, 1941a(Kolmogorov, , 1941b. Based on Kolmogorov's theory, this flatter energy spectrum suggests that the energy is not cascading constantly in the wave number space, but there exist input and output of energy at various scales. We can speculate that the cumuli population is maintaining the turbulence because they also have a self-similar size distribution as shown in Figure 3(a). More detailed analyses to understand the mechanism will be conducted in the future.

Concluding remarks
In this study, the influence of spatial resolution was investigated through a grid refinement approach using three-dimensional LES simulations. For the cloud microphysics, we used the SDM, a stochastic LCM, to reduce the uncertainty of cloud microphysics and numerical diffusion as much as possible. The results of the sensitivity experiment indicated that a numerical convergence for cloud cover of trade wind cumuli was achieved with a horizontal and vertical grid spacing of 12.5 and 10 m. The detailed structure of turbulence, which can only be resolved with fine spatial resolution, played a notable role in simulating cloud cover.
As suggested by previous studies (Matheou et al., 2011), the grid size required for accurate simulation can differ for precipitating cumulus and stratocumulus. Further investigation targeting the properties of turbulence when using LCMs is required to improve our   understanding of the impact of spatial resolution on the detailed structure of turbulence. These analyses also enable a more comprehensive understanding of the origin of the flatter turbulence energy spectrum (Figure 4). Comparing the grid size sensitivity between LCMs and ECMs is another important issue for future studies. Investigating all these issues should provide useful information to improve the cloud models used in mesoscale and global-scale models. KAKENHI

Supporting information
The following supporting information is available: Video S1. Movie of the spatial distribution of the liquid water path (LWP) from t = 4 h to t = 6 h simulated with Δx = 6.25 m.
Video S2. Movie of the spatial distribution of the liquid water path (LWP) from t = 4 h to t = 6 h simulated with Δx = 12.5 m.
Video S3. Movie of the spatial distribution of the liquid water path (LWP) from t = 4 h to t = 6 h simulated with Δx = 25 m.
Video S4. Movie of the spatial distribution of the liquid water path (LWP) from t = 4 h to t = 6 h simulated with Δx = 50 m.
Video S5. Movie of the spatial distribution of the liquid water path (LWP) from t = 4 h to t = 6 h simulated with Δx = 100 m.
Video S6. Movie of the spatial distribution of vertical velocity at z = 50 m from t = 4 h to t = 6 h simulated with Δx = 6.25 m.
Video S7. Movie of the spatial distribution of vertical velocity at z = 50 m from t = 4 h to t = 6 h simulated with Δx = 12.5 m.
Video S8. Movie of the spatial distribution of vertical velocity at z = 50 m from t = 4 h to t = 6 h simulated with Δx = 25 m.
Video S9. Movie of the spatial distribution of vertical velocity at z = 50 m from t = 4 h to t = 6 h simulated with Δx = 50 m.
Video S10. Movie of the spatial distribution of vertical velocity at z = 50 m from t = 4 h to t = 6 h simulated with Δx = 100 m.