Abstract

To better understand the processes involved in tropical cyclone development, the authors simulate an axisymmetric tropical-cyclone-like vortex using a two-dimensional model based on nonhydrostatic dynamics, equilibrium thermodynamics, and bulk microphysics. The potential vorticity principle for this nonhydrostatic, moist, precipitating atmosphere is derived. The appropriate generalization of the dry potential vorticity is found to be P = ρ−1 {(−∂υ/∂z) (∂θρ/∂r) + [ f + ∂()/rr] (∂θρ/∂z)}, where ρ is the total density, υ is the azimuthal component of velocity, and θρ is the virtual potential temperature. It is shown that P carries all the essential dynamical information about the balanced wind and mass fields. In the fully developed, quasi-steady-state cyclone, the P field and the θ̇ρ field become locked together, with each field having an outward sloping region of peak values on the inside edge of the eyewall cloud. In this remarkable structure, the P field consists of a narrow, leaning tower in which the value of P can reach several hundred potential vorticity (PV) units.

Sensitivity experiments reveal that the simulated cyclones are sensitive to the effects of ice, primarily through the reduced fall velocity of precipitation above the freezing level rather than through the latent heat of fusion, and to the effects of vertical entropy transport by precipitation.

1. Introduction

The potential vorticity conservation principle provides a basis for understanding midlatitude weather systems, both through balanced models and primitive equation models. In many tropical weather systems, such as the ITCZ and tropical cyclones, the release of latent heat plays a crucial role, so that potential vorticity is not materially conserved. However, even in these cases, the relevant dynamics is “slow manifold dynamics,” with adjusted wind and mass fields intimately related to the evolving potential vorticity field. Thus, the analysis of tropical flows in terms of potential vorticity dynamics yields insights into such processes as ITCZ breakdown, the formation of easterly waves, and into the extreme inner core structures of tropical cyclones. The primary purpose of the present paper is to study the inner core potential vorticity structure of tropical cyclones simulated with a high resolution, axisymmetric, nonhydrostatic tropical cyclone model whose thermodynamic/dynamic foundations and discretizations are based on the work of Ooyama (1990, 2001, 2002, hereafter O90, O1, O2, respectively). The model has several unique features: 1) a very accurate treatment of moist processes within the context of equilibrium thermodynamics, with a simple switch to include or exclude the effects of ice; 2) an associated potential vorticity (PV) principle and an invertibility principle, both exactly derivable from the original model equations; 3) the inclusion of the thermodynamic and dynamic effects of precipitation, in particular the vertical transport of entropy and momentum by precipitation; 4) spatial numerics based on the cubic spline transform (CST) method, which results in small computational dispersion errors and noise-free nesting.

We begin in section 2 by presenting a concise, self-contained description of the model. In section 3 and appendix B we derive the potential vorticity and invertibility principles associated with our cloudy, precipitating model atmosphere. The moist generalization of the dry Ertel potential vorticity turns out to be P = ρ−1ζ · θρ, where ρ is the total density of moist air, ζ is the absolute vorticity vector, and θρ is the virtual potential temperature. In the axisymmetric case, the vorticity vector has components in the vertical, radial, and azimuthal directions, but θρ has components in the vertical and radial directions only. As a consequence, the azimuthal component of ζ is lost in performing the product ζ · θρ, and the potential vorticity simplifies to P = ρ−1 {(−∂υ/∂z) (∂θρ/∂r) + [ f + ∂()/rr] (∂θρ/∂z)}, where υ is the azimuthal component of the wind. Since the evolving hurricane remains close to a state of gradient and hydrostatic balance, υ and ρ can be expressed in terms of the pressure p. Then, since θρ is a function of ρ and p only, the potential vorticity can be expressed solely in terms of the pressure. In other words, the P field contains all the required information about the balanced part of the wind and mass fields, and the invertibility principle, which determines the pressure from the potential vorticity, is an elliptic partial differential equation in the radial-vertical plane. Understanding the evolution of the P field is thus a crucial part of understanding the whole intensification problem.

In section 4 we present the results of the control experiment, paying particular attention to the potential vorticity dynamics. With 500-m horizontal and vertical resolution in the inner core, a remarkable, hollow tower potential vorticity structure emerges in the quasi–steady state, with values of potential vorticity as high as 275 PV units (PVU; where 1 PVU = 1.0 × 10−6 m2 s−1 K kg−1) in the eyewall. This is in sharp contrast to the “zero PV picture” that emerges when equivalent potential temperature or saturation equivalent potential temperature is used as the thermodynamic variable in the definition of potential vorticity.

Section 5 contains a discussion of selected sensitivity experiments, including the important effects of ice and the effects of precipitation terms in the entropy and momentum budgets. Concluding remarks are given in section 6.

2. Model

Consider atmospheric matter to consist of dry air, airborne moisture, and precipitation, with respective densities ρa, ρm, ρr, so that the total density is given by ρ = ρa + ρm + ρr.1 The mass density of airborne moisture is the sum of the densities of water vapor and airborne condensate, so that ρm = ρυ + ρc. However, the partition of ρm into ρυ and ρc is not considered in the prognostic stage, but only later in the diagnostic stage. In cylindrical coordinates, with the assumption of axisymmetry, the mass conservation law for dry air is ∂ρa/∂t + ∂(ρaru)/rr + ∂(ρaw)/∂z = 0, the advective form of which is (2.1). Similarly, the conservation law for ρr is ∂ρr/∂t + ∂(ρrru)/rr + ∂[ρr(w + W)]/∂z = Qr, where w + W is the vertical velocity of the precipitation (so that W is the fall velocity of the precipitation relative to the dry air and airborne moisture) and Qr is the rate of conversion from airborne moisture to precipitation. Defining the precipitation mixing ratio as μr = ρr/ρa, the conservation law for ρr can be combined with (2.1) and written in the advective form (2.3). Finally, the conservation law for the airborne moisture is ∂ρm/∂t + ∂(ρmru)/rr + ∂(ρmw + Fm)/∂z = −Qr, where Fm is the boundary layer turbulent flux of water vapor. Adding this conservation law for ρm to the conservation law for ρr and converting the result into an advective form for the total water mixing ratio μ = (ρm + ρr)/ρa, we obtain (2.2). Note that, although there are four types of matter (with densities ρa, ρυ, ρc, ρr), there are only three prognostic equations, (2.1)(2.3), for the distribution of mass. As we shall see, the separation of the airborne moisture density ρm into the vapor density ρυ and the airborne condensate density ρc will be accomplished diagnostically in the two alternatives of (2.14).

The total entropy density is σ = σa + σm + σr, consisting of the sum of the entropy densities of dry air, airborne moisture, and precipitation. Since the vertical entropy flux is given by σaw + σmw + σr(w + W) + ρaFs = σw + σrW + ρaFs, we can write the flux form of the entropy conservation principle as ∂σ/∂t + ∂(σru)/rr + ∂(σw + σrW + ρaFs)/∂z = 0, where Fs denotes the turbulent eddy flux in the atmospheric boundary layer, and where radiative effects have been neglected. Defining s = σ/ρa as the “dry-air-specific” entropy of moist air, we can combine (2.1) with the above entropy conservation principle to obtain (2.4).

We next consider the momentum equations. Defining the absolute angular momentum per unit mass as m = + ½fr2, we can write the absolute angular momentum budget as ∂(ρm)/∂t + ∂(ρmru)/rr + ∂[(ρw + ρrW + Fm)m]/∂z = −∂(ρarFυ)/∂z, where Fυ denotes the turbulent eddy flux of υ in the boundary layer. With the aid of the continuity equation for total density, this absolute angular momentum budget can be written as (2.6). In a similar fashion we can derive the radial wind Eq. (2.5), where p = pa + pυ is the sum of the partial pressures of dry air and water vapor. The derivation of the vertical equation of motion is somewhat more complicated. While all the matter moves with the same horizontal velocity, the vertical velocity of dry air and airborne moisture is w, while the vertical velocity of precipitation is w + W. The prediction of both w and w + W is equivalent to predicting W, which is inconsistent with the diagnostic treatment of W through parameterized cloud microphysics. A solution to this problem, proposed in O1, is to write a single budget equation for the total vertical momentum ρw + ρrW, and then approximate this equation by neglecting the material derivative of W along the precipitation path, that is, by neglecting ∂W/∂t + u(∂W/∂r) + (w + W)(∂W/∂z). With this approximation our vertical momentum equation becomes (2.7). The neglect of the material derivative of W along the precipitation path is consistent with the parameterization Eq. (2.18), which gives a slow variation of W because of the small fractional power of ρr and the inverse square root of ρa. The largest errors due to this assumption are expected in a small region near the melting level, where the ice factor fice results in an increase of the fall velocity of the precipitation.

Consolidating our results so far, the prognostic equations for the mass density of dry air ρa, the total water mixing ratio μ, the precipitation mixing ratio μr, the specific entropy s, and the velocity components u, υ, w are

 
formula
 
formula
 
formula
 
formula
 
formula
 
formula
 
formula

Associated with these seven prognostic equations is a set of diagnostic equations that is required to determine ρ, σr, p, W, Qr, Fm, Fs, Fu, Fυ. The diagnostic set of equations can be divided into three subsets: thermodynamic diagnosis to determine ρ, σr, p; precipitation microphysics diagnosis to determine W and Qr; and air–sea interaction diagnosis to determine Fm, Fs, Fu, Fυ. The equations for the thermodynamic diagnosis are

 
formula
 
formula
 
formula
 
formula
 
formula
 
formula
 
formula
 
formula

where the known functions S1, S2, E, C, ρ*υ are defined in appendix A. The thermodynamic diagnosis can be interpreted conceptually as follows: Input {ρa, μ, μr, s} ⇒ Output{ρ, ρr, ρm, ρυ, ρc, T1, T2, T, σr, pa, pυ, p}. Specifically, starting with the values of the prognostic, thermodynamic state variables ρa, μ, μr, s, the thermodynamic diagnosis proceeds in the order given, that is, determination of the densities ρ, ρr, ρm from (2.8), the thermodynamically possible temperatures T1, T2, and the entropy density of precipitation σr from (2.9) to (2.11), the actual temperature T of the gaseous matter from (2.12), the partial pressure of dry air from (2.13), the water vapor density ρυ, cloud condensate density ρc, and water vapor partial pressure pυ from the appropriate alternative in (2.14), and the total pressure from (2.15). It should be noted that all the other required thermodynamic functions, such as C(T), ρ*υ(T), etc., can be determined from E(T), once it is specified. If, for all T, E(T) is specified as Ew(T), the saturation vapor pressure over liquid water, then the effects of the latent heat of fusion are not included. If a synthesized E(T) is obtained from Ew(T) and Ei(T), the saturation vapor pressure over ice, then the model includes the effects of the latent heat of fusion. As discussed in O90, the synthesized E(T) involves two specified parameters: the center of the freezing zone Tf and the width of the freezing zone ΔTf . For our experiments with ice, we use Tf = 273.15 K and ΔTf = 1.0 K, with Ew(T) and Ei(T) obtained from the new Goff formulas (World Meteorological Organization 1979, appendix A).

The purpose of the precipitation microphysics diagnosis is to determine the precipitation fall velocity W and the precipitation source/sink term Qr. Our microphysics parameterization is identical to that used by Ooyama (2001), who modified the Klemp and Wilhelmson (1978) formulation to include ice. The equations for the precipitation microphysics diagnosis are

 
formula
 
formula
 
formula
 
formula
 
formula
 
formula
 
formula

where the constants ρr0, W0, τauto, τcol, τevap are given in appendix A. Note that the ice factor, defined as the ratio of ice to rain terminal velocity, is assumed to have the temperature dependence given by (2.16). When used in (2.18), this ice factor provides for a fivefold increase in the terminal velocity of precipitation upon descent through the melting layer. Note also that the precipitation source/sink term Qr is the sum of three terms: the autoconversion term Qauto, which initiates the precipitation process when the cloud condensate density ρc exceeds 0.001ρa; the collection term Qcol, which is a function of the cloud condensate and precipitation densities ρc and ρr; and the evaporation term Qevap, which is based on the analysis of the evaporation of a stationary drop and then augmented by the dimensionless ventilation factor fvent (with typical magnitude ∼10).

The purpose of the air–sea interaction diagnosis is to determine the boundary layer turbulent fluxes of water vapor, entropy, and radial and tangential momentum. We make the simple assumption that these fluxes vary linearly with height in the boundary layer, vanish at and above z = 1000 m, and have surface values determined by the bulk aerodynamic formulas.2 For example, for the momentum fluxes Fu and Fυ appearing in (2.5) and (2.6), the surface values are

 
formula

where the subscript zero denotes a surface value and V0 = (u20 + υ20)1/2 is the surface wind speed. Similarly, the surface values of the water vapor and entropy fluxes are

 
formula

where μ*υ0 and s(2)m0 are the saturation water vapor mixing ratio and the saturation entropy at the sea surface temperature and pressure. In all experiments reported here we have assumed a sea surface temperature of 28°C. We have also assumed that the drag and exchange coefficients are equal and depend on wind speed through Deacon’s formula

 
formula

where k = 4.0 × 105 m−1 s. It is well known (e.g., Malkus and Riehl 1960; Ooyama 1969; Emanuel 1995; Braun and Tao 2000) that the intensity of simulated tropical cyclones is sensitive to the relative magnitude of CD and CH. For example, if, in place of (2.25), one chooses CD = 0.0011 + kDV0 and CH = 0.0011 + kHV0, very intense storms are obtained when kD = 0 and kH = k, while much weaker storms are obtained when kD = k and kH = 0. Although such variations of kD and kH are extreme, it is an uncomfortable fact that the behavior of CD and CH at high wind speeds is one of the most uncertain aspects of all tropical cyclone models. The subject is one of active research (e.g., see Lighthill 1999; Emanuel 2003). Since our emphasis in this paper is on other aspects of the tropical cyclone problem, we use (2.25) in all the experiments reported here.

To summarize, the model consists of the seven prognostic Eqs. (2.1)(2.7) and the eighteen diagnostic Eqs. (2.8)(2.25). However, before discretization we make two further transformations. The first is motivated by the fact that direct use of the mixing ratio Eqs. (2.2) and (2.3) can lead to the problem of “negative water.” To prevent this, the model actually predicts ν = bhyp(μ) and νr = bhyp(μr), and then, for the purpose of thermodynamic diagnosis only, inverts these via μ = ahyp(ν) and μr = ahyp(νr), where bhyp and ahyp denote the biased hyperbolic transform and its quasi inverse, defined by

 
formula
 
formula

where the bias constant μ0 is set equal to 10−7. This procedure yields nonnegative values of μ and μr, which can be interpreted without difficulty by the thermodynamic diagnosis (2.8)(2.15). In the second modification made before discretization, the prognostic equations for ρa, ν, νr, s are rewritten in terms of deviations from a resting background state that depends on z only. For example, (2.1) becomes

 
formula

where a′ = aâ, with a = ln(ρa/ρa0), â = ln(ρ̂a/ρa0), and ρ̂a(z) the specified background dry air density. This procedure of predicting deviations from a background state enhances the numerical accuracy of the model.

The domain extends from the sea surface to a height of 24 km and from the vortex center to a radius of 1536 km. The outer boundary is open, and the solutions are assumed to decay exponentially with an e-folding distance of 1400 km. The top and bottom boundaries are assumed rigid, with the boundary condition w = 0. To reduce the effect of the reflection of gravity waves off the top boundary, we have also included Rayleigh-type damping terms in (2.4)(2.7). These terms damp the winds to zero and the specific entropy to its background value. The damping coefficient vanishes for z ≤ 18 km and increases linearly with height from 18 km to a maximum value of 0.015 s−1 at 24 km. Thus, the region between 18 km and the lid at 24 km is a sponge layer that effectively damps vertically propagating gravity waves that would otherwise reenter the region below 18 km after falsely reflecting off the lid. In all the cross sections of the present paper, only the region below 18 km is displayed.

The model spatial numerics are based on the CST method described in O2. As the name implies, the CST method uses the cubic B spline as the basis function. Because the first two derivatives of the B spline are continuous, the CST method has small computational dispersion errors, similar to the Fourier spectral method. Yet, because the B spline is locally defined, the CST method allows flexibility with regard to boundary conditions. With reduced dispersion errors and flexible boundary conditions, the CST method provides for noise-free nesting. To take advantage of this, the domain is discretized into a series of nested grids. As illustrated in Fig. 1, the horizontal domain consists of six grids. The horizontal grid spacing Δr within each grid increases by a factor of 2 from the finest grid at 0.5 km to the coarsest grid at 16.0 km. The vertical domain, in contrast, consists of a single grid with 48 grid intervals and a grid spacing Δz of 0.5 km. Time integration is accomplished in a two-stage process. In the first stage all the prognostic variables are advanced explicitly using the leapfrog scheme with a small enough time step for stability of gravity waves. In the second stage these explicit predictions are implicitly adjusted via a second-order elliptic equation. This particular implementation of the semi-implicit method allows a tenfold increase of the time step. The time steps are 2.5 s on the finest grid, 5 s on the next two grids, and 10 s on the coarsest three grids.

Fig. 1.

Vertical cross section of the initial tangential wind. Headings above this and the following figures indicate (left) the experiment, (center) the contoured variable, including units and contour interval Δ, and (right) the time in units of hh:mm:ss.s. Perturbation variables are identified by the (-bg) to the right of the variable name, indicating that a background state has been subtracted. The vertically oriented dashed lines mark the interfaces between nested grids, while the distances straddling these lines indicate the horizontal grid spacing to either side.

Fig. 1.

Vertical cross section of the initial tangential wind. Headings above this and the following figures indicate (left) the experiment, (center) the contoured variable, including units and contour interval Δ, and (right) the time in units of hh:mm:ss.s. Perturbation variables are identified by the (-bg) to the right of the variable name, indicating that a background state has been subtracted. The vertically oriented dashed lines mark the interfaces between nested grids, while the distances straddling these lines indicate the horizontal grid spacing to either side.

3. Potential vorticity equation and invertibility principle

Under the assumption of axisymmetry, the radial and vertical components of the absolute vorticity vector are (ξ, ζ) = (−∂m/rz, ∂m/rr) = [−∂υ/∂z, f + ∂()/rr], where m = + ½fr2 is the absolute angular momentum per unit mass. These are the only vorticity components that contribute to the potential vorticity. Since υ is the only velocity component appearing in the definitions of ξ and ζ, it follows that the only momentum equation involved in the PV derivation is (2.6). The derivation of the PV equation from (2.6) and the continuity equation for total density is given in appendix B. The result is

 
formula

where θρ is the virtual potential temperature and P is defined below in (3.5). The effects of θ̇ρ and appear in the first and second terms on the right-hand side of (3.1), with the former being interpreted as a measure of the variation of θ̇ρ along the m surfaces [or equivalently the variation of θ̇ρ along the projection of the vorticity vector onto the (r, z) plane], and with the latter being interpreted as a measure of the variation of along the θρ surfaces. The last term in (3.1) is generally small and describes the effects of precipitation and boundary layer water vapor flux on the potential vorticity.

The variable P has two important properties that contribute to its fundamental importance (Schubert 2004): (i) it reduces to the classical Ertel PV in the limit of a completely dry atmosphere; (ii) it is invertible, that is, it carries all the dynamical information about the balanced part of the wind and mass fields. Other choices for the scalar field in the definition of PV are lacking in this regard. For example, if saturation equivalent potential temperature is used in place of θρ, property (i) is lost, while if equivalent potential temperature is used in place of θρ, property (ii) is lost.

To see that an invertibility principle exists for P, consider the situation in which the radial momentum Eq. (2.5) and the vertical momentum Eq. (2.7) reduce to the gradient wind Eq. (3.2) and the hydrostatic Eq. (3.3), respectively. Then, including the definition of virtual potential temperature (3.4) and the definition of PV (3.5), we have the system

 
formula
 
formula
 
formula
 
formula

This is a system of four equations for the four unknowns υ, ρ, p, θρ, with given P. By eliminating υ, ρ, and θρ we can obtain a single partial differential equation relating the pressure p to the potential vorticity P. Although this equation is complicated in the z coordinate, simpler equivalent forms can be obtained by using either (p/p0)κ or θρ as the vertical coordinate (Schubert et al. 2001). Note that the solution of the invertibility problem gives us the total density ρ, the total pressure p, and the virtual potential temperature θρ, that is, the parts of the mass field that are of direct dynamical significance. It can be shown (Schubert 2004) that the moist invertibility principle (3.2)(3.5) is isomorphic with the dry invertibility principle solved in previous studies (e.g., Schubert and Alworth 1987) under the interchanges ρaρ, pap, and θρθ. Thus, previous studies of the dry invertibility principle are easily interpreted in terms of the moist invertibility principle.

In this paper we will not be concerned with actually solving the invertibility principle (3.2)(3.5). Rather, we simply use the existence of the principle as justification for studying model output fields of P. In section 4 we will use a transformed version of (3.1) to understand the PV dynamics of a quasi-steady tropical cyclone. When attempting to obtain physical understanding of certain phenomena, it has been said (Stommel 1995) that complex models such as (2.1)(2.27) are “not much help: like vegetable soup, they have too many ingredients to reveal which one imparts the flavor.” In the case of a tropical cyclone, it is definitely the PV that imparts the flavor.

4. Control experiment

The initial condition for all experiments is a purely azimuthal vortex in gradient and hydrostatic balance, with a horizontally uniform relative humidity field. At the surface the radial distribution of initial azimuthal velocity is 2υm(r/rm)/[1 + (r/rm)2], with the maximum wind specified as υm = 12 m s−1 and the radius of maximum wind as rm = 100 km. The initial azimuthal velocity is assumed to decrease linearly with height to zero at z = 18 km, and to be zero between 18 km and the model top at 24 km. Figure 1 shows this initial azimuthal wind field. The associated pressure anomaly field has a minimum of −7.1 hPa at the surface in the vortex core. The far-field temperature and humidity is taken from Jordan’s (1958) mean hurricane season sounding and is shown in Fig. 2. At the initial time this sounding is approximately valid at all radii, since the initial horizontal temperature gradient is weak, with the lower tropospheric vortex core approximately 1.1 K warmer than the far field.

Fig. 2.

Initial far-field soundings of temperature (solid) and dewpoint temperature (dashed).

Fig. 2.

Initial far-field soundings of temperature (solid) and dewpoint temperature (dashed).

To summarize the temporal and spatial evolution of this control (CNTL) experiment, Fig. 3 depicts the 3-h running mean surface tangential wind speed as a function of radius and time. Between 40 and 50 h, a very small, intense central vortex develops. As discussed in appendix C, this small, central vortex is not the result of a problem with CST numerics. In fact such small, central vortices also appear in previous high resolution axisymmetric tropical cyclone simulations (e.g., Yamasaki 1983; Willoughby et al. 1984) that use entirely different numerics. Although strong cumulus convection (referred to as a hub cloud) is sometimes observed in the center of a hurricane eye, the convection appearing near r = 0 in our CNTL experiment is unrealistically deep and intense. Beyond 80 h, however, this central vortex weakens substantially, while a secondary tangential wind maximum forms at about 30 km and propagates inward, intensifying along the way. This wind maximum is collocated with a ring of convection that eventually contracts to become the new tropical cyclone eyewall. From 100 to 180 h, the tropical cyclone gradually intensifies. As shown in Fig. 3, this overall increase of intensity is not steady but highly variable. For instance, during the 12-h period following 118 h, the maximum tangential wind oscillates over 20 m s−1. This variability, which is very similar to that obtained by Willoughby et al. (1984), is caused by disturbances of the boundary layer inflow and the formation of secondary rings of convection that propagate inward. The secondary circulation of a new, inward-propagating ring and the restriction of high θe inflow into the core, results in the dissipation of the primary ring, and its ultimate replacement by the secondary ring.

Fig. 3.

Isolines of the 3-h-averaged surface tangential wind speed as a function of (r, t) from 0 to 100 km and 0 to 240 h for the CNTL experiment. The contour interval is 15 m s−1, starting from the 10 m s−1 contour.

Fig. 3.

Isolines of the 3-h-averaged surface tangential wind speed as a function of (r, t) from 0 to 100 km and 0 to 240 h for the CNTL experiment. The contour interval is 15 m s−1, starting from the 10 m s−1 contour.

Beyond 180 h, the vortex settles into a quasi steady state. Figures 4a–c depict the 180–240-h average primary and secondary circulations. The model steady-state cyclone is similar in structure to observed storms but is more intense. The tangential wind has a maximum of over 90 m s−1, with the peak located beneath the eyewall at a radius of 14 km. This intense vortex is produced by the advection of high angular momentum air toward the center. This same air is later advected away from the center within the outflow branch of the secondary circulation, producing an anticyclonic vortex with a wind speed of −30 m s−1 at a radius of 1250 km and a height of 14 km. The horizontal branches of the secondary circulation are concentrated into shallow layers at the surface and near the tropopause. As shown in Fig. 4b, the inflow branch is confined3 below 2 km, with a maximum inward radial flow of 31.8 m s−1 at r = 21 km. The maximum inflow lies just outside the radius of maximum tangential wind. The most intense outflow is confined to a layer between 12 and 18 km. The maximum outflow of 29.5 m s−1 is located at a radius of about 75 km, which is more than 50 km away from the eyewall. The two vertical branches of the secondary circulation are very different in terms of intensity and horizontal scale. Immediately above and sloping away from the region of maximum surface convergence is the ascending branch of the secondary circulation, embedded within the eyewall cloud. In the ascending branch, vertical velocities reach 5.5 m s−1 just above the boundary layer and in the upper troposphere. The width of the ascending branch increases with height. In contrast, the compensating subsidence in the descending branch of the secondary circulation is very weak and extends into the far field of the domain. The transverse circulation depicted in Fig. 4 provides a deformation field in the (r, z) plane, with especially large deformation in the lower troposphere at radii between 10 and 20 km. This frontogenetic effect (Eliassen 1959; Emanuel 1997) acts to crowd together the absolute angular momentum surfaces, resulting in large vorticity.

Fig. 4.

Time-averaged (180–240 h) cross sections of the (a) tangential, (b) radial, and (c) vertical wind speeds (m s−1) for the CNTL experiment. Solid curves indicate positive values of u, υ, or w, while the dotted curves indicate zero or negative values. Shading in this figure and the following two figures denotes the region where the airborne condensate mixing ratio is greater than 0.1 g kg−1.

Fig. 4.

Time-averaged (180–240 h) cross sections of the (a) tangential, (b) radial, and (c) vertical wind speeds (m s−1) for the CNTL experiment. Solid curves indicate positive values of u, υ, or w, while the dotted curves indicate zero or negative values. Shading in this figure and the following two figures denotes the region where the airborne condensate mixing ratio is greater than 0.1 g kg−1.

Figure 5 shows the 180–240-h mean cross sections of T and the water vapor mixing ratio μυ. The maximum T ′, disregarding the effect of the very small central vortex, is approximately 17 K and is located at a height of 12 km. The axis of maximum T ′ extends both outward into the stratiform precipitation and downward along the inner edge of the eyewall. The horizontal extent of these temperature anomalies depends on the local radius of deformation, defined by (gH)1/2{[ f + ∂()/rr]( f + 2υ/r)}−1/2, where H is the equivalent depth and where the factor in the denominator is a measure of the inertial stability or “stiffness” of the vortex. When the radius of deformation is small, the vortex is stiffened such that the horizontal extent of the secondary circulation is restricted. Assuming (gH)1/2 ≈ 60 m s−1, the radius of deformation is less than 5 km along the inner edge of the eyewall; however, within the outflow, it is larger than 300 km. Therefore, we observe a narrow secondary circulation and adiabatic warming on the inner edge of the eyewall, as compared to a broad secondary circulation outside the eyewall.

Fig. 5.

Time-averaged (180–240 h) cross section of the (a) temperature deviation (temperature minus height-dependent background temperature) and (b) water vapor mixing ratio for the CNTL experiment.

Fig. 5.

Time-averaged (180–240 h) cross section of the (a) temperature deviation (temperature minus height-dependent background temperature) and (b) water vapor mixing ratio for the CNTL experiment.

The dynamics of the eye and eyewall also have a distinct influence on the distribution of water vapor, as shown in Fig. 5b. Radially, the water vapor mixing ratio μυ is a maximum in the eyewall, because of vertical advection within the updraft. Because of subsidence within the eye, μυ is a minimum along the inside of the eyewall. At the surface, the large flux of water vapor increases μυ outside 25 km to 23.9 g kg−1, which is 5.5 g kg−1 greater than the initial value. The corresponding surface relative humidity is nearly 100%. Beneath the eyewall, downdrafts produce a local surface minimum of 21.6 g kg−1 at a radius of 14 km. Above the 0°C isotherm, the terminal velocity of precipitation is approximately 1.5–2 m s−1, whereas below this level, it increases to as much as 8 m s−1 within the eyewall. Because of the relatively small terminal velocity aloft and the intense outflow, precipitation is advected far from the eyewall, resulting in surface precipitation rates of about 10–40 mm h−1. However, beneath the eyewall, W and μr are both large, producing a precipitation rate of over 200 mm h−1.

Figure 6 shows 180–240-h mean cross sections of the vertical component of the absolute vorticity ζ, the virtual potential temperature θρ, and the potential vorticity anomaly P′ = PP, where the far-field potential vorticity is defined by P = ( f/)(dρ/dz). Since the horizontal shear across the inner edge of the eyewall is very large in the middle troposphere, with tangential wind speeds increasing by 60 m s−1 over only 5 km radius, there is a midtropospheric maximum of ζ and an associated maximum of P′ ≈ 275 PVU. In contrast, for the upper-tropospheric maximum of P′ ≈ 275 PVU the (∂υ/∂z) (∂θρ/∂r) term in (3.5) becomes more important because of the large horizontal θρ gradient resulting from the warm core of the vortex and the relatively large vertical shear of the tangential wind beneath the tropopause. In addition there is enhanced PV along the central axis, a remnant of the small central vortex. If the r axis in Fig. 6c were stretched to correspond to the z axis, the outward slope of the large PV on the inner edge of the eyewall cloud would become more apparent, and the structure would resemble a leaning tower of high PV. In three dimensions it might be viewed as a bowl of high PV.

Fig. 6.

Time-averaged (180–240 h) cross sections of (a) vertical component of absolute vorticity, (b) virtual potential temperature, and (c) potential vorticity anomaly for the CNTL experiment.

Fig. 6.

Time-averaged (180–240 h) cross sections of (a) vertical component of absolute vorticity, (b) virtual potential temperature, and (c) potential vorticity anomaly for the CNTL experiment.

To understand the origin of this remarkable PV structure, let us return to (3.1), neglecting the generally small last term. To rewrite this equation in a more physically revealing form, first define the potential radius R by ½fR2 = m = + ½fr2, so that fRṘ = where = DR/Dt. Transforming from (r, z, t) to (R, θρ, 𝒯), where 𝒯 = t, the first two terms on the right-hand side of (3.1) can be written as

 
formula
 
formula

and the material derivative operator as

 
formula

Note that ∂/∂R implies fixed θρ, ∂/∂θρ implies fixed R, and ∂/∂𝒯 implies fixed R, θρ. Using (4.1)(4.3) in (3.1), we obtain

 
formula

In the steady state and above the frictional boundary layer we have ∂P/∂𝒯 = 0 and = 0, so that the potential vorticity Eq. (4.4) reduces to θ̇ρ(∂P/∂θρ) = P(∂θ̇ρ/∂θρ), which can also be written as

 
formula

Thus, in the steady state above the boundary layer, the ratio of θ̇ρ to P is constant along an absolute angular momentum surface. For an angular momentum surface erupting from the boundary layer into the eyewall updraft we can assume that the integration constant required by (4.5) is set by P(R, θρB) and θ̇ρ(R, θρB), the values of P and θ̇ρ at the top of the boundary layer. Then, integration of (4.5) yields

 
formula

Thus, in the steady state there is an intimate connection between the P field and the θ̇ρ field. In essence, the steady-state P field has become locked to the θ̇ρ field.

Since the P field determines the primary circulation and the θ̇ρ field, along with frictional effects, determine the secondary circulation, (4.6) can be interpreted as the fundamental relation between the primary and secondary parts of the steady-state circulation.

In the steady state, a parcel of eyewall air erupting from the boundary layer on a given angular momentum surface stays on this same outward-sloping surface as it spirals upward, crossing θρ surfaces and changing its PV at the rate P(∂θ̇ρ/∂θρ), which is positive below the level of maximum θ̇ρ and negative above this level. Note that the level of maximum θ̇ρ is also the level of maximum P through the coupling described by (4.6). All parcels erupting from the boundary layer on the same angular momentum surface have the same Lagrangian history, but the Lagrangian histories are generally distinct on different angular momentum surfaces. Unfortunately, when applied to the CNTL experiment, the above argument is compromised by the unsteadiness illustrated in Fig. 3. However, as we shall see, the argument more accurately holds for the “no ice experiment” discussed in section 5.

It is interesting to note that the generalized moist PV, defined by (3.5), is only a modest modification of the dry Ertel PV, which is obtained by replacing ρ with ρa and θρ with θ in (3.5). This modest modification is in sharp contrast to other proposed PV generalizations that involve use of the equivalent potential temperature or the saturation equivalent potential temperature, both of which result in drastic differences with the dry Ertel PV. To confirm that (3.5) yields PV fields that are very similar to dry Ertel PV fields, we have produced a 180–240-h mean cross section of the dry Ertel PV. This dry PV cross section (not shown) is nearly identical to Fig. 6c. Thus, the dry Ertel PV is useful for diagnostic analysis of moist models (e.g., Yau et al. 2004; Braun et al. 2006). In addition, although a hurricane is definitely a phenomenon involving moist physics, a reasonable approximation to the balanced dynamics can be constructed (e.g., Schubert and Alworth 1987; Möller and Smith 1994) using the dry PV invertibility principle and the dry PV evolution equation, as long as diabatic and frictional effects are properly included in the latter equation.

5. Sensitivity experiments

a. Sensitivity to ice microphysics

In our representation of thermodynamics and microphysics there are two effects of ice: (i) the latent heat effects of ice are incorporated by computing E(T), from which L(T), C(T), D(T), ρ*υ(T) are determined, as a temperature-dependent interpolation of the saturation vapor pressure over a plane surface of water and a plane surface of ice; (ii) the reduced settling speed of geometrically complex ice particles is included in the model through the temperature-dependent ice factor defined by (2.16). Remembering that our CNTL experiment included both these effects of ice, we now perform a sensitivity experiment with no latent heat of fusion and with the ice factor set to unity.4

Figure 7 shows the 3-h-averaged surface tangential wind as a function of r, t for this “no ice” (NICE) experiment, while Figs. 8 –10 show the quasi-steady-state structure in the (r, z) plane. The formats of Figs. 8 –10 are identical to Figs. 4 –6 except that the fields are averaged from 120 to 240 h. Comparing Fig. 7 with Fig. 3, we see that in the NICE experiment the simulated tropical cyclone develops more rapidly and attains a more intense steady state than in the CNTL experiment. In addition the NICE experiment is distinctly less variable. Comparison of Figs. 8 –10 with Figs. 4 –6 indicates that the NICE experiment has stronger tangential winds, a smaller radius of maximum tangential winds, a low-level inflow penetrating further inward, a narrower and more intense updraft, a warmer and dryer central core, a thinner annular ring of eyewall cloud, and a narrower and more intense outward tilting region of maximum ζ and P. In particular, note that the peak values of P in the NICE experiment are approximately 400 PVU, which is even more extreme than those found in the CNTL experiment. Note also that our 0.5 km by 0.5 km resolution on the inner grid is required to accurately capture these extreme features of the ζ and P fields. These differences between the CNTL and NICE experiments are caused by the development and contraction of secondary eyewalls in the CNTL experiment, features that do not develop in the NICE experiment. The reason that secondary eyewalls do not form in the NICE experiment is that extensive stratiform precipitation does not develop. Above 5 km the downward terminal velocity of precipitation reaches 10 m s−1 in the NICE experiment and only 2 m s−1 in the CNTL experiment. Even though the CNTL experiment has a weaker time-averaged secondary circulation than the NICE experiment, the small terminal velocities in the CNTL experiment allow the precipitation to be lofted into the outflow and advected far from the eyewall. However, in the NICE experiment, much of the precipitation falls from the sloping updraft without being ejected into the outflow. As a result, the surface precipitation extends to almost 100 km in the CNTL experiment but only about 45 km in the NICE experiment. Furthermore, since the precipitation is distributed over a smaller area in the NICE experiment, the precipitation rate is a factor of four greater than in the CNTL experiment. Without the stratiform precipitation to provide the mesoscale downdrafts that induce surface convergence and ascent, secondary eyewalls do not form in the NICE experiment.

Fig. 7.

Same as Fig. 3, but for the NICE experiment.

Fig. 7.

Same as Fig. 3, but for the NICE experiment.

Fig. 8.

Same as Fig. 4 but for 120–240-h time average of the NICE experiment.

Fig. 8.

Same as Fig. 4 but for 120–240-h time average of the NICE experiment.

Fig. 10.

Same as Fig. 6 but for 120–240-h time average of the NICE experiment.

Fig. 10.

Same as Fig. 6 but for 120–240-h time average of the NICE experiment.

These differences between the CNTL and NICE experiments are consistent with the results reported by Willoughby et al. (1984), Lord et al. (1984), and Lord and Lord (1988) using a nonhydrostatic moist model with 1-km vertical grid spacing and 2-km radial grid spacing in the inner core, but with a more elaborate microphysical parameterization. To compare our model with the Willoughby–Lord model, it is useful to measure the complexity of a nonhydrostatic moist model by the number of prognostic equations over and above the five required for a dry model. For the nonhydrostatic moist model (2.1)(2.25) the count is two, that is, (2.2) for the total water mixing ratio μ and (2.3) for the precipitation mixing ratio μr. In fact, this nonhydrostatic moist model can be considered the model of maximum simplicity among the class of models that explicitly calculate the movement of solid and liquid precipitation from its formation to its impact with the earth’s surface. Contributing to the simplicity of this model is the fact that the distribution of precipitation is described by a single scalar field μr, even though the precipitation can be solid or liquid. In contrast, the count for the Willoughby–Lord model is six, one for each of the mixing ratios of water vapor, cloud water, cloud ice, rain, snow, and graupel. With the increased number of prognostic equations comes an increased number of conversion processes. For example, with the model of maximum simplicity, there are only three parameterized conversions, Qauto, Qcol, Qevap, as given by (2.19)(2.21). In the Willoughby–Lord model the number of parameterized conversions increases to twenty-four, which requires a more extensive theoretical and observational basis for the microphysical parameterization. Thus, it is interesting that very similar results to those of obtained by the Willoughby–Lord model can be obtained with a microphysical model of maximum simplicity.

In a second sensitivity experiment (not shown) we included the latent heat effect of ice but set fice to unity. This yields results qualitatively similar to the NICE experiment. In both of these experiments with fice = 1, the precipitation is not advected away from the eyewall but rapidly falls from the sloping updraft. Without extensive stratiform precipitation and the resulting horizontal dipole of freezing and melting, secondary eyewalls do not form; thus, the latent effects of ice alone cannot explain the significant differences between the NICE and CNTL experiments. In a third sensitivity experiment (not shown) we excluded the latent heat effect of ice but left the ice factor as defined by (2.16). This experiment yields results qualitatively similar to the CNTL experiment. Thus, our results indicate that the reduced terminal fall velocity associated with frozen precipitation is the primary factor responsible for the qualitative differences between the CNTL and NICE experiments. The latent effects of ice are a secondary effect. Thus, using a simple microphysical scheme, we have obtained results consistent with those obtained using more sophisticated parameterizations (Willoughby et al. 1984; Lord et al. 1984; Lord and Lord 1988).

b. Sensitivity to precipitation effects

One distinctive feature of our model is the inclusion in (2.4)(2.7) of terms representing the vertical fluxes of entropy and momentum by precipitation. To evaluate the sensitivity of simulated tropical cyclones to these vertical fluxes by precipitation, we have performed two additional experiments, the first of which includes the W term in (2.4) but excludes the W terms in (2.5)(2.7). The results of this experiment (not shown) demonstrate that the precipitation terms in the momentum equations have relatively little impact on tropical cyclone development. This result is not unexpected considering that the momentum of the precipitation, ρrW, is typically much smaller than ρw.

The final sensitivity experiment excludes the W term in (2.4) but includes the W terms in (2.5)(2.7). The results of this experiment (not shown) differ substantially from the CNTL. For instance, the tropical cyclone intensifies more rapidly and attains a more intense steady state than the CNTL experiment. These differences result from changes in the boundary layer entropy, with the θe of inflowing air being 10 K warmer beneath the eyewall compared to the CNTL experiment. With higher θe in the boundary layer, the development of secondary eyewalls is inhibited. We conclude that the vertical flux of entropy by precipitation is an important effect that should be included in tropical cyclone models. Perhaps this is not surprising, since it is well known that accurate simulation of the boundary layer moist entropy budget is an important part of the tropical cyclone forecasting problem.

6. Concluding remarks

We have presented results from an axisymmetric, nonhydrostatic, full-physics model of the tropical cyclone and have used the moist potential vorticity principle as a diagnostic for interpreting the quasi-steady-state structure. This analysis demonstrates how the P and θ̇ρ fields become locked together in a thin leaning tower on the inner edge of the eyewall cloud. In addition, our sensitivity experiments indicate that accurate simulations of tropical cyclones require 1) the inclusion of ice and, in particular, the slow terminal velocities associated with frozen precipitation above the melting level; 2) the inclusion of the vertical transport of moist entropy by precipitation.

Since the numerical model used in this study is based on the equation set (2.1)(2.25), it is interesting to comment on the usefulness of these equations as a physical model of a tropical cyclone. We certainly have the most confidence in the prognostic Eqs. (2.1)(2.7) and the equilibrium thermodynamics (2.8)(2.15) as being part of an accurate description of nature. Although precipitation microphysics has been parameterized as a bulk process, (2.16)(2.18) and (2.21) have considerable observational and laboratory support (e.g., see Kessler 1969). However, (2.19) and (2.20) must be regarded as a crude parameterization of the intricate process of precipitation formation. As for the air–sea interaction parameterization (2.23)(2.25), the most uncertain aspect is (2.25), especially the values of CD and CH at high wind speeds. Finally, a major limitation of the model (2.1)(2.25) is the assumption of axisymmetry.

Consistent with the results of Persing and Montgomery (2003), our results show that the quasi-steady-state intensity for an axisymmetric storm may be much greater than predicted by energetically based, steady-state maximum potential intensity (MPI) theories. To explain this discrepancy, Persing and Montgomery have used the nonhydrostatic model developed by Rotunno and Emanuel (1987) to show that past model simulations supporting such MPI theories may have been limited by their coarse resolution and large diffusion—problems that have been overcome in both the Persing–Montgomery simulations and those presented here. At the same time it should be realized that high resolution, low diffusion, axisymmetric models probably produce overly intense vortices because of the suppression of variations in the azimuthal direction. In fact, the PV distributions produced in such models would be unstable in a combined baroclinic–barotropic sense if the models allowed variations in the azimuthal direction. Over their life cycle such instabilities would radially mix the PV, thereby reducing the maximum tangential wind (Schubert et al. 1999). In addition, although difficult to observe and thus inadequately documented, Kelvin–Helmholtz instability along the sloping inner edge of the eyewall may also play a role in limiting the maximum tangential winds.

Although the PV mixing process can be crudely parameterized in an axisymmetric model by using radial diffusion or hyperdiffusion, there are unrealistic aspects associated with this type of parameterization (Kossin and Schubert 2003). In this regard it should be recalled that our model has no explicit frictional effects above the boundary layer, although the CST method does include a third-order derivative constraint in the least squares minimization that defines the transform of spatial fields to nodal amplitudes. This third order derivative constraint is equivalent, in wavenumber space, to a sharp, sixth-order, low-pass filter that effectively eliminates small-scale errors at the resolution limit. This filter is weak compared to the usual diffusion or hyperdiffusion used in most finite difference models, so that the model flow above the boundary layer should be regarded as quite inviscid. Thus, the assumption of axisymmetry, in conjunction with the nearly inviscid nature of the flow above the boundary layer, probably leads to the overly intense vortices produced by the model.

Even with the above limitations in mind, our results help answer the question, “What is a quasi-steady-state hurricane?” It is no doubt a complicated structure involving all three flow components and many moist thermodynamic fields. But, at its core, it is an extreme structure in which the P field and the θ̇ρ field have become intimately coupled in such a way that they vary in a similar fashion along a tightly packed group of absolute angular momentum surfaces. To capture the formation and asymmetric evolution of such extreme PV structures in full-physics 3D models is obviously a very challenging problem, but it may be necessary for accurate intensity forecasts. Perhaps this is part of a sobering realization that hurricane intensity forecasting is fundamentally much more difficult than hurricane track forecasting.

Fig. 9.

Same as Fig. 5 but for 120–240-h time average of the NICE experiment.

Fig. 9.

Same as Fig. 5 but for 120–240-h time average of the NICE experiment.

Fig. C1. Results of shallow water experiments on the Rayleigh recoil problem under four different dynamics: axisymmetric, nonlinear; axisymmetric, linear; slab symmetric, nonlinear; slab symmetric, linear. (upper) The fluid depth h (m) at r = 0 as a function of time (in seconds), and (lower) h as function of r (in kilometers) at the time of maximum fluid depth at r = 0.

Fig. C1. Results of shallow water experiments on the Rayleigh recoil problem under four different dynamics: axisymmetric, nonlinear; axisymmetric, linear; slab symmetric, nonlinear; slab symmetric, linear. (upper) The fluid depth h (m) at r = 0 as a function of time (in seconds), and (lower) h as function of r (in kilometers) at the time of maximum fluid depth at r = 0.

Acknowledgments

The authors thank Paul Ciesielski, William Cotton, Rob Fleishauer, Matthew Garcia, Richard Johnson, James Kossin, John McGinley, Brian McNoldy, Michael Montgomery, John Persing, Richard Taft, Jonathan Vigh, and three anonymous reviewers for helpful comments. This work was supported by NASA/CAMEX Grant NAG5-11010, NASA/TCSP Grant 04-0007-0031, NSF Grant ATM-0332197, and NOAA Grant NA17RJ1228.

REFERENCES

REFERENCES
Braun
,
S. A.
,
2002
:
A cloud-resolving simulation of Hurricane Bob (1991): Storm structure and eyewall buoyancy.
Mon. Wea. Rev.
,
130
,
1573
1592
.
Braun
,
S. A.
, and
W-K.
Tao
,
2000
:
Sensitivity of high-resolution simulations of Hurricane Bob (1991) to planetary boundary layer parameterizations.
Mon. Wea. Rev.
,
128
,
3941
3961
.
Braun
,
S. A.
,
M. T.
Montgomery
, and
Z.
Pu
,
2006
:
High-resolution simulation of Hurricane Bonnie (1998). Part I: The organization of eyewall vertical motion.
J. Atmos. Sci.
,
63
,
19
42
.
Eliassen
,
A.
,
1959
:
On the formation of fronts in the atmosphere.
The Rossby Memorial Volume, Rockefeller Institute and Oxford University, 277–287
.
Emanuel
,
K. A.
,
1995
:
Sensitivity of tropical cyclones to surface exchange coefficients and a revised steady-state model incorporating eye dynamics.
J. Atmos. Sci.
,
52
,
3969
3976
.
Emanuel
,
K. A.
,
1997
:
Some aspects of hurricane inner-core dynamics and energetics.
J. Atmos. Sci.
,
54
,
1014
1026
.
Emanuel
,
K. A.
,
2003
:
A similarity hypothesis for air–sea exchange at extreme wind speeds.
J. Atmos. Sci.
,
60
,
1420
1428
.
Hausman
,
S. A.
,
2001
:
Formulation and sensitivity analysis of a nonhydrostatic, axisymmetric tropical cyclone model. Atmospheric Science Paper 701, Colorado State University, 210 pp
.
Jordan
,
C. L.
,
1958
:
Mean soundings for the West Indies area.
J. Meteor.
,
15
,
91
97
.
Kessler
,
E.
,
1969
:
On the Distribution and Continuity of Water Substance in Atmospheric Circulations. Meteor. Monogr., No. 32, Amer. Meteor. Soc., 84 pp
.
Klemp
,
J. B.
, and
R.
Wilhelmson
,
1978
:
The simulation of three-dimensional convective storm dynamics.
J. Atmos. Sci.
,
35
,
1070
1096
.
Kossin
,
J. P.
, and
W. H.
Schubert
,
2003
:
Diffusion versus advective rearrangement of a circular vortex sheet.
J. Atmos. Sci.
,
60
,
586
589
.
Lighthill
,
J.
,
1999
:
Ocean spray and the thermodynamics of tropical cyclones.
J. Eng. Math.
,
35
,
11
42
.
Lord
,
S. J.
, and
J. M.
Lord
,
1988
:
Vertical velocity in an axisymmetric, nonhydrostatic tropical cyclone model.
J. Atmos. Sci.
,
45
,
1453
1461
.
Lord
,
S. J.
,
H. E.
Willoughby
, and
J. M.
Piotrowicz
,
1984
:
Role of a parameterized ice phase microphysics in an axisymmetric, nonhydrostatic tropical cyclone model.
J. Atmos. Sci.
,
41
,
2836
2848
.
Malkus
,
J. S.
, and
H.
Riehl
,
1960
:
On the dynamics and energy transformations in steady state hurricanes.
Tellus
,
12
,
1
20
.
Manzello
,
S. L.
, and
J. C.
Yang
,
2002
:
An experimental study of a water droplet impinging on a liquid surface.
Exp. Fluids
,
32
,
580
589
.
Möller
,
J. D.
, and
R. K.
Smith
,
1994
:
The development of potential vorticity in a hurricane-like vortex.
Quart. J. Roy. Meteor. Soc.
,
120
,
1255
1265
.
Ooyama
,
K. V.
,
1969
:
Numerical simulation of the life cycle of tropical cyclones.
J. Atmos. Sci.
,
26
,
3
40
.
Ooyama
,
K. V.
,
1990
:
A thermodynamic foundation for modeling the moist atmosphere.
J. Atmos. Sci.
,
47
,
2580
2593
.
Ooyama
,
K. V.
,
2001
:
A dynamic and thermodynamic foundation for modeling the moist atmosphere with parameterized microphysics.
J. Atmos. Sci.
,
58
,
2073
2102
.
Ooyama
,
K. V.
,
2002
:
The cubic-spline transform method: Basic definitions and tests in a 1D single domain.
Mon. Wea. Rev.
,
130
,
2392
2415
.
Persing
,
J.
, and
M. T.
Montgomery
,
2003
:
Hurricane superintensity.
J. Atmos. Sci.
,
60
,
2349
2371
.
Rein
,
M.
,
1996
:
The transitional regime between coalescing and splashing drops.
J. Fluid Mech.
,
306
,
145
165
.
Rogers
,
R.
,
S.
Chen
,
J.
Tenerelli
, and
H.
Willoughby
,
2003
:
A numerical study of the impact of vertical shear on the distribution of rainfall in Hurricane Bonnie (1998).
Mon. Wea. Rev.
,
131
,
1577
1599
.
Rotunno
,
R.
, and
K. A.
Emanuel
,
1987
:
An air–sea interaction theory for tropical cyclones. Part II: Evolutionary study using a nonhydrostatic axisymmetric numerical model.
J. Atmos. Sci.
,
44
,
542
561
.
Schubert
,
W. H.
,
2004
:
A generalization of Ertel’s potential vorticity to a cloudy, precipitating atmosphere.
Meteor. Z.
,
13
,
465
471
.
Schubert
,
W. H.
, and
B. T.
Alworth
,
1987
:
Evolution of potential vorticity in tropical cyclones.
Quart. J. Roy. Meteor. Soc.
,
113
,
147
162
.
Schubert
,
W. H.
,
M. T.
Montgomery
,
R. K.
Taft
,
T. A.
Guinn
,
S. R.
Fulton
,
J. P.
Kossin
, and
J. P.
Edwards
,
1999
:
Polygonal eyewalls, asymmetric eye contraction, and potential vorticity mixing in hurricanes.
J. Atmos. Sci.
,
56
,
1197
1223
.
Schubert
,
W. H.
,
S. A.
Hausman
,
M.
Garcia
,
K. V.
Ooyama
, and
H-C.
Kuo
,
2001
:
Potential vorticity in a moist atmosphere.
J. Atmos. Sci.
,
58
,
3148
3157
.
Stommel
,
H. M.
,
1995
:
Autobiography.
Collected Works of Henry M. Stommel, N. G. Hogg and R. N. Huang, Eds., Vol. 1, Amer. Meteor. Soc. 380 pp
.
Willoughby
,
H. E.
,
H-L.
Jin
,
S. J.
Lord
, and
J. M.
Piotrowicz
,
1984
:
Hurricane structure and evolution as simulated by an axisymmetric nonhydrostatic numerical model.
J. Atmos. Sci.
,
41
,
1169
1186
.
World Meteorological Organization
,
1979
:
General. Vol. I, Technical Regulations, WMO 49, 3 pp
.
Yamasaki
,
M.
,
1983
:
A further study of the tropical cyclone without parameterizing the effects of cumulus convection.
Pap. Meteor. Geophys.
,
34
,
221
260
.
Yau
,
M. K.
,
Y.
Liu
,
D-L.
Zhang
, and
Y.
Chen
,
2004
:
A multiscale numerical study of Hurricane Andrew (1992). Part VI: Small-scale inner-core structures and wind streaks.
Mon. Wea. Rev.
,
132
,
1410
1433
.
Zhang
,
D-L.
, and
X.
Wang
,
2003
:
Dependence of hurricane intensity and structures on vertical resolution and time-step size.
Adv. Atmos. Sci.
,
20
,
711
725
.
Zhang
,
D-L.
,
Y.
Liu
, and
M. K.
Yau
,
2000
:
A multiscale numerical study of Hurricane Andrew (1992). Part III: Dynamically induced vertical motion.
Mon. Wea. Rev.
,
128
,
3772
3788
.

APPENDIX A

List of Symbols

Mass densities, mixing ratios, temperatures, pressures, velocities
  • ρa  mass density of dry air

  • ρυ  mass density of water vapor

  • ρc  mass density of airborne condensate

  • ρr  mass density of precipitating water substance

  • ρm = ρυ + ρc  mass density of airborne moisture

  • ρ = ρa + ρm + ρr  total mass density

  • μm = ρm/ρa  mixing ratio of airborne moisture

  • μr = ρr/ρa  mixing ratio of precipitating water substance

  • μ = μm + μr  mixing ratio of total water substance

  • T1  temperature for thermodynamic state 1

  • T2  temperature for thermodynamic state 2

  • T = max(T1, T2)  temperature

  • Tρ = p/(ρRa)  virtual temperature

  • θρ = Tρ(p0/p)κ  virtual potential temperature

  • θ = T(p0/pa)κ  dry potential temperature

  • θe = T0 exp[σ/(ρacpa)]  equivalent potential temperature

  • pa  partial pressure of dry air

  • pυ  partial pressure of water vapor

  • p = pa + pυ  total pressure of moist air

  • u, υ, w  radial, azimuthal, and vertical components of the velocity of dry air and airborne moisture

  • W  vertical velocity of precipitation (relative to air)

  • w = w + ρ−1 (ρrW + Fm)  density-weighted-mean vertical velocity

Specific entropies (J kg−1 K−1) and entropy densities (J m−3 K−1)
  • sa(ρa, T)  specific entropy of dry air, defined by sa(ρa, T) = cυa ln(T/T0) − Ra ln(ρa/ρa0)

  • s(1)m(ρm, T)  specific entropy of airborne moisture in state 1: s(1)m(ρm, T) = cυυ ln(T/T0) − Rυ ln(ρm/ρ*υ0) + L(T0)/T0

  • s(2)m(ρm, T)  specific entropy of airborne moisture in state 2: s(2)m(ρm, T) = C(T) + D(T)/ρm

  • sr = C(T2)  specific entropy of condensed water

  • s = σ/ρa  dry-air-specific entropy of moist air

  • σa = ρasa  entropy density of dry air

  • σm = ρmsm  entropy density of airborne water substance

  • σr = ρrsr  entropy density of precipitating water substance

  • σ = σa + σm + σr  total entropy density

  • S1(ρa, ρm, T)  entropy density function for state 1, defined by S1(ρa, ρm, T) = ρasa(ρa, T) + ρms(1)m (ρm, T)

  • S2(ρa, ρm, T)  entropy density function for state 2, defined by S2(ρa, ρm, T) = ρasa(ρa, T) + ρms(2)m (ρm, T)

Constants and defined functions of temperature
  • f = 5.0 × 10−5 s−1  Coriolis parameter

  • g = 9.80665 m s2  acceleration of gravity

  • Ra = 287.05 J kg−1 K−1  gas constant of dry air

  • Rυ = 461.51 J kg−1 K−1  gas constant of water vapor

  • cpa = 1004.675 J kg−1 K−1  specific heat of dry air at constant pressure

  • c = 1850.0 J kg−1 K−1  specific heat of water vapor at constant pressure

  • cυa = cpaRa  specific heat of dry air at constant volume

  • cυυ = cRυ  specific heat of water vapor at constant volume

  • κ = Ra/cpa

  • p0 = 100 kPa  Reference pressure

  • T0 = 273.15 K  Reference temperature

  • ρa0 = p0/(RaT0)  reference density for dry air

  • ρ*υ0 = ρ*υ(T0)  mass density of saturated vapor at T0

  • ρr0 = 1.0 × 10−3kg m−3  reference density for precipitation

  • W0 = 5.5206 m s−1  reference fall velocity

  • τauto = 1000 s  Autoconversion time scale

  • τcol = 455 s  Collection time scale

  • τevap = 763 s  Evaporation time scale

  • E(T)  saturation vapor pressure; E(T) is synthesized from the saturation vapor pressures over water and ice

  • ρ*υ(T) = E(T)/(RυT)  mass density of saturated vapor

  • L(T) = RυT2[d ln E(T)/dT]  specific latent heat for vaporizing condensate at T

  • C(T)  entropy of a unit mass of condensate at T, given by C(T) = cυυ ln(T/T0) − Rυ ln[ρ*υ(T)/ρ*υ0] + L(T0)/T0L(T)/T

  • D(T) = L(T)ρ*υ(T)/T  gain of entropy per unit volume by evaporating a sufficient amount of water, ρ*υ(T), to saturate the volume at T

Others
  • Fm, Fs, Fu, Fυ  boundary layer turbulent fluxes of water vapor, entropy, radial, and azimuthal momentum

  • Qr  conversion rate of ρm to ρr

  • (ξ, ζ) = (−∂m/rz, ∂m/rr)  radial and vertical components of vorticity

  • D/DT  operator ∂/∂t + u∂/∂r + w∂/∂z

  • D/Dt  operator ∂/∂t + u∂/∂r + w∂/∂z

  • m = + ½fr2  absolute angular momentum per unit mass

  • = Dm/Dt  source term for absolute angular momentum

  • = Dm/Dt  source term for absolute angular momentum: = + ρ−1(ρrW + Fm)(∂m/∂z)

  • θ̇ρ = ρ/Dt  source term for virtual potential temperature

  • ρ = Dθρ/Dt  source term for virtual potential temperature: ρ = θ̇ρ + ρ−1(ρrW + Fm)(∂θρ/∂z)

  • P  potential vorticity, as defined as (3.5), or by P = ρ−1∂(m, θρ)/r∂(r, z)

APPENDIX B

Potential Vorticity Equation

The only two prognostic equations involved in the PV derivation are the continuity equation for total density and the tangential momentum Eq. (2.6). The continuity equation for total density is obtained by adding (1 + μ) times (2.1) to ρa times (2.2), which results in

 
formula

where w = w + ρ−1(ρrW + Fm) and D/Dt = ∂/∂t + u∂/∂r + w∂/∂z. The tangential momentum Eq. (2.6) can be written more compactly as

 
formula

where = −ρ−1∂(ρarFυ)/∂z. By differentiation of (B.2) we obtain

 
formula
 
formula

where (ξ, ζ) = (−∂m/rz, ∂m/rr) are the radial and vertical components of vorticity. Next, we define the virtual temperature by Tρ = p/(ρRa) and the virtual potential temperature by θρ = Tρ(p0/p)κ. Taking ∂/∂r and ∂/∂z of Dθρ/Dt = ρ we obtain

 
formula
 
formula

Forming the sum (∂θρ/∂r) (B.3) + ξ(B.5) + (∂θρ/∂z) (B.4) + ζ(B.6), we obtain

 
formula

Using the continuity Eq. (B.1) we can eliminate the divergence from (B.7) to obtain the potential vorticity equation

 
formula

where P is defined in (3.5). Although (B.8) has a compact form, it is desirable to transfer the irreversible processes embedded in the operator D/Dt to the right-hand side of (B.8). Thus, using the relations between D/Dt and D/Dt, ρ and θ̇ρ, and and , we can transform (B.8) to the equivalent form (3.1). It should be noted that the difference between the density-weighted-mean vertical velocity w and dry air vertical velocity w tends to be small. The two are identical in nonprecipitating regions above the boundary layer (where ρr = 0 and Fm = 0), and, even in heavily precipitating regions with rainfall rates of 36 mm h−1, the magnitude of (ρr/ρ)W is only 0.01 m s−1.

APPENDIX C

Axisymmetric Shallow Water Equations

To examine possible problems with the CST numerics under the axisymmetric assumption, we have performed some simple, idealized experiments with the shallow water equations. In the nonrotating, axisymmetric case, the equations for the radial velocity u and the fluid depth h are

 
formula
 
formula

The intent of the numerical experiments is to simulate Rayleigh’s recoil column or the Worthington jet, that is, the column of water that erupts when a drop of water (or milk) impacts the surface after falling from above (Manzello and Yang 2002). The details of the entire process are intricate and involve surface tension effects, which are not included in (C.1) and (C.2). In fact, Rayleigh’s interest concerned the breakup of the recoil column into droplets due to such capillary effects. Our interest is limited to part of the entire process. We simply wish to simulate the column that forms at the center when water rushes into the initially assumed well that is supposed to have been excavated by the impact of the drop in Rayleigh’s experiment. The initial condition is u = 0 and h = h0h1sech[(r/r1)2]. While actual laboratory experiments have spatial scales of centimeters and time scales of milliseconds, our numerical experiment uses the scaled-up dimensions h0 = 1000 m, h1 = 900 m, g = 10 m s−2, r1 = 8 km, Δr = 1 km, Δt = 1 s.

In the top panel of Fig. C1, the thick solid line (labeled Axi.Nonlin) shows the time history of the surface height h at r = 0 for the axisymmetric, nonlinear case, that is, the case when the flow is governed by (C.1) and (C.2). The thin solid line (labeled Axi.Linear) shows the corresponding time history for the axisymmetric, linear case, that is, the case when the flow is governed by the linearized versions of (C.1) and (C.2). Corresponding results for the cases in which axisymmetry is replaced by slab symmetry are shown by the curves labeled Slab.Nonlin and Slab.Linear. The bottom panel of Fig. 11 shows the radial profiles of h at the time of maximum h at r = 0. A prominent recoil column occurs only in the axisymmetric, nonlinear case (with the derivative constraint filter parameter, described in O2, chosen as lc = 2). Although there is no “truth solution” with which to compare Fig. 11, the recoil column produced in the axisymmetric nonlinear case is qualitatively similar to such columns produced in laboratory experiments (e.g., Rein 1996).

We have run similar experiments with different resolutions, with different values of lc, and with different initial conditions. An examination of all these results does not reveal anything peculiar or unexpected in the work of the derivative constraint filter under the axisymmetric assumption. We have also run similar experiments with the Coriolis force and with an azimuthal component υ, but the only effect of rotation is a progressive reduction of the recoil height with increasing f. From these experiments we conclude that there is no fundamental problem with the CST numerics applied to axisymmetric flows.

Footnotes

Corresponding author address: Wayne H. Schubert, Dept. of Atmospheric Science, Colorado State University, Fort Collins, CO 80523-1371. Email: waynes@atmos.colostate.edu

1

A list of symbols is provided in appendix A.

2

Although the assumed boundary layer depth of 1000 m is consistent with many observed tropical flows, the inner core hurricane boundary layer may be shallower than 1000 m, which requires high vertical resolution in the lower troposphere. For a detailed study of the effect of lower tropospheric vertical resolution on hurricane intensity, see Zhang and Wang (2003).

3

The vertical confinement of the boundary layer radial inflow is even more pronounced in models with high vertical resolution in the lowest 2 km (Zhang et al. 2000; Braun 2002; Rogers et al. 2003).

4

More extensive discussions of these and other numerical experiments can be found in Hausman (2001).