The magazine publisher is the copyright holder of this article and it is reproduced with permission. Further reproduction of this article in violation of the copyright is prohibited. ABSTRACT 1. INTRODUCTION 2. MODEL FORMULATION A. PRIMITIVE EQUATION FORMULATION Equations (1)-(4) form a system of six equations for the following explicitly calculated parameters: u, v, w, p[subM], theta[subM], and q[subM]. Other unknowns, such as turbulent flux terms (introduced by Reynolds averaging), and the diabatic heat source term Q in the energy equation, must be parameterized. B. VORTICITY-STREAMFUNCTION FORMULATION Prognostic vorticity equations are obtained by application of the curl operator to (1), to produce which by definition satisfies [universal quantifier] • rho[sub0]V = 0. Vorticity is then related to Psi* by with arbitrary scalar function lambda. These transformations leave the velocity field unchanged and can be used to specify special conditions on Psi. We therefore choose lambda so that the third streamfunction component vanishes; that is, By substitution of (10) in (11), the following condition on lambda must be fulfilled: The resulting stream vector Psi = (Psi[sub1], Psi[sub2], 0) fulfills both rho[sub0]V = [universal quantifier] × Psi and (11) and may be used without loss of generality. The major advantage of this choice of lambda is that only two streamfunction components are necessary to totally describe an arbitrary three-dimensional flow filed. In consequence, only the following two components (horizontal ones) of (9) need to be solved: C. COORDINATE TRANSFORMATION D. THERMAL ENERGY, SPECIFIC HUMIDITY, AND TKE EQUATIONS E. BOUNDARY CONDITIONS 1) SURFACE BOUNDARY CONDITIONS 3) For thermally forced circulations, temperature and specific humidity in TVM are calculated either by a modified force restore model of Deardorff (Schayes et al. 1996) or by the Institut d' Astronomie et de Géophysique Georges Lema^itre land surface model (De Ridder and Schayes 1997). For dynamically forced circulations, temperature and specific humidity are fixed. 2) UPPER BOUNDARY CONDITIONS 3) LATERAL BOUNDARY CONDITIONS F. NUMERICS 3. EVALUATION OF THE NONHYDROSTATIC ANELASTIC VORTICITY FORMULATION A. HYDROSTATIC LINEAR MOUNTAIN WAVES Solutions are governed by two dimensionless quantities. The first, F[subx] = al, is the ratio between the horizontal obstacle length scale a and the intrinsic vertical flow length scale (inverse Scorer parameter) l[sup-1], and the second, F[subz] = h[subm]l, is the ratio between the vertical obstacle length scale h[subm] and l[sup-1]. While F[subx] defines the degree of nonhydrostaticity of the flow, F[subz] gives information about its nonlinearity. If F[subx] >> 1, the flow is virtually hydrostatic, while nonhydrostatic effects become important when F[subx] ˜ 1 (Schumann et al. 1987). If F[subz] << 1, the flow is essentially linear, while nonlinear effects become important when F[subz] ˜ 1. B. NONLINEAR MOUNTAIN WAVES C. THE BOULDER WINDSTORM 4. EVALUATION OF THE HYDROSTATIC AND INCOMPRESSIBLE APPROXIMATIONS FOR THERMALLY INDUCED CIRCULATIONS A. SEA BREEZES CONVERGENCE OVER AN ISLAND in which alpha = 8 and tau[subd] = 86 400 s. Water temperature is kept constant at theta[sub0] during the entire simulation. All parameters are kept similar to those of Yang, that is, an horizontal grid resolution of 1 km, a potential temperature gradient of 4 K km[sup-1], and a zero geostrophic wind, except for the vertical grid, which is stretched differently in this study. Note that ensuring an identical grid configuration will not ensure better comparisons since implicit and explicit diffusion (both horizontally and vertically) could be significantly different and modify the effective grid resolution used in the two studies. B. SEA BREEZES OVER COMPLEX TERRAIN 5. CONCLUSIONS FIG. 1. (a), (b) Analytical and (c), (d) TVM nonhydrostatic solutions for vertical wind velocity (cm s[sup-1]) for the hydrostatic linear mountain wave (MW) case. The left (a)-(c) and right (b)-(d) columns illustrate the incompressible and anelastic solutions, respectively. REFERENCES APPENDIX LIST OF SYMBOLS A. ROMAN B. GREEK C. SUBSCRIPTS D. SPECIALAUTHOR: P. THUNIS; A. CLAPPIER TITLE: Formulation and Evaluation of a Nonhydrostatic Mesoscale Vorticity Model (TVM) SOURCE: Monthly Weather Review 128 no9 3236-51 S 2000
This paper describes the formulation and the application of the nonhydrostatic anelastic vorticity model (TVM). This model is constructed using a method involving two horizontal streamfunctions and two horizontal vorticity components. The evaluation of this formulation is performed by simulating various bidimensional hydrostatic and nonhydrostatic mountain wave cases. Results are compared with analytical solutions and in particular with those developed by Laprise and Peltier for nonlinear forcings. The anelastic formulation is also validated with respect to the highly nonlinear 1972 Boulder windstorm. TVM is shown to accurately reproduce these mountain wave test cases in both its incompressible and anelastic formulations. In the second part of this paper, the adequacy of the hydrostatic and anelastic assumptions in simulating thermally induced circulations is investigated and compared to previous works. For a moderate surface thermal forcing, typical geographical setups are used and show slight differences between hydrostatic and nonhydrostatic horizontal wind speeds. For vertical wind speeds, differences are shown to be much larger and more sensitive to changes in grid resolution. For more stringent thermal surface forcing, differences remain low for horizontal wind speeds but increase considerably for the vertical wind component. The comparison between anelastic and incompressible solutions for the same cases shows the adequacy of the incompressible assumption when circulations are forced by the surface and are characterized by a relatively shallow vertical extent. In such conditions, virtually no differences are observed between the two formulations.
Although the vorticity approach in fluid dynamics modeling is not new, it has not been frequently used in mesoscale atmospheric flow applications. One horizontal vorticity component, however, has been used to simulate circulations in vertical cross sections or for cumulus cloud simulations (Orville and Sloan 1970; Pearson 1973; Moeng and Arakawa 1980).
The first attempt to incorporate all three space dimensions into a vorticity model was by Thyer (1966). Following ideas put forth by Malkus and Witt (1959), he successfully used a nonhydrostatic (NH) incompressible model with the two horizontal vorticity and two streamfunction components to reproduce mountain-valley wind circulations. Due to contemporary computer capabilities and to numerical scheme instabilities, however, the equations were successfully integrated for only 120 s.
This approach was also followed by Bornstein et al. (1986), Sievers and Zdunkowski (1986), and Schayes et al. (1996). While the second authors applied their nonhydrostatic incompressible model mainly to urban canopy layer flows, the other two developed a hydrostatic (H) and incompressible formulation for mesoscale flows (e.g., sea breezes and mountain winds). A different approach, which retains all three vorticity and streamfunction components, was used by Saitoh et al. (1996) to resolve nonhydrostatic urban circulations.
One of the characteristics of the vorticity approach is that both H and NH formulations are based on the same set of prognostic conservation equations; that is, the number of equations (prognostic or diagnostic) and numerical algorithms used to solve them are strictly identical. This is not strictly true for primitive-equation formulations for which the number of prognostic momentum equations differs. Indeed, an additional prognostic equation is required in the NH formulation for the vertical momentum component. Resolution of this additional equation may add numerical bias to the physical H-NH differences, preventing strict isolation of the impact of the hydrostatic assumption. This property is used in the second part of this work to carry out an intercomparison of H and NH solutions for different atmospheric circulations.
Section 2 describes the formulation of a nonhydrostatic anelastic version of the hydrostatic incompressible Thermal Vorticity Model (TVM) of Schayes et al. (1996) and Bornstein et al. (1996) with a special focus on the new developments (i.e., anelasticity and nonhydrostaticity). Section 3 is devoted to the validation of this new formulation. This is realized by simulating mountain wave cases with different degrees of nonlinearity, nonhydrostaticity, and anelasticity using analytical solutions as reference. The 1972 Boulder windstorm is also simulated and comparisons are made to both measurements and previous model results. In section 4 comparisons of H versus NH simulations of thermally induced circulations are presented and compared to previous studies (e.g., Yang 1991). Finally, conclusions are presented in section 5.
One of the main characteristics of the vorticity formulation is the total elimination of pressure from the prognostic equations (Pielke 1984). As vorticity equations are constructed by taking the curl of the momentum equation, the pressure gradient term drops by definition, but mesoscale pressure remains in the buoyancy term and in the thermal energy equation. Total elimination of mesoscale pressure from the prognostic vorticity equations implies one supplementary assumption, namely that thermal forcing is dominant in the atmospheric flows of interest. Indeed, for such cases (subclass of deep Boussinesq flows), referred to as deep thermal convection (Thunis and Bornstein 1996), it may be shown that mesoscale pressure contributions can be neglected in all terms.
The corresponding set of primitive equations is thus selected as the starting point for derivation of the TVM vorticity equations. The Reynolds-averaged equations (within a rectangular coordinate system) for momentum, thermal energy, specific humidity, and mass are expressed, respectively, as follows:
[df]rho[sub0]V / [df]t = -[universal quantifier] • (rho[sub0]VV) - [universal quantifier] • (rho[sub0]V'V') - [universal quantifier][subH]p[sub0] - [universal quantifier]p[subM] - theta[subM] / theta[sub0] rho[sub0]g - 2Omega X rho[sub0]V, (1)
[df]theta[subM] / [df]t = -1 / rho[sub0][universal quantifier] • (Vrho[sub0]theta[subM]) - 1 / rho[sub0][universal quantifier] • (rho[sub0]V'theta') - 1 / rho[sub0] [df]wrho[sub0]theta[sub0] / [df]z - (1 / rho[sub0]c[subp] [df]Q / [df]z), (2)
[df]q[subM] / [df]t = -1 / rho[sub0][universal quantifier] • (Vrho[sub0]q[subM] - 1 / rho[sub0][universal quantifier] • (rho[sub0]V'q') - 1 / rho[sub0] [df]wrho[sub0]q[sub0] / [df]z, and (3)
[universal quantifier] • rho[sub0]V = 0, (4)
where Reynolds-averaged quantities are unbarred for convenience and where all symbols are defined in the appendix. As seen from (3), this version of TVM does not include water vapor source/sink terms (no phase change). Also, basic-state thermodynamic variables (phi[sub0] = {Begin Greek}f - f{End Greek}[subM]) are assumed stationary and in hydrostatic balance; that is,[df]p[sub0] / [df]z = -rho[sub0]g.
For deep Boussinesq flows, the vorticity vector zeta is defined aszeta = [universal quantifier] X rho[sub0]V. (5)
[df]zeta / [df]t = -[universal quantifier] × ([universal quantifier] • rho[sub0]VV) - [universal quantifier] × [[universal quantifier] • rho[sub0](V'V')] - [universal quantifier] × theta[subM] / theta[sub0]rho[sub0]g - 2[universal quantifier] × (Omega × rho[sub0]V), (6)
in which the pressure gradient term disappears by virtue of the vector identity [universal quantifier] × [universal quantifier]p = 0, and where Omega = (0, Omega cos{Begin Greek}f, Q {End Greek}sinphi) = (0, f*/2, f/2). Equation (6) may be rearranged as
[df]zeta / [df]t = -[universal quantifier] • (Vzeta) - [universal quantifier] • [rho[sub0] V([universal quantifier] × V)] - [universal quantifier] × [[universal quantifier] • (rho[sub0] V'V')] - [universal quantifier] × theta[subM] / theta[sub0 rho[sub0]g + 2(omega • [universal quantifier])rho[sub0]V. (7)
To recover velocities from the vorticities, a streamfunction Psi* is defined asrho[sub0]V = [universal quantifier] X Psi*, (8)
zeta = [universal quantifier] × rho[sub0]V = [universal quantifier] × [universal quantifier] × Psi* = -[universal quantifier][sup2]Psi* + [universal quantifier]([universal quantifier] • Psi*). (9)
As shown by Sievers and Zdunkowski (1986), Psi* in (8) is not uniquely determined, since a gauge transformation exists, that is,Psi --> Psi* + [universal quantifier]lambda (10)
Psi • k = 0. (11)
lambda = [integral][supz[sub[sub0] Psi* • k dz'.
omega[sub1] = [Graphic Character Omitted] + [Graphic Character Omitted] and omega[sub2] = [Graphic Character Omitted] + [Graphic Character Omitted], (12)
where for notational convenience, new vorticity and streamfunction fields were defined as follows:
omega = (-zeta[sub1], zeta[sub2], zeta[sub3]) and F = ({Begin Greek}W, -W{End Greek}[sub2], 0).
Note, however that each horizontal vorticity in (12) is a function of both streamfunction components. A further advantage of the gauge transformation (11) is that only the two horizontal vorticity components must be numerically integrated:
[df]omega[sub1] / [df]t = -[df]uomega[sub1] / [df]x - [df]vomega[sub1] / [df]y - [df]womega[sub1] / [df]z - [df] / [df]x(rho[sub0]v[df]u / [df]z - rho[sub0]w[df]u / [df]y) - [df] / [df]y(rho[sub0]v[df]v / [df]z - rho[sub0]w[df]v / [df]y) - [df] / [df]z(rho[sub0]v[df]w / [df]z - rho[sub0]w[df]w / [df]y) + K[subH][df][sup2]omega[sub1] / [df]x[sup2] + K[subH][df][sup2]omega[sub1] / [df]y[sup2] + [df][sup2] / [df]z[sup2](K[subm]omega[sub1]) - g / theta[sub0] [df](rho[sub0]theta[subM]) / [df]y - f[df](rho[sub0]u) / [df]z - f*[df](rho[sub0]u) / [df]y and
[df]omega[sub2] / [df]t = -[df]uomega[sup2] / [df]x - [df]vomega[sup2] / [df]y - [df]womega[sup2] / [df]z - [df] / [df]x(rho[sub0]u[df]u / [df]z - rho[sub0]w[df]u / [df]x) -[df] / [df]y(rho[sub0]u[df]v / [df]z - rho[sub0]w[df]v / [df]x) - [df] / [df]z(rho[sub0]u[df]w / [df]z - rho[sub0]w[df]w / [df]x) +K[subH][df][sup2]omega[sub2] / [df]x[sup2] + K[subH][df][sup2]omega[sub2] / [df]y[sup2] + [df][sup2] / [df]z[sup2](K[subm]omega[sub2]) - g / theta[sub0] [df](rho[sub0]theta[subM]) / [df]x + f[df](rho[sub0]v) / [df]z - f*[df](rho[sub0]u) / [df]x - f*[df](rho[sub0]w) / [df]z. (13)
The H and NH contributions in (12) were underlined whereas (13) remains unchanged for both formulations. The horizontal diffusion coefficient K[subH] was assumed constant in space, an approximation commonly made in Reynolds-averaged mesoscale models since it normally serves only as a numerical filter.
In mesoscale atmospheric applications, it has become standard practice to transform only the vertical coordinate. The functional form of the resulting generalized contravariant vertical coordinate system in terms of the original Cartesian system is
~x[sup1] = x[sup1], ~x[sup2] = x[sup2], and ~x[sup3] = sigma = z[subt]z - z[subg](x, y) / z[subt] - z[subg](x, y), (14)
in which z[subt] and z[subg] are the model top and surface orography, respectively. This transformation, first derived by Gal-Chen and Somerville (1975), has the advantages of invariance with time and application of lower boundary conditions at the first sigma level. Equations (12) and (13) become, after regrouping terms,
[df]õmega[sup1] / [df]t = -[df]~u[sup1]õmega[sup1] / [df]~x - [df]~u[sup2]õmega[sup1] / [df]~y - [df]~u[sup3]õmega[sup1] / [df]sigma - [df] / [df]~x(rho[sub0]~u[sub2]~u[sup1[sub[sub;3] - rho[sub0]~u[sub3]~u[sup1[sub[sub;2]) - [df] / [df]~y(rho[sub0]~u[sub2]~u[sup2[sub[sub;3] - rho[sub0]~u[sub3]~u[sup2[sub[sub;2]) - [df] / [df]{Begin Greek}s(r{End Greek}[sub0]~u[sub2]~u[sup3[sub[sub;3] - rho[sub0]~u[sub3]~u[sup3[sub[sub;2]) + f*[df] / [df]sigma([df]z / [df]~yrho[sub0]~u[sup1]) - f[df]rho[sub0]~u[sup1] / [df]sigma - f*[df] / [df]~y([df]z / [df]{Begin Greek}sr{End Greek}[sub0]~u[sup1]) + g[df] / [df]{Begin Greek}s(j{End Greek}[subM] / theta[sub0] [df]z / [df]~yrho[sub0]) - g[df] / [df]~y(theta[subM] / theta[sub0] [df]z / [df]{Begin Greek}sr{End Greek}[sub0]) + [df][sup2](K[subH]õmega[sup1]) / [df]~x[sup2] + [df][sup2](K[subH]õmega[sup1]) / [df]~y[sup2] + [df][sup2](K[subm]õmega[sup1]) / [df]sigma[sup2] - [df][sup2] / [df]sigma[sup2](rho[sub0]~u[sub3] [df]K[subm] / [df]~y) + [df][sup2] / [df]sigma[df]~y(rho[sub0]~u[sub3][df]K[subm] / [df]{Begin Greek}s) - G{End Greek}[sup3[sub[sub23][df] / [df]sigma(K[subm][df]rho[sub0]~u[sub3] / [df]sigma) and [df]õmega[sup2] / [df]t = - [df]~u[sup1]õmega[sup2] / [df]~x - [df]~u[sup2]õmega[sup2] / [df]~y - [df]~u[sup3]õmega[sup2] / [df]sigma - [df] / [df]~x(rho[sub0]~u[sub1]~u[sup1[sub[sub;3] - rho[sub0]~u[sub3]~u[sup1[sub[sub;1]) - [df] / [df]~y(rho[sub0]~u[sub1]~u[sup2[sub[sub;3] - rho[sub0]~u[sub3]~u[sup2[sub[sub;1]) - [df] / [df]{Begin Greek}s(r{End Greek}[sub0]~u[sub1]~u[sup3[sub[sub;3] - rho[sub0]~u[sub3]~u[sup3[sub[sub;1]) - [df]z / [df]sigmaf*[df]rho[sub0]~u[sup3] / [df]sigma - f*[df] / [df]sigma([df]z / [df]~yrho[sub0]~u[sup2]) + f[df]rho[sub0]~u[sup2] / [df]sigma - f*[df] / [df]~x([df]z / [df]{Begin Greek}sr{End Greek}[sub0]~u[sup1]) + g[df] / [df]{Begin Greek}s(j{End Greek}[subM] / theta[sub0] [df]z / [df]~xrho[sub0]) - g[df] / [df]~x(theta[subM] / theta[sub0] [df]z / [df]{Begin Greek}sr{End Greek}[sub0]) + [df][sup2](K[subH]õmega[sup2]) / [df]~x[sup2] + [df][sup2](K[subH]õmega[sup2]) / [df]~y[sup2] + [df][sup2](K[subm]õmega[sup2]) / [df]sigma[sup2] - [df][sup2] / [df]sigma[sup2](rho[sub0]~u[sub3][df]K[subm] / [df]~x) + [df][sup2] / [df]sigma[df]~x(rho[sub0]~u[sub3][df]K[subm] / [df]{Begin Greek}s) - G{End Greek}[sup3[sub[sub13][df] / [df]sigma(K[subm][df]rho[sub0]~u[sub3] / [df]sigma),
in which the only approximation made has consisted of neglecting topographical second-order derivatives in estimating the turbulent diffusion contribution. Elliptic equations (12) linking streamfunctions to vorticities become
[df][sup2]~Phi[sub1] / [df]~y[sup2] + [([df]sigma / [df]z)[sup2] + ([df]sigma / [df]y)[sup2]] [df][sup2]~Phi[sub1] / [df]sigma[sup2] + 2[df]sigma / [df]y [df][sup2]~Phi[sub1] / [df]sigma[df]y + [df]~Phi[sub1] / [df]sigma[[df] / [df]sigma([df]sigma / [df]y)[sup2] - [df]sigma / [df]z [df] / [df]~y([df]z / [df]~y)] = [df]sigma / [df]zõmega[sup1] - [df]sigma / [df]y [df]sigma / [df]x [df][sup2]~Phi[sub2] / [df]sigma[sup2] - [df][sup2]~Phi[sub2] / [df]~y[df]~x - [df]sigma / [df]y [df][sup2]~Phi[sub2] / [df]sigma[df]~x - [df]sigma / [df]x [df][sup2]~Phi[sub2] / [df]sigma[df]~y - [df]~Phi[sub2] / [df]sigma[[df] / [df]sigma([df]sigma / [df]y [df]sigma / [df]x) - [df]sigma / [df]z [df] / [df]~y([df]z / [df]~x)] and (15) [df][sup2]~Phi[sub2] / [df]~x[sup2] + [([df]sigma / [df]z)[sup2] + ([df]sigma / [df]x)[sup2]][df][sup2]~Phi[sub2] / [df]sigma[sup2] + 2[df]sigma / [df]x [df][sup2]~Phi[sub2] / [df]sigma[df]x + [df]~Phi[sub2] / [df]sigma[[df] / [df]sigma([df]sigma / [df]x)[sup2] - [df]sigma / [df]z [df] / [df]~x([df]z / [df]~x)] = [df]sigma / [df]zõmega[sup2] - [df]sigma / [df]y [df]sigma / [df]x [df][sup2]~Phi[sub1] / [df]sigma[sup2] - [df][sup2]~Phi[sub1] / [df]~y[df]~x - [df]sigma / [df]y [df][sup2]~Phi[sub1] / [df]sigma[df]~x - [df]sigma / [df]x [df][sup2]~Phi[sub1] / [df]sigma[df]~y - [df]~Phi[sub1] / [df]sigma[[df] / [df]sigma([df]sigma / [df]x [df]sigma / [df]y) - [df]sigma / [df]z [df] / [df]~x ([df]z / [df]~y)]. (16)
Following Pielke (1984), covariant and contravariant tensors are represented by subscripts and superscripts, respectively. Note that only the underlined quantities in Eqs. (15) and (16) contribute in hydrostatic flows.
Velocities are finally recovered from the streamfunctions, via (8) now applied to ~Phi:
rho[sub0]~u[sup1] = [df]sigma / [df]z [df]~Phi[sub2] / [df]{Begin Greek}s, r{End Greek}[sub0]~u[sup2] = [df]sigma / [df]z [df]~Phi[sub1] / [df]sigma, and rho[sub0]~u[sup3] = [df]sigma / [df]z([df]~Phi[sub1] / [df]~y + [df]~Phi[sub2] / [df]~x). (17)
The thermal energy and specific humidity equations (2) and (3), respectively, can be transformed in a similar manner:
[df]theta[subM] / [df]t = - 1 / rho[sub0] [df]sigma / [df]z [df] / [df]~x(rho[sub0]~u[sup1]theta[subM][df]z / [df] sigma) - 1 / rho[sub0] [df]sigma / [df]z [df] / [df]~y(rho[sub0]~u[sup2]theta[subM][df]z / [df]sigma) - 1 / rho[sub0] [df]sigma / [df]z [df] / [df]{Begin Greek}s(r{End Greek}[sub0]~u[sup3]theta[subM][df]z / [df]sigma) + [df]sigma / [df]z [df] / [df]~x([df]z / [df]sigma K[subH] [df]theta[subM] / [df]~x) + [df]sigma / [df]z [df] / [df]~y([df]z / [df]sigmaK[subH][df]theta[subM] / [df]~y) + [df] / [df]sigma(K[subh][df]theta[subM] / [df]sigma) - 1 / rho[sub0] [df]~u[sup3]rho[sub0]theta[sub0] / [df]sigma - ([df]sigma / [df]z)(1 / rho[sub0]c[subp] [df]Q / [df]sigma) and
[df]q[subM] / [df]t = - 1 / rho[sub0] [df]sigma / [df]z [df] / [df]~x(rho[sub0]~u[sup1]q[subM][df]z / [df]sigma) - 1 / rho[sub0] [df]sigma / [df]z [df] / [df]~y(rho[sub0]~u[sup2]q[subM][df]z / [df]sigma) - 1 / rho[sub0] [df]sigma / [df]z [df] / [df]{Begin Greek}s(r{End Greek}[sub0]~u[sup3]q[subM][df]z / [df]sigma) + [df]sigma / [df]z [df] / [df]~x([df]z / [df]sigmaK[subH][df]q[subM] / [df]~x) + [df]sigma / [df]z [df] / [df]~y([df]z / [df]sigmaK[subH] [df]q[subM] / [df]~y) + [df] / [df]sigma(K[subh][df]q[subM] / [df]sigma) - 1 / rho[sub0] [df]~u[sup3]rho[sub0]q[sub0] / [df]sigma.
TVM applies a 1.5-order closure scheme (Therry and Lacarerre 1983) to predict vertical diffusion coefficients (Schayes et al. 1996), which leads to a supplementary prognostic equation for turbulent kinetic energy (TKE).
1) ~Phi[sub1] = ~Phi[sub2] = constant (zero for convenience).
2) For thermally forced circulations, a no-slip boundary is chosen (friction). In the vicinity of the surface, vertical velocities disappear and vorticities simplify into vertical gradient of horizontal velocities, that is,
õmega[sup1] ~ [df]~u[sub2] / [df]sigma, õmega[sup2] ~ [df]~u[sub1] / [df]sigma, (18)
for which an analytical SBL expression is obtained:
[df]V[subH](z) / [df]z = u* / k[sub0]z[sub0]Phi[subm].
Vorticity boundary conditions (18) therefore become
õmega[sup1] = ~u[sub2] / V[subH] u* / k[sub0]z[sub0]Phi[subm], õmega[sup2] = ~u[sub1] / V[subH] u* / k[sub0]z[sub0]Phi[subm].
For dynamically forced circulations, a free-slip boundary condition (no friction) is necessary for comparison with analytical solutions. The resulting boundary conditions for vorticities transform into õmega[sup1] = õmega[sup2] = 0.
4) For thermally forced circulations, TKE is fixed at 4u[sup2[sub[sub*], following Therry and Lacarrère (1983), whereas TKE is not used for dynamically forced circulations.
1) Wind is assumed to be in geostrophic equilibrium:
[df]~Phi[sub1] / [df]sigma = [df]z / [df]{Begin Greek}sr{End Greek}[sub0]v[subg], [df]~Phi[sub2] / [df]sigma = [df]z / [df]{Begin Greek}sr{End Greek}[sub0]u[subg].
2) Vorticities are fixed to zero.
3) Temperature and specific humidity are fixed at largescale values (kept constant for the cases presented hereafter).
4) For vertically propagating perturbations, a damping Rayleigh friction layer similar to that of Tripoli and Cotton (1982) is used (for temperature only) in the highest several layers to avoid reflection.
1) Topography is assumed to be constant in directions perpendicular to boundaries.
2) Zero-gradient boundary conditions were selected for all variables.
Note that this set of boundary conditions is only representative of the idealized bidimensional applications considered in the present study.
The procedure by which TVM solves the PBL hydrodynamic and thermodynamic transport equations at each time step in the calculation is described in Schayes et al. (1996). For each prognostic equation, first each one-dimensional advection term is successively solved, then the three one-dimensional diffusion equations, and finally the remaining body force or source/sink terms. Finite-difference calculations are performed on the non-uniformly spaced, three-dimensional staggered Arakawa C grid, which locates TKE values at the center of horizontal cell faces.
Advection terms are approximated by the piecewise parabolic method (Collela and Woodward 1984; Carpenter et al. 1990), which may be classified as a "finite volume" scheme that possesses third-order accuracy in both time and space. As it was originally developed for astrophysical flows with strong shocks and discontinuities, it is especially designed to reduce explicit numerical diffusion without creating numerical instabilities, that is, over- and undershooting. This characteristic that was valid for 1D applications has been extended to multidimensional applications by Clappier (1998). The diffusion, vorticity tilting, and buoyancy terms are approximated by classical forward in time, centered in space differencing, which possesses second-order accuracy in space and first in time.
The elliptic streamfunction equations (12) are solved through an iterative smoothly converging variant of the biconjugate gradient method (van der Vorst 1992), using the hydrostatic solution as the initial guess at each time step.
While the change from an H to an NH formulation in primitive equation mesoscale models requires one additional prognostic equation for w, the current vorticity formulation contains the same equations [there are only a few additional terms in the elliptic equation (12)], solved by identical numerical methods. In addition, the current formulation only requires solution of two bidimensional elliptic equations for streamfunctions, as opposed to a similar three-dimensional equation for pressure in primitive equation formulations. The required storage is reduced approximately by a factor N/2, where N is the number of vertical grid levels. Additional advantages, as well as disadvantages, of the vorticity approach are provided in Schayes et al. (1996).
The objective of this section is to evaluate the ability of the TVM formulation to reproduce different types of atmospheric circulations that differ by their degree of nonhydrostaticity, nonlinearity, and anelasticity. Mountain wave circulations constitute a good evaluation test for mesoscale models since analytical solutions are available. A brief summary of the linear mountain wave theory is first given. Comparison to linear hydrostatic mountain waves is then used to evaluate the anelastic and incompressible formulations. The behavior of the incompressible version of TVM to nonlinear types of forcings is then investigated by use of the nonlinear analytical solutions proposed by Laprise and Peltier (1989). Finally, the Boulder windstorm is used to test the anelastic formulation to a highly nonlinear forcing.
Fundamental properties of mountain flows have been mainly studied with an assumption of a two-dimensional, nonviscous, adiabatic flow. Solutions for such linear flows are reviewed by Smith (1979) and Durran (1990). Combination of the five basic governing equations (u and w momentum, continuity, the ideal gas law, and the first law of thermodynamics) yields the following deep Boussinesq equation for perturbation vertical velocity w[subM]:
[df][sup2]([square root]rho[sub0](z) / rho[sub0](0) w[subM]) / [df]x[sup2] + [df][sup2]([square root]rho[sub0](z) / rho[sub0](0) w[subM]) / [df]z[sup2] + l[sup2](z)([square root]rho[sub0](z) / rho[sub0](0) w[subM]) = 0. (19)
The Scorer parameter l, defined by
l[sup2](z) = N[sup2] / U[sup2[sub[sub0] - 1 / U[sub0] d[sup2] U[sub0] / dz[sup2] - 1 / U[sub0] d lnrho / dz dU[sub0] / dz - 1 / 4(d lnrho / dz)[sup2] - 1 / 2 d[sup2] lnrho / dz[sup2],
represents an inverse intrinsic vertical scale and determines whether wave modes can propagate vertically. Note that only the first two terms remain when the incompressible approximation is invoked.
Standard practice is to solve (19) for flows with uniform wind speed [U[sub0](z) = U[sub0]] and stability [N(z) = N] profiles over a bell-shaped mountain described byz[subg](x) = h[subm]a[sup2] / x[sup2] + a[sup2]. (20)
Hydrostatic mountain wave numerical simulations constitute a standard evaluation case of mesoscale models (Durran and Klemp 1983; Yang 1993). For such cases, F[subx] << 1, the solution of (19) simplifies to
w[subM](x, z) = [square root]rho[sub0](O) / rho[sub0](z)U[sub0]h[subm]a (x[sup2] - a[sup2]) sinlz - 2xa coslz / (a[sup2] + x[sup2])[sup2] (21)
Relation (21) is illustrated in Fig. 1a for a mountain and atmosphere having the following characteristics: a = 10 km, h[subm] = 50 m, U[sub0] = 10 m s[sup-1], and N = 0.01 s[sup-1]. These parameters yield F[subx] = 10 and F[subz] = 0.05, with the latter ratio indicating that topography height h[subm] is small enough to ensure mostly linear flow. For clarity, both analytical and numerical results are amplified by a factor of 10. The figures clearly show internal gravity waves propagating vertically and affecting the atmosphere through its entire depth. When density variations are accounted for (Fig. 1b), the anelastic analytical solution shows an intensification of the vertical velocities with height.
For the following simulations, identical (to the analytical) values of stability, wind speed, and topography aspect ratio were used to initialize TVM. The simulation domain (81 × 28 grid points) covers 121 km horizontally with a constant grid resolution of 1500 m, whereas the model top is fixed at 25 km with a resolution ranging from 100 m near the surface to 1000 m near the model top.
TVM results reach a quasi-stationary state after a dimensionless time of t[sup*] = tU[sub0]/a = 14 (required for adjustment between thermal and dynamic fields). Wave amplitude is well represented for both the incompressible (Fig. 1c) and anelastic (Fig. 1d) solutions as seen from the vertical velocity cells and the simulated vertical wavelength is close to the analytical value (i.e., 6.3 km).
In summary, TVM produces hydrostatic linear wave results that are in close agreement with the analytical, incompressible, and anelastic solutions.
Both H and NH analytical solutions to the incompressible mountain wave problem have been presented by Laprise and Peltier (1989) for nonlinear types of forcing. In their work, the correct nonlinear surface boundary condition is introduced iteratively through a numerical procedure. One of the parameters they looked at in their study is the maximum steepness of the streamlines. The steepness predominantly reflects the amplitude of the internal waves launched by flow over the obstacle and, as such, it scales almost linearly with obstacle height. Figure 2a (from Laprise and Peltier) displays the maximum steepness as a function of the obstacle height (F[subz]) and half-width (F[subx]), in nondimensional terms for different types of flows. As seen from this figure, some attenuation occurs in NH conditions, as lee waves start to develop and progress horizontally ahead of the vertically propagating H waves discussed earlier. An example of such NH type of circulation (F[subx] = 1, F[subz] = 0.5) is given in Fig. 3a for the vertical wind speed fields. The nonhydrostaticity of such flow is clear when comparing with the solution obtained with the H formulation (Fig. 3b). Comparisons of H versus NH solutions for various mountain wave cases are also discussed by Keller (1994), Elkhaffi and Carissimo (1993), Sharman et al. (1988), or Yang (1993).
In this section, different simulations corresponding to hydrostatic (F[subx] = 10) and nonhydrostatic (F[subx] = 2) circulation types will be made for different degrees of nonlinearity (F[subz] = 0.25, F[subz] = 0.5, and F[subz] = 0.75) and the maximum steepness of the streamlines for each of those six cases is compared to the analytical values. For the NH forcing (F[subx] = 2), the parameters a = 2000 m, N = 0.01 s[sup-1], and U = 10 m s[sup-1] have been selected for three mountain heights [h[subm] = 250 m (F[subz] = 0.25), 500 m (F[subz] = 0.5), and 750 m (F[subz] = 0.75)]. The horizontal and vertical grids are stretched and have as highest resolution 1 km and 200 m (in the first 7 km), respectively. The grid domain covers 85 km by 24 km with 61 by 72 grid points, respectively. For the H cases (F[subx] = 10), the parameters a = 5000 m, N = 0.01 s[sup-1], and U = 5 m s[sup-1] have been selected for three mountain heights [h[subm] 125 m (F[subz] = 0.25), 250 m (F[subz] = 0.5), and 375 m (F[subz] = 0.75)]. The horizontal domain covers 80 km with a uniform resolution of 1 km whereas the vertical grid has a maximum resolution of 100 m (90 grid points). All cases were simulated without any explicit horizontal and vertical diffusion, the only extra filter being the diffusive layer at the top of the model to damp vertically propagating waves. Values of the maximum steepness are obtained at a height corresponding to 3/4 of the vertical wavelength of the flow, that is, at 4700 and 2350 m for the NH and H forcings, respectively. The time evolution of the nondimensional steepness values for those six cases may be seen in Fig. 2b. As can be seen from this figure, convergence is achieved much faster for moderately linear forcings for both NH and H forcings. For the same level of nonlinearity, convergence is achieved in a smoother manner for NH simulations. The resulting convergence values are overplotted on the analytical curves of Laprise and Peltier (1989) in Fig. 2a.
As can be seen, TVM results agree well with the analytical values indicating the good behavior of the code for both H and NH nonlinear circulations.
In January 1972, Boulder, Colorado, on the eastern slope of the Rocky Mountains, experienced a severe chinook windstorm. Although similar winds have been observed in the lee of other mountain ridges, this storm is unique for the richness of its experimental in situ dataset (Lilly and Zipser 1972; Lilly 1978). During the storm, a train of lee waves propagated in the lower part of the atmosphere east of the Rockies, and even more spectacularly, a strong stratospheric air descent developed in the lee of the Continental Divide. The corresponding wind field (Fig. 4a) shows a strong maxima of over 60 m s[sup-1] on the downslope edge, and near-zero wind speed regions above and farther downwind.
Initial conditions in previous simulations (e.g., Durran and Klemp 1983) were selected for the TVM simulation, for example, a simplified bell-shape mountain with an elevation change of 2000 m and a half-width of 10 000 m. Initial temperature and wind soundings (Figs. 4b, c) were obtained at Grand Junction, 300 km west (upstream) of Boulder. The temperature sounding indicates two stable layers: the first from 0 to 10 km, which includes a shallow weak inversion layer near 3000 m and the isothermal lower stratosphere above 10 km. The simulation domain (54 × 25 grid points) has a constant grid resolution of 1500 m horizontally, whereas the vertical resolution ranges from 100 m near the surface to 1000 m near the model top (fixed at 21 km).
At the time of fully developed windstorm, the NH-anelastic version of TVM produces horizontal wind speeds on the lee slope that reach 70 m s[sup-1] (Fig. 5a), the stratospheric air descent goes down to 6 km, and the train of waves propagating in the lee of the mountain is well developed. These results exhibit similar qualitative and quantitative patterns to the observations as well as to the results obtained by Durran and Klemp (1983) and Buty et al. (1988), both of whom used nonhydrostatic models, although slight small-scale differences exist.
This test case was mainly performed to test the behavior of TVM in highly nonlinear conditions but also to check its anelastic formulation. Indeed, anelasticity capability is required to correctly simulate the Boulder case due to the vertical extent of its circulation. To illustrate this fact, the incompressible version has been run on the same case (Fig. 5b). The results show a pattern that is qualitatively similar to the anelastic one but underestimated. In this configuration, the model only produces a maximum horizontal velocity of 55 m s[sup-1], the negative speed area in altitude is not present, and the stratospheric air descent is limited here to 9 km. It should be pointed out that the differences among different models (all being anelastic) for a similar test case (although initial conditions were slightly different) can also be quite significant, as illustrated by Doyle et al. (2000) who presented an intercomparison among 11 different anelastic models for the 1972 Boulder windstorm. In their work, differences are shown to be relatively large and to arise from factors such as subgrid-scale diffusion or lateral boundary conditions. Note also that the different levels of refinement that can be used in the same model to approximate anelasticity may also lead to significant differences, as shown by Nance and Durran (1994).
In conclusion, TVM has been shown capable of accurately simulating one of the most nonlinear observed atmospheric situations. Results were similar to observed values and to those from comparable primitive equation models.
As introduced in the previous section, one of the characteristics of the vorticity formulation is to keep a quasi-identical formulation for both the H and NH sets of equations. Indeed, the two prognostic vorticity equations remain unchanged. The only differences appear in the elliptic equations that link streamfunctions and vorticities where extra terms are present for the NH formulation. In such a configuration, the numerical schemes used in the two formulations (H and NH) are identical. This feature of the vorticity formulation ensures us that only the physical aspects of the approximation are isolated without any interference that could result from numerical treatment. In primitive equation models, this is not so straightforward since the set of equations differs for the H and NH formulations (one supplementary prognostic equation is required for w).
Throughout this work, all mentioned relative errors between NH and H solutions are obtained through the following relation (Martin and Pielke 1983):
P(%) = 200 (|A[subNH]| - |A[subH]|) / (|A[subNH]| + |A[subH]|),
where |A[subNH]| and |A[subH]| represent, respectively, the maximum NH and H absolute amplitudes (over space) for this variable.
Previous studies have investigated the importance of NH effects in thermally induced circulations (e.g., Martin and Pielke 1983; Yang 1991). Yang (1991) considered the special case of two bidimensional sea breezes converging over an island and tested the validity of the H assumption for several situations by varying critical parameters, that is, the geostrophic wind, the thermal stability, or the Coriolis forcing. This section is devoted to repeating one of the previous study test cases and comparing the vorticity H-NH differences for this specific case. Note that the objective here is not to present an exhaustive study that would investigate the NH impact of a series of significant parameters. The first selected test case is one of Yang's (1991) in which he investigated the importance of the horizontal heating scale (island extension) on the circulation. A second series of tests is then carried out to check the validity of the hydrostatic assumption in the presence of complex terrain. In this second configuration, sea breezes interact with slope winds and the effect of mountain height and slope on the validity of the H assumption is investigated. As the validity of the incompressible approximation is often not known for thermally induced mesoscale circulations, differences between incompressible and anelastic solutions are also shown hereafter (for the two geographical setups), although the vorticity formulation is comparable to a primitive equation approach to establish the validity of this approximation. Note that all results presented below are valid only for dry conditions. An investigation of the impact of the hydrostatic assumption on moist convection can be found in Kato (1997).
This case is designed as follows: a 16-km-wide island is bordered on each side by a 20-km-wide water surface (this case is denoted SB). For island surface grid points, a sinusoidal wave temperature variation is calculated via theta(t) = alpha sin(2pit/tau[subd]) + theta[sub0], (22)
The time evolution of the generated circulations is the following: at 0300 LST, shallow land breezes are initiated on the two sides of the island with maximum wind speeds of about 1 m s[sup-1]. By 0900 LST those land breezes have turned into sea breezes of the same intensity, which converge around 1200 LST. Sea breezes then reach their maximum intensity and progressively damp afterward as the surface heat supply decreases. The sea breezes at maximum intensity (13:00 LST) are shown in Fig. 6a. Note that left and right breeze cells are symmetric with respect to the middle of the island as the large-scale wind is zero.
A time comparison of (NH-anelastic), (H-anelastic), and (NH-Incompressible) maximum horizontal (Fig. 6c) and vertical (Fig. 6e) velocities indicates that NH circulations are weaker than H ones. This observation agrees with results obtained by both Martin and Pielke (1983) and Yang (1991). For this spatial resolution (1000 m), however, no major difference is observed in the horizontal or vertical wind components. Indeed, the first are virtually identical whereas vertical speeds exhibit differences that do not exceed 10% during the maximum convergence period (from 1200 to 1400 LST).
If the dimensions of the island are now reduced by a factor of 2 (Fig. 6b), that is, the island and water surfaces now have widths of 8 and 10 km, respectively, and if the grid resolution is refined to 500 m [this corresponds exactly to a second case presented in Yang (1991)], the comparison of H and NH results leads to the following comments. 1) The NH breezes are still weaker than H ones, but despite the relatively fine resolution, horizontal wind fields look quantitatively and qualitatively similar (P < 10%). 2) Significant differences are now observed for vertical speeds for which peak values of P ˜ 20% are reached between 1100 and 1400 LST when breezes are at their maximum intensities.
Qualitatively similar results were presented in Yang (1991) although differences between H and NH fields were much larger. Indeed, he obtained at maximum seabreeze development the following values: for the 1-km resolution case, P values of 6% for horizontal velocity, and 28% and 10% for downward and upward velocities, respectively, and the corresponding TVM values are 1%, 5%, and 7%; and for the 500-m resolution case, P values of 38% for horizontal velocity, and 76% and 28% for downward and upward velocities, respectively, and the corresponding TVM differences are 2%, 5%, and 14%.
The large discrepancies between H and NH simulations shown by Yang probably arise partially from different numerical treatment in the H and NH formulations.
In Figs. 6c-f, a comparison with the incompressible version of TVM is also provided. Almost no difference is observed when density variations are neglected as could be expected since the vertical extent of the circulations is limited.
Although the two test cases presented above, corresponding to alpha = 8 in (22), were among the most NH cases proposed by Yang (they indeed range from alpha = 4 to alpha = 10), H-NH differences obtained with the vorticity formulation remain small. To test the behavior of the H approximation in more stringent cases, the same tests (two island cases) with alpha = 16 are realized, and all other parameters remaining unchanged. This doubling of the thermal forcing still leads to differences below 10% for the horizontal wind component for both the 500- and 1000-m horizontal resolutions, whereas differences in the vertical wind field now reach 21% and 52% for the 1000- and 500-m horizontal resolutions, respectively. Differences between the incompressible and anelastic formulations do not exceed 10% for all fields and resolutions. Table 1 summarizes the quantitative differences for the different tests presented above (denoted as SB).
In conclusion, although differences in vertical velocity fields are well marked, they do not significantly influence other variable fields as the grid resolution remains large enough. The hydrostatic assumption seems to remain valid even for relatively small scales when thermal forcing is not too strong. Similar conclusions were given by Martin and Pielke (1983) who, by use of NH and H models, stated that the hydrostatic assumption was valid for sea breeze when the heating scale is as small as 6 km over flat terrain. For strong forcings, however, significant errors in the calculation of the vertical wind field are present and could be of importance in air pollution modeling.
Another possible way to generate NH effects is the presence of complex terrain. In this section, we construct an idealized case that consists of a mountain (h[subm] = 500 m, a = 10 km) boarded on its left by a 25-km-wide water surface (this case is denoted MSB). A sinusoidal wave temperature identical to (22) with alpha = 8 is prescribed for land grid points while water temperature remains constant. The grid configuration (75 × 24 grid points) has a resolution of 1 km horizontally and is stretched along the vertical direction. The initial lapse rate is 4 K km[sup-1].
At 0400 LST, slope winds develop on each side of the mountain and converge on its western side with the land breeze. At 0900 LST, both a sea-breeze cell on the coastline and two anabatic cells on the mountain slopes have formed. Their respective intensities are 1.5 and 0.9 m s[sup-1]. At 1500 LST, the sea-breeze cell, which has then reached an intensity of 3.5 m s[sup-1], merges with the left upslope wind cell. At 1800 LST, the sea-breeze core has passed the ridge (Fig. 7a) and is then subject to a strong downslope acceleration to finally reach a maximum velocity of 4.8 m s[sup-1].
Differences between NH and H simulations in the along-slope velocity (Fig. 7c) component are surprisingly weak. The H winds generally overestimate NH ones but errors remain under 5%. For vertical wind speeds (Fig. 7e), the error also remains below 10% except when the sea breeze overcomes the ridge and is subject to a strong downslope acceleration (it then reaches 15%). If the mountain has a height of 1000 m (Figs. 7b,d,f), differences remain very weak and, for vertical wind speed, they are even weaker since the sea breeze then never overcomes the ridge.
As in the previous section for sea breezes over flat terrain, the two tests with different mountain heights are repeated for alpha = 16, that is, when the thermal surface forcing is doubled. In those two cases, differences remain below 10% for the horizontal wind whereas they reach 45% and 25% for the vertical wind speed for mountain heights of 500 and 1000 m, respectively. Note that in this case the two sea breezes have enough energy to pass over the ridge but do not prevent the maximum differences occurring for the lowest hill.
Differences between incompressible and anelastic simulations are below 10% for the horizontal components whereas they reach 15% and 30% for the vertical wind speed for the 500- and 1000-m heights, respectively
Note that all results presented in sections 3 and 4 were obtained by neglecting Coriolis forcing. No radiative lateral boundary conditions were used but a damping layer was present in upper vertical layers for the mountain waves and Boulder windstorm cases. The horizontal and vertical diffusion coefficients were fixed at zero for all dynamically forced circulations. For thermally induced circulation, the horizontal diffusion coefficient was fixed at a value proportional to the grid size (-200 m s[sup-2]) while vertically, they were obtained through the TKE closure.
This paper described the formulation and application of the NH anelastic vorticity model (TVM). Evaluation of the new formulation was realized by simulations of 2D H and NH mountain wave cases with different degrees of nonlinearity. Various types of topography-induced waves, including the well-known 1972 Boulder windstorm, were numerically simulated, and results compared to available analytical solutions, measurements, and/or results from other primitive equation mesoscale models.
The change from an H to an NH formulation in primitive equation mesoscale models requires one additional prognostic equation for the vertical velocity, whereas the current vorticity formulation contains the same prognostic and diagnostic equations (additional terms are present in the streamfunction NH elliptic equation only), solved by identical numerical schemes. The model was constructed for deep thermal convection flows using a method involving only two horizontal streamfunctions and two horizontal vorticity components, shown to be sufficient to totally describe three-dimensional wind fields. In addition, the current anelastic formulation only requires solution of two bidimensional elliptic equations for streamfunctions, as opposed to a similar three-dimensional equation for pressure in primitive equation formulations. The decrease in required storage in TVM over primitive equation models for the same 3D grid is thus approximately 2/N, where N is the number of vertical grid levels. A generalized vertical coordinate system was introduced to account for complex terrain.
TVM in both its NH and H formulations was shown to accurately (qualitatively and quantitatively) reproduce H and NH mountain wave behavior for both linear and nonlinear forcings. As noted, results were obtained without any explicit diffusion.
The adequacy of the hydrostatic assumption in simulating thermally induced circulations has also been investigated. For a moderate surface thermal forcing, typical geographical setups were used and showed very slight differences between H and NH horizontal wind speeds. Indeed, a decrease in grid resolution (down to 500 m) or an increase in mountain slopes (up to 15%) did not generate differences over 10%. For vertical wind speeds, differences were shown to be much larger and more sensitive to changes in grid resolution. In fact, an increase in mountain slope even reduced the differences between H and NH simulations since in the steeper case, the slope prohibited the passage of the core of the sea breeze over the ridge and prevented the successive downslope acceleration in which significant NH effects are present. When a more stringent thermal surface forcing was imposed, differences remained low for horizontal wind speeds but considerably increased for the vertical wind component (above 20% for all cases) illustrating the failure of the hydrostatic assumption in reproducing those circulations.
The comparison between anelastic and incompressible solutions for the same cases showed the adequacy of the incompressible assumption in simulating thermally driven circulations that are forced by the surface and that are characterized by a relatively shallow vertical extent. In such conditions, virtually no differences are observed between the two formulations. The only exception being the 1-km-high mountain with a strong thermal surface forcing for which differences reached 30% in the calculation of the vertical wind speed.
The comparison of this work with previous numerical studies related to the impact of assumptions on atmospheric conservation equations incites prudence in interpreting the results. Indeed, the use of different numerical treatments for different formulations could lead to an overestimation of the impact of a given assumption. As long as an identical numerical treatment is not used for the comparison of different formulations, the impact of a given approximation will remain doubtful.
Although the different tests presented here were all bidimensional in order to better evaluate the importance of the hydrostatic assumption, validation studies of this vorticity formulation on three-dimensional coastal, urban, and/or complex terrain conditions have also been performed. Applications over the Athens basin (Greece) and over the Valencia area (Spain) incorporating photochemistry are presented in Clappier et al. (2000), Grossi et al. (2000), and Thunis and Cuvelier (1999).
ADDED MATERIAL
P. THUNIS
Environment Institute, Joint Research Centre, Ispra, Italy
A. CLAPPIER
Laboratoire d'Etude de la Pollution de l'Air, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
Corresponding author address: Dr. Philippe Thunis, Environment Institute, Joint Research Centre, TP 280, 21020 Ispra, Italy.
E-mail: philippe.thunis@jrc.it
Acknowledgments. The help of Prof. G. Schayes, Dr. G. Graziani, Prof. R. Bornstein, and Dr. B. Piraux is deeply acknowledged.
TABLE 1. Quantitative summary table of the H - NH and anelastic -- incompressible (A - I) differences (factor P) for the eight thermally induced circulation cases presented in section 4. Values of P are representative of the fully developed sea-breeze period. Case U[subH-NH] W[subH-NH] U[subA-I] W[subA-I]
SB1000 (alpha = 8) <10% <10% <10% <10%
SB1000 (alpha = 16) ˜10% ˜20% ˜10% <10%
SB500 (alpha = 8) <10% ˜20% <10% <10%
SB500 (alpha = 16) <10% ˜50% <10% <10%
MSB500 (alpha = 8) <10% ˜15% <10% <10%
MSB500 (alpha = 16) <10% ˜45% <10% ˜15%
MSB1000 (alpha = 8) <10% <10% <10% ˜15%
MSB1000 (alpha = 16) <10% ˜25% ˜10% ˜30%
FIG. 2. (a) Time evolution of the maximum nondimensional steepness (eta[subz]) as a function of the obstacle width (F[subx]) and of the obstacle height (F[subz]). (b) Analytical maximum steepness (from Laprise and Peltier 1989) in terms of F[subx] and F[subz]. TVM numerical values for the six cases in (a) are indicated by solid circles.
FIG. 3. (a) TVM nonhydrostatic and (b) hydrostatic incompressible solutions for vertical wind velocity (m s[sup-1]) for conditions characterized by F[subz] = 0.5 and F[subx] = 1, corresponding to a nonhydrostatic, moderately nonlinear type of surface forcing.
FIG. 4. (a) West-east aircraft cross section through 1972 Boulder windstorm for horizontal speed (m s[sup-1]), and Grand Junction rawinsonde soundings for 1972 Boulder windstorm for (b) temperature (K) and (c) for horizontal speed (m s[sup-1]) (from Peltier and Clark 1979).
FIG. 5. (a) TVM nonhydrostatic anelastic and (b) nonhydrostatic incompressible Boulder windstorm solutions for horizontal speed (m s[sup-1]).
FIG. 6. (a), (b) Horizontal wind speed obtained with TVM nonhydrostatic at full sea-breeze development for the (a), (c), (e) 1-km and (b), (d), (f) 500-m grid spatial resolutions. The shaded part of the horizontal axis is indicative of the land area. (c), (d) The time evolution of the maximum horizontal wind speeds and (e), (f) the time evolution of the maximum vertical wind speeds. The nonhydrostatic anelastic, hydrostatic anelastic, and nonhydrostatic incompressible results are represented by the continuous line, dashed line, and squares, respectively.
FIG. 7. Same as Fig. 6 except for the case MSB. (a), (b) The wind field for the (a) 500- and (b) 1000-m mountains, at 1800 LST.
Bornstein, R. D., S. Klotz, U. Pechinger, R. Salvador, R. Street, L. J. Shieh, F. Ludwig, and R. Miller, 1986: Application of linked three-dimensional PBL and dispersion models to New York City. Air Pollution Modeling and Its Application, Vol. V, Plenum Press, 543-564.
Bornstein, P. Thunis, P. Grossi, and G. Schayes, 1996: Topographic-vorticity mode mesoscale-beta (TVM) model. Part II: Evaluation. J. Appl. Meteor., 35, 1824-1834.
Buty, D., J. Caneill, and B. Carissimo, 1988: Simulation numérique de la couche limite atmosphérique en terrain complexe au moyen d'un modèle mésométéorologique nonhydrostatique: Le code MERCURE. J. Theor. Appl. Mech., 7, 35-52.
Carpenter, R. L., K. Droegemeier, P. R. Woodward, and C. E. Hane, 1990: Application of the Piecewise Parabolic Method (PPM) to meteorological modeling. Mon. Wea. Rev., 118, 586-612.
Clappier, A., 1998: A correction method for use in multidimensional time-splitting advection algorithms: Application to two- and three-dimensional transport. Mon. Wea. Rev., 126, 232-242.
Clappier, A., and Coauthors, 2000: Effect of sea breeze on air pollution in the greater Athens area. Part I: Numerical simulations and field observations. J. Appl. Meteor., 39, 546-562.
Collela, P., and P. Woodward, 1984: The Piecewise Parabolic Method (PPM) for gas dynamical simulations. J. Comput. Phys., 54, 174-201.
De Ridder, K., and G. Schayes, 1997: The IAGL land surface model. J. Appl. Meteor., 36, 167-182.
Doyle, J. D., and Coauthors, 2000: An intercomparison of model-predicted wave breaking for the 11 January 1972 Boulder windstorm. Mon. Wea. Rev., 128, 901-914.
Durran, D. R., 1990: Mountain waves and downslope winds. Atmospheric Processes over Complex Terrain, Meteor. Monogr., No. 45, Amer. Meteor. Soc., 59-81.
Durran, D. R., and J. B. Klemp, 1983: A compressible model for the simulation of moist mountain waves. Mon. Wea. Rev., 111, 2341-2361.
Elkhafi, A., and B. Carissimo, 1993: Numerical simulations of a mountain wave observed during the "Pyrenees Experiment": Hydrostatic/non-hydrostatic comparison and time evolution. Beitr. Phys. Atmos., 66, 183-200.
Gal-Chen, T., and R. Somerville, 1975: On the use of a coordinate transformation for the solution of the Navier-Stokes equations. J. Comput. Phys., 17, 209-228.
Grossi, P., P. Thunis, A. Martilli, and A. Clappier, 2000: Effect of sea breeze on air pollution in the greater Athens area. Part II: Analysis of different emission scenarios. J. Appl. Meteor., 39, 563-575.
Kato, T., 1997: Hydrostatic and nonhydrostatic simulations of moist convection: Review and further study. Meteor. Atmos. Phys., 63, 39-51.
Keller, T. L., 1994: Implications of the hydrostatic assumption on atmospheric gravity waves. J. Atmos. Sci., 51, 1915-1929.
Laprise, R., and W. R. Peltier, 1989: On the structural characteristics of steady finite-amplitude mountain waves over bell-shaped topography. J. Atmos. Sci., 46, 586-595.
Lilly, D. K., 1978: A severe downslope windstorm and aircraft turbulence event induced by a mountain wave. J. Atmos. Sci., 35, 59-77.
Lilly, D. K., and E. J. Zipser, 1972: The front range windstorm of 11 January 1972--A meteorological narrative. Weatherwise, 25, 56-63.
Malkus, J. S., and G. Witt, 1959: The evolution of a convective element: A numerical calculation. The Atmosphere and the Sea in Motion: Scientific Contributions to the Rossby Memorial Volume, B. Bolin, Ed., Oxford University Press, 425-439.
Martin, L. C., and R. A. Pielke, 1983: The adequacy of the hydrostatic assumption in sea breeze modeling over flat terrain. J. Atmos. Sci., 40, 1472-1481.
Moeng, C., and A. Arakawa, 1980: A numerical study of a marine subtropical stratus cloud layer and its stability. J. Atmos. Sci., 37, 2661-2676.
Nance, L. B., and D. R. Durran, 1994: A comparison of the accuracy of three anelastic systems and the pseudo-incompressible system. J. Atmos. Sci., 51, 3549-3565.
Orville, H. D., and L. J. Sloan, 1970: A numerical simulation of the life history of a rainstorm. J. Atmos. Sci., 27, 1148-1159.
Pearson, R. A., 1973: Properties of the sea breeze front as shown by a numerical model. J. Atmos. Sci., 30, 1050-1060.
Peltier, W. R., and T. L. Clark, 1979: The evolution and stability of finite-amplitude mountain waves. Part II: Surface wave drag and severe downslope windstorms. J. Atmos. Sci., 36, 1498-1529.
Pielke, R., 1984: Mesoscale Meteorological Modeling. Academic Press, 468 pp.
Saitoh, T. S., T. Shimada, and H. Hoshi, 1996: Modeling and simulation of the Tokyo urban heat island. Atmos. Environ., 30, 3431-3442.
Schayes, G., P. Thunis, and R. Bornstein, 1996: Topographic vorticity-mode mesoscale-beta (TVM) model. Part I: Formulation. J. Appl. Meteor., 35, 1815-1823.
Schumann, U., T. Hauf, H. Holler, H. Schmidt, and H. Volkert, 1987: A mesoscale model for the simulation of turbulence, clouds and flow over mountains: Formulation and validation examples. Beitr. Phys. Atmos., 60, 413-446.
Sharman, R. D., T. L. Keller, and M. G. Wurtele, 1988: Incompressible and anelastic flow simulations on numerically generated grids. Mon. Wea. Rev., 116, 1124-1136.
Sievers, U., and W. G. Zdunkowski, 1986: A microscale urban climate model. Beitr. Phys. Atmos., 59, 13-40.
Smith, R. B., 1979: The influence of mountains on the atmosphere. Advances in Geophysics, Vol. 21, Academic Press, 87-230.
Therry, G., and P. Lacarrère, 1983: Improving the eddy kinetic energy model for planetary boundary layer description. Bound.-Layer Meteor., 25, 63-88.
Thunis, P., and R. Bornstein, 1996: Hierarchy of mesoscale flow assumptions and equations. J. Atmos. Sci., 53, 380-397.
Thunis, P., and C. Cuvelier, 1999: Impact of biogenic emissions on ozone formation in the Mediterranean area--A BEMA modelling study. Atmos. Environ., 34, 467-481.
Thyer, N. H., 1966: A theoretical explanation of mountain and valley winds by a numerical model. Arch. Meteor. Geophys. Bioklimatol., 26, 39-50.
Tripoli, G. J., and W. R. Cotton, 1982: Colorado State University three-dimensional cloud/mesoscale model, 1982: Part I, General theoretical framework and sensitivity experiments. J. Rech. Atmos., 16, 185-219.
van der Vorst, H. A., 1992: BI-CGSTAB: A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems. J. Sci. Stat. Comput., 13, 631-644.
Yang, X., 1991: A study of nonhydrostatic effects in idealized sea breeze systems. Bound.-Layer Meteor., 54, 183-208.
Yang, X., 1993: A nonhydrostatic model for simulation of airflow over mesoscale bell-shaped ridges. Bound.-Layer Meteor., 65, 401-424.
a = Mountain half-width
c[subp] = Specific heat at constant pressure (= 1005 J kg[sup-1] K[sup-1])
f, f* = Coriolis parameters
F[subx] = Nondimensional mountain width
F[subz] = Nondimensional mountain height
g = Vector acceleration of gravity (magnitude = 9.81 m s[sup-2])
h[subm] = Maximum mountain height
K[subH] = Horizontal turbulent diffusion coefficient
K[subh], K[subm], = Heat, momentum, and specific humidity turbulent
K[subq] = diffusivity
k = Vertical unit vector
k[sub0] = Von Kármán constant (=0.40)
l = Scorer parameter
L[subH] = Mountain horizontal length scale
N = Brunt-Väisälä frequency
p = Pressure
q = Specific humidity
Q = Radiative flux
t = Time
t* = Dimensionless time
T = Temperature
u* = Friction velocity
u, v, w = East, north, and upwind component
V = Three-dimensional wind vectors
x, y, z = Eastward, northward, and upward coordinates
z[subg] = Topography height
z[sub0] = Roughness length
z[subt] = Model-top height
FTCS = Forward in time, centered in space
PBL = Planetary boundary layer
PPM = Piecewise parabolic method
SBL = Surface boundary layer
TKE = Turbulent kinetic energy
TVM = Thermal Vorticity Model
eta[subz] = Maximum streamline steepness
theta = Potential temperature
{Begin Greek}x, q {End Greek}= Vorticity vectors
rho = Density
sigma = Topography transformed vertical coordinate
phi = Latitude
Phi[subm] = SBL stability function for momentum
{Begin Greek}F, F{End Greek}* = Vector streamfunctions
Omega = Earth rotation vector (magnitude = 7.292 × 10[sup-5] rad s[sup-1])
g = Geostrophic
M = Mesoscale perturbation
m = Momentum
0 = Static reference state
() = Reynolds decomposed mean value (no oversymbol)
()' = Turbulent perturbation
(˜) = Topography coordinate system variable
[universal quantifier] = 3D del operator
[universal quantifier][subH] = Horizontal del operator
()[subi] = Covariant tensor component (i = 1-3, for x, y, z)
()[supi] = Contravariant tensor component (i = 1-3, for x, y, z)
()[sub;i] = Covariant derivative in the ith direction
Formulation and evaluation of a nonhydrostatic mesoscale vorticity model (TVM)
In: Monthly weather review, Jg. 128 (2000), Heft 9, S. 3236-3251
Online
academicJournal
- print, 41 ref
Zugriff:
This paper describes the formulation and the application of the nonhydrostatic anelastic vorticity model (TVM). This model is constructed using a method involving two horizontal streamfunctions and two horizontal vorticity components. The evaluation of this formulation is performed by simulating various bidimensional hydrostatic and nonhydrostatic mountain wave cases. Results are compared with analytical solutions and in particular with those developed by Laprise and Peltier for nonlinear forcings. The anelastic formulation is also validated with respect to the highly nonlinear 1972 Boulder windstorm. TVM is shown to accurately reproduce these mountain wave test cases in both its incompressible and anelastic formulations. In the second part of this paper, the adequacy of the hydrostatic and anelastic assumptions in simulating thermally induced circulations is investigated and compared to previous works. For a moderate surface thermal forcing, typical geographical setups are used and show slight differences between hydrostatic and nonhydrostatic horizontal wind speeds. For vertical wind speeds, differences are shown to be much larger and more sensitive to changes in grid resolution. For more stringent thermal surface forcing, differences remain low for horizontal wind speeds but increase considerably for the vertical wind component. The comparison between anelastic and incompressible solutions for the same cases shows the adequacy of the incompressible assumption when circulations are forced by the surface and are characterized by a relatively shallow vertical extent. In such conditions, virtually no differences are observed between the two formulations.
| Titel: |
Formulation and evaluation of a nonhydrostatic mesoscale vorticity model (TVM)
|
|---|---|
| Autor/in / Beteiligte Person: | THUNIS, P ; CLAPPIER, A |
| Link: | |
| Zeitschrift: | Monthly weather review, Jg. 128 (2000), Heft 9, S. 3236-3251 |
| Veröffentlichung: | Boston, MA: American Meteorological Society, 2000 |
| Medientyp: | academicJournal |
| Umfang: | print, 41 ref |
| ISSN: | 0027-0644 (print) |
| Schlagwort: |
|
| Sonstiges: |
|