Influence of grid aspect ratio on planetary boundary layer turbulence in large-eddy simulations

We examine the influence of the grid aspect ratio of horizontal to vertical grid spacing on turbulence in the planetary boundary layer (PBL) in a large-eddy simulation (LES). In order to clarify and distinguish them from other artificial effects caused by numerical schemes, we used a fully compressible meteorological LES model with a fully explicit scheme of temporal integration. The influences are investigated with a series of sensitivity tests with parameter sweeps of spatial resolution and grid aspect ratio. We confirmed that the mixing length of the eddy viscosity and diffusion due to sub-grid-scale turbulence plays an essential role in reproducing the theoretical − 5/3 slope of the energy spectrum. If we define the filter length in LES modeling based on consideration of the numerical scheme, and introduce a corrective factor for the grid aspect ratio into the mixing length, the theoretical slope of the energy spectrum can be obtained; otherwise, spurious energy piling appears at high wave numbers. We also found that the grid aspect ratio has influence on the turbulent statistics, especially the skewness of the vertical velocity near the top of the PBL, which becomes spuriously large with large aspect ratio, even if a reasonable spectrum is obtained.


Introduction
In meteorological simulations, the grid aspect ratio a of horizontal ∆x to vertical grid spacing ∆z 15 is generally much larger than that in other fluid dynamics fields. Here, we define the aspect ratio as a = ∆x/∆z. The use of such large aspect ratios has been validated based on the phenomena studied in this field. So far, most of the phenomena of interest are considerably affected by the rotation of the earth and vertical stratification; both lead to differences between the horizontal and vertical scales of atmospheric phenomena. Further, since it is natural to treat the atmospheric motions 20 separately in the horizontal and vertical directions in order to understand the atmospheric dynamics, a large aspect ratio is used for grids in meteorological simulations. We will have a brief review of the atmospheric phenomena of various scales and the grid aspect ratio used in numerical simulations for the phenomena in the previous studies. This background would help us to understand the grid configuration in current meteorological large-eddy simulations (LESs) and the importance of the 25 investigation of the effect of the grid aspect ratio on the turbulence in LES.
Atmospheric phenomena have a wide range in terms of both time and space. One of the initial targets of the studies involving meteorological numerical simulations was the planetary wave. In this case, the horizontal scale, at O(1000)-O(10 000) km, is much larger than the vertical, and is dominated by the effect of rotation. The vertical scale is mainly determined by atmospheric thick- If we turn our attention to mesoscale phenomena, whose horizontal scales are O(10)-O(100) km, a more detailed expression of vertical structure is required. Since moist convection plays an important role in the vertical motion of such phenomena, hydrostatic balance, which is a good approxi-40 mation for planetary and synoptic-scale phenomena, is no longer satisfied. The spatial scale of their vertical motion is roughly O(1) km. Therefore, the grid aspect ratio in this case should be reasonably determined with O(10)-O(100). In the beginning of the 1970s, one of the most significant paradigm changes in meteorological numerical studies was the appearance of non-hydrostatic models (e.g., Klemp and Wilhelmson, 1978). Cumulus convection is explicitly represented in these models (e.g., 45 Skamarock et al., 2008;Hodur, 1997;Tripoli, 1982;Xue et al., 2000;Black, 1994;Saito et al., 2006;Tsuboki and Sakakibara, 2002). Recently, global convection system-resolving models, such as NICAM (Tomita and Satoh, 2004;Satoh et al., 2014), have been developed to investigate the interactions between synoptic-scale phenomena, cloud clusters, and individual cumulus clouds. To resolve vertical motion precisely, the number of vertical layers in such models is larger than in gen-50 eral circulation models (GCMs), but the grid aspect ratio is still large. For example, it is five in Miyamoto et al. (2013Miyamoto et al. ( , 2015.
In the atmosphere, there are many smaller scale variabilities than associated with cumulus clouds.
For example, the spatial scale of the dominant motion in the unstable planetary boundary layer (PBL) is smaller than 1 km. Clearly, turbulence in the PBL is an important issue for all meteorological 55 models, because heat and mass transport by turbulence in both the horizontal and vertical directions strongly affect the atmospheric mass and energy balances. In particular, vertical transport is important for atmospheric variability over the PBL. Its effect is parameterized as turbulence models in GCMs and cloud-resolving or permitting models (e.g., Mellor and Yamada, 1982;Nakanishi and Niino, 2004). At this scale, the effects of rotation and stratification are smaller compared with those schemes in Appendix B, and the validity of the dynamical core in Appendix C. Symbols used in this paper are summarized in Table 1.

100
Since the purpose of this study is to clarify the impact of the grid aspect ratio on the turbulent aspects, the influences of the approximations to the governing equations should be reduced as much as possible. For this reason, we employ the set of fully compressible non-hydrostatic equations as the governing equations. The three-dimensional momentum (ρu, ρv, and ρw, for the x, y, and z directions, respectively), total density (ρ), mass-weighted potential temperature (ρθ), and mass 105 concentration of tracers (ρq x ; x ∈ {v, l, s}) are used as the prognostic variables. See Appendix A1 for the detailed formulation of the governing equations.
The central difference schemes are used for spatial discretization; the second-order scheme is applied to the pressure gradient terms in the momentum equations and divergence terms of mass flux in the continuity equation, while the fourth-order scheme is applied to the advection terms in the 110 momentum and thermodynamics equations. The reason the advection terms are discretized by the high-order scheme is based on the accuracy of the eddy viscosity and the diffusion terms representing effect of SGS turbulence. The coefficients of the viscosity and diffusion terms are proportional to the square of the grid spacing, so that the magnitude of the terms would be comparable to the truncation error of the advection terms, in terms of order of accuracy, if the second-order scheme 115 is employed. Additionally, the higher-order treatment for the advection terms is necessary from a different viewpoint as well. Since the advective term is a non-linear convolution, it requires a higherorder treatment to resolve the additional modes. The use of the lower-order scheme is justified by the scale-separation of the fast modes (acoustic and fast gravity waves) and slow modes (advection). In the meteorological phenomena, the terms of the pressure gradient in the momentum equations and 120 the divergence in the continuity equation are dominant for the fast modes, while the advection term is dominant for the slow modes. The interaction between the fast and slow modes is not significant generally. If we consider SGS mixing in a local field such as a several-grid scale, the fast waves would pass over this field soon before completing the SGS mixing. This means that fast waves do not participate much in the local mixing, compared with the mixing process itself.

125
The even-order schemes require the explicit numerical filter for numerical stability. However, they have the advantage that it is easier to evaluate the effect of the numerical diffusivity using the explicit numerical filter than that introduced implicitly by an odd-order scheme. See Appendices A3 and A4 for details of the discretization and a discussion of the numerical filter, respectively.
A fully explicit scheme, i.e., HE-VE (horizontally explicit and vertically explicit) scheme, is used 130 for temporal integration in this study. This scheme generally has less implicit numerical diffusion than implicit schemes such as HI-VI (horizontally implicit and vertically implicit) and HE-VI (horizontally explicit and vertically implicit) schemes. To focus on the influences of the grid aspect ratio, an explicit scheme is more suitable than implicit schemes. Additionally, we do not use a time-split scheme to avoid numerical damping of the time-splitting. See Appendix A2 for a detailed description 135 of the temporal scheme used in this study.
In order to validate the dynamical core of this model, we performed a density current experiment (Straka et al., 1993) as a standard test case. The detailed results are described in Appendix C; we posit that the dynamical core of the model shows reasonable performance, and is reliable enough for the investigations performed in this study. As well as validating the dynamical core, Sato et al. 140 (2014) validated this model from the viewpoint of moist physical processes.

Overview of LES modeling
The model for effect of SGS turbulence used in this model is a Smagorinsky-Lilly-type model (Smagorinsky, 1963;Lilly, 1962). In the governing equations shown in Appendix A1, the effect of 145 the SGS turbulence can be written as ∂ρτij ∂xj in the momentum equations, and ∂ρτ * i ∂xi in the thermodynamics and tracer equations. In this subsection, the subscripts imply summation over the set {1, 2, 3}.
The τ ij and τ * i are parameterized as

150
where ν SGS and ν * SGS are the coefficients of the SGS eddy viscosity and diffusion, respectively, δ ij is the Kronecker delta, and φ represents scalar quantities such as θ and q x . S is the strain tensor given by S ij = 1 2 TKE SGS is the turbulence kinetic energy (TKE) of the SGS turbulence: The second term on the right-hand side of Eq. (1) arose because of the consistency with Eq. (4) in the compressive flow (e.g., Moin et al., 1991).
The role of the sub-grid model is to parameterize the effect of SGS turbulence based on the energy cascade theory of three-dimensional isotropic turbulence. The eddy viscosity and diffusion model is 160 employed as a sub-grid model to represent the effect. For the determination of the amount of the energy cascade, the mixing length of the eddy viscosity and diffusion is the most important factor.
The coefficient of the SGS eddy viscosity, ν SGS in Eq. (1), is proportional to the square of the mixing length (L mix ) in Prandtl's mixing length theory;

165
The mixing length depends on what type of spatial filter we employ on the variables in the equations, and the length scale of the spatial filter is the essential factor for the mixing length. The spatial filtering is inevitable in the discretization of the equations. The spatial filter in a numerical model is implicitly determined by the grid spacing and discretization schemes. The artificial length characterizing the spatial filter owing to discretization is defined as the filter length. Besides the filter length,

170
the shape of the grid should be also considered in the mixing length determination. Its effect can be represented by the grid aspect ratio. Scotti et al. (1993) proposed an equation for the effect of the grid aspect ratio on the mixing length. Using the filter length ∆ and grid aspect ratio a, the mixing length can be represented as

175
where f (a) is a function of the grid aspect ratio, and represents the effect of the grid aspect ratio on the mixing length. C s is the Smagorinsky constant, and F S represents the effect of static stability.
The effect of the static stability on the mixing length is introduced in this model according to Brown et al. (1994). They extended the Lilly model to improve the performance of the simulation by considering the dependency of the Prandtl number, Pr = ν SGS /ν * SGS . For the unstable condition 180 (Ri < 0), where Ri is the local (pointwise) gradient Richardson number, defined as

185
Pr N is the Prandtl number in neutral conditions; Pr N = 0.7. The constants c and b are 16 and 40, respectively, following Brown et al. (1994). They are empirically given to fit surface-layer observations. For the weakly stable condition, where Ri is smaller than the critical Richardson number Ri c (= 0.25), i.e., 0 ≤ Ri < Ri c , For the strongly stable condition (Ri ≥ Ri c ), The TKE SGS in Eq. (1) is where C k is an SGS constant and assumed to be 0.1, following Deardorff (1980) and Moeng and Wyngaard (1988).

Problems of filter length and grid aspect ratio
There are two ambiguous factors in the configuration for determining the mixing length in recent 200 meteorological LES models. One problem is the configuration of the filter length in Eq. (6) in the models. The value of ∆ is usually set as the grid spacing for simplicity, not by considering the numerical filter inherited by an individual model scheme. The other problem is that the effect of the grid aspect ratio on the mixing length is not considered in determining the mixing length; correspondingly, f (a) = 1 in Eq. (6).

205
Such rough treatment of the filter length and no consideration of the grid aspect ratio lead incorrect turbulent properties. A usual remedy to such effects is a posteriori tuning of the Smagorinsky constant, C s in Eq. (6), so as to reproduce realistic results (e.g., Deardorff, 1971 , 1991). They are advantageous in the case in which the assumption of the isotropic turbulence is not justified, e.g., very close to a boundary. However, in our opinion, dynamic SGS models are associated with some practical problems, e.g., numerical stability; more importantly, they seem to be a type of a mathematical procedure and do not seem to have a physical basis.
Here, we should note that each theory is based on its own basic concepts. Although this means that 215 different theoretical concepts lead to different values for the constants, we should keep in mind that a certain constant is uniquely determined by a particular model on theoretical basis. The Smagorinsky constant is derived by integration of the kinetic energy filtered out by the spatial filter. If the cutoff filter is employed as the spatial filter, the integration can be performed in the cubic B in Fourier space: where k i is the wavenumber in the i-direction. Lilly (1967) obtained the constant (C s = 0.16) by integrating in the sphere inscribed in the cubic B, while Scotti et al. (1993) obtained it more precisely by integrating in the entire cubic. In this study, we follow Scotti et al. (1993), and use C s = 0.13.
Let us return to the first problem, i.e., the determination of the spatial filter length ∆. Most models 225 in previous studies assume that the filter length is determined as the geometric mean of the grid spacings in the three directions (e.g., Deardorff, 1980): It is generally difficult to strictly define the spatial filter in a model because the filtering effect implicitly introduced by discretization is very complicated. This difficulty is especially significant if 230 we apply the implicit method of time integration. On the other hand, the effect of the explicitly introduced numerical filter can be more easily estimated than that of an implicitly introduced filter. In our model, we introduce hyperviscosity and diffusion explicitly as the numerical filter (see Appendix A4). The explicit numerical filter is thought to be the dominant spatial filter because we employ the explicit temporal integration method and even-order spatial difference scheme. This type 235 of filter removes the two-grid scale noise selectively. In other words, the two-grid scale variability is filtered out from the resolved variability. Consequently, we can say that in our model, the two-grid scale is preferred as the filter length rather than the grid spacing defined in Eq. (16): In the next section, we investigate the validity of this configuration using simulation results.

240
The second problem is the treatment of the grid aspect ratio in LES in actual grid systems. In the equilibrium condition with the universal Kolmogorov spectrum, the energy flux cascaded into the SGS variability is equal to the SGS dissipation. Since the dissipation does not depend on the artificial grid configuration but on the physical configuration, energy cascaded to the SGS turbulence should not depend on grid configuration, including the grid aspect ratio. The energy cascade flux, which is 245 equal to the dissipation, can be written as a function of the strain tensor S ij defined in Eq.
(3) and the mixing length in Eq. (6) in Smagorinsky-type models. S ij depends on the grid configuration, including aspect ratio. The mixing length should be thus dependent on the aspect ratio to cancel out the dependency of S ij on the aspect ratio, otherwise the energy cascade flux and dissipation should be dependent on the grid configuration. Scotti et al. (1993) derived an approximate function of the where b 1 = arctan(1/a), b 2 = arctan(a) = π/2 − b 1 , and 255 P 1 (z) = 2.5P 2 (z) − 1.5(cos(z)) 2/3 sin(z), Note that they considered two grid aspect ratios, while here we consider only one for simplicity, i.e., ∆x = ∆y > ∆z. Note that f (a) is a monotonic function for a, for instance, f (2) = 1.036, f (5) = 260 1.231, f (10) = 1.469, and f (20) = 1.790. Roughly speaking, there can exist larger scale (or lower wavenumber) variability in filtered out variabilities for a larger grid aspect ratio. This implies larger turbulence kinetic energy for larger aspect ratios, and larger TKE results in a larger mixing length.
In our study, the effect of the grid aspect ratio on the mixing length is introduced in Eq. (18) by following Scotti et al. (1993).

Numerical experiments
We performed three PBL experiments to examine the influences of the grid aspect ratio and filter length on simulated turbulence, which are summarized in Table 2. In the control experiment, the filter length is double the grid spacing, as in Eq. (17), and the mixing length in the sub-grid model is modified by the grid aspect ratio, as in Eq. (18). The second experiment is a small filter length exper-270 iment in which the filter length is set as in Eq. (16). The last is a fixed mixing length experiment, in which the mixing length in the SGS model is not modified by the grid aspect ratio, that is, f (a) = 1.
In the most current meteorological simulations, the configuration of the filter length follows Eq. (16) and f (a) = 1 in Eq. (6).
In the three experiments above, a systematic parameter sweep of resolution and grid aspect ratios 275 was conducted. The spatial resolution of each run performed in the control experiment is shown in Table 3. The parameter set in the small filter length experiment is the same as the control experiment, except for the 5mAR10 and 5mAR20 runs, while in the fixed mixing length experiment it is the same as the control, except for the 10mAR3 run.
In all the runs, the domain size is 9.6 km × 9.6 km in the horizontal direction, and 3 km in the  Ito et al. (2010). In this study, we consider dry conditions without moist and radiation processes.
Different temporal intervals are used for the dynamical and physical processes. We define the dynamical process as that related to fluid dynamics; the advection, pressure gradient, and gravita-295 tional force terms in the governing equations are treated as dynamical processes in this model. The other processes are physical processes in this model. In this study, the physical processes are only the surface flux for the momentum and the eddy viscosity and diffusion for SGS turbulence. Note that we treat the eddy viscosity and diffusion terms, which originated from the advection term, as the physical process for correspondence with the sub-grid turbulence model in RANS mode in this 300 model. The intervals for the dynamical process ∆t dyn are 0.006, 0.012, and 0.03 s for ∆z = 5, 10, and 30 m, respectively. The interval for the physical processes ∆t phys is ten times ∆t dyn to reduce computational resource. One concern is that use of too large a factor of ∆t phys /∆t dyn has the nonnegligible effect of an artificial acoustic wave excited by intermittent forcing added by the physical processes. In this model, the artificial wave can be reduced as described in Appendix A2.2. In a pre-305 liminary experiment, we confirmed that simulated phenomena in the runs with ∆t phys = ∆t dyn and ∆t phys = 10∆t dyn do not have essential differences in the context of this study.
We used an eighth-order numerical filter, i.e., n = 8 in Eq. (A125), with the non-dimensional coefficient γ of 2 × 10 −4 in this study. The higher order is preferred in order to limit the numerical filter to smaller scale variability, but the higher-order filter requires more expensive computational 310 costs. As discussed in Appendix A4, we confirmed that the effect of the numerical filter of order eight is less than that of the eddy diffusion and viscosity. The smaller coefficient γ is preferred as long as two-grid scale noise is prevented. The coefficient was determined by another preliminary parameter sweep experiment of the coefficient without the SGS model for the numerical filter to prevent two-grid noise.

315
4 Results and discussions

Energy spectrum
Most sub-grid models are based on the idea of an energy cascade due to three-dimensional isotropic turbulence so that the energy spectrum has the slope of k −5/3 , where k is the wavenumber. the spectrum of the horizontal velocity at the lowest layer shows a slope of almost −5/3, while that of the vertical velocity does not. These differences, due to the existence of the PBL top and bottom 340 boundaries, cause anisotropy as well as static stability (Yakhot et al., 1989;Horiuti, 1993).
In most LES models, the filter length for the sub-grid model is set to be equal to the grid spacing itself. In our model, the filter length is set to be twice the grid spacing, as described in Sect. 2. In order to examine the influence of the filter length, we performed the small filter length experiment, in which we use a conventional filter length, i.e., as in Eq. (16). Figure 1b shows the energy spectrum in 345 this experiment (only runs with ∆z = 10 m are shown). The spectra show a spurious energy profile, in which energy piles up at the higher wavenumbers. The slopes of the spectra become gentler than −5/3 because of the spurious energy. The spurious energy piling is considered to be a result of the lack of energy cascade because of SGS variability, caused by the use of a small mixing length.
We found that the situation is more significant for larger grid aspect ratios and coarser horizontal 350 resolution, as shown Fig. 1b, although the effect of the grid aspect ratio represented by Eq. (18) is taken into account in this experiment, as well as in the control experiment.
The deviation from the theoretical −5/3 slope of the energy spectrum of resolved variability is more obvious at the energy-dissipative range. However, the reduced energy is compensated by the SGS energy. On the other hand, the energy pile cannot be compensated by any SGS components.

355
Therefore, the spurious energy pile is the most important point in terms of the reproducibility of the spectrum. In order to estimate the degree of the energy piling quantitatively, we introduce an index denoted by SEP: where E(k) is the energy spectrum in each run, and the fitting coefficient A is calculated with 360 a spectrum ranging from 1/1000 to 1/100 m −1 using the highest resolution run with a grid aspect ratio of unity (10mAR1 run) in the control experiment as a reference solution. Thus, SEP indicates the maximum ratio of the energy spectrum in each run to the fitting slope of −5/3 power in the reference run. Note that the index is not influenced by the lower energy than the theoretical slope at the energy-containing and energy-dissipative scales. Figure 4 shows the dependency of the SEP 365 index on the grid aspect ratio. Since the energy spectrum is not perfectly a power law with the −5/3 slope, the index is not one, even in the 10mAR1 run. The index in the control experiment is smaller than 1.2 for all the runs. On the other hand, the indices in the small filter length experiment are larger than those in the control experiment. It tends to be larger for the larger aspect ratio. This tendency can be found for both the vertical grid resolutions of 10 and 30 m. The magnitude of the SEP is larger 370 than 1.2 for all the runs, except when the grid aspect ratio is unity. This shows that, in this model, twice the grid spacing is more suitable as the filter length rather than the grid spacing itself, and suggests that the filter length should be chosen according to the numerical schemes used in a model.
In addition, the dependency on the grid aspect ratio can clearly be seen. These results suggest that an appropriate filter length is necessary, especially in the case of a large grid aspect ratio.

375
Next, we consider the influence of f (a) in the mixing length, as in Eq. (6), on the spectrum. The mixing length is modified by the grid aspect ratio in the control experiment according to Scotti et al. (1993), as described in Sect. 2. In the fixed mixing length experiment, we fixed f (a) in Eq. (6) equal to one in order to ignore the effect of the grid aspect ratio. Figure 1c shows the spectra of the runs in the fixed mixing length experiment; only runs with larger grid aspect ratios are shown. The spectra 380 for the large aspect ratio, especially greater than ten, show the spurious energy pile, as in the case of the small filter length experiment. The SEP index is greater than 1.3 in these cases, as shown in In all the experiments, the magnitude of the SEP tends to be smaller for the larger-vertical resolution. This tendency is more apparent for the larger-aspect ratio. It is possible that the amount of energy dissipation depends on the grid configuration, although it should be identical. The larger dissipation could result in smaller total energy and consequently a smaller SEP in the coarser resolution 390 runs.

Turbulence statistics
In this subsection, we show the influence of the grid configuration of the runs on the turbulence statistics in the PBL in the control experiment. Figure 5 shows the dependency of the vertical profiles of several turbulence statistics on the spatial resolution and grid aspect ratio in the control experiment: 395 the horizontal mean potential temperature, the vertical heat flux, the variance of the vertical velocity, and the skewness of the vertical velocity. These values are averaged over the last half hour (t = 3.5-4 h). The change in the 30-min averaged kinetic energy at 500 m in t = 3-3.5 h to that in t = 3.5-4,h is small enough: only 3.4%. Consequently, we assume that the state in the last half hour is appropriate enough for the analysis. In the calculation of these values, the horizontal mean is defined as 400 the mass weighted average, ρφ ρ , where the over-line represents the mathematical horizontal average. These profiles are calculated from the resolved variability, except for the heat flux. Figure 5a gives the vertical profile of the potential temperature for all the runs. The potential temperature is almost uniform, and the atmosphere is well mixed in the PBL below 1.2 km height; the difference is about 0.05 K. On the other hand, the profile in the surface layer indicates some 405 difference according to the resolution and grid aspect ratio. Here, we define the surface layer as the layer that keeps the potential temperature gradient less than −10 −3 K m −1 . The conclusions in this paragraph do not change qualitatively for a range of the threshold. In the 10mAR1 run, this steep potential temperature gradient is reproduced in the unstable surface layer below 75 m depth. Figure 6a shows the dependency of the depth of the surface layer on the horizontal resolution. We found that the depth becomes larger for coarser horizontal resolution, even if the vertical resolution is higher. Typically, in the 5mAR20 run, the depth of the surface layer is about 140 m. The variance of the resolved variability is smaller in the 5mAR20 run than in the 10mAR1 run, as shown in Fig. 7, and the strength of mixing by the resolved motion is not a cause of this tendency. Figure 7 shows that the potential temperature and the variance of the horizontal velocity have similar vertical profiles 415 against z=z/∆ in terms of the vertical gradient. This suggests that the vertical gradient of the profile in the surface layer is determined by the filter length. Theoretically, the variance of vertical velocity should converge to a certain value with grid refine-435 ment for the following reason. The variance is approximately equal to the sum of squares of each wavenumber component as where k max andŵ(k) are the maximum horizontal wavenumber and amplitude of the vertical velocity of wavenumber k, respectively. Note that here we ignore the density variability for simplicity.

440
Under the condition that the energy spectrum |ŵ| 2 has an exponential decay of k −5/3 , the accumulated energy, which is equal to half of the variance, can be obtained analytically as where A and B are constants. Nevertheless, as shown in Fig. 6b, the variance has not yet converged in the range of resolutions. The convergence point is still a debatable issue at higher resolutions.

445
Here, we try to estimate the convergence point. We estimate the total energy for k max = ∞ from Eq. 24 in the 10mAR1 run, and it is 104% of the observed total energy in the simulation. This suggests that the convergence point of the variance would be about 1.82 m 2 s −2 , i.e., 104% of 1.75 m 2 s −2 .
The profile of skewness shows an almost linear slope in the PBL except near the surface, as shown in Fig. 5d. The skewness has a similar value around 500 m height for all the runs, while it  (2011) tried to explain the dependency of the skewness on the horizontal resolution by the SGS moment, but they showed that its effect is quantitatively not sufficient to explain the difference. We suppose 460 that localized strong upward plumes due to the larger eddy viscosity is a possible cause of the larger skewness based on the following explanation. Coarser resolution corresponds to larger filter length and consequently, larger eddy viscosity. The larger viscosity prevents small scale motions. However, the amount of heat transferred vertically should be almost the same as in other runs, because static instability becomes strong if the vertical heat flux is smaller because motion is prevented. In fact, the 465 horizontal mean vertical heat flux is almost identical in all the runs (Fig. 5b). Individual convective plumes in coarser resolution runs could be stronger, transferring more heat than in higher resolution runs. It is possible that such localized stronger upward plumes are the cause of the larger skewness in coarser resolution runs.
Residuals from logarithmic linear regression of the skewness in the horizontal resolution are rel-470 atively larger for coarser horizontal resolution: e.g., the 5mAR20 (cyan square) and the 30mAR5 (green circle) runs. Dependency of the skewness on the grid aspect ratio is one of the reasons for the residuals. The dependency on the aspect ratio can be seen in Fig. 6d. The skewness tends to be larger for larger aspect ratios. Positive residuals from the 5mAR20 run and negative ones from the 30mAR5 run in the regression would be due to large and small aspect ratios, respectively. Also, the 475 skewness of the 5mAR20 run is larger than that of the 10mAR10 run (blue × symbol) at the same horizontal resolutions (100 m). Their difference could be due to the difference in aspect ratio. The skewness in the 10mAR10 run is closer to that of the 10mAR1 run (black × symbol), than that of the 5mAR20 run, even if the vertical resolution is higher in the 5mAR20 run than in the 10mAR10 run.
This suggests that a large aspect ratio produces spurious large skewness in vertical velocity, that is, artificially strong upwind. The dependence on the grid aspect ratio seems somewhat complex, and the reason for the dependence is not clear. We assume that it is related to static stability. For higher grid aspect ratios, i.e., relatively higher vertical resolution or lower horizontal resolution, static stability could vary vertically rather than horizontally. In that case, a strong shallow stable layer could cover a relatively wider region horizontally, and strong static stability and small horizontal contrast 485 might prevent vertical motion and suppress the vertical eddy viscosity and diffusion. Under this condition, only a strong parcel can go upward against the thin strong stable layer. For the same amount of vertical heat transportation independent of grid configuration, a larger grid aspect ratio might require stronger upwind, indicating larger skewness.

490
We conducted a series of planetary boundary layer experiments to examine the influences of the aspect ratio of horizontal to vertical grid spacing on the atmospheric turbulence in a large-eddy simulation. In order to focus on the influences, we tried to avoid artificial effects as much as possible by employing the fully compressible governing equations. A fully explicit (i.e., HE-VE) temporal integration scheme and a high-order central difference scheme for spatial differentials are adopted 495 to reduce implicitly introduced numerical viscosity and diffusion by discretization.
In the model used in this paper, we considered the effect of spatial filter length and grid aspect ratio on the mixing length of the eddy viscosity and diffusion, which is a parameterization of the energy cascade to SGS variability. Explicit numerical hyperviscosity and diffusion is introduced to reduce the two-grid scale noise in this model. This can be considered a spatial filter in LES. As 500 a result, the filter length in this model is double the grid spacing, while the grid spacing is used as the length in most LESs. The effect of the grid aspect ratio on the mixing length proposed by Scotti et al. (1993) was also introduced.
If we use a reasonable filter length and introduce the grid aspect ratio effect, the energy spectra of all the runs show good correspondence with the theoretical k −5/3 slope, regardless of spatial 505 resolution and grid aspect ratio. On the other hand, the theoretical slope in the spectrum is no longer obtained, and a spurious energy pile is found if the spatial filter length is set to be equal to the grid spacing as in most meteorological LESs. The spurious energy piling is more significant for larger grid aspect ratios. A shorter filter length results in a shorter mixing length and consequently, a smaller energy cascade to SGS variability. Thus, an inappropriate filter length causes spurious energy.

510
The effect of the grid aspect ratio on the mixing length of the eddy viscosity and diffusion cannot be ignored. If the mixing length is not modified by the grid aspect ratio, spurious energy also piles at higher wavenumbers because of an insufficient energy cascade to SGS turbulence at large grid aspect ratios. In previous studies, the Smagorinsky constant was often modified as a tuning parameter to obtain the expected energy spectrum. However, the constant should be determined theoretically and 515 should not be tuned to make the energy cascade large instead of considering filter length and grid aspect ratio.
The horizontal resolution and grid aspect ratio also influence the turbulence statistics. The vertical profiles of several turbulence statistics depend on the grid configuration as follows. The depth of the unstable surface layer is larger for coarser horizontal resolution (but not vertical resolution). The In this subsection, we introduce the governing equations for the prognostic variables (ρ, ρu, ρv, ρw, ρθ, and ρq x s). The gas constant and specific heat are those for total (moist) air in the 540 thermodynamics equation and equation of state. The θ in this model is not the conventional potential temperature for dry air, but the corresponding value for total air, considering water content. These variables are spatially filtered quantities, and the Favre filter (Favre, 1983) is used for u, v, w, θ, and q x s.

A1.1 Continuity equations 545
The continuity equations for each material can be described in flux form:

550
where ∇ is the gradient vector, and DIFF [φ] is the diffusion term by SGS turbulence (see Sect. 2.2).
The sum of the ratios of their mass should be unity: The source terms for water substances should satisfy the following relationship: The sum of Eqs. (A1)-(A4) gives the continuity equation of total density: For this derivation, we use the fact that the operator DIFF [φ] is a linear operator;

A1.2 Momentum equations
The momentum equations for the gas, liquid, and solid materials are described as

570
where ⊗ represents the tensor product. The pressure is derived from the equation of state as The sum of Eqs. (A10)-(A12) gives the total momentum equation as Note that the drag forces from loading do not appear in Eq. (A14), because those terms are canceled out by the summation.

A1.3 Thermodynamics equations
The equations of internal energies are given as ∂ρq l e l ∂t + ∇ · (ρq l e l u) + ∂ρq l e l w l ∂z where e x are the internal energies, and defined as The DIFF operator represents the mixing by SGS turbulence. The SGS eddy viscosity should be used for conserved quantities following the motion. The internal energies are not conserved quantities in 590 the Lagrangian sense, and the diffusion term in these equations and T * are only conceptual, and should be determined to be consistent with the diffusion term in an equation for a conserved quantity along the flow trajectory (i.e., potential temperature), which is discussed later.
The sum of Eqs. (A15)-(A17) gives the total internal energy equations: and the total diabatic heating is described as Equations (

A1.4 Conservation of thermodynamics in the dynamical process
Equation (A22) is not a complete flux form, because the internal energy itself is not conserved both in the Euler sense and the Lagrangian sense. In this section, we consider the conserved quantity for the thermodynamics equation.

605
The potential temperature for dry air, which is defined as is used as a conserved quantity in traditional models. Although it is conserved along a Lagrangian trajectory, it is no longer conserved when the water substances are included. We introduce a new conserved quantity following the motion in moist conditions.

625
where θ v is the potential temperature for water vapor, defined as Thus, The quantity is conserved along the flow trajectory, where Here we define a new potential temperature This θ satisfies and θ is a conserved quantity along the flow trajectory, even in moist conditions. We employ ρθ for the prognostic variable.

640
The pressure expression is derived diagnostically as: (A38) Figure A1a gives the vertical profile of the temperature in the US control atmosphere, and Fig. A1b shows the vertical profiles of θ/θ d under this temperature condition when we assume that q v is the specific humidity at saturation, and q l + q s gives 0.0, 0.01, 0.02, and 0.04. The difference between θ 645 and θ d becomes larger with height and may not be negligible.

A1.5 Diabatic heating in the physical process
Changing the prognostic variable for thermodynamics from the internal energy to the newly defined potential temperature θ should modify the diabatic heating and sedimentation terms in Eq. (A22).

A1.6 Summary of governing equations
The governing equations are summarized as follows: at the top and bottom boundaries.

A2 Temporal integration scheme 680
We conceptually separate the complete set of governing equations into dynamical and physical parts: The diabatic heating process, diffusion by SGS turbulence, sedimentation process of liquid and solid waters, and the source and sink of water substances are treated as physical processes, and the others are treated as dynamical processes.

685
According to this scheme, the dynamical processes can be written as On the other hand, the physical processes are as follows. The governing equations for the physical processes being

A2.1 Dynamical processes
A Runge-Kutta (RK) scheme with three steps is used as the temporal integration scheme for the dynamical processes. The RK scheme with three steps used in this model is defined as where f (φ) ≡ ∂φ ∂t . Taylor expansion of the φ t+∆t calculated by Eq. (A72) around φ t yields Theoretically, the Taylor expansion of φ(t + ∆t) is This shows that the three-step scheme used in this model only has second-order accuracy, despite the three steps. However, it has third-order accuracy for linear equations, because is a linear function. This scheme can calculate small-amplitude perturbations with third-order accuracy, and finite amplitude waves with second-order accuracy. In terms of numerical stability, this 720 scheme has a stability region, which is almost the same as in standard explicit RK schemes with three steps for small perturbations. The stability region is wider than for RK schemes with secondorder accuracy, and includes the imaginary axis, which corresponds to neutral rotating modes, while that for second-order schemes does not. This is why we choose the three-step method instead of the two-step method.

725
The advantage of this scheme compared with ordinary three-step schemes is the reduction of memory load and storage, which is one of the most expensive components of recent computers with low byte per flop (B/F) ratios, along with the benefit of higher numerical stability.

A2.2 Physical processes
The acoustic wave is the fastest mode in the dynamical processes, and the temporal interval for 730 dynamical processes must be less than the grid spacing divided by the speed of the acoustic wave, to satisfy the Courant-Friedrichs-Lewy (CFL) condition. However, the time scale of the physical processes is usually much longer than the interval, so the temporal interval for the physical processes can be much longer than for the dynamical processes. We use a larger temporal interval to calculate the tendencies of the physical processes than of the dynamical processes. We call the time step for 735 the physical processes a large time step, and for dynamical processes, a small time step.
Traditionally, tendencies in physical processes are calculated with large time steps, and the prognostic variables are updated with the Euler scheme with the tendency for large time steps. This causes an artificial acoustic wave, as described below. In this model, the tendency in some physical processes, such as surface flux (sensible heat flux and latent heat flux) and eddy viscosity and dif-740 fusion, are calculated with large time steps in the same way, but they are added to the prognostic variables with the tendency of dynamical processes in small time steps. Figure A2 shows the horizontal averaged vertical velocity in the planetary boundary layer experiment, whose details are described in Sect. 3. In the control run, in which the temporal interval for the large time step (∆t phys ) is the same as that for the small time step (∆t dyn ), the velocity is almost 745 zero, as shown in Fig. A2a. However, in the case where ∆t phys > ∆t dyn and the tendency of the surface flux is added to the prognostic variables at the large time step as in the traditional manner, a vertically propagating artificial wave can be seen, as in Fig. A2b.
In addition, such a wave radiating periodically from a fixed location can cause a spurious stationary wave in the simulation output. This could lead to misinterpretation of the simulation results, 750 although this problem is not directly based on physics or modeling. In most practical cases, the temporal interval for historical output (∆t output ) is a multiple of ∆t large , and snapshots of physical quantities are the output. In such cases, every snapshot has artificial acoustic waves with the same phase, because the sampling frequency is not sufficient. This results in a spurious stationary wave in the historical output, as we can see in Fig. A2c, that is, forcing by physical processes added at small 755 frequencies results in a spurious stationary wave in the historical data.
The artificial acoustic wave and resulting spurious stationary wave can be avoided if the prognostic variables are updated with the tendency calculated in the physical processes with a small time step, although the tendency is calculated with a large time step. Figure A2d shows that the artificial acoustic wave in (b) and (c) does not appear with this method. This shows that the tendency of the 760 physical processes can be calculated with a large time step (i.e., not every small time step), while the time integration must be done with the dynamical process at a small time step.

A3 Spatial discretization for the dynamical processes
We employ the Arakawa-C staggered grid. Central difference schemes are used for the spatial differential for the dynamical processes, because waves such as gravity waves and acoustic waves 765 generally cannot keep their spatial symmetry with odd-order schemes. Based on a consideration of numerical stability (see Appendix B), we choose the fourth-order central difference scheme for the advection (or convection) terms, and the second-order central difference scheme for the pressure gradient term in the momentum equations and divergence term in the continuity equation.
Before the discretization of the differential equations, we diagnose several quantities from the 770 prognostic variables; Full-level pressure and potential temperature Half-level density Half-level velocity

A3.1 Continuity equation
Divergence in the continuity equation is calculated with the second-order central difference scheme.
The continuity equation is discretized as

A3.2 Momentum equations
The advection terms and pressure gradient term are calculated with the fourth and second-order central difference schemes, respectively. The momentum equation is discretized as where 800 (ρu) i,j,k = −(ρu) i+ 3 2 ,j,k + 7(ρu) i+ 1 2 ,j,k 12 and the velocities at the cell wall for the staggered control volume in the x direction are defined as v i+ 1 2 ,j+ 1 In this form, the fourth-order accuracy in the advection term is guaranteed on the condition of constant velocity. The momentum equations in the y and z directions are discretized in the same way.
In the vertical equation, ρ i,j,k+ 1 2 g is added.

A3.4 Tracer advection
The tracer advection process is done with the Euler scheme after the temporal integration of the dynamical variables (ρ, ρu, ρv, ρw, and ρθ). We impose the consistency with continuity (CWC; 825 Gross et al., 2002) and monotonicity of the tracer advection following Niwa et al. (2011).
With the condition of no source/sink, the ratio of the mass of tracers to the total mass in the advection process should be conserved along the trajectory. It is, at the least, necessary that the spatially constant mass concentration should be kept in any motion of fluid. In order to satisfy this condition, we use the same mass flux at the last Runge-Kutta process (k 3 /∆t in Eq. (A72)) for 830 integration of tracers: In order to satisfy the monotonicity of tracer advection, we employ the Flux Corrected Transport 835 (FCT) scheme, which is a hybrid scheme with a higher-order difference scheme and first-order upwind scheme (Zalesak, 1979). In this model, we use the fourth-order central difference scheme as a higher-order scheme.
If the fourth-order central difference is applied, q is discretized as 840 and q high i,j+ 1 2 ,k and q high i,j,k+ 1 2 are defined in the same way. On the other hand, in the first-order upwind scheme, q is described as and q low i,j+ 1 2 ,k and q low i,j,k+ 1 2 are defined in the same way. The actual q is described as 845 and q i,j+ 1 2 ,k and q i,j,k+ 1 2 are defined in the same way. Equation (A95) can be written as where 855 F high,low i+ 1 2 ,j,k = ∆t∆y j ∆z k (ρu) i+ 1 2 ,j,k q high,low and F high,low i,j+ 1 2 ,k and F high,low i,j,k+ 1 2 are defined in the same way. The anti-diffusive flux is defined as and A i,j+ 1 2 ,k and A i,j,k+ 1 2 are defined in the same way. Equation (A99) can be rewritten as In practice, we calculate Eq. (A102) by the following steps: 1. The tentative values are calculated using the low-order flux: 2. Allowable maximum and minimum values are calculated: 3. Several values for the flux limiter are calculated: 4. The flux limiters at the cell wall are calculated: and C i,j+ 1 2 ,k and C i,j,k+ 1 2 are defined as the same way.

A3.5 Boundary condition 905
The boundary condition at the top and bottom boundaries is This leads to the boundary condition of the vertical momentum: For other prognostic variables, the vertical fluxes at the top and bottom boundaries are zero, except 920 those from physical processes.

A4 Numerical filter
We impose an explicit numerical filter using the numerical viscosity and diffusion. Although the filter is necessary for numerical stability, too strong a filter could dampen any physically meaningful variability. In this subsection, we describe the numerical filters used in this model, and discuss the 925 strength of the filter.
In order to damp the higher wavenumber component selectively, we adopt the hyperviscosity and diffusion in the traditional way. The hyperviscosity and diffusion of the nth order is defined as where f is an arbitrary variable (f ∈ ρ, u, v, w, θ, q).

930
The Laplacian of f is discretized as and Here we consider spatially dependent grid interval in calculating the Laplacian. If it is calculated with constant ∆x i as non-negligible numerical noise appears where the grid spacing varies (e.g., stretching layer near the 940 top boundary).
The hyperviscosity and diffusion can be discretized as where The coefficient, ν, is written as where γ is a non-dimensional coefficient. One-dimensional sinusoidal two-grid noise will decay to 1/e with 1/γ time steps. Note that the theoretical e-folding time is 2 n π n ∆t γ . However, it is ∆t γ with the fourth-order central scheme used in this model.
For the numerical stability of the numerical filter itself, it should satisfy for the one-dimensional two-grid noise, and for the three-dimensional two-grid noise. The conditions might be stricter for other types of noise.
The flux, F , for the numerical filter is added to the advective flux as where the first term of the right-hand side is the flux calculated by the advection scheme. In the present model, the advection scheme is the fourth-order central difference scheme (see Appendix A3). This concept is very important for the CWC condition in the tracer equations (see Ap-960 pendix A3.4). The modified mass flux of the numerical filter should be used in the tracer advection, otherwise the CWC condition is violated.
The numerical viscosity and diffusion in the y and z directions are formulated in the same way as in the x direction, although a special treatment for the z direction is needed. At the top and bottom boundaries, the flux must be zero, F kmax+ 1 2 = F kmin− 1 2 = 0. In order to calculate the F kmax− 1 2 and 965 F kmin+ 1 2 , values beyond the boundaries, f kmax+1 and f kmin−1 , are required, then the mirror boundary condition is assumed; f kmax+1 = −f kmax and f kmin−1 = −f kmin . This condition is appropriate to cause the decay the vertical two-grid noise.
Vertical profiles of density, potential temperature, and water vapor usually have significant (e.g., logarithmic) dependencies on height. Eq. (A130) has a non-zero value even for the steady state, 970 and the numerical filter produces artificial motion. To reduce this artificial motion, we introduce a reference profile which is a function of height, and deviation from the reference is used as f instead of ρ, θ, and q v in calculating the numerical filter. The reference profile can be chosen arbitrarily, but a profile under hydrostatic balance is usually chosen.
Determination of the value of the non-dimensional coefficient γ is an important issue. If it is too 975 small, the simulation could be numerically unstable or numerical noise could violate the physical variability, while variability could be dampened too much if γ is too large. At least, the effect of the numerical filter on the phenomena of the scale we want to simulate must be reasonably smaller than the effect of the physically oriented viscosity and diffusion, i.e., the eddy viscosity and diffusion representing the effect of the SGS turbulence.

980
The e-folding time of an eddy whose spatial scale is L is 2 n L n γ∆x n ∆t according to the numerical filter, and is L 2 C 2 s (2∆x) 2 |S| by eddy viscosity and diffusion for an isotropic grid under neutral stability conditions. The e-folding time of an eddy whose scale is larger than the effective resolution of the numerical filter should be longer than that for eddy viscosity and diffusion, where the effective resolution is the smallest scale of physically meaningful phenomena that a simulation can repre-985 sent. It requires that γ < 2 n+2 C 2 s |S|∆t L n−2 ∆x n−2 . The magnitude of the characteristic velocity, U , is assumed to be O(1)-O(10) m s −1 for motions whose spatial scale is approximately equal to the effective resolution. Then, |S|∆t(∼ U ∆x ∆t) is estimated to be O(10 −3 )-O(10 −2 ), because ∆x/∆t is O(350)-O(1000) m s −1 according to the CFL condition. Roughly speaking, the effective resolution is usually several times the grid spacing. Here, we assume L/∆x ∼ 2 − 10; Skamarock (2004) suggested that the effective resolution is about six times the grid spacing. It is found that γ should be This condition is stricter than the stability condition Eq. (A134) for n = 4, so this condition must be considered when the value of γ is determined if n = 4 is chosen. On the other hand, it is always 995 satisfied if the stability condition is satisfied for n ≥ 8. This means that in this case, the effect of the numerical filter on phenomena whose spatial scale is larger than the effective resolution is always smaller than that of the eddy viscosity and diffusion, and this condition for γ does not have to be considered if n ≥ 8 is chosen. The advantage of an explicit numerical filter is that its strength can be controlled or tuned based on such physical considerations.

1000
Odd-order advection schemes are generally more stable numerically than even-order schemes. In general, odd-order schemes can be divided into a central difference term and a filter term (and sometime other additional terms) conceptually. This implicit numerical filter stabilizes the calculations.
We estimate the corresponding value of γ for the implicit filter. The spatial differential of a scalar quantity φ with the ordinary third-order upwind difference scheme is written as The first term on the right-hand side of the equation represents the central difference with fourthorder accuracy. The second term, divided by |U |, is written as 1010 and corresponds to fourth-order hyperviscosity and diffusion. This can be considered as implicit numerical filter introduced to the scheme. Their dimensional coefficient is |U | ∆x 3 12 , and the nondimensional coefficient is Considering the CFL condition, The third-order scheme by Kawamura and Kuwahara (1984) is a modification of the third-order scheme to enable changing the strength of the numerical filter, and is written as fect is not always smaller than that of the eddy viscosity and diffusion representing the effect of SGS turbulence. As we discuss here, it is important to estimate the effect of the numerical filter used in a model, and to control or confirm that the artificial filter does not violate the viscosity and diffusion as a physical parameterization. This point is one of the most important factors in choosing a numerical scheme in this type of study.

A5 Physical processes
Currently, the following processes are implemented as physical processes in SCALE-LES.
-Cloud microphysics -A one-moment three-category bulk scheme (Kessler, 1969) -A one-moment six-category bulk scheme (Tomita, 2008) 1035 -A two-moment six-category bulk scheme (Seiki and Nakajima, 2014) -A bin method (Suzuki et al., 2010) -Sub-grid turbulence -A Smagorinsky-Lilly-type scheme including stability effect developed by Brown et al. (1994) 1040 -A RANS (Reynolds Averaged Navier-Stokes Simulation) turbulence model (Nakanishi and Niino, 2004, level 2.5) -Radiation -A parallel plane radiation model (MstrnX; Sekiguchi and Nakajima, 2008) -Surface flux 1045 -A Louis-type bulk model (Louis, 1979;Uno et al., 1995) -A Beljaars-type bulk model (Beljaars and Holtslag, 1991;Wilson, 2001) -Urban Canopy -A single-layer urban canopy model (Kusaka et al., 2001;Kusaka and Kimura, 2004) In this model, we use the second-order central difference scheme for the spatial differential terms of the pressure gradient and the divergence of mass flux, and the fourth-order scheme for advection terms. In this section, we investigate the numerical instability of the terms and show why we use the second-order scheme for these terms.
For simplicity, we assume a case in which the initial potential temperature is constant θ 0 , and the initial momentum is also constant. In addition, dry conditions are assumed. The potential temperature is to be constant at all times.
The governing equations are To investigate the stability, the equations are linearized. The density is divided into the basic value and the deviation from the basic: where Pressure is written as where c is the speed of the acoustic wave (= c pd c vd p00 ρ0 ) in the reference state of p = p 00 and ρ = ρ 0 . The time derivatives are written as For the second-order central difference, a = 0 and b = 1, while a = 1 and b = 9 for the fourth order.
Now we consider the numerical stability of two-grid noise. An eigen analysis showed that the two-grid noise is the most unstable eigen mode in all the cases we tested. The density with the noise can be written as ρ t i,j,k = ρ 0 + A sin(π(i + j + k)). Using the third-step RK scheme, the density after 1080 one time step, ρ t+∆t i,j,k , is ρ t+∆t i,j,k − ρ 0 = 1 − 6 (a + b) 2 (b − 3a) 2 c 2 ∆t 2 ∆x 2 ρ t i,j,k − ρ 0 .
The necessary condition for the two-grid noise to decrease in time is for the second-order spatial difference, and for the fourth order, where ∆t CFL = ∆x/c. If the fourth-order scheme is adopted, we need to make ∆t 0.6 times (or number of time steps multiplied by 1.67) larger than that in the second-order 1090 scheme.
The accuracy of the pressure gradient term and divergence term could especially affect the highfrequency modes of acoustic and gravity waves. The high-frequency modes seem to be less significant meteorologically. We choose the second-order spatial scheme for terms of the pressure gradient and divergence terms to increase ∆t in spite of the decreasing accuracy of the terms, and the fourth-1095 order scheme for the other terms.
Appendix C: Density current experiment Straka et al. (1993) proposed a standard benchmark test of a non-linear density current problem, and we performed a test experiment with the same settings. A cold bubble whose size is 4 × 2 km (horizontal and vertical, respectively) and minimum thermal perturbation is −15 K, is placed in the 1100 domain's center horizontally, and at 3 km height in a basic resting state whose potential temperature is constant at 300 K. The domain is two-dimensional, and its size is 51.2 km × 6.4 km. Runs are done with various spatial resolutions of 200, 100, 50, and 25 m. The viscosity and diffusion are given by a coefficient of 75 m 2 s −1 for the velocities and scalar quantities. Figure A3 shows a plot of the potential temperature at t =900 s for the 100 m resolution run. The structure, such as the position of the front of the density current, rolls caused by the K-H instability, and the magnitude at the local maxima, is reasonably similar with results of various models shown in Straka et al. (1993). The dependency of the L 2 norm of the potential temperature perturbation as a reference solution of the 25 m run of spatial resolution shows higher-order convergence than the second order between 100 and 50 m runs (Fig. A4). This model uses a fourth-order central difference scheme for advection, which has fourth-order accuracy with constant velocity, and second-order central difference scheme for the acoustic wave. The higher than second-order convergence seems to be due to the fourth-order difference scheme, although the velocity is not constant, and the scheme does not have fourth-order accuracy. The convergence is almost first-order between the 200 and 100 m runs. As Straka et al. (1993) mentioned, 200 m is too coarse to resolved this density current. Only spectra for ∆z = 10 m runs in the small filter length experiment and for large aspect ratio runs (a ≥ 5) in the fixed mixing length experiment are shown.  Table 3. Names of runs we performed in the control experiment. Columns and rows correspond to vertical and horizontal resolutions, respectively. The number following character "AR" represents the grid aspect ratio.