## 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*+ ∂(

*rυ*)/

*r*∂

*r*] (∂

*θ*/∂

_{ρ}*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*+ ∂(

*rυ*)/

*r*∂

*r*] (∂

*θ*/∂

_{ρ}*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} m^{2} 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.

## 2. Model

Consider atmospheric matter to consist of dry air, airborne moisture, and precipitation, with respective densities *ρ _{a}*,

*ρ*,

_{m}*ρ*, so that the total density is given by

_{r}*ρ*=

*ρ*+

_{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}*ρ*+

_{υ}*ρ*. However, the partition of

_{c}*ρ*into

_{m}*ρ*and

_{υ}*ρ*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 ∂

_{c}*ρ*/∂

_{a}*t*+ ∂(

*ρ*)/

_{a}ru*r*∂

*r*+ ∂(

*ρ*)/∂

_{a}w*z*= 0, the advective form of which is (2.1). Similarly, the conservation law for

*ρ*is ∂

_{r}*ρ*/∂

_{r}*t*+ ∂(

*ρ*)/

_{r}ru*r*∂

*r*+ ∂[

*ρ*(

_{r}*w + W*)]/∂

*z*=

*Q*, where

_{r}*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

*Q*is the rate of conversion from airborne moisture to precipitation. Defining the precipitation mixing ratio as

_{r}*μ*=

_{r}*ρ*/

_{r}*ρ*, the conservation law for

_{a}*ρ*can be combined with (2.1) and written in the advective form (2.3). Finally, the conservation law for the airborne moisture is ∂

_{r}*ρ*/∂

_{m}*t*+ ∂(

*ρ*)/

_{m}ru*r*∂

*r*+ ∂(

*ρ*+

_{m}w*F*)/∂

_{m}*z*= −

*Q*, where

_{r}*F*is the boundary layer turbulent flux of water vapor. Adding this conservation law for

_{m}*ρ*to the conservation law for

_{m}*ρ*and converting the result into an advective form for the total water mixing ratio

_{r}*μ*= (

*ρ*+

_{m}*ρ*)/

_{r}*ρ*, we obtain (2.2). Note that, although there are four types of matter (with densities

_{a}*ρ*,

_{a}*ρ*,

_{υ}*ρ*,

_{c}*ρ*), 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

_{r}*ρ*into the vapor density

_{m}*ρ*and the airborne condensate density

_{υ}*ρ*will be accomplished diagnostically in the two alternatives of (2.14).

_{c}The total entropy density is *σ* = *σ _{a}* +

*σ*+

_{m}*σ*, consisting of the sum of the entropy densities of dry air, airborne moisture, and precipitation. Since the vertical entropy flux is given by

_{r}*σ*+

_{a}w*σ*+

_{m}w*σ*(

_{r}*w*+

*W*) +

*ρ*=

_{a}F_{s}*σw*+

*σ*+

_{r}W*ρ*, we can write the flux form of the entropy conservation principle as ∂

_{a}F_{s}*σ*/∂

*t*+ ∂(

*σru*)/

*r*∂

*r*+ ∂(

*σw*+

*σ*+

_{r}W*ρ*)/∂

_{a}F_{s}*z*= 0, where

*F*denotes the turbulent eddy flux in the atmospheric boundary layer, and where radiative effects have been neglected. Defining

_{s}*s*=

*σ*/

*ρ*as the “dry-air-specific” entropy of moist air, we can combine (2.1) with the above entropy conservation principle to obtain (2.4).

_{a}We next consider the momentum equations. Defining the absolute angular momentum per unit mass as *m* = *rυ* + ½*fr* ^{2}, we can write the absolute angular momentum budget as ∂(*ρm*)/∂*t* + ∂(*ρmru*)/*r*∂*r* + ∂[(*ρw* + *ρ _{r}W* +

*F*)

_{m}*m*]/∂

*z*= −∂(

*ρ*)/∂

_{a}rF_{υ}*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*=

*p*+

_{a}*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*+

*ρ*, and then approximate this equation by neglecting the material derivative of

_{r}W*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

*ρ*and the inverse square root of

_{r}*ρ*. The largest errors due to this assumption are expected in a small region near the melting level, where the ice factor

_{a}*f*

_{ice}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

*μ*, the specific entropy

_{r}*s*, and the velocity components

*u*,

*υ*,

*w*are

Associated with these seven prognostic equations is a set of diagnostic equations that is required to determine *ρ*, *σ _{r}*,

*p*,

*W*,

*Q*,

_{r}*F*,

_{m}*F*,

_{s}*F*,

_{u}*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

*Q*; and air–sea interaction diagnosis to determine

_{r}*F*,

_{m}*F*,

_{s}*F*,

_{u}*F*. The equations for the thermodynamic diagnosis are

_{υ}where the known functions *S*_{1}, *S*_{2}, *E*, *C*, *ρ**_{υ} are defined in appendix A. The thermodynamic diagnosis can be interpreted conceptually as follows: Input {*ρ _{a}*,

*μ*,

*μ*,

_{r}*s*} ⇒ Output{

*ρ*,

*ρ*,

_{r}*ρ*,

_{m}*ρ*,

_{υ}*ρ*,

_{c}*T*

_{1},

*T*

_{2},

*T*,

*σ*,

_{r}*p*,

_{a}*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}*ρ*from (2.8), the thermodynamically possible temperatures

_{m}*T*

_{1},

*T*

_{2}, and the entropy density of precipitation

*σ*from (2.9) to (2.11), the actual temperature

_{r}*T*of the gaseous matter from (2.12), the partial pressure of dry air from (2.13), the water vapor density

*ρ*, cloud condensate density

_{υ}*ρ*, and water vapor partial pressure

_{c}*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

*E*(

_{w}*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

*E*(

_{w}*T*) and

*E*(

_{i}*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

*T*and the width of the freezing zone Δ

_{f}*T*. For our experiments with ice, we use

_{f}*T*= 273.15 K and Δ

_{f}*T*= 1.0 K, with

_{f}*E*(

_{w}*T*) and

*E*(

_{i}*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 *Q _{r}*. 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

where the constants *ρ _{r}*

_{0},

*W*

_{0},

*τ*

_{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

*Q*is the sum of three terms: the autoconversion term

_{r}*Q*

_{auto}, which initiates the precipitation process when the cloud condensate density

*ρ*exceeds 0.001

_{c}*ρ*; the collection term

_{a}*Q*

_{col}, which is a function of the cloud condensate and precipitation densities

*ρ*and

_{c}*ρ*; and the evaporation term

_{r}*Q*

_{evap}, which is based on the analysis of the evaporation of a stationary drop and then augmented by the dimensionless ventilation factor

*f*

_{vent}(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 *F _{u}* and

*F*appearing in (2.5) and (2.6), the surface values are

_{υ}where the subscript zero denotes a surface value and *V*_{0} = (*u*^{2}_{0} + *υ*^{2}_{0})^{1/2} is the surface wind speed. Similarly, the surface values of the water vapor and entropy fluxes are

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

where *k* = 4.0 × 10^{5} 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 *C _{D}* and

*C*. For example, if, in place of (2.25), one chooses

_{H}*C*= 0.0011 +

_{D}*k*

_{D}V_{0}and

*C*= 0.0011 +

_{H}*k*

_{H}V_{0}, very intense storms are obtained when

*k*= 0 and

_{D}*k*=

_{H}*k*, while much weaker storms are obtained when

*k*=

_{D}*k*and

*k*= 0. Although such variations of

_{H}*k*and

_{D}*k*are extreme, it is an uncomfortable fact that the behavior of

_{H}*C*and

_{D}*C*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.

_{H}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(

*μ*), and then, for the purpose of thermodynamic diagnosis only, inverts these via

_{r}*μ*= ahyp(

*ν*) and

*μ*= ahyp(

_{r}*ν*), where bhyp and ahyp denote the biased hyperbolic transform and its quasi inverse, defined by

_{r}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

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.

## 3. Potential vorticity equation and invertibility principle

Under the assumption of axisymmetry, the radial and vertical components of the absolute vorticity vector are (*ξ*, *ζ*) = (−∂*m*/*r*∂*z*, ∂*m*/*r*∂*r*) = [−∂*υ*/∂*z*, *f* + ∂(*rυ*)/*r*∂*r*], where *m* = *rυ* + ½*fr* ^{2} 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

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

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*/

*p*

_{0})

*r*

^{κ}o*θ*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}*ρ*,

*p*↔

_{a}*p*, 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*/

*r*)/[1 + (

_{m}*r*/

*r*)

_{m}^{2}], with the maximum wind specified as

*υ*= 12 m s

_{m}^{−1}and the radius of maximum wind as

*r*= 100 km. The initial azimuthal velocity is assumed to decrease linearly with height to zero at

_{m}*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.

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.

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 confined^{3} 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.

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*+ ∂(

*rυ*)/

*r*∂

*r*](

*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.

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

*μ*are both large, producing a precipitation rate of over 200 mm h

_{r}^{−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*′ =

*P*−

*P*, 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.

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 ½*fR*^{2} = *m* = *rυ* + ½*fr* ^{2}, 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

and the material derivative operator as

Note that ∂/∂*R* implies fixed *θ _{ρ}*, ∂/∂

*θ*implies fixed

_{ρ}*R*, and ∂/∂𝒯 implies fixed

*R*,

*θ*. Using (4.1)–(4.3) in (3.1), we obtain

_{ρ}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

_{ρ}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*,

*θ*) and

_{ρB}*θ̇*(

_{ρ}*R*,

*θ*), the values of

_{ρB}*P*and

*θ̇*at the top of the boundary layer. Then, integration of (4.5) yields

_{ρ}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.

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

*μ*, 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,

_{r}*Q*

_{auto},

*Q*

_{col},

*Q*

_{evap}, 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 *f*_{ice} to unity. This yields results qualitatively similar to the NICE experiment. In both of these experiments with *f*_{ice} = 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, *ρ _{r}W*, 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

*θ*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.

_{e}## 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 *C _{D}* and

*C*at high wind speeds. Finally, a major limitation of the model (2.1)–(2.25) is the assumption of axisymmetry.

_{H}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.

## 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

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### List of Symbols

##### Mass densities, mixing ratios, temperatures, pressures, velocities

*ρ*mass density of dry air_{a}*ρ*mass density of water vapor_{υ}*ρ*mass density of airborne condensate_{c}*ρ*mass density of precipitating water substance_{r}*ρ*=_{m}*ρ*+_{υ}*ρ*mass density of airborne moisture_{c}*ρ*=*ρ*+_{a}*ρ*+_{m}*ρ*total mass density_{r}*μ*=_{m}*ρ*/_{m}*ρ*mixing ratio of airborne moisture_{a}*μ*=_{r}*ρ*/_{r}*ρ*mixing ratio of precipitating water substance_{a}*μ*=*μ*+_{m}*μ*mixing ratio of total water substance_{r}*T*_{1}temperature for thermodynamic state 1*T*_{2}temperature for thermodynamic state 2*T*= max(*T*_{1},*T*_{2}) temperature*T*=_{ρ}*p*/(*ρR*) virtual temperature_{a}*θ*=_{ρ}*T*(_{ρ}*p*_{0}/*p*)virtual potential temperature^{κ}*θ*=*T*(*p*_{0}/*p*)_{a}dry potential temperature^{κ}*θ*=_{e}*T*_{0}exp[*σ*/(*ρ*)] equivalent potential temperature_{a}c_{pa}*p*partial pressure of dry air_{a}*p*partial pressure of water vapor_{υ}*p*=*p*+_{a}*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}(*ρ*+_{r}W*F*) density-weighted-mean vertical velocity_{m}

##### Specific entropies (J kg^{−1} K^{−1}) and entropy densities (J m^{−3} K^{−1})

*s*(_{a}*ρ*,_{a}*T*) specific entropy of dry air, defined by*s*(_{a}*ρ*,_{a}*T*) =*c*ln(_{υa}*T*/*T*_{0}) −*R*ln(_{a}*ρ*/_{a}*ρ*_{a}_{0})*s*^{(1)}_{m}(*ρ*,_{m}*T*) specific entropy of airborne moisture in state 1:*s*^{(1)}_{m}(*ρ*,_{m}*T*) =*c*ln(_{υυ}*T*/*T*_{0}) −*R*ln(_{υ}*ρ*/_{m}*ρ**_{υ0}) +*L*(*T*_{0})/*T*_{0}*s*^{(2)}_{m}(*ρ*,_{m}*T*) specific entropy of airborne moisture in state 2:*s*^{(2)}_{m}(*ρ*,_{m}*T*) =*C*(*T*) +*D*(*T*)/*ρ*_{m}*s*=_{r}*C*(*T*_{2}) specific entropy of condensed water*s*=*σ*/*ρ*dry-air-specific entropy of moist air_{a}*σ*=_{a}*ρ*entropy density of dry air_{a}s_{a}*σ*=_{m}*ρ*entropy density of airborne water substance_{m}s_{m}*σ*=_{r}*ρ*entropy density of precipitating water substance_{r}s_{r}*σ*=*σ*+_{a}*σ*+_{m}*σ*total entropy density_{r}*S*_{1}(*ρ*,_{a}*ρ*,_{m}*T*) entropy density function for state 1, defined by*S*_{1}(*ρ*,_{a}*ρ*,_{m}*T*) =*ρ*(_{a}s_{a}*ρ*,_{a}*T*) +*ρ*_{m}s^{(1)}_{m}(*ρ*,_{m}*T*)*S*_{2}(*ρ*,_{a}*ρ*,_{m}*T*) entropy density function for state 2, defined by*S*_{2}(*ρ*,_{a}*ρ*,_{m}*T*) =*ρ*(_{a}s_{a}*ρ*,_{a}*T*) +*ρ*_{m}s^{(2)}_{m}(*ρ*,_{m}*T*)

##### Constants and defined functions of temperature

*f*= 5.0 × 10^{−5}s^{−1}Coriolis parameter*g*= 9.80665 m s^{2}acceleration of gravity*R*= 287.05 J kg_{a}^{−1 }K^{−1}gas constant of dry air*R*= 461.51 J kg_{υ}^{−1 }K^{−1}gas constant of water vapor*c*= 1004.675 J kg_{pa}^{−1}K^{−1}specific heat of dry air at constant pressure*c*= 1850.0 J kg_{pυ}^{−1}K^{−1}specific heat of water vapor at constant pressure*c*=_{υa}*c*−_{pa}*R*specific heat of dry air at constant volume_{a}*c*=_{υυ}*c*−_{pυ}*R*specific heat of water vapor at constant volume_{υ}*κ*=*R*/_{a}*c*_{pa}*p*_{0}= 100 kPa Reference pressure*T*_{0}= 273.15 K Reference temperature*ρ*_{a}_{0}=*p*_{0}/(*R*_{a}T_{0}) reference density for dry air*ρ**_{υ0}=*ρ**_{υ}(*T*_{0}) mass density of saturated vapor at*T*_{0}*ρ*_{r}_{0}= 1.0 × 10^{−3}kg m^{−3}reference density for precipitation*W*_{0}= 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*) mass density of saturated vapor_{υ}T*L*(*T*) =*R*_{υ}T^{2}[*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*/*T*_{0}) −*R*ln[_{υ}*ρ**_{υ}(*T*)/*ρ**_{υ0}] +*L*(*T*_{0})/*T*_{0}−*L*(*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

*F*,_{m}*F*,_{s}*F*,_{u}*F*boundary layer turbulent fluxes of water vapor, entropy, radial, and azimuthal momentum_{υ}*Q*conversion rate of_{r}*ρ*to_{m}*ρ*_{r}(

*ξ*,*ζ*) = (−∂*m*/*r*∂*z*, ∂*m*/*r*∂*r*) radial and vertical components of vorticity*D*/*DT*operator ∂/∂*t*+*u*∂/∂*r*+*w*∂/∂*z**D*/*Dt*operator ∂/∂*t*+*u*∂/∂*r*+*w*∂/∂*z**m*=*rυ*+ ½*fr*^{2}absolute angular momentum per unit mass*ṁ*=*Dm*/*Dt*source term for absolute angular momentum=

*Dm*/*Dt*source term for absolute angular momentum: =*ṁ*+*ρ*^{−1}(*ρ*+_{r}W*F*)(∂_{m}*m*/∂*z*)*θ̇*=_{ρ}*Dθ*/_{ρ}*Dt*source term for virtual potential temperature_{ρ}=*Dθ*/_{ρ}*Dt*source term for virtual potential temperature:_{ρ}=*θ̇*+_{ρ}*ρ*^{−1}(*ρ*+_{r}W*F*)(∂_{m}*θ*/∂_{ρ}*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

where *w* = *w* + *ρ*^{−1}(*ρ _{r}W* +

*F*) and

_{m}*D*/

*Dt*= ∂/∂

*t*+

*u*∂/∂

*r*+

*w*∂/∂

*z*. The tangential momentum Eq. (2.6) can be written more compactly as

where = −*ρ*^{−1}∂(*ρ _{a}rF_{υ}*)/∂

*z*. By differentiation of (B.2) we obtain

where (*ξ*, *ζ*) = (−∂*m*/*r*∂*z*, ∂*m*/*r*∂*r*) are the radial and vertical components of vorticity. Next, we define the virtual temperature by *T _{ρ}* =

*p*/(

*ρR*) and the virtual potential temperature by

_{a}*θ*=

_{ρ}*T*(

_{ρ}*p*

_{0}/

*p*)

*. Taking ∂/∂*

^{κ}*and ∂/∂*

_{r}*of*

_{z}*Dθ*/

_{ρ}*Dt*=

_{ρ}we obtain

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

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

*ρ*= 0 and

_{r}*F*= 0), and, even in heavily precipitating regions with rainfall rates of 36 mm h

_{m}^{−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

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* = *h*_{0} − *h*_{1}sech[(*r*/*r*_{1})^{2}]. While actual laboratory experiments have spatial scales of centimeters and time scales of milliseconds, our numerical experiment uses the scaled-up dimensions *h*_{0} = 1000 m, *h*_{1} = 900 m, *g* = 10 m s^{−2}, *r*_{1} = 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 *l _{c}* = 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 *l _{c}*, 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).