Expanding Grids for E cient Cloud Dynamics Simulations Across Scales

With large eddy simulations (LES) and/or cloud-resolving models (CRMs), it is now possible to simultaneously simulate shallowanddeep convection.However, using traditionalmethods, the computational expense is typically very large, due to the small grid spacingsneeded to resolve shallowclouds.Here, themain purpose is to present a method that is computationally less expensive by a factor of roughly 10 to 50. Unlike traditional grid stretching of only the vertical z grid spacing, the present method involves expansion of the grid spacing in all coordinate directions (x,y,z) and time t. A ne grid spacing of O(10)-O(100) m can be used near the surface to resolve boundary layer turbulence, and the grid spacing expands to beO(1000)mat higher altitudes, which reduces computational cost while still resolving deep convection. Example simulations are conducted with a simpli ed LES/CRM in 2D to verify the theoretical cost savings.


Introduction
Clouds and convection occur across a range of scales, and the grid spacing needed in numerical simulations can be di erent for di erent cloud phenomena. For instance, a grid spacing of roughly 1 km is needed to simulate deep convection [5,6], whereas a grid spacing of roughly 10 to 100 m is needed to simulate shallow clouds [25,27].
To simulate both shallow and deep convection together, the computational expense is very large. Not only is a small grid spacing needed to resolve shallow clouds, but a large domain is also needed to encompass the scales of deep convective clouds. It is only somewhat recently that computational power has increased to a point where shallow and deep convection could be simulated together in large eddy simulations (LES) [8,13,18].
Here we present a setup that could bring a signi cant computational savings. The basic idea is illustrated in Fig. 1, and it involves the use of a ne grid spacing of roughly 10 to 100 m only in the lowest portions of the atmosphere where such high resolution is needed to resolve shallow clouds. At higher altitudes, the grid expands to a coarser grid spacing of roughly 1 km, thereby bringing a large computational savings, while still allowing deep convective clouds to be resolved. In contrast to the common approach of stretching of a grid in only the vertical z direction, the present method involves an expansion of the grid spacing in all coordinate directions x, y, z and time t, which brings additional computational savings.
The setup of the expanding grid could be useful in any setting where shallow clouds and deep convection are simulated together. For instance, while superparameterization techniques for climate models were  The expanding grid shown is a schematic diagram and has 5 grid sections, and each section has a height of 4 grid cells (except for the lowest grid section, which has 8 grid cells). In the numerical experiments, the expanding grid uses 6 grid sections, and each section has a height of 10 grid cells (except for the lowest grid section, which has 20 grid cells); and the nest and coarsest grid sections have grid spacings of 50 and 1600 meters, respectively. originally targeted at resolving deep convection without resolving shallow clouds [4,21], new con gurations are currently being explored to resolve shallow clouds as well [3,11,19,24]. Also, while it is perhaps too computationally expensive at the present time, in the future it is feasible to consider the possibility of resolving shallow clouds in weather prediction and global CRMs [22,23,32]. At the present time, for LES of shallow and deep convection over limited areas [8,13,18], grid nesting is often used to combine horizontal regions of different horizontal grid spacing; the expanding grid could possibly allow for additional computational savings that could be reinvested elsewhere, such as in cloud microphysics and subgrid-scale processes [7,26,30], or in larger domain sizes to include additional mesoscale or synoptic-scale variability, or in larger ensembles of simulations.
This paper is organized as follows. In section 2, we describe the expanding grid and the setup of the numerical simulations. In sections 3, 4, and 5, we compare numerical simulations of the expanding grid versus uniform grids with ne or coarse grid spacings. A setup with 2D simulations is used to allow comparison with ne-grid-spacing simulations at reasonably achievable computational cost. The numerical simulations also provide a veri cation of the theoretical estimates of computational cost savings, which are described in detail in the appendix.

Numerical Setup
The expanding grid is set up in the following way. The domain is decomposed into M grid sections, which we denote as G i , for ≤ i ≤ M. See Fig. 1 for a schematic illustration of the setup (but note that the actual grid setup di ers from Fig. 1 in some details, due to space limitations for a single page of the journal; these di erences are described in the caption of Fig. 1). Each grid section has uniform grid spacing of h i in all spatial directions (so ∆x i = ∆y i = ∆z i = h i ), with time step of ∆t i . The grids are re ned in the sense that h i = h i− / , and ∆t i = ∆t i− / . For the speci c tests of the present paper, we use 6 grid sections, where grid section G extends from an altitude of kilometers to kilometers with a grid spacing of meters and a time step of seconds. The thickness of grid G i is then half the thickness of G i− for i > , so that G uses h = m and covers the lower most kilometer in the vertical direction. Other con gurations are also possible by, e.g., changing the number of grid sections M, allowing di erent numbers of vertical grid cells per grid section, etc. We choose the present setup because it can be described easily in terms of powers of 2 and because it meets the basic requirements on grid spacings for simulating both shallow and deep convection (see section 1). For additional tests, we have also used setups with and grid sections, so that the cost of inserting or removing a ne grid adjacent to the bottom of the domain can be determined.
To implement an expanding grid setup, two options are as follows. First, one could use existing adaptive mesh re nement (AMR) software in a static con guration that remains the same at each time step, which would essentially amount to a nested grid setup [2,29,31]. One advantage of this is its exibility, as AMR techniques can generally re ne any region of the domain and can even use dynamic/adaptive mesh re nement criteria. A potential disadvantage would be the di culty in incorporating AMR software into existing codes, so this approach may be best implemented in new codes that are designed from the start with AMR capabilities in mind, such as those used in [14][15][16]28], and [10], where the focus has traditionally been on horizontal rather than vertical re nement. In AMR, each grid cell is often represented as the node of a tree, the depth of which is equal to the number of levels of re nement at a given point in space. So to incorporate AMR into an existing code, one would need to change the array data structures into these tree-based ones. Besides this additional di culty, one would still expect an AMR code with the expanding grid de ned in a statically re ned con guration would o er the same speedup as the second approach.
As a second approach, which is taken in the present paper, modi cations could be made to an existing code. In making the modi cations, the basic idea is to leverage the fact that an existing code can solve the dynamics over an individual grid section, which is arranged as a regular array. What is needed in addition, then, is a wrapper code that calls the existing code to update each grid section, as needed. This calling procedure is then implemented as in AMR and grid nesting methods to move forward all grid sections one full time step [1,20]. For example, since grid section G uses ∆t = s and grid section G uses ∆t = s, grid section G must advance two ∆t time steps for every one ∆t time step.
The data structures for an expanding grid require some special treatment, but can be accomplished in a simple way using native Fortran capabilities. In particular, the numerically expanding grid is not arranged in space as a regularly spaced array. Nevertheless, each of the individual grid sections is a regularly spaced array. Therefore, the many regularly spaced arrays can be collected together and stored using a Fortran derived data type, and individual grid sections can be accessed e ciently using pointers. Note that this setup could possibly be extended so that grid sections could be inserted or removed based on dynamic criteria by using Fortran's dynamic memory allocation, although the simpler static con guration is used for the present initial study.
As the standard case, the domain size is kilometers in the vertical and 16 kilometers in the horizontal direction. To provide some additional tests, we also consider cases where the horizontal length of the domain is increased to be 24 km or 32 km. The simulations presented here are all 2D in order to allow for comparison with simulations with the ne uniform grid at a reasonable computational cost. Theoretical cost estimates will be presented below for both the 2D and 3D cases, in order to provide both a comparison with the 2D numerical simulations here and an estimate for expectations for 3D simulations. As comparison cases, in addition to the expanding grid, we also consider a uniform ne grid with a grid spacing of meters and a time step of . seconds, and a uniform coarse grid with a grid spacing of meters and a time step of seconds.
To test a model's ability to simulate both shallow and deep convection, numerical experiments are set up to simulate the transition from shallow to deep convection, similar to, e.g., [12]. We solve the equations for moist, non-precipitating atmospheric dynamics with phase changes of water between vapor and cloud water, with some details described in the appendix. To remain as simple as possible for these tests, a Boussinesq rather than anelastic setup is used [17]. The simple setup allows us to use a simple code, and the code is available on Zenodo with the DOI 10.5281/zenodo.3857736 to provide an illustration of how to modify an existing code in order to allow the capability of the expanding grid. To initialize the simulations, small amplitude random temperature perturbations are inserted near the surface. Prescribed large scale moistening and cooling are applied as forcing terms to the prognostic equations for the thermodynamic variables [5,9], and a sponge layer is added to all variables in the upper part of the domain.

Shallow and deep convection together: What is cost of ne vs. expanding grid?
Fig. 2 shows a comparison of numerical simulations for the expanding grid and the uniform ne grid. Specically, the horizontally-averaged liquid water mixing ratio is plotted, as a function of height z and time t. Both the expanding grid and ne grid show a gradual deepening in a transition from shallow convection to deep convection. This suggests that the expanding grid can simulate the same types of convective features as the ne grid. While it is clear from Fig. 1 that the expanding grid is less expensive than the uniform ne grid, it is useful to have quantitative estimates of the cost savings. As theoretical estimates of cost savings, we nd where N f and Nex are the total number of grid points in the ne grid and expanding grid, respectively, and T f and Tex are the total time to carry out the simulations for the ne grid and expanding grid, respectively. Details of the estimates are shown in the appendix. The approximate estimate of M− holds for both the 2D and 3D cases. Note that the cost savings depends only on M, which is the number of grid sections (i.e., the number of times the grid is re ned, plus one), and does not depend on the number of vertical grid cells per grid section. For the standard setup here with M = , the theoretical estimate of cost savings is then a factor of 32 (or, based on a more re ned estimate in the appendix, a factor of ). What we nd numerically is that it is times faster than the ne grid, which suggests that the cost savings is even greater in practice than in the crude theoretical estimates. To provide some additional tests, we also compared the sensitivity of the time savings to di erent domains by performing additional simulations on domains with horizontal extents of and kilometers; we observe that the expanding grid is and times faster than the ne grid, respectively.
Based on the considerations above, we estimate that the expanding grid is 10 to 50 times faster than traditional ne grids. The lower estimate represents an expanding grid with M = grid sections, as in the standard case considered here, with a grid spacing that expands from h = m to h = m. While such a setup was seen here to o er a speedup of roughly 45 times, that comparison was relative to a uniform grid, so we decrease the estimate by a factor of 4 to o er a comparison with a vertically stretched grid. The upper estimate of 50 times, on the other hand, represents an expanding grid with M = sections, which could be set up with a grid spacing that expands from h = . m to h = m. While the increase from M = to M = should come with an increased speedup from roughly 45 times to roughly 180 times, we again decrease the estimate (now to 50 times) to remain conservative and to o er a comparison to a stretched grid rather than a uniform grid. The choice of the number of grid sections may be somewhat empirical, and problem dependent. For instance, a domain with a large vertical extent, and therefore a wide range of scales, would bene t from a large number of grid sections, perhaps as many as . However, for domains with relatively small vertical extents, or grid sections may be su cient. In any case, the expanding grid o ers a substantial cost savings.

Re ning a deep-convection grid to resolve shallow convection: How much added expense?
Another question of interest is, How much extra expense is needed to re ne a deep-convection grid in order to resolve shallow convection? Theoretically, we expect the expanding grid to come with additional expense, beyond the coarse grid expense, by a factor of about . What we observe in numerical simulations is somewhat close to this: the expanding grid requires more computation time than the coarse grid by a factor of O( ) for the three di erent domain sizes. However, this factor of 700-1000 is small in a relative sense; if one instead re nes the uniform coarse grid to the uniform ne grid, the additional expense is theoretically a factor of roughly 33,000.
The value of the expanding grid is in its additional realism, compared with the coarse grid, and in achieving that realism for the least added expense possible. In particular, in the uniform coarse grid, cloud formation does not occur until about hours, and there is less of a gradual transition from shallow clouds to deep. See Figs. 2 and 3. Instead, there is a rapid build-up of liquid water starting at approximately kilometers in altitude, associated with a maximum vertical velocity of m/s. The cloud top then increases up until hours, when it reaches its greatest height of km. Beyond this time, the solution remains nearly constant in that the maximum vertical velocity is nearly zero and the clouds cover a constant % of the domain, with the maximum horizontally averaged liquid water staying at about . g/kg. The coarse grid is evidently very di erent than the expanding and ne grids, as the coarse grid does not resolve shallow clouds.

Which grid gives the most realistic simulation, given available resources?
In this section, we assume that there is a given amount of resources, and we ask for the grid that gives the most realistic simulations. To that end, we compare the expanding grid to the "equivalent uniform grid," or just the "equivalent grid" for short, which is de ned as the uniform grid with the same number of grid cells as  the expanding grid. This results in a uniform grid spacing of meters. As the resolution is approximately four times as coarse as the ne grid, a time step of second is used. Through this comparison, we test whether the grid arrangement of the expanding grid o ers more realism in comparison to a uniformly arranged grid.
The horizontally averaged liquid water for the equivalent grid is shown in gure 2. Broadly speaking, the equivalent grid lacks the gradual transition from shallow to deep convection that was seen in the ne grid and expanding grid. Instead, in the equivalent grid, the formation of deep convection occurs somewhat abruptly. Cloud formation rst occurs at a height of approximately km at . hours, and the transition to deep convection occurs shortly after, with the cloud top at a height of about km. Compared to the coarse grid, the equivalent grid features the same abrupt transition to deep convection, with an extremely strong updraft.

Conclusions
An expanding grid setup has been presented as a way of allowing e cient simulation of both shallow and deep convection. In particular, the expanding grid was estimated to be 10 to 50 times faster than traditional ne grids. Both theoretical and numerical estimates were presented for the cost savings, and details were described in section 3 and the appendix. In comparing the expanding grid and a ne uniform grid, a rough theoretical estimate of the cost savings is a factor of M− , where M is the number of grid sections in the expanding grid (i.e., the number of levels of re nement, plus 1). In the standard con guration tested here with M = grid sections, the cost savings is then roughly estimated to be a factor of 32. The rough theoretical estimates of cost savings are the same for either 2D or 3D and for either number of grid points or total computation time. In practice, numerical examples showed an even greater speedup of a factor of 40 to 50. Implementation of the expanding grid could be done in a number of ways. For instance, one could use AMR in either a static or dynamic/adaptive con guration. Here, as another option, we described how one can modify an existing code without needing full AMR capabilities. The expanding grid can be implemented in a straightforward way using a Fortran derived data type to store the expanding grid as a collection of regular arrays.
Some rst comparisons were shown between the expanding grid and a uniform ne grid. Due to the large cost of the uniform ne grid, a simpli ed LES/CRM in 2D was explored as a feasible setup. The comparisons show similar gradual deepening from shallow convection to deep convection in both the expanding grid and the uniform ne grid, and the maximum vertical velocity in the two cases shows similar magnitudes and variability. It would be interesting in the future to make further detailed comparisons and in 3D to better understand the similarities and di erences in variability on the expanding grid versus a uniform ne grid. It would also be interesting to consider other con gurations of the expanding grid, by trying, e.g., di erent re nement altitudes or di erent numbers of vertical grid cells for each grid section.

A Appendix
In this section, we derive theoretical estimates for the computational speedup o ered by the expanding grid. The derivation is split into two parts: rst the number of grid cells is calculated, and second the computation time is calculated by also estimating the number of time steps.
First we derive an expression for the total number of grid cells in the expanding grid. It will be a function of the total number of grid sections, M, and the number of grid cells on the coarsest section in the vertical and horizontal direction, which we denote N z ex , and N x ex . The 2D case will be considered rst for simplicity, and the 3D case is presented later below. Here, we assume that in the expanding grid, the vertical extent of each grid section is half that of the previous grid section, except the height of the lowest grid section is the same as that of the previous grid section. Grid section i then has i− N x ex N z ex grid cells for i < M, and grid section M has M N x ex N z ex cells. Summing over the total number of grid sections yields the expression where Nex is the total number of grid cells of the expanding grid. The number of grid cells in the uniform coarse grid, Nc, and uniform ne grid, N f , can be determined in terms of N x ex . If we denote the width of the domain by Lx, and the height by Lz, and assume that ∆x = ∆z, then the area of a single coarse grid cell is (Lx /N x ex ) , and the area of a single ne grid cell is [Lx /( M− N x ex )] . Dividing the total area of the domain by the area of a single grid cell yields Now we determine estimates for the amount of time each grid will take for a single time step. We assume that the coarse grid has a time step of ∆t, and the ne grid has a time step of ∆t/ M− . For a regularly spaced grid, we compute time estimates according to the formula Time = T × (number of grid cells) × (number of time steps per coarse time step), where T is a dimensional constant with units of time. (In a more sophisticated calculation, we could replace T with the function T(i), where i is the grid section number, and T(i) is the time complexity for grid section i. So theoretically, we could have di erent results if there is signi cant variation in T(i) for di erent grid sections.) With equation 5, we nd that Tc = T Nc , where Nc and N f are given by equations (3), and (4), respectively. For the expanding grid, we apply equation (5) to each grid section individually and then sum over the number of grid sections. Note that grid section i takes i− time steps per coarse time step. With this in mind, we have This simpli es to the expression In comparisons of computational cost, it is ratios that are needed. Based on the formulas above, one can determine the ratios and In arriving at these ratio formulas, one must use the fact that ∆x = ∆z within each grid section and that (It is ( / )Lz that arises, not Lz, since the height of the coarsest grid section is one-half the height of the entire domain.) Two noteworthy items are the following. First, notice that the cost speedup is roughly the same (≈ M− ) for both the ratio of the number of grid cells, N f /Nex, and the ratio of the computation times, T f /Tex. Second, notice that these ratios depend only on the number of grid sections, M, and not on the number of vertical grid cells per grid section.
The choice ∆x = ∆z was made for simplicity, but our results are actually the same for anisotropic grids with the same assumptions. We illustrate this with an example by calculating the number of grid cells on expanding and ne grids. The grid spacings for each grid section on a two-dimensional expanding grid can be written as Lx / M− N x ex and Lz / M− N z ex , and the total number of grid cells is still given by equation (2). The total number of grid points on the ne grid is then M− N x ex N z ex . Comparing this with equation (2) indicates that cost savings in the anisotropic case only depends on the number of grid sections.
We can easily generalize to three dimensions. In this case, grid section i has i− N x ex N y ex N z ex grid boxes for i < M, and grid section M has M− N x ex N y ex N z ex , where N y ex is the number of grid boxes in the y direction on the coarsest grid section. One can show that The time estimates are similar, and we arrive at the equations Based on the formulas above, one can determine the 3D ratios and The derivation is similar to the 2D case. Notice that the cost speedups are roughly the same (≈ M− ) and roughly the same as in the 2D case. In fact, the 2D ratio T f /Tex is exactly the same as the 3D ratio N d f /N d ex , which is reasonable because they both represent three-dimensional grids: (x, z, t) in the former case and (x, y, z) in the latter.
We end by noting that, thus far, we have compared the expanding grid to uniform grids. While some LES and CRM setups do use uniform grid spacings [18,25], it is also common to use a vertically stretched grid, where the horizontal grid spacings remain xed but the vertical grid spacing is stretched (e.g., it may be 50 m near the surface and stretched to be 500 m near the tropopause). So we remark brie y about the savings we could expect relative to a vertically stretched grid. We compare to the vertical structure of the grid in [13], where vertical levels are used, and the smallest vertical grid spacing is meters. Assuming that the stretched grid has a uniform horizontal resolution of meters, we nd that the three dimensional expanding grid is times faster, which is a reduction in the speedup by a factor of roughly 3 relative to the estimate of M− = that arises when comparing against a uniform grid.

B Appendix
In this section, we describe our model in more detail, and provide a discussion of some aspects of its implementation. We solve the moist, non-precipitating Boussinesq equations, which are where p is the pressure, ρ is a constant background density, ϕ = p /ρ , u(x, t) is the velocity vector, θe is the equivalent potential temperature anomaly, and q t is the anomalous total water mixing ratio. The buoyancy, b, is de ned as where θ = K is a constant background potential temperature, g = . m s − is the acceleration due to gravity and R vd = Rv /R d − = . , where R d is the gas constant for dry air and Rv is the gas constant for water vapor.
Simple parameterizations were chosen for the initial tests shown here. For example, subgrid-scale turbulence was parameterized using eddy viscosity and eddy di usion in the vertical and hyperviscosity and hyperdi usion in the horizontal. Note that the parameterizations should be adapted to the grid spacing. For example, if a di erent subgrid-scale turbulence parameterization were used, such as a Smagorinsky or turbulent kinetic energy closure, the mixing length scale depends on the grid spacing, so one would need make sure the appropriate length scales are used for each grid section. One set of physical processes that was not included is rainfall, ice, and associated cloud microphysics; instead, only phase changes between water vapor and cloud liquid water were included. Another physical process that was not included was radiative transfer. We point out, however, that a variety of approaches could potentially be used, depending on how much detail is desired for the radiation model. For instance, one could use all available grid information with a radiative transfer code that can be used with mesh re nement. As a less expensive alternative, one could average over the ne grid data and use a conventional radiative transfer model on coarsened columns. Our reason for using such a simpli ed setup was to demonstrate that basic features of multiscale convection could be reproduced in a more computationally e cient manner. It will be interesting in the future to investigate the e ects of additional physical complexity.
For the numerical methods, we used a third order upwind discretization of the advection terms, and the pressure gradient and incompessibility constraint are treated using a projection method. The Poisson equation for the pressure was solved using an iterative method (conjugate gradient method), although one could also use Fourier transforms in the horizontal direction to yield a one-dimensional di erential equation in the vertical direction, as is often done.
For passing information between di erent grid sections, the array for each variable has extra rows and columns to store ghost cells values to be used at physical boundaries and boundaries between grid sections. This is needed, for instance, for the advection terms. At grid section boundaries, the coarse grid variables are interpolated to provide these ghost cell values for the adjacent ne grid. For the Poisson equation for the pressure, rather than solving the Poisson equation over each grid section individually, we solve it over the desired grid section and all lower-altitude grid sections (albeit using a uniform grid spacing across all of those grid sections) in order to avoid the need for a pressure boundary condition at the grid section's lower boundary.
To store the data output, NetCDF or similar formats can be used, even for the expanding grid, since the expanding grid can be treated as a collection of regular arrays, and regular arrays can be stored e ciently and easily. In the data le, each variable is labelled with the grid section to which it corresponds. For example, in the expanding grid used in the numerical simulation in the paper, we have six horizontal velocity variables, each corresponding to a di erent grid section. These variables are denoted by "ui, " where i is the grid section number.