Next Article in Journal
Thermal Comfort Assessment in University Classrooms: A Discriminant Analysis for Categorizing Individuals According to Gender and Thermal Preferences
Previous Article in Journal
Recent Advances in Ni-Based Catalysts for CH4-CO2 Reforming (2013–2023)
Previous Article in Special Issue
The 50th Anniversary of the Metaphorical Butterfly Effect since Lorenz (1972): Multistability, Multiscale Predictability, and Sensitivity in Numerical Models
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

Challenges and Progress in Computational Geophysical Fluid Dynamics in Recent Decades

1
Department of Earth, Atmospheric, and Planetary Sciences, Purdue University, West Lafayette, IN 47907, USA
2
Department of Atmospheric Sciences, National Central University, Zhongli 320, Taiwan
3
Institute for Space-Earth Environmental Research, Nagoya University, Nagoya 464-8601, Japan
Atmosphere 2023, 14(9), 1324; https://doi.org/10.3390/atmos14091324
Submission received: 17 May 2023 / Revised: 27 July 2023 / Accepted: 28 July 2023 / Published: 22 August 2023

Abstract

:
Here we present the numerical methods, applications, and comparisons with observations and previous studies. It includes numerical analyses of shallow water equations, Sun’s scheme, and nonlinear model simulations of a dam break, solitary Rossby wave, and hydraulic jump without smoothing. We reproduce the longitude and transverse cloud bands in the Equator; two-day mesoscale waves in Brazil; Ekman spirals in the atmosphere and oceans, and a resonance instability at 30° from the linearized equations. The Purdue Regional Climate Model (PRCM) reproduces the explosive severe winter storms in the Western USA; lee-vortices in Taiwan; deformation of the cold front by mountains in Taiwan; flooding and drought in the USA; flooding in Asia; and the Southeast Asia monsoons. The model can correct the small-scale errors if the synoptic systems are correct. Usually, large-scale systems are more important than small-scale disturbances, and the predictability of NWP is better than the simplified dynamics models. We discuss the difference between Boussinesq fluid and the compressible fluid. The Bernoulli function in compressible atmosphere conserving the total energy, is better than the convective available potential energy (CAPE) or the Froude number, because storms can develop without CAPE, and downslope wind can form against a positive buoyancy. We also present a new terrain following coordinate, a turbulence-diffusion model in the convective boundary layer (CBL), and a new backward-integration model including turbulence mixing in the atmosphere.

Contents:
Abstract
Keywords
1: Introduction:
2: Shallow Water Equations and Numerical Schemes
    2.1: Instability Forward–Backward (FB) and Leapfrog (LF) Schemes
    2.2: Modified Leapfrog Scheme for Shallow Water Equations
    2.3: New Scheme for Shallow Water Equations
    2.4: Nonlinear Shallow Water Simulations
         2.4.1: Simulation of Dam Break:
         2.4.2: Solitary Rossby Wave
         2.4.3: Hydraulic Jumps
         2.4.4: Vortex Moving over Central Mountains Range in Taiwan
         2.4.5: Inviscid Vortex
         2.4.6: Trajectories of PV with Surface Friction
         2.4.7: Vortices Merge
3: Linearized Equations
    3.1: Coexist of Two Different Types of Cloud Bands
    3.2: Solve Linearized Equations as Initial Value Problems
         3.2.1: Rainbands and Symmetric Instability
    3.3: Cloud Bands over Tropical Continent
    3.4: Ekman layers in Atmosphere and Ocean
         3.4.1: Simulation at 30° N with Diurnal Variation in Atmosphere and Ocean
         3.4.2.: Simulation at 60° S with Diurnal Variation in Atmosphere and Ocean
4: PURDUE REGIONAL CLIMATE MODEL (PRCM)-a hydrostatic model
    4.1: Basic Equations
    4.2: Inland Sea Breeze and Dryline
    4.3: Effect of mountain
         4.3.1: Lee Vortices and Mountains
         4.3.2: Front Deforms around the Central Mountain Range in Taiwan
    4.4: Cyclogenesis and Winter Severe Storms in the USA
5: Regional Climate
    5.1: 1993 Flood in the USA
    5.2: 1988 Drought in the USA
    5.3: Snow, Land Surface, and Regional Climate
         5.3.1: One-Dimensional Snow–Soil Model
         5.3.2: Equations for Snow Layer
    5.4.: Flooding due to Snow Melt
    5.5: Southeast and East Asian Monsoon
    5.6: Simulations of Dust Storm
6: Nonhydrostatic Models
    6.1. Boussinesq Fluid Versus Compressible Fluid
    6.2. Modified Forward–Backward Scheme with Smoothing (MFBS)
    6.3: Kelvin–Helmholtz Wave
    6.4: Equations of NTU/Purdue Nonhydrostatic Model
    6.5: Numerical Scheme for Advection Equations
    6.6: Nonhydrostatic Model—Parcel Method, Froude Number, Bernoulli, Downslope Wind, and Waves
    6.7: Boulder Severe Downslope Windstorm
    6.8: Convective Available Potential Energy (CAPE)
         6.8.1: Dry Plume
         6.8.2: Moist Plume
         6.8.3: Nonhydrostatic Pressure Inside a Cloud Model
    6.9.: Lee Vortices and Hydraulic Jump in White Sand Missile Range
7: Terrain Following Coordinate in Atmospheric Model
    7.1: Equations
    7.2: New Terrain Following Coordinate
    7.3: Gal-Chen and Somerville Terrain Following Coordinate
    7.4: Numerical Simulations
         7.4.1: Gradient
         7.4.2: Divergence
         7.4.3: Curl along
         7.4.4: The Navier–Stokes Equations in New Terrain-Following Coordinate
8: Pollution Model
    8.1: Pollution in Convective Boundary Layer (CBL)
    8.2: Backward Integration of Diffusion Equation
         8.2.1: Forward-in-Time Integration
         8.2.2: Reverse-in-Time Integration
9: Summary
Acknowledgements:

1. Introduction

The foundation of fluid dynamics includes conservations of mass, momentum, and energy, as well as the thermodynamic equation of state. Fluids are assumed to obey the continuum assumption. At infinitesimal space, it contains infinite numbers of molecules and the property of the fluid (i.e., temperature, density, velocity) inside is uniform and well-defined. Fluid dynamics consists of a set of nonlinear differential equations to present the interactions among multi-length and temporal scales under different mechanical and thermodynamic environments. The computational geofluid dynamics (CGFD) has been applied to study problems ranging from small turbulence, pollution, storm, weather, to global climate, oceanic circulations, weather on other planets, and convection in the Sun. The Navier–Stokes equations (NS) do not have a well-defined size or analytical solutions. Usually, the partial differential equations are converted to the finite difference or spectrum equations at discrete point/mode, which are difficult to keep the original property of the equations. The numerical equations should be consistent and converge to the original differential equations, but it is difficult to prove that the numerical solution is unique and converges to the true solution, except a few special cases. It is also difficult to calculate the sub-grid scale turbulences/variables which cannot be resolved. Furthermore, those equations have frequently been simplified and solved by numerical methods with the help of high-performance computers. The results depend upon the detailed formulation of the forcings, the numerical schemes applied, and the initial and boundary conditions.
Bjerknes (1904) [1] formulated the primitive equations in numerical weather prediction (NWP) and climate modeling. Richardson (1922) [2] took six weeks to produce a six-hour forecast at two points in Central Europe based on the equations proposed by Bjerknes. Charney et al. (1950) [3] used a computer to simulate a barotropic model. Phillips (1956) [4] developed the first successful weather/climate model. Now, the weekly weather forecasting generated from the NWP models are shown in TV and newspapers everywhere. People are concerned about the accuracy of the forecast, especially for severe weather and how long the forecasting is still valuable.
Poincare (1882, 1890, 1892) [5,6,7], discovered that a small change in one state of a deterministic nonlinear system can result in large differences in a later state in his research of the three-body problem, became the first person to discover a chaotic deterministic system which laid the foundations of modern chaos theory. He also introduced phase space (Poincare 1881), Poincare section (map) [7,8], and other concepts to work on dynamic systems. Lorenz (1963) [9] had rediscovered the chaotic behavior and attractors of a nonlinear system when he worked on simplified Benard–Rayleigh convection. Both Poincare and Lorenz emphasized the chaos could come from the uncertainty of initial state. A predictability limit of two weeks (Lorenz 1969) [10] has been suggested and highly cited. Meanwhile, Smagorinsky (1969) [11] documented that “In general even at 21 days, the perturbed map is still far from randomly related to the control—another way of saying that the deterministic limit has not yet been reached”. Arakawa (1966) [12] believed that the predictability limit is not necessarily a fixed number. Pagel et al. (1997) [13] found that regional models are not sensitive to the small variation of initial conditions. Wiin-Nelson (1998) [14] suggested that the practical limit is determined by the uncertainty of the initial state due to the accuracy and distributions of observations and the uncertainty in describing the forcing of a given system. Charney and Shukla (1981) [15], and Hoskins (2012) [16] pointed out that the improvement of the initial data, modeling, parameterization, and better understanding of the large-scale systems (blockings, equatorial waves, etc.), as well as the influence of soil and ocean, which has much longer predictable timescale, increase our predictability significantly. The idea of a predictability limit of two weeks is based on a highly simplified model, which is sensitive to the parameters and model chosen (Hand and Finch 2008) [17]. The dissipation in the atmosphere and oceans also makes the systems differ from the conservative systems discussed by Poincare and Lorenz. The sensitivity of the numerical model was also found in Beltrami (1987) [18] and Wiin-Nielsen (1998) [14]. Recently, Shen (2020) [19], and Shen et al. (2022a, b) [20,21] added more modes in Lorenz models and obtained both point attractors and chaotic attractors due to the aggregated negative feedback of small-scale convective process. Poincare (1892) [7] noted that to study the evolution of a physical system over time, one must model based on a choice of laws of physics and the necessary and sufficient parameters that characterize the system. The equations of the motions in the atmosphere and oceans are much more complicated than those studied by Poincare or Lorenz. It may be impossible to predict the long-term behavior in the chaotic weather/climate systems, but the improvement of the NWP; intensive research and test on prediction systems in many institutes/organizations; and our experiences working on various models indicate that we will be able to predict the change in weather/climate near the saddle points and extend the short-term forecast from a couple weeks to a season or longer.
Even if the dynamic system is stable, there may be numerical instabilities resulting from the method of applying (Greenwood 2003) [22], Sun (2010, 2011) [23,24]. Instead of working on the initial condition, we used the ECMWF reanalysis as the initial and lateral boundary conditions for our real data simulations. We have focused on the equations and numerical modeling, which correspond to the parameters and the mathematics models discussed by Poincare, Lorenz, and Wiin-Nielsen, etc. We also used the observations data to verify the numerical results and improve the numerical model. We are aware of uncertainty and difficulty in CGFD, and the nature is always more complicated than those presented by CGFD. But CGFD is still the best tool to study and predict the motions in the atmosphere and oceans.
This article is limited to the author’s learning and working experiences with the teachers, former students, and colleagues, which only covers a small portion of CGFD. Although majority of the material comes from the previous publications, it may be different from the traditional texts: (Mesinger and Arakawa 1976 [25], Gill 1982 [26], Haltiner and Williams 1983 [27], Pielke Sr. 2002 [28], Lin 2007 [29], Holton and Hakim 2012 [30], Versteeg and Malalasekera, 1995 [31], Vallis 2013 [32], etc.) or the conventional theories. This paper tries to link dynamics instability and weather/climate research. Each section includes the basic equations, approach, and results, although the detailed derivations are incomplete. The details should be referred to the original papers. The reader with numerical analysis, computational fluid dynamics, dynamic meteorology or oceanography should be able to follow. Comments, suggestions, and corrections will be appreciated. We hope readers will discover the beauty, uncertainty, and challenge in CGFD and linkage between nature and sciences.

2. Shallow Water Equations and Numerical Schemes

The shallow water equations (SWE) are simple but capable of capturing the major characteristics of the large-scale motions in the atmosphere and oceans, because of a small ratio of the vertical to the horizontal length scale and in the hydrostatic equilibrium. They have been applied to simulate the Rossby wave, Kevin wave, Poincare wave, geostrophic adjustment, and surface gravity waves, etc. (Pedlosky 1979 [33], Gill 1982 [26], Vallis 2013 [32]), as well as the flow in lakes and rivers. It is noted that the barotropic model is a non-divergence shallow water system. Lorenz (1980) [34] added diffusion in the continuity equation to represent thermal dissipation, and derive a low-order primitive equation by dropping the time derivatives in the divergence equations to study attractor and quasi-geostrophic flow in SWE. His results showed that gravity wave decays with time and the quasi-geostrophic mode remains as expected. It is likely that Kelvin waves may be affected too. The dispersive relation of an inertial surface wave is
σ 2 = f 2 + ( k 2 + l 2 ) g H
where σ is frequency, f is the Coriolis parameter, k and l are the wave number, H is the mean water depth. Another important equation is the conservation of potential vorticity PV:
D D t ς + f h = 0
where ς is the relative vorticity, and h is the depth of water. The corresponding PV in the compressible atmosphere:
D D t ω a · θ ρ = 0
where ω a is the absolute vorticity, θ is potential temperature, and ρ is atmosphere density.

2.1. Instability of Forward–Backward (FB) and Leapfrog (LP) Schemes (Sun 2010) [23]

The 1D linearized shallow water equations are
h t = H u x + δ ν 2 h x 2
u t = g h x + ν 2 u x 2
where H is the mean depth, g is gravity, h and u are depth perturbation and velocity, ν is viscosity. The control parameter δ = 0 or 1. In the staggered C-grids, the finite difference equations in the forward–backward (FB) scheme become:
( h p n + 1 h p n ) Δ t = H u p + 1 / 2 n u p 1 / 2 n Δ x + δ ν h p + 1 n 2 h p n + h p 1 n Δ x 2 ( u q n + 1 u q n ) Δ t = g h q + 1 / 2 n + 1 h q 1 / 2 n + 1 Δ x + ν u q + 1 n 2 u q n + u q 1 n Δ x 2
where q = p ± 1/2, p is the grid for h, and Δx is the spatial interval. A wave-type (normal mode) solutions at the nth time step are:
h p n = h ^ n exp [ i k p Δ x ] = h ^ 0 λ n exp [ i k p Δ x ] u q n = u ^ n exp [ i k q Δ x ] = u ^ 0 λ n exp [ i k q Δ x ]
where i = 1 , h ^ n and u ^ n are the amplitude of h and u at the nth time step; k is the wave number. We obtain the difference equations
h ^ n + 1 u ^ n + 1 = 1 δ S 2 i Δ t H X i g Δ t X ( 1 δ S 2 ) 1 S 2 R 2 h ^ n u ^ n
and the eigenvalue
λ = 2 ( 1 + δ ) S 2 R 2 ± [ 2 ( 1 + δ ) S 2 R 2 ] 2 4 ( 1 S 2 ) ( 1 δ S 2 ) / 2
where X = sin [ k Δ x / 2 ] Δ x / 2 ,   R 2 = g H Δ t 2 X 2 = ( C X Δ t ) 2 = 4 C o 2 sin 2 ( k Δ x / 2 ) , S 2 = ν Δ t X 2 , C = g H , and Courant number C o = C Δ x Δ t .
For simplicity, we set g = 1, H = 1, and Δx = 1.
The difference equations for the Leapfrog (LF) scheme are:
( h p n + 1 h p n 1 ) 2 Δ t = H u p + 1 / 2 n u p 1 / 2 n Δ x + δ ν h p + 1 n 1 2 h p n 1 + h p 1 n 1 Δ x 2 ( u q n + 1 u q n 1 ) 2 Δ t = g h q + 1 / 2 n h q 1 / 2 n Δ x + ν u p + 1 n 1 2 u p n 1 + u p 1 n 1 Δ x 2
and the corresponding eigenvalue
λ 2 = 1 ( 1 + δ ) S 2 2 R 2 ± 1 ( 1 + δ ) S 2 2 R 2 2 ( 1 2 δ S 2 ) ( 1 2 S 2 )
Without viscosity (i.e., ν = 0), the FB is neutrally stable, |λ| = 1, when the Courant number C o = C Δ t Δ x < 1 ; and the LF is also neutrally stable when Co < 0.5 (Mesinger and Arakawa 1976 [25], Haltiner and Williams, 1980 [27]). At Co = 1, the FB has a repeated eigenvalue. λ 1 , 2 = 2 R 2 ± [ 2 R 2 ] 2 4 / 2 = 1 for 2Δx waves. Equation (5) becomes
x n + 1 = h ^ n + 1 u ^ n + 1 = 1 2 i 2 i 3 h ^ n u ^ n = A x n
An eigenvector corresponding to λ1 = −1 is x 1 = 1 i . The generated eigenvector x 2 can be found from
A x 2 = λ 1 x 2 + x 1
and
x 2 = 0.5 0
If we define a matrix P = ( x 1 , x 2 ) = 1 0.5 i 0 , then, we can have
A B = P 1 A P = 0 i 2 2 i 1 2 i 2 i 3 1 0.5 i 0 = 1 1 0 1 = λ 1 0 λ
Which is in Jordan block form (Nobel 1969) [35], and since
λ 1 0 λ n = λ n n λ n 1 0 λ n
We obtain
A n = P A B n P 1 = 1 0.5 i 0 λ n n λ n 1 0 λ n 0 i 2 2 i
= 1 0.5 i 0 ( 1 ) n n ( 1 ) n 1 0 ( 1 ) n 0 i 2 2 i
Red dashed line in Figure 1 shows the magnitude of 2Δx wave linearly increases with time and with 2Δt oscillation, which can be controlled by the fourth-order Shuman smoothing (blue solid line) (Sun 2010) [23]. When Co = 0.5, (i.e., R = 1), and S 2 = 0 , (i.e., ν = 0), the LF scheme has the repeated eigenvalues too, λ 1 = λ 2 = i and λ 3 = λ 4 = i . Hence, it becomes weakly unstable, the same as the FB scheme.
Sun (1982) [36] applied viscosity to the burger equation to suppress instability. But in SWE, viscous terms create an unstable 2Δx and 3Δx waves at Co = 1, and an unstable 2Δx wave at Co = 0.8 in the FB. The growth rate becomes less than one for longer waves and the maximum damping occurs at 4Δx or 3Δx waves with δ = 1 (Figure 2). The LF also becomes unstable at R2 = 1, (i.e., 2Δx-wave at Co = 0.5). The amplification factor for 2Δx-wave equals 2 for ν = 0.25 (i.e., S2 = 0.5), as shown in Figure 3, because viscosity enhances the 2-Δt oscillation for the short waves (wavelength ≤ 3Δx). Instability also shows up for 2Δx-wave when ν > 0.13 at Co = 0.4 in LF. Viscosity has the largest damping at 4Δx waves when Co = 1 in FB and Co = 0.5 in LF; and at 3Δx waves when Co = 0.8 in FB and Co = 0.4 in LF. The instability becomes weaker with δ = 0.

2.2. Modified Leapfrog Scheme for Shallow Water Equations

Sun and Sun (2011) [37] proposed to replace h q + 1 / 2 n in (7) by h ¯ q + 1 / 2 n = h q + 1 / 2 n + 1 + 2 h q + 1 / 2 n + h q + 1 / 2 n 1 4 .
The finite different equations become
( h p n + 1 h p n 1 ) 2 Δ t = H u p + 1 / 2 n u p 1 / 2 n Δ x
( u q n + 1 u q n 1 ) 2 Δ t = g h ¯ q + 1 / 2 n h ¯ q 1 / 2 n Δ x
After solving (15), we obtain h q + 1 / 2 n + 1 in (16). No iteration is required. It is a second-order accuracy in both space and time. The eigenvalues are:
λ = 2 R 2 ± [ 2 R 2 ] 2 4 / 2
|λ| = 1, if Co < 1 which is twice allowed in the original LF scheme. The eigenvalue of (17) is identical to (6) of the FB scheme without viscosity (S = 0) (Sun 1984a, 2010) [23,38]. This simple method enables the LF users to save one half of the computing time.

2.3. New Scheme for Shallow Water Equations, Sun (2011) [24]

Typically, finite difference schemes produce oscillations behind the shock in simulating the dam break. Therefore, more sophisticated methods have been proposed (Roe 1981 [39], Toro 1999 [40]), including characteristic based method (Wang and Yeh 2005 [41], Oh 2007 [42], etc.) and various Riemann solvers (MacCormack 1969 [43], etc.). Unfortunately, it is impossible to diagonalize the matrices on both x- and y-directions of characteristic equations simultaneously, those methods require complicated algorithms with some uncertainty (Oh 2007) [42]. Sun (2011) [24] proposed a new method by using the time average on the right-hand sides of (3):
( h p n + 1 h p n ) Δ t = H u p + 1 / 2 n + 1 / 2 u p 1 / 2 n + 1 / 2 Δ x = H u + Δ t 2 u t x p n = H u Δ t 2 g h x x p n = H u x p n + g H 2 2 h x 2 p n
h p n + 1 = h p n + H Δ t Δ x ( u p 1 / 2 n u p + 1 / 2 n ) + g H Δ t 2 2 Δ x 2 ( h p 1 n 2 h p n + h p + 1 n )
( u q n + 1 u q n ) Δ t = g h q + 1 / 2 n + 1 / 2 h q 1 / 2 n + 1 / 2 Δ x = g 2 h n + h n + 1 x q
u q n + 1 = u q n + Δ t g 2 Δ x h q 1 / 2 n h q + 1 / 2 n + h q 1 / 2 n + 1 h q + 1 / 2 n + 1
The last term in (19) is like the diffusion term in the Lax–Wendroff scheme (Mesinger and Arakawa 1976) [25]. Using the normal mode solutions of (4), we obtain the eigenvalues
λ = 1 2 C o 2 sin 2 ( k Δ x / 2 ) ± 2 C o sin ( k Δ x / 2 ) C o 2 sin 2 ( k Δ x / 2 ) 1
where C o = g H Δ t / Δ x is the Courant number. (22) is identical to (17) or (6) with S = 0. Equation (16) shows that|λ| = 1 if Co < 1. Same as the FB discussed previously (Sun 2010) [23]. The phase angle ϕ and the phase speed Ĉ of the numerical scheme are
λ = λ exp ( i ϕ ) = λ exp ( i k C Δ t )
and
C p h = φ / ( k Δ t ) = sin 1 4 2 C o 2 sin 2 ( k Δ x / 2 ) 2 2 / ( k Δ t )
The phase speed Cph is the same as the FB (Sun 1984a; 2010) [23,38]. It is noted that the diffusion term in (19) does not cause damping in the solutions because the amplification factor |λ| = 1 when Co < 1. On the other hand, the diffusion term in the Lax–Wendroff scheme can severely damp the short waves (Mesinger and Arakawa, 1976) [25]. The ratio of the computed phase speed to the analytical phase speed C p h / C is presented in Table 1.

2.4. Nonlinear Shallow Water Simulations

The flux form of the nonlinear SWE can be written as:
h t = ψ x ϕ y
ψ t + u ψ x + v ψ y f ϕ = g h ( h + h s ) x k ψ
φ t + u φ x + v φ y + f ψ = g h ( h + h s ) y k φ
where ψ = hu and ϕ = hv, and hs is bottom elevation, k is friction

2.4.1. Simulation of Dam Break: Water Depth Is 10 m within a Radius of 11 m, and 1 m Outside

It is also assumed f = hs = k = 0. Following the procedure of Section 2.3, we can write the difference–differential equations of (25)–(27):
h p , q n + 1 h p , q n Δ t = ψ n + 1 / 2 x p , q ϕ n + 1 / 2 y p , q = ψ + Δ t 2 ψ t x p , q n ϕ + Δ t 2 ϕ t y p , q n = ψ n x p , q ϕ n y p , q + Δ t 2 u ψ x + v ψ y + g 2 h 2 x x p , q n + Δ t 2 u ϕ x + v φ y + g 2 h 2 y y p , q n
ψ p , q n + 1 ψ p , q n Δ t = u ψ n + 1 / 2 x p , q v ψ n + 1 / 2 y p , q g 2 h n + 1 / 2 2 x p , q   = u ψ Δ t 2 u u ψ x + v ψ y + g 2 h 2 x + ψ u u x + v u y + g h x + g 2 h n + h n + 1 2 2 x p , q n v ψ Δ t 2 v u ψ x + v ψ y + g 2 h 2 x + ψ u v x + v v y + g h y y p , q n
and
ϕ p , q * n + 1 ϕ p , q * n Δ t = u ϕ n + 1 / 2 x p , q * v ϕ n + 1 / 2 y p , q * g 2 h n + 1 / 2 2 y p , q * = u ϕ Δ t 2 u u ϕ x + v ϕ y + g 2 h 2 x + ϕ u u x + v u y + g h x x p , q * n v ϕ Δ t 2 v u ϕ x + v ϕ y + g 2 h 2 x + ϕ u v x + v v y + g h y n + g 2 h n + h n + 1 2 2 y p , q *  
where p* = p ± 1/2, and q* = q ± 1/2.
This scheme is simple, accurate, and conserves the mass and other properties in an inviscid fluid. Most important is that it can be applied to the nonlinear SWE without any diffusion or numerical smoothing. Hence, it can reproduce the solitary Rossby wave, hydraulic jump, the flow moving over mountains, and vortex-merging, etc. Without smoothing/diffusion, our numerical results are usually different from the previous studies that were generated by the models with smoothing. Hence, the physical interpretations can be different. Solutions are usually diverse among elliptic, hyperbolic, and parabolic partial differential equations.
The domain consists of 400 × 400 grids with Δx = Δy = 1 m. The fourth-order Shuman smoothing with coefficient α = 0, 0.1, or 0.2 is applied to the prognostic variables each time step on both second-order and fourth-order schemes (Sun 2011) [24]. Most previous publications show the results at t = 0.69 s. The simulated vertical cross section at t = 0.69 s for both the second-order (α = 0, and 0.2) and fourth-order schemes (α = 0.2) are shown in Figure 4. The differences among the second-, third- and fourth-order schemes are insignificant. The results show that the water level h remains at 10 m near the center; the shock reaches at r = 18 m from the center with a crest of h ≈ 3.7 m at r = 16 m; and a band of low water (h ≈ 3.2 m) forms around r = 13 m. They are consistent with the second-order Godunov-type numerical model (Fujihara and Borthwick 2000) [44]; the unstructured triangular meshes and Roe’s flux function solver (Anastasiou and Chan 1997) [45]; and the characteristic Lagrangian method (Oh 2007) [42], etc.
As pointed out by Tseng and Chu (2000) [46], the traditional MacCormack scheme breaks down at t = 0.22 s for this case. The current model running to 4 s without smoothing indicates that the diffusion term in (19) may stabilize the system. The water moves away and creates a deep hole near the center by t = 3 s. The hole deepens and widens while the shock propagates away. After t = 22 s, water starts returning to the center, eventually a high-water column with a sharp gradient forms near the center, while the original shock continuously propagates away, as shown at t = 25 s in Figure 5. The simulations agree excellently with Toro (1999) [40].
The total mass is conserved. The total energy slightly decreases with time due to the truncation errors and smoothing applied. The total energy is given by
E = p q h p , q u p , q 2 + v p , q 2 2 + g ( h p , q ) 2 2 Δ x Δ y
The ratio of the total energy to its initial value Ra(t) = E(t)/E(t = 0) is about 0.96 at t = 1 s and 0.91 at t = 3 s. Using the Lagrangian method, Oh (2007) [42] obtains Ra = 0.85 at t = 1.5 s. The current scheme manages both mass and kinetic energy well, which is a major challenge to dam-break simulations (Liang et al., 2004) [47].

2.4.2. Solitary Rossby Wave (Sun and Sun 2013) [48]

The flux form upstream scheme, the center difference scheme, and the Lax–Wendroff scheme, etc., failed to simulate the Rossby soliton (Chu and Fan, 2010) [49]. The major difference between the Lax–Wendroff scheme and Sun (2011) [24] is that the diffusion term in the Lax–Wendroff scheme can severely damp the short waves, but no damping in Sun scheme as discussed in Section 2.3. Boyd (1980, 1985) [50,51] has shown that an asymptotic solution exists to the inviscid, nonlinear shallow water equations for the Rossby soliton on an equatorial beta plane, which are governed by Korteweg–de Vries (KDV) equation (1895) [52].
u t 6 u u x + 3 u x 3 = 0 :   solution :   u ( x , t ) = c 2 sech 2 c 2 ( x c t )
The soliton propagates to the west at a fixed phase speed of linear wave c, without change in shape. After a coordinate transfer s = x c t . The equations become
u 0 ( s , y , t ) = ζ ( s , t ) 9 + 6 y 2 4 e 0.5 y 2
v 0 ( s , y , t ) = 2 y ζ ( s , t ) s e 0.5 y 2
h 0 ( s , y , t ) = ζ ( s , t ) 3 + 6 y 2 4 e 0.5 y 2
and the variable ϛ(s, t) satisfies
ς t 1.5366 ς ς s 0.098765 3 ς s 3 = 0
Which is the Korteweg–de Vries (KDV) equation with the exact solution
ζ ( s , t ) = 0.772 B 2 sech 2 ( B ( s + 0.395 B 2 t ) )
Equations (32)–(36) were also referred as the zeroth-order solution (Boyd 1980, 1985) [50,51]. The initial condition is given by (32)–(36) with t = 0.
The model equations including the Coriolis parameter f and the terrain height hs are shown in (25)–(27) without friction (k = 0). Using Sun (2011) [24] scheme, Sun and Sun (2013) [48] simulated solitary waves which are comparable with solutions of Boyd (1980, 1985) [50,51]. The dispersive waves from the linearized equations are also comparable with Boyd’s results as shown in Figure 6 and Figure 7. The amplitude difference in soliton may come from the initial condition not being resolved by the grids. The phase lag is due to model phase speed is slightly slower than analytic solution, the lag becomes larger for the short waves in Figure 7.
The results also confirm that solitary waves of large amplitude containing a recirculating fluid trapped within the moving disturbance (Figure 8). However, the simulated recirculation for B = 0.6 is more complicated than Boyd’s, which may indicate the pattern may not be well represented by Boyd’s asymptotic series. The simulations also show that terrain can slow down the speed of the soliton but do not significantly change the general pattern of the waves farther downstream (Sun and Sun 2013) [48]. Solitary wave solutions are also related to the non-dissipative equations of the Lorenz 1963 model [9]:
d X / d t = σ ( X + Y ) ,   d Y / d t = X Z + r X Y ,   and   d Z / d t = X Y b Z  
with   σ = 10 ,   r = 28 ,   and   b = 8 / 3 .
Combining those equations into a single equation of Z, Shen (2020) [19] obtained:
d 2 Z d ς 2 + 3 Z 2 4 r Z = 0 , which is identical to the integral of Korteweg–de Vries equation (i.e., d 2 f d ξ 2 + 3 f 2 c f = 0 ). It will be interested to see the modification in Figure 1 of Shen (2020) [19] if he keeps more terms in his model.

2.4.3. Hydraulic Jumps

The inviscid, one-dimensional SWE without f is
u t + x u 2 2 + g ( h + h s ) = 0
And B e = u 2 2 + g h s + h is the Bernoulli function, and h f = h s + h . The jump is defined where h and u are discontinuous. The flux from becomes
h u t + x h u 2 + g h 2 2 + g h h s x = 0
Let us define the potential of mass flux J = J 1 + J 2 where
J 1 = h u 2 + g h 2 2   and   J 2 = g h h s x
The equation of total energy density per unit length E becomes:
t E = t h u 2 2 + g h h s + h 2 = x u h u 2 2 + u h g h s + h = x u h · B e
It is assumed that mass flux remains constant; ∂h/∂x = 0 at the inflow; and open boundary condition at the outflow, they are more realistic than the periodic boundary conditions used by Houghton and Kasahara (1968) [53] (will be referred as HK) and others.
Sun (2018) [54] has presented six different cases, only Case B in (Sun 2018) [54] will be discussed here. The terrain is:
h s = h c 1 x x c x L 2   for   0 x x c x L   and   h s = 0   otherwise .
where h0 = 20 cm, hc = 10 cm, Δx = 1.0 cm, xL= 40 Δx, xc = 600 or 400 cm (600 or 400 grid), g = 980 cm s−2, and the domain L = 1200 Δx, and u0 = 42 cm s−1. Combining the conservation of mass flux and the steady Bernoulli function, HK obtained:
F o 2 2 U 3 + U ( M F o 2 2 1 ) + 1 = 0
where U = u u 0 , M = h s / h 0 , and the Froude number F 0 = u 0 c 0 = u 0 g h 0 , where u 0 and h0 are velocity and depth at inflow. There exist three real roots of (41) when
M M * F o 2 2 3 2 F o 2 / 3 + 1
Figure 9a,b show the time sequence of the simulated hf and u at t = 0 (black), 2.4 (red), 4.8 (green), 7.2 (dark blue), and 9.6 (purple). At beginning, the shock waves propagate both directions. The simulated Φc = hc/h0 and F c = u c / c 0 at crest at =4.8 s agree with asymptotic solution of HK. There are no real solutions of Equation (41), because M * = F 0 2 2 3 2 F 0 2 / 3 + 1 = 0.3728 .
Therefore, the flow with singularity originally has changed to a critical flow by decreasing velocity and increasing flow height at the inflow. At the crest, the simulated velocity U c s is identical to U * s , which has a physical solution (i.e., a real number instead of a complex number) of Equation (41). Equation (37) is also held at the crest, but not at the jump at x = 627 with u =175.7 cm s−1, where the velocity calculated from Equation (41) is only 157.58 cm s−1. The Bernoulli function (black line) drops from 22,948 to 20,485 cm2 s−2 across the jump at t = 96 s, but J remains almost the same (284,582 vs. 284,581 cm3 s−2), except small kinks at the jump and obstacle foots, because it is derived from an inviscid Equation (38) of mass flux in Figure 9.
The shockwave reflected from the inflow boundary changes the flow property from the region II of HK initially to the critical flow at the final steady state, which occurs to all cases having hydraulic jumps. It may imply that the critical flow is a stable state according to Poincare or an attractor to Lorenz.

2.4.4. Vortex Moving over Central Mountains Range in Taiwan (Sun 2016) [55]

Previously, the surface friction was believed secondary and not included in Lin et al. (2005) [56], Huang and Lin (2008) [57]. The surface friction in momentum equations in (26) and (27) is assumed proportional to terrain elevation:
k = α h s
The potential vorticity in the SWE becomes
d d t ln ζ + f h = 1 ζ + f ( k u ) y ( k v ) x
Equation (44) is used to simulate the flow in the troposphere (with a depth ~10 km), which is topped by an undisturbed stratosphere (Matsuno 1966) [58]. In water, we obtain: the reduced gravity g′ = 4 m s−2; depth H = 100 m; the horizontal length scale R = 0.2 R*, where R* is the radius of the vortex in the real atmosphere; the Coriolis parameter f = f23 (ϕ = 23° N) = 5.7 × 10−5 s−1; and the peak of mountain in the model hm = 10 m at the grid (100, 85) corresponds to the mountain peak of h m * = 3658 m with a 4 km horizontal resolution. Table 2a,b show the radius of vortex (R); fluid depth (H), the height of mountain peak (hm), Rossby number (Ro), Froude number of the mountain ( F r = U g h m ) (Sun and Chern 1994 [59], Sun and Sun 2015 [60]), Froude number far away from topography ( F r = U g H ), and LR/R, where LR the radius of deformation. They are comparable to the corresponding parameters in the atmosphere (Sun 2016) [55] indicated by superscript *. The initial surface height depression of the vortex is
h = h o / 1 + x x o R 2 + y y o R 2 = h o / 1 + x i o Δ x R 2 + y j o Δ y R 2
and geostrophic wind relation is used to represent the wind around the vortex initially.

2.4.5. Inviscid Vortex with α = 0 and U = −1.2 m s−1 (U* = −6 m s−1)

Figure 10a shows the trajectories of vorticities A1, A2, and A3, initially located at grids LN(280, 170), LC(280, 150), and LS(280,130), respectively. The amplitudes of initial surface depression ( h o ) are 2, 3, and 4 m; the initial radii R are 16, 24, and 40 km (R* = 80, 120, and 200 km). The movement of vortices is independent of size or intensity. The streamlines, and wind speed (in red) in Figure 10b show the maximum wind of 27 m s−1 at r = 60 km in the atmosphere for A2 when the vortex is located over the Central Mountain Range (CMR) at 96,000 s (26.7 h). Figure 10c–e show
D Π D t D ( PV ) D t = d d t ζ + f h = t ζ + f h + v · ζ + f h = 0 .
It confirms that the PV is conserved when an inviscid vortex passes over the mountains.

2.4.6. Trajectories of PV with Surface Friction

The trajectories of the vortex with h o = 3 m and R = 24 km (R* = 120 km in the real atmosphere) are shown in Figure 11; in which the solid lines represent the trajectories with U = −1.2 m s−1 (U* = −6 m s−1); dash-dot-lines with U = −0.8 m s−1 (U* = −4 m s−1); the line with * for U = −1.0 m s−1 (U* = −5 m s−1); the line with Δ for U = −1.6 m s−1 (U* = −8 m s−1), except the black solid line for F4 (with U = −1.6 m s−1 (U* = −8 m s−1), α = 1.5 × 1.0−4 m−1/2s−1, starting at LN); the lines with o including the beta (β = df/dy) effect for B4 and B5; the yellow line with black dots for α = 0 (no friction) for A2 and B4; magenta for α = 1 × 10−5 m−1/2 s−1; cyan for α = 3 × 10−5 m−1/2s−1; green for α = 5 × 10−5 m−1/2s−1; blue for α = 8 × 10−5 m−1/2s−1; red for α = 1 × 10−4 m−1/2s−1; and black for α = 1.5 × 1.0−4 m−1/2s−1.
Since there is no surface friction over the ocean, at beginning vortices move westward. The vortices with beta effect ( β = d f / d y ), deflect slightly northward initially, because a higher value of f is carried southward by cyclonic circulation ahead of the vortex.
With moderate and large surface friction, when a vortex approaching Taiwan, surface friction creates a strong northerly channel flow westside of the vortex. A strong southerly flow also forms eastside of the vortex because the slowdown vortex is pushed by the easterly flow from behind. Meanwhile, the vortex deforms. The asymmetric vortex and the accompanying two high-wind zones rotate cyclonically while the vortex moves westward. At first, the northerly channel flow is stronger than the southerly flow on the east of the vortex, causing the vortex to deflect southward. The zone of high PV and strong wind continues rotating around the vortex center. The movement of the center depends on the summation of the velocity vector on the vortex. It can continue moving westward or deflect southward or northward. It can also retreat and form a loop. The mechanism is different from the tradition theory of channel flow or the inner core asymmetric flow coming from mid-troposphere (Wu et al., 2015) [61]. If a large vortex moves around the northern coast, it can induce a strong southwesterly flow to create flooding on the Western Taiwan. On the other hand, the secondary circulation is weak when the vortex directly passed over the CMR. Although SWE may explain the effects of mountains and surface friction, the dynamics of a real typhoon are far more complicated than the SWE.

2.4.7. Vortices Merge (Sun and Oh 2022) [62]

The domain consists of 450 × 450 grids, with Δx = Δy = 0.025 m. Following Waugh (1992) [63] and Oh (2007) [42], we set g = 1 m s−2, and the initial surface perturbation is
h = h o / x x o r o 2 + y y o r o 2 + 1 1.5
where yo = 0, xo = ±0.5 × do (±40 Δx), ro. = 1.2 m (48 Δx), and geostrophic wind at t = 0.
Case A (H = 1 m, f = 1 ms−2, do = 2 m and ho = −0.2 m at t = 0)
If we use the length scale L = do/2 and the maximum velocity U, we obtain the Rossby number, Ro = U/fL = 0.16, and radius of deformation R a = g H / f = 1 m at t = 0. Figure 12a shows the initial PV (thick black), surface heigh h (thin white), streamlines (thin black), and speed (shaded color). Figure 12b reveals that the strong wind zones (red zone) advect downstream at t = 20 s. The PV rotates slower at point B, but moves faster at point A; hence, development of the filaments is due to differential velocity across the PV contours. The streamlines show that the vorticity of filaments does not advect or ejected from the vortex cores as proposed by Kimura and Herring (2001) [64], Cerretelli and Williamson (2003) [65], and Melander et al., 1987 [66], etc. Our simulations show that the filaments are shattered with little effect to vortex merging in the inner core.
The asymmetric structure of quadrants does exist in the simulated PV. Advection generates positive ∂∏/∂t on the lee of the PV ridge, where contours become more convex and the angle between PV contours and streamlines increases. Hence, positive PV can be effectively transported along two parallel but oppositive paths, from A to B or from C to D shown in Figure 12c. Hence, the vortices may move closer while rotating around each other. The vortices never merge when negative ∂∏/∂t prevails between two vortices. Hence, the advance of positive ∂∏/∂t provides a simple yet eloquent mechanism for vortex merging, which is superior to the following mechanisms: the axisymmetrization principle (Melander et al., 1987 [66], Melander et al., 1988 [67]); the angle between vortex major axis and vortex center line (Huang 2005) [68]; interactions between vorticity and rate of strain (Moore and Saffman, 1975 [69]; Brandt and Nomura 2007 [70]); or production of palinstrophy in Kimura and Herring (2001) [64] and Cerretelli and Williamson (2003) [65], etc.
The contours of h and PV at t = 60 s in Figure 12 are in good agreement with the Lagrangian simulations (Oh 2007) [42]. The maximum speed, 0.18 m s−1, occurs along the inner core boundary. The velocity at the vortex centers, 0.09 m s−1, is identical to the angular velocity, ω = v/r = 0.09 m s−1 /0.36 m = 0.25 = 14° s−1. They move along the constant h-contour and rotate around the center of system, like the observed cold front moving along the eastern coastal mountain ranges of Taiwan and China with the observed wind U instead of a Kelvin wave of U + g h (Sun and Chern 2006) [71]. As time increases, filaments stretch longer and wrap around the shrinking inner core at 100 s (Figure 12e). The shape of inner core gradually changes from elliptic to circular. The filaments become long, thin strips around the center (Figure 12f at t = 120 s). The vorticities are still distinguishable at the end of integration, t = 130 s and total kinetic energy, mass, and PV remain as initial values (Sun and Oh 2022) [62].
Two-dimensional, inviscid turbulence has both kinetic energy and vorticity square (enstrophy) in the inertial range satisfy, E ( k ) ε 2 / 3 k 5 / 3 and Ζ ( k ) η 2 / 3 k 3 , where ε is the rate of cascade of kinetic energy per unit mass, η is the rate of cascade of mean-square vorticity. The backward energy cascade from higher to lower wavenumbers k. The −3 range gives an upward vorticity flow (Fjortoft 1953 [72], Kraichnan [73]). They are different from the energy cascade from lower wave number to higher wave number in 3D flow (Kolmogorov 1941 [74], Obukhov 1941 [75]). Vortex tubes/filaments can effectively transfer the energy from the resolvable scale to small and/or subgrid scales. Meanwhile, the unstable disturbance can grow from a small size to large one. The weather/climate models still require rigorous turbulence parameterizations to handle the cascade of energy and enstrophy.

3. Linearized Equations

In fluid dynamics, the dispersive relationship of waves is derived from the linearized equations, including acoustic waves, gravity waves, inertia internal gravity waves, Poincare waves, Kelvin waves, and Rossby waves, etc. Stability analysis of Benard–Rayleigh convection, Eady baroclinic wave (Eady 1949 [76], Pedlosky 1979 [33], Gill 1982 [26], etc.), and numerical schemes are also derived from the linearized equations, as shown in Section 2.1 and Section 2.3.

3.1. Coexist of Two Different Types of Cloud Bands (Sun 1978) [77]

Figure 13a shows that both longitudinal and transverse cloud bands coexist to the east of disturbed wave trough near the Equator (Malkus and Riehl 1964) [78]. The rows of longitudinal mode composed of medium or small clouds are parallel to the wind shear with a distance between two rows about 25 km. The large clouds line up at right angles to the wind shear (i.e., transverse mode) with a wavelength of 100 km. Sun (1978) [77] applied the normal mode ( ϕ = Φ ( z ) exp ( i ( k x + l y σ t ) ) ) for perturbation in the linearized equations including a wave-CISK type of latent heat to study those modes. His results show the coexistence of the longitudinal mode consisting of the small–mid size clouds in the lower atmosphere (lest or near 5–6 km) where the buoyancy is greater than effect of stratification, and the transverse model with the cloud tops reaching 9.6 km, in which the total cooling from stratification is slightly larger than buoyancy. Hence, the transverse mode chooses the strongest wind shear direction (Figure 13b). Kuo (1963) [79] pointed out that in a stable stratified fluid, the transverse mode is more unstable than the longitudinal mode, but the reverse is true for an unstable stratified flow.

3.2. Solve Linearized Equation as Initial Value Problems

The normal mode ϕ = Φ ( z ) exp ( i ( k x + l y σ t ) ) or ϕ = Φ exp ( i ( k x + l y + m z σ t ) ) cannot be applied to latent heating (positive warming inside cloud and negative cooling outside) or time dependent basic state. Hence, we treated the linearized equations as an initial value problem.

Rainbands and Symmetric Instability (Sun, W. Y., 1984b) [80]

We set the y-axis to the direction of wind shear coinciding with the direction of cloud band when the buoyancy is dominated, as discuss in Section 3.1. Because the variation along y-direction is much smaller than that along x-direction, the linearized equations become:
u t f v + 1 ρ 0   p   x = 0 ,
v t + u V ˜ x + w V ˜ z + f u = 0
w t + 1 ρ 0 p z + g θ θ 0 = 0
θ t + u θ ˙ ˜ x + w θ ˙ ˜ z = θ ˙ ˜ Q c p T ˜ ( heating )
u x + w z = 0
The environment is assumed in thermal wind balance.
g θ 0 θ ˜ x = f V ˜ z
The equations can be combined into a single vorticity equation:
2 ς t 2 = 2 t 2 w x u z = 2 t 2 2 ψ z 2 + 2 ψ x 2 = 2 S 2 2 ψ x z N 2 2 ψ z 2 F 2 2 ψ z 2 + g θ o θ ˙ ˜ c p T ˜ Q x
where
N 2 = g θ o θ ˜ z ,   S 2 = g θ o θ ˜ x = f V ˜ z ,   F 2 = f ( f + V ˜ x ) ,   w = ψ x ,   u = ψ z
If Q = 0, and assume a normal mode solution:
ψ = Ψ exp ( i ( k x + m z σ t ) )
We obtain
σ 2 = N 2 k 2 + F 2 m 2 2 S 2 k m k 2 + m 2
Equation (55) is identical to inertia internal gravity waves if V ˜ x = V ˜ z = 0 . Fjortoft (1950) [81], Ooyama (1966) [82], Hoskins (1978) [83] and others found symmetric instability in (55) if
R i = g θ o θ ˜ z V ˜ z 2 < f 2 F 2 = f 2 f ( f + V ˜ x )
Including latent heat by assuming Q > 0 in ascending and Q = 0 in descending motion, Sun (1984b) [80] showed the stream function at t = 2.2743 × 104 s in Figure 14a, and heating rate in Figure 14b. The temperature in the entire cloud layer increases: area inside clouds due to condensation, and cloud-free area by subsidence warming. Cooling occurs right above the cloud layer due to penetrative convection (Sun 1975 [84], Kuo and Sun 1976 [85]). The cooling in the planetary boundary layer (PBL) is not exactly beneath the condensation warming; but slightly toward the cool side because of cold advection coming from the left. It generates a local circulation and induces the rolls to move toward the warm side, as shown in Figure 14c at t = 2.653 × 104 s.
Meanwhile, the small bands are suppressed by the subsidence of the larger bands, as observed squall lines in Midwest of the USA on 3 April 1974 (Figure 15). The downward motion also carries the high wind from the upper layer down to PBL, which enhances the low-level convergence and may cause the wide spread of tornados (Agee et al., 1975 [86]). The descending flow carrying warm, dry air may play a significant role in the development of severe storms too.

3.3. Cloud Bands over Tropical Continent

Sun and Orlanski (1981a) [87] solved the linearized equations and found the coexisting 1-day and 2-days mesoscale waves (Figure 16a,b) excited by resonance of trapezoidal instability and the land sea breeze. The energy of those inertia internal gravity waves propagates upward (Figure 16b) as simulated by Kuo and Sun (1976) [85]. They spread inland too, but there was no convection over the ocean (Figure 17) because the variation of the sea surface temperature is small. Those mesoscale convections were also confirmed by the nonlinear simulations (Sun and Orlanski, 1981b) [88] and the observed cloudiness in tropical regions (Orlanski and Polinsky 1977) [89]. Including the beta effect, Sun (1995b) [90] showed that the growth rate of instability is unsymmetric with respect to latitude.

3.4. Ekman Layers in Atmosphere and Ocean (Sun and Sun 2020) [91]

The equations of 1D Ekman layer are:
u ρ t f v ρ = f v g ρ + z μ u ρ z v ρ t + f u ρ = f u g ρ + z μ v ρ z
They   can   be   combined   into   a   sin gle   equation : W t + i f W = i f W g + z μ W z
where W = u ρ + i v ρ , and W g = u g ρ + i v g ρ ( u g   and   v g ) are geostrophic wind in x and y direction, respectively, and i = 1 . It is assumed ρ = 1000 kg m−3 in the ocean, and ρ = 1 kg m−3 in the atmosphere. The domain is from −400 m to 4 km with a spatial interval of 2 m in the ocean and 4 m in the atmosphere, respectively. The stress ρ μ W z is calculated at the mid-point. The velocity (u, v) is continuous at the interface, z = 0, but the surface stresses can be different across it. Viscosity μ is function of both space and time in both atmosphere and ocean (Figure 18), which consists of 1 day, half-day, and other modes.

3.4.1. Simulation at 30° N with Diurnal Variation in Atmosphere and Ocean (with δ z o _ max = 40   m and μ o _ max = 0.04 m2 s−1) (Case C of Sun and SUN 2020)

The time sequences at z = 0 (blue, multiplied by 10) and 600 m (red) in Figure 19a show that the parcels rotate clockwise, and the amplitudes grow with time. The blue Δ shows a weak wind at 12L12 (12 p.m. day 12), and blue ο reveals a strong wind at 24L12 (midnight), corresponding to a nocturnal low-level jet (LLJ) (Blackadar 1957 [92], Palmer and Newton 1969 [93], etc.) in the atmosphere. Meanwhile a jet develops in the upper oceanic mixed layer during daytime and week current at night (Price et al., 1986 [94]). The growth rate of velocity is about 1.00043/h at z = 600 m in Figure 19b.
Equation (58) is identical to the equation of parametric resonances (Kalashnik 2018 [95]). The velocities at 03L12 (dotted dashed green),12L12 (thick dashed red), 15L12 (dotted dashed magenta), 24L12 (thick dashed black), and hourly velocity oscillations at various heights in the atmosphere (Figure 19c) show that the individual spiral significantly departs from the daily average (thick blue), and the shape is close to the analytic solution (thick green). The velocity grows from a small elliptical shape near surface to larger and more circular circulations at 200 m (blue), and 600 m (red). Large amplitude oscillations also exist at the upper mixed layer and beyond: at 1 km (green); 1.6 km (blue), because of resonance with inertia oscillation of the Coriolis force. Velocities with a large magnitude also develop in the oceanic mixed layer and beyond. Ekman spirals rotate a complete cycle in 24 h (Figure 19d). Figure 19e shows that in the morning, the velocity at surface is parallel to the stress, and the magnitude reaches ~0.1 m s−1 as observed at 30° N (Price et al., 1986 [94]). The average velocity (thick blue) at −4 m is ~60° to the right of the averaged wind stress, the averaged surface stress is 0.1 Pa and the angle ψ between the velocity at −2 m and surface stress, is 58°. They agree with 0.08 Pa and ψ = 60°, respectively, observed at 30.9° N, 123.5° W (Price et al., 1986 [94]). The simulations also show a compacted spiral at night (Figure 19d,e). The daily averaged spiral of D12 (thick blue in Figure 19e) is between 12L and 24L, but quite different from either of them. It is close to the analytic solution (thick green). The averaged spiral is also comparable with the observations shown in Figure 19f when the mean stress is aligned to the direction of the observed stress. If we re-align the averaged stress as the y-axis, we obtain the hourly stress and mass flux in Figure 19g. The integrated averaged spirals in the atmosphere and ocean are not as flat as suggested by Brown (1970) [96], Price et al. (1986) [94] and others. Polton et al. (2013) [97] believed that observed spiral is not flatten if the geostrophic shear can be removed from the data.
The velocity in Figure 19h shows a weak velocity at night and before sunrise, a strong diurnal jet trapped in the upper mixed layer due to a strong surface stress and weak viscosity in the ocean in the daytime. The Coriolis force becomes dominated as the PBL force decays away from the boundary. The rotation of velocity changes 180° in 12 h (half of inertia period); the small velocities cancel out each other near 30 m ( δ E o = 33   m ) and below. Hence, the magnitude decreases from 8.9 cm s−1 to 1 cm s−1, but the angle changes is no more than 90 degrees from surface to −32 m. It usually took from weeks to months to collect observed data (Chereskin 1995 [98], Price and Sundermeyer 1999 [99], and others), which cancels out oscillations with the period less than weeks or so. Hence, many scientists proposed a large viscosity for rotation but small value for the magnitude of the mean velocity (Lenn et al., 2009 [100], etc.). Although the magnitude of averaged spiral is small, the parcels rotate with large amplitude near and beyond the Ekman depth in both atmosphere (blue and green lines with * in Figure 19c) and ocean (thin black at −40 m in Figure 19e). Therefore, the change in magnitude and rotation of velocity of the mean Ekman spiral cannot be interpreted by viscosity alone because the Coriolis can be important when the time to take the average is comparable or longer than the inertia period. Furthermore, viscosity should not be a complex number either, as it is used by many scientists.
The resonance of f and PBL forcing explains the frequent development of the nocturnal low-level jet (LLJ), dryline, and severe storms in the lee of the Rocky Mountains around 30° N (Fujita 1958 [101], Wexler 1961 [102], Palmer and Newton 1969 [93], Sun and Ogura 1979 [103], and Sun and Wu 1992 [104]).

3.4.2. Simulation at 60° S with Diurnal Variation in Atmosphere and Ocean (Case G of Sun and Sun 2020 [91])

We set δzo_ max = 60 m, μo_ max = 0.5 m2 s−1, and δEa = 2νa_anl/f = 252 m in the atmosphere and the depth of Ekman thickness, δEo = 2νo_anl/f = 89 m at 60° S. When the direction of daily mean stress is aligned with the observed stress (Figure 20a,b), the average oceanic spiral in Figure 20b is comparable with the observation at Drake Passage between z = −26 and z = −90 m (Lenn and Chereskin, 2009 [100]), and Polton et al. (2013) [97].
As discussed before, the time-averaged Ekman spiral is different from the steady state solution of (58). It is unlikely to find a viscosity to fit the observed average profile. Hence, “Diffusion theory is a convection and not a testable scientific theory” (Price and Sundermeyer 1999) [99] is incorrect.

4. Purdue Regional Climate Model (PRCM)—A Hydrostatic Model

The linear equations cannot predict the amplitude of disturbances, unless they are driven by a given force and the equations become inhomogeneous, as in the sea breeze circulation (Sun and Orlanski 1981a [87]), or Ekman layer simulation (Sun and Sun 2020 [91]). Hence, we need the nonlinear atmospheric and oceanic models.

4.1. Basic Equations

The PRCM, a hydrostatic primitive equation model, utilizes the terrain-following coordinate (σ) in the vertical direction:
σ = p p t p * = p p t p s p t
where p, pt, and ps are pressure, pressure at domain top, and surface pressure, respectively.
u t = Adv ( u ) + fv m ρ p x + m p * p σ φ x + Diff ( u )
v t = Adv ( v ) + fu m ρ p y + m p * p σ φ y + Diff ( v )
where   φ   ( geopotential ) = g z , f = 2 Ω sin φ , and m is a mapping factor.
Adv = mu x mv y σ · σ
The ice equivalent potential temperature θ ei is
θ ei = θ + θ T L v c p q v L f c p q i
and
θ e i t = A d v θ e i + L v c p q v L f c p q i d d t θ T + D i f f θ e i + R a d θ e i + C o n v θ e i
The change in the surface pressure is calculated by
p * t = 0 1 σ p * V d σ
where
σ p * V = m 2 up * m x + vp * m y
The vertical velocity in the sigma coordinate is
σ · = 1 p * 0 σ σ · ( p * V ) d σ + σ p * 0 1 σ · ( p * V ) d σ
The hydrostatic equation is
φ ln   p = R d T 1 + 0 . 61 q v q l q i
The total water content qw, liquid water ql and ice qi are given by
q w t = Adv q w + Diff q w + Conv q w + Micro q w
q l t = Adv q l + Diff q l + Conv q l + Micro q l
q i t = Adv q i + Diff q i + Conv q i + Micro q i
The sub-grid scale turbulence kinetic energy, E ¯ , is governed by:
E ¯ t = x i u i ¯ E ¯ u i u j ¯ u i ¯ x i + g θ o w θ v ¯ u i e + p / ρ o ¯ x i ε
where e 1 2 u 2 + v 2 + w 2 and E ¯ e ¯ . The sub-grid fluxes were parameterized by
u i u j ¯ = K m u i ¯ x j + u j ¯ x j ,
w i θ i e ¯ = K h θ e ¯ z γ c ,
w q w ¯ = K h q w ¯ z ,
w θ v ¯ = A 1 w θ e i ¯ + B 1 w q w ¯ ,
w e + p / ρ o ¯ = 3 K m E ¯ / z ,
where K m is the sub-grid scale eddy coefficient for momentum, K h the sub-grid eddy coefficient for scalar quantities. The coefficients, A and B in (76) for unsaturated air are given by:
A 1 = 1 + 0.61 q w
B 1 = 0.61 θ L v c p θ T 1 + 0.61 q w
And for saturated air,
A 1 = R v T 2 1 + 1.61 q s q w + 1.61 T L v q s R v T 2 + L v 2 c p q s
B 1 = α i L f c p θ T A θ
Cloud water and cloud ice can coexist in this model. The averaged saturation specific humidity is given by
q s = α i q s i + 1 α i q s w
where
α i = q i q i + q l
The average latent heat of saturated water vapor can be written as
L q s = α i L s q s i + 1 α i L v q s w
The cloud water flux is diagnosed from the equivalent potential temperature flux and total water flux by
w q l ¯ = A 2 w θ e i ¯ + B 2 w q w ¯
where
A 2 = ( 1 α i ) T θ L q s R v T 2 + L v c p + α i L f c p L q s
B 2 = ( 1 α i ) ( R v T 2 + L v c p L q s ) R v T 2 + L v c p + α i L f c p L q s
Similarly, the cloud ice flux can be expressed as
A 3 = α i T θ L q s R v T 2 + L v c p + α i L f c p L q s
B 3 = α i ( R v T 2 + L v c p L q s ) R v T 2 + L v c p + α i L f c p L q s
The covariance of the vertical velocity is parameterized as
w 2 ¯ = 2 3 E ¯ ,
And the eddy viscosities are
K m = 0.10 l E ¯ 1 / 2 ,
K h = ( 1 + 2 l / Δ z ) K m 3 K m if   l < Δ z otherwise
where l is the mixing length scale of the ensemble average turbulence and Δ z is vertical grid space. Following Sun and Ogura (1980) [105], Deardorff (1980) [106], the length scale in the stable atmosphere is assumed:
l = 0.6 E ¯ 1 / 2 g θ o w θ v ¯ K h 1 / 2 if   l < Δ z Δ z otherwise .
In an unstable PBL, according to Caughey and Palmer (1979) [107], it is
l = 0.25 1.8 Z i 1 exp ( 4 z / Z i ) 0.003 exp ( 8 z / Z i ) .
And the upper limit of Zi ≤ 530 m. Dissipation rate ε is given by
ε = 0.41 E ¯ 3 / 2 l
The model has the prognostic equations for the momentum, surface pressure, ice-equivalent potential temperature θei, turbulent kinetic energy (TKE), vapor, cloud, ice, rain, snow, and graupel, etc. (Sun and Chang 1986a [108], Sun 1986 [109], Chern 1994 [110], Chen and Sun 2002 [111]). It consists of multi-layers of snow, soil temperature, and wetness, etc. (Sun and Wu 1992 [104], Bosilovich and Sun 1995 [112], Chern and Sun 1998 [113], Sun and Chern 2005 [114]). It also includes radiation parameterizations (Chou and Suarez 1994 [115], Chou et al., 2001 [116]), and cumulus parameterizations (Kuo 1965, 1974 [117,118], Anthes 1977 [119], Molinari 1982 [120]). The forward–backward scheme is applied in the Arakawa C grid to permit more computational accuracy and efficiency (Sun 1980, 1984a [38,121]). A fourth-order advection scheme (Sun 1993c [122]) or mass-conserved semi-Lagrangian scheme (Sun et al., 1996 [123], Sun and Yeh 1997 [124], and Sun and Sun 2004 [125]) is used to calculate the advection terms. A local reference is applied to calculate the pressure gradient force, which significantly reduces the error of the pressure gradient terms in the σp coordinate over steep topography (Sun 1995a [126]). Yang (2004a, b) [127,128] has added the transport of dusts and the atmospheric chemistry module, as shown in Figure 21.

4.2. Inland Sea Breeze and Dryline

Sun and Ogura (1979) [103] proposed the “Inland-sea-breeze” hypothesis that the pre-storm convergence, vertical circulation, and the nocturnal low-level jet along dryline (Figure 22) on the lee of the Rocky Mountains are triggered by a strong horizontal soil moisture and temperature gradients. Sun and Wu (1992) [104] successfully simulated the formation and diurnal movement of the dryline. In the late spring and early summer, on a slope terrain with a strong soil moisture gradient under a vertical wind shear, a dryline can form within 12 h in Oklahoma and Texas. They reproduced a deeper intrusion of the moisture field far above the potential temperature inversion (Figure 23a,b), as observed by Schaefer (1974) [129]. In the daytime, due to the strong mixing with the warm, dry air descending from the Rocky Mountains, the moisture on the west boundary of dryline dilutes quickly and results in an eastward movement of dryline. At night without mixing, the cool, moist air on the eastside thrusts into the deep, warm, dry air on the westside (Figure 23a–d), and it also gradually turns northward due to the Coriolis force and forms a nocturnal LLJ on the slope of the Rocky Mountains (Figure 23c), as proposed by Sun and Ogura (1979) [103]. The southerly LLJ can effectively transport warm, moisture northward. Hence, the vicinity of the dryline becomes a favorable zone for convections. A weak instability caused by the resonance of a diurnal oscillation of the PBL, and inertial oscillation of the Coriolis force may enhance the disturbances around 30° N, as discussed in Section 3.4.1 (Sun and Sun 2020) [91].

4.3. Effect of Mountain

4.3.1. Lee Vortices and Mountains

With detailed physics on the non-slip surface, Sun and Chern (1993) [130] confirmed the importance of vertical stretching on the formation of vortices proposed by Sun et al. (1991) [131]. They also reproduced the observed lee vortices, meso-low, and downslope wind on the east of the CMR under a southwesterly monsoon during the TAMEX-IOP-2 (Taiwan Area Mesoscale Experiment intensive observation period 2). The lee vortex formed near the southeast coast at 2000LST 16 May 1987 (Figure 24a,b) moved to the east of Taiwan at 0200 LST 17 (Figure 25a,b). The air moved across the CMR, then sank adiabatically and formed a warm, meso-low to the south of the lee vortex as observed in Figure 26a,b. Hence, it should not have an artificial meso-β front as shown in Figure 26b by Kuo and Chen (1990) [132]. The PRCM also replicated the observed severe front at 00Z 17 May 1987 (Figure 27). This well-developed front extended from Japan and the Japan Sea through Taiwan into the South China Sea and Vietnam, and had the property of the typical Mei-yu front discussed by Hsu and Sun (1994) [133] and others.
The clouds with super-cooled water observed by aircraft and Denver cyclone observed by radar on 13 February 1990 during Winter Icing and Storms Project’s (WISP) Valentine’s Day storm (VDS) were also reproduced by PRCM as shown in Figure 28 (Haines et al., 1997 [134]).

4.3.2. Front Deforms around the Central Mountain Range (CMR) in Taiwan

Sun and Chern (2006) [71] replicated the deformation of a weak cold front by the CMR after 30 h integration. The large-scale patterns are comparable with the ECWMF analysis at 06Z15, 1997 (Figure 29a–c) except that the deformation of the front was missing in the analysis. As the cold air approached Northern Taiwan from the north, a relative high pressure with an anticyclonic circulation built up on the windward side. More oncoming flow was diverted to the northeast of Taiwan than the northwest and caused the difference in the frontal speed between the eastern and western parts of the front. The wind speed also increased as the northeasterly flow was blocked by the mountain range. The speed of the front was comparable to U, the wind behind the front, instead of (U+ g h ), the speed of a Kelvin wave.

4.4. Cyclogenesis and Winter Severe Storms in the USA

Chern (1994) [110] has simulated the severe winter storms in the US in March 1985. The unique features of the event are: (1) started from a fair weather, two surface low pressure centers formed near the Great Basin region, one was initiated at 1200 UTC 1 March near the border of Oregon, California, and Nevada; and the other took shape at 0000 UTC 2 March over Eastern Idaho; (2) an explosive development occurred over Nevada with a deepening rate of 18.5 mb in 24 h; (3) the cyclone over Eastern Idaho also exposed with a pressure fall of 13.7 mb in 24 h; (4) a Pacific block with a Rex-type high–low dipole dominated the upstream flow patterns during the entire event; (5) an upper-level short-wave trough grew rapidly from a relatively undisturbed stage and propagated southeastward to the Great Basin region; (6) the upper-level jet stream underwent a dramatic change in wind direction as the upper-level trough developed; and (7) the storm brought heavy snow in excess of two feet over Montana, South Dakota, and southwest of Minnesota.
A control simulation was conducted with full model physics on a 45 km grid to study the rapid cyclogenesis during 1–3 March 1985. The simulation replicated the development and movement of two surface lows in the mountainous areas (Figure 30 and Figure 31). The central pressures of the simulated lows were within 3 mb of the observed values and the predicted cyclone paths were within 135 km of the analysis. The sharp transition in direction of the 300 mb jet stream during the early development of the upper-level ridge/trough system was well captured by the model. The root mean square errors of sea level pressure were 1.90 mb and 3.47 mb for 24 and 48 h predictions, respectively. The correlation coefficients of the change in sea level pressure between the simulation and ECMWF/TOGA analysis were very high, 0.984 for the 24 h prediction and 0.972 for the 48 h prediction (Figure 30 and Figure 31).
Six sensitivity tests were performed to identify the impacts of atmospheric radiation, surface fluxes, surface friction, latent heat release, and the effect of topography on the event. Results showed: (a) complex terrain was quite important to the structure, location, and movement of the cyclones. (b) Surface friction and latent heating were more important than the surface fluxes and atmospheric radiation over the entire domain, and (c) latent heating was minor, and it only contributed 5.4 mb increase in surface pressure over South Dakota, to the northeast of the surface low center.
The development of the Great Basin cyclones was controlled mainly by adiabatic dynamics, and the upper-level forcing. The positive PV anomaly coincided with a tropopause folding through which air originating from the lower stratosphere penetrating deeply into the mid- and lower troposphere as shown in ECMWF and simulated geopotential at 500 mb (Figure 32). The absolute vorticity increased due to the decrease in static stability as the stratospheric high PV air parcels moved into the mid troposphere. The ageostrophic circulations show that the development of the cyclone over Nevada coincided with the rising branch of the indirect circulation in the left quadrant exit region of the upper-level jet streaks upstream from the developing trough. The subsiding branch of the secondary circulation in the warm side of the jet frontal system was important to the tropopause folding. The heavy snowfall and strong upward motion over South Dakota were due to the interactions of the ascending branch of a thermally direct circulation in the entrance of a polar jet streak, the effect of latent heating, and the lifting of air over sloping terrains and a surface front. The system was far more complicated than the idealized baroclinic instability or cyclogenesis presented in the paper (Eady 1949 [76]) or idealized model simulations (Yildirim 1994 [135]).
It was a fair weather before 1 March and the NMC model failed to predict the explosive development of the storms, which severely affected the areas from Wisconsin to Nevada longer than a week. Meanwhile, the PRCM simulated the development of storms (unstable mode) well. It indicates that a comprehensive model can predict the development of the system near the saddle point as the dynamic system changes drastically.

5. Regional Climate

5.1. 1993 Flood in the USA

Bosilovich and Sun (1999a, b) [136,137] reproduced the large-scale atmospheric features in the summer of 1993, when heavy precipitation caused a long-lived, catastrophic flood in the Midwestern United States (Figure 33). The upper-level jet stream and trough were over the Northwestern United States, the Great Plains LLJ supported the moisture transport, and heavy precipitation occurred in the midwest. The model also reproduced the observed precipitation pattern: synoptic wave-type in June and mesoscale convective-system (MCS) in July, over the heaviest rainfall area with an excessive (≥20%) simulated precipitation. The sensitivity tests show: the increased surface heating caused by a strong dry anomaly induces a large-scale surface pressure perturbation and weakens the low-level jet and moisture convergence within the flood region. Both wet and dry regional anomalies in the Southern Great Plains cause less precipitation in the flood region. The uniform soil moisture of both anomalies leads to a reduction in differential heating, surface pressure gradient, and the low-level jet. The model also indicates that the long-lived, 1993 catastrophic flood may be related to the large-scale system.

5.2. 1988 Drought in the USA

Many scientists have applied numerical models to study the cause(s) of the 1988 drought in the USA shown in Figure 34. However, discrepancies exist among the different hypotheses. Oglesby and Erikson (1989) [138], and Oglesby et al. (1994) [139], using the NCAR-CCM1, suggested that the draught was caused by the persistence of a soil moisture anomaly in the US. However, using the same CCM1, Sun et al. (1997) [140] carried out four experiments, which included those using the climatological SSTs (control case); 1988 SSTs; dry soil moisture anomaly (25% of soil moisture in May in the normal year between 35–50° N in the US); and 1988 SSTs and dry soil moisture anomaly. Our results show that the 1988 SST did not cause the simulated weather pattern over the US to be drier than the control case with the climatological SST. The ensemble mean with the perturbed soil moisture experiment did show less precipitation than the control; however, due to the large variance in the data, the reduction in precipitation was not statistically significant. In the experiment with the soil moisture anomaly and 1988 SST, the most significant differences were obtained, specifically in the reduction of precipitation in the US. This may indicate that the 1988 severe drought may be caused by a combination of dry soil, and a special large scale weather pattern related to the SST, instead of by either SST or dry soil alone.
Using the PRCM with the observed SST and the ECMWF analysis as initial and lateral boundary conditions, Sun et al. (2004) [141] reproduced a strong warm ridge at 500 hPa in North America over a hot, dry land. The monthly precipitation and soil moisture were far below the normal values in the Midwest and the Gulf States, but more than normal in the Rocky Mountains as observed (Figure 34), corresponding to a stronger N. American monsoon than normal. A sensitivity test shows that the monthly precipitation could significantly increase using the saturated soil moisture as the initial condition. The soil would dry up eventually, because the wet soil does not provide positive feedback with the low-level jet. Hence, the change in the large-scale weather pattern might trigger the dry episode by forming a ridge in North America, which gradually cut down precipitation in the Gulf States and midwestern region. The soil became dry and hot, which further intensified the blocking of the warm ridge in the US. Hence, dry soil had positive feedback for the severe drought in the summer of 1988. But the dry soil alone was not the cause. Soil moisture became much drier in July than June, but precipitation amounts returned to normal after mid-July 1988, which was also well simulated. It indicates that PRCM can reproduce a rapid change from an abnormal status to a normal situation, i.e., a jump from one phase space to another phase space in the dynamics system (Lorenz 1963 [9]).

5.3. Snow, Land Surface, and Regional Climate

5.3.1. One-Dimensional Snow–Soil Model (Sun and Chern, 116 [114])

The land surface module of the PRCM without snow/ice or frozen soil has been described by Bosilovich and Sun (1995, 1998) [112,142], Sun and Bosilovich (1996) [143], and Sun (2001) [144]. Here we present the model with snow/ice. According to Jordan (1991) [145], the heat (enthalpy) equation for a potentially freezing becomes:
( H e n Δ z ) t k = ( F k + 1 / 2 F k 1 / 2 )
where
H e n = { ( T 273.15 ) ( c d ρ d Δ ζ d + c w ρ w Δ ζ w + c i ρ i Δ ζ i ) + L e i ρ w Δ ζ w } / Δ z = { ( T 273.15 ) ( C d Δ ζ d + C w Δ ζ w + C i Δ ζ i ) + L e i ρ w Δ ζ w } / Δ z ,
L e i ( = L f ) is latent heat of fusion for ice (3.335 × 105 J/kg); Δz is the thickness of the layer; ρ d ,   ρ w ,   and   ρ i are density of dry soil, water, and solid ice within the layer; c d ,   c w ,   and   c i as well as Δ ζ d / Δ z ,   Δ ζ w / Δ z ,   and   Δ ζ i / Δ z are the specific heat capacity and component for dry soil, water, and ice, respectively. The mean heat capacity of soil is C s o i l = ( C d Δ ζ d + C w Δ ζ w + C i Δ ζ i ) / Δ z . F k ± 1 / 2 is the heat flux at the interface between k and k ± 1 layer
F k ± 1 / 2 = K s o i l ( T s z ) + H f l o w k ± 1 / 2
where the effective thermal conductivity K s o i l is calculated following Johansen (1975) [146], reported by Farouki (1982) [147]; H f l o w is the heat transfer due to water passing through layer k. Without the snow coverage, the heat flux at the soil top is
F 1 / 2 = R n H s 0 H L 0
where R n is net radiation. The downward shortwave radiation R s absorbed at surface is
R s = ( 1 A ) R s s f c t o p
where A is surface albedo, and R s s f c t o p is short wave radiation reaching the surface. The net long-wave radiation at surface is
R l = ε R a ε σ T s f c 4
where σ   ( Stefan Boltzmann   constant ) = 5.67 × 10 8   W   m 2 K 4 . The values of A, R s s f c t o p , and surface emissivity ε (ε = 1 for soil) are provided by field experiment. Both sensible heat flux H s 0 and latent heat flux H L 0 are calculated by similarity equations (Bosilovich and Sun 1995 [112]). When the soil is covered by snow, Fk±1/2, the heat flux inside the snow may include the downward shortwave radiation.
The equation for liquid in a frozen or unfrozen soil is
w t = z K ψ z K z + S W
where S w is the source term, including melt of ice/snow, precipitation, evaporation, and runoff within the layer, K is hydraulic conductivity, and ψ is matrix potential. Following Clapp and Hornberger (1978) [148], soil hydraulic property functions can be defined as:
ψ = ψ a e w s w b K K s = w w s 3 + 2 b
The subscript s denotes the value at saturation, and b is an empirical constant which depends upon the texture of the soil.
ψ a e is the air entry tension (Dingman 2002 [149]). Combining (102) and (103), we obtain:
w b ψ ψ t = z K ψ z K z + S W

5.3.2. Equations for Snow Layer

Based on energy and mass balances, the layer’s number and thickness of one-dimensional snow model can change with time. The temperature, density, water content, and thickness of each layer are calculated separately. Each snow layer may consist of solid snow and liquid water. Hence, the mass of a snow layer ( M s n ) is
M s n = M s + M w
where M s is the mass of solid snow and M w is the mass of liquid water content. The thickness and density of a snow layer are given by
Δ z s n = Δ z s + Δ z w ,
and the mean density is
ρ s n = ρ s Δ z s + ρ w Δ z w / Δ z s n
where subscripts “sn”, “s”, “w” are related to snow layer, solid snow, and liquid water, respectively. The temperature ( T s n ) is derived from equation of enthalpy within snow
( H e n Δ z s n ) t k = ( F k + 1 / 2 s n F k 1 / 2 s n ) ,
where enthalpy is defined as:
H e n = ( C w + C s ) ( T s n 273.15 ) + ρ w Δ z w L l i
The heat capacity of a snow layer ( C s n ) is computed from the heat capacities of solid snow ( C s ) and liquid water ( C w ), weighted by their respective volume fractions:
C s n = C s Δ z s + C w Δ z w / Δ z s n ,
where C w = 4.218 × 10 6   J   m 3   K 1 . According to Verseghy (1991) [150], the heat capacity of solid snow is
C s = 1.9 × 10 6   ρ s / ρ i 0   J   m 3   K 1 ,
where ρ i 0 = 920   kg   m 3 is the density of pure ice. The heat flux at the snow top is
F 1 / 2 = R n H s 0 H L 0
where R n is net radiation. The surface energy balance Equations (99)–(101) are also applied to the snow surface. The infrared emissivity is 0.99 for fresh snow and 0.82 for old snow (Oke 1987 [151]). Snow emissivity = 0.95 is used in this study. The snow albedo is calculated by the equation proposed by Bohren and Barkstrom (1974) [152]:
A s = 1.0 0.206   C ν d s 1 / 2 ,
where mean snow grain diameter ds (in mm) is derived from an empirical formula of snow density ρs:
d s = G 1 + G 2 ρ s 2 + G 3 ( 10 3 ρ s ) 4 ,   ( in   mm )
where C ν = 1.2, G1 = 0.16 mm, G2 = 0.0, and G3 = 110 mm g−4 cm12, and the unit of density of solid snow ( ρ s ) is kg   m 3 (Anderson 1976) [153].
The heat fluxes at the top and bottom of layer k is defined as
F k ± 1 / 2 s n = K s n ( T s n z ) R s s n H f l o w k ± 1 / 2 ,
where K s n is the effective thermal conductivity; R s s n is the downward solar radiation inside snow; H f l o w is the heat transfer due to water passing through the layer. The effective thermal conductivity of a snow layer is calculated from volume weighted thermal conductivities of solid snow ( K s ) and liquid water ( K w ):
K s n = K s Δ z s + K w Δ z w / Δ z s n ,
where the thermal conductivity of liquid water K w is set to 0.6   W m 1 K 1 . The effective conductivity of solid snow represents the heat transport through both conduction and vapor diffusion. The empirical formulation of Yen (1965) [154]
K s = 3.221 × 10 6 ρ s 2 ,
Equation (117) is used for the effective conductivity of solid snow. Here, the unit of density of solid snow ( ρ s ) is kg   m 3 .
The change in solar radiation inside kth layer is
R s s n k + 1 / 2 = R s s n k 1 / 2 exp β Δ z s n k ,
where β is the bulk extinction coefficient. Following Jordan (1991) [145], the incident solar radiation is divided into near-infrared and visible components. Infrared radiation is assumed to be totally absorbed within the top layer. The asymptotic bulk extinction coefficient is used for the attenuation of visible radiation, which is computed from the formula proposed by Bohren and Barkstrom (1974) [152], and Anderson (1976) [155].
β = 100 C v ( 1.0 3 ρ s ) d s ,   in   ( m 1 )
The unit of ds is mm in (119). The melting and freezing processes in the model are based on the conservation of mass and heat. Should the calculated new snow temperature from Equations (108) and (109) be greater than the melting point, the excess energy goes toward melting the snow, so that T s n remains at 273.15 K until snow completely melts. Freezing occurs when liquid water exists while snow temperature is below 273.15 K. The changes in density and thickness due to the melting and freezing processes are also considered. The change in mass within each layer is
( ρ s n Δ z s n ) t = { ρ p r e c i p i t a t i o n Pr E e v a p } s f c + { net   water   flow }
where { ρ p r e c i p i t a t i o n Pr E e v a p } s f c is the net change in the mass at the surface due to evaporation Eevap and precipitation ρ p r e c i p i t a t i o n Pr at surface. If the wet-bulb temperature Tw is less than 0 °C, precipitation is assumed to be snow, otherwise is rain (Anderson 1976) [155]. The density of new snow is calculated by
ρ n s = 1000 [ 0.09 + 0.0017 ( T w 258.16 ) 3 / 2 ]   ( in   kg   m 3 )
It is noted that 0.09 is used in (121) instead of 0.05 proposed by LaChapelle (1969) [156].
Following Lynch-Stieglitz (1994) [157], we assume that the maximum liquid water holding capacity is 5.5% by height of a compacted snow layer thickness. Should the liquid water content exceed the capacity, the excess water flows to the next lower layer. In the lower layer, the liquid water may freeze, remain in the layer, or flow through to the next lower layer. Liquid water leaving the bottom of the snow either infiltrates into ground or leaves the system as surface runoff. The heat transport due to the water flowing through a snow layer is also considered as shown in Equation (115).
The mechanical compaction is adapted from Pitman et al. (1991) [158], and Lynch-Stieglitz (1994) [5]:
ρ s n + 1 = ρ s n + 0.5 × 10 7 ρ s n g N k × exp 14.643 4000 min ( T , 273.16 ) 0.02 ρ s n   Δ t
where g N k is the weight of the snowpack above the midpoint of layer k.
The number of snow layers varies from 0 (no snow) to 6. The depth of each layer changes with time according to the conservation of mass and heat. Furthermore, the thickness of the top and bottom layer is always less than or equal to 6 cm to simulate the drastic change in snow property near the interfaces. The 5-year simulations of snow depth, temperature inside snow and beneath are compared excellently with observations shown in Figure 35, and the predictable time scale of the stiff nonlinear soil–snow model can be longer than a year. The snow–ice model was also coupled with a 1D mixed layer ocean model to simulate the sea ice at the SHEBA Experiment in 1997–1998 (Sun and Chern 2010 [159]).

5.4. Flooding Due to Snow Melt

(Min 2005 [160]), Min and Sun (2015) [161] applied the PRCM with the multilayer of snow–land–surface (SLS) module of Sun and Chern (2005) [114] to study the 1997 spring flooding in Minnesota and the Dakotas due to snowmelt. In the winter of 1996–1997, blizzard from the second half of November through January built up an enormous snowpack; many areas had more than 250 cm of snowfall, 2–3 times the normal annual amount. The frigid conditions were throughout much of the winter and early March in 1997, delaying the onset of snowmelt. Significant melt of the deep snow started with particularly warm conditions at the end of March and early April. Many rivers in South Dakota, Southern Minnesota, and Southern North Dakota were rising well above flood stage. The left panels of Figure 36 show the observed snow depth (top), simulation with a single snow layer (middle, referred as CTL), and with comprehensive snow–soil physics (bottom, referred as EXP 1) in March of 1997. The right panels show the observed wind at 200 hPa (top), the simulation from single snow layer (middle), and from comprehensive snow–soil model (bottom). The sea level pressure simulated with comprehensive snow–soil model (EXP1—solid black line) is also more accurate than the simple snow–soil model (red circle) compared with the ECMWF analysis (green circle) in Figure 37. Including cryosphere model physics lowers the surface temperatures due in part by the initial frozen soil conditions and the reduction of incoming solar radiation at the snow surface because of a higher albedo. These affect the horizontal temperature gradients, and in turn change the location of synoptic weather systems over the cold land region. Hence, an accurate snow–land model is crucial to short term weather forecasting or long-term climate study. Charney and Shukla (1981) [15] showed that the tropical sea surface temperatures or land surface soil moisture may evolve slowly or in predictable ways. Consequently, they can give a positive bias to the subsequent behavior of the atmosphere, and therefore provide the basis for some predictive power.

5.5. Southeast and East Asian Monsoon

Yu et al. (2004a, b) [162,163] applied the PRCM with 60 km resolution and comprehensive physics to study the regional climate in 10 summers (MJJA) during 1991–2000. Compared with ECMWF reanalysis, the simulated mean sea level pressure is excellent with bias = 0.15 hPa, the RMSE = 0.60 hPa, and the pattern correlation is 95%; the mean temperature has bias = 0.47 K, RMSE = 0.72 K, and correlation = 0.99 (Figure 38a,b). The observed and simulated precipitation shows correlation coefficient for mean precipitation = 35% (Figure 38c). Yu et al. suggested that the excessive precipitation may be caused by the misrepresented terrain or the southwestern lateral boundary of the model, and after the systematic error is removed, the errors reduce significantly, and the correlation becomes 0.84 (Table 3). The vorticity at 850 hPa shows excellent intraseasonal variability with the correlation coefficient as high as 0.9 or more (Figure 39).
The results obtained by Hsu et al. (2004) [164] indicate that the PRCM could simulate the overall characteristics of the 1998 East Asian summer monsoon (EASM) on onset and the daily, seasonal, intraseasonal time scales. On the seasonal time scale, the model tended to produce more precipitation over the land and less precipitation over the ocean. This land–sea contrast in the model was stronger than observed Pacific anticyclone. The over simulated anticyclone by the PRCM, a fundamental problem, was also found in many models.
The northward-propagating rain bands and intraseasonal oscillation events were well reproduced (Figure 40). The model can produce realistic sub-seasonal variability in the inner domain once the proper lateral boundary forcing is provided. It correctly simulated the onset timing and dramatic changes before and after the onset. But the model incorrectly simulated the circulation and precipitation during the onset because it failed to simulate the rapid development of a weak trough in the Northern South China Sea and missed some extreme events that resulted in flooding in China. After the onset period, the model performance became reliable again. This indicates that the model can fix the existing large biases with the proper lateral boundary, and the dynamic system with the strange attractor being locally unstable yet globally stable.
Sun et al. (2011) [165] compared the PRCM simulation with ECMWF on heavy precipitation over Korea and China from 30 July to 18 August 1998. They found that heavy rainfall along the Mei-Yu (Changma) front was maintained by: (1) an anomalous 850 hPa subtropical high, which intruded westward and enhanced the LLJ over Southeastern Asia, (2) a stronger baroclinicity around 40° N over Eastern Asia and a low pressure located over Japan, Korean, and Mongolia, which not only enhanced southwesterly wind along the frontal zone but also blocked the northward movement of the rain band, and (3) excessive evaporation from abnormally wet, warm land. The simulated precipitation ended by 18 August when the subtropical high retreated and the low-pressure in Mongolia also moved away from the Asian continent, which was consistent with the ECWMF reanalysis and observed precipitation estimates from GPCP. High correlations were found for both surface and upper-level atmospheric fields when compared to ECMWF reanalysis for the 20-day means, as shown in Table 4. It is noted that the slower propagation of the simulated low pressure in Mongolia compared to the observed data may be responsible for the deepening of the low and excessive precipitation in that area. It is also suspected that the model simulates too strong a convection caused by the modified Kuo-type cumulus parameterization scheme (Molinari and Dudek, 1992 [166]), which may enhance a cyclonic circulation and strengthen the low pressure. This may lead to stronger convergence of moisture and heavier precipitation.
The well simulations of seasonable variations of 850 vorticity over S. China Sea, the northward propagation of the Monsoon front in East Asia, as well as 1988 drought, 1993 flood, 1997 snow-melted flood in the USA, etc. show that the PRCM can perform the monthly ~ seasonal simulations. It also indicated that the lateral boundary conditions are more important than the initial small-scale disturbances. They are consistent with Paegle et al. (1997) [13] and Hoskins (2012) [16] statement: “Smaller scale phenomena that cannot be represented explicitly may be partially “slave” to the retained scales, like general regions of convection to a front, in which case aspects of their feedback on the retained scales may be well determined by retained scales”. Those large-scale systems with slow time variations usually carry more energy, which can transfer to the smaller scales with high frequencies due to nonlinear interactions. The simplified low-order dynamic models contain the nonlinear interaction among the spectral components in a rudimentary form only. Therefore, they do not have cascade processes linking the spectral region of the large-scale forcing with the dissipation range (Wiin-Nielsen 1998 [14]). They also miss the feedback from the smaller waves to the large waves, which can be important in meteorology. Hence, Shen (2019) [167] includes more modes to Lorenz model and found the stable point in addition to the original strange attractors.

5.6. Simulations of Dust Storm

Yang (2004a) [127] has coupled the PRCM and the dust module as illustrated in Figure 41 to simulate the severe dust storms in Asia in April 1998. The dust module consists of: (a) Dust source function derived by Ginoux et al. (2001) [168] based on 1° by 1° terrain and vegetation data set derived from the advanced very high-resolution radiometer (AVHRR) data (DeFries and Townshend 1994 [169]); (b) Particle size, which is a function of the source region’s soil properties (Tegen and Fung 1995) [170]). Seven size bins (0.1–0.18 μm, 0.18–0.3 μm, 0.3–0.6 μm, 0.6–1 μm, 1–1.8 μm, 1.8–3 μm, and 3–6 μm, with corresponding effective radii of 0.15, 0.25, 0.4, 0.8, 1.5, 2.5, and 4 μm) are applied in the model; (c) The threshold friction velocity, defined as the horizontal wind velocity required to lift dust particles from the surface, is a function of particle diameter (Marticorena and Bergametti 1995 [171]) and soil wetness; (d) Dust emission, which depends on source function, surface wind speed, threshold velocity, and the fraction of each size class; and (e) Transport and removal processes. Dust aerosols in the PRCM-Dust are transported by advection, dispersion, sub-grid cumulus convection, wet and dry depositions. The dry deposition of dust aerosols is assessed by the gravitational settling for each model vertical layer and surface deposition velocity. The removal of dust aerosols by wet deposition is calculated using the model precipitation rate for both stratified and convective clouds.
With a resolution of 60 km × 60 km, Yang (2004a) [127] has integrated the model continuously from 8–24 April 1998 without nudging or restart. The control run is to provide the environment and weather conditions without dusts, and a reference to compare with the case including dusts run (Sun et al., 2013a, b [172,173]). The ECMWF mean geopotential (m) at 850 hPa and wind at 700 hPa during April 8–24 (Figure 42) are well simulated (Figure 43) (Sun et al., 2013a, b) [172,173]. The ending on April 23 is in good agreement (Figure 44). The uplifted dusts reached around 800 hPa or higher over the source region and remained at 3–5 km or higher in the downwind regions. Dusts were transported southeastward with a cutoff low, then moved northeastward before 18 April and then were transported eastward after. They are consistent with the observed trajectory shown in Figure 45. The horizontal distributions of the dust aerosols were consistent with the satellite images, the TOMS Aerosol Index maps, lidars, surface network reports, and the isentropic trajectory derived from the weather map. The radiative forcing of dust aerosols induced a warming in Northern Asia and cooling in Southern Asia. They are similar with KHSS97 (Koepke et al., 1997) [174] and TF95 (Tegen and Fung 1995) [170] data. Figure 46a shows the difference in the simulated surface temperature with and without dust. The difference is +1.81 K in Northeast Asia and −3.03 K in Southeast Asia due to the radiative effects induced by the large amount of dust aerosols in North and Northeast Asia. This is consistent with the results of Myhre and Stordal (2001) [175] as well as the bias of April 1998 from the 10-year mean of the ECMWF reanalysis (Figure 46b). The influence of dust was far beyond the polluted areas and reached all Asia and India. The warming was over the areas with a high concentration of dust in April 1998, because most of the dust stayed above the low clouds or even coexisted with the cirrus clouds. This was also validated by satellite images that showed the dust layer resided on top of low-lying clouds and significantly reduced the cloud albedo, particularly near the ultraviolet (Husar et al., 2001) [176]. The dust layer also trapped longwave radiation from below. The results are different from the fifth generation Penn State/NCAR Mesoscale Model (MM5) and other results. Their simulated mean height of the dusts was lower than observed and resulted in a southward transport of dusts (In and Park, 2002) [177]. It is also noted that many models used nudging, data assimilation to adjust the meteorological field, or restarted the model daily or every other day, and they cannot be used to assess the effect of aerosols on weather or climate.
The PRCM was also used by Wu et al. (2003) [178] to study the radiative aerosols released from Chinshan Nuclear Power Plant in Taipei County, Northern Taiwan during winter.
In addition to the initial condition, a regional model requires the lateral boundary conditions from the general circulation models (GCMs) or reanalysis. The PRCM can replicate the weather systems and show the physics and dynamics involved for a month or longer integrations, because the ECMWF reanalysis provides excellent lateral boundaries. If a regional model is applied to do forecasting, the performance depends on the GCMs which provide the lateral boundary (Sun 2006) [179]. In addition to numerical methods, we should improve our understanding of the physical processes among aerosols, clouds, precipitation, radiation, and the turbulent fluxes from sea, land, and mountains to improve the predictability of a model. Meanwhile, if we cannot duplicate the previous and current weather/climate, it is hard to claim that we can predict the climate in the future.

6. Nonhydrostatic Models

6.1. Boussinesq Fluid Versus Compressible Fluid

The vertical momentum equation is:
δ d   w d   t = 1 ρ   p   z g
where δ = 0 for a hydrostatic model, and δ = 1 for a nonhydrostatic model. In Boussinesq fluid, the frequency of an internal gravity wave becomes:
σ 2 = k 2 N 2 δ k 2 + m 2
It shows that if k 2 m 2 , nonhydrostatic simulations can be represented by hydrostatic models. The nonhydrostatic models need lots of computing time to solve either the time-consuming Poisson’s equation in an anelastic (or incompressible) atmosphere, or acoustic waves with a very small-time interval in the compressible atmosphere. Hence, Hsu and Sun (1991, 1992) [180,181] added the selected diffusion in the hydrostatic model to simulate the hexagon convective clouds over ocean. Laprice (1992) [182] and Yeh et al. (2002) [183] used a hydrostatic pressure as the basis for the vertical coordinate, then incorporated the nonhydrostatic departure from the hydrostatic system as a perturbation. The resulting set of nonlinear equations is solved iteratively. Janjic (2003) [184] used a mass-based vertical coordinate in his nonhydrostatic Meso Model developed at NCEP, and the nonhydrostatic dynamics has been introduced through an add-on module. It is noted that the speed of sound is infinite in an incompressible fluid.
The atmospheric nonhydrostatic models are frequently solved either as a compressible or anelastic fluid. An anelastic fluid that eliminates acoustic waves is required to solve the time-consuming Poisson equation (Kuo and Sun 1976 [85], Sun and Orlanski 1981b [88], and Chen and Sun 2001 [185]). Ogura and Phillips (1962) [186] showed that if the percentage range in potential temperature is small, and the vertical scale is small compared to the depth of an adiabatic atmosphere, the anelastic system reduces to the Boussinesq equations. The Boussinesq approximations include: the density variation is ignored except in the buoyancy term; and the density variation in buoyancy can be represented by
d ln ρ = d ln T ,   or   d ln ρ = d ln θ .
The equation of state of the air is
p = ρ R T ,   or   d ln p = d ln ρ + d ln T
If   d ln p d ln ρ ,   we have
d ln ρ = d ln T ,   or   d ln ρ = d ln θ ,
which is one of the Boussinesq assumptions. But in an adiabatic process, we also have
c v / c p d ln p = d ln ρ .
Hence, the Boussinesq assumption is inconsistent when it is applied to the atmosphere. But Boussinesq approximations simplify the equations and become popular in meteorology and oceanography. The traditional parcel method, the Froude number in mountain dynamics, and CAPE are based on the Boussinesq fluid, their limitations will be discussed in Section 6.5, Section 6.6 and Section 6.7. The dispersive relation of the inertia internal gravity wave is
σ 2 = k 2 + l 2 N 2 + f 2 m 2 ( k 2 + l 2 + m 2 ) .
Equation (129) is also derived from the Boussinesq fluid. Property of anelastic fluid ( ρ t = · ( ρ v ) = 0 ) is similar to the incompressible fluid ( d ρ ρ d t = · v = 0 ). But ρ obtained from the Poisson equation can change due to the change in temperature and velocity field, which may disrupt the conservation of mass. Hence, the compressible fluid may be a better choice. The frequency in an irrotational, compressible atmosphere is
σ 4 σ 2 C 2 g S 0 C 2 + m 2 + k 2 + l 2 + Γ 0 2 + k 2 + l 2 g S 0 C 2 = 0
where Γ 0 = g C 2 R 0 2 = S 0 + R 0 2 , and subscript 0 indicates the environment value. Equation (130) contains two acoustic waves and two internal gravity waves, the latter is corresponding to (124).
The compressible nonhydrostatic model requires a small-step to handle acoustic waves according to the Courant–Friedrichs–Lewy (CFL) criterion. MacDonald et al. (2000) [187] proposed a quasi-nonhydrostatic model by multiplying α (typically the square of the vertical to horizontal aspect ratio) to the hydrostatic terms in the vertical momentum equation to decrease both the frequency and amplitude of gravity waves. It enables us to calculate the vertical equations explicitly in a longer time step but slow down the hydrostatic adjustment process. Klemp and Wilhemson (1978) [188] proposed a horizontal explicit and a vertical implicit scheme (HE-VI). Sun and Hsu (2005) [189] used a parameter on the mass divergence term in the continuity equation to damp the high frequency acoustic waves in the NTU/Purdue nonhydrostatic model (Hsu and Sun 2001 [190], Hsu et al., 2004 [191], Sun et al., 2009 [192], and Hsieh et al., 2006, 2010 [193,194], etc.). Using the Cloud Resolving Storm Simulator (CReSS, Tsuboki and Sakakibara 2007 [195]), Sun et al. (2012) [196] proved the method can significantly increase the time interval without affecting the gravity waves or thermal convections. Sun et al. (2012) [196] also demonstrated that the horizontal explicit and a vertical implicit scheme (HE-VI) can introduce big errors if the horizontal length scale is comparable with the vertical length scale.

6.2. Modified Forward–Backward Scheme with Smoothing (MFBS)

The non-hydrostatic equations are:
u i t + u j u i x j = 1 ρ p x i + g i + F r r i = F i ρ
The buoyancy
g i = g   if   i = 3 0 otherwise
θ t + u j θ x j = 1 c p θ T d q d t + D θ
δ ρ t + ρ u j x j = 0
  p = ρ R T
  θ = T P O P R C P
The 2D linearized equations in (x, z) directions become:
u t + 1 ρ o p x = 0
  w t + 1 ρ o p z + g ρ ρ o = 0
  θ t + w θ o z = 0
δ ρ t + ρ o u x + w z + w ρ o z = 0
  ρ = p C 2 ρ o θ θ o
where primes are deviations from background (with subscript “0”), which are function of height only, and C = ( c p / c v ) R T 0 . The background wind is zero. δ > 1 was introduced in the continuity equation to filter the high frequency acoustic waves. If we define: p ^ = p / ρ 0 ,   θ ^ = ρ 0 θ / θ 0 ,   u ^ = u ρ 0 ,   w ^ = w ρ 0 ,   S 0 = θ 0 θ 0 z ,   and   R 0 = ρ 0 ρ 0 z , we have:
u ^ t + p ^ x = 0
w ^ t + p ^ z + ( g C 2 + R 0 2 ) p ^ g θ ^ = 0
θ ^ t + w ^ S 0 = 0
p ^ t = C 2 δ u ^ x + w ^ z + ( S 0 + R 0 2 ) w ^
Sun et al. (2013c) [197] showed that the 2Δx and/or 2Δz waves require the shortest Δt to satisfy the CFL criterion. To relax the time step limitation, Sun et al. (2013c) [197] apply smoothing on the right-hand side of (145). Thus, equation (145) is replaced by
p ^ t = C 2 δ f ¯ x z
where
f = u ^ x + w ^ z + ( S 0 + R 0 2 ) w ^
and
f ¯ x z = ( f ¯ x ) ¯ z = f p 1 , q + 2 f p , q + f p + 1 , q 4 ¯ z = cos 2 ( k Δ x 2 ) cos 2 ( m Δ z 2 ) f
where k and m are the wave numbers in the x- and z-directions, respectively. Equation (145) becomes
p ^ t = C 2 δ f ¯ x z = C 2 cos 2 ( k Δ x 2 ) cos 2 ( m Δ z 2 ) δ u ^ x + w ^ z + ( S 0 + R 0 2 ) w ^ = C 2 δ ˜ u ^ x + w ^ z + ( S 0 + R 0 2 ) w ^
and
δ ˜ δ cos 2 ( k Δ x 2 ) cos 2 ( m Δ z 2 )
Equation (148) indicates that 2Δx and 2Δz waves are removed from p ^ t . The solutions of (142)–(146) can be assumed:
θ ^ p ^ u ^ w ^ ( p Δ x , q Δ z , n Δ t ) = λ ˜ n Θ P U W e i ( k p Δ x + m q Δ z )
where λ ˜ is the complex amplification factor. A forward in time scheme is applied to the momentum fields and a backward to the pressure and temperature, they become:
λ ˜ 1 Δ t U + i X P = 0
λ ˜ 1 Δ t W + i Z Γ 0 P g Θ = 0
λ ˜ 1 Δ t δ ˜ P C 2 + i λ ˜ X U + λ ˜ i Z + Γ 0 W = 0
λ ˜ 1 Δ t Θ + λ ˜ S 0 W = 0
where Γ 0 = g C 2 R 0 2 = S 0 + R 0 2 , R 0 = ρ 0 ρ 0 z ,   and   S 0 = θ 0 θ 0 z . The eigenvalue for the modified forward–backward scheme with smoothing (MFBS) becomes
δ ˜ C 2 λ ˜ 1 Δ t 4 + λ ˜ 1 Δ t 2 λ ˜ X 2 + Z 2 + Γ 0 2 + δ ˜ g S 0 C 2 + λ ˜ 2 g S 0 X 2 = 0
or
δ C 2 cos 2 ( k Δ x 2 ) cos 2 ( m Δ z 2 ) λ ˜ 1 Δ t 4 + λ ˜ 1 Δ t 2 λ ˜ X 2 + Z 2 + Γ 0 2 + δ g S 0 C 2 cos 2 ( k Δ x 2 ) cos 2 ( m Δ z 2 ) + λ ˜ 2 g S 0 X 2 = 0
Which consists of two acoustic waves and two gravity waves, as discussed in Equation (130). Let us define:
Δ t ˜ δ = 2 δ ˜ / C 2 X 2 + Z 2 + Γ 0 2 + δ ˜ g S 0 C 2
The amplification factor λ ˜ is then unity if Δt < Δ t ˜ δ . In that case, the frequency ω ˜ δ is defined by
λ ˜ = λ ˜ r + i λ ˜ i = | λ ˜ | exp ( i ω ˜ δ Δ t ) = cos ( ω ˜ δ Δ t ) i sin ( ω ˜ δ Δ t )
and
ω ˜ δ = tan 1 ( λ ˜ i / λ ˜ r ) / Δ t
Without smoothing, the values of λ ˜ ,   ω ˜ δ ,   and   Δ t ˜ δ of Equations (153)–(156) are equal to λδ, ωδ, and Δtδ of MFB, the solution without smoothing (Sun et al., 2012) [196]. The detailed amplification factor λ ˜ and phase speed can be found in Sun et al. (2013c) [197].

6.3. Kelvin–Helmholtz Wave with Δx = 10 m and Δz = 5 m from CReSS

Figure 47 shows the simulated potential temperatures θ for Kelvin–Helmholtz instability at t = 320 s from the traditional FB with Δts = 0.01 s; modified FB without smoothing (MFB) with δ = 16, Δts = 0.04 s; modified FB with smoothing (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, and their accuracies are comparable. However, the Δts = 0.08 s for (MFBS) is twice that for (MFB) with the same δ, and four times HE-VI. Sun et al. (2012) [196] also showed that even Co being less than 1, HE-VI result can depart from the FB, and the HE-VI model becomes unstable when Δts ≥ 0.04 s.

6.4. Equations of NTU/Purdue Nonhydrostatic Model

The vertical coordinate of the NTU/Purdue nonhydrostatic model is defined as:
σ = p 0 ( z ) p 0 ( z t o p ) p 0 ( z s u r f a c e ) p 0 ( z t o p )
where p 0 , pressure of a reference atmosphere, is strictly a function of height. Following Kasahara (1974) [198], we convert (131)–(134) to the equations in (x, y, σ) coordinate. The cloud physics, PBL, and turbulence parameterizations are the same as those in the hydrostatic model (Sun 1993a [199]). The advection terms were solved based on Sun (1993c) [122], or the semi-Lagrangian schemes (Sun et al., 1996 [123], Sun and Yeh 1997 [124], Sun 2023 [200]).

6.5. Numerical Scheme for Advection Equations

Sun (1993c) [122] combined Crowley fourth-order scheme and Gadd scheme to reduce their phase error, by integrating Crowley scheme and Gadd scheme alternately with time step, because the ratio of phase speed to real phase is slightly less than 1 for Growly scheme, but slightly larger than one for Gadd scheme. The initial square step function and simulations after 800 time steps with the Courant = 0.2 for Crowley, Gadd, and Sun’s modified schemes are shown in Figure 48a–d (Sun 1993c) [122], Sun’s is better than the original ones. However, phase error and dispersion still exist. Therefore, the semi-Lagrangian schemes is also used to solve the advection equation.
In Figure 49, the fluid moves from O i ( ξ , η ) to P i ( X , Y ) (indicated by red circle) within Δt in the regular coordinates. We connect those red circles to form the Lagrangian curve (η = yj), which intercepts with the vertical line X = Xk at blue diamond point Q. Because X and Y are well functions of ξ and η, F(Xk) or Y(Xk) can be derived from the coordinates of the neighboring reds or blues by the Lagrangian polynomials (Sun 2023) [200]:
F ( X ) = F i 2 , j ( X X i 1 , j ) ( X X i , j ) ( X X i + 1 , j ) ( X i 2 , j X i 1 , j ) ( X i 2 , j X i , j ) ( X i 2 , j X i + 1 , j ) + F i 1 , j ( X X i 2 , j ) ( X X i , j ) ( X X i + 1 , j ) ( X i 1 , j X i 2 , j ) ( X i 1 , j X i , j ) ( X i 2 , j X i + 1 , j ) + F i , j ( X X i 2 , j ) ( X X i 1 , j ) ( X X i + 1 , j ) ( X i , j X i 2 , j ) ( X i , j X i 1 , j ) ( X i , j X i + 1 , j ) + F i + 1 , j ( X X i 2 , j ) ( X X i 1 , j ) ( X X i , j ) ( X i + 1 , j X i 2 , j ) ( X i + 1 , j X i 1 , j ) ( X i + 1 , j X i , j )
Y ( X ) = Y i 2 , j ( X X i 1 , j ) ( X X i , j ) ( X X i + 1 , j ) ( X i 2 , j X i 1 , j ) ( X i 2 , j X i , j ) ( X i 2 , j X i + 1 , j ) + Y i 1 , j ( X X i 2 , j ) ( X X i , j ) ( X X i + 1 , j ) ( X i 1 , j X i 2 , j ) ( X i 1 , j X i , j ) ( X i 2 , j X i + 1 , j ) + Y i , j ( X X i 2 , j ) ( X X i 1 , j ) ( X X i + 1 , j ) ( X i , j X i 2 , j ) ( X i , j X i 1 , j ) ( X i , j X i + 1 , j ) + Y i + 1 , j ( X X i 2 , j ) ( X X i 1 , j ) ( X X i , j ) ( X i + 1 , j X i 2 , j ) ( X i + 1 , j X i 1 , j ) ( X i + 1 , j X i , j )
After obtaining F and Y of all intercepts (blue diamonds) on vertical line at X = Xk, we apply Equation (161) again but on the vertical direction by replacing X by Y to obtain the value F at the Euler grid C. The same procedure is applied to the entire domain. It is noted that following the same procedure, the vertical/(horizontal) Lagrangian curves can be interpolated from the fluids initially coming from the vertical lines, xi/(horizontal lines yj). Interpolation from the Lagrangian curves projected by fluids coming from either a horizontal or vertical curve to the regular grids is referred as “Economic Internet Interpolation”. The interpolation using both vertical and horizontal lines is called “Complete Internet Interpolation” (Sun and Yeh 1997 [124]). They showed that complete interpolation is slightly more accurate than economic interpolation, but the complete version requires twice computing resources than the economic one. It is also noted that Equations (160) and (161) used in this model can be replaced by the fifth or higher order polynomials (Sun 2023 [200]).
The scheme has been applied to simulate an idealized cyclogenesis of Doswell (1984) [201], which has a tangential velocity of the circular vortex:
V (r) = Asech2 (r) tanh(r),
where r is the radius of the vortex; and A = 2.598 so that the maximum value of V equals 1. The initial condition is
F (x, y, 0) = −tanh[(yyc)/δ],
where δ is the characteristic width of the frontal zone, i.e., one half is F = 1, the other half F = −1 initially. The analytic solution is
F ( x , y , t ) = tanh ( y y c ) δ cos ( ω t ) ( x x c ) δ sin ( ω t )
where (xc, yc) is the center of rotation and ω = V/r is the angular velocity. The domain is 10 units long with resolutions of 129 × 129 grid points. An integration time of five units is chosen so that the analytic solution is still resolvable on the low-resolution grid. The numerical simulation for an idealized cyclogenesis after 16-time steps (with δ = 0.05, Co = 4, t = 5, and Error = 0.076) shown in Figure 50a, is in good agreement with the analytic solution shown in Figure 50b (Sun 2023) [200]. The semi-Lagrangian scheme also avoids the phase error or nonlinear instability which may exist in the finite difference or finite volume schemes (Mesinger and Arakawa 1976) [25].

6.6. Nonhydrostatic Model—Parcel Method, Froude Number, Bernoulli, Downslope Wind, and Waves (Sun and Sun 2015, 2019) [60,202]

The parcel method, which has been cited in many texts, assumes: (a) the air parcel rises or sinks adiabatically; (b) the pressure of the parcel is always the same as the ambient pressure at the same level; and (c) the atmosphere is in hydrostatic equilibrium. Unfortunately, (b) and (c) cannot hold at the same time (Sun and Sun 2015, 2019) [60,202].
Queney (1948) [203] and Queney et al. (1960) [204] first derived the mountain waves and downslope wind based on linearized equations. Long (1953a, b) [205,206], Smith (1985) [207] and others proposed hydraulic jump theory to explain the wave breaking and downslope wind. Clark and Peltier (1984) [208], Peltier and Clark (1979) [209] proposed the amplification of resonant lee waves as well as the trapping and subsequent amplification of internal wave beneath the critical layer. More discussions can be found in Lin (2007) [29]. Both theories propose that the high winds occur in the layer below the dividing streamline in hydraulic jump theory, or the critical layer in the resonant amplification. But simulated large amplitude waves develop far beyond the critical layer (where the wind reverses) as shown in Figure 51, because the density decreases with height. It indicates that energy may not be trapped between the critical layer and ground surface. Long’s model states that the increase in kinetic energy (KE) comes from the decrease in potential energy (PE). The simulations show an increase in both KE and PE when the air parcel climbs up the mountain. Then, air descends against buoyancy and produces a warm, low-pressure center on the lee due to adiabatic warming, which enhances the downslope wind and generates a strong updraft.
Sun and Sun (2015) [60] proposed a hypothesis based on conservations of Bernoulli function that the change in kinetic energy can also come from enthalpy (EN) in compressible fluid. The momentum equation can be written as:
V t + V · V + 2 Ω × V = p ρ + g f r
where V = (u, w). If we define dr = (dx, dz) = Vdt and integrate (165) from i to f:
i V t · d r + i V · V 2 + ( 2 Ω + × V ) × V + p ρ + g z · d r = i f r · d r
Because ( × V ) × V is orthogonal to dr, ( 2 Ω + × V ) × V · d r = 0 and c p d T d p ρ = 0 in an adiabatic process, Bernoulli function remains as a constant, for a steady, inviscid flow:
V · V 2 f + g z f + c p T f = V · V 2 i + g z i + c p T i = B .
or
d V · V 2 + c p T + g z = d K E + d ( E N ) + g d z ( P E ) = 0  
or
d V · V 2 + c p d T d z + g c p d z = 0
where EN (specific enthalpy) = cpT and d EN = cpdT. Bernoulli function is the total energy, i.e., the summation of the kinetic energy (KE), enthalpy (EN), and potential energy (PE).
If we define the hydrostatic–adiabatic process as an air parcel moves adiabatically from zi to zf and follows the hydrostatic equation (i.e., dw/dt = 0), i.e., according to the parcel method:
d p = ρ g d z   and   c p d T α d p = 0
We obtain the hydrostatic–adiabatic lapse rate β h a :
β h a = d T d z = g c p
β h a is traditionally called “dry adiabatic lapse rate”. Let us define p = p h + p , p″ is nonhydrostatic pressure (or dynamic pressure), and hydrostatic pressure satisfy p h ρ z g . From (136) for a dry adiabatic process, we obtain:
d ln θ = d ln T R c p d ln p = 0
θ θ z = T T z R c p p h + p p z = 1 T T z + g c p 1 c p p ρ z = 0
and
  d w d t = p h + p ρ z g = p ρ z = 0       If d T d z = g c p p ρ z 0       If d T d z g c p
If the parcel moves adiabatically with the hydrostatic–adiabatic lapse rate, the dynamic pressure is disregarded and w remains constant, that is consistent with V · V 2 f = V · V 2 i in (169).
Sun and Sun (2015) [60] also proved that conventional Froude number in mountain meteorology
F r = U i N h
Equation (175) is based on both hydrostatic and adiabatic approximations. The Froude number can also be derived from a Boussinesq fluid (ρ′/ρ = −θ′/θ) having a hydraulic jump over an obstacle of height h:
F r = U i g h = U i g θ θ h = U i g h θ i + d θ ¯ d z h θ i θ = U i g θ d θ ¯ d z h U i N h
Since an air parcel cannot follow the hydrostatic adiabatic lapse rate, the Froude number is invalid in the atmosphere, but it can be applied in the incompressible water.
Figure 52a shows the simulated x-component wind (shaded color), potential temperature θ (dashed black lines), Bernoulli function B (thick write line), streamlines (thin black line), and pressure p (thin black line) at t = 6 hr. Figure 52b is the same as Figure 52a, except that the shaded colors display temperature. The contours of B, θ, and streamlines are parallel before flows reach the turbulent areas on the lee side, indicating that the flow is steady, adiabatic, and inviscid, and the contours of B or θ can be interpreted as trajectories outside turbulence areas.
We can follow the trajectory of an air parcel along B = 280,000 m2 s−2, which passes 1A (x = 160 km, z = 1099.14 m) in the upstream region, where θ = 291.1 K, T = 268.01 K, U = 11.24 m s−1, and ρ = 0.9733 kg m−3; 1B (x = 320 km, z = 2391 m) over the mountain peak where θ = 291.1 K, T = 254.91 K, U = 31.89 m s−1, and ρ = 0.8588 kg m−3; and 1C (x = 348 km, z = 378 m) with strong downslope wind, U = 74.44 m s−1, θ = 290.7 K, T = 272.30 K, and ρ = 1.0178 kg m−3, as shown in Table 1 of Sun and Sun (2015) [60]. Figure 52a,b also show that near the mountain top (d PE > 0), the wind becomes stronger (d KE > 0), but the temperature is lower than surrounding. Both PE and KE increase but EN decreases; and the parcel does not follow the hydrostatic adiabatic lapse rate β h a .

6.7. Boulder Severe Downslope Windstorm

Severe downslope windstorm occurred during 11 and 12 January 1972. Here, the observed soundings at Grand Junction, CO are used at inflow. On a free-slip surface without the effect of the Earth’s rotation, the 3 h simulations from 11-hydrostatci models (Doyle et al., 2000) [210] show that all models produce a significant leeside windstorm, but the wind maximum at the 3 h time varies from 86 m s−1 to 43 m s−1 and substantial variations exist in the movement of the jump. However, on non-slip surface without rotation, the simulated downslope wind decreases with time after 2–3 h integration. At z~25 m, u = 35 m s−1 at t = 2 h and u = 27 m s−1 after t = 6 h (Sun and Hsu 2005) [189]. Furthermore, the updraft moves farther downstream as time increases (Sun 2013) [211], and Figure 45 in Lin 2007 [29]). As Lilly (1978) [212] pointed out that the air flow in the upstream changed little during 11 and 12 January, and a surge of warmer air moved across the Northwestern US at 500 hPa, accompanied by rapid pressure falls and strong surface cyclogenesis in the lee of the Rocky Mountains, which triggered outbreak of the downslope windstorm. Hence, Sun (2013) [211] proposed that the northwesterly upper-level jet is unbalanced by the large-scale pressure gradient, the Coriolis force decreases the westerly wind and increases the northerly wind of the upper-level jet as it approaches mountain. Convergence forms in the upper layer and forces air to descend and form a severe, long-lasting downslope windstorm. The 2D equations of Sun (2013) [211] are:
u t + A d v ( u ) = f ( v δ 2 v g ) + D i f f ( u )
v t + A d v ( v ) = f ( u δ 2 u g ) + D i f f ( v )
w t + A d v ( w ) = 1 ρ p z + b u g + D i f f ( w )
Here f = 2 Ω sin ϕ ; b = 2 Ω cos ϕ ; ϕ = 36° N; Other equations are the same as in Section 6.1. It is assumed that f u g = 1 ρ p y b a c k g r o u n d and   f v g = 1 ρ p x b a c k g r o u n d ; and δ2 is set to be 0, or 0.5, or 1.0. The background θeg and (ug, vg) satisfy the thermal wind equation initially. At t > 0 and away from inflow, the geostrophic adjustment occurs, if δ2 ≠ 1. The detailed turbulent kinetic energy, parameterizations, and numerical schemes can be found in (Hsu and Sun, 2001 [190], Sun and Hsu, 2005 [189], Sun et al., 2012 [196]).
The Case D of Sun (2023) [211] is presented. The time-independent observed wind is used at inflow (x = 0, black line in Figure 53a) on non-slip surface with f ≠ 0, b ≠ 0. The red line D in Figure 53a shows the profiles of westerly winds at t = 4 h on the windward side at the mountain foot (x = 300 km). With δ2 = 0, at x = 300 km, the westerly wind of the upper-level northwesterly jet decreases with time as flow approaches the mountain, because d u d t f ( v δ 2 v g ) = f v , and d v d t f ( δ 2 u g u ) = f u with u > 0 and v < 0. Convergence in the upper layer enhances the downward motion over the mountain. Consequently, a strong downslope wind develops on the lee side, as shown in Figure 53b,c, which reaches about 55 ms−1 at t = 2 h, and about 60 m s−1 at 6 h. The maximum downdraft reaches −18 m s−1 at t = 2 h, and −27 m s−1 at t = 6 h; and the peak updraft reaches 30 m s−1 at t = 2 h and 33 m s−1 at t = 6 h. They are comparable with the magnitude of 30 m s−1 of the vertical motion observed over the turbulent zone at z~6 km. The simulated westerly wind varies from −2 m s−1 to 52 m s−1 in 320 km < x < 360 km at z~6 km is also comparable with observations (−5 ms−1 < u < 55 ms−1) presented in Figure 13 of Lilly (1978) [212]. A strong/weak westerly wind corresponding to a cold/warm potential temperature is consistent with observations as well. The vertical velocities and downslope wind are much stronger than the cases with f = 0 or geostrophic balance cases. The hydraulic jump locates at x~340 km. Waves, rotors, and turbulences appear in downstream regions. A zone of the reversed and weak winds develops above the strong downslope wind and on the left-hand side of a deep hydraulic jump. The existence of the reversed flow (u < 0) near z~6 km is consistent with observation, and also turbulences on the lee side and upper layers of the mountain. As the downslope wind intensifies with increasing time, the zone of the reversed/weak winds broadens and descends. On the lee side, the waves, rotors, and turbulences gradually develop over a deep layer from surface to 10 km, and in the upper layers (z > 16 km); meanwhile, strong wind zones exist near θv~310 K (z~4 km) and θv~400 K (z~14 km) at t = 6 h (Figure 53b). The westerly wind at z ~ 25 m is about 43 m s−1 at t = 2 h, and 48 ms−1 at t = 6 h. They are comparable with observations (~50 m s−1). Furthermore, the intensive downslope wind and strong vertical motions last for the entire 9 h integration, and the hydraulic jump does not propagate downstream. They are different from the case with a free-slip surface.

6.8. Convective Available Potential Energy (CAPE) (Sun and Sun 2019) [202]

Sun and Sun (2015, 2019) [60,202] proposed that an increase in KE and PE can come from a decrease in EN according to conservation of Bernoulli function in dry atmospheres. Equations with moisture condensation are:
c p d T d t α d p d t = L d q s d t D h ,
c p d ln θ d t = 1 T L d q s d t + D h ,   and   θ e = θ + θ T q s c p ,
where θ e is the equivalent potential temperature.
V t + V · V + 2 Ω × V = p ρ + g f r
Integrating along the parcel trajectory, we obtain
i V t r · d r + i d V · V 2 + g z + c p T + L q s = i D h + f r · d r
The Bernoulli function with condensation in an inviscid, steady flow becomes:
Bernoulli   function   ( B ) = V · V 2 + g z + c p T + L q s = constant ,
i.e., d B = 0 = d (KE + PE + EN) = 0 for an inviscid, steady flow, where the moist EN (enthalpy) = c p T + L q s . Following (172)–(174) and putting (180) into the vertical momentum equation, we obtain
d w d t = p ρ z g = L q s T T z c p T z g = c p T z 1 + L c p q s T g = c p 1 + L c p q s T T z + β h s = p ρ z
where p″ is nonhydrostatic pressure and the hydrostatic moisture adiabatic lapse rate
β h s g c p / 1 + L c p q s T
Both (184) and (185) are identical to the dry case, except the dry EN = c p T being replaced by the moist EN = c p T + L q s . It also confirms that vertical velocity does not change if it follows hydrostatic saturate adiabatic lapse rate. It is noted that (174) or (185) does not require a steady state condition (i.e., t = 0 ). The change in momentum is driven by the nonhydrostatic pressure, which comes from the parcel departing from hydrostatic adiabatic lapse rate or hydrostatic moisture adiabatic lapse rate. More discussion can be found in Sun and Sun (2019) [202].
The convective available potential energy (CAPE), one of the most popular indicators in predicting severe weather, is the integration of the buoyancy (Buoy = g θ / θ ¯ ) of a parcel from free convection (LFC) to the equilibrium level (EL): In steady Boussinesq fluid, it becomes:
C A P E = z L F C z E L g ρ ρ ¯ d z = z L F C z E L g θ θ ¯ d z = z L F C z E L g T T ¯ d z = i d w 2 2
It shows that “the kinetic energy per unit mass of the outflow is equal to the energy per unit mass of the inflow plus the work conducted by buoyancy forcing during adiabatic ascent” (Moncrief and Green 1972) [213]. But CAPE of the hot towers observed in the Equator trough is negative (Riehl and Malkus 1958) [214]. Strong updrafts are also frequently observed in tropical cyclones without CAPE (Ebert and Holland 1992) [215]. Strong warm downdrafts without evaporation cooling but with positive buoyancy are also observed in the eyewall of hurricane Emily (1987) (Black et al., 1994 [216], etc.), and strong warm, downslope windstorms on the lee of mountains, as discussed previously.
The NTU/Purdue nonhydrostatic model (Hsu and Sun, 2001 [133], Sun and Hsu 2005 [189]), which includes warm cloud physics (Chen and Sun, 2002) [111] and turbulence parameterization (Sun and Chang 1986a, b [108,217]; Sun 1988, 1989 [218,219]), but no radiation with dx = 500 m, and dz ∼75 m, except dz = 25 m for the first layer above the surface was employed. A moist-stable (CAPE = 0) sounding of (98646 RPMT (Mactan Observations) at 00Z 18 June 2017 (Figure 54) was used as the initial temperature and moisture fields, unless the moisture was specified differently. We imposed a weak, vertical velocity wgr at the surface:
w g r ( x , t ) = w a m p tanh ( t / 702 s ) x x c m Δ x 2 + 1

6.8.1. Dry Plume with wamp = 3 m s−1 and m = 30, Case C of Sun and Sun (2019) [202]

Figure 55a shows the B u o y (shaded colors), w (purple line), θ (white), and streamlines (black) at t = 1.5 h. At x = 335 grid, Buoy = −0.04 m s−2 near the surface decreases to −0.62 m s−2 at 4.6 km, but w~3 m s−1 at surface increases to 12.4 m s−1 at 2 km, then deceases to −0.26 m s−1 at 4.9 km (Figure 55a,b). Figure 55c shows that θ (red) and B (green) remain constant from the surface to 4.4 km. Turbulence mixing decreases the values of θ and B at the top of the updraft. Part of EN is used to support the increase in KE and PE from surface to 2 km, as discussed in Sun and Sun (2015) [60]. The maximum adiabatic cooling (negative buoyancy) coincides with the maximum updraft along x = 335 grid in Figure 55a.

6.8.2. Moist Plume with wamp = 1 m s−1 and m = 30, Case D of Sun and Sun (2019) [202]

The observed temperature and humidity of the sounding in Figure 54 are used. Figure 56a shows the Buoy (shaded colors), which includes the loading of cloud water (ql; thick white) and rain (qr; thin white), θe (nearly horizontal white thin line), w (purple), and streamlines at t = 0.56 h. The cloud forms (thick white contour) between 0.4 and 3.7 km, with ql > 0.001. Strong adiabatic cooling occurs in the lower updraft around the center. Inside the cloud, adiabatic cooling is partially cancelled by the latent heat released. Strong negative buoyancy also occurs at cloud top and above. Although the buoyancy is negative, the vertical velocity increases from 1 m s−1 at surface to the maximum of 3.3 m s−1 at 0.8 km at x = 335 grid, where a stronger updraft releases more latent heat according to w ( q s / z ) , and the buoyancy becomes positive in 2.1 km ≤ z ≤ 3 km at t = 0.67 h (red area in Figure 56b). Meanwhile, loading and evaporative cooling from raindrops near and beneath cloud base strengthen negative buoyancy. Consequently, a red, pear-shaped positive buoyancy forms in the middle of cloud (Figure 56b). The updraft increases from 1 m s−1 at surface to 10.1 m s−1 at 1.5 km against negative buoyancy (Figure 56c). At x = 335 grid, B and θe are near constant from surface to 2.5 km (Figure 56c,d). They are not as straight as the dry case shown in Figure 55c due to stronger turbulent mixing as well as evaporation and loading of cloud water and rain.
Figure 56e shows that at t = 1.16 h, negative buoyancy still exists in the updraft below 2.6 km. Above it, positive buoyancy enhances the updraft. The vertical velocity reaches 18.3 m s−1 at 9.3 km (Figure 56e,f). A strong cooling occurs at the dome-shaped cloud top near the tropopause. The property of the updraft is modified above 6 km because of entrainment of dry air, as shown in streamline, θe and B at t = 1.16 h (Figure 56e,f). Figure 6e,f show adiabatic warming on downward motions outside the cloud at z~3 km at x~285 grid and 385 grid; and z~7 km at x~310 grid and x~360 grid, and gaps between clouds (at z~3 km and x~320 grid and x~350 grid. The downdraft is not due to evaporative cooling or loading of rain. More cases can be found in Sun and Sun (2019) [202]. Conservation of Bernoulli function shows that the enthalpy EN can be converted to KE and/or PE so the plume can ascend or descend vertically. Consequently, deep convection develops without CAPE, and downslope can form against positive buoyancy due to the dynamic pressure in (174) and (185).
Hot towers of tall cumulonimbus clouds in the equatorial trough have a horizontally span about 2–4 km across. They can reach altitudes as high as 12–18 km and exhibit high reflectivity. Hot towers are effectively undilute; as a result, the equivalent potential temperature within a hot tower remains nearly constant throughout their entire vertical extent. The release of latent heat caused by condensation of those hot towers supplied the energy necessary to maintain Hadley cells and the trade winds according to Riehl and Malkus (1958) [214]. They found that convections can rise against the stable environment with a temperature lower than that of the surroundings, where CAPE is negative according to Equation (187). The Bernoulli function provides plausible explanations that the cumulonimbus chimneys in the Equator, hurricanes and other severe storms can develop in an environment without CAPE, reach 12–18 km, and transport moisture and aerosols effectively into the stratosphere.
The observations show that weak thunderstorms may have vertical updraft speeds of 6–12 m s−1. At such a speed, a thunderstorm will vertically develop by 10 km in about 15 min. A severe thunderstorm may have a vertical updraft going up at 30–35 m s−1, and an overshooting top lasting for 10 min or longer. If the updraft at 20 m s−1, the parcel takes about 10 min reaching the cloud top at 12 km from the ground, which is shorter than the life span of storms (about an hour or longer). The convection and moisture experiment in 1998 also showed that the equivalent potential temperature within hot towers was virtually constant across their entire vertical extent, confirming the lack of entrainment. Hence, the Bernoulli function, which assumed a steady state and inviscid fluid, should be applicable to storms as well as the flow passing over mountains.
The hot towers in the Equator or deep convections of typhoons can interact with the large-scale systems and influence the large-scale systems in the atmosphere and the ocean. The individual cumulonimbus or storms are difficult to predict, but the development of typhoons or cluster of hot towers depends on the large-scale environment, including sea surface temperature, wind, and the supply of heat and moisture, etc., which are easier than the individual convection to predict.

6.8.3. Nonhydrostatic Pressure Inside a Cloud Model

The importance of nonhydrostatic pressure in an air parcel was discussed in Equations (174) and (185) for both dry and moist air parcels. In a one-dimensional time dependent cloud model including detailed cloud microphysics discussed in Section 4.1, entrainment, and nonhydrostatic (i.e., dynamic) pressure, Chen and Sun (2002) [111] found that the inclusion of the nonhydrostatic pressure can (1) affect vertical velocities, (2) help the cloud develop sooner, (3) help maintain a longer mature stage, (4) produce a stronger overshooting cooling, and (5) approximately double the precipitation amount as shown in Figure 57. It is noted that the dynamics pressure is solved from the continuity equation of an anelastic atmosphere. They also show that the one-dimensional cloud model can reproduce the time evolution of the area averaged of three-dimensional WRF simulations shown in Figure 58.

6.9. Lee Vortices and Hydraulic Jump in White Sand Missile Range (Haines et al., 2019) [220]

The National Taiwan University/Purdue University (NTU/PU) nonhydrostatic model using a modified forward–backward integration scheme (Sun et al., 2013c) [197] to retains internal gravity waves but suppresses unwanted sound waves was use to study the field experiment conducted at White Sand Missile Range (WSMR), New Mexico, USA on 25 January 2004, which was designed to study a downslope windstorm, lee vortices, and hydraulic jumps along the lee of the Organ and San Andres Mountains (Haines et al., 2006). It was set up with a vertical grid spacing of 300 m, and a horizontal resolution of 1 km on a domain of 201 × 201 km with the sounding at EL Paso, Texas at 1200 UTC on 25 January 2004 as the initial condition. After 4 h integration, the simulated flows including hydraulic jumps and vortex shedding in the lee of the Organ Mountains agree well with observations as show in Figure 59. The simulations from the WRF used the Gal-Chen and Somerville coordinate (1975) [221] are presented in Figure 60, which failed to reproduce the vortex shedding as observed (Haines et al., 2006) [222]. Figure 61 shows that the observed wind pattern on January 25 can be reproduced very well by an initial westerly wind of 5 m s−1 with Δx = Δy = 2 km and Δz = 300 after 4 h integration. Furthermore, the train of lee waves is also revealed clearly in this larger domain simulation, which has been reported frequently in this area (Haines et al., 2006) [222].

7. Terrain Following Coordinate in Atmospheric Model (Sun 2021) [223]

Numerical atmospheric and oceanic models usually require terrain following coordinates to handle the complex terrain (Kasahara 1974 [198], Sun and Chern1993 [130], Pacanowski 1995 [224], Saito et al., 2001 [225], Staniforth et al., 2004 [226], Lin 2007 [29], Lin et al., 2018 [227], Sun and Sun 2019 [202], etc.) The coordinates proposed by Gal-Chen and Somerville (1975) [221] (will be referred as GS) have been applied to the Regional Atmospheric Modeling System (RAMS, Pielke et al., 2002) [228], Geesthacht Simulation Model of the Atmosphere (GESIMA, Kapitza and Eppel 1992) [229], Cloud Resolving Strom Simulator (CReSS, Tsuboki and Sakakibara 2007) [195], Japan Meteorological Research Institute (MRI-model, Saito et al., 2001) [225], Weather Research and Forecasting Model (WRF) (Skamarock et al., 2008) [230], and other models. The relationship between their coordinates ( x ¯ ,   y ¯ ,   z ¯ ) and the Cartesian coordinates (x, y, z) satisfies: x ¯ = x ,   y ¯ = y , and   z ¯ = z t z z b z t z b , where z t   and   z b are the height of the domain and the terrain elevation, respectively. This popular coordinate works well when it is applied to the gradient ∇ψ but fails to produce the accurate divergence or the curl over a sloped mountain. It cannot be applied to the Navier–Stokes equations either. We propose a new terrain-following coordinate, in which the spatial intervals are not constant.

7.1. Equations

In a curvilinear coordinate (Wikipedia), a position vector can be presented by
r = x ^ i g i ( = 1 3 x ^ i g i ) = x ^ i g i ( = 1 3 x ^ i g i ) = x i e i
where x i e i is in Cartesian coordinate, g i is the tangent (contravariant) basis vector along the curvilinear coordinate, x ^ i (a contravariant component), x ^ i is covariant component, and covector (or dual basis vector) g i are shown in Figure 62.
d r = r x ^ i d x ^ i = g i d x ^ i = r x ^ i d x ^ i = g i d x ^ i = e i d x i = e x d x + e y d y + e z d z
g i r x ^ i = e j x j x ^ i = e x x x ^ i + e y y x ^ i + e z z x ^ i
and
r x ^ i = g i
They satisfy:
g i · g j = δ i j
x ^ i = r g i = x ^ j g j g i ;   x ^ i = r g i = x ^ j g j g i
and
g 1 = g 2 × g 3 J ,   g 2 = g 3 × g 1 J ,   g 3 = g 1 × g 2 J ,
Jacobin is
J = g 1 g 2 × g 3 = g 1 × g 2 g 3 = g 2 × g 3 g 1
J = g 1 , g 2 , g 3 = det r x ^ 1 , r x ^ 2 , r x ^ 3 = det x x ^ 1 x x ^ 2 x x ^ 3 y x ^ 1 y x ^ 2 y x ^ 3 z x ^ 1 z x ^ 2 z x ^ 3
and
d f = f x ^ i d x ^ i = g i f x ^ i · g k d x ^ k = g i f x ^ i · d r = f · d r
Hence, the gradient operator is:
g i x ^ i = e i x i
and
x ^ i = g j x ^ i x ^ j = g i = e j x ^ i x j
According to the divergence theorem, V d v = V n d A , we can calculate the divergence in the volume of d v = g 1 d x ^ 1 × g 2 d x ^ 2 g 3 d x ^ 3 = J d x ^ 1 d x ^ 2 d x ^ 3 in Figure 63.
The flux of a velocity V normal to the lower surface is
V n d A = V g 1 d x ^ 1 × g 2 d x ^ 2 = V g 1 × g 2 d x ^ 1 d x ^ 2 = J V g 3 d x ^ 1 d x ^ 2 = J v ^ 3 d x ^ 1 d x ^ 2
where v ^ 3 = V g 3 , and the flux on the upper surfaces can be estimated by: J v ^ 3 d x ^ 1 d x ^ 2 + ( J v ^ 3 d x ^ 1 d x ^ 2 ) x ^ 3 d x ^ 3 . Hence, the net flux along g 3 is ~ ( J v ^ 3 ) x ^ 3 d x ^ 1 d x ^ 2 d x ^ 3 . It is also applied to two other directions. The summation of the net fluxes becomes
V n d A = i = 1 3 ( J v ^ i ) x ^ i d x ^ 1 d x ^ 2 d x ^ 3
Therefore:
V d v = V d x ^ 1 g 1 × d x ^ 2 g 2 d x ^ 3 g 3 = V g 1 × g 2 g 3 d x ^ 1 d x ^ 2 d x ^ 3 = V J d x ^ 1 d x ^ 2 d x ^ 3 = V n d A = i = 1 3 ( J v ^ i ) x ^ i d x ^ 1 d x ^ 2 d x ^ 3
and
V = 1 J ( J v ^ 1 ) x ^ 1 + ( J v ^ 2 ) x ^ 2 + ( J v ^ 3 ) x ^ 3
The curl of vector V becomes:
× V = g 1 × x ^ 1 V + g 2 × x ^ 2 V + g 3 × x ^ 3 V = J 1 g 2 × g 3 × x ^ 1 V + g 3 × g 1 × x ^ 2 V + g 1 × g 2 × x ^ 3 V = J 1 g 3 x ^ 1 g 2 V g 3 g 2 x ^ 1 V g 2 x ^ 1 g 3 V + g 2 g 3 x ^ 1 V + g 1 x ^ 2 g 3 V g 1 g 3 x ^ 2 V g 3 x ^ 2 g 1 V + g 3 g 1 x ^ 2 V + g 2 x ^ 3 ( g 1 V ) g 2 g 1 x ^ 3 V g 1 x ^ 3 ( g 2 V ) + g 1 g 2 x ^ 3 V
But
g 1 x ^ 2 = x ^ 2 r x ^ 1 = x ^ 1 r x ^ 2 = g 2 x ^ 1 .
Therefore
× V = J 1 g 3 x ^ 1 ( g 2 V ) g 2 x ^ 1 ( g 3 V ) + g 1 x ^ 2 ( g 3 V ) g 3 x ^ 2 ( g 1 V ) + g 2 x ^ 3 ( g 1 V ) g 1 x ^ 3 ( g 2 V )
Finally, we obtain:
× V = 1 J g 1 g 2 g 3 x ^ 1 x ^ 2 x ^ 3 g 1 V g 2 V g 3 V = 1 J g 1 g 2 g 3 x ^ 1 x ^ 2 x ^ 3 v ^ 1 v ^ 2 v ^ 3
where v ^ i = V g i = v ^ j g j g i .

7.2. New Terrain Following Coordinate

Figure 64 shows the 2D diagram of the curvilinear coordinate ( x ^ 1 , x ^ 3 ) , the inclination angle θ of the z ^ -surface. For convenience, we also define the vertical coordinate z ^ = x ^ 3 and:
z ^ = x ^ 3 = z z b z t z b
where z b , is the terrain elevation and z t is the domain height. From (204), we obtain:
z ^ z x , y = z z b z t z b z x , y = 1 z t z b
and
z x ^ 3 x ^ 1 , x ^ 2 = z z ^ x ^ 1 , x ^ 2 = z b + z ^ z t z b z ^ x ^ 1 , x ^ 2 = z t z b
The inclinations of z ^ along x and y–directions are
z x z ^ = z b + z ^ z t z b x z ^ = ( 1 z ^ ) z b x
and
z y z ^ = ( 1 z ^ ) z b y
The changes z ^ with respect to x and y are:
z ^ x z = z z b z t z b x = z t z b z z b x z z b z t z b x z t z b 2 = z b x z ^ 1 z t z b
z ^ y z = z z b z t z b y = z b y z ^ 1 z t z b
The tangent (contravariant) basis vectors in the terrain coordinate system become:
g 1 = r x ^ 1 x ^ 3 = x e x + y e y + z e z x ^ 1 x ^ 3 = x x ^ 1 e x + e z ( 1 z ^ ) z b x
g 2 = r x ^ 2 x ^ 3 = x e x + y e y + z e z x ^ 2 x ^ 3 = y x ^ 2 e y + e z ( 1 z ^ ) z b y
and
g 3 = r x ^ 3 x ^ 1 , x ^ 2 = x e x + y e y + z e z x ^ 3 x ^ 1 , x ^ 2 = ( z t z b ) e z
At surface ( z ^ = 0 ) , g 1 and g 2 are tangent to the terrain surface, while g 3 always points to the vertical direction. The Jacobian becomes:
J = g 1 g 2 × g 3 = det x x ^ 1 x x ^ 2 x x ^ 3 y x ^ 1 y x ^ 2 y x ^ 3 z x ^ 1 z x ^ 2 z x ^ 3 = det x x ^ 1 0 0 0 y x ^ 2 0 ( 1 z ^ ) z b x x x ^ 1 ( 1 z ^ ) z b y y x ^ 2 z t z b = x x ^ 1 y x ^ 2 z t z b
From (210)–(214), the covector (dual) basis vectors become:
g 1 = g 2 × g 3 J = e x y x ^ 2 ( z t z b ) / x x ^ 1 y x ^ 2 z t z b = e x / x x ^ 1 = x ^ 1 x j e j
g 2 = g 3 × g 1 J = e y x x ^ 1 ( z t z b ) / x x ^ 1 y x ^ 2 z t z b = e y / y x ^ 2 = x ^ 2 x j e j
and
g 3 = g 1 × g 2 J = z b x z ^ 1 e x ( z t z b ) + z b y z ^ 1 e y ( z t z b ) + e z ( z t z b ) = x ^ 3 x j e j
The transformation between the velocity (u, v, w) in the Cartesian coordinate (x, y, z) and the velocity ( u ^ , v ^ , w ^ ) in the new coordinate is
u v w = d x d t d y d t d z d t = x x ^ 1 d x ^ 1 d t + x x ^ 2 d x ^ 2 d t + x x ^ 3 d x ^ 3 d t y x ^ 1 d x ^ 1 d t + y x ^ 2 d x ^ 2 d t + y x ^ 3 d x ^ 3 d t z x ^ 1 d x ^ 1 d t + z x ^ 2 d x ^ 2 d t + z x ^ 3 d x ^ 3 d t = x x ^ 1 0 0 0 y x ^ 2 0 ( 1 z ^ ) z b x x x ^ 1 ( 1 z ^ ) z b y y y ^ 1 z t z b u ^ 1 u ^ 2 u ^ 3 = J u ^ 1 u ^ 2 u ^ 3
and
u ^ 1 u ^ 2 u ^ 3 = d x ^ 1 d t d x ^ 2 d t d x ^ 3 d t = x ^ 1 x x ^ 1 y x ^ 1 z x ^ 2 x x ^ 2 y x ^ 2 z x ^ 3 x x ^ 3 y x ^ 3 z u v w = x ^ 1 x 0 0 0 x ^ 2 y 0 z b x z ^ 1 ( z t z b ) z b x z ^ 1 ( z t z b ) 1 ( z t z b ) u v w = J 1 u v w
and
J J 1 = x x ^ 1 0 0 0 y x ^ 2 0 ( 1 z ^ ) z b x x x ^ 1 ( 1 z ^ ) z b y y x ^ 2 z t z b x ^ 1 x 0 0 0 x ^ 2 y 0 z b x z ^ 1 ( z t z b ) z b x z ^ 1 ( z t z b ) 1 ( z t z b ) = 1 0 0 0 1 0 0 0 1 = J 1 J
It shows that g 3 parallels to ez, g 1 parallels ex, and g 2 parallels ey. In this new system, g i is different from g i .

7.3. Gal-Chen and Somerville Terrain Following Coordinate

Gal-Chen and Somerville (1975) [221] assumed:
x ¯ = x ;   y ¯ = y ;   and   z ¯ = z t z z b z t z b
And the inverse transformation:
x = x ¯ ;   y = y ¯ ;   z = z ¯ z t z b z t + z b
The relationships of the basis vectors and other variables between the systems in Section 7.2 and Section 7.3 are:
z ^ = z z b / z t z b = z ¯ / z t
and
x x ¯ 1 = 1   and   y x ¯ 2 = 1 ,   but   x x ^ 1   or y x ^ 2   can   be   different   from   1 .
Their basic vector g ¯ i and the relationship with g i are:
g ¯ 1 = r x ¯ 1 = x e x + y e y + z e z x ¯ 1 x ¯ 3 = x x ¯ 1 e x + e z x ¯ 1 ( z ¯ z t z b z t + z b ) x ¯ 3 = x x ¯ 1 e x + e z z b x ¯ 1 ( 1 z ¯ z t ) x ¯ 3 = x x ¯ 1 e x + e z ( 1 z ^ ) z b x x ¯ 3 = x ^ 1 x ¯ 1 g 1
g ¯ 2 = r x ¯ 2 = x e x + y e y + z e z x ¯ 2 x ¯ 3 = y x ¯ 2 e y + e z ( z t z ¯ ) z t z b y = y x ¯ 2 e y + e z ( 1 z ^ ) z b y = r x ¯ 2 = x ^ 2 x ¯ 2 g 2
g ¯ 3 = r x ¯ 3 x ¯ 1 , x ¯ 2 = x e x + y e y + z e z x ¯ 3 x ¯ 1 , x ¯ 2 = e z ( z t z b ) z t = r x ¯ 3 = g 3 z t
And the Jacobian
J ¯ = g ¯ 1 g ¯ 2 × g ¯ 3 = det r x ¯ 1 r x ¯ 3 = det x x ¯ 1 x x ¯ 2 x x ¯ 3 y x ¯ 1 y x ¯ 2 y x ¯ 3 z x ¯ 1 z x ¯ 2 z x ¯ 3 = det x x ¯ 1 0 0 0 y x ¯ 2 0 z b x ( z t z ¯ ) z t z b y ( z t z ¯ ) z t ( z t z b ) z t = ( z t z b ) z t
As well as
g ¯ 1 = g ¯ 2 × g ¯ 3 / J ¯ = ( z t z b ) z t e y × z t e z ( z t z b ) = e x = g 1 x ^ 1 x
g ¯ 2 = g ¯ 3 × g ¯ 1 / J ¯ = e y = g 2 x ^ 2 y
g ¯ 3 = g ¯ 1 × g ¯ 2 / J ¯ = x x ¯ 1 e x + e z z b x ¯ 1 ( 1 z ¯ z t ) x ¯ 3 × y x ¯ 2 e y + e z ( z t z ¯ ) z t z b x ¯ 2 / ( z t z b ) z t = z t ( z t z b ) e z ( z t z ¯ ) ( z t z b ) z b y e y ( z t z ¯ ) ( z t z b ) z b x e x = z t z ^ 1 ( z t z b ) z b x e x + z ^ 1 ( z t z b ) z b y e y + 1 ( z t z b ) e z = z t g 3
The velocity turns out:
u ¯ v ¯ w ¯ = d x ¯ 1 d t d x ¯ 2 d t d x ¯ 3 d t = x ¯ 1 x x ¯ 1 y x ¯ 1 z x ¯ 2 x x ¯ 2 y x ¯ 2 z ( z ¯ z t ) z b x ( z t z b ) ( z ¯ z t ) z b y ( z t z b ) z t ( z t z b ) u v w = 1 0 0 0 1 0 ( z ¯ z t ) z b x ( z t z b ) ( z ¯ z t ) z b y ( z t z b ) z t ( z t z b ) u v w
GS-derived Equation (232). It is noted that:
u ¯ = d x ¯ 1 d t = d x d t = u = u ^ x x ^ 1 ,
v ¯ = d x ¯ 2 d t = d y d t = v = v ^ y x ^ 2 ,
Hence, their system becomes quite simple, and
w ¯ = d x ¯ 3 d t = ( z ¯ z t ) z b x ( z t z b ) u + ( z ¯ z t ) z b y ( z t z b ) v + z t ( z t z b ) w = z t d z ^ d t = z t ( 1 z ^ ) z b x ( z t z b ) u ( 1 z ^ ) z b y ( z t z b ) v + 1 ( z t z b ) w = z t u ^ 3
The corresponding gradient, the divergence, and curl are
S = g ¯ i S x ¯ i ,
V = 1 J ¯ ( J ¯ v ¯ 1 ) x ¯ 1 + ( J ¯ v ¯ 2 ) x ¯ 2 ( J ¯ v ¯ 3 ) x ¯ 3 ,
and
× V = 1 J ¯ g ¯ 1 g ¯ 2 g ¯ 3 x ¯ 1 x ¯ 2 x ¯ 3 g ¯ 1 V g ¯ 2 V g ¯ 3 V = 1 J ¯ g ¯ 1 g ¯ 2 g ¯ 3 x ¯ 1 x ¯ 2 x ¯ 3 v ¯ 1 v ¯ 2 v ¯ 3

7.4. Numerical Simulations

The 2D model consists of two sets of grids: The first one is the C-grids in the Cartesian coordinate without mountain. The second set is derived from the current and GS coordinates. The Cartesian 2D model includes (20 × 20) uniform grids with Δx = 1 km, zt = 12 km, and Δz = 12 km/20. The elevation of a bell-shaped mountain is given by:
z b = z b m / x x c w a 2 + 1.0
where zbm = 2000 m, wa = 4000 m, and xc is located at the mid-point in the x-axis. The coordinate z ^ is given by (204) and z ¯ by (221).

7.4.1. Gradient

A scalar variable,
S = C s ( x x c ) z
With Cs = 3.0. The second-order centered difference in the gradient in the Cartesian coordinate is:
S a = S x e x + S z e z = S j + 1 / 2 , k S j 1 / 2 , k Δ x e x + S j , k + 1 / 2 S j , k 1 / 2 Δ z e z
The numerical scheme in the current terrain following coordinate is:
S b = g j S x ^ j = e x x / x ^ 1 S x ^ 1 + z b x z ^ 1 e x ( z t z b ) + e z ( z t z b ) S x ^ 3
The scheme in GS is
S c = g ¯ j S x ¯ j = e x x / x ¯ 1 S x ¯ 1 + z ¯ z t ( z t z b ) z b x e x + z t ( z t z b ) e z S x ¯ 3
The continuous black lines in Figure 65 show the gradient in the Cartesian coordinate of ∇Sa along e x   and   e z directions, respectively. They are completely covered by the dashed red lines from ∇Sb and the blue dots of ∇Sc, i.e., ∇Sb = ∇Sc.

7.4.2. Divergence

A vector
V = D c ( x x c ) z 2 e x + D c ( x x c ) 2 z e z
With Dc = 6 × 10−8. We solve the divergences numerically:
D i v a = V = v x + w z ,
D i v b = V = g j x ^ j ( v ^ i g i ) = 1 J x ^ j ( J v ^ j )
and
D i v c = V = g ¯ j x ¯ j ( v ¯ i g ¯ i ) = 1 J ¯ x ¯ j ( J ¯ v ¯ j )
Figure 66 shows that Divb (dashed red lines) agrees with Diva (black lines). But Divc (blue dots) departs from Diva (black lines) over the slope of the mountain, where the difference between δ x ¯ 1 and   δ x ^ 1 is large.

7.4.3. Curl along y-Direction

The vector V is given by
V = 2.526 × 10 6 ( x x c ) z e x + 2.526 × 10 3 ( x x c ) ( z 10 Δ z ) 2 e z
The y-component of curls becomes:
C u r l a = u z w x ,
C u r l b = 1 J v ^ 1 x ^ 3 v ^ 3 x ^ 1 ,
and
C u r l c = 1 J ¯ v ¯ 1 x ¯ 3 v ¯ 3 x ¯ 1
The results of Curla (black), Curlb (dashed red), and Curlc (blue dots) are shown in Figure 67.
Overall, Curlb resembles Curla, except some errors in the shortwaves. But Curlc can be quite different from Curla.

7.4.4. The Navier–Stokes Equations in New Terrain-Following Coordinate

As discuss Sect 6, the advection form of nonlinear equations has been applied to the NTU/Purdue nonhydrostatic model (Hsu and Sun 2001 [190], MacCall et al., 2015 [231]). It is simple and can be solved numerically with high accuracy (Sun and Yeh 1997 [124], Sun and Sun 2004 [125]), but it does not conserve the total energy or momentum. Combining those equations, we obtain the flux forms:
ρ t + ρ u j x j = 0
u i ρ t + u i ρ u j x j = p x j + g i + F r i = F i
and
ρ θ t + ρ θ u j x j = 1 c p θ ρ T d q d t + D θ ρ = Q
Without forcing terms on the right-hand side, Equation (131) in the curvilinear coordinate becomes
v ^ k g k t = v ^ j g j g i x ^ i ( v ^ m g m ) = v ^ i x ^ i ( v ^ m g m ) = v ^ i v ^ m g m x ^ i v ^ i v ^ m x ^ i g m or v ^ k t = v ^ i v ^ k x ^ i v ^ i v ^ m g m x ^ i g k .
With terrain given by (239), and the velocity: v = 4 sin ( π x / x L ) sin ( 2 π z / z L ) e x + sin ( 2 π x / x L ) sin ( π z / z L ) e z , where x L = 20 Δx, is the length of the horizontal domain, z L is the height of the domain, we calculate the local change in velocities ( u ^ , v ^ , w ^ ) and ( u ¯ , v ¯ , w ¯ ) in the terrain following coordinates, and compare with (u, v, w) calculated from the Cartesian coordinate. The numerical results of time-derivative of x-component velocities in Figure 68 show that results in ( x ^ 1 , x ^ 2 , x ^ 3 ) coordinate is better than that in ( x ¯ 1 , x ¯ 2 , x ¯ 3 ) coordinate as expected.
The conventional flux form of momentum in the curvilinear is
ρ v ^ j t + 1 J x ^ i ( J ρ v ^ i v ^ j ) g k v ^ i v ^ k g j x ^ i = g j g i p x ^ i +
In the Cartesian coordinate, the component momentum over the entire domain ( u i ρ d v ) in (253) is conserved with Fi = 0. ρ θ d v in (254) is also conserved if Q = 0. But (256) is not conserved if Fi = 0, it is also quite complicated. In the curvilinear systems, we may use the divergence operator to present the flux form of the equations, i.e.,
ψ t = ( ψ v ) = 1 J x ^ j ( J ψ v ^ j ) = 1 J x ^ j ( J v ˜ j )
where v ˜ j = ψ v ^ j , ψ = ρ in (252), ψ = ρ θ in (254), and ψ = ρ u i in (253) where ui is the i-component velocity in the Cartesian coordinate. Equation (257) is much simpler than (256) and can be solved easily. Equation (257) can also be written:
ρ J t = ρ J u ^ i x ^ i
ρ θ J t = ρ θ J u ^ i x ^ i + Q
and
ρ u j J t = ρ u j J u ^ i x ^ i + F j
The total mass is conserved, total energy is conserved for adiabatic flow, and total momentum along each component (uj, in Cartesian coordinate) is also conserved when Fj = 0 in (260). The Navier–Stokes equations including forcing terms can be solved by the finite volume method (Sun 2011) [24]. It is noted that our new coordinate is consistent with Kasahara (1974) [198]. The latter is not in the flux form nor conserves the mass or momentum. The differences between simulations from the WFR based on GS coordinate and NTU/Purdue based on Kasahara are discussed in Section 6.8 and will be discussed more in the Summary.

8. Pollution Model

8.1. Pollution in Convective Boundary Layer (CBL)

Sun and Ogura (1980) [105] included the equation of the turbulent potential temperature–humidity covariance ( θ q ¯ ) to simulate the observed ( θ q ¯ ) in their turbulence parameterization. Sun and Chang (1986a, b) [108,217] and Sun (1988, 1989) [218,219] extended Sun and Ogura’s work by including the temperature–concentration covariance ( c θ ¯ ) to simulate the dispersion of plumes in convective boundary layer (CBL). The 1D equations for wind, potential temperature Θ , and humidity are:
U t = f V f V g u w ¯ z
V t = f U + f U g v w ¯ z
Θ t = θ w ¯ z
Q t = q w ¯ z
The equations of subgrid-turbulences are shown in Section 4 of this article. The model does not include radiation or cloud physics, but they were included in Sun (1993a, b) [199,232] to study the nocturnal precipitation triggered by radiative cooling at stratocumulus cloud top. The equations for the crosswind integration concentration Cy, the concentration flux w c y ¯ and the covariance of concentration and temperature c θ y ¯ are
C y t = U C y x w c y ¯ z
w c y ¯ t = U w c y ¯ x w 2 ¯ C y ¯ z + g Θ ( 1 α 2 ) c θ y ¯ + z 2 E λ w c y ¯ z α 1 2 E λ w c y ¯
c θ y ¯ t = U c θ y ¯ x w θ ¯ C y ¯ z w c y ¯ Θ z + z 2 E λ c θ y ¯ z α 3 2 E λ c θ y ¯
The results shown in Figure 69 are quite comparable with the laboratory experiments shown in Figure 70 (Deardorff and Willis 1975 [233], Willis and Deardorff 1978, 1981 [234,235]), in which the plume released from surface ascends to upper mixed layer, then gradually descends. If the sources are far away from ground, the plumes descend, touch the ground and bounce to the middle and upper mixed layer. The simulated surface concentration is also in good agreement with field experiments (Sun 1988, 1989) [218,219].

8.2. Backward Integration of Diffusion Equation (Sun and Sun 2017) [236]

In geoscience, Lagrangian backward trajectory without mixing is widely used to estimate the property at the source. With mixing, the plume can mix with surrounding and become quite different form the source. Sun and Sun (2017) [236] proposed reverse-in-time integrating to assess the concentrations at the source. The model consists of five different sizes of concentric cylinders with different species. We integrate it forward in time, then use data collected at downwind to integrate the model backward in time. The Leibniz integral rule for density ρ m j for the m-th specie in the j-th cell with velocity V j is:
D D t V j ρ m j d V = V j ρ m j t d V + ρ m j V j · n d A = V j ρ m j S m j d V ,
where S m j is sink, let us define V j = V j ¯ ( cell - sized   mean ) + V j   ( perturbation ) , ρ m j = ρ m j ¯ + ρ m j . It is also assumed · V j ¯ = 0, and sizes of each cells remain constant. (268) becomes:
V j ρ m j ¯ t d V + ρ m j ¯ + ρ m j V j ¯ + V j ¯ · n d A = V j ρ m j ¯ S m j d V , V j D ρ m j ¯ D t d V = D D t V j ρ m j ¯ d V = ρ m j V j ¯ · n d A V j ρ m j ¯ S m j d V .
The bar of cell-sized mean will be dropped hereafter. Following the Lagrangian frame, the change in the density ρ m j for the m-th specie in the j-th cell with volume Vj is
D ρ m j D t = 1 V j ρ m j i j A i , j u i , j + i j ρ m i A i , j u i , j ρ m j V j S m j   for   i = 1 , 2 . , I ;   j = 1 , 2 , . , J ;   and   m = 1 , 2 , M
where A i , j is the surface area between the i-th and j-th cells, and u i , j is the turbulent velocity on A i , j . It is noted that (270) is equivalent to
ρ ¯ t = · ( ρ + ρ ) ( V + V ¯ ) = · ρ ¯ V ¯ + ρ V ¯ = V ¯ · ρ ¯ · ρ V ¯ ;   D ρ ¯ D t = · ρ V ¯
Traditionally, · ρ V ¯ is ignored in the continuity equation for air density, but the subgrid turbulence fluxes for moisture and aerosols are kept in the equations. For simplicity, Vj, A i , j , and u i , j are time independent. The finite difference form of Equation (270) becomes:
ρ m j n + 1 ρ m j n Δ t = 1 V j p ρ m j n + q ρ m j n + 1 i j A i , j u i , j + i j p ρ m i n + q ρ m i n + 1 A i , j u i , j p ρ m j n + q ρ m j n + 1 V j S m j for   j = 1 , 2 , , J ;   i = 1 , 2 , , J ;   and   m = 1 , 2 , , M .
An implicit time step is applied on the right hand, where q = 1 − p, and 0 ≤ p ≤ 1. (272) Becomes:
ρ m j n + 1 1 + q i j A i , j V j u i , j Δ t + S m j q i j A i , j V j u i , j Δ t ρ m i n + 1 = ρ m j n 1 p i j A i , j V j u i , j Δ t + S m j + p i j A i , j V j u i , j Δ t ρ m i n   for   j = 1 , . . , J ;   i = 1 , . , I ;   and   m = 1 , . , M .
We assume that p = 0.5 and S m j = 0. Equation (273) can be integrated forward-in-time from an initial condition or reverse-in-time from that collected at downwind region. The model consists of two layers of cylinders as shown in Figure 71. The depth of the lower layer is 2000 m, which is topped by the cylinder E (light yellow) with radius of 30,000 m and thickness of 2222.2 m. The lower layer includes four concentric cylinders (that will be referred as cylinder or cell hereafter): A (black), B (red), C (green), and D (blue), with radii of RAB (between cylinders A and B) = 5000 m, RBC (between cylinders B and C) = 10,000 m, RCD (between cylinders C and D) = 20,000 m, and RD (outside of D) = 30,000 m. The volume is 157,079 × 106 m3 for cylinder A; 471,239 × 106 m3 for cylinder B; 188,496 × 107 m3 for cylinder C; 314,159 × 107 m3 for cylinder D; and 628,256 × 107 m3 for cylinder E (light blue).
Severe pollutions usually occur in a calm, stable condition. Hence, the value of horizontal u′ is set equal to 2.5 × 10 2   m   s 1 . The vertical mixing u′ is set equal to 0.1 of the horizontal u′ in the lower layer. The time interval Δt = 50 s is applied to both the forward and reverse integrations.

8.2.1. Forward-in-Time Integration

Each cylinder has its unique specie with one unit amount of density initially. The time evolutions of different species in each cylinder are shown in Figure 72a–e. Initially, ρ 1 A (specie 1 in cylinder A, black line) decreases with time exponentially, which is replaced by the material coming from cylinders B ( ρ 2 A specie 2 in cylinder A, red line), C ( ρ 3 A green), D ( ρ 4 A blue), and E ( ρ 5 A light blue). The density of those invaders starts from zero and increases with time. The decrease in ρ 1 A slows down as time increases, at the end of integration, tf = 4 × 105 s (~111.1 h), ρ 1 A ≅ 0.07, while ρ 2 A = 0.144, ρ 3 A = 0.26, ρ 4 A = 0.20, and ρ 5 A = 0.32. Most of original material in cylinder A is replaced by the substances from outside, only 7% of the original material remains at t = 4 × 105 s.
The distribution of ρ 2 B (species 2 in cylinder B), shown in Figure 72b is quickly replaced by the material coming from other cylinders. Only about 11% remains at tf = 4 × 105 s. Figure 72c–e show that the larger the volume, the more the original material can maintain. It is also noted that Cylinder E (Figure 72e) can keep up more than 70% of the original material because of large volume as well as a small u′ between it and the cylinders in the lower layer. Figure 72e also shows that the contributions from the lower cylinders are proportionally to their contact areas with Cylinder E.
The forward-in-time integration (marketed with +) of each specie among different cylinders are indicated by black (in cylinder A), red (B), green (C), blue (D), and light blue (E) lines, and shown in Figure 73a–e.

8.2.2. Reverse-in-Time Integration

Using the density distribution of species at t = 4 × 105 s from forward integration in Figure 73a–e as initial condition, we integrate Equation (273) backward-in-time, i.e., we calculate ρ m j n from a given ρ m j n + 1 by Equation (273). The integrated density from tf = 4 × 105 s to the end at ti = 0 are shown by black, red, green, blue, and light blue lines marked with “ ” signs in Figure 73a–e for each species. They coincide with the forward-in-time integrations (with +). The results end up almost identical to the initial condition in the forward-in-time integration. ρ 5 A (Figure 73) has an error of 1.15% after 2 × 8000 times of forward and backward integrations. Figure 73 shows that ρ 2 A (=0.145) is larger than ρ 2 B (=0.113) at tf = 4 × 105 s, but at the end of the reverse integration ρ 2 A = 0 and ρ 2 B = 1. The values of ρ 3 A , ρ 3 B , and ρ 3 C are around 0.26 at tf = 4 × 105 s, after the reverse integration ρ 3 A = ρ 3 B = 0 and ρ 3 C = 1 at t = 0. It indicates that the model can handle the counter gradient flux during backward integration, a difficult task in diffusion modeling (Sun and Chang 1986a, b, Sun 1988, 1989).
The simulations show the importance of mixing. Because of mixing, the material of each cell at the downwind region is significantly different from its beginning. Hence, the conventional forward or backward trajectory method without mixing can lead to a significant error to interpret the property at downwind regions or sources.
Equation (273) can be written as
x n = A x n 1 = = A n x 0 for t = 0 to tf. in forward integration. Then y n = A 1 y n 1 = = A n y 0 = A n x n in the backward integration starts from t = tf to t = 0 with y 0 = x n as initial condition. The period is 2 tf, and the system after a complete cycle can be determined by matrix: y n = A n x n = A n A n x 0 = I x 0 according to Floquet theorem (Hand and Finch 1998) [17]. The system starts from a simple system that the species mix by perturbation velocity and become quite different from the original situation in the forward integration. After backward integration, the species return to their initial position, which is like the Arnold cat map (Arnold 1967) [237]. It is also noted that in the forward integration, each species approaches a uniform distribution among different cells as time approaches infinity, then the system becomes a fixed point, f(x) = x, no matter forward or backward integration (Strogatz 2015) [238]. Because it is a first-order system, stability is much simpler than the nonlinear dynamics (Poincare, Lyapunnov [239], Lorenz [240], Hand and Finch 1998 [17], Greenwood 2003 [22], and others), but most inverse value problems are difficult to solve.

9. Summary

We discuss the equations of CGFD, numerical methods, applications, and comparisons with observations and earlier studies. In Section 2 of shallow water, in addition to stability analysis of the forward–backward and leapfrog schemes with or without diffusion, we successfully applied Sun’s scheme (2011) [24] without smoothing or nudging to simulate the nonlinear dam break, solitary Rossby wave, hydraulic jump, and vortices merging, etc. Solving the linearized equations (in Section 3), we obtained the coexistence of the longitudinal, shallow cloud bands and transverse, deep cloud bands developed in the tropical region. We also reproduced the observed mesoscale convection system in Brazil due to the resonance of trapezoidal instability and the PBL diurnal oscillation. The coupled atmosphere–ocean Ekman layer shows the resonance instability due to the frequency of the PBL oscillation coinciding with the inertial frequency of the Coriolis force at 30°. It also proves that the weekly or monthly mean of Ekman profiles should not be obtained by solving the steady diffusion equation using different diffusion coefficients for phase angle and amplitude of the velocity.
In Section 4, the numerical schemes and the physical components of the PRCM, except radiation, are presented. The PRCM successfully captured the development of explosive severe winter storms in Western US from a fair weather. It reproduced the formation and advection of lee vortices in Taiwan, the winter storms, and a small observed vortex over the Rocky Mountains, the deformation of a weak cold front by the Central Mountain Range in Taiwan during TAMEX. The results also show that the vertical stretching is most important to the formation of lee vortex. The regional climate simulations in Section 5 show that snow–land module can affect the wind and geopotential at 200 mb. The detailed formulation of snow-soil-vegetation is crucial to the short-term weather forecasting or long-term climate simulations. The PRCM reproduced the 1993 flooding, and 1988 drought and hot dome in the US, it also well simulated the 1998 flooding in China and Korea, as well as the weather and northward movement of East Asian Monsoon. It predicted the onset and decay of the summer Monsoon. The model also replicated the potential vorticity at 850 mb in the South China Sea in the summers of 1991–1998. The PRCM can simulate the transition from a stable state to an unstable state and vice versa. It is also noted the PRCM has an excessive rainfall in the flooding areas. The model seems able to correct the small-scale errors if the synoptic systems are correct. The result confirms that the large-scale systems are more important, which can guard the small-scale disturbances. It implies that the predictability from the NWP model seems longer than the simplified dynamics models in chaos theory.
Despite the chaotic nature of the atmosphere, ECMWF also believes that long-term and seasonal (up to 13 months) forecasting are possible because ENSO shows variations on long time scales (seasons and years) and, to a certain extent, are predictable. The slow variations of snow and soil wetness also affect the long-term forecast. (https://effis.jrc.ec.europa.eu/about-effis/technical-background/seasonal-forecast-explained (accessed on 25 July 2023)). Lorenz (1982) [241] stated “Predictions at least ten days ahead as skillful as predictions now made seven days ahead appear to be possible. Additional improvements at extended range may be realized if the one-day forecast is capable of being improved significantly”.
In Section 6, we discuss Boussinesq fluid versus compressible fluid and different approaches to cut down the computing time in the atmospheric nonhydrostatic models. We also discuss the mistakes of the parcel method, Froude number, and CAPE, because they ignored the dynamic pressure based on the hydrostatic system, in which the flow cannot move against buoyancy; and the KE remains constant if the parcel moves vertically following the hydrostatic–adiabatic lapse rate. On the other hand, the Bernoulli function conserved the summation of KE, PE, and EN based on the nonhydrostatic, compressible atmosphere. According to the Bernoulli function, the plumes can rise against buoyancy and reach 12–18 km as the hot towers of cumulonimbus clouds in the Equator, or the deep convections in hurricanes and severe storms. They can also effectively pump moisture and pollutants deeply into the stratosphere. The simulations also show an upper layer convergence is needed to sustain the downslope windstorm on the lee of mountain on a non-slip surface.
We present a new terrain-following coordinate and the comparisons with the Gal-Chen and Somerville (GS) model in Section 7. GS model introduces errors in divergence, curl, and the advection terms of the equations because they ignored the variation of the (horizontal) spatial interval over the mountain. Using the divergence formula in curvilinear coordinate, we also provide an efficient flux form of the Navier–Stokes equations which conserved the mass, energy, and momentum.
Section 8 shows the model and numerical results from a turbulence–pollution model, which replicates the ascending and descending plumes in a convective boundary layer. The backward integration scheme shows that the data collected from the downwind region can be used to assess the property of pollutants at the source. Because of including mixing, the results are more accurate than those based on the traditional Lagrangian backward trajectory method.
The numerical method is applied to integrate the Navier–Stokes equations from the initial condition (at time step n = 0) to the next time step, n = 1. Then, the results at n = 1 are used to calculate results at n = 2, 3, … The solutions at each time step depend upon not only the previous condition but also the equations and numerical method applied. With the same initial condition, Figure 74a,b show the simulated flows over an elliptic mountain at 10 and 20 h integration from NTU/Purdue and WRF. The vortices at t = 10 h are similar. The vortices developed earlier and moved to the farther downstream are still quite similar at t = 20 h in Figure 74b, but those developed later become quite different near the mountain. It is also noted that the stagnation exists on the windward side of the mountain in the NTU/Purdue simulation, and flow is also wider on the lee compared with the WRF simulation (Sun et al., 2010) [242]. The differences grow with time even for smooth, well resolved cases.
The solutions of the Navier–Stokes equations depend upon the initial and boundary conditions (Keiss and Lorenz 2004 [243]), they gradually deteriorate with time and eventually lose the memory of the initial conditions even without any strange attractor. Beyond that, the NWP cannot be applied (a predictability limit?). On the other hand, the numerical results are continuously influenced by the forcing and lateral boundary conditions, examples include the seasonable variations, land-sea contrast, terrain effect and ground/sea surface, etc. Advance in supercomputer does help the CGFD and NWP, and scientists will continue to face the challenge of formulating the physics, sub-grid fluxes, and handing the boundaries. Sometimes, a model may not be able to predict a well-resolved large system either, for example, it is difficulty to accurately replicate the strength and extent of the Western Pacific subtropical high, which affects the regional flooding in Southeast Asia (Hsu et al., 2004 [164], Sun et al., 2011 [165], Sun et al. [192]).

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The author is indebted to H. L. Kuo and other teachers for their teaching and stimulating. He also appreciates his family for their support, and the opportunities working with the former students and colleagues, who have contributed significantly to the research and model development. He is grateful to N. H. Lin, H. H. Hsu, K. Tsuboki, B. Shen, Mandić and the reviewers for their helps, comments, and suggestions. He also thanks Purdue University and P. L. Lin at National Central University for providing computing facility. The author withdrew Sun 2023 [200] at the final stage of the gallery proof, because OJFD broke the agreement and refused including Taiwan in author’s affiliation.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Bjerknes, V. Das Problem der Wettervorhersage, Betrachtet vom Stanpunkt der Mechanik und der Physik. Meteor. Zeits 1904, 21, 1–7. [Google Scholar]
  2. Richardson, L.F. Weather Prediction by Numerical Processes; Cambridge University Press: Boston, MA, USA, 1922; p. 66. ISBN 9780511618291. [Google Scholar]
  3. Charney, J.; Fjørtoft, R.; von Neumann, J. Numerical Integration of the Barotropic Vorticity Equation. Tellus 1950, 2, 237–254. [Google Scholar] [CrossRef]
  4. Phillips, N.A. The general circulation of the atmosphere: A numerical experiment. Q. J. R. Meteorol. Soc. 1956, 82, 123–154. [Google Scholar] [CrossRef]
  5. Poincaré, H. Sur les Fonctions Fuchsiennes. Comptes Rendus Hebd. L’académie Sci. Paris 1882, 94, 1166–1167. [Google Scholar]
  6. Poincaré, H. Sur le problème des trois corps et les équations de la dynamique. Acta Math. 1890, 13, 1–270. (In French) [Google Scholar]
  7. Poincaré, H. Solutions Periodiques, Non-Existence des Integrales Uniformes, Solutions Asymptotiques; Gauthier-Villars: Paris, France, 1892; Volume 1. (In French) [Google Scholar]
  8. Poincare, H. Memoire sur Jes courbes definies par une equation differentielle. J. Math. 1881, 7, 375–442. [Google Scholar]
  9. Lorenz, E. Deterministic nonperiodic flow. J. Atmos. Sci. 1963, 20, 130–141. [Google Scholar] [CrossRef]
  10. Lorenz, E. The predictability of a flow which possess many scales of motion. Tellus 1969, 21, 289–307. [Google Scholar] [CrossRef]
  11. Smagorinsky, J. Problems and promises of deterministic extended range forecasting. Bull. Am. Meteorol. Soc. 1969, 50, 286–311. [Google Scholar] [CrossRef]
  12. Arakawa, A. Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow. Part I. J. Comput. Phys. 1966, 1, 119–143. [Google Scholar] [CrossRef]
  13. Paegle, P.; Yang, Q.; Wang, M. Predictability in limited area and global models. Meteorol. Atmos. Phys. 1997, 63, 53–69. [Google Scholar] [CrossRef]
  14. Wiin-Nielsen, A. On the stable part of the Lorenz attractor. Atmôsfera 1998, 11, 61–73. [Google Scholar]
  15. Charney, J.G.; Shukla, J. Predictability of monsoons. In Monsoon Dynamics; Lighthill, S.J., Pearce, R.P., Eds.; Cambridge University Press: Cambridge, UK, 1981; pp. 99–109. [Google Scholar]
  16. Hoskins, B.J. Predictability Beyond the Deterministic Limit. WMO Bull. 2012, 61, 33–36. [Google Scholar]
  17. Hand, L.N.; Finch, J.D. Analytical Mechanics; Cambridge University Press: Cambridge, UK, 2008; 575p. [Google Scholar]
  18. Beltrami, E. Mathematic for Dynamic Modeling; Academic Press: Cambridge, MA, USA, 1987; 277p. [Google Scholar]
  19. Shen, B.W. Homoclinic Orbits and Solitary Waves within the Nondissipative Lorenz Model and KdV Equation. Int. J. Bifurc. Chaos 2020, 30, 15. [Google Scholar] [CrossRef]
  20. Shen, B.-W.; Pielke, S.; Pielke, R.; Cui, J.; Faghih-Naini, S.; Paxson, W.; Kesarkar, A.; Atlas, R. The dual nature of chaos and order in the atmosphere. Atmosphere 2022, 13, 1892. [Google Scholar] [CrossRef]
  21. Shen, B.-W.; Pielke, R.A.; Zeng, X.; Cui, J.; Faghih-Naini, S.; Paxson, W.; Atlas, R. Three kinds of butterfly effects within Lorenz Models. Encyclopedia 2022, 2, 1250–1259. [Google Scholar] [CrossRef]
  22. Greenwood, D.T. Advanced Dynamics; Cambridge Press: New York, NY, USA, 2003; 425p. [Google Scholar]
  23. Sun, W.Y. Instability in leapfrog and forward-backward Schemes. Mon. Wea. Rev. 2010, 138, 1497–1501. [Google Scholar] [CrossRef]
  24. Sun, W.Y. Instability in Leapfrog and Forward-Backward Schemes: Part II: Numerical Simulation of Dam Break. J. Comput. Fluids 2011, 45, 70–76. [Google Scholar] [CrossRef]
  25. Mesinger, F.; Arakawa, A. Numerical Methods Used in Atmospheric Models. GARP Glob. Atmos. Res. Program 1976, 17, 64. [Google Scholar]
  26. Gill, A.E. Atmospheric-Ocean Dynamics; Academic Press: San Diego, CA, USA, 1982; 661p. [Google Scholar]
  27. Haltiner, G.J.; Williams, R.T. Numerical Prediction and Dynamic Meteorology, 2nd ed.; Wiley: Hoboken, NJ, USA, 1980; 477p. [Google Scholar]
  28. Pielke, S.R.A. Mesoscale Meteorological Modeling, 2nd ed.; Academic Press: Cambridge, MA, USA, 2002; 676p. [Google Scholar]
  29. Lin, Y.L. Mesoscale Dynamics; Cambridge University Press: New York, NY, USA, 2007; 633p. [Google Scholar]
  30. Holton, J.R.; Hakim, G.J. An Introduction to Dynamic Meteorology, 5th ed.; Academic Press: Cambridge, MA, USA, 2012; 531p. [Google Scholar]
  31. Versteeg, H.K.; Malalasekera, W. An Introduction to Computational Fluid Dynamics: The Finite Volume Method; Longman Scientific & Technical: Harlow, UK, 1995; 257p. [Google Scholar]
  32. Pielke, R.A. Mesoscale Metorological Modeling; Elsevier Inc.: Amsterdam, The Netherlands, 1984; 612p. [Google Scholar]
  33. Pedlosky, J. Geophysical Fluid Dynamics; Springer: Berlin/Heidelberg, Germany, 1979; 626p. [Google Scholar]
  34. Lorenz, E.N. Attractor sets and quasi-geostrophic equilibrium. J. Atmos. Sci. 1980, 37, 1685–1699. [Google Scholar] [CrossRef]
  35. Nobel, B. Applied Linear Algebra; Prentice-Hall: Upper Saddle River, NJ, USA, 1969; 523p. [Google Scholar]
  36. Sun, W.Y. A comparison of two explicit time integration schemes applied to the transient heat equation. Mon. Wea. Rev. 1982, 110, 1645–1652. [Google Scholar] [CrossRef]
  37. Sun, W.Y.; Sun, O.M.T. A modified leapfrog scheme for shallow water equations. Comput. Fluids 2011, 52, 69–72. [Google Scholar] [CrossRef]
  38. Sun, W.Y. Numerical analysis for hydrostatic and nonhydrostatic equations of inertial-internal gravity waves. Mon. Wea. Rev. 1984, 112, 259–268. [Google Scholar] [CrossRef]
  39. Roe, P.L. Approximate Riemann Solvers, Parameter Vectors and Difference Schemes. J. Compu. Phys. 1981, 43, 357–372. [Google Scholar] [CrossRef]
  40. Toro, E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics; Springer: Berlin/Heidelberg, Germany, 1999; 309p. [Google Scholar]
  41. Wang, H.; Yeh, G.T. A characteristic-based semi-Lagrangian method for hyperbolic systems of conservation laws. Chin. J. Atmos. Sci. 2005, 29, 21–42. [Google Scholar]
  42. Oh, T.J. Development and Testing of Characteristic-Based Semi-Lagrangian Two-Dimensional Shallow Water Equations Development and Testing of Characteristic-Based Semi-Lagrangian Two-Dimensional Shallow Water Equations Model. Ph.D. Thesis, Purdue University, W. Lafayette, IN, USA, 2007; 148p. Available online: https://docs.lib.purdue.edu/dissertations/AAI3278686/ (accessed on 1 January 2008).
  43. MacCormack, R.W. The Effect of Viscosity in Hyper-Velocity Scattering; AIAA: Washington, DC, USA, 1969; pp. 69–354. [Google Scholar]
  44. Fujihara, M.; Borthwick, A.G.L. Godunov-type solution of curvilinear shallow-water equations. J. Hydraul. Eng. ASCE 2000, 126, 827–836. [Google Scholar] [CrossRef]
  45. Anastasiou, K.; Chan, C.T. Solution of the 2D shallow water equations using the finite volume method on unstructured triangular meshes. Int. J. Numer. Methods Fluids 1997, 24, 1225–1245. [Google Scholar] [CrossRef]
  46. Tseng, M.H.; Chu, C.R. Two-dimensional shallow water flows simulation using TVD-MacCormack scheme. J. Hydraul. Res. 2000, 38, 123–131. [Google Scholar]
  47. Liang, Q.; Borthwick, A.G.L.; Stelling, G. Simulation of dam- and dyke-break hydrodynamics on dynamically adaptive quadtree grids. Int. J. Numer. Methods Fluids 2004, 46, 127–162. [Google Scholar] [CrossRef]
  48. Sun, W.Y.; Sun, O.M.T. Numerical Simulation of Rossby Wave in Shallow Water. Comput. Fluids 2013, 76, 116–127. [Google Scholar] [CrossRef]
  49. Chu, P.C.; Fan, C. Space-Time Transformation in Flux-form Semi-Lagrangian Schemes. Terr. Atmos. Ocean. Sci. 2010, 21, 17–26. [Google Scholar] [CrossRef]
  50. Boyd, J.P. Equatorial solitary waves. Part-1: Rossby solitons. J. Phys. Oceanogr. 1980, 10, 1699–1717. [Google Scholar] [CrossRef]
  51. Boyd, J.P. Equatorial solitary waves. Part 3: Westward-traveling Motions. J. Phys. Oceanogr. 1985, 15, 46–54. [Google Scholar] [CrossRef]
  52. Korteweg, D.J.; de Vries, G. On the change of long waves advancing in a rectangular canal, and on a new type of long stationary waves. Phil. Mag. 1895, 39, 422–433. [Google Scholar] [CrossRef]
  53. Houghton, D.D.; Kasahara, A. Nonlinear shallow fluid flow over an isolated ridge. Commun. Pure Appl. Math 1968, 11, 1–23. [Google Scholar] [CrossRef]
  54. Sun, W.Y. Hydraulic jump and Bernoulli equation in nonlinear shallow water model. Dyn. Atmos. Ocean. 2018, 82, 37–53. [Google Scholar] [CrossRef]
  55. Sun, W.Y. The Vortex Moving toward Taiwan and the influence of the Central Mountain Range. Geosci. Lett. 2016, 3, 21. [Google Scholar] [CrossRef]
  56. Lin, Y.-L.; Chen, S.-Y.; Hill, C.M.; Huang, C.-Y. Control parameters for tropical cyclones passing over mesoscale mountains. J. Atmos. Sci. 2005, 62, 1849–1866. [Google Scholar] [CrossRef]
  57. Huang, C.-Y.; Lin, Y.-L. The influence of mesoscale mountains on vortex tracks: Shallow-water modeling study. Meteor. Atmos. Phys. 2008, 101, 1–20. [Google Scholar] [CrossRef]
  58. Mastsuno, T. Quasi-geosprophic motions in the equator area. J. Meteor. Soc. Jpn. 1966, 44, 25–43. [Google Scholar] [CrossRef]
  59. Sun, W.Y.; Chern, J.D. Numerical experiments of vortices in the wake of idealized large mountains. J. Atmos. Sci. 1994, 51, 191–209. [Google Scholar] [CrossRef]
  60. Sun, W.Y.; Sun, O.M.T. Bernoulli equation and flow over a mountain. Geosci. Lett. 2015, 2, 7. [Google Scholar] [CrossRef]
  61. Wu, C.C.; Li, T.-H.; Huang, Y.-H. Influence of mesoscale topography on tropical cyclone tracks: Further examination of the channeling effect. J. Atmos. Sci. 2015, 72, 3032–3050. [Google Scholar] [CrossRef]
  62. Sun, W.Y.; Oh, T.J. Vortex Merger in Shallow Water Model. Asia-Pac. J. Atmos. Sci. 2022, 58, 533–547. [Google Scholar] [CrossRef]
  63. Waugh, D.W. The efficiency of symmetric vortex merger. Phys. Fluid Dyn. 1992, 4, 1745–1758. [Google Scholar] [CrossRef]
  64. Kimura, Y.; Herring, J.R. Gradient enhancement and filament ejection for a non-uniform elliptic vortex in two-dimensional turbulence. J. Fluid Mech. 2001, 439, 43–56. [Google Scholar] [CrossRef]
  65. Cerretelli, C.; Williamson, C.H.K. The physical mechanism for vortex merging. J. Fluid Mech. 2003, 475, 41–77. [Google Scholar] [CrossRef]
  66. Melander, M.V.; McWilliams, J.C.; Zabusky, N.J. Axisymmetrization and vorticity-gradient intensification of an isolated two-dimensional vortex through filamentation. J. Fluid Mech. 1987, 178, 137–159. [Google Scholar] [CrossRef]
  67. Melander, M.V.; Zabusky, N.J.; McWilliams, J.C. Symmetric vortex merger in two dimensions: Causes and conditions. J. Fluid Mech. 1988, 195, 305–340. [Google Scholar] [CrossRef]
  68. Huang, M.J. The physical mechanism of symmetric vortex merger: A new viewpoint. Phys. Fluids 2005, 17, 1–7. [Google Scholar] [CrossRef]
  69. Moore, D.W.; Saffman, P.G. The instability of a straight vortex filament in a strain field. Proc. R. Soc. Lond. A 1975, 346, 413–425. [Google Scholar]
  70. Brandt, L.K.; Nomura, K.K. The physics of vortex merger: Further insight. Phys. Fluids 2006, 18, 051701. [Google Scholar] [CrossRef]
  71. Sun, W.Y.; Chern, J.D. Numerical study of influence of mountain ranges on Mei-yu front. J. Meteorol. Soc. Jpn. 2006, 84, 27–46. [Google Scholar] [CrossRef]
  72. Fjortoft, R. On the changes in the spectral distribution of kinetic energy for two dimensional nondivergent flow. Tellus 1953, 5, 225–230. [Google Scholar] [CrossRef]
  73. Kraichnan, R.H. Tnertial renges in two-dimensional turbulence. Phys. Fluid 1967, 10, 1417–1423. [Google Scholar] [CrossRef]
  74. Kolmogorov, A.N. Local structure of turbulence in an incompressible fluid at very high Reynolds numbers. Dokl. Akad. Nauk. SSSR 1941, 31, 99–101. [Google Scholar] [CrossRef]
  75. Obukhov, A.M. Spectral energy distribution in a turbulent flow. Dokl. Akad. Nauk. SSSR 1941, 32, 22–24. [Google Scholar]
  76. Eady, E.T. Long wave and cyclone waves. Tellus 1949, 1, 33–52. [Google Scholar] [CrossRef]
  77. Sun, W.Y. Stability of deep cloud streets. J. Atmos. Sci. 1978, 35, 466–483. [Google Scholar] [CrossRef]
  78. Malkus, J.S.; Riehl, H. Cloud Structure and Distribuions over the Tropical Pacific Ocean; University of California Press: Berkeley, CA, USA; Los Angeles, CA, USA, 1964; 229p. [Google Scholar]
  79. Kuo, H.L. Perturbations of plane Couette flow in a stratified fluid and origin of cloud streets. Phys. Fluid 1963, 6, 195–211. [Google Scholar] [CrossRef]
  80. Sun, W.Y. Rainbands and symmetric instability. J. Atmos. Sci. 1984, 41, 3412–3426. [Google Scholar] [CrossRef]
  81. Fjortoft, R. Application of integral theorems in deriving criteria of stability for laminar flows and for the baroclinic circular vortex. Geofys. Pub. 1950, 17, 1–52. [Google Scholar]
  82. Ooyama, K. On the stability of baroclinic circular vortex: A sufficient criterion for instability. J. Atmos. Sci. 1966, 23, 43–53. [Google Scholar] [CrossRef]
  83. Hoskins, B.J. Baroclinic Instability and Frontogensis in Rotating Fluids in Geophysics; Roberts, P.H., Ed.; Academic Press: Cambridge, MA, USA, 1978; pp. 171–204. [Google Scholar]
  84. Sun, W.Y. Convection in the Surface Layer of the Atmosphere. Ph.D. Thesis, The University of Chicago, Chicago, IL, USA, 1975; 130p. [Google Scholar]
  85. Kuo, H.L.; Sun, W.Y. Convection in the lower atmosphere and its effects. J. Atmos. Sci. 1976, 33, 21–40. [Google Scholar] [CrossRef]
  86. Agee, E.; Church, C.; Morris, C.; Snow, J. Some synoptic aspects and dynamic features of vortices associated with the tornado outbreak of 3 April 1974. Mon. Wea. Rev. 1975, 103, 318–353. [Google Scholar] [CrossRef]
  87. Sun, W.Y.; Orlanski, I. Large mesoscale convection and sea breeze circulation. Part I: Linear stability analysis. J. Atmos. Sci. 1981, 38, 1675–1693. [Google Scholar] [CrossRef]
  88. Sun, W.Y.; Orlanski, I. Large mesoscale convection and sea breeze circulation. Part 2: Non-Linear numerical model. J. Atmos. Sci. 1981, 38, 1694–1706. [Google Scholar] [CrossRef]
  89. Orlanski, I.; Polinsky, L.J. Spectral distribution of cloud cover over Africa. J. Meteor. Soc. 1977, 55, 483–494. [Google Scholar] [CrossRef]
  90. Sun, W.Y. Unsymmetrical symmetric instability. Q. J. R. Meteorol. Soc. 1995, 121, 419–431. [Google Scholar] [CrossRef]
  91. Sun, W.Y.; Sun, O.M. Inertia and diurnal oscillations of Ekman layers in atmosphere and ocean. Dyn. Atmos. Ocean. 2020, 90, 289. [Google Scholar] [CrossRef]
  92. Blackadar, A.K. Boundary layer wind maxima and their significance for the growth of nocturnal inversion. Bull. Am. Meteorol. Soc. 1957, 38, 283–290. [Google Scholar] [CrossRef]
  93. Palmer, E.; Newton, C.W. Atmospheric Circulation Systems; Academic Press: London, UK, 1969; 606p. [Google Scholar]
  94. Price, J.F.; Weller, R.A.; Pinkel, R. Diurnal cycling: Observations and models of the upper ocean response to heating, cooling, and wind mixing. J. Geophys. Res. 1986, 91, 8411–8427. [Google Scholar] [CrossRef]
  95. Kalashnik, M.V. Resonant excitation of baroclinic waves in the presence of Ekman friction. Izvestiya. Atmos. Ocean. Phys. 2018, 54, 109–113. [Google Scholar] [CrossRef]
  96. Brown, R.A. A secondary flow model for the planetary boundary layer. J. Atmos. Sci. 1970, 27, 742–757. [Google Scholar] [CrossRef]
  97. Polton, J.A.; Lenn, Y.D.; Elipot, S.; Chereskin, T.K.; Sprintall, J. Can Drake Passage observations match Ekman’s classic theory? J. Phys. Oceanogr. 2013, 43, 1733–1740. [Google Scholar] [CrossRef]
  98. Chereskin, T.K. Direct evidence for an Ekman balance in the California Current. J. Geo. Res. 1995, 100, 261–269. [Google Scholar] [CrossRef]
  99. Price, J.F.; Sundermeyer, M.A. Stratified Ekman layers. J. Geoph. Res. 1999, 104, 20467–20494. [Google Scholar] [CrossRef]
  100. Lenn, Y.-D.; Chereskin, T.K. Observations of Ekman currents in the Southern Ocean. J. Phys. Ocean. 2009, 39, 768–799. [Google Scholar] [CrossRef]
  101. Fujita, T.T. Structure and movement of a dry front. Bull. Am. Meteor. Soc. 1958, 39, 574–582. [Google Scholar] [CrossRef]
  102. Wexler, H. A boundary layer interpretation of low-level jet. Tellus 1961, 12, 268–378. [Google Scholar]
  103. Sun, W.Y.; Ogura, Y. Boundary layer forcing as a possible trigger to a squall-line formation. J. Atmos. Sci. 1979, 36, 235–254. [Google Scholar] [CrossRef]
  104. Sun, W.Y.; Wu, C.C. Formation and diurnal oscillation of dryline. J. Atmos. Sci. 1992, 49, 1606–1619. [Google Scholar] [CrossRef]
  105. Sun, W.Y.; Ogura, Y. Modeling the evolution of the convective planetary boundary layer. J. Atmos. Sci. 1980, 37, 1558–1572. [Google Scholar] [CrossRef]
  106. Deardorff, J.W. Stratocumulus-capped mixed layers derived from a three-dimensional model. Bound-Layer Meteor. 1980, 18, 495–527. [Google Scholar] [CrossRef]
  107. Caughey, S.J.; Palmer, S.G. Some aspects of turbulence structure through the depth of the convective layer. Quart. J. Roy. Meteor. Soc. 1979, 105, 811–827. [Google Scholar] [CrossRef]
  108. Sun, W.Y.; Chang, C.Z. Diffusion model for a convective layer: Part 1. Numerical simulation of a convective boundary layer. J. Clim. Appl. Meteor. 1986, 25, 1445–1453. [Google Scholar] [CrossRef]
  109. Sun, W.Y. Air pollution in a convective boundary. Atmos. Environ. 1986, 20, 1877–1886. [Google Scholar] [CrossRef]
  110. Chern, J.D. Numerical Simulation of Cyclogenesis over the Western United States. Ph.D. Thesis, Purdue University, W. Lafayette, IN, USA, 1994; 178p. [Google Scholar]
  111. Chen, S.H.; Sun, W.Y. A one-dimensional time dependent cloud model. J. Meteor. Soc. Jpn. 2002, 80, 99–118. [Google Scholar] [CrossRef]
  112. Bosilovich, M.G.; Sun, W.Y. Formation and Verification of a land surface parameterization for atmospheric models. Bound.-Layer Meteorol. 1995, 73, 321–341. [Google Scholar] [CrossRef]
  113. Chern, J.D.; Sun, W.Y. Formulation and validation of a snow model. In Proceedings of the A11A-06, 1998 Fall Meeting, AGU, Supplement to EOS, AGU, 79, 45, F96, Washingto, DC, USA, 10 November 1998. [Google Scholar]
  114. Sun, W.Y.; Chern, J.D. One dimensional snow-land surface model applied to Sleepers experiment. Bound.-Layer Meteorol. 2005, 116, 95–115. [Google Scholar] [CrossRef]
  115. Chou, M.-D.; Suarez, M.J. A solar radiation parameterization for atmospheric studies. NASA Tech. Memo 1994, 15, 40. [Google Scholar]
  116. Chou, M.-D.; Suarez, M.J.; Liang, X.Z.; Yan, M.H.H. A thermal infrared radiation parameterization for atmospheric sciences. NASA Tech. Memo 2001, 19, 68. [Google Scholar]
  117. Kuo, H.L. On formation and intensification of tropical cyclones through latent heat release by cumulus convection. J. Atmos. Sci. 1965, 22, 40–63. [Google Scholar] [CrossRef]
  118. Kuo, H.L. Further studies of the parameterization of the influence of cumulus convection on large-scale flow. J. Atmos. Sci. 1974, 31, 1232–1240. [Google Scholar] [CrossRef]
  119. Anthes, R.A. A cumulus parameterization scheme utilizing a one-dimensional cloud model. Mon. Wea. Rev. 1977, 105, 270–286. [Google Scholar] [CrossRef]
  120. Molinari, J. A method for calculating the effects of deep cumulus convection in numerical models. Mon. Wea. Rev. 1982, 110, 1527–1534. [Google Scholar] [CrossRef]
  121. Sun, W.Y. A forward-backward time integration scheme to treat internal gravity waves. Mon. Wea. Rev. 1980, 108, 402–407. [Google Scholar] [CrossRef]
  122. Sun, W.Y. Numerical experiments for advection equation. J. Comput. Phys. 1993, 108, 264–271. [Google Scholar] [CrossRef]
  123. Sun, W.Y.; Yeh, K.S.; Sun, R.Y. A simple Semi-Lagrangian Scheme for advection equation. Q. J. R. Meteorol. Soc. 1996, 122, 1211–1226. [Google Scholar] [CrossRef]
  124. Sun, W.Y.; Yeh, K.S. A general semi-Lagrangian advection scheme employing forward trajectories. Quart. J. Roy. Meteor. Soc. 1997, 123, 2463–2476. [Google Scholar] [CrossRef]
  125. Sun, W.Y.; Sun, O.M.T. Mass Correction Applied to Semi-Lagrangian Advection Scheme. Mon. Wea. Rev. 2004, 132, 975–984. [Google Scholar] [CrossRef]
  126. Sun, W.Y. Pressure gradient in a sigma coordinate. Terr. Atmos. Ocean. Sci. 1995, 4, 579–590. [Google Scholar] [CrossRef] [PubMed]
  127. Yang, K.J.-S. Numerical Simulations of Asian Aerosols and Their Regional Climatic Impacts: Dust Storms and Bio-Mass Burning. Master’s Thesis, Purdue University, W. Lafayette, IN, USA, 2004; 141p. [Google Scholar]
  128. Yang, K.J.-S. Regional Climate-Chemistry Model Simulations of Ozone in the Lower Troposphere and Its Climatic Impacts. Ph.D. Thesis, Purdue University, W. Lafayette, IN, USA, 2004; 157p. [Google Scholar]
  129. Schaefer, J.T. The life cycle of the dryline. J. Appl. Meteorol. Climatol. 1974, 13, 444–449. [Google Scholar] [CrossRef]
  130. Sun, W.Y.; Chern, J.D. Diurnal variation of lee-vortexes in Taiwan and surrounding area. J. Atmos. Sci. 1993, 50, 3404–3430. [Google Scholar] [CrossRef]
  131. Sun, W.Y.; Chern, J.D.; Wu, C.C.; Hsu, W.R. Numerical simulation of mesoscale circulation in Taiwan and surrounding area. Mon. Wea. Rev. 1991, 119, 2558–2573. [Google Scholar] [CrossRef]
  132. Kuo, Y.; Chen, G. The Taiwan Area Mesoscale Experiment (TAMEX): An overview. Bull. Am. Meteor. Sci. 1990, 71, 488–503. [Google Scholar] [CrossRef]
  133. Hsu, W.R.; Sun, W.Y. A numerical study of a low-level jet and its accompanying secondary circulation in a Mei-Yu system. Mon. Wea. Rev. 1994, 122, 324–340. [Google Scholar] [CrossRef]
  134. Haines, P.A.; Chern, J.D.; Sun, W.Y. Numerical simulation of the valentine’s Day storm during WISP 1990. Tellus 1997, 49, 595–612. [Google Scholar] [CrossRef]
  135. Yildirim, A. Numerical Study of an Idealized Cyclone Evolution and Its Subsynoptic Features. Ph.D. Thesis, Purdue University, W. Lafayette, IN, USA, 1994; 145p. [Google Scholar]
  136. Bosilovich, M.G.; Sun, W.-Y. Numerical simulation of the 1993 Midwestern Flood: Local and remote sources of water. J. Geophys. Res. 1999, 104, 415–423. [Google Scholar] [CrossRef]
  137. Bosilovich, M.G.; Sun, W.-Y. Numerical simulation of the 1993 Midwestern flood: Land-Atmosphere Interactions. J. Clim. 1999, 12, 1490–1505. [Google Scholar] [CrossRef]
  138. Oglesby, R.J.; Erikson, D. Soil moisture and persistence of North American drought. J. Clim. 1989, 2, 1362–1380. [Google Scholar] [CrossRef]
  139. Oglesby, R.J.; Leap, D.I.; Sun, W.Y.; Agee, E.M. Modeling the effects of climatic changes on availability and quality of water. Comput. Mech. 1994, 2, 9–17. [Google Scholar]
  140. Sun, W.Y.; Bosilovich, M.G.; Chern, J.D. CCM1 regional response to anomalous surface properties. Terr. Atmos. Ocean. Sci. 1997, 8, 271–288. [Google Scholar] [CrossRef]
  141. Sun, W.Y.; Chern, J.D.; Bosilovich, M. Numerical Study of the 1988 Drought in the United States. J. Meteorol. Soc. Jpn. 2004, 82, 1667–1678. [Google Scholar] [CrossRef]
  142. Bosilovich, M.B.; Sun, W.Y. Simulation of soil moisture and temperature. J. Atmos. Sci. 1998, 55, 1170–1183. [Google Scholar] [CrossRef]
  143. Sun, W.Y.; Bosilovich, M.G. Planetary boundary layer and surface layer sensitivity to land surface parameters. Bound.-Layer Meteorol. 1996, 77, 353–378. [Google Scholar] [CrossRef]
  144. Sun, W.Y. Chapter 6. Interactions among atmosphere, soil, and vegetation. In Applications and Development of Agrometeorology; Yang, L., Ed.; Taiwan Agricultural Research Institute and Chinese Society of Agrometeorology: Taipei, Taiwan, 2001; pp. 69–92. [Google Scholar]
  145. Jordan, R. A One-Dimensional Temperature Model for a Snow Cover; Special Report 91-16; U.S. Army Corps of Engineers, Cold Regions Research and Engineering Laboratory: Hanover, NH, USA, 1991; 49p. [Google Scholar]
  146. Johansen, O. Thermal Conductivity of Soils. Ph.D. Thesis, Trondheim University, Trondheim, Norway, 1975. [Google Scholar]
  147. Farouki, O. Evaluation of Methods for Calculating Soil Thermal Conductivity; CRREL Rep. 82-8; CRREL: Hanover, NH, USA, 1982; 90p. [Google Scholar]
  148. Clapp, R.B.; Hornberg, G.M. Empirical equations for some soil hydraulic properties. Water Resour. Res. 1978, 13, 601–604. [Google Scholar] [CrossRef]
  149. Dingman, S.L. Physical Hydrology, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, USA, 2002; 646p. [Google Scholar]
  150. Verseghy, D.L. CLASS- A Canadian land surface scheme for GCMS. I: Soil model. Int. J. Climatol. 1991, 11, 111–133. [Google Scholar] [CrossRef]
  151. Oke, T.R. Boundary Layer Climates; Routledge: London, UK, 1988; 464p. [Google Scholar]
  152. Bohren, C.; Barkstrom, B. Theory of the optical properties of snow. J. Geophys. Res. 1974, 79, 4527–4535. [Google Scholar] [CrossRef]
  153. Anderson, E.A.; Greenan, H.J.; Whipkey, R.Z.; Machell, C.T. 1977: NOAA-ARS cooperative snow research project—Watershed hydro-climatology and data for water year 1960–1974. NOAA Tech. Rep. 1977, 316. Available online: https://onlinebooks.library.upenn.edu/webbin/book/lookupid?key=olbp53867 (accessed on 25 July 2023).
  154. Yen, Y.C. Heat Transfer Characteristics of Naturally Compacted Snow; Research Report 166; U.S. Army Cold Regions Research and Engineering Laboratory: Hanover, NH, USA, 1965; 9p. [Google Scholar]
  155. Anderson, E.A. A Point Energy Balance Model for a Snow Cover; NOAA Tech. Rep. NWS 19; Office of Hydrology, National Weather Service: Chanhassen, MN, USA, 1976.
  156. LaChapelle, E.R. Field Guide to Snow Crystals; University of Washington Press: Seattle, WA, USA; London, UK, 1969; 101p. [Google Scholar]
  157. Lynch-Stieglitz, M. The development and validation of a simple snow model for GISS GCM. J. Clim. 1994, 7, 1842–1855. [Google Scholar] [CrossRef]
  158. Pitman, A.J.; Yang, Z.-L.; Cogley, J.G.; Henderson-Sellers, A. Description of Bare Essentials of Surface Transfer for the Bureau of Meteorological Research Centre AGCM; BMRC Research Report No. 32; BRMC: Ballarat, Australia, 1991; 117p. [Google Scholar]
  159. Sun, W.Y.; Chern, J.D. One-Dimensional Sea Ice-Ocean Model Applied to SHEBA Experiment in 1997–1998 Winter: Terr. Atmos. Ocean. Sci. 2010, 21, 1–15. [Google Scholar] [CrossRef]
  160. Min, K.H. Coupled Cryosphere Model Development for Regional Climate Study and the Initialization of Purdue Regional Model with the Land Data Assimilation System. Ph.D. Thesis, Purdue University, W. Lafayette, IN, USA, 2005; 149p. [Google Scholar]
  161. Min, K.-H.; Sun, W.Y. Atmosphere-Cryosphere Coupled Model Development and Its Application for Regional Climate Studies. Adv. Meteorol. 2015, 2015, 764970. [Google Scholar] [CrossRef]
  162. Yu, Y.C.; Hsu, H.H.; Kau, W.S.; Tsou, T.H.; Hsu, W.R.; Sun, W.Y. An evaluation of the East Asian summer monsoon simulation by the Purdue Regional Model. J. Environ. Prot. Soc. Taiwan 2004, 27, 41–56. [Google Scholar]
  163. Yu, Y.C.; Hsu, H.-H.; Kau, W.S.; Tsou, T.H.; Hsu, W.R.; Sun, W.Y.; Yu, W.N. An evaluation of the east Asian summer monsoon simulation using the Purdue Regional Model. In Proceedings of the Third Workshop on Regional Climate Modeling for Monsoon System, Honolulu, HI, USA, 17 February 2004. 23p. [Google Scholar]
  164. Hsu, W.-R.; Hou, J.P.; Wu, C.C.; Sun, W.Y.; Tcheng, S.C.; Chang, H.Y. Large-eddy simulation of cloud streets over the East China Sea during cold-air outbreak events, pm.2.4. In Proceedings of the 16th Symposium on Boundary Layers and Turbulence, Portland, ME, USA, 9–13 August 2004; pp. 1–7. [Google Scholar]
  165. Sun, W.Y.; Min, K.H.; Chern, J.-D. Numerical study of 1998 flooding in East Asia. Asia-Pac. J. Atmos. Sci. 2011, 47, 123–135. [Google Scholar] [CrossRef]
  166. Molinari, J.; Dudek, M. Parameterization of Convective Precipitation in Mesoscale Numerical Models: A Critical Review. Mon. Wea. Rev. 1992, 120, 326–344. [Google Scholar] [CrossRef]
  167. Shen, B.W. Aggregated negative feedback in a generalized Lorenz model. Int. J. Bifurc. Chaos 2019, 29, 1950037–1950091. [Google Scholar] [CrossRef]
  168. Ginoux, P.; Chin, M.; Tegen, I.; Prospero, J.M.; Holben, B.; Dubovik, O.; Lin, S.-J. Sources and distributions of dust aerosols simulated with the GOCART model. J. Geophys. Res. 2001, 106, 255–273. [Google Scholar] [CrossRef]
  169. DeFries, R.S.; Townshend, J.R.G. NDVI-derived land cover classification a global scale. Int. J. Remote Sens. 1994, 15, 3567–3586. [Google Scholar] [CrossRef]
  170. Tegen, I.; Fung, I. Contribution to the atmospheric mineral aerosol load from land surface modification. J. Geophys. Res. 1995, 100, 707–726. [Google Scholar] [CrossRef]
  171. Marticorena, B.; Bergametti, G. Modeling the atmospheric dust cycle, I, Design of a soil-derived dust emission scheme. J. Geophys. Res. 1995, 100, 16415–16430. [Google Scholar] [CrossRef]
  172. Sun, W.Y.; Yang, K.J.-S.; Lin, N.-H. Numerical Simulations of Asian Dust-Aerosols and Regional Impacts on Weather and Climate-Part I: Control Case-PRCM Simulation without Dust-Aerosols. Aerosol Air Qual. Res. 2013, 13, 1630–1640. [Google Scholar] [CrossRef]
  173. Sun, W.Y.; Yang, K.J.-S.; Lin, N.-H. Numerical Simulations of Asian Dust-Aerosols and Regional Impacts on Weather and Climate-Part II: PRCM-Dust Model Simulation. Aerosol Air Qual. Res. 2013, 13, 1641–1654. [Google Scholar] [CrossRef]
  174. Koepke, P.; Hess, M.; Schult, I.; Shettle, E.P. Global Aerosol Data Set; Tech. Rcp. 243; Max Planck Inst. for Meteorol: Hamburg, Germany, 1997; 44p. [Google Scholar]
  175. Myhre, G.; Stordal, F. Global Sensitivity Experiments of the Radiative Forcing due to Mineral Aerosols. J. Geophys. Res. 2001, 106, 18193–18204. [Google Scholar] [CrossRef]
  176. Husar, R.B.; Tratt, D.M.; Schichtel, B.A.; Falke, S.R.; Li, F.; Jaffe, D.; Gasso, S.; Gill, T.; Laulainen, N.S.; Lu, F.; et al. Asian Dust Events of April 1998. J. Geophys. Res. 2001, 106, 18317–18330. [Google Scholar] [CrossRef]
  177. In, H.J.; Park, S.U. A Simulation of Long-range Transport of Yellow Sand Observed in April 1998 in Korea. Atmos. Environ. 2002, 36, 4173–4187. [Google Scholar] [CrossRef]
  178. Wu, C.C.; Yu, Y.C.; Hsu, W.R.; Hsu, K.J.; Sun, W.Y. Numerical Study on the Wind Fields and Atmospheric Transports of a Typical Winter Case in Taiwan and Surrounding Area. Atmos. Sci. 2003, 31, 29–54. (In Chinese) [Google Scholar]
  179. Sun, W.Y. Uncertainties in Climate Models, 10th ed.; McGraw-Hill Encyclopedia of Sciences & Technology: New York, NY, USA, 2006. [Google Scholar]
  180. Hsu, W.R.; Sun, W.Y. Numerical study of mesoscale cellular convection. Bound.-Layer Meteor. 1991, 57, 167–186. [Google Scholar] [CrossRef]
  181. Hsu, W.-R.; Sun, W.Y. Reply to “Comment on the numerical study of mesoscale cellular convection”. Bound.-Layer Meteorol. 1992, 60, 405–406. [Google Scholar] [CrossRef]
  182. Laprise, R. The Euler equations of motion with hydrostatic pressure as an independent variable. Mon. Wea. Rev. 1992, 120, 197–207. [Google Scholar] [CrossRef]
  183. Yeh, K.S.; Cote, J.; Gravel, S.; Methot, A.; Patoine, A.; Roch, M.; Staniforth, A. The CMC–MRB Global Environmental Multiscale (GEM) Model. Part III: Nonhydrostatic Formulation. Mon. Wea. Rev. 2002, 130, 339–356. [Google Scholar] [CrossRef]
  184. Janjic, Z.I. Nonhydrostatic model based on a new approach. Meteorol. Atmos. Phys. 2003, 82, 271–285. [Google Scholar] [CrossRef]
  185. Chen, S.H.; Sun, W.Y. The applications of the multigrid method and a flexible hybrid coordinate in a nonhydrostatic model. Mon. Wea. Rev. 2001, 129, 2660–2676. [Google Scholar] [CrossRef]
  186. Ogura, Y.; Phillips, N.A. Scale analysis of deep and shallow convection in the atmosphere. J. Atmos. Sci. 1962, 19, 173–179. [Google Scholar] [CrossRef]
  187. MacDonald, A.E.; Lee, J.L.; Sun, S. QNH: Design and test of a quasi-nonhydrostatic model for mesoscale weather prediction. Mon. Wea. Rev. 2000, 128, 1016–1036. [Google Scholar] [CrossRef]
  188. Klemp, J.; Wiehelmson, R. The Simulation of Three-Dimensional Convective Storm Dynamics. J. Atmos. Sci. 1978, 35, 1070–1096. [Google Scholar] [CrossRef]
  189. Sun, W.Y.; Hsu, W.R. Effect of surface friction on downslope wind and mountain waves. Terr. Atmos. Ocean. Sci. 2005, 16, 393–418. [Google Scholar] [CrossRef]
  190. Hsu, W.R.; Sun, W.Y. A time-split, forward-backward numerical model for solving a nonhydrostatic and compressible system of equations. Tellus 2001, 53, 279–299. [Google Scholar] [CrossRef]
  191. Hsu, H.H.; Yu, Y.C.; Kau, W.S.; Hsu, W.R.; Sun, W.Y. Simulation of the 1998 East Asian summer monsoon using Purdue regional model. J. Meteor. Soc. Jpn. 2004, 82, 1715–1733. [Google Scholar] [CrossRef]
  192. Sun, W.Y.; Hsu, W.R.; Chern, J.-D.; Chen, S.H.; Yang, K.J.-S.; Yeh, K.; Bosilovich, M.G.; Haines, P.A.; Min, K.-H.; Oh, T.-J.; et al. Purdue atmospheric models and applications. In Recent Progress in Atmospheric Sciences: Applications to the Asia-Pacific Region; Liao, K.N., Chou, M.D., Eds.; World Scientific Publishing: Singapore, 2009; 496p, pp. 200–230. [Google Scholar]
  193. Hsieh, M.N. A Monotone Semi-Lagrangian Advection Scheme in a 3D Nonhydrostatic Model and Applications in Squall Line Simulations. Ph.D. Thesis, NTU, Taipei, Taiwan, 2006; 146p. (In Chinese). [Google Scholar]
  194. Hsieh, M.-E.; Hsu, W.R.; Sun, W.Y. Applications of a three-dimensional nonhydrostatic atmospheric model on uniform flows over an idealized mountain. In Proceedings of the 2010 Computational Fluid Dynamics Taiwan, Chung-Li, Taiwan, 29–31 July 2010; pp. 1–12. [Google Scholar]
  195. Tsuboki, K.; Sakakibara, A. Numerical Prediction of High-Impact Weather Systems. In The Seventeenth IHP Training Course (International Hydrological Program); Nagoya University and UNESCO HYARC: Nagoya, Japan, 2007; 273p. [Google Scholar]
  196. Sun, W.Y.; Sun, O.M.T.; Tsuboki, K. A Modified Atmospheric Nonhydrostatic Model on Low Aspect Ratio Grids. Tellus A 2012, 64, 17516. [Google Scholar] [CrossRef]
  197. Sun, W.Y.; OSun, M.; Tsuboki, K. A Modified Atmospheric Nonhydrostatic Model on Low Aspect Ratio Grids-Part II. Tellus A 2013, 65, 19681. [Google Scholar] [CrossRef]
  198. Kasahara, A. Various vertical coordinate systems used for numerical weather prediction. Mon. Wea. Rev. 1974, 102, 509–522. [Google Scholar] [CrossRef]
  199. Sun, W.Y. Numerical simulation of a planetary boundary layer: Part I. Cloud free case. Beitr. Phys. Atmos. 1993, 66, 3–16. [Google Scholar]
  200. Sun, W.Y. An Efficient Forward Semi-Lagrangian Model. Open J. Fluid Dyn. (OJFD) 2023. accepted. [Google Scholar]
  201. Doswell, C.A., III. A kinematic analysis associated with a nondivergent vortex. J. Atmos. Sci. 1984, 41, 1242–1248. [Google Scholar] [CrossRef]
  202. Sun, W.Y.; Sun, O.M.T. Revisiting the parcel method and CAPE. Dyn. Atmos. Ocean. 2019, 86, 134–152. [Google Scholar] [CrossRef]
  203. Queney, P. The problem of airflow over mountains: A summary of theoretical studies. Bull. Am. Meteor. Soc. 1948, 29, 16–26. [Google Scholar] [CrossRef]
  204. Queney, P.; Corby, G.A.; Gerbier, N.; Koschmieder, H.; Zierep, J. The airflow over mountains. WMO Tech. 1960, 34, 135. [Google Scholar]
  205. Long, R.R. Some aspects of the flow of stratified fluids. Part I. A theoretical investigation. Tellus 1953, 5, 42–58. [Google Scholar] [CrossRef]
  206. Long, R.R. A laboratory model resembling the Bishop-waves phenomena. Bull. Am. Meteorol. Soc. 1953, 34, 205–211. [Google Scholar] [CrossRef]
  207. Smith, R.B. On severe downslope winds. J. Atmos. Sci. 1985, 42, 2597–2603. [Google Scholar] [CrossRef]
  208. Clark, T.L.; Peltier, W.R. Critical level reflection and the resonant growth of nonlinear mountain waves. J. Atmos. Sci. 1984, 41, 3122–3134. [Google Scholar] [CrossRef]
  209. Peltier, W.R.; Clark, T.L. The evolution and stability of finite-amplitude mountain waves. Part II: Surface drag and severe downslope windstorms. J. Atmos. Sci. 1979, 36, 1498–1529. [Google Scholar] [CrossRef]
  210. Doyle, J.D.; Durran, D.R.; Chen, C.; Colle, B.A.; Georgelin, M.; Grubisic, V.; Hsu, W.R.; Huang, C.Y.; Landau, D.; Lin, Y.L.; et al. An intercomparison of model predicted wave breaking for the 11 January 1972 Boulder windstorm. Mon. Wea. Rev. 2000, 128, 901–914. [Google Scholar] [CrossRef]
  211. Sun, W.Y. Numerical Study of Severe Downslope Storm. Weather Clim. Extrem. 2013, 2, 22–30. [Google Scholar] [CrossRef]
  212. Lilly, O.K. A severe downslope windstorm and aircraft turbulence induced by a mountain wave. J. Atmos. Sci. 1978, 35, 59–77. [Google Scholar] [CrossRef]
  213. Moncrieff, M.W.; Green, J.S.A. The propagation and transfer properties of steady convective overturning in shear. Q. J. R. Meteorol. Soc. 1972, 98, 336–352. [Google Scholar] [CrossRef]
  214. Riehl, H.; Malkus, J.S. On the heat balance in the equatorial trough zone. Gephysica 1958, 6, 503–538. [Google Scholar]
  215. Ebert, E.; Holland, G.J. Observations of record cold cloud-top temperatures in tropical cyclone Hilda (1990). Mon. Wea. Rev. 1992, 120, 2240–2251. [Google Scholar] [CrossRef]
  216. Black, R.A.; Bluestein, H.B.; Black, M.L. Unusually strong vertical motions in a Caribbean hurricane. Mon. Wea. Rev. 1994, 122, 2722–3739. [Google Scholar] [CrossRef]
  217. Sun, W.Y.; Chang, C.Z. Diffusion model for a convective layer: Part II: Plume released from a continuous point source. J. Clim. Appl. Meteor. 1986, 25, 1454–1463. [Google Scholar] [CrossRef]
  218. Sun, W.Y. Air pollution in a convective atmosphere. In Library of Environmental Control Technology; Cheremisinoff, Ed.; Gulf Publishing: Houston, TX, USA, 1988; pp. 515–546. [Google Scholar]
  219. Sun, W.Y. Diffusion modeling in a convective boundary layer. Atmos. Environ. 1989, 23, 1205–1217. [Google Scholar]
  220. Haines, P.A.; Sun, W.Y.; Chen, S.H.; Hsu, W.R.; Hsieh, M.E. NTU/PU model simulations and observed flow over mountain. Terr. Atmos. Ocean. Sci. 2019, 30, 171–184. [Google Scholar] [CrossRef]
  221. Gal-Chen, T.; Somerville, R.C.J. Numerical solution of the Navier-Stokes equations with topography. J. Comput. Phys. 1975, 17, 276–309. [Google Scholar] [CrossRef]
  222. Haines, P.A.; Grove, D.J.; Sun, W.Y.; Hsu, W.R. Observations and Numerical Modeling of Sub and Super-critical Flow at White Sands Missile Range. In Proceedings of the 12th Conference on Mountain Meteorology, the American Meteorological Society, Santa Fe, NM, USA, 28 August–2 September 2006. [Google Scholar]
  223. Sun, W.Y. Coordinates over complex terrain in atmospheric model. J. Atmos. Sci. Res. 2021, 4, 3–49. [Google Scholar] [CrossRef]
  224. Pacanowski, R.C. MOM 2 Documentation User’s Guide and Reference Manual Version 1.0; GFDL Ocean Technical Report #3; GFDL, NOAA: Princeton, NJ, USA, 1995; 233p. [Google Scholar]
  225. Saito, K.; Kato, K.T.; Eito, H.; Muroi, C. Numerical Prediction Division Unified Nonhydrostatic Model; MRI Report 42; Documentation of the Meteorological Research Institute: Singapore, 2001; 133p.
  226. Staniforth, A.; White, W.; Wood, N.; Thuburn, J.; Zerroukat, M.; Corderro, E. Unified Model Documentation Paper. No. 15. Joy of U.M. 6.0-Model Formulation. Dynamics Research Numerical Weather Prediction; Met Office: Devan, UK, 2004.
  227. Lin, M.-Y.; Sun, W.Y.; Chiou, M.-D.; Chen, C.-Y.; Cheng, H.-Y.; Chen, C.-H. Development and evaluation of a storm surge warning system in Taiwan. Ocean. Dyn. 2018, 68, 1025–1049. [Google Scholar] [CrossRef]
  228. Pielke, R.A.; Cotton, W.R.; Walko, R.L.; Tremback, C.J.; Lyons, W.A.; Grasso, L.D.; Nicholls, M.E.; Moran, M.D.; Wesley, D.A.; Lee, T.J.; et al. A comprehensive meteorological modeling system-RAMS. Meteor. Atmos. Phys. 1992, 49, 69–91. [Google Scholar] [CrossRef]
  229. Kapitza, H.; Eppel, D.P. The non-hydrostatic mesoscale GESIMA. Part I: Dynamical equations and tests. Beitr. Phys. Atmos. 1992, 65, 129–145. [Google Scholar]
  230. Skamarock, W.C.; Klemp, J.B.; Dudhia, J.; Gill, D.O.; Barker, D.M.; Duda, M.G.; Huang, X.-Y.; Wang, W.; Powers, J.G. A Description of Advanced Research WRF Version 3; NCAR/TN–475+STR; NCAR: Boulder, CO, USA, 2008; 113p. [Google Scholar]
  231. MacCall, B.-T.; Wang, Y.; Sun, W.Y. A New Semi-Implicit Time Integration Scheme for the Time-Dependent Atmospheric Boundary Layer Environment (ABLE) Model; ARL-TR-7428; ARL: Aberdeen Proving Ground, MD, USA, 2015. [Google Scholar]
  232. Sun, W.Y. Numerical simulation of a planetary boundary layer: Part II. Cloudy case. Beitr. Phys. Atmos. 1993, 66, 17–30. [Google Scholar]
  233. Deardorff, J.W.; Willis, G.E. A parameterization of diffusion into the mixed layer. J. Appl. Met. 1975, 14, 1451–1458. [Google Scholar] [CrossRef]
  234. Willis, G.E.; Deardorff, J.W. A laboratory study of dispersion from an elevated source within a modeled convective planetary boundary layer. Atmos. Environ. 1978, 12, 1305–1311. [Google Scholar] [CrossRef]
  235. Willis, G.E.; Deardorff, J.W. A laboratory study of dispersion from a source in the middle of the convective mixed layer. Atmos. Environ. 1981, 15, 109–117. [Google Scholar] [CrossRef]
  236. Sun, W.Y.; Sun, O.M.T. Backward integration of diffusion equation. Aerosol Air Qual. Res. 2017, 17, 278–289. [Google Scholar] [CrossRef]
  237. Arnold, V.I.; Avez, A. Problèmes Ergodiques de la Mécanique Classique; Gauthier-Villars: Paris, France, 1967. (In French) [Google Scholar]
  238. Strogatz, S.H. Nonlinear Dynamics and Chaos. With Applications to Physics, Biology, Chemistry, and Engineering; Westpress View: Boulder, CO, USA, 2015; 513p. [Google Scholar]
  239. Lyapunov, A.M. Problème général de la stabilité du mouvement. Ann. Fac. Sci. Toulouse Mathématiques Série 2 1947, 9, 203–474. [Google Scholar]
  240. Lorenz, E.N. Predictability: Does the flap of a butterfly’s wings in Brazil set off a tornado in Texas? In Proceedings of the 139th Meeting of AAAS Section on Environmental Sciences, New Approaches to Global Weather: GARP, AAAS, Cambridge, MA, USA, 29 December 1972. 5p. [Google Scholar]
  241. Lorenz, E.N. Atmospheric predictability experiments with a large numerical model. Tellus 1982, 34, 505–513. [Google Scholar] [CrossRef]
  242. Sun, W.Y.; Hsu, W.R.; Haines, P.K.; Hsieh, M.E.; Chen, S.H.; Grove, D. Numerical simulations over mountains. In Proceedings of the Third International Workshop on Next- Generation NWP Models: Bridging Parameterization, Explicit Clouds and Large Eddies, Jeju Island, Republic of Korea, 29 August–1 September 2010. [Google Scholar]
  243. Kreiss, H.-O.; Lorenz, J. Initial-Boundary Value Problems and the Navier-Stokes Equations; SIAM, Publications Library: Philadelphia, PA, USA, 2004; 488p. [Google Scholar]
Figure 1. Magnitude of n(−1)n−1: h or u without smoothing, (red dashed line), and n(−0.91)n−1 with smoothing (blue solid line) as function of nth time step, where λ′ = Shλ = 0.91 for η = 0.15. (Sun 2010) [23].
Figure 1. Magnitude of n(−1)n−1: h or u without smoothing, (red dashed line), and n(−0.91)n−1 with smoothing (blue solid line) as function of nth time step, where λ′ = Shλ = 0.91 for η = 0.15. (Sun 2010) [23].
Atmosphere 14 01324 g001
Figure 2. Growth rate for forward–backward scheme as function of wavelength (Δx) and viscosity (ν) at: (a) Courant number Co = 1; (b) Co = 0.8. (Sun 2010) [23].
Figure 2. Growth rate for forward–backward scheme as function of wavelength (Δx) and viscosity (ν) at: (a) Courant number Co = 1; (b) Co = 0.8. (Sun 2010) [23].
Atmosphere 14 01324 g002
Figure 3. Growth rate for Leapfrog scheme as function of wavelength (Δx) and viscosity (ν) at: (a) Courant number Co = 0.5; (b) Co = 0.4. (Sun 2010) [23].
Figure 3. Growth rate for Leapfrog scheme as function of wavelength (Δx) and viscosity (ν) at: (a) Courant number Co = 0.5; (b) Co = 0.4. (Sun 2010) [23].
Atmosphere 14 01324 g003
Figure 4. Vertical cross section of water depth h at distance R from the center at t = 0.69 s from 2nd-order scheme with Shuman smoothing coefficient α = 0 (×) and α = 0.2 (▲), and 4th-order scheme with α = 0.2 (□). (Sun 2011) [24].
Figure 4. Vertical cross section of water depth h at distance R from the center at t = 0.69 s from 2nd-order scheme with Shuman smoothing coefficient α = 0 (×) and α = 0.2 (▲), and 4th-order scheme with α = 0.2 (□). (Sun 2011) [24].
Atmosphere 14 01324 g004
Figure 5. Vertical cross sections of water depth at t = 1, 3, 10, 15 and 25 s from the 2nd-order scheme with smoothing coefficient α = 0.2. Sun (2011) [24].
Figure 5. Vertical cross sections of water depth at t = 1, 3, 10, 15 and 25 s from the 2nd-order scheme with smoothing coefficient α = 0.2. Sun (2011) [24].
Atmosphere 14 01324 g005
Figure 6. Height of Rossby wave propagated westward (a): analytic solution at t = 120 (dashed line) and simulated height at t = 0 (at x = 0) and t = 30, 60, 90, 120 (solid lines), the contours from 1.0 to 1.14 with interval 0.02; (b) simulated height (solid line) at t = 120 from linearized wave, contours are from 0.95 to 1.11 with interval of 0.01, and initial condition at t = 0 (dashed line) (Sun and Sun 2013) [48].
Figure 6. Height of Rossby wave propagated westward (a): analytic solution at t = 120 (dashed line) and simulated height at t = 0 (at x = 0) and t = 30, 60, 90, 120 (solid lines), the contours from 1.0 to 1.14 with interval 0.02; (b) simulated height (solid line) at t = 120 from linearized wave, contours are from 0.95 to 1.11 with interval of 0.01, and initial condition at t = 0 (dashed line) (Sun and Sun 2013) [48].
Atmosphere 14 01324 g006
Figure 7. Simulated solution of linearized equations (purple solid line), soliton of nonlinear equations (blue dashed line), and Boyd’s solutions from nonlinear equations (bold solid line) and linearized equations (black dotted line) at t = 100. (Sun and Sun 2013) [48].
Figure 7. Simulated solution of linearized equations (purple solid line), soliton of nonlinear equations (blue dashed line), and Boyd’s solutions from nonlinear equations (bold solid line) and linearized equations (black dotted line) at t = 100. (Sun and Sun 2013) [48].
Atmosphere 14 01324 g007
Figure 8. Simulated heights and streamlines in coordinate moving with wave from sum of zeroth-and first-order solutions as initial conditions with Δx = Δy = 0.2 and smooth coefficient = 0.1: (a) at t = 121.5 for B = 0.395; and (b) at t = 100.9389 for B = 0.6. (Sun and Sun 2013) [48].
Figure 8. Simulated heights and streamlines in coordinate moving with wave from sum of zeroth-and first-order solutions as initial conditions with Δx = Δy = 0.2 and smooth coefficient = 0.1: (a) at t = 121.5 for B = 0.395; and (b) at t = 100.9389 for B = 0.6. (Sun and Sun 2013) [48].
Atmosphere 14 01324 g008
Figure 9. (a): Time sequences of water surface height hf and velocity u (velocity at inflow u0 = 42 cm s−1 and Froude number at inflow Fo = 0.3), red line at t = 2.4 s, green at t = 4.8 s, dark blue at t = 7.2 s, light blue at t = 9.6 s, purple at t = 96 s; (b): Left: mass flux without terrain J1 (red), mass flux with terrain J2 (blue), and J (green) = J1 + J2; Right: h × u (purple), Bernoulli function Be (black); dashed lines at t = 0, solid lines at t = 96 s. (Sun 2018) [54].
Figure 9. (a): Time sequences of water surface height hf and velocity u (velocity at inflow u0 = 42 cm s−1 and Froude number at inflow Fo = 0.3), red line at t = 2.4 s, green at t = 4.8 s, dark blue at t = 7.2 s, light blue at t = 9.6 s, purple at t = 96 s; (b): Left: mass flux without terrain J1 (red), mass flux with terrain J2 (blue), and J (green) = J1 + J2; Right: h × u (purple), Bernoulli function Be (black); dashed lines at t = 0, solid lines at t = 96 s. (Sun 2018) [54].
Atmosphere 14 01324 g009aAtmosphere 14 01324 g009b
Figure 10. (a) Trajectories for A1 (contours from 1.2 × 10−5 to 2.2 × 10−5 m−1s−1), A2 (contours from 1.2 × 10−5 to 1.5 × 10−5 m−1s−1), and A3 (contour is 6 × 10−6 m−1 s−1) with α = 0, black line is hs = 5 m; (b) streamline, wind V (red) and y-component mass flux, hv (green); PV budget for: (c) local rate of change, (d) advection, and (e) total derivative o for A2 at t = 9.6 × 104 s (Sun, 2016) [55].
Figure 10. (a) Trajectories for A1 (contours from 1.2 × 10−5 to 2.2 × 10−5 m−1s−1), A2 (contours from 1.2 × 10−5 to 1.5 × 10−5 m−1s−1), and A3 (contour is 6 × 10−6 m−1 s−1) with α = 0, black line is hs = 5 m; (b) streamline, wind V (red) and y-component mass flux, hv (green); PV budget for: (c) local rate of change, (d) advection, and (e) total derivative o for A2 at t = 9.6 × 104 s (Sun, 2016) [55].
Atmosphere 14 01324 g010
Figure 11. Trajectories of PV with initial surface depression h o = 3 m, initial radius R = 24 km (Sun, 2016) [55].
Figure 11. Trajectories of PV with initial surface depression h o = 3 m, initial radius R = 24 km (Sun, 2016) [55].
Atmosphere 14 01324 g011
Figure 12. (a) PV (thick black), h (white), streamline (thin back), and speed (shaded color) at t = 0; (b) Same as (a) but at t = 20 s (Sun and Oh 2022) [62]. (c) ∏ ≡ PV (thick black), h (white), streamline (thin back), and ∂∏/∂t (shaded color) at t = 20 s; (d) Same as (c) but at t = 60 s; (e) Same as (c) but at t = 100 s, (f) Same as (c) but at t = 120 s (Sun and Oh 2022) [62].
Figure 12. (a) PV (thick black), h (white), streamline (thin back), and speed (shaded color) at t = 0; (b) Same as (a) but at t = 20 s (Sun and Oh 2022) [62]. (c) ∏ ≡ PV (thick black), h (white), streamline (thin back), and ∂∏/∂t (shaded color) at t = 20 s; (d) Same as (c) but at t = 60 s; (e) Same as (c) but at t = 100 s, (f) Same as (c) but at t = 120 s (Sun and Oh 2022) [62].
Atmosphere 14 01324 g012aAtmosphere 14 01324 g012b
Figure 13. (a): Schematic distribution of cumuliform clouds near a wave trough in the easterlies; the orientation of cloud rows and the spacing between them was measured exactly from the film. The direction of the wind in the lower layer is indicated by an arrow. (Malkus and Riehl, 1964) [78]. (b): Schematic distribution of cumuliform clouds near a wave trough in the easterlies; the orientation of cloud rows and the spacing between them was measured from the film. The wind in the lower layer is indicated by VB, wind at upper layer is VT, (Sun 1978) [77].
Figure 13. (a): Schematic distribution of cumuliform clouds near a wave trough in the easterlies; the orientation of cloud rows and the spacing between them was measured exactly from the film. The direction of the wind in the lower layer is indicated by an arrow. (Malkus and Riehl, 1964) [78]. (b): Schematic distribution of cumuliform clouds near a wave trough in the easterlies; the orientation of cloud rows and the spacing between them was measured from the film. The wind in the lower layer is indicated by VB, wind at upper layer is VT, (Sun 1978) [77].
Atmosphere 14 01324 g013
Figure 14. (a) Stream function at 2.274 × 104 s for Case D. Contours are from 0.16 to +0.18 with interval 2.0 × 10−2 m2 s−1; (b) Vertical cross section of the heating rate ( Q t ) at 2.27 × 104 s for Case D. Contours are from −1.800 to 3.9 (×10−5 K s−1) with interval 3.0 × 10−6 Ks−1 and labels scaled by 10−7 (c). Stream function at 2.653 × 104 s for Case D. Contours are from −0.56 to +0.64 with interval of 8.00 × 10−2 m2 s−1. (Sun 1984b) [80].
Figure 14. (a) Stream function at 2.274 × 104 s for Case D. Contours are from 0.16 to +0.18 with interval 2.0 × 10−2 m2 s−1; (b) Vertical cross section of the heating rate ( Q t ) at 2.27 × 104 s for Case D. Contours are from −1.800 to 3.9 (×10−5 K s−1) with interval 3.0 × 10−6 Ks−1 and labels scaled by 10−7 (c). Stream function at 2.653 × 104 s for Case D. Contours are from −0.56 to +0.64 with interval of 8.00 × 10−2 m2 s−1. (Sun 1984b) [80].
Atmosphere 14 01324 g014aAtmosphere 14 01324 g014b
Figure 15. Surface weather depiction at 2000 GMT (1500 EST) 3 April 1974. Three active squall lines were all producing tornadoes in Illinois, Indiana, and Tennessee. (After Agee et al., 1975) [86].
Figure 15. Surface weather depiction at 2000 GMT (1500 EST) 3 April 1974. Three active squall lines were all producing tornadoes in Illinois, Indiana, and Tennessee. (After Agee et al., 1975) [86].
Atmosphere 14 01324 g015
Figure 16. (a) A composite picture of the vertical velocity in the lower atmosphere from 1400 LST of day 14 to 1400 LST of day 16 at 12 h intervals for case B. The maximum is 20.1 cm s−1 and the minimum −32.4 cm s−1 at 1400 LST of day 14; the maximum 26.5 cm s−1 and minimum −24.4 cm s−1 at 0200 LST of day 15; the maximum 45.6 cm s−1 and minimum −21.6 cm s−1 at 1400 LST of day 15; the maximum 52.6 cm s−1 and minimum −42.4 cm s−1 at 0200 LST of day 16; the maximum 41.4 cm s−1 and minimum −51.8 cm s−1 at 1400 LST of day 16. (b) Temporal variation of z velocity at x = 768 km from 00 LST of day 0 to 06 LST of day 8 for case C. The maximum is 3.1 cm s−1, the minimum −7.53 cm s−1, and the interval is 1.06 cm s−1 (Sun and Orlanski 1981a) [87].
Figure 16. (a) A composite picture of the vertical velocity in the lower atmosphere from 1400 LST of day 14 to 1400 LST of day 16 at 12 h intervals for case B. The maximum is 20.1 cm s−1 and the minimum −32.4 cm s−1 at 1400 LST of day 14; the maximum 26.5 cm s−1 and minimum −24.4 cm s−1 at 0200 LST of day 15; the maximum 45.6 cm s−1 and minimum −21.6 cm s−1 at 1400 LST of day 15; the maximum 52.6 cm s−1 and minimum −42.4 cm s−1 at 0200 LST of day 16; the maximum 41.4 cm s−1 and minimum −51.8 cm s−1 at 1400 LST of day 16. (b) Temporal variation of z velocity at x = 768 km from 00 LST of day 0 to 06 LST of day 8 for case C. The maximum is 3.1 cm s−1, the minimum −7.53 cm s−1, and the interval is 1.06 cm s−1 (Sun and Orlanski 1981a) [87].
Atmosphere 14 01324 g016
Figure 17. GOES East satellite picture at 2130 GMT 5 September 1974, showing three cloud bands aligned parallel to the coast of South America. Their horizontal wavelength is of a few hundred kilometers. There is a clear sky off the coast (Sun and Orlanski 1981a) [87].
Figure 17. GOES East satellite picture at 2130 GMT 5 September 1974, showing three cloud bands aligned parallel to the coast of South America. Their horizontal wavelength is of a few hundred kilometers. There is a clear sky off the coast (Sun and Orlanski 1981a) [87].
Atmosphere 14 01324 g017
Figure 18. Diurnal variation of log10 (μ) in atmosphere and ocean (Sun and Sun 2020) [91].
Figure 18. Diurnal variation of log10 (μ) in atmosphere and ocean (Sun and Sun 2020) [91].
Atmosphere 14 01324 g018
Figure 19. (a) Time evolutions of velocity at z = 600 m (red), and z = 0 (bule, multiplied by 100), box indicates at t = 0, × at t = 1 s, o at 12L12, Δ at 24L12. (b) Time sequences of u (red line with ×), v (red line) at 600 m, and 10 × u (blue line with ×), 10 × v (blue line) at z = 0 for Case C (30° N). (c) Simulated velocities (u, v) in atmosphere on D12: at 03L (thin dotted–dashed green); 12L (thick dashed red); 15 L (thin dotted–dashed magenta); and 24L (dashed black), and hourly velocity at z = 200 m (blue with *), z = 600 m (red with *). (d) Oceanic Ekman spiral at every 3 h and 0.4 × stress at 24L12 (dashed black). (e) Simulated (u, v) in ocean at 12L (thick dashed red) and 24L (thick dashed black), hourly (u, v) at z = 0 (thin magenta with *), −4 m (thin green), −10 m (thin blue with *), −20 m (thin red), −40 m (thin black); analytic velocity (thick green), daily average (thick blue); and 0.4 × stress at 12L (thick red), at 24L (thick black) on Day 12. (f) Simulated daily mean Ekman spiral (blue, red), and observation from Price et al. (1986) at 30° N. Hodograph starting from 2 m down to 70 m with 4 m interval in solid black curve of observed and red line of simulation. (g) Hourly surface stress ( τ ˜ x (blue *), τ ˜ y (blue o) and mass flux f l u x ( u ˜ ) (red *), f l u x ( v ˜ ) (red o) and their daily averages (at 25 h). (h) Simulated (u, v) in water on Day 12, the top row shows the surface stress, and the last column is daily mean. (Sun and Sun 2020) [91].
Figure 19. (a) Time evolutions of velocity at z = 600 m (red), and z = 0 (bule, multiplied by 100), box indicates at t = 0, × at t = 1 s, o at 12L12, Δ at 24L12. (b) Time sequences of u (red line with ×), v (red line) at 600 m, and 10 × u (blue line with ×), 10 × v (blue line) at z = 0 for Case C (30° N). (c) Simulated velocities (u, v) in atmosphere on D12: at 03L (thin dotted–dashed green); 12L (thick dashed red); 15 L (thin dotted–dashed magenta); and 24L (dashed black), and hourly velocity at z = 200 m (blue with *), z = 600 m (red with *). (d) Oceanic Ekman spiral at every 3 h and 0.4 × stress at 24L12 (dashed black). (e) Simulated (u, v) in ocean at 12L (thick dashed red) and 24L (thick dashed black), hourly (u, v) at z = 0 (thin magenta with *), −4 m (thin green), −10 m (thin blue with *), −20 m (thin red), −40 m (thin black); analytic velocity (thick green), daily average (thick blue); and 0.4 × stress at 12L (thick red), at 24L (thick black) on Day 12. (f) Simulated daily mean Ekman spiral (blue, red), and observation from Price et al. (1986) at 30° N. Hodograph starting from 2 m down to 70 m with 4 m interval in solid black curve of observed and red line of simulation. (g) Hourly surface stress ( τ ˜ x (blue *), τ ˜ y (blue o) and mass flux f l u x ( u ˜ ) (red *), f l u x ( v ˜ ) (red o) and their daily averages (at 25 h). (h) Simulated (u, v) in water on Day 12, the top row shows the surface stress, and the last column is daily mean. (Sun and Sun 2020) [91].
Atmosphere 14 01324 g019aAtmosphere 14 01324 g019b
Figure 20. (a) Model simulations, (b) observed surface stress and velocity vectors near 60° S (Lenn and Chereskin 2009) [100], (Sun and Sun 2020) [91].
Figure 20. (a) Model simulations, (b) observed surface stress and velocity vectors near 60° S (Lenn and Chereskin 2009) [100], (Sun and Sun 2020) [91].
Atmosphere 14 01324 g020
Figure 21. Schematic diagram of PRCM (Purdue Regional Climate Model).
Figure 21. Schematic diagram of PRCM (Purdue Regional Climate Model).
Atmosphere 14 01324 g021
Figure 22. Schematic vertical cross section across the simulated dryline at 1630 LT 8 June 1966. J, U, and D denoted the location of the cores of the jet in the y-direction, the upward motion and downward motion, respectively. The dashed lines are the contour lines of 2 cm s−1 for the vertical velocity (Sun and Ogura 1979) [103].
Figure 22. Schematic vertical cross section across the simulated dryline at 1630 LT 8 June 1966. J, U, and D denoted the location of the cores of the jet in the y-direction, the upward motion and downward motion, respectively. The dashed lines are the contour lines of 2 cm s−1 for the vertical velocity (Sun and Ogura 1979) [103].
Atmosphere 14 01324 g022
Figure 23. (a) θv,; (b) qw (total water content) at 1800 CST on day 1; (c) southerly wind v and low-lvel-jet (LLJ) at 0200 CST on day 2; (d) time–space variation of qw at z~10 m, with a contour interval of 1 g Kg−1, labels scale by 10. (Sun and Wu 1992) [104].
Figure 23. (a) θv,; (b) qw (total water content) at 1800 CST on day 1; (c) southerly wind v and low-lvel-jet (LLJ) at 0200 CST on day 2; (d) time–space variation of qw at z~10 m, with a contour interval of 1 g Kg−1, labels scale by 10. (Sun and Wu 1992) [104].
Atmosphere 14 01324 g023aAtmosphere 14 01324 g023b
Figure 24. (a) Streamlines at z = 1000 m at 2000 LST 16 May (Sun and Chern 1993). (b) Observed surface data at 2000 LST 16 May 1987 [130].
Figure 24. (a) Streamlines at z = 1000 m at 2000 LST 16 May (Sun and Chern 1993). (b) Observed surface data at 2000 LST 16 May 1987 [130].
Atmosphere 14 01324 g024
Figure 25. (a) Simulated surface pressure deviation from −240 to 120 Pa with interval of 20 Pa. (b) Streamline at z = 1000 at 0200 LST 17 May (Sun and Chern 1993) [130].
Figure 25. (a) Simulated surface pressure deviation from −240 to 120 Pa with interval of 20 Pa. (b) Streamline at z = 1000 at 0200 LST 17 May (Sun and Chern 1993) [130].
Atmosphere 14 01324 g025
Figure 26. (a) Observed sea surface pressure, (b) streamline at 900 mb observed by airplane at 0200 LST 17 (1800 UTC 16) May 1989 (Kuo and Chen 1990) [132].
Figure 26. (a) Observed sea surface pressure, (b) streamline at 900 mb observed by airplane at 0200 LST 17 (1800 UTC 16) May 1989 (Kuo and Chen 1990) [132].
Atmosphere 14 01324 g026
Figure 27. (Left) simulated cloud water at 700 mb at 0000 UTC May 17 (after 48 h integration), and (right) satellite image at 2137 UTC May 16.
Figure 27. (Left) simulated cloud water at 700 mb at 0000 UTC May 17 (after 48 h integration), and (right) satellite image at 2137 UTC May 16.
Atmosphere 14 01324 g027
Figure 28. (Left) radar observed wind and vortex, and (Right) Simulated wind and Denver Cyclone (red) on 13 February 1990 during WISP, color bar indicates elevation height in m (Haines et al., 1997) [134].
Figure 28. (Left) radar observed wind and vortex, and (Right) Simulated wind and Denver Cyclone (red) on 13 February 1990 during WISP, color bar indicates elevation height in m (Haines et al., 1997) [134].
Atmosphere 14 01324 g028
Figure 29. (a) Isochrones of leading edge of the cold air from the surface and satellite data. hhZdd is to denote the time and date. (b) Simulated streamline at 06Z15, 1997 at z ≈ 25 m (Sun and Chern 2004). (c) ECMWF reanalysis streamline at 06Z15 at z ≈ 25 m [71].
Figure 29. (a) Isochrones of leading edge of the cold air from the surface and satellite data. hhZdd is to denote the time and date. (b) Simulated streamline at 06Z15, 1997 at z ≈ 25 m (Sun and Chern 2004). (c) ECMWF reanalysis streamline at 06Z15 at z ≈ 25 m [71].
Atmosphere 14 01324 g029
Figure 30. Locations of the cyclones at 6 h interval. Solid circles represent observation from ECMWF/TOGA analysis. Open from PRCM. The cyclones over Nevada (eastern Idaho) first appeared at 1200 UTC1 (0000 UTC 2) March 1985 (Chern 1994) [110].
Figure 30. Locations of the cyclones at 6 h interval. Solid circles represent observation from ECMWF/TOGA analysis. Open from PRCM. The cyclones over Nevada (eastern Idaho) first appeared at 1200 UTC1 (0000 UTC 2) March 1985 (Chern 1994) [110].
Atmosphere 14 01324 g030
Figure 31. Observed (left) and simulated (right) surface pressure at low pressure centers in Nevada (filled bar) and Idaho (hatched bar) (Chern 1994) [110].
Figure 31. Observed (left) and simulated (right) surface pressure at low pressure centers in Nevada (filled bar) and Idaho (hatched bar) (Chern 1994) [110].
Atmosphere 14 01324 g031
Figure 32. (a) ECMWF geopotential Z; (b) model simulated Z at 500 m at 00Z/03/1985 after 2 days integration (Chern 1994) [110].
Figure 32. (a) ECMWF geopotential Z; (b) model simulated Z at 500 m at 00Z/03/1985 after 2 days integration (Chern 1994) [110].
Atmosphere 14 01324 g032
Figure 33. Simulated from PRCM and NCDC observation of daily integrated precipitation over flooding area in the Midwestern US for: (a) June associated with synoptic waves and (b) July 1993 associated with mesoscale convective system (Bosilovich and Sun 1999a) [136].
Figure 33. Simulated from PRCM and NCDC observation of daily integrated precipitation over flooding area in the Midwestern US for: (a) June associated with synoptic waves and (b) July 1993 associated with mesoscale convective system (Bosilovich and Sun 1999a) [136].
Atmosphere 14 01324 g033
Figure 34. Observed and simulation of temperature, geopotential, moisture, and wind at 500 mb for June 1988 drought in USA (Sun et al., 2004) [141].
Figure 34. Observed and simulation of temperature, geopotential, moisture, and wind at 500 mb for June 1988 drought in USA (Sun et al., 2004) [141].
Atmosphere 14 01324 g034aAtmosphere 14 01324 g034b
Figure 35. Observed (green circle) and simulated (red x for Case A, blue line for Case B) during November 1970 to May 1974 at Sleepers Watershed. (Case A: Fresh snow density is given by observed mean value of each winter. Case B: Fresh snow density = 150 kg m−3): (a) snow depth, (b) ground temperature at layer 1–5850 and (c) snowpack temperatures at 6, 12, and 24 inches above ground (Sun and Chern 2005) [114].
Figure 35. Observed (green circle) and simulated (red x for Case A, blue line for Case B) during November 1970 to May 1974 at Sleepers Watershed. (Case A: Fresh snow density is given by observed mean value of each winter. Case B: Fresh snow density = 150 kg m−3): (a) snow depth, (b) ground temperature at layer 1–5850 and (c) snowpack temperatures at 6, 12, and 24 inches above ground (Sun and Chern 2005) [114].
Atmosphere 14 01324 g035aAtmosphere 14 01324 g035bAtmosphere 14 01324 g035c
Figure 36. Left panels: the observed snow depth (top), simulation with a single snow layer (middle, CTL), simulated snow depth with comprehensive snow–soil physics (bottom, EXP 1) in March of 1997. Right panels: observed wind at 200 hPa (top), the simulation from single snow layer (middle), and simulation from comprehensive snow–soil (bottom). (Min and Sun 2015) [161].
Figure 36. Left panels: the observed snow depth (top), simulation with a single snow layer (middle, CTL), simulated snow depth with comprehensive snow–soil physics (bottom, EXP 1) in March of 1997. Right panels: observed wind at 200 hPa (top), the simulation from single snow layer (middle), and simulation from comprehensive snow–soil (bottom). (Min and Sun 2015) [161].
Atmosphere 14 01324 g036aAtmosphere 14 01324 g036b
Figure 37. Time series of mean sea level pressure [hPa] from ECMWF (green circle), CTL (red cross), and EXP1 (solid line) simulations in Red River Region in March of 1997 (Min and Sun 2015) [161].
Figure 37. Time series of mean sea level pressure [hPa] from ECMWF (green circle), CTL (red cross), and EXP1 (solid line) simulations in Red River Region in March of 1997 (Min and Sun 2015) [161].
Atmosphere 14 01324 g037
Figure 38. (a) ECMWF, PRCM in mean sea level pressures and difference, bias = −0.15 Pa, RMSE = 0.60 Pa, pattern correlation coefficient = 0.95 Yu et al. (2004a, b) [162,163]. (b) ECMWF, PRCM surface temperature and difference, bias 0.47 K, RMSE = 0.72 K, pattern correlation coefficient 0.99 Yu et al. (2004a, b) [162,163]. (c) ECMWF, PRCM precipitation and difference, terrain height indicated by contours in, pattern correlation is 0.35 (Yu et al., 2004a, b) [162,163].
Figure 38. (a) ECMWF, PRCM in mean sea level pressures and difference, bias = −0.15 Pa, RMSE = 0.60 Pa, pattern correlation coefficient = 0.95 Yu et al. (2004a, b) [162,163]. (b) ECMWF, PRCM surface temperature and difference, bias 0.47 K, RMSE = 0.72 K, pattern correlation coefficient 0.99 Yu et al. (2004a, b) [162,163]. (c) ECMWF, PRCM precipitation and difference, terrain height indicated by contours in, pattern correlation is 0.35 (Yu et al., 2004a, b) [162,163].
Atmosphere 14 01324 g038aAtmosphere 14 01324 g038bAtmosphere 14 01324 g038c
Figure 39. Vorticity at 850 hPa (red for ECMWF, blue for PRCM) and correlation for summers (May, June, July and August) from in 1991–1998 (Yu et al., 2004a, b) [162,163].
Figure 39. Vorticity at 850 hPa (red for ECMWF, blue for PRCM) and correlation for summers (May, June, July and August) from in 1991–1998 (Yu et al., 2004a, b) [162,163].
Atmosphere 14 01324 g039aAtmosphere 14 01324 g039b
Figure 40. The 30–60-day band-pass filtered GPCP (left) and PRM (right) precipitation for the 5-day periods: (a,e) 21–25 May, (b,f) 31 May−4 June, (c,g) 10–14 June, and (d,h) 20–24 June. In 1998. The contour interval is 2 mm/day, and the shading indicates the positive anomalies (Hsu et al., 2004) [164].
Figure 40. The 30–60-day band-pass filtered GPCP (left) and PRM (right) precipitation for the 5-day periods: (a,e) 21–25 May, (b,f) 31 May−4 June, (c,g) 10–14 June, and (d,h) 20–24 June. In 1998. The contour interval is 2 mm/day, and the shading indicates the positive anomalies (Hsu et al., 2004) [164].
Atmosphere 14 01324 g040
Figure 41. The schematic illustration of the components in the integrated PRCM-Dust model (Yang 2004a) [127].
Figure 41. The schematic illustration of the components in the integrated PRCM-Dust model (Yang 2004a) [127].
Atmosphere 14 01324 g041
Figure 42. (a) ECMWF Mean 850 hPa geopotential (m) during April 8–24; (b) wind at 700 hPa (unit: m/s). (Yang 2004a) [127].
Figure 42. (a) ECMWF Mean 850 hPa geopotential (m) during April 8–24; (b) wind at 700 hPa (unit: m/s). (Yang 2004a) [127].
Atmosphere 14 01324 g042
Figure 43. (a) Simulated mean 850 hPa geopotential (m) during 8–24 April; (b) wind at 700 hPa (unit: m/s). (Sun et al., 2013a, b) [172,173].
Figure 43. (a) Simulated mean 850 hPa geopotential (m) during 8–24 April; (b) wind at 700 hPa (unit: m/s). (Sun et al., 2013a, b) [172,173].
Atmosphere 14 01324 g043
Figure 44. Total column aerosol optical thickness at 500 nm from PRCM-Dust model and AERONET measurement site at Dalanzadgad, Mongolia. (Sun et al., 2013a, b) [172,173].
Figure 44. Total column aerosol optical thickness at 500 nm from PRCM-Dust model and AERONET measurement site at Dalanzadgad, Mongolia. (Sun et al., 2013a, b) [172,173].
Atmosphere 14 01324 g044
Figure 45. (a) Simulated mixing ratio of dust near 550 hPa and (b) TOM AI maps at 00Z 22 April 1998. (Sun et al., 2013a, b) [172,173].
Figure 45. (a) Simulated mixing ratio of dust near 550 hPa and (b) TOM AI maps at 00Z 22 April 1998. (Sun et al., 2013a, b) [172,173].
Atmosphere 14 01324 g045
Figure 46. (a) PRCM (Dust)–PRCM (without Dust) daily mean surface temperature at 00Z for 8–24 April, (b) ECMWF April 1998 monthly surface mean temperature 10-years mean. (Sun et al., 2013a, b) [172,173].
Figure 46. (a) PRCM (Dust)–PRCM (without Dust) daily mean surface temperature at 00Z for 8–24 April, (b) ECMWF April 1998 monthly surface mean temperature 10-years mean. (Sun et al., 2013a, b) [172,173].
Atmosphere 14 01324 g046
Figure 47. Simulations Kelvin–Helmholtz instability at t = 320 s from: potential temperature from FB: θFB (multiplication of density δ = 1, Δts = 0.01 s); modified FB: MFB (δ = 16, Δts = 0.04 s); modified FB with smoothing: MFBS (δ = 4, Δts = 0.04 s); modified FB with smoothing: MFBS (δ = 16, Δts = 0.08 s), and horizontal explicit–vertical implicit: HE-VI (Δts = 0.02 s). (Sun et al., 2013c) [197].
Figure 47. Simulations Kelvin–Helmholtz instability at t = 320 s from: potential temperature from FB: θFB (multiplication of density δ = 1, Δts = 0.01 s); modified FB: MFB (δ = 16, Δts = 0.04 s); modified FB with smoothing: MFBS (δ = 4, Δts = 0.04 s); modified FB with smoothing: MFBS (δ = 16, Δts = 0.08 s), and horizontal explicit–vertical implicit: HE-VI (Δts = 0.02 s). (Sun et al., 2013c) [197].
Atmosphere 14 01324 g047
Figure 48. (a) Initial step function, and 800 time step integrations with Co = 0.2 for: (b) Crowley scheme, (c) Gadd scheme, and (d) Sun’s modified scheme (Sun, 1993c) [122].
Figure 48. (a) Initial step function, and 800 time step integrations with Co = 0.2 for: (b) Crowley scheme, (c) Gadd scheme, and (d) Sun’s modified scheme (Sun, 1993c) [122].
Atmosphere 14 01324 g048
Figure 49. A particle O in Euler grids (dashed lines) moves to P in Lagrangian grids (full lines), and the corresponding Eulerian coordinate lines x = xi and y = yj are transformed to the Lagrangian coordinate lines ξij and ηij, respectively. the fluid moves from O i ( ξ , η ) to P i ( X , Y ) (indicated by red circle) within Δt in the regular coordinates. We connect those red circles to form the Lagrangian curve (η = yj), which intercepts with the vertical line X = Xk at blue diamond point Q. (Sun 2023) [200].
Figure 49. A particle O in Euler grids (dashed lines) moves to P in Lagrangian grids (full lines), and the corresponding Eulerian coordinate lines x = xi and y = yj are transformed to the Lagrangian coordinate lines ξij and ηij, respectively. the fluid moves from O i ( ξ , η ) to P i ( X , Y ) (indicated by red circle) within Δt in the regular coordinates. We connect those red circles to form the Lagrangian curve (η = yj), which intercepts with the vertical line X = Xk at blue diamond point Q. (Sun 2023) [200].
Atmosphere 14 01324 g049
Figure 50. The idealized cyclogenesis test on 129 × 129 grids, with δ = 0.05. (a) Numerical simulation after 16 time steps with Co = 4, error = 0.076, (b) analytic solution. (Sun 2023) [200].
Figure 50. The idealized cyclogenesis test on 129 × 129 grids, with δ = 0.05. (a) Numerical simulation after 16 time steps with Co = 4, error = 0.076, (b) analytic solution. (Sun 2023) [200].
Atmosphere 14 01324 g050
Figure 51. Simulated streamline and θ (dashed black lines) at t = 6 h integration for case with width ɑ = 10 km, and height hm = 2 km, mean wind U = 20 m s−1, and temperature lapse rate β = 3.5 K km−1 for z < 12 km, and β = 0.7 K km−1 for z > 12 km at inflow. Vertical velocity is multiple by 10 (Sun and Sun 2015) [60].
Figure 51. Simulated streamline and θ (dashed black lines) at t = 6 h integration for case with width ɑ = 10 km, and height hm = 2 km, mean wind U = 20 m s−1, and temperature lapse rate β = 3.5 K km−1 for z < 12 km, and β = 0.7 K km−1 for z > 12 km at inflow. Vertical velocity is multiple by 10 (Sun and Sun 2015) [60].
Atmosphere 14 01324 g051
Figure 52. (a) Simulated wind U (shaded color), Bernoulli function B (thick white lines), potential temperature θ (dashed black lines), pressure p (thin white lines), and streamlines at t = 6 h integration for case; (b) Same as (a) except that shaded color represents temperature T with different scale in x and z. (Sun and Sun 2015) [60].
Figure 52. (a) Simulated wind U (shaded color), Bernoulli function B (thick white lines), potential temperature θ (dashed black lines), pressure p (thin white lines), and streamlines at t = 6 h integration for case; (b) Same as (a) except that shaded color represents temperature T with different scale in x and z. (Sun and Sun 2015) [60].
Atmosphere 14 01324 g052aAtmosphere 14 01324 g052b
Figure 53. (a) Fix westerly wind profiles at X = 0 (black), and at X = 300 km (red line D) at t = 4 h for Case D: (b) Simulated u (shaded color with red contour and magnitude (m s−1)), streamlines (thin black), and θv (thick black and in (K) at 2 h: (c) same as (b) except at t = 6 h for Case D (Sun 2013) [211].
Figure 53. (a) Fix westerly wind profiles at X = 0 (black), and at X = 300 km (red line D) at t = 4 h for Case D: (b) Simulated u (shaded color with red contour and magnitude (m s−1)), streamlines (thin black), and θv (thick black and in (K) at 2 h: (c) same as (b) except at t = 6 h for Case D (Sun 2013) [211].
Atmosphere 14 01324 g053aAtmosphere 14 01324 g053bAtmosphere 14 01324 g053c
Figure 54. Sounding at RPMT (98646, at Mactan Observation Station) at 00Z18 Jun 2017 [202].
Figure 54. Sounding at RPMT (98646, at Mactan Observation Station) at 00Z18 Jun 2017 [202].
Atmosphere 14 01324 g054
Figure 55. (a) Bernoulli function: B (shaded color), θ (white), w (purple), and streamline at 1.5 h. (b) Buoy (red) and w (black) at x = 335 grid at 1.5 h. (c) B (green) and θ (red) at x = 335 grid at 1.5 h (Sun and Sun 2019) [202].
Figure 55. (a) Bernoulli function: B (shaded color), θ (white), w (purple), and streamline at 1.5 h. (b) Buoy (red) and w (black) at x = 335 grid at 1.5 h. (c) B (green) and θ (red) at x = 335 grid at 1.5 h (Sun and Sun 2019) [202].
Atmosphere 14 01324 g055aAtmosphere 14 01324 g055b
Figure 56. (a): Bernoulli function B (shaded color), ql (thick white), θe (thin white), w (purple), and streamline for Case D at 0.56 h, (b): buoy (shaded color), ql (thick white), w (purple) and streamline at 0.67 h [202], (c): same as (a) but at 0.67 h, (d) θe (black), B (red) at 0.67 h, (e): Buoyancy, buoy (shaded color), ql (thick white), θe (thin white), w (purple), and streamline for Case D. at 1.16 h. (f): B (shaded color), ql (thick white), qr (blue), θe (thin white), and w (black) at 1.16 h (Sun and Sun 2019) [202].
Figure 56. (a): Bernoulli function B (shaded color), ql (thick white), θe (thin white), w (purple), and streamline for Case D at 0.56 h, (b): buoy (shaded color), ql (thick white), w (purple) and streamline at 0.67 h [202], (c): same as (a) but at 0.67 h, (d) θe (black), B (red) at 0.67 h, (e): Buoyancy, buoy (shaded color), ql (thick white), θe (thin white), w (purple), and streamline for Case D. at 1.16 h. (f): B (shaded color), ql (thick white), qr (blue), θe (thin white), and w (black) at 1.16 h (Sun and Sun 2019) [202].
Atmosphere 14 01324 g056aAtmosphere 14 01324 g056b
Figure 57. Accumulated rain (mm) without nonhydrostatic pressure (line R3) and with nonhydrostatic pressure (line R4) (Chen and Sun 2002) [111].
Figure 57. Accumulated rain (mm) without nonhydrostatic pressure (line R3) and with nonhydrostatic pressure (line R4) (Chen and Sun 2002) [111].
Atmosphere 14 01324 g057
Figure 58. Time evolutions of (a) vertical velocity (m s−1), (c) potential temperature anomaly (K), and (e) moisture anomaly (kg kg−1) from the one-dimensional cloud model, and the averaged (b) vertical velocity (m s−1), (d) potential temperature anomaly (K), and (f) moisture anomaly (kg kg−1) within the radius of 5000 m cloud from the WRF model. The values in (e,f) are multiplied by 103 (Chen and Sun 2002), [111].
Figure 58. Time evolutions of (a) vertical velocity (m s−1), (c) potential temperature anomaly (K), and (e) moisture anomaly (kg kg−1) from the one-dimensional cloud model, and the averaged (b) vertical velocity (m s−1), (d) potential temperature anomaly (K), and (f) moisture anomaly (kg kg−1) within the radius of 5000 m cloud from the WRF model. The values in (e,f) are multiplied by 103 (Chen and Sun 2002), [111].
Atmosphere 14 01324 g058
Figure 59. The NTU/PU simulation of January 25 downslope windstorm case at WSMR over a no-slip surface at 32 N after 4 h integration with a mountain height of 1.5 km. The vertical grid interval is 300 m and horizontal grid spacing is 1 km. The model’s 10 m winds are shown by the wind barbs. Observed 10 m winds are shown by the thick arrows (Haines et al., 2019) [220].
Figure 59. The NTU/PU simulation of January 25 downslope windstorm case at WSMR over a no-slip surface at 32 N after 4 h integration with a mountain height of 1.5 km. The vertical grid interval is 300 m and horizontal grid spacing is 1 km. The model’s 10 m winds are shown by the wind barbs. Observed 10 m winds are shown by the thick arrows (Haines et al., 2019) [220].
Atmosphere 14 01324 g059
Figure 60. Same as Figure 59, except for WRF (Haines et al., 2006 [222]).
Figure 60. Same as Figure 59, except for WRF (Haines et al., 2006 [222]).
Atmosphere 14 01324 g060
Figure 61. NTU/PU simulated wind vector at z = 10 m above ground after 4 h integration from a constant westerly wind of 5 m s−1 with Δx = Δy = 2 km and Δz = 300 m. The observed wind at 0800L25 25 January 2004, is indicated by yellow arrows. The arrow at right-lower corner indicates 20 m s−1. The elevation is shaded by color. Yellow arrows are the observed wind at the time (Haines et al., 2019) [220].
Figure 61. NTU/PU simulated wind vector at z = 10 m above ground after 4 h integration from a constant westerly wind of 5 m s−1 with Δx = Δy = 2 km and Δz = 300 m. The observed wind at 0800L25 25 January 2004, is indicated by yellow arrows. The arrow at right-lower corner indicates 20 m s−1. The elevation is shaded by color. Yellow arrows are the observed wind at the time (Haines et al., 2019) [220].
Atmosphere 14 01324 g061
Figure 62. Contravariant (tangent) basis vectors g i , covector (dual) basis g i , and a 2D vector r = x ^ i g i ( = 1 2 x ^ i g i ) = x ^ i g i ( = 1 2 x ^ i g i ) (Sun 2021) [223].
Figure 62. Contravariant (tangent) basis vectors g i , covector (dual) basis g i , and a 2D vector r = x ^ i g i ( = 1 2 x ^ i g i ) = x ^ i g i ( = 1 2 x ^ i g i ) (Sun 2021) [223].
Atmosphere 14 01324 g062
Figure 63. 3D curvilinear coordinates and volume dv = d x ^ 1 g 1 × d x ^ 2 g 2 d x ^ 3 g 3 (Sun 2021) [223].
Figure 63. 3D curvilinear coordinates and volume dv = d x ^ 1 g 1 × d x ^ 2 g 2 d x ^ 3 g 3 (Sun 2021) [223].
Atmosphere 14 01324 g063
Figure 64. Basis vectors of 2D Cartesian coordinate and terrain following coordinates. (Sun 2021) [223].
Figure 64. Basis vectors of 2D Cartesian coordinate and terrain following coordinates. (Sun 2021) [223].
Atmosphere 14 01324 g064
Figure 65. Function S in (240) is shown by yellow lines, ∇Sa by color vector; x and y components of gradient in Cartesian coordinate by black lines, which are covered by dashed red lines of ∇Sb and blue dots of ∇Sc in terrain following coordinates. Green lines show terrain and depth beneath (in meters.) (Sun 2021) [223].
Figure 65. Function S in (240) is shown by yellow lines, ∇Sa by color vector; x and y components of gradient in Cartesian coordinate by black lines, which are covered by dashed red lines of ∇Sb and blue dots of ∇Sc in terrain following coordinates. Green lines show terrain and depth beneath (in meters.) (Sun 2021) [223].
Atmosphere 14 01324 g065
Figure 66. Black lines for divergence of V of (244) calculated from Diva, dashed red lines from Divb, and blue dots from Divc. (Sun 2021) [223].
Figure 66. Black lines for divergence of V of (244) calculated from Diva, dashed red lines from Divb, and blue dots from Divc. (Sun 2021) [223].
Atmosphere 14 01324 g066
Figure 67. y-component curl: Black lines of V of (248) calculated from Curla, dashed red lines from Curlb and blue dots from Curlc, and green lines for terrain. (Sun 2021) [223].
Figure 67. y-component curl: Black lines of V of (248) calculated from Curla, dashed red lines from Curlb and blue dots from Curlc, and green lines for terrain. (Sun 2021) [223].
Atmosphere 14 01324 g067
Figure 68. Local time derivative of u: black lines calculated from Cartesian coordinate, red from (253) based on new coordinate x ^ , y ^ , z ^ , and blue dots based on Gal-Chen and Somerville coordinate x ¯ , y ¯ , z ¯ . (Sun 2021) [223].
Figure 68. Local time derivative of u: black lines calculated from Cartesian coordinate, red from (253) based on new coordinate x ^ , y ^ , z ^ , and blue dots based on Gal-Chen and Somerville coordinate x ¯ , y ¯ , z ¯ . (Sun 2021) [223].
Atmosphere 14 01324 g068
Figure 69. Simulated non-dimensional crosswind integrated concentration C y * as a function of non-dimensional distance, x * = [ w * / ( U m Z i ) ] C y * = ( U m Z i / Q ) C y measured form point sources at (a) zs = 0.75 Zi; (b) zs = 0.49 Zi; (c) zs = 0.24 Zi; and (d) zs = 0.067 Zi (Zi is depth of mixed layer, Um is mean wind, w* is convective vertical velocity at surface) (Sun 1988, 1989) [218,219].
Figure 69. Simulated non-dimensional crosswind integrated concentration C y * as a function of non-dimensional distance, x * = [ w * / ( U m Z i ) ] C y * = ( U m Z i / Q ) C y measured form point sources at (a) zs = 0.75 Zi; (b) zs = 0.49 Zi; (c) zs = 0.24 Zi; and (d) zs = 0.067 Zi (Zi is depth of mixed layer, Um is mean wind, w* is convective vertical velocity at surface) (Sun 1988, 1989) [218,219].
Atmosphere 14 01324 g069
Figure 70. Crosswind integrated concentration C y * vs. x* from laboratory experiments for point sources of zs = 0.49, 0.24, and 0.067 Zi (Deardorff and Willis 1975, Willis and Deardorff 1978, 1981) [233,234,235].
Figure 70. Crosswind integrated concentration C y * vs. x* from laboratory experiments for point sources of zs = 0.49, 0.24, and 0.067 Zi (Deardorff and Willis 1975, Willis and Deardorff 1978, 1981) [233,234,235].
Atmosphere 14 01324 g070
Figure 71. The size and arrangement of cylinders. (Sun and Sun 2017) [236].
Figure 71. The size and arrangement of cylinders. (Sun and Sun 2017) [236].
Atmosphere 14 01324 g071
Figure 72. Different species (black, red, green, blue, and light blue for specie 1, 2, 3, 4, and 5, respectively) at each cell (A, B, C, D, E): (a) ρ 1 A ,   ρ 2 A , ρ 3 A ,   ρ 4 A , and ρ 5 A ; (b) ρ 1 B ,   ρ 2 B ,   ρ 3 B ,   ρ 4 B , and ρ 5 B ; (c) ρ 1 C ,   ρ 2 C ,   ρ 3 C ,   ρ 4 C , and ρ 5 C ; (d) ρ 1 D ,   ρ 2 D ,   ρ 3 D ,   ρ 4 D , and ρ 5 D ; and (e) ρ 1 E ,   ρ 2 E ,   ρ 3 E ,   ρ 4 E , and ρ 5 E (Sun and Sun 2017) [236].
Figure 72. Different species (black, red, green, blue, and light blue for specie 1, 2, 3, 4, and 5, respectively) at each cell (A, B, C, D, E): (a) ρ 1 A ,   ρ 2 A , ρ 3 A ,   ρ 4 A , and ρ 5 A ; (b) ρ 1 B ,   ρ 2 B ,   ρ 3 B ,   ρ 4 B , and ρ 5 B ; (c) ρ 1 C ,   ρ 2 C ,   ρ 3 C ,   ρ 4 C , and ρ 5 C ; (d) ρ 1 D ,   ρ 2 D ,   ρ 3 D ,   ρ 4 D , and ρ 5 D ; and (e) ρ 1 E ,   ρ 2 E ,   ρ 3 E ,   ρ 4 E , and ρ 5 E (Sun and Sun 2017) [236].
Atmosphere 14 01324 g072aAtmosphere 14 01324 g072bAtmosphere 14 01324 g072c
Figure 73. Forward-integrations (+) and backward integration ⊗ for species (1, 2, 3, 4, 5) among different cells (A, B, C, D. E), (a) ρ 1 A ,   ρ 1 B ,   ρ 1 C ,   ρ 1 D , and ρ 1 E ; (b) ρ 2 A ,   ρ 2 B ,   ρ 2 C ,   ρ 2 D , and ρ 2 E ; (c) ρ 3 A ,   ρ 3 B ,   ρ 3 C ,   ρ 3 D , and ρ 3 E ; (d) ρ 4 A ,   ρ 4 B ,   ρ 4 C ,   ρ 4 D , and ρ 4 E ; and (e) ρ 5 A ,   ρ 5 B ,   ρ 5 C ,   ρ 5 D , and ρ 5 E (Sun and Sun 2017) [236].
Figure 73. Forward-integrations (+) and backward integration ⊗ for species (1, 2, 3, 4, 5) among different cells (A, B, C, D. E), (a) ρ 1 A ,   ρ 1 B ,   ρ 1 C ,   ρ 1 D , and ρ 1 E ; (b) ρ 2 A ,   ρ 2 B ,   ρ 2 C ,   ρ 2 D , and ρ 2 E ; (c) ρ 3 A ,   ρ 3 B ,   ρ 3 C ,   ρ 3 D , and ρ 3 E ; (d) ρ 4 A ,   ρ 4 B ,   ρ 4 C ,   ρ 4 D , and ρ 4 E ; and (e) ρ 5 A ,   ρ 5 B ,   ρ 5 C ,   ρ 5 D , and ρ 5 E (Sun and Sun 2017) [236].
Atmosphere 14 01324 g073aAtmosphere 14 01324 g073bAtmosphere 14 01324 g073c
Figure 74. Numerical simulations from the NTU/Purdue model and WRF over an idealized elliptical mountain: They contain 600 × 600 × 50 in x, y, and z directions with space interval Δx = Δy = 1 km; Δz~300 m on a free-slip surface; initial wind U = 4 m s−1 and the Brunt–Vaisala frequency N = 10−2 s−1, mountain peak = 2 km; half width lengths: 5 and 10 km (tilted by 30 degrees), and the Froude number Fr = U/Nh = 0.2, (a) after 10 h integration, (b) after 20 h integration (Sun et al., 2010) [242].
Figure 74. Numerical simulations from the NTU/Purdue model and WRF over an idealized elliptical mountain: They contain 600 × 600 × 50 in x, y, and z directions with space interval Δx = Δy = 1 km; Δz~300 m on a free-slip surface; initial wind U = 4 m s−1 and the Brunt–Vaisala frequency N = 10−2 s−1, mountain peak = 2 km; half width lengths: 5 and 10 km (tilted by 30 degrees), and the Froude number Fr = U/Nh = 0.2, (a) after 10 h integration, (b) after 20 h integration (Sun et al., 2010) [242].
Atmosphere 14 01324 g074aAtmosphere 14 01324 g074b
Table 1. Ratio of computed phase speed to analytic phase speed C p h / C as function of wavelength (λ) and Courant number Co for 2nd-order scheme according to (24) (Sun and Sun 2011) [37].
Table 1. Ratio of computed phase speed to analytic phase speed C p h / C as function of wavelength (λ) and Courant number Co for 2nd-order scheme according to (24) (Sun and Sun 2011) [37].
λ2Δx3Δx4Δx5Δx6Δx7Δx8Δx9Δx10Δx
Co
0.2500.6430.8340.9050.9390.9570.9690.9760.9810.985
0.5000.6670.8550.9200.9500.9650.9750.9810.9850.988
0.7500.6130.9000.9490.9690.9790.9850.9880.9910.993
1.000.0000.5001.001.001.001.001.001.001.00
Table 2. Parameters in atmosphere and in shallow-water model with (a) R = 24 km and U = 0.8 ms−1, and (b) R = 40 km and U = 1.2 m s−1, (b) Same as Table 2a except R = 40 km and U = 1.2 m s−1. (Sun 2016) [55].
Table 2. Parameters in atmosphere and in shallow-water model with (a) R = 24 km and U = 0.8 ms−1, and (b) R = 40 km and U = 1.2 m s−1, (b) Same as Table 2a except R = 40 km and U = 1.2 m s−1. (Sun 2016) [55].
(a)R (km)H (m)f
(10−5 s−1)
hm(m)vel (m/s)N
or. g′
RoLRFrmFr
AtmR* = 120H* = 1045.783 h m * = 3500U* = 4N = 10−2 s−1U*/f*R* = 0. 58NH*/(f*R*) = 102/(f*R*) = 14.41U*/(Nhm) = U/(35 ms−1) = 0.11U*/(NH*) = 0.04
WaterR = 24H = 1005.783hm = 10U = 0.8g′ = 4 m s−2U/fR = 0.58 g H / ( fR ) = 20/(fR) = 14.41 U / g h m = 0.2 U*/6.3 = U/(32 ms−1) = 0.12 U / g H = 0.04
(b)R (km)H (m)f
(10−5 s−1)
hm(m)vel (m/s)N
or. g′
RoLRFrmFr
AtmR* = 2001045.783 h m * = 3500U* = 6N = 10−2 s−1U*/f*R* = 0.52NH*/(f*R*) = 8.65U*/(Nhm) = 0.17U*/(NH*) = 0.06
WaterR = 40H* = 1005.783hm = 10U = 1.2g′ = 4 m s−2U/fR = 0.53 g H / ( fR ) = 8.65 U / g h m = 0.19 U / g H = 0.06
Table 3. Climatological Pattern Correlation of PRCM (Yu et al., 1994b) [162,163].
Table 3. Climatological Pattern Correlation of PRCM (Yu et al., 1994b) [162,163].
Climatological Pattern Correlation of PRCM
MJJAMayJunJulAug
MSLP0.980.960.980.970.96
TEM2M0.990.990.990.970.97
850T0.870.970.900.730.59
200T0.880.830.830.940.94
850Z0.970.960.990.980.97
200Z0.960.990.980.940.77
850Q0.960.980.960.950.91
200Q0.840.930.800.700.71
PRECIP0.350.370.520.270.29
PRECIP
(bias removed) 0.720.890.900.85
Table 4. 20-day mean statistics pf PRCM compared to reanalysis data of ECMWF except for precipitation which is derived from GPCP.
Table 4. 20-day mean statistics pf PRCM compared to reanalysis data of ECMWF except for precipitation which is derived from GPCP.
FIELDLEVELBIASRMSCOR
Mean Sea Level Pressure (hPa) Air Qv (kg/kg) Temperature (K)
Wind speed (ms−1)
Precipitation (mmday−1)
SFC−0.481.570.95
2.45 × 10−32.81 × 10−30.96
0.731.750.97
0.781.260.84
0.951.380.73
Height (m)200 hPa3745.10.99
500 hPa12.517.70.98
850 hPa0.2510.30.97
Temperature (K)200 hPa0.741.170.9
500 hPa0.630.860.98
850 hPa0.131.220.94
Qv (kg/kg)500 hPa−1.34 × 10−46.59 × 10−40.83
700 hPa3.01 × 10−41.22 × 10−30.82
850 hPa1.13 × 10−31.61 × 10−30.92
Wind speed (ms−1)200 hPa−1.603.980.91
500 hPa−0.722.140.8
850 hPa1.282.160.83
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Sun, W.-Y. Challenges and Progress in Computational Geophysical Fluid Dynamics in Recent Decades. Atmosphere 2023, 14, 1324. https://doi.org/10.3390/atmos14091324

AMA Style

Sun W-Y. Challenges and Progress in Computational Geophysical Fluid Dynamics in Recent Decades. Atmosphere. 2023; 14(9):1324. https://doi.org/10.3390/atmos14091324

Chicago/Turabian Style

Sun, Wen-Yih. 2023. "Challenges and Progress in Computational Geophysical Fluid Dynamics in Recent Decades" Atmosphere 14, no. 9: 1324. https://doi.org/10.3390/atmos14091324

APA Style

Sun, W. -Y. (2023). Challenges and Progress in Computational Geophysical Fluid Dynamics in Recent Decades. Atmosphere, 14(9), 1324. https://doi.org/10.3390/atmos14091324

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop