A Simple Family of Tropical Cyclone Models

,


Introduction
An extensive hierarchy of numerical models allows the tropical cyclone research community to simulate the intensification and movement of tropical cyclones. At the top of this hierarchy, non-hydrostatic, three-dimensional, full-physics models are the ones of most interest for accurate tropical cyclone prediction. At the bottom of this hierarchy, simple, axisymmetric models remain useful for providing a conceptual understanding of certain aspects of the tropical cyclone problem. Here, the purpose is to discuss the equations that comprise a simple family of dynamical models capable of simulating tropical cyclone life cycles, including intensification, the formation of boundary layer shocks, and the development of an eyewall. In particular, we will consider the following four models: All four models are axisymmetric and hydrostatic, with only three layers and with moisture predicted only in the lowest layer. All have the same simple parameterizations of convective mass fluxes and air-sea interaction. The most general model considered here is the primitive equation model (PE), which retains the primitive form of the momentum equations in all three layers. It is the only one of the four models that does not filter propagating inertia-gravity waves. The first hybrid balanced model (HB1) assumes gradient balance in the upper two layers but retains the primitive form of the momentum equations in the boundary layer. The second hybrid balanced model (HB2) is slightly simpler than HB1 in that diagnostic forms of the boundary layer momentum equations are used, so that the boundary layer pumping is diagnostic but more accurate than in GB. HB2 is identical to the model briefly studied by Ooyama [1]. Finally, the gradient balanced model (GB) assumes gradient balance in all three layers, with the boundary layer pumping becoming both diagnostic and local. GB is identical to the model extensively studied by Ooyama [2].
This review is organized in the following way. Section 2 and Table 1 give a summary of the three-layer primitive equation model. Then, Sections 3 and 4, along with Tables 2 and 3, provide corresponding summaries of the hybrid balanced models HB1 and HB2, while Section 5 and Table 4 give a summary of the fully balanced model GB. Section 6 discusses how the prognostic boundary layer dynamics of the PE and HB1 models can be isolated from the dynamics of the overlying layers and then studied as a one-way interaction problem. The potential vorticity aspects of the four models are then discussed in Section 7.

Primitive Equation Model (PE)
As illustrated in Figure 1, all four models consider axisymmetric motions in three layers of an incompressible fluid on an f -plane with subscripts denoting the layers from 0 for the lowest, 1 for the middle, and 2 for the top layer. The constant density in the lowest two layers is ρ 0 = ρ 1 = ρ, and that in the upper layer is ρ 2 = ρ, where the constant satisfies < 1. The lowest layer has constant thickness h 0 , while the upper two layers have variable thicknesses h 1 (r, t) and h 2 (r, t). The radial velocity components in the three layers are u 0 (r, t), u 1 (r, t), and u 2 (r, t), and the corresponding radial mass fluxes are ψ 0 (r, t) = −ru 0 h 0 , ψ 1 (r, t) = −ru 1 h 1 , and ψ 2 (r, t) = −ru 2 h 2 . The tangential velocity components in the three layers are v 0 (r, t), v 1 (r, t), and v 2 (r, t).
The hydrostatic pressure at height z, which depends on the weight of the overlying fluid, is given by where the first line applies in layers 0 and 1 (the lower two layers), while the second line applies in layer 2 (the upper layer). Assume that, in the far-field, the upper two layer thicknesses approach the constant valuesh 1 andh 2 , so that the corresponding hydrostatic pressures in the far-field arep 1 (z) = gρ(h 0 +h 1 − z) + g ρh 2 , p 2 (z) = g ρ(h 0 +h 1 +h 2 − z). ( Defining the geopotential anomalies as φ 1 = (p 1 −p 1 )/ρ and φ 2 = (p 2 −p 2 )/ ρ, it is easily shown from (1) and (2) that so that the radial pressure gradient force in the lower two layers is (1/ρ)(∂p 1 /∂r) = (∂φ 1 /∂r), while that in the upper layer is (1/ ρ)(∂p 2 /∂r) = (∂φ 2 /∂r). Figure 1. Schematic diagram of the axisymmetric, three-layer, f -plane family of models. The constant density in the lower two layers is ρ 0 = ρ 1 = ρ, while the constant density in the upper layer is ρ 2 = ρ, with < 1. The thicknesses of the layers are h 0 , h 1 , and h 2 , with h 0 a specified constant, and with h 1 and h 2 functions of (r, t). The radial velocity components are u 0 , u 1 , and u 2 , and the radial mass fluxes are ψ 0 = −ru 0 h 0 , ψ 1 = −ru 1 h 1 , and ψ 2 = −ru 2 h 2 . The tangential velocity components are v 0 , v 1 , and v 2 . The surface stress τ s drives the boundary layer radial inflow, which results in boundary layer pumping (w > 0) in the inner region and boundary layer suction (w < 0) in the outer region. Diabatic processes are parameterized through the mass transport terms Q + and Q − . Adapted from ([2], Figure 1).

The mass continuity equations are
where w is the vertical velocity at the top of the boundary layer and Q = Q + − Q − is the diabatic mass flux between the upper two layers with Q + being the upward diabatic mass flux (heating) from layer 1 to layer 2 and Q − being the downward diabatic mass flux (cooling) from layer 2 to layer 1. Adding the three continuity Equations (4)-(6), noting the cancellation of the terms on the right-hand side, and then integrating over the area and specifying the outer boundary condition to be we conclude that i.e., the total mass is conserved.
The radial momentum equations for the three layers are given by Equations (1.11)-(1.13) of Table 1, where f is the constant Coriolis parameter, the right-hand side terms are given by w ± = 1 2 (|w| ± w), so w = w + − w − , µ is the constant coefficient for the frictional stress between the layers, c D is the drag coefficient, 1/2 is the boundary layer wind speed, and the three Λ (u) j terms are horizontal diffusive fluxes for the u-equations as given by (A2) in Appendix A. Similarly, the tangential momentum equations for the three layers are given by Equations (1.14)-(1.16) of Table 1, where the right-hand side terms are given by where the three Λ (m) j terms are horizontal diffusive fluxes for the m-equations as given by (A4) in Appendix A. Since the momentum Equations (1.11)-(1.16) may have a somewhat unfamiliar form, a derivation from first principles is given in Appendix A.
Rather than predict h 1 (r, t) and h 2 (r, t) via (5) and (6), it is more convenient to predict φ 1 (r, t) and φ 2 (r, t). The prognostic equations for φ 1 (r, t) and φ 2 (r, t) are easily obtained from (3), (5), and (6), with the result given in Equations (1.17) and (1.18) of Table 1. After φ 1 (r, t) and φ 2 (r, t) are predicted from (1.17) and (1.18), the layer depths can be recovered from (1.1) and (1.2), which follow directly from (3). Note that σ = 1 − is a measure of static stability. Now, consider the parameterization of the convective mass flux Q + . The equivalent potential temperature budget in the convection is ηθ ec = θ e0 + (η − 1)θ e1 , where η is the entrainment parameter, θ ec is the equivalent potential temperature in the convection in the upper layer, θ e0 is the equivalent potential temperature in the boundary layer, and θ e1 is the equivalent potential temperature in layer 1. We assume that the convection has neutral buoyancy in the upper layer, so that θ ec = θ * e2 , where θ * e2 is the saturation equivalent potential temperature in the upper layer. Combining these last two relations, we obtain Equation (1.9) of Table 1. Note that conditional instability occurs in a column when θ e0 > θ * e2 , which means η > 1 since θ * e2 > θ e1 . The equivalent potential temperature of the middle layer is assumed to be a specified parameter, while the boundary layer equivalent potential temperature, θ e0 (a prognostic variable), and the upper layer saturation equivalent potential temperature, θ * e2 (a diagnostic variable), depend on radius and time. Finally, Q + is given by Q + = ηw + , so it then follows that Q = ηw + − Q − as indicated by Equation (1.10) of Table 1, where Q − is a specified parameter (typically zero).  (9) and (10). Because the radial components u 1 (r, t) and u 2 (r, t) are predicted, the PE model allows propagating inertia-gravity waves. Because the boundary layer radial and tangential components u 0 (r, t) and v 0 (r, t) are predicted, and because large inflow often occurs in the boundary layer, the PE model can develop boundary layer shocks. The diabatic mass flux Q + depends on boundary layer convergence (w > 0) and the entrainment parameter η, which in turn depends on the boundary layer equivalent potential temperature θ e0 , predicted by (1.19), and the upper tropospheric saturation equivalent potential temperature θ * e2 , diagnosed from (1.4).

PROGNOSTIC EQUATIONS:
Radial wind components: Tangential wind components: Geopotential anomalies: Boundary layer equivalent potential temperature: To include the effects of the warm core on convection, we need to interpret θ * e2 in terms of potential temperature, even though potential temperature does not appear explicitly in the three-layer model. This can be accomplished by writing the hydrostatic equation for a compressible fluid in the differential form d(φ −φ) = −c p (θ −θ) d(p/p 00 ) κ , where φ −φ is the deviation of the geopotential φ from its far-field valueφ, θ −θ is the deviation of the potential temperature θ from its far-field valueθ, p 00 = 1000 hPa is the reference pressure, and κ = R/c p . Based on a discrete form of this hydrostatic equation, it is reasonable to define the middle-tropospheric potential temperature anomaly θ m −θ m , whereθ m is the middle-tropospheric potential temperature in the far field, as where we use 700 and 300 hPa as representative values of the pressure in the upper two layers. Thus, in the context of the three-layer model, φ 2 − φ 1 can be interpreted as proportional to the middle-tropospheric potential temperature anomaly. Using a Taylor series approximation (or simple inspection of the upper tropospheric spacing of θ-lines and θ * e -lines on a thermodynamic diagram), we can justify the approximation θ * e2 =θ * e2 + 2.0(θ m −θ m ), whereθ * e2 is the upper layer saturation equivalent potential temperature in the far field. Combining this approximation with (11) yields Equation (1.4) of Table 1. The spatial and temporal dependence of θ * e2 can serve as an important negative feedback on tropical cyclone intensification, since an evolving warm core can reduce θ e0 − θ * e2 and hence reduce η via Equation (1.9).
The equivalent potential temperature is predicted only in the lowest layer. In flux form, the boundary layer equivalent potential temperature equation is where c E is the exchange coefficient, U = (u 2 0 + v 2 0 ) 1/2 is the boundary layer wind speed, and θ * es is the saturation equivalent potential temperature at the sea-surface temperature and pressure. When (12) is converted to advective form, we obtain Equation (1.19) of Table 1. Important controls on moist convection occur in the terms of (1.19). In the tropical cyclone inner region, where there is boundary layer pumping (w > 0 so that w − = 0), the term w − (θ e0 − θ e1 )/h 0 vanishes. However, in the tropical cyclone outer region, where there is boundary layer suction (w < 0 so that w − > 0), the term w − (θ e0 − θ e1 )/h 0 is positive and tends to produce a decrease of boundary layer equivalent potential temperature through the downward transport of dry middle-tropospheric air (θ e1 < θ e0 ). Air-sea interaction plays a crucial role in maintaining high values of θ e0 through the surface flux term c E U(θ * es − θ e0 )/h 0 . The factor c E U becomes large in regions of high surface winds, while the saturation equivalent potential temperature at the sea surface (θ * es ) becomes large in regions where the surface pressure is low, an effect that is incorporated through Table 1. Both effects enhance the surface flux c E U(θ * es − θ e0 ) and can maintain θ e0 > θ * e2 , and hence η slightly larger than unity, even though an upper tropospheric warm core is being produced and θ * e2 is increasing through Equation (1.4). Since θ * es depends on both the sea surface temperature and pressure, and since the sea surface pressure depends on φ 1 , we should regard θ * es as a dependent variable that increases as the sea surface pressure decreases. Using a Taylor series approximation (or simple inspection of the lower tropospheric spacing of p-lines and θ * e -lines on a thermodynamic diagram) the dependent variable θ * es can be parameterized as given in Equation (1.3) of Table 1. Thus, in the far field, where the geopotential anomaly φ 1 is very small, the value of θ * es isθ * es = 372 K for the typical case examined in Section 8. However, if φ 1 in the inner core is approximately −5000 m 2 s −2 (i.e., a −50 hPa sea surface pressure anomaly), then the value of θ * es is locally increased to 382 K. This dependence of θ * es on φ 1 plays an important role in enhancing the surface flux of θ e in the inner core of the vortex. Finally, we list in Equation (1.6) the assumed c D and c E formulas used by Ooyama for his standard case. The results of some very interesting experiments with other parameterizations for c D and c E are discussed by Ooyama ([2], Section 14).
The nine prognostic variables of the primitive equation model are u 0 (r, t), u 1 (r, t), , and θ e0 (r, t), so all the other variables are diagnostic. The steps of the numerical integration procedure start with knowledge of these nine prognostic variables, either from the initial conditions or from the previous time step. The procedure for advancing one time level is as follows: 1.
Specify desired values for the various fixed parameters: g, c p , f , or where K is the constant kinematic coefficient of eddy viscosity appearing in the equations for the horizontal diffusive fluxes given in (A2) and (A4) in Appendix A.

3.
Compute the boundary layer wind speed U(r, t) from (1.5), and then the exchange coefficient c E (r, t) and the drag coefficient c D (r, t) from (1.6).
Using the definitions of F j (r, t) and G j (r, t) given in (9) and (10), predict new values of the radial and tangential wind components from (1.11)-(1.16). 7.
Using the known boundary layer radial velocity u 0 (r, t), the diagnosed boundary layer suction w − (r, t), and the diagnosed θ * es (r, t), predict the boundary layer equivalent potential temperature θ e0 (r, t) via (1. 19), and then return to step 2 for the beginning of the next time step.
The dynamical core of this model has no explicit thermodynamics. Since there is a discontinuity in the density at the interface between layers 1 and 2, we can regard this interface as containing many constant-density surfaces packed closely together. Through the analogy between density surfaces and potential temperature surfaces (see [3]), we can also regard this interface as containing many closely packed θ surfaces. The mass flux Q + can therefore be regarded as being due to latent heat release in the moist convective updrafts that carry cloudy air across θ surfaces into the upper troposphere. The four models discussed here connect Q + to boundary layer convergence, boundary layer θ e , and upper tropospheric temperature in a very simple way. It should be noted that nonhydrostatic, full-physics models of tropical cyclones explicitly model the details of water and ice microphysics, the different categories of ice, the different fall velocities of liquid and ice, the role of aerosols in the size distributions of condensate, etc. The simplified three-layer models presented here neglect all such microphysical details, but do capture the most fundamental effect of precipitating moist convection, i.e., the convective mass flux across θ surfaces from the lower troposphere to the upper troposphere. The dynamical adjustments to this mass flux are what produce the lower tropospheric cyclone and the upper tropospheric anticyclone.
Equations (1.1)-(1.19) constitute one of the simplest PE models capable of simulating tropical cyclone life cycles, including intensification, the formation of boundary layer shocks, and the development of an eyewall.

Hybrid Balanced Model 1 (HB1)
The equations for the HB1 model are listed in Table 2. The most obvious difference between Tables 1 and 2 is that for HB1 the number of diagnostic equations has increased from ten to fourteen, while the number of prognostic equations has decreased from nine to five. The ten diagnostic equations in Table 1 carry over directly to Table 2, and four additional diagnostic equations appear. In the HB1 model, the prognostic equations for u 1 and u 2 are discarded and replaced by the gradient wind relations (2.11) and (2.12). With gradient balance in the middle and upper layers, the predictions of v 1 and v 2 are redundant with the predictions of φ 1 and φ 2 . In other words, taking the local time derivative of the gradient wind equation yields Using the tangential wind equation to eliminate (∂v j /∂t) and the mass continuity principle to eliminate (∂φ j /∂t), we obtain the diagnostic Equations (2.13) and (2.14) of Table 2, where the local inertial stability S j in layer j is defined by where ζ j = ∂(rv j )/r∂r is the relative vorticity in layer j and P j = ( f + ζ j )h j /h j is the potential vorticity in layer j. Equations (2.13) and (2.14) constitute a coupled pair of second order diagnostic equations (i.e., the Eliassen transverse circulation equations given in reference [4]) for the radial mass fluxes ψ 1 = −ru 1 h 1 and ψ 2 = −ru 2 h 2 . These radial mass fluxes are forced by the diabatic and frictional effects that appear on the right-hand sides of (2.13) and (2.14), while the spatial distribution of these fluxes is shaped by the inertial stability factors on the left-hand sides. Since (2.13) and (2.14) are diagnostic equations, an instantaneous change in the diabatic and frictional forcing produces an immediate response in the entire fluid. This "action at a distance" property does not occur in the PE model since fluctuations of the inner core forcing excite inertia-gravity waves, which require a finite time to spread information laterally and thereby adjust the radial mass fluxes in the surrounding fluid. The forcing terms on the right-hand side of the transverse circulation Equations (2.13) and (2.14) involve G 1 , G 2 , (∂w/∂r), and (∂Q/∂r). For the idealized situation in which the absolute angular momentum is materially conserved in the upper two layers, we have G 1 = G 2 = 0, so the forcing of the transverse circulation then involves only the radial derivatives of the boundary layer pumping w and the diabatic mass flux Q. If the diabatic mass flux Q were to vanish and w were to be largest at the center of the vortex, then (∂w/∂r) < 0 over a core region, which results in u 1 > 0 and u 2 > 0, with the relative magnitudes of u 1 and u 2 determined by the relative magnitudes of the inertial stability factors S 1 and S 2 . Such radial outflow leads to the spin-down (at all radii) of the cyclonic vortex in layers 1 and 2. As discussed by Greenspan and Howard [5] and Greenspan ([6], Sections 2.3 and 2.4), this is what occurs in a laminar Ekman layer with a no-slip lower boundary condition. However, as discussed by Eliassen and Lystad [7], the situation is different in a turbulent Ekman layer with a thin Prandtl surface layer and stress condition at the lower boundary. In that case, the Ekman pumping does not maximize at the vortex axis, but rather at some finite radius. This effect seems crucial in understanding the structure of the eye and eyewall. All four models discussed here can simulate this important effect.  The upper half of the table lists the diagnostic equations, while the lower half lists the five prognostic equations. As indicated in (2.11) and (2.12), gradient balance is assumed in the upper two layers, but, as indicated in (2.15) and (2.16), the primitive forms of the momentum equations are retained in the boundary layer, thus allowing the formation of boundary layer shocks. Because the prognostic equations for u 1 (r, t) and u 2 (r, t) have been discarded and replaced by the coupled Eliassen Equations (2.13) and (2.14), propagating inertia-gravity waves are filtered in the HB1 model. A sufficient (but not necessary) condition for the uniqueness of the solutions ψ 1 , ψ 2 of the Eliassen equations is that S 1 ≥ 0 and S 2 ≥ 0. Numerical integrations of the Eliassen equations reveal that S 2 can become negative at some distance from the cyclone center, but no difficulty is encountered in the solution of the discretized Eliassen equations if a direct method of solution (such as Gaussian elimination) is used.

PROGNOSTIC EQUATIONS:
Boundary layer radial and tangential wind components: Geopotential anomalies: Boundary layer equivalent potential temperature: The five prognostic variables of the HB1 model are u 0 (r, t), v 0 (r, t), φ 1 (r, t), φ 2 (r, t), and θ e0 (r, t), so all the other variables are diagnostic. HB1 is "hybrid" in the sense that gradient balance is used in the upper two layers but not in the boundary layer. The steps of the numerical integration procedure start with knowledge of these five prognostic variables, either from the initial conditions or from the previous time step. The procedure for advancing one time level is very similar to that outlined in the previous section for the PE model.

Hybrid Balanced Model 2 (HB2)
The equations for the HB2 model are listed in Table 3. The most obvious difference between Tables 2 and 3 is that for HB2 the number of prognostic equations has decreased from five to three. The diagnostic equations in Table 2 carry over directly to Table 3, and two additional diagnostic equations appear. In the HB2 model the prognostic equations for u 0 and v 0 have been discarded and replaced by the diagnostic relations (3.15) and (3.16), which are essentially quasi-steady state versions of (2.15) and (2.16). Ooyama [1] argued that the system of boundary layer diagnostic primitive Equations (3.15) and (3.16) can be numerically integrated (radially inward) under appropriate boundary conditions on u 0 and v 0 at a large radius (1000 km, say). He did not impose boundary conditions at r = 0, but found that u 0 and v 0 approach zero as r goes to zero, which is obviously a relief but not an entirely satisfactory result. Frisius and Lee [8] provided a different, and quite successful, approach to this issue. Since the equations contain both a forcing term (the imposed horizontal pressure gradient of the overlying layer) and surface friction terms, the solutions u 0 and v 0 are determined more or less locally. Here, we use the term "semi-local" as synonymous with Ooyama's "more or less local"'. Note that Ooyama's research with the HB2 model [1] was published in the Proceedings of the WMO/IUGG Symposium on NWP, held in Tokyo in 1968. That work is an extension of his modeling work with the GB model on the life cycle of tropical cyclones [2] and shows the importance of the u 0 (∂u 0 /∂r) term in shaping the boundary layer pumping and thus influencing storm size and intensification rate. A comparison between HB1 and GB will be provided in Section 8. Although Ooyama's HB2 model work [1] is less well-known than his GB model work [2], it is a fundamental contribution to our understanding of tropical cyclones. Finally, it might be said that, without due caution, solving for the boundary layer winds u 0 and v 0 in HB2 can be fraught with danger, in the sense that a particular numerical method may become erratic as it tries to produce multivalued solutions [9,10].

Gradient Balanced Model (GB)
The equations for the gradient balanced model are listed in Table 4. The only differences between Tables 3 and 4 are that for the GB model, the gradient balance assumption has also been used for the boundary layer and the boundary layer tangential wind equation has been replaced by the local relation (4.15), which is essentially a high Rossby number version of Ekman balance. In spite of the simplicity of the GB model, it is fairly successful at simulating the life cycle of tropical cyclones. However, as discussed by Williams et al. [11], one deficiency of the GB model is its inability to simulate the formation of boundary layer shock-like structures in the radial wind and hence the very concentrated boundary layer pumping associated with eyewalls. This deficiency is primarily due to the neglect of the u 0 (∂u 0 /∂r) term in the boundary layer dynamics of the GB model.
Ooyama extensively studied the GB model [2] and concluded that this model is capable of simulating many aspects of developing and mature tropical cyclones, including the importance of the oceanic latent heat supply. However, one obvious defect of the GB model is the production of a radius of maximum wind that is too large and that expands too rapidly with time. It was hypothesized that this defect is related to the use of the gradient balance approximation in the boundary layer. Therefore, Ooyama proposed the HB2 model [1], which differs from the GB model only in the sense that the boundary layer wind components u 0 and v 0 are computed using the quasi-steady state form of the boundary layer primitive equations, i.e., using the boundary layer equations listed as (3.15) and (3.16) of Table 3. Table 3. Summary of the HB2 model. The upper half of the table lists the diagnostic equations, while the lower half lists the three prognostic equations. Some results from this model were given by Ooyama [1], who discussed how the results differ from the GB model shown in Table 4. These differences between HB2 and GB are discussed in Section 8, along with the observation that the numerical solution of the nonlocal diagnostic Equations (3.15) and (3.16) can lead to mathematical difficulties.

PROGNOSTIC EQUATIONS:
Geopotential anomalies: Boundary layer equivalent potential temperature:  The upper half of the table lists the diagnostic equations, while the lower half lists the three prognostic equations. Note that the boundary layer radial inflow u 0 (r, t) and the boundary layer pumping w(r, t) are computed from (4.15) and (4.7), which are local diagnostic relations. This differs from the other three models in which the boundary layer dynamics is nonlocal. Although the boundary layer radial flow u 0 (r, t), determined by (4.15), can vary rapidly in r because ζ 1 (r, t) varies rapidly in r, true boundary layer shocks are not produced in the GB model because of the local nature of the simplified boundary layer dynamics. An extensive set of numerical integrations using this model is given by Ooyama [2].

PROGNOSTIC EQUATIONS:
Geopotential anomalies: Boundary layer equivalent potential temperature:

Isolating the Slab Boundary Layer Model from the PE and HB1 Models
To understand the important role of boundary layer dynamics in shaping the spatial distribution of convective mass flux, it is useful to isolate the slab boundary layer equations from the PE model equations of Table 1 and the HB1 model equations of Table 2. This can be accomplished by assuming that µ = 0, u 1 is negligible, and v 1 is the specified gradient wind v gr (r, t), where With these assumptions, the boundary layer radial, tangential, and continuity equations become We can regard (16) as a closed system in u 0 (r, t), v 0 (r, t), w(r, t), and w − (r, t), with the imposed pressure gradient force specified in terms of v gr (r, t). As discussed by Williams et al. [11], in the absence of horizontal diffusion terms, the slab boundary layer Equations (16) form a hyperbolic system, which means that it can be written in characteristic form. The derivation of this characteristic form was given by Slocum et al. [12], whose Equations (A.13) and (A.14) show that in regions where w < 0 there are two families of characteristics, one given by (dr/dt) = 2u and one given by (dr/dt) = u, while in regions where w > 0 there is only one family of characteristics, given by (dr/dt) = u. Since these characteristics can intersect, there is the possibility of shocks developing, resulting in discontinuities in u and v and singularities in w and ζ.
The high-resolution (∆r = 100 m) numerical solutions of (16) presented by Williams et al. [11] and Slocum et al. [12] do contain near-singularities in the boundary layer pumping w(r, t) and the vorticity ζ(r, t). Obviously, true singularities do not occur in nature; their mathematical existence reflects the simplicity of the physics included in (16). In a three-dimensional, non-hydrostatic, full-physics hurricane model, spikes in the radial distribution of boundary layer pumping might be expected to collapse to the spatial scale of an individual cumulonimbus cloud, within which the vertical velocity would be limited by non-hydrostatic moist physics. Such full-physics simulations of Supertyphoon Haiyan (2013) have been performed by Tsujino and Kuo [13]. In their paper, Figures 10 and 12 illustrate the development of a near-singularity in the vertical motion field (w ≈ 20 m s −1 ) at z ≈ 434 m and r ≈ 20 km, if the horizontal resolution of their three-dimensional CReSS model is reduced to ∆x = ∆y = 500 m.

Potential Vorticity Aspects of the Models
Since it can lead to additional insights, it is interesting to interpret the four models in terms of potential vorticity (PV) dynamics (see [14][15][16][17]). Although the four models differ in their treatment of the boundary layer, they have identical continuity equations and angular momentum equations in the two main layers. Thus, all four models have the same PV dynamics in the main layers. Note that because h 0 is a constant, there is no difference between vorticity and potential vorticity in the boundary layer. The PV equations for the two main layers are where D j /Dt = (∂/∂t) + u j (∂/∂r) and P j = ( f + ζ j )(h j /h j ) for j = 1, 2. Since the HB1, HB2, and GB models use gradient wind balance, these three models possess an invertibility principle relating v 1 , v 2 , P 1 , P 2 . This invertibility principle is obtained by taking (∂/∂r) ofh j ( f + ζ j ) = h j P j to obtainh j (∂ζ j /∂r) − P j (∂h j /∂r) = h j (∂P j /∂r). Then, expressing (∂h j /∂r) in terms of (∂φ j /∂r), and making use of the gradient wind relation ( f + v j /r)v j = (∂φ j /∂r), we obtain ∂ ∂r This is a pair of coupled, second order equations for v 1 (r, t) and v 2 (r, t), knowing the potential vorticity fields P 1 (r, t) and P 2 (r, t) predicted from (17). Neglecting the G 1 and G 2 terms in (17) and considering the conditionally unstable vortex core region where η > 1 and w > 0, the PV equations take the form (D 1 P 1 /Dt) = (P 1 /h 1 )(η − 1)w and (D 2 P 2 /Dt) = −(P 2 / h 2 )ηw, so that (D 1 P 1 /Dt) > 0 and (D 2 P 2 /Dt) < 0 there. Note that, if the development of the upper tropospheric warm core reduces η to unity, this serves as a brake on the otherwise exponential growth of P 1 . Very large values of P 1 in the vortex core result not only in large cyclonic flow in layer 1 but also in substantial cyclonic flow in layer 2, a result that can be attributed to the coupling between the two equations in (18).
Ooyama performed a linear stability analysis of the GB model ( [2], Section 8). In this simplified linear analysis, curvature effects are neglected, so axisymmetry is replaced by line-symmetry and gradient balance is replaced by geostrophic balance. The problem is reduced to the coupled equations for φ 1 and φ 2 listed in his entry (8.8). Although he does not interpret these equations as PV equations, they can easily be rearranged into the PV forms ∂ ∂t whereh is the basic state depth of each layer,η is the basic state entrainment parameter, and the constant k s corresponds to c D U. The two quantities in large parentheses in (19) are proportional to the perturbation PV in each layer, with (∂ 2 φ 1 /∂r 2 ) and (∂ 2 φ 2 /∂r 2 ) proportional to the relative vorticity in the two layers. The two terms on the right-hand sides of (19), which are proportional to k s (∂ 2 φ 1 /∂r 2 ), represent the effects of Ekman pumping. When coupled with convective mass fluxes and withη > 1, these terms result in an inner core increase of PV in layer 1 and a decrease in layer 2. During such a linear "CISK" process, both the primary and secondary circulations grow exponentially with time. However, in the nonlinear context, described for layer 1 by (D 1 P 1 /Dt) = (P 1 /h 1 )(η − 1)w, the PV can grow exponentially even when the secondary circulation and (η − 1)w are nearly fixed in time. Since the assumptions of linear dynamics are violated even in tropical depressions (see the caption of Table 5), it follows that it is the P 1 factor in (P 1 /h 1 )(η − 1)w that is most relevant to tropical cyclone rapid intensification, with the factor (η − 1) serving as a convective limiter on the intensification process. Recently, Hendricks et al. [18] have developed a shallow water axisymmetric model for intensification (SWAMI) and have shown how a PV limiter is important for the practical use of their model.

Concluding Remarks
A comparison of the four models is given in Figure 2. The dependent variables of the four models are listed in the left column. The variables h 1 , h 2 , θ * es , θ * e2 , U, c D , c E , w, w ± , η, and Q are diagnostic in all models, while the variables φ 1 , φ 2 , and θ e0 are prognostic in all models. The boundary layer wind components u 0 and v 0 are prognostic in the PE and HB1 models, but are diagnostic in the HB2 and GB models; this is related to the appearance of boundary layer shocks, as indicated by the inner arrow on the right. The radial and tangential wind components u 1 , u 2 , v 1 , and v 2 are prognostic in the PE model but are diagnostic in the HB1, HB2 and GB models; this is related to the occurrence of propagating inertia-gravity waves, as indicated by the outer arrow on the right.  Through a comparison of the GB model Equations (Table 4) with the HB2 model Equations (Table 3), it is apparent that the only difference in the model equations is the more general, nonlocal boundary layer dynamics of the HB2 model. However, this difference leads to significant changes to the tropical cyclone life cycles simulated by the two models. These result from different Ekman pumping distributions. To illustrate these differences, Ooyama made comparison runs [1] using the constants listed in Table 5, the initial condition φ 2 (r, 0) = 0, and with the initial φ 1 (r, 0) in gradient balance with v 1 (r, 0) =v 2(r/r)

Model
where the constantsr andv are also given in Table 5. The initial condition on θ e0 (r, 0) is determined from η(r, 0) = 2. Ooyama found in a numerical test that the values of u 0 and v 0 within 50 km from the vortex center are not influenced very much by their values outside 100 km. Solutions from HB2 are plotted as the solid curves in Figure 3 below (adapted from [1], Figure 2) and show three features found in later modeling studies [9][10][11][12]19] and also in observations of Hurricane Hugo [20]: (1) Subgradient boundary layer tangential flow outside the radius of maximum wind; (2) Supergradient boundary layer tangential flow near the radius of maximum wind; (3) Maximum Ekman pumping at a much smaller radius than in the GB model, resulting in a tighter vortex that intensifies more rapidly. As shown in Figure 3, the GB model produces an Ekman pumping with a 2 m s −1 maximum at r ≈ 60 km with a sharp edge on the outside of the eyewall, while the HB2 model produces an Ekman pumping that maximizes at r ≈ 15 km with a sharp edge on the inside of the eyewall. This difference is due to the u 0 (∂u 0 /∂r) term in the HB2 model, i.e., to the nonlocal character of the boundary layer dynamics in HB2. This term leads to a nearly discontinuous behavior of the boundary layer radial inflow near r ≈ 15 km, and hence the maximum Ekman pumping there. This behavior of u 0 in the HB2 model indicates the danger of formulating the nonlocal boundary layer dynamics without horizontal diffusion terms. If, without horizontal diffusion, the u 0 (∂u 0 /∂r) term leads to the formation of a shock, the derivatives in the boundary layer equations can cease to exist, and the equations are no longer a useful description of the physics near the shock. This situation can be addressed by fitting in a shock condition or using horizontal diffusion to prevent multivalued radial inflow. The latter was the approach originally suggested by von Neumann and Richtmyer [21]. Table 5. Specified constants for the GB and HB2 model runs shown in Figure 3. Note that the large value of the initial average vorticity insider, given by 2v/r = 8 f , makes linear CISK theory (see Section 8 and references [2,22]) of limited use. A particular limitation of the Ooyama (1969) cumulus parameterization scheme is the assumption that, if the column is conditionally unstable, the total cloud base vertical mass flux due to deep convection is equal to the vertical mass flux associated with the larger-scale horizontal convergence in the frictional boundary layer. Later studies showed that this assumption is not consistent with observations and with non-hydrostatic, full-physics models. For example, Yanai et al. [23] performed a heat and moisture budget analysis using an array of five radiosonde stations in the Marshall Islands, which are located in the tropical Pacific ITCZ. From the apparent heat source Q 1 and the apparent moisture sink Q 2 , and using a simple bulk cloud model, Yanai et al. found that for cloud clusters in the ITCZ region, the total cloud base vertical mass flux is much larger than the vertical mass flux associated with the horizontal convergence in the frictional boundary layer. This conclusion about the vertical convective mass flux at the top of the boundary layer in weak tropical disturbances has been confirmed by several other observational studies, for example, [24][25][26]. However, as found by Smith et al. [27] using the full-physics CM1 model [28], this may not be the case for hurricanes and typhoons in their mature and decaying phases. The vortex in HB2 is tighter and has gone through a more rapid intensification than the vortex in GB. This difference in structure and intensification rate is due to the striking differences in Ekman pumping, with the GB model having its 2 m s −1 maximum w at r ≈ 60 km and the HB2 model having a similar maximum but at r ≈ 15 km. This illustrates the fundamental importance of the u 0 (∂u 0 /∂r) term in the radial equation of boundary layer motion in the HB2 model. Note that the boundary layer tangential flow v 0 in the HB2 model is slightly subgradient for r > 30 km and slightly supergradient for r < 30 km. This subgradient/supergradient effect can be more apparent in strong storms such as Hurricane Hugo (see [11,20]). Adapted from ([1], Figure 2).

Parameter
In order to remove some of the limitations of the Ooyama model, Zhu et al. [29] developed a three-layer, σ-coordinate tropical cyclone model with the option of explicit moist physics or different types of parameterized moist physics. The parameterized moist physics is primarily based on generalizations of the parameterization scheme designed by Arakawa [30] for the original 3-layer version of the UCLA general circulation model. An axisymmetric version of the Zhu et al. model was explored by Nguyen et al. [31], with follow-up studies carried out by Zhu and Smith [32,33], and Zhu et al. [34]. Together, these papers clarify the roles of shallow convection, precipitation-cooled downdrafts, and the vertical transport of momentum by deep convection on the dynamics of tropical cyclone intensification. Another family of simple hurricane models, not reviewed here, is the Maximum Potential Intensity (MPI) family first proposed by Emanuel [35][36][37] to help understand the air-sea interaction aspects and the finite amplitude and steady-state aspects of tropical cyclones. An extensive review of the MPI models and their use in understanding the basic intensification process and longer-term evolution of tropical cyclones can be found in recent reviews by Smith and Montgomery [38,39].
It is interesting to note that the family of models defined in Tables 1-4 and summarized in Figure 2 has been organized in an unconventional manner. In the conventional organization (for example with community models such as WRF) all members of the family have the same dry dynamical core but have a variety of choices for the parameterized moist physics and air-sea interaction. In contrast, the four models in Tables 1-4 and Figure 2 all have the same parameterized moist physics and air-sea interaction but have four different choices for the dry dynamical core. These two approaches to the organization of model families seem complementary and useful for helping understand the physics of tropical cyclones.
Other extensions of the four models listed in Tables 1-4 have been described by Bliss [40] and DeMaria and Schubert [41], who relaxed the axisymmetric and balanced assumptions and formulated the three-dimensional (spectral) problem on a middle latitude β-plane to simulate storm motion. As another extension, DeMaria and Schubert [42,43] have formulated a generalized version of the axisymmetric PE model (Table 1) in order to apply the concepts of nonlinear normal mode analysis and initialization to the hurricane problem. The four models listed in Tables 1-4 have no explicit thermodynamics, which at first sight seems problematic for a model of a phenomenon that involves latent heat release. However, DeMaria and Pickle [3] have shown that a model with explicit thermodynamics, if formulated in terms of isentropic layers, is almost mathematically equivalent to the incompressible layer model. The most complete comparisons of the four models of Tables 1-4, along with several additional refinements, were given by Frisius and Lee [8], who performed numerical integrations that were, in the core region, twenty times the resolution (0.25 km vs. 5 km) used by Ooyama. Their results show additional fine structure and generally confirm the coarser resolution results shown in Figure 3.
The PE model described in Table 1 is not the only one of the four models that can be generalized to three dimensions. However, the generalizations of HB1, HB2, and GB to three dimensions can become somewhat complicated. The reason is that these involve generalizations of the classic ω-equation of quasi-geostrophic theory, which is based on the selective use of the simplest relation between the wind and mass fields. This quasi-geostrophic ω-equation is not valid for the inner core regions of intense tropical cyclones because of the large Rossby numbers found for the highly curved flows in the core. More general ω-equations that are valid for the high Rossby number flows in tropical cyclones have been derived by Arakawa [44], Thompson [45,46], and Shapiro and Montgomery [47]. These arguments are based on generalizations of gradient wind balance, e.g., on the nonlinear balance equation in the derivations of Arakawa and Thompson and on asymmetric balance in the derivation of Shapiro and Montgomery. These generalized ωequations yield fully three-dimensional formulas for local inertial stability and local Rossby length. These second-order, elliptic, ω-equations have the potential to provide diagnostic tools for a better understanding of such problems as the role of large-scale vertical wind shear on tropical cyclone structure and intensity change.
In closing, it is interesting to note that simple axisymmetric models have greatly expanded our understanding of tropical cyclone dynamics. For additional research that illustrates the rich history of these (and related models), see the studies by Emanuel [48,49], Shapiro [50], Camp and Montgomery [51], Smith and Vogl [10], Smith and Montgomery [52], Schecter and Dunkerton [53], Schecter [54,55], and Frisius et al. [56]. In many scientific fields, there is a tendency for the simplest models to experience a phase in which they yield new, remarkable physical insights (e.g., the Bohr model of the hydrogen atom, the Eady model of baroclinic instability, and the Ooyama model of tropical cyclones). Later, these simple models are used less in the research literature but are still very useful for introducing students to basic concepts. Although the simple family of tropical cyclone mod-els described in Tables 1-4 and Figure 2 may be entering their pedagogic phase, they may still prove useful in providing an increased understanding of the remarkable inner-core, lower-tropospheric structures observed in intense storms such as Hurricanes Hugo [20] and Patricia [57]. as a product of ru 0 h 0 and u 0 and then make use of the continuity Equation (6) to obtain