Computational Domain Size Effects on Large-Eddy Simulations of Precipitating Shallow Cumulus Convection

Abstract

Idealized large-eddy simulations of shallow convection often utilize horizontally periodic computational domains. The development of precipitation in shallow cumulus convection changes the spatial structure of convection and creates large-scale organization. However, the limited periodic domain constrains the horizontal variability of the atmospheric boundary layer. Small computational domains cannot capture the mesoscale boundary layer organization and artificially constrain the horizontal convection structure. The effects of the horizontal domain size on large-eddy simulations of shallow precipitating cumulus convection are investigated using four computational domains, ranging from $40\times 40{\mathrm{km}}^2$ to $320\times 320{\mathrm{km}}^2$ and fine grid resolution (40 m). The horizontal variability of the boundary layer is captured in computational domains of $160\times 160{\mathrm{km}}^2$. Small LES domains (≤40 km) cannot reproduce the mesoscale flow features, which are about $100\mathrm{km}$ long, but the boundary layer mean profiles are similar to those of the larger domains. Turbulent fluxes, temperature and moisture variances, and horizontal length scales are converged with respect to domain size for domains equal to or larger than $160\times 160{\mathrm{km}}^2$. Vertical velocity flow statistics, such as variance and spectra, are essentially identical in all domains and show minor dependence on domain size. Characteristic horizontal length scales (i.e., those relating to the mesoscale organization) of horizontal wind components, temperature and moisture reach an equilibrium after about hour 30.

1. Introduction

Clouds forming in the atmospheric boundary layer play a crucial role in the Earth’s energy balance as they influence radiative transfer, surface energy fluxes and the hydrological cycle [1,2,3]. Shallow clouds, in particular, have been identified as one of the main sources of uncertainty in climate projections [1]. Precipitation alters the boundary layer spatial structure, leading to large-scale organization [4]. Accurate simulations of shallow cumulus convection are thus essential for improving climate models and reducing uncertainties in future projections.

Large-eddy simulation (LES) is a major research tool in atmospheric boundary layer studies. LES is a “turbulence resolving” model that yields a physically consistent description of the spatial variability over a large range of scales [5,6]. For example, in investigations of processes in shallow cumulus convection [7,8,9,10,11,12]. The domain size in LES can significantly impact the simulated cloud fields and boundary layer structure [13,14]. Particularly in idealized process-level LES investigations where periodic boundary conditions in the horizontal directions are used. A growing body of literature has explored the influence of domain size on LES of shallow cumulus convection [6,13,14,15,16,17,18]. For instance, refs. [12,14] demonstrated that increasing domain size could lead to the formation of larger convective structures.

Several studies focused on the effects of precipitation on shallow cumulus convection [19,20]. Precipitation can alter the spatial organization of clouds, leading to the formation of large persistent structures [1]. Even in non-precipitating cumulus convection, variance accumulation at the largest scales and generation of large horizontal structures is observed [21]. The interplay between domain size and precipitation remains relatively unexplored. Lamaakel and Matheou [12] showed that larger domains were required to accurately capture the cloud organization in precipitating shallow convection. Whereas [12] studied the development of shallow convective organization in relatively large LES computational domains, up to $160\times 160{\mathrm{km}}^2$, presently we focus on the boundary layer steady state, which occurs after mesoscale organization has been formed.

This study examines the effects of computational domain size on the equilibrium state of a typical case of trade-wind precipitating shallow cumulus convection. The present investigation is the extension of previous work [12] to steady state boundary layers. One important open question from that study [12] is whether the horizontal length scales will continue to grow or reach a steady state. Previous LES investigations suggested that feedbacks can influence mesoscale variability of domain-mean properties and “very large computations are required to obtain meaningful cloud statistics” [15]. To help achieve this goal, an additional run to those of [12] with an unprecedented combination of domain size ($320\times 320{\mathrm{km}}^2$ in the horizontal directions) and fine grid resolution ($\Delta x=40\mathrm{m}$) is carried out.

Section 2 briefly describes the LES model and the setup of the simulations. Results are discussed in Section 3. The LES output for different domain sizes is compared with respect to domain-averaged statistics versus time, vertical profiles at the end of the run, spectra and horizontal length scales. First, the evolution of cloud organization is presented, followed by physical-space turbulent flow statistics. In the second part of the results, spectra and length scales are discussed. Finally, conclusions are summarized in Section 4.

2. Methods

The LES model of [22] is used to simulate shallow convection in a trade-wind atmospheric boundary layer. The LES model numerically integrates the anelastic approximation of the equations of motion on an f-plane. Warm rain processes are parameterized with the two-moment bulk warm-rain parameterization of Seifert and Beheng [23] and turbulent transport is parameterized with the buoyancy-adjusted stretched-vortex subgrid scale (SGS) model [24]. The buoyancy-adjusted stretched-vortex model is a state of the art structural turbulence closure that accounts for both the anisotropy of the SGS flow and density-stratification effects. Momentum, liquid water potential temperature ${\theta }_l$ and total water mixing ratio $q_t$ advection terms are approximated with the fourth-order centered fully conservative scheme of [25]. Microphysical-variable advection terms are discretized with a monotone flux-limited scheme to ensure water-mass conservation. The third-order Runge–Kutta method of [26] is used for time integration. To limit gravity wave reflections, a Rayleigh damping layer is used at the top $0.5\mathrm{km}$ of the domain.

The idealized case corresponding to the Rain In shallow Cumulus over the Ocean (RICO) observational campaign is simulated [27,28]. The initial soundings have a typical trade-wind cumulus-topped boundary layer structure. Large-scale forcings include a constant in time geostrophic wind, subsidence, advective tendencies and clear air radiative cooling. Surface fluxes are dynamically computed using a constant sea surface temperature (SST) bulk parameterization following [28].

Simulations with progressively larger domain sizes are carried out. All domains are square in the horizontal directions, i.e., equal domain size lengths $L_x=L_y$. The boundary conditions are doubly periodic in the horizontal. The x and y coordinates are along the zonal and meridional directions, respectively. The smallest domain size has $L_x=40.96\mathrm{km}$, which is twice as large in each direction compared to the LES domains used in the model intercomparison [28]. The domain horizontal length doubles in each direction as the domain sizes increase in the present runs: $L_x=L_y=40.96$, 81.92, 163.84 and $327.68\mathrm{km}$. The setup of the runs is shown in

Table 1. Details of the large-eddy simulations. The number of grid points in the meridional x and zonal y directions is $N_x=N_y$ and the number of vertical grid points is $N_z$. The grid spacing is isotropic $\Delta x=\Delta y=\Delta z$. All domains are square in the horizontal directions with extents $L_x=L_y$. The height of the computational domain is $L_z$.

Run	${\mathit{N}}_{\mathit{x}}$	${\mathit{N}}_{\mathit{z}}$	${\mathit{L}}_{\mathit{x}}\left(\mathbf{km}\right)$	${\mathit{L}}_{\mathit{z}}\left(\mathbf{km}\right)$	$\Delta \mathit{x}\left(\mathbf{m}\right)$
A	1024	125	40.96	5	40
B	2048	125	80.92	5	40
C	4096	125	163.84	5	40
D	8192	125	327.68	5	40

. The simulations are labeled A to D in the order of increasing computational domain size. A uniform and isotropic grid spacing $\Delta x=\Delta y=\Delta z=40\mathrm{m}$ is used. The simulations are performed in a horizontally translating frame of reference with approximately the domain-mean wind $[-6,-4]\mathrm{m}{\mathrm{s}}^{-1}$ [11]. The simulations are run for $36\mathrm{h}$ to reach a statistically stationary state. A statistically steady state is expected after $t>30\mathrm{h}$ [8,10]. The case setup does not include any diurnal cycle effects and the forcings are constant through the duration of the run.

The two largest runs, C and D, are fairly large computations. The largest simulation consists of 8.3 billion grid cells. Large-scale organization develops after $t>16\mathrm{h}$, see [7,8,10,12,28], whereas during the first part of the simulations ($t<16\mathrm{h}$) homogeneous scattered cumulus develop and progressively deepen. As discussed in [12], smaller domains sufficiently capture the organization and flow statistics for $t<16\mathrm{h}$. Therefore, run D is initialized based on run C output at $t=16\mathrm{h}$, to save computing time. The output of run C is repeated in a $2\times 2$ pattern to span the twice larger in each direction domain D. Random temperature and moisture perturbations are added in the lowest $200\mathrm{m}$ to help randomize the periodic pattern. After $t>30\mathrm{h}$ run D becomes de-correlated from run C.

Model validation and sensitivity to model parameters are discussed in [22]. In addition, the present simulations yield cloud patterns, including their horizontal dimensions, similar to the observations in [27,29,30] and other model results [19,28,31]. At the present grid resolution $\Delta x=40\mathrm{m}$, the LES generates grid-independent results [22]. Further validation and assessment of the LES model for a diverse set of meteorological conditions is discussed in [7,32,33,34,35,36,37].

3. Results

3.1. Convection Organization

The organization of convection is visualized in Figure 1 and Figure 2 by plotting contours of cloud-liquid water path (LWP) and rain water path (RWP) fields at $t=36\mathrm{h}$. All four LES domains are shown in each figure directly displaying the differences in area coverage of each domain. The characteristic structure of shallow precipitating convection with cold pools, cloud arcs and occasional stratiform cloud anvils, similar to observations and previous modeling [8,10] is clear in the largest domains. A comparison between the LWP and RWP fields shows the locations of stratiform cloud anvils, which mostly register in the LWP field but not in RWP, because precipitation primarily develops in convective clouds. Most rain develops in the long cumulus clusters that form the cloud arcs at the leading boundary of the cold pools.

Structures longer than $100\mathrm{km}$ form in the boundary layer. The spacing between the cloud arcs is similar in the two largest domains C and D. Smaller domains, particularly domain A, are very small compared to the cloud arcs and cold pool structures of runs C and D. Smaller domains cannot realistically represent the size and spacing of cloud arcs and cold pool structures.

3.2. Domain-Averaged Turbulent-Flow Statistics

Figure 3 shows boundary layer statistics averaged over the computational domain volume. Panels correspond to the time evolution of the vertical integral of turbulent kinetic energy, $\mathrm{VTKE}\left(t\right)=\frac{1}{2}{\int }_0^{L_z}\rho \langle u^{\prime 2}+v^{\prime 2}+w^{\prime 2}\rangle \mathrm{d}z$, cloud-liquid (LWP, suspended condensate) and rain water (RWP) paths, cloud statistics (cloud top $z_c$, base $z_b$ and cloud fraction $cc$), inversion height $z_i$ and surface precipitation rate. The angle brackets $\langle •\rangle$ denote the instantaneous horizontal average. The inversion height $z_i$ is defined as the height of the maximum potential temperature, $\theta$, gradient. Cloud base and cloud top are defined as the minimum and maximum heights with a horizontal-mean cloud-liquid water mixing ratio $\langle q_l\rangle >{10}^{-5}\mathrm{kg}{\mathrm{kg}}^{-1}$. Lines for run D ($L_x=320\mathrm{km}$) start at $t=17\mathrm{h}$, an hour after initialization, to remove the initial impulse caused by adding random perturbations when domain D is initialized. Traces between runs C and D show visible differences by hour 25 and larger differences after hour 30. Data for simulations A–C for $t<30\mathrm{h}$ are the same as those of [12].

Runs C and D show good agreement with respect to domain-averaged statistics. The differences in cloud-top height between C and D are expected, because $z_c$ tracks the domain-maximum of the cloud top. As the domain area increases, the instances of locally higher cloud tops are more likely. As discussed in [12], the area fraction of cloud above the inversion is very small. Thus, high $z_c$ values are not representative of the entire area or the depth of the turbulent layer. See also TKE and $q_l$ profiles in the next section.

3.3. Vertical Profiles

Figure 4 shows profiles at $t=36\mathrm{h}$. Profiles are instantaneous horizontal averages without any time averaging. The total water vertical flux ${\overline{wq}}_t$ is additionally time-averaged in a ten-minute interval to create smoother curves. Turbulent fluxes include both the resolved-scale and subgrid-scale contributions. The time-variability of the profiles is small because all computational domains are relatively large and the horizontal averages include most of the horizontal variability of the boundary layer. The variation with respect to time of the vertical profiles is discussed in Appendix A.

Only the horizontal component of TKE exhibits large differences with respect to domain size. The differences in the TKE profile are caused by the horizontal component $\frac{1}{2}(u^{\prime 2}+v^{\prime 2})$. The profiles of Figure 4e,f are essentially identical. Note that the x-axis scale of Figure 4g, the vertical component of TKE, is different than Figure 4e,f. Further, most of the differences in the horizontal component between runs A and B are because of the $v^{\prime 2}$ component (not shown in Figure 4).

With the exception of the moisture flux and $q_l$, the smallest domain, run A, has the largest differences with respect to the other runs in Figure 4. The two largest domains, runs C and D, are in good agreement for u, $\theta$, $q_t$ and TKE. Most of the differences between all runs are observed near the inversion layer, $z\approx 2\mathrm{km}$. As discussed in the previous section, even though clouds occasionally rise up to $4\mathrm{km}$, TKE, $q_l$ and the moisture flux are small above the inversion, indicating that very deep clouds have very small area fractions. This is in agreement with previous run C results [12] before the boundary layer reached steady state.

The profiles of wind, $\theta$ and $q_t$ are essentially identical for all domains. Relatively small differences are observed near the inversion layer ($z\approx 2\mathrm{km}$) between the smallest domain (run A, $40\times 40{\mathrm{km}}^2$) and the rest of the simulations. Larger differences are observed in second-order flow statistics, covariances and TKE. In addition, TKE (Figure 4e) and temperature variance (Figure 5b) are in good agreement only in the two largest-domain runs, C and D. Interestingly, the cloud liquid at steady state (Figure 4d) does not exhibit large differences but temperature and moisture variance (Figure 5) profiles show larger differences, particularly between runs A–B and runs C–D.

3.4. Inversion Strength

The small differences in the mean $\theta$ and $q_t$ profiles near the inversion layer (Figure 4) and larger differences with respect to their variance (Figure 5) motivate the study of the inversion characteristics. Figure 6a shows the inversion strength, defined as the maximum vertical gradient of the mean potential temperature $\mathrm{d}\langle \theta \rangle /\mathrm{d}z$ as a function of time. After large-scale organization develops, $t>16\mathrm{h}$, the inversion strength evolution significantly differs between run A and the rest of the simulations. The smallest domain run sustains a stronger inversion strength. The inversion strength is a characteristic of the simulated temperature profile, which can be useful in understanding the dynamics of cloud development. For example, a strong inversion layer can limit the upward growth of clouds and restrict their vertical extent, while a weaker inversion layer may allow for more vigorous or deeper cloud development.

Figure 6b shows another diagnostic to help understand the dynamics of the inversion. At $t=36\mathrm{h}$ the largest computational domain (run D) is partitioned into $8\times 8$ subdomains in the horizontal. Each of the subdomains has the same horizontal area as domain A. In each of the run D subdomains we compute the inversion strength and the inversion height. These pairs are plotted in Figure 6b as open circles and quantify the inversion strength at the “local” level, whereas the lines of Figure 6a correspond to the entire-domain average of $\theta$. The domain average is plotted in Figure 6b as a filled circle, which is the same value as the run D line at $t=36\mathrm{h}$ in Figure 6a. Figure 6b shows that there is large within-domain variability of the inversion strength and height. Values of inversion strength in the D-run subdomains are similar to the inversion strength value of run A. The distribution of values in the y-axis of Figure 6b quantifies the undulations of the inversion in run D, which span about half a km within the area of the simulation. Because of the feedback between the rising cumulus-topped updrafts and the evolution of the trade-wind inversion, a causal relationship between inversion undulations and boundary layer depth cannot be inferred from Figure 6b.

3.5. Spectra and Length Scales

Lamaakel and Matheou [12] studied the growth rate of the horizontal organization in the RICO case. The growth of characteristic horizontal length scales was quantified during the transition from scattered cumulus to the mesoscale-organized cloud clusters shown in Figure 1. Between $20<t<30\mathrm{h}$ the horizontal characteristic length scales of u, v, ${\theta }_l$ and $q_t$ grow at a rate of 3–$4\mathrm{km}{\mathrm{h}}^{-1}$, which is in agreement with the growth rate of the developing cold pools. As discussed in [12], the multiple cold pools developing in domain C interact resulting in a boundary layer composed of both scattered cumulus and cloud arcs (Figure 1 and Figure 2). Lamaakel and Matheou [12] studied the first generation of cold pools developing in a fairly homogeneous scattered cumulus field. Presently, we focus on the boundary layer steady state where cold pools are present in various stages of their life cycle and they are interacting. The main question is if the characteristic length scales reach an equilibrium value.

The horizontal characteristic length scales are defined based on the one-dimensional spectra along the longitudinal direction. Similar results are obtained when considering the transverse direction [12]. For u, w, ${\theta }_l$ and $q_t$ the x-direction (zonal) is used. For v the meridional direction is used. The one-dimensional premultiplied spectra are shown in Figure 7. The spectra are normalized by the corresponding variance, e.g., for u, the pre-multiplied spectrum is $k_x{\mathrm{\Phi }}_{uu}\left(k_x\right){\langle u^{\prime }{u}^{\prime }\rangle }^{-1}$, where ${\mathrm{\Phi }}_{uu}\left(k_x\right)$ is the one-dimensional spectral function, which is computed by performing a two-dimensional discrete Fourier transform on a horizontal plane and then averaging along the $k_y$ wavenumbers. All spectra are computed in the subcloud layer ($z=360\mathrm{m}$) where the flow is fully turbulent, i.e., three-dimensional turbulence spans the entire layer. In Figure 7 the x-axis is converted from wavenumber to length scale, $l_x=2\pi /k_x$, to assist the physical interpretation.

The local maxima of the premultiplied spectra are used to identify energy-containing structures. The length scales correspond to narrow wavenumber bands where spectral energy is higher than the neighboring wavenumbers. In the spectra of Figure 7 we are interested in the change of the location of the maxima as well as the shape of the overall curves. Following [12] the spectra are smoothed by a Gaussian filter.

As observed in [12], the vertical velocity spectra do not depend on the domain size and the premultiplied spectrum always has a single peak. The present data extend the observation of [12] regarding the vertical velocity premultiplied spectra to all stages of the boundary layer evolution, from organization development to steady state. Moreover, the characteristic w-velocity scale is small, about $1\mathrm{km}$, and it is not influenced by the size of the computational domains presently used.

In general, the spectra of runs C and D are in good agreement compared to those of A and B runs. The spectra of runs A and B (except w-velocity) contain more energy at the smaller scales ($l_x<10\mathrm{km}$). Considering that the larger domains result in a more realistic representation of the boundary layer, the smaller domain sizes cause an artificial accumulation of fluctuations at small scales. Further, in some cases (e.g., $k_x{\mathrm{\Phi }}_{uu}$ and $k_x{\mathrm{\Phi }}_{qq}$, $k_x{\mathrm{\Phi }}_{\theta \theta }$ at $t=30\mathrm{h}$) the maxima for runs C and D are at scales larger than the A and B domains.

Figure 8 shows the time evolution of characteristic length scales $l_u$, $l_v$, $l_w$, $l_q$ and $l_{\theta }$, where the subscript corresponds to the variable (i.e., $l_q$ corresponds to $q_t$ length scales and $l_{\theta }$ to ${\theta }_l$). Only length scales computed from the output of runs C and D are shown in Figure 8 because their spectra are not sensitive to the domain size. At a given time, i.e., x-coordinate in the plot, multiple symbols for are plotted when the premultiplied spectra have multiple peaks. In Figure 8d, the radii of the first two cold pools $r_{cp1}$ and $r_{cp2}$ are also shown using the data of [12]. The growth rates of $r_{cp1}$ and $r_{cp2}$ agree well with $l_q$. Note that in Figure 8 the y axis of the $l_u$, $l_v$, $l_q$ and $l_{\theta }$ plots is shown in logarithmic scale to better show the large variability during the evolution of the length scales. Figure 8e corresponds to the vertical velocity length scale. A single length scale $l_w$ is always observed, which grows slowly with time, thus a linear y-axis is used.

Some spectra have sharp peaks between wavenumbers 2 and 4, e.g., Figure 7a,c,h. These peaks do not vary in time and were not included in the data plotted in Figure 8. These peaks are likely artifacts of the periodic domain and do not correspond to physical length scales. Data for domain D are shown after hour 24 to allow sufficient time for run D to de-correlate from the initial condition that was based on run C.

As discussed in [12], precipitation-induced large-scale organization develops after $t>16\mathrm{h}$. Horizontal large scales grow relatively fast. All variables, including w, have a small scale which is observed to grow in time, but at a much smaller rate compared to the larger scale during $20<t<30\mathrm{h}$. The large-scale growth terminates at about $t=30\mathrm{h}$ when the boundary layer is expected to reach steady state [10]. In steady state, the premultiplied spectra have peaks in 10–$100\mathrm{km}$ with most of the peaks concentrating in the higher values of this range. Moreover, in Figure 8 the length scales diagnosed from runs C and D show the same characteristics, suggesting that the domain size of run C is sufficiently large.

4. Conclusions

The effects of computational domain size are investigated in large-eddy simulations (LES) of a typical case of a trade-wind shallow precipitating cumulus boundary layer [27,28]. The LES domain size varies from 40 km to 320 km. The simulations focus on the equilibrium state of the boundary layer, after 36 h of shallow convection evolution. For most flow statistics the $40\times 40{\mathrm{km}}^2$ domain is too small to realistically capture the flow. The $80\times 80{\mathrm{km}}^2$ domain adequately captures some flow statistics but not the cloud mesoscale structure. The flow statistics in the two largest domains with sizes of $160\times 160{\mathrm{km}}^2$ and $320\times 320{\mathrm{km}}^2$ are in good agreement, suggesting that a domain equal or larger than $160\times 160{\mathrm{km}}^2$ is required to represent the mesoscale variability of the boundary layer.

The smallest domain LES ($40\mathrm{km}$) cannot reproduce the mesoscale flow features that are about $100\mathrm{km}$ long, but boundary layer mean profiles are similar to the results of the larger domains. Vertical velocity flow statistics, such as variance and spectra, are essentially identical in all domains and show minor dependence on the domain size. A central question of the present study is whether the horizontal length scales will reach an equilibrium in the larger domain or continue to grow as the computational domain area increases. The present simulations show that the large length scales (i.e., those relating to the mesoscale organization) of horizontal wind components, temperature and moisture reach an equilibrium after about hour 30. Figure 8 shows the evolution of the length scales and it contains the main conclusion. Based on previous studies [14], length scales are expected to continuously grow. However this type of growth rate is expected to be slower than the rapid mesoscale organization observed between hours 20 and 30 in the present simulations. In fact, the smallest characteristic length scale for all variables continuously increases throughout the duration of the present simulations.