The eigenvalues and non-linear simulations of the non-hydrostatic model (NH) using the conventional forward–backward scheme (FB) (Mesinger and Arakawa, 1976; Sun, 1980, 1984), the horizontal explicit and vertical implicit scheme (HE–VI) (Klemp and Wilhelmson, 1978; Saito, 2007) and the modified non-hydrostatic model (MNH) have been discussed in Sun et al. (2012), hereafter referred to as Part I. The MNH is designed to suppress high-frequency acoustic waves by multiplying the left-hand side of the continuity equation by a parameter δ, where 16≥δ≥4. When δ=1, the MNH and NH are identical. Since the FB is applied to solve NH and MNH, they called their approach a modified forward–backward scheme (MFB) for δ>1; MFB is the same as the conventional FB when δ=1.
In Part I, frequency of acoustic waves is greatly reduced by the MFB, but much less so for gravity waves. Also, the MFB is more accurate for smaller values of δ or finer space intervals. The non-linear non-hydrostatic Cloud Resolving Storm Simulator (CReSS; Tsuboki and Sakakibara, 2002, 2007) showed that Kelvin-Helmholtz instabilities, thermal bubbles and mountain waves simulated by a conventional FB can be reproduced using an MFB with a time interval four times longer. It was also found that the shortest waves (i.e. 2Δx and/or 2Δz) require the shortest time interval according to the Courant-Fredrich-Lewy (CFL) criterion. Since the accuracy of gravity waves decreases as δ increases, it seems preferable to filter 2Δx and 2Δz waves rather than increase the value of δ. Hence, in addition to multiplying δ (δ≤16 in Part I) in the continuity equation of MNH, we propose to apply smoothing on the right hand of the continuity equation to damp the short waves. We refer to the resulting model as the modified non-hydrostatic model with smoothing (MNHS). The eigenvalue and non-linear model simulations show that, for the same accuracy, MNHS allows twice the time interval of MNH. It is also noted that other numerical schemes, for example: the modified leapfrog scheme (Sun and Sun, 2011) and/or the new semi-implicit scheme (Sun, 2011) can be applied to the wave-related terms of the MNHS to achieve the same efficiency as discussed in this paper.
The 2-D non-hydrostatic equations for the dry atmosphere can be written as:
where u and w are the x and z components of wind; p is pressure; θ is potential temperature; T is temperature; ρ is density; R is gas constant; C_{p} is specific heat at constant pressure; D_{u} and D_{w} are momentum diffusions along the x and z directions; D_{θ} is the heat diffusion. The linearised equations become:
where primes indicate deviations from the basic state variables (with subscript ‘0’), which are functions of height only, and $C=\sqrt{\left(\stackrel{{c}_{p}}{\phantom{\rule{0ex}{0.6}}}/\underset{{c}_{v}}{\phantom{\rule[--0.8]{0ex}{0ex}}}\right)R{T}_{0}}$ The basic state wind is assumed to be zero. In Part I, δ>1 was introduced in the continuity equation in order to filter high-frequency acoustic waves. Following the procedure in Part I, eqs. (7)–(11) can be formulated as:
$\stackrel{\u02c6}{p}=p\text{'}/\sqrt{{\rho}_{0}},\stackrel{\u02c6}{\theta}=\sqrt{{\rho}_{0}}\theta \text{'}/{\theta}_{0},\stackrel{\u02c6}{u}=u\text{'}\sqrt{{\rho}_{0}},$$\stackrel{\u02c6}{w}=w\text{'}\sqrt{{\rho}_{0}},{S}_{0}=\frac{\partial {\theta}_{0}}{{\theta}_{0}\partial z},\text{and}{R}_{0}=\frac{\partial {\rho}_{0}}{{\rho}_{0}\partial z}$
The linearised equations show that the 2Δx and/or 2Δz waves require the shortest Δt to satisfy the CFL criterion. Therefore, to relax the time step limitation, we apply smoothing on the right-hand side of eq. (15). Thus, eq. (15) is replaced by
where
and
is the two-dimensionally smoothed version of f. Where k and m are the wave numbers in the x and z directions, respectively. Equation (15) can now be written
and
Equations (17) and (18) indicate that 2Δx and 2Δz waves are removed from convergence/divergence. The 3Δx and 3Δz waves, which can create non-linear instability, are also significantly reduced by the second-order smoothing. The effect of smoothing decreases with increasing wavelength. The value of cos(kΔx/2)=0.809 for a 5Δx-wave seems to imply that the errors can be quite large for 5Δx and shorter waves according to eq. (17). However, linear stability analysis shows that smoothing does not significantly affect the results for wavelengths longer than 3Δx. Non-linear model simulations also reveal that the accuracy of the MNH (without smoothing) is comparable to that of the MNHS (with smooth) for the same value of δ, even when the time interval for MHNS is twice that for MNH. In Section 3, we will show the eigenvalues of the MNH, MNHS and HE–VI models, as well as their comparisons with , the eigenvalue of gravity waves derived from the differential in time and difference in space [eq. (2. 18) in Part I]:
where ${\mathrm{\Gamma}}_{0}=-\frac{g}{{C}^{2}}-\frac{{R}_{0}}{2}={S}_{0}+\frac{{R}_{0}}{2};$$X=\frac{sin(k\mathrm{\Delta}x/2)}{\mathrm{\Delta}x/2};$ and $Z=\frac{sin(m\mathrm{\Delta}z/2)}{\mathrm{\Delta}z/2}.$
The solutions of eqs. (12)–(15) can be assumed:
where ˜λ is the complex amplification factor.
Following Part I, the Arakawa C grids are used in the eigenvalue analysis; a forward in time scheme is applied to the momentum fields and a backward to the pressure and temperature, producing finite difference equations for eqs. (12)–(15),
The eigenvalues for the modified forward–backward scheme with smoothing (MFBS) can now be determined from
Without smoothing, eqs. (25) and (26) are identical to eq. (3.6) of Part I, which consists of two acoustic waves and two gravity waves. Let us define:
The amplification ∣˜λ∣ factor is then unity if Δt≤Δ˜t_{δ}. In that case, the frequency ˜ω_{δ} is defined by
and
Without smoothing, the values of ˜λ, ˜ω_{δ}, and Δ˜t_{δ} of eqs. (25)–(29) are also equal to λ_{δ}, ω_{δ} , and Δt_{δ} of MFB discussed in Part I.
With g=10 m s^{−2}, C=300 m s^{−1}, Δx=400 m, Δz=125 m, S_{o}=1.0×1.0^{−5} m^{−1} and δ=16, we obtain Δt_{δ=16}=1.59 s according to eq. (27) without smoothing. Figure 1 shows that: (a) the frequency of ω_{δ=16} (from 0.002 near the right-bottom corner to 0.009 s ^{−1} with interval of 0.001 s^{−1}) and ω_{δ=16}−σ_{g} (from −3.5×10^{−6} to −5.0×10^{−7} s^{−1}); (b) ${\stackrel{\u0303}{\omega}}_{\delta =4}$and ${\stackrel{\u0303}{\omega}}_{\delta =4}-{\sigma}_{g}$ (from –8×10^{−7} to −1.0×10^{−7} s^{−1}); and (c) ${\stackrel{\u0303}{\omega}}_{\delta =16}$ and ${\stackrel{\u0303}{\omega}}_{\delta =16}-{\sigma}_{g}$ (from −4.×10^{−6} to −5.0×10^{−7} s^{−1}). Smoothing of the continuity equation creates a large error for the 2Δx (and/or 2Δz) waves in Fig. 1b and c. Because those short-waves cannot be calculated accurately using the finite difference schemes, the eigenvalues of the unwanted 2Δx and 2Δz waves are not presented here. The error associated with smoothing decreases drastically with increasing wavelength. The errors for 3Δx waves, of –8×10^{−7} s^{−1} in Fig. 1b and −3.5×10^{−6} s^{−1} in Fig. 1c, are comparable with the errors at longer wavelengths. With approximately the same Δt (~1.6 s), the maximum value of ${\stackrel{\u0303}{\omega}}_{\delta =4}-{\sigma}_{g}$(=∣8×10^{−7}∣ s^{−1}) is much less than the maximum of ${\omega}_{\delta =16}-{\sigma}_{g}$(=∣3.5×10^{−6}∣ s^{−1}). The error of ${\stackrel{\u0303}{\omega}}_{\delta =16}$(=∣4×10^{−6}∣ s^{−1}) is also comparable with the error of ω_{δ=16}(=∣3.5×10^{−6}∣ s^{−1}), where the time interval used for ${\stackrel{\u0303}{\omega}}_{\delta =16}$(Δt=3.08s) is twice that for ω_{δ=16}(Δt=1.59s).
When the basic state remains the same and Δx=5 m and Δz=5 m, Δt_{δ=16}=4.71×10^{−2} s according to eq. (27) without smoothing. Figure 2a shows that the frequency of ω_{δ=16}(Δt=4.71×10^{−2}s) is from 0.002 to 0.009 s^{−1}, and ω_{δ=16}−σ_{g} is from −2.5×10^{−8} to 2.0×10^{−8} s^{−1}. With smoothing, δ=16, and Δt=9.42×10^{−2} s_{,}${\stackrel{\u0303}{\omega}}_{\delta =16}$ and ${\stackrel{\u0303}{\omega}}_{\delta =16}-{\sigma}_{g}$(from −3.5×10^{−9} to 5.0×10^{−10} s^{−1}) are shown in Fig. 2b. The accuracy of ${\stackrel{\u0303}{\omega}}_{\delta =16}$ is better than ω_{δ=16}, which indicates that smoothing of the continuity equation does not deteriorate the accuracy of eigenvalues in the fine resolution case. Part I showed that the MNH (with δ=16) is capable of reproducing the FB results with a time interval four times longer. When spatial smoothing is added to MNH, the allowable time interval becomes eight times that of FB, according to linearised equations, without significantly reducing the accuracy. Note also that Part I showed that in a fine-resolution model with small aspect ratio, the MFB is more accurate and efficient than HE–VI.
The non-linear non-hydrostatic CReSS (Tsuboki and Sakakibara, 2002, 2007), has been applied to simulate shear instability, thermal bubbles, and mountain waves with the FB, MFB, MFBS and HE–VI. CReSS uses Arakawa C and Lorenz staggered grids in the horizontal and vertical, respectively. Prognostic variables consist of the 3-D velocity components, perturbations of pressure and potential temperature, water vapour mixing ratio, subgrid scale turbulent kinetic energy (TKE) and cloud physical variables. In the FB version of CReSS, u, v and w are calculated using forward-differences, while the pressure perturbation p’ uses the backward-difference with a small time interval (Δts). The non-linear and other forcing terms, which are calculated at a larger time interval (Δtb=nΔts), remain constant when integrating sound waves. The leapfrog scheme is applied to the advection terms. The model also includes viscosity and divergence damping of the pressure gradient force in the momentum equations. In the HE–VI, the forward difference is applied to calculate u and v. Then, the implicit scheme is applied to solve w and p’ using a tri-diagonal matrix solver, as discussed in Saito (2007). The detailed equations and numerical schemes are referred to Tsuboki and Sakakibara (2007). There is no analytical solution to the non-linear eqs. (1)–(6). Hence, the results of the conventional FB (δ=1) will be used as the references to compare with the simulations from the HE–VI, MFB and the new MFBS.
The initial x-component wind u is shown in Fig. 3a. The initial background potential temperature and elliptic-type perturbation with the amplitude of –0.4 K, and the horizontal and vertical radii of 180 m and 60 m are shown in Fig. 3b. Figure 4a shows the simulated potential temperatures θ at t=320 s with Δtb=0.4 s from: FB with Δts=0.01 s; MFB with δ=16, Δts=0.04 s; MFBS with δ=4, Δts=0.04 s; MFBS with δ=16, Δts=0.08 s; and HE–VI with Δts=0.02 s. The five simulations nearly coincide. The differences between θ_{MFBS}(δ=4, Δts=0.04 s) and θ_{FB} (δ=1, Δts=0.01 s) range from –0.03 to 0.03 K, (Fig. 4b). Meanwhile, θ_{MFBS}(δ=16, Δts=0.08 s)–θ_{FB} (δ=1, Δts =0.01 s) is from –0.12 to 0.08 K, (Fig. 4c). Fig. 5b–d of Part I show the differences between the MFB (δ=16, Δts=0.04 s) and FB, which range from –0.12 K to 0.1 K; between the HE–VI (Δts=0.02 s) and FB they are from –0.02 to 0.03 K. The accuracies are comparable between the MFBS with δ=4, Δts=0.04 s and HE–VI with Δts=0.02 s, as well as between the MFB with δ=16, Δts=0.04 s and MFBS with δ=16, Δts =0.08 s. However, the Δts for the MFBS is twice that for the MFB with the same δ. It is also noted that the HE–VI model becomes unstable when Δts=0.04 s.
The initial x-component wind and the background potential temperature are the same as in Section 4.1. The amplitude of initial potential temperature perturbation θ’ is –0.4K, which consists of sine-waves in the horizontal with a wavelength of 40 m and an elliptic-type perturbation in the vertical with a radius of 120 m, as shown in Fig. 5a. Figure 5b shows the simulated θ’ at t =240 s with Δtb=0.12 s from FB with δ=1 and Δts=0.005 s (solid line) and from MFBS with δ=16 and Δts=0.04 s (dashed line). The temperature contours are virtually indistinguishable. The same holds for the velocity fields (not shown). However, the simulated temperature and velocity from HE–VI (Δts=0.005 s) depart significantly from the FB with the same time interval, as discussed in Part I.
The initial potential temperature of a thermal bubble with a Gaussian profile is shown in Fig. 8a of Part I, following eqs. (38) and (39) of Robert (1993). Simulations were performed with the same Δx=Δz=5 m and Δtb=0.12 s as in Part I. The simulated θ and velocity vector at t=6, and 12 min of MFBS (δ=16 and Δts=0.06 s) are almost identical to those of FB (δ=1 and Δts=0.008 s) shown in Fig. 8b and c in Part I. Figure 6a shows the simulated θ_{MFBS} (Δts=0.06 s with black contours) and θ_{FB} (Δts=0.008 s with colour shaded plot) at t=18 min for comparison. The difference between θ_{MFBS} and θ_{FB} ranges from –0.03 to 0.04 (Fig. 6b). Between θ_{MFB} (δ= 16, Δts=0.03 s) and θ_{FB} the difference ranges from –0.03 to 0.06 K, and for θ_{HE–VI} (Δts=0.008 s) and θ_{FB, } from –0.3 to 0.2 K, as shown in Fig. 9c and d in Part I. The time interval Δts for the MFBS is twice that for MFB and almost eight times that for HE–VI. Yet, the simulations using MFBS are as accurate as MFB and better than HE–VI. The simulations are also comparable to those of Robert (1993), Hsu and Sun (2001) and Chen and Sun (2001), etc., although the details depend upon the particular numerical schemes and turbulence parameterisations chosen for each model.
The initial conditions are identical to those described in Section 4.3, except for the addition of a 10 ms^{−1} prevailing mean wind. This allows for an interaction between thermal convection and the mean wind, which is a more realistic scenario than the previous case. The simulated potential temperature perturbation and velocity perturbation V’ from the mean wind for the MFBS (δ=16 and Δts=0.06 s) at t=12 and 24 minutes are shown in Fig. 7a and b with a periodic boundary condition. The pattern is almost symmetric at t=12 min and becomes more asymmetric with increasing time, because the x-component wind is slightly stronger on the windward side of the thermal. Figure 7b shows that the warm area shifts to the right-hand side of the front at t=24 minutes. Differences from the FB (Δts=0.008 s) simulations remain small: between –0.004 and 0.006 K at t=12 minutes, and between –0.15 and 0.2 K at t=24 minutes (Fig. 7c and d).
A uniform 10 m s^{−1} wind flows over a bell-shaped mountain with a peak height of 500 m and a half-width of 2000 m. The background buoyancy frequency is 0.01 s ^{−1}. The simulated w with Δx=400 m, Δz=125 m and Δtb=10 s, at t=9000 s from: FB with δ=1, Δts=0.25 s; MFB with δ=16, Δts=1.0 s; MFBS with smoothing, δ=16, Δts=2.0 s; and HE–VI with Δts=1.0 s are shown in Fig. 8a. They coincide among each other. Figurw 8b shows that the difference in simulated w between the FB (Δts=0.25 s) and MFBS (δ=16, Δts=2.0 s) at t=9000 s ranges from –0.1 to 0.08 m s^{−1} (solid lines), and the difference between the FB and MFB (δ=16, Δts=1.0 s) is also between –0.1 and 0.08 m s^{−1} (dashed lines). Similarly, the difference between the FB and HE–VI (Δts=1.0 s) is between –0.1 m s^{−1} and 0.08 m s^{−1} (dashed line) in Fig. 8c. Hence, the accuracy is comparable among the MFB (δ=16, Δts=1.0 s), MFBS (δ=16, Δts=2.0 s with smoothing) and HE–VI (Δts=1.0 s).
In Part I, it was demonstrated that the MNH with a parameter δ (between 4 and 16) applied to the continuity equation can suppress high-frequency acoustic waves effectively without a significant impact on gravity waves, enabling the use of a longer time step. Here, we have added a smoothing to the right-hand side of the continuity equation in the MNH, allowing us to double once more the time interval which was used for the MFB in Part I. Both eigenvalue analyses and non-linear model simulations show that the smoothing does not significantly degrade the accuracy of the results. The MNHS approach is simple and can be easily incorporated into many different numerical schemes. Furthermore, the results do not depart from the conventional FB result, even when the time interval used for the MNHS is eight times that for FB. However, HE–VI simulations can be significantly different from FB even when using an identical time interval that is allowed by the CFL criterion, as discussed in Part I. When the aspect ratio is around 6, the same Δts can be applied to the MFBS and HE–VI, but the HE–VI still needs more computing time to solve a tri-diagonal matrix. Besides, the HE–VI results may depart from that of FB, as discussed in Part I. We conclude that the MNHS is more accurate and efficient than the HE–VI when the aspect ratio is small, as in the cases presented here.
We would like to thank Profs. Uyeda and Ohigashi, Kato San and Tsujino San at Nagoya University for their help when the senior author was a Visiting Professor at the Hydrospheric Atmospheric Research Center (HyARC), Nagoya University. We also like to thank Mr. Seaman at Purdue Information Technology Center for his assistance during this study. The reviewers’ comments are greatly appreciated. The computing facilities provided by Nagoya University, Purdue University and Taiwan High Performance Computing Center are highly appreciated as well.
ChenS. H, SunW. Y. The applications of the multigrid method and a flexible hybrid coordinate in a nonhydrostatic model. Mon. Weather Rev. 2001; 129: 2660–2676.
HsuW.-R, SunW. Y. A time-split, forward-backward numerical model for solving a nonhydrostatic and compressible system of equations. Tellus. 2001; 53A: 279–299.
KlempJ. B, WilhelmsonR. B. The simulation of three-dimensional convective storm dynamics. J. Atmos. Sci. 1978; 35: 1070–1096.
MesingerF, ArakawaA. Numerical Methods used in Atmospheric Models. 1976; , Paris, France64. GARP Publication Series No. 14, WMO/ICSU Joint Organizing Committee.
RobertA. Bubble convection experiment with a semi-implicit formation of the Euler equations. J. Atmos. Sci. 1993; 50: 1865–1973.
SaitoK. Nonhydrostatic atmospheric models and operational development at JMA. J. Meteo. Soc. Japan. 2007; 85B: 271–304.
SunW. Y. A forward-backward time integration scheme to treat internal gravity waves. Mon. Weather Rev. 1980; 108: 402–407.
SunW. Y. Numerical analysis for hydrostatic and nonhydrostatic equations of inertial-internal gravity waves. Mon. Weather Rev. 1984; 112: 259–268.
SunW. Y. Instability in leapfrog and forward-backward schemes: part II: numerical simulation of dam break. J. Comput. Fluids. 2011; 45: 70–76.
SunW. Y, SunO. M. T. A modified leapfrog scheme for shallow water equations. J. Comput. Fluids. 2011; 52: 69–72.
SunW. Y, SunO. M. T, TsubokiK. A modified atmospheric nonhydrostatic model on low aspect ratio grids. Tellus A. 2012; 64: 17516.
TsubokiK, SakakibaraA, ZimaH. P, JoeK, SatoM, SeoY, ShimasakiM. Large-scale parallel computing of Cloud Resolving Storm Simulator. High Performance Computing. 2002; Berlin, Germany: Springer Verlag. 243–259.
TsubokiK, SakakibaraA. Numerical prediction of high-impact weather systems. The Seventeenth IHP Training Course (International Hydrological Program). 2007; December. 2-15Nagoya University and UNESCO, Nagoya, Japan273.