On the instabilities of tropical cyclones generated by cloud resolving models

Abstract An approximate method is developed for finding and analysing the main instability modes of a tropical cyclone whose basic state is obtained from a cloud resolving numerical simulation. The method is based on a linearised model of the perturbation dynamics that distinctly incorporates the overturning secondary circulation of the vortex, spatially inhomogeneous eddy diffusivities, and diabatic forcing associated with disturbances of moist convection. Although a general formula is provided for the latter, only parameterisations of diabatic forcing proportional to the local vertical velocity perturbation and modulated by local cloudiness of the basic state are implemented herein. The instability analysis is primarily illustrated for a mature tropical cyclone representative of a category 4 hurricane. For eddy diffusivities consistent with the fairly conventional configuration of the simulation that generates the basic state, perturbation growth is dominated by a low azimuthal wavenumber instability having greatest asymmetric kinetic energy density in the lower tropospheric region of the inner core of the vortex. The characteristics of the instability mode are inadequately explained by nondivergent 2D dynamics. Moreover, the growth rate and modal structure are sensitive to reasonable variations of the diabatic forcing. A second instability analysis is conducted for a mature tropical cyclone generated under conditions of much weaker horizontal diffusion. In this case, the linear model predicts a relatively fast high-wavenumber instability that is insensitive to the parameterisation of diabatic forcing. The prediction is in very good quantitative agreement with a previously published analysis of how the instability develops in a cloud resolving model on the way to creating mesovortices slightly inward of the central part of the eyewall.


Introduction
Satellite and radar images of mature tropical cyclones commonly reveal deformed eyewalls and mesovortices along the periphery of the eye. There has been longstanding interest in understanding how such features develop and whether the process appreciably affects the temporal trend of vortex intensity. One plausible explanation for the emergence of prominent waves and mesovortices involves an instability of the local circular shear flow. Although such an explanation is prevalent in the literature, there has been limited progress in advancing an instability theory for realistically modelled tropical cyclones.
Basic insights have been gained through the study of idealised two-dimensional (2D) models. Such models show that a vorticity annulus similar to that on the inward edge of an eyewall is usually unstable. The onset of perturbation growth may involve the mutual amplification of counter-propagating vortex Rossby waves or a destabilising wave-critical layer interaction. Depending on specifics, the instability may generate robust arrays of mesovortices or engender transient turbulence that thoroughly redistributes inner core vorticity into a centralised monopole (Schubert et al., 1999;Kossin and Schubert, 2001). The latter transformation may appreciably deepen the central pressure deficit while diminishing the maximum azimuthally averaged wind speed (ibid). Adding simplified parameterisations of diabatic forcing (moist convection) to a nondivergent 2D model or a shallow-water system generally modifies the development of an instability and the coinciding change of vortex intensity. Details depend on the parameterisation, and published results on the topic (Rozoff et al., 2009;Hendricks et al., 2014;Lahaye and Zeitlin, 2016) await rigorous comparison to more realistic theories and numerical simulations.
Additional insights have been gained from the study of three-dimensional (3D) stratified vortices whose basic states do not possess secondary circulations. The dominant modes of instability often resemble their 2D counterparts but differ in quantitative details [Nolan and Montgomery, 2002 (NM02)]. The qualitative similarities can extend beyond wave growth to nonlinear mesovortex formation and potential vorticity mixing (Hendricks and Schubert, 2010). On the other hand, adding vertical structure to the vortex introduces the possibility of baroclinic instability (Kwon and Frank, 2005). Moreover, instability mechanisms involving the interactions of vortex Rossby waves and inertia-gravity waves become potentially relevant in the parameter regime of a major hurricane Montgomery, 2003, 2004;Hodyss and Nolan, 2008;Menelaou et al., 2016;Schecter and Menelaou, 2017).
The final step toward a realistic perturbation theory is to generalise a 3D model to incorporate moisture and secondary circulation into the basic state of the vortex. The inclusion of cloud coverage alone has the effect of substantially reducing static stability (Durran and Klemp, 1982). In principle, such reduction can alter the structure and growth rate of the linear eigenmode associated with an instability (Schecter and Montgomery, 2007 (SM07) ;Menelaou et al., 2016). The importance of the secondary circulation to the prevailing mechanism of perturbation growth is presently unclear. Although secondary circulations are known to significantly influence the inner core instabilities of tornado-like vortices (Rotunno, 1978;Gall, 1983;Walko and Gall, 1984;Nolan, 2013), tropical cyclones are distinct atmospheric systems.
Needless to say, cloud coverage in a mature tropical cyclone is largely linked to the secondary circulation. Therefore, including one without the other in a model could yield misleading results. Naylor and Schecter (2014) (NS14) recently examined the consequences of having both. They found only subtle differences between perturbation growth in realistically simulated (moist convective) tropical cyclones and the instabilities of analogous dry (nonconvective) vortices. However, there is no firm reason to believe that the results of NS14 are general. A more extensive investigation is necessary.
NM02 contains the underpinnings of an appropriate linear model for investigating perturbation dynamics in a moist convective tropical cyclone. The NM02 model accommodates the incorporation of an adequately resolved boundary layer and the complete overturning secondary circulation of the basic state, but does not close the book on the thermodynamics. Proper parameterisation of the perturbation to diabatic forcing as a function of the prognostic fluid variables is necessary for a realistic instability analysis and remains an open issue. A separate challenge pertinent to analysing instabilities is to move beyond the conventional but questionable simplification of using constant eddy diffusivities.
Section 2 of this paper presents a somewhat distinct linear model of the perturbation dynamics that includes tuneable formulas for diabatic forcing and subgrid turbulent transport with inhomogeneous eddy diffusivities. The parameterisation of diabatic forcing does not provide a definitive closure of the thermodynamic equation, but facilitates assessment of how an instability mode might change with plausible variation in the treatment of cloud processes. Section 3 outlines a numerical method for finding the main instability modes of a tropical cyclone and the second-order response of symmetric fields to the growth of an asymmetric mode. Section 4 describes the basic state of a mature tropical cyclone generated by an axisymmetric model with explicit cloud microphysics. Section 5 analyzes the 3D instability of the aforementioned system and illustrates sensitivities to the representations of diabatic forcing and subgrid turbulence in the perturbation equations. Results of the analysis are compared to those of an ostensibly analogous 2D (barotropic) model. Section 6 presents an additional instability analysis for one of the tropical cyclones examined in NS14. The adequacy of the analysis is evaluated by direct comparison to the initial stage of perturbation growth simulated (in NS14) with a three-dimensional cloud resolving model. Section 7 summarises our main findings. The appendices provide some technical details excluded from the main text.

The perturbation equations
The present study is based on a compressible nonhydrostatic model of a tropical cyclone. The equations of motion are expressed in a cylindrical coordinate system whose central axis corresponds to that of the vortex. The radial, azimuthal and vertical coordinates are respectively denoted by r, u and z. As usual, time is denoted by t. The prognostic fluid variables are the radial velocity u, the azimuthal velocity v, the vertical velocity w, the density potential temperature h q and the total density q. Tendency equations for the mixing ratios of water vapour and hydrometeors are not explicitly considered. The influence of cloud processes on the perturbation dynamics is parameterized as explained in Section 2.2.

Basic formulation of the model
The nonlinear equations of motion governing the tropical cyclone are given by in which v is the three-dimensional velocity vector, f is the (constant) Coriolis parameter, g is the gravitational acceleration, and c pd is the specific heat of dry air at constant pressure. The Exner function satisfies the relation in which p is total pressure, p a 10 5 hPa, R d is the gas constant of dry air, and c vd is the specific heat of dry air at constant volume. Each D a represents a tendency (of field-a) induced by surface fluxes and unresolved turbulence within the vortex. S h represents the tendency of h q induced by cloud processes, radiative transfer and dissipative heating. S q is the density tendency attributable to mass changes of water content. Standard notations have been used for the gradient operator r ro r þûr À1 o u þ zo z and the divergence divðhÞ used interchangeably with o ox in this paper to denote a partial derivative with respect to any variable x.
A generic field F may be written as follows: F F b ðr; zÞ þ F 0 ðr; u; z; tÞ, in which the subscript b denotes the component associated with a suitably defined basic state of the vortex. The preceding decomposition may be applied to both the fluid variables fu; v; w; h q ; q; Pg and the forcing functions fD u ; D v ; D w ; D h ; S h ; S q g in the nonlinear model [Equations (1a)-(1f)]. The result is a perturbation equation for each prognostic fluid variable of the form in which L F consists of terms linear in F 0 and all other perturbation fields, N L F represents nonlinear terms of higher order in the perturbation amplitude, and B F accounts for residual terms involving only basic state variables along with -g in the vertical velocity equation. Ideally, the basic state is chosen to be sufficiently close to equilibrium that the magnitude of B F is no greater than second-order in the perturbation amplitude. Neglecting the relatively small terms B F and N L F reduces the dynamics to o t F 0 % L F . The azimuthal symmetry of the basic state facilitates an azimuthal Fourier decomposition of the reduced system. Letting F 0 ¼ P 1 n¼À1 F n ðr; z; tÞe inu for all F yields To stay within the realm of standard practice, the arbitrary function h qba is equated to h qb ðr B ; zÞ, in which r B is the outer boundary radius of the model. To prevent artificial trapping of radiated waves in a finite domain, we have letD un D un Àcu n and likewise for all otherD-functions. Similarly, we have letS qn S qn Àcq n . In the preceding definitions, the terms proportional to cðr; zÞ represent sponge-damping near r B and near the upper vertical boundary of the model. By design, the positive function c is negligible inside the tropical cyclone. To simplify matters, the parameterisations utilised for this study restrict D an and S an (for all applicable a) to be linear functions of the wavenumber-n components of the prognostic fluid variables. It follows that Equations (3a)-(3e) constitute an autonomous linear system. Note that the reality condition F Àn ¼ F Ã n eliminates the need to explicitly solve for the negative-n Fourier components. Here and elsewhere, the superscript 'Ã' denotes the complex conjugate of the dressed variable.
The feedback of an asymmetric linear perturbation on the mean vortex is essentially a second-order 'eddy forcing' of the symmetric (n ¼ 0) fluid variables. The full nonlinear equation of motion for a symmetric perturbation field (F 0 ) is schematically given by in which L F 0 is the right-hand side of the linear equation for F 0 [see Equations (3a)-(3e)], and the leading order contribution to N LA F (N LS F ) is quadratic in the asymmetric (symmetric) component of the perturbation. The primary quadratic part of the asymmetry term is conveniently written as follows: The summands are given specifically by Þ=c vd is the linear approximation of P m in terms of prognostic fluid variables, The operators R and I in Equations (5b)-(5f) respectively yield the real and imaginary parts of their operands. If every F 0 is initially subdominant to the asymmetric perturbation, N LS F will be negligible for an extended period of time. Forthcoming analysis of wave-mean flow interaction will set both N LS F and B F to zero in Equation (4). The latter approximation goes beyond that made in the reduced linear model for symmetric perturbations [(3a)-(3e) with n ¼ 0] by assuming that B F is much smaller than a second-order correction to the dynamics.

Parameterisation of the influence of moisture
The definition of our chosen thermodynamic variable [h q p=ðqR d PÞ] implies that in which e R d =R v , R v is the gas constant of water vapour, s d ¼ c pd ln TÀR d ln p d is the specific entropy of dry air, T is absolute temperature, p d ¼ p=½1 þ q v =e is the partial pressure of dry air, q v (q t ) is the mixing ratio of water vapour (total water content), and the overdot represents a material derivative minus any tendency directly connected to small-scale turbulence. To facilitate discussion hereafter, S h will be referred to as diabatic forcing. Equation (6) indicates that S h involves more than a term proportional to the dry-air heating rate. Nevertheless, in cloudy regions of a tropical cyclone, the reasonable assumption that _ s d is of order jL v=s _ q v =Tj suggests that _ s d will largely control the sign of the sum in parentheses. Here, the symbol L v=s has been used to denote the latent heat of vapourisation/ sublimation.
To devise a parameterisation for S 0 h , one might first consider an idealised cloudy vortex governed by reversible moist-adiabatic thermodynamics with ice-only or liquidonly condensate. The diabatic forcing in such a system satisfies an equation of the form in which q vÃ is the saturation vapour mixing ratio with respect to ice or liquid. The step function H(x) is formally defined to equal unity (zero) when x is positive (negative). The subscripts on the partial derivatives with respect to pressure indicate that the specific moist-entropy (s m ) and total water mixing ratio (q t ) are held constant. The symbol h s q (h u q ) represents the functional form of h q in terms of p, s m and q t under the assumption that the air is saturated (unsaturated) and q v equals q vÃ (q t ). Appendix A provides practical formulas for both partial derivatives that appear in Equation (7).
In the preceding reversible moist-adiabatic vortex model, the perturbation to diabatic forcing can be written as follows: The rightmost term in Equation (8) involving the product of two perturbation fields presumably has minimal effect on the weak disturbances of interest (see SM07 for caveats). Furthermore, the middle term would be negligible in a cloudy vortex whose basic state had no secondary circulation. Keeping only the first term in Equation (8), assuming _ p 0 % Àq b gw 0 , and lettingṽ b Àq b gv b =o z h qb would yield There is no firm reason to believe that a parameterisation of the diabatic forcing anomaly based on Equation (9a) would be quantitatively accurate for realistic tropical cyclones that have pronounced secondary circulations with precipitating clouds of both liquid and solid hydrometeors. On the other hand, for the class of parameterisations proportional to w 0 , Equation (9a) provides a reasonable starting point for a process of systematic adjustment toward a decent fit with experimental data. Sensitivity of results will be assessed by using the more flexible formula and letting e v vary between 0 and 1.1. For the majority of calculations in this paper, v b will be evaluated assuming ice/liquid condensate above/below the freezing level in the troposphere. The reader is referred to Appendix A for details on how v b is extracted from a numerically simulated tropical cyclone, and for further commentary on the relation S 0 h / w 0 . A more general linearised parameterisation of the perturbation to S h may have the form in which F denotes a prognostic fluid variable,L F j is the j th member of a generic set of linear operators (including differential operators) acting on F 0 ; G F j is an integration kernel paired with that operator, and the volume integral is taken over the entire domain of the system. Equation (9a) can be obtained from (10) by lettingL w 1 ½w 0 ¼ w 0 , G w 1 ¼ṽ b ðr;zÞo z h qb ðr;zÞdðrÀrÞdðuÀũÞdðzÀzÞ=r, and G F j ¼ 0 for F 6 ¼ w or j 6 ¼ 1. As usual, the symbol d has been used to represent the Dirac distribution. Note that Equation (10) is somewhat restrictive; neither the integrals nor kernels involve time. On the other hand, Equation (10) includes parameterisations that relate the perturbation of diabatic forcing at a point ðr; u; zÞ in the free troposphere to the perturbation of vertical velocity at a point ½r c ðr; u; zÞ; u c ðr; u; zÞ; z c at the top of the boundary layer (z ¼ z c ). A simple example that maintains the dynamical independence of the azimuthal Fourier transforms of the perturbation fields (in linear theory) might have an integration kernel of the form paired with the operatorL w 1 ½w 0 ¼ w 0 , while G F j ¼ 0 for all other combinations of F and j. Note that we have let u c ¼ u c þ u. Exploration of the preceding type of parameterisation will be deferred to future study.
One potential deficiency of the foregoing parameterisations [Equation (9b); Equation (11)] is their neglect of any direct response of moist convection to small enhancements or reductions of surface enthalpy fluxes coinciding with surface wind speed perturbations. Such a response could be incorporated into Equation (10), but the importance of such inclusion to mature tropical cyclone instabilities is presently unclear. Note also that the parameterisation used for this study [Equation (9b)] is not designed for high frequency perturbations exemplified by ordinary acoustic oscillations. It so happens that such rapidly oscillating modes have either subdominant or negative growth rates in our applications of the linear model. Purists might reasonably argue that the fast modes should be filtered out of the dynamical system for consistency. However, filtering out the acoustic modes alone is somewhat complicated and seems to have negligible effect on the main tropical cyclone instabilities that are investigated in this paper (Appendix B).

Parameterisation of small-scale turbulence
The influence of small-scale turbulence on the velocity perturbation is parameterized with a linear eddy viscosity scheme that incorporates a modification of the oceanic surface drag. The velocity tendencies associated with turbulence can be expressed as follows: in which the momentum eddy diffusivities (K m h and K m v ) are assumed to be functions of only r and z. Azimuthal and temporal dependencies of the diffusivities are neglected for simplicity. The rz and uz components of the stress tensor appearing in Equations (12a) and (12b) are given by in which juj ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi u 2 þ v 2 p . Unless stated otherwise, the drag coefficient is given by Note that the velocity fields in all formulas pertaining to the surface stress [Equations (13a)-(13b) at z ¼ 0; Equation (13c)] are evaluated at the lowest active grid level above the ocean in our numerical version of the linear model. Specifications of K m h ; K m v , C d0 and C d1 are forthcoming.
Several remarks are in order regarding the preceding representation of turbulent transport in the velocity equations. To begin with, Equations (12a)-(12c) above the surface are equivalent to a parameterisation of the form in which D 0 i is the tendency associated with turbulence in the prognostic equation for the ith component of the velocity perturbation (v 0 i ) in a Cartesian coordinate system Bear in mind that such a parameterisation does not follow from direct linearisation of a typical nonlinear model. Direct linearisation would produce additional terms accounting for perturbative variations of the eddy diffusivities. Note also that the usual density factors have been neglected. Despite such imperfections, Equations (12a)-(12c) are believed to provide a reasonable framework for estimating how inhomogeneous anisotropic turbulent viscosity should influence the perturbation dynamics.
Moving on to the thermodynamic equation, the effect of small-scale turbulence on h 0 q is parameterized by in which K h h=v depends only on r and z. For simplicity, the perturbation to the surface flux of h q is set to zero (see Section 2.5). The loose application of a simple diffusion scheme to the density potential temperature perturbation is deemed adequate for the present study. It is provisionally assumed that any subtle imprecision in formally representing D 0 h by Equation (15) does not affect an instability analysis more than moderate variation of the e k parameter defined below.
Several remaining formulas are required to complete the turbulence parameterisation in the linear system. To begin with, the momentum eddy diffusivities are given by in which 'max' returns the greater of its two arguments at each point in the r-z plane. The variables K m h;sm ðr; zÞ and K m v;sm ðr; zÞ in Equation (16) are obtained directly from the simulation (sm) that generates the tropical cyclone under consideration. In particular, they correspond to the horizontal and vertical momentum eddy diffusivities averaged over u (if the simulation is 3D) and over the time period that is used to define the basic state. The multiplier e k is allowed to deviate from unity for sensitivity tests. The constants K m h;min and K m v;min are lower limits of the diffusivities to be specified in due course. The previously unspecified parameters associated with the drag coefficient are given by C d0 ¼ 0:001e k and C d1 ¼ 0:0024e k for the primary instability analysis in Section 5 of this paper. The preceding formulas permit consistency with the simulation that generates the basic state when e k ¼ 1 (see Section 4). For further consistency, the thermal eddy diffusivities are given by K h h=v ¼ K m h=v , so that the Prandtl number is unity throughout the domain of the linear model.

Additional simplifications
Perturbations to radiative transfer and dissipative heating are neglected in forthcoming sections of this paper. The potential impact of radiation on the development of instabilities has been examined to some extent by adding Newtonian relaxation of the form to the perturbation of diabatic forcing in several sensitivity tests. The dominant instabilities considered herein normally have shorter time scales than a typical 12-h value of the radiative adjustment time s r . Accordingly, S r hn is normally found to have negligible consequence.
The perturbation to S q in the mass continuity equation is difficult to properly model without explicit moisture equations. The present study simply lets under the provisional assumption that it is of minor consequence to the main instabilities of a tropical cyclone. Equation (18) reducesS qn to the artificial damping term activated near the upper and outer boundaries of the domain of the dynamical system.

Boundary conditions
The linear model employs a standard set of boundary conditions for a fluid in a rigid cylindrical enclosure. At r ¼ 0, in which d nm equals 1 for n ¼ m and is otherwise 0. At the outer boundary radius r B , u n ¼ 0, o r ðv n =rÞ ¼ 0, and o r F n ¼ 0 for F n 2 fw n ; h qn g. At the surface and upper boundary (z ¼ 0 and z B ), w n ¼ 0 and o z F n ¼ 0 for F n 2 fu n ; v n ; h qn g. Consistent boundary conditions for q n are implicit in the linear model. Note that all constraints imposed on the perturbation fields at r ¼ r B and z ¼ z B are incidental when sponge damping is activated. Note further that the free-slip conditions ðo z u n ¼ o z v n ¼ 0Þ at z ¼ 0 are replaced with the surface drag parameterisation when C d > 0. As a final remark, the velocity fields of the basic state are assumed to satisfy

General theory
Let x n denote the state vector of the linearised system, with each element representing the value of one of the prognostic perturbation fields (u n , v n , w n , h qn , or q n ) at a specific point in the r-z plane. In practice, each field F n is represented on a grid with N rF points in r and N zF points in z. It follows that x n has a total of N t P F N rF N zF elements, in which the sum is over all 5 prognostic variables. The preceding discretization transforms the continuous linear model [Equations (3a)-(3e)] into a system of the form in which M n is an N t Â N t non-Hermitian matrix of complex coefficients. The eigenmodes of the discretized linear system are solutions to Equation (20) of the form in which X k is the time-independent complex eigenvector associated with the complex eigenfrequency k k R þ ik I , and a k is an arbitrary complex amplitude. Substituting Equation (21) into Equation (20) and switching the left and right sides yields Under ordinary circumstances, there are N t independent solutions to Equation (22) composing a complete eigenbasis of the wavenumber-n linear system. It follows that the solution to a generic initial value problem can be written The symbol x i (y Ãi ) in Equation (23c) denotes the i th element of x (y Ã ). The symbol M † n in Equation (23d) represents the conjugate transpose of the matrix M n . The eigenmode associated with the greatest positive value of k R will eventually dominate the right-hand side of Equation (23a). Should there exist no eigenmodes with positive k R , transient or sustained algebraic growth of the perturbation may still occur (Smith and Rosenbluth, 1990;Nolan and Farrell, 1999;Antkowiak and Brancher, 2004). Examination of such nonexponential growth in the linear model at hand is deferred to future study.
So as not to be lost in abstraction, it is worth remarking that the physical perturbation corresponding to a complex eigenmode is usually given by 2R½a k X k e inuþkt . In other words, if a k X k ðu nk ; v nk ; w nk ; h nk ; q nk Þ, the physical perturbation has the form u 0 ¼ 2ju nk ðr; zÞj cos fnu þ k I t þ arg½u nk ðr; zÞge kRt and likewise for all other fields. The coefficient 2 is replaced by 1 if n ¼ 0 and both k and a k X k are real. Suppose that the system is initially perturbed with a single asymmetric (n 6 ¼ 0) eigenmode. Consideration of Equation (4) suggests that the discretized symmetric component of the disturbance will be governed by in which b k / ja k j 2 is the time-independent part of a forcing vector obtained by evaluating the right-hand side of Equation (5a) with the eigenmode solution x m ¼ a k X k e kt for m ¼ n and x m ¼ 0 otherwise. In addition to neglecting N LS F and B F , the foregoing simplification of Equation (4) assumes that all asymmetric modes initialised to zero remain subdominant over the time period of interest. Equation (24a) is readily solved by the method of Laplace transforms and the calculus of residues after expanding both sides in the eigenvectors fX m g of M 0 . The result for x 0 ¼ 0 at t ¼ 0 is given below: in which The second term on the right-hand side of Equation (24b) will eventually dominate if m R < 2k R for all eigenfrequencies fmg of the linear symmetric system. The second term is merely the particular solution to Equation (24a) given by The particular solution is considered herein to be the intrinsic response of the mean vortex to an asymmetric instability mode. It is reasonable to consider the intrinsic response to be an essential part of the mode itself.

Computation of the main instability modes
Each fluid variable in the linearised model is discretized on a rectangular grid in the r-z plane with nonuniform spacing in both coordinates, as in earlier studies such as NM02. Finer resolution generally exists near the surface and within the core of the tropical cyclone. The discretized representations of v n , h qn and q n share the same grid. The representation of u n (w n ) is radially (vertically) staggered with respect to v n . The basic state variables and eddy diffusivities are defined on all of the staggered grids. Boundary values are not explicitly stored but are incorporated into computations where necessary.
The following simple formulas are normally used for finite differencing and linear interpolation of a generic field F n : in which ' represents either r or z, ' 0 denotes ' at the evaluation point, d' 6 is the distance from ' 0 to the nearest staggered grid point in the positive (þ) or negative (À) direction, and F 6 n is the value of F n at ' 6 ' 0 6d' 6 . For example, if ' represents r (z) and ' 0 is on the v-grid, then F 6 n and ' 6 are on the u-grid (w-grid). Formulas for secondorder derivatives and bilinear interpolations are generally obtained through repeated applications of Equations (25a) and (25b). Implementation of more accurate discretization techniques will be explored at a future time.
The computation of the complete eigenbasis of a finelystructured tropical cyclone is usually too expensive to achieve with confidence of correct results. Although M n is sparse and has a storage requirement proportional to N t , the eigenbasis fX k g has a storage requirement proportional to N 2 t . The consequent demand on memory becomes difficult to handle for grids comparable to those used in modern tropical cyclone simulations. Furthermore, the time required to compute a complete eigenbasis on a modern simulation grid is excessive. Grids of lower resolution should be avoided, because they are prone to introduce spurious eigenmodes with dominant growth rates. Moreover, grids of higher resolution are desirable to check for convergence of the numerics.
The present study employs a less ambitious approach that begins by extracting the dominant eigenmode from a solution of the initial value problem. The discretized linear model [Equation (20)] is set up on a dense mesh [see Section 5.1] and integrated forward in time with a 4thorder Runge-Kutta algorithm. The initialisation involves assigning small random values to the real and imaginary parts of h qn at each grid point; all other fields contained in x n are initialised to zero. It is provisionally assumed that the preceding disturbance excites the main instability modes of a tropical cyclone and eventually evolves into a state dominated by the most unstable member of the group. The real and imaginary parts of the eigenfrequency k of the most unstable eigenmode are readily obtained from the late time series of a selected element of x n . The right-hand eigenvector X k is very well approximated by the late spatial structure of x n . The validity of the mode is generally cross-checked against the output of a standard sparse-matrix eigensolver (eigs) packed into Scientific Python (SciPy). Validation is efficiently completed by searching exclusively for the eigenmode of M n with k closest to that obtained from the initial value problem. The SciPy eigensolver is also used to find the corresponding left-hand eigenvector X L k . Because the restricted searches are fast, they are usually repeated on a grid with twice the original resolution (in both r and z) to slightly improve the accuracy of presented results.
Suppose that the eigenfrequencies fk ðaÞ g are ordered such that k R . . . In principle, if all eigenmodes with a < b are known, a minor variant of the foregoing procedure can be repeated to obtain eigenmode b. The variant involves filtering out all eigenmodes with a < b from the initial condition of the state vector that is integrated forward in time; that is, letting in which Y is an arbitrary vector. One may reasonably assume that the time asymptotic solution of x n will be dominated by the eigenmode labelled b. All eigenmodes of interest can thus be found iteratively. Note that the unfiltered initialisation vector Y need not be random after the first iteration; the approach taken here is to let Y equal the end-state of x n from the preceding time integration used to find the eigenmode labelled bÀ1.

The basic state of a mature tropical cyclone
The primary basic state considered herein corresponds to a mature tropical cyclone simulated with Cloud Model 1 (CM1-r19.4) in an energy-conserving axisymmetric mode of operation [Bryan and Fritsch, 2002;Bryan and Rotunno, 2009 (BR09)]. The model is configured with a variant of the two-moment Morrison microphysics parameterisation (Morrison et al., 2005(Morrison et al., , 2009, having graupel as the large icy-hydrometeor category and a constant cloud-droplet concentration of 100 cm -3 . Radiative transfer is not explicitly calculated, but potential temperature (h) is relaxed toward its ambient value on a 12-h time scale with a rate not to exceed 2 K d -1 in magnitude. The influence of subgrid-turbulence above the surface is represented by an anisotropic Smagorinsky-type scheme resembling that described in BR09. The nominal mixing lengths are given by CM1-formulas tailored for tropical cyclones in an axisymmetric framework or on grids that are deemed insufficiently dense for a standard large-eddysimulation scheme. The horizontal mixing length increases from 100 m to 1 km as the underlying surface pressure decreases from 1015 to 900 hPa. The vertical mixing length increases asymptotically to 100 m with increasing z. The resulting eddy diffusivities will be discussed in due course. Heating associated with frictional dissipation is activated. Surface fluxes are parameterized with bulk-aerodynamic formulas. The drag coefficient conforms to Equation (13c) with C d0 ¼ 0:001 and C d1 ¼ 0:0024, based roughly on the findings of Fairall et al. (2003) and Donelan et al. (2004). The enthalpy exchange coefficient is given by C e ¼ 0:0012 based on Drennan et al. (2007). The radial vorticity (contours) and azimuthal vorticity (colour) distributions. The solid yellow and dashed white curves in (b)-(e) correspond to the principal AM isoline, which is commonly shown for spatial reference in the contour plots of this paper. Note that v bm is located where the radius of the isoline is minimised.
The computational domain extends radially to r B ¼ 1061:25 km and vertically to z B ¼ 29:5 km. Free-slip boundary conditions are imposed at the upper and outer boundaries, but they are largely incidental. Rayleigh damping is activated above z ¼ 25 km and within 100 km of r B . Such damping not only minimises the reflection of upward and outward propagating waves, but also prevents the circulation of the tropical cyclone from expanding to the outer wall. The spatial resolution is fairly typical for modern tropical cyclone simulations. The radial grid spacing is 250 m for r 87:5 km, and gradually stretches to 10 km as r approaches the boundary radius r B . The vertical grid spacing increases from 50 to 500 m between the sea-surface and z ¼ 5.5 km, whereupon it remains uniform up to z B .
The mature tropical cyclone is generated with a standard spinup procedure on the oceanic f-plane. The model is initialised with a surface-concentrated mesoscale cyclone in gradient-wind and hydrostatic balance as in NS14. The ambient atmosphere is initialised with the Dunion (2011) moist tropical sounding. The sea-surface temperature T s is 27 C, and the Coriolis parameter f is 5 Â 10 À5 s À1 . After approximately 7 days of intensification, the maximum azimuthal velocity of the tropical cyclone remains steady over an extended period of time. The basic state variables Þ and principal eddy diffusivities ðK m h;sm ; K m v;sm Þ appearing in the linear model are obtained by averaging 25 consecutive hourly snapshots starting 2 days into the aforementioned period of steady intensity. Figure 1a shows that the basic state of the tropical cyclone exhibits typical warm-core baroclinic structure throughout most (but not all) of the troposphere. The absolute maximum azimuthal wind speed (v bm ) of 84.9 m s À1 is located 36 km from the centre of the vortex and nearly 1 km above the surface. The maximum azimuthal wind speed of 61.2 m s À1 at the lowest grid level (z ¼ 25 m) is indicative of a category-4 hurricane. It is worth remarking that the primary circulation does not robustly satisfy gradient balance (Fig. 1b). The fractional error defined by is most pronounced (66%) in the vicinity of v bm . Figure 1c shows that throughout the lower and middle troposphere, the potential vorticity distribution is generally peaked off centre within the area bounded by the principal angular momentum (AM rv þ fr 2 =2) isoline. Here, the potential vorticity is defined by PV f a Á rh=q, in which f a fẑ þ r Â v is the absolute vorticity vector. The principal AM isoline is defined so as to pass through v bm . By analogy to the behaviour of dry vortices in gradient and hydrostatic balance, the radially nonmonotonic variation of PV suggests that the tropical cyclone is susceptible to vortex Rossby wave instability mechanisms (see Montgomery and Shapiro, 1995). The pocket of negative PV extending up to 6 km above sea level slightly outward of the principal AM isoline suggests that (neglecting viscous dissipation) the vortex may also be susceptible to inertial instability mechanisms in the lower-to-middle tropospheric region of its core (see Eliassen, 1951). Figure 1d demonstrates that the distribution of relative vertical vorticity (f zb ) in the lower and middle troposphere basically resembles that of PV. The most notable deviation is seen where the PV distribution is thermally enhanced at the top of the boundary layer in the eye of the storm. Note that the primary circulation also possesses appreciable radial vorticity associated with its vertical shear (Fig. 1e). Evidently, the radial vorticity achieves magnitudes greater than f zb near the surface.
of the surface inflow is relatively strong but not much greater than typical observations pertaining to major hurricanes (Zhang et al., 2011). The inflow intensity is greatest slightly outward of the corner flow region, where the streamlines rapidly turn upward into the eyewall cloud. Note that the azimuthal vorticity associated with the secondary circulation in the vicinity of the corner flow (Fig. 1e) is comparable in magnitude to the peak vertical and radial vorticities associated with the primary circulation. Note also that the streamline associated with the deep updraft and outflow passing through the location of v bm is virtually congruent with the principal AM isoline. Such a condition is to be expected for a nearly equilibrated axisymmetric vortex. As usual, the secondary circulation in the eye is dominated by weak subsidence. The streamlines are somewhat less coherent at larger radii between the surface inflow and upper tropospheric outflow. Concerns that such incoherence may indicate a significant departure from equilibrium are alleviated by noting that the regional wind speeds are minute compared to peak values. Figure 3 illustrates the moist-thermodynamic structure of the basic state. Contours of saturated pseudoadiabatic entropy (s pÃ ) are shown superimposed on the relative humidity distribution. Relative humidity is calculated with respect to liquid water if the absolute temperature satisfies T > T 0 273:15 K, but is otherwise calculated with respect to ice. It should come as no surprise that the eyewall and outflow regions are predominantly saturated or slightly supersaturated. The dashed black-and-white contours correspond to s pÃ for liquid-only condensate (Bryan, 2008). It is seen that the angular momentum and liquid-only s pÃ contours passing through the location of v bm are congruent as they ascend along the eyewall up to the freezing level. At higher altitudes, the angular momentum contour appears to stay closer to the dotted-blue s pÃ contour calculated under the assumption of ice-only condensate (Hakim, 2011). The preceding observations suggest that the eyewall updraft region is in a state of approximate slantwise convective neutrality with respect to appropriately defined pseudoadiabatic thermodynamics. Such a state of affairs is consistent with the classical steady state theory expounded by Emanuel (1986).
The cloud structure of the basic state is important to the linear model insofar as it determines the proportionality between S 0 h and w 0 in the local parameterisation of diabatic processes given by Equation (9b). When the aforementioned parameterisation is activated, there are two terms proportional to w n on the right-hand side of Equation (3d) that may be unified as follows: in whichÑ A typical value of e vṽb between zero and one reducesÑ 2 and thereby diminishes the negative/positive Eulerian change in h q associated with a perturbative updraft/downdraft. While not precisely the conventionally defined moist static stability,Ñ 2 has a similar significance. Figure   4 compares the distribution ofÑ 2 in the approximate dry limit (e v ¼ 0) to the moist variant with e v ¼ 1. It is seen that incorporating the cloud coverage of the basic state reducesÑ 2 up to an order of magnitude in the eyewall updraft. Significant reduction is also found over much of the depicted area within and underneath the upper outflow of the tropical cyclone. By contrast,Ñ 2 exhibits minimal change in the virtually cloud-free region of the eye situated above the boundary layer. As explained earlier, the eddy diffusivities used by the linear model are linked to those regulating the basic state. Figure 5 shows the eddy diffusivities that are defined by Equation (16) with e k ¼ 1. Choosing e k ¼ 1 lets K m h=v ¼ K m h=v;sm throughout much of the inner core. The maximum values are somewhat large but have orders of magnitude consistent with those inferred from observations (Zhang and Montgomery, 2012;Rogers et al., 2013). To reduce the potential for spurious or uninteresting small-scale instabilities where the eddy diffusivities in the CM1 simulation are

Linear instability analysis of a mature tropical cyclone
The present section of this article examines the instability of the tropical cyclone described in Section 4. The primary objective is to elucidate the dependence of the dominant instability mode on the parameterisation of the perturbation of diabatic forcing. Sensitivity to the parameterisation of small-scale turbulence is also addressed. The analysis concludes with an assessment of the relevance of 2D instability theory. A few preliminary remarks are warranted. Henceforth, the meaning of F 0 is subtly changed from the exact difference FÀF b to the first-order perturbation of the generic field F obtained from the linear model [Equations (3a)-(3e); Equation (20)]. The new meaning of F 0 applies to both figures and text. Moreover, the amplitudes of displayed instability modes are invariably chosen to render the maximum value of v 0 (2jv n j) equal to one-tenth of v bm . The preceding convention amounts to letting ja k j ¼ 0:1v bm =j2V k e kRts j, in which V k is the azimuthal velocity element of X k with the greatest magnitude, and t s is the time of the snapshot. In some cases, the second-order change to the mean vortex (x 0 ) that will have attended the creation of such a state from a weaker disturbance by way of Equation (24a) is found to have winds moderately stronger than v 0 in certain areas of the flow. Such a result indicates that the arbitrarily chosen mode amplitude is slightly beyond the threshold for the quantitative accuracy of Equation (24a). Choosing a smaller amplitude for rigorous compliance with the assumptions of our theoretical framework would not change forthcoming depictions of the spatial structure of x 0 or the dependent kinetic energy perturbation dKE defined later.
Finally, although the physics parameterisations are varied, the domain size and peripheral sponge-layer of the linear model used to find the instability modes do not change from one calculation to the next. As in the CM1 simulation used to generate the basic state, the invariant domain of the linear model extends radially to r B ¼ 1061:25 km and vertically to z B ¼ 29:5 km. The sponge damping coefficient is given by c ¼ f2 þ tanh½ðrÀr c Þ=dr c þ tanh½ðzÀz c Þ=dz c g=2s c , in which r c ¼ 961:25 km, dr c ¼ z c ¼ 25 km, dz c ¼ 0:75 km and s c ¼ 300 s. Further computational details are provided in due course.

Sensitivity to the parameterisation of diabatic forcing
The dominant instability of the tropical cyclone under consideration is sensitive to the degree of diabatic forcing allowed in the linear model. The sensitivity is illustrated below by adjusting e v in Equation (9b) for S hn while keeping turbulent transport consistently parameterized with e k ¼ 1. A value of the diabatic forcing parameter (e v ) in the neighbourhood of unity has some basic credibility (Section 2.2) but may not coincide with the best representation of reality. A smaller value between 0 and 1 seems plausible if, say, the eyewall were to become nonuniformly saturated around an azimuthal circuit. Values of e v very close to 0 or appreciably greater than 1 seem difficult to justify on physical grounds, but are of theoretical interest. The present method for computing the primary instability modes of the vortex follows the general procedure outlined in Section 3.2. The most unstable eigenmode (MUM) for a given e v and azimuthal wavenumber n is provisionally equated to that which dominates a perturbation within 1 day of initialising the linear model [Equations (3a)-(3e)] with random noise in h 0 q . The absolute MUM (AMUM) is defined to be that which possesses the largest growth rate for all n in the closed interval between 0 and 8. The time integration is conducted on a grid (set of staggered grids) with double the resolution of the CM1 simulation that generated the basic state. The aforementioned grid is denoted G2 and holds N t ¼ 733; 184 values of the prognostic perturbation fields. All MUMs are confirmed to be solutions of the eigenproblem on a second grid (G4) with quadruple the resolution of the CM1 grid (G1). All eigenfrequencies and eigenfunctions shown in Section 5 of this paper are taken from the G4 solutions. For those interested, Appendix C discusses convergence of numerical results with increasing resolution.
Extensive computations reveal that the AMUM corresponds to n ¼ 2 for e v 2 f0; 0:5; 1g. Despite its common dominance, the n ¼ 2 MUM varies considerably with the allowed degree of diabatic forcing measured by e v . Figure  6 shows the variation of the complex eigenfrequency k. The growth rate (k R ) gradually decays with increasing e v until apparently vanishing at 0.9. By contrast, the oscillation frequency (k I ) changes little. The preceding behaviour is similar to that reported by SM07 for the n ¼ 3 MUM of a cloudy vortex resembling a category-3 hurricane with no mean secondary circulation. On the other hand, increasing e v from 0.9 to 1 introduces a new mode of instability that oscillates slower and grows faster than any of its predecessors. Further amplification of e v to 1.1 substantially increases both k R and jk I j. One might reasonably speculate that high sensitivity to variation of e v in the neighbourhood of unity is related to approximate slantwise convective neutrality with respect to pseudoadiabatic thermodynamics in the eyewall (Fig. 3). Figure 7 shows the basic inner-core structure of the n ¼ 2 MUM for values of the diabatic forcing parameter below (e v ¼ 0:5) and above (e v ¼ 1) the apparent stability point. The left column shows selected views of the asymmetric velocity perturbation. The middle column illustrates the thermal structure of each mode in terms of h 0 q and P 0 . The right column shows the distributions of diabatic forcing. The velocity perturbations of the two modes are qualitatively similar near the surface but clearly differ aloft. Whereas the pressure perturbations seem only subtly distinct, disparities in h 0 q are pronounced. Marked distinctions in the perturbations of the secondary circulation and potential temperature in the middle and upper troposphere coincide with substantial differences in S 0 h . Not only does S 0 h have a greater amplitude in the MUM corresponding to e v ¼ 1, but the two spatial patterns diverge considerably above 4 km in the eyewall updraft region of the vortex. Figure 8 elaborates on the inner-core structure of each MUM. The left column shows the intensity of the vertical vorticity perturbation f 0 z , as measured by its maximum solid contours: K (m 2 s 1 ) Fig. 5. Horizontal (colour) and vertical (solid contours) momentum eddy diffusivities in the middle-to-lower tropospheric core of the simulated tropical cyclone. The dashed curve is the principal AM isoline. Fig. 6. Complex eigenfrequency k of the n ¼ 2 MUM of the simulated tropical cyclone of Section 4 versus the diabatic forcing parameter e v . The real (blue) and imaginary (red) parts of the eigenfrequency are normalised to their respective values (k RÃ ¼ 7:89 Â 10 À5 s À1 and k IÃ ¼ À1:30 Â 10 À3 s À1 ) obtained for e v ¼ 1. The absence of a discernible instability precludes the plotting of data for e v ¼ 0:9. Note that the positive ratio k I =k IÃ represents a nondimensional magnitude of the oscillation frequency; the actual value of k I is negative. All results depicted here and in Figs. 7-12 are for systems in which turbulent transport is parameterized with e k ¼ 1. value over u. In each MUM, the intensity peaks of f 0 z roughly coincide with a subset of regions where the radial gradient of basic state potential vorticity is locally enhanced (see Fig. 1c). The amplitudes of the peaks differ considerably between the two modes, especially in the middle-to-upper troposphere. The middle column depicts the maximum magnitude of the horizontal vorticity perturbation f 0 h along an azimuthal circuit. In both MUMs, jf 0 h j broadly exceeds the vertical vorticity perturbation. As before, differences between the two MUMs are mainly seen in the amplitudes of various peaks of the plotted field. The right column shows the circuit-maximum of the horizontal divergence, defined by r 0 ½o r ðru 0 Þ þ o u v 0 =r. In both MUMs, r 0 is broadly smaller than the vertical vorticity perturbation near the surface, but is far from negligible. In the middle tropospheric region of the eyewall cloud, the amplitudes of r 0 and f 0 z are comparable to each other. The MUM corresponding to e v ¼ 1 is distinguished by having a middle tropospheric peak of r 0 that slightly exceeds the inner core maximum of f 0 z . Figure 9 illustrates for each MUM how the azimuthal phase velocity minus the local angular velocity of the primary circulation (c u Àk I =n À X b ) varies over the core of the tropical cyclone. The dashed green curves representing the zero contours of c u correspond to where the mode corotates with the mean flow. Negative/positive values of c u indicate locally retrograde/prograde wave propagation in the azimuthal direction. The superimposed vertical vorticity distribution (solid black contours) of the MUM corresponding to e v ¼ 1 is concentrated in the region of retrograde propagation. On the other hand, the MUM with weaker diabatic forcing has a middle-toupper tropospheric swath of intense f 0 z that extends well into the region of prograde propagation. In both cases, the magnitude of the intrinsic frequency of the mode (nc u ) is less than the nominal inertial frequency ( ffiffiffiffiffiffiffiffiffiffi g b n b p ) where the vorticity anomalies are peaked. While notable, such local slowness does not necessarily indicate that traditional asymmetric balance theory (Shapiro and Montgomery, 1993) would provide an accurate description of the wave dynamics. Bear in mind that the issue is complicated by the moist secondary circulation and the vertical shear in v b . Moreover, even small deviations from balanced dynamics are potentially important to the instability mechanism.
Moving outward to where r exceeds 100 km, the MUMs acquire intrinsic frequencies that broadly satisfy g b n b ( ðnc u Þ 2 (Ñ 2 (not shown). The preceding condition suggests that the intrinsic frequency lies comfortably within the regime of inertia-gravity waves. Consistent with such waves, one finds that jf zn j ( jr n j beyond the core of the vortex, barring sporadic pockets of violation. The right panels in Fig. 10 convey the basic structure of the outer waves as represented by w 0 in the two MUMs under consideration. Although both modes are normalised to have the same inner core maximum value of v 0 , the outer waves have appreciably stronger vertical velocities for the case in which e v ¼ 0:5. Whether such a distinction is relevant to the mechanism of modal growth is a question left for future analysis. In theory, seemingly weak inertia-gravity wave radiation may contribute significantly to the prevailing low-n instability of an intense tropical cyclone (Menelaou et al., 2016;Schecter and Menelaou, 2017). However, the author is unaware of any existing method for assessing the importance of inertia- gravity wave emission to the growth of a multifaceted instability mode of a convective vortex with the geometrical complexity of a realistic hurricane.
Figures 11a,d show changes to the mean flow that attend the growth of each instability mode from an asymptotically small disturbance. The symmetric component of the perturbation is given by x 0 ¼ X p e 2kRts , with X p given by Equation (24d). The growth of either instability mode modestly reduces the u-averaged azimuthal wind speed at the initial location of maximal intensity while accelerating the cyclonic rotation of the inner eye, at least in the lower troposphere. The middle tropospheric patterns of symmetric azimuthal acceleration and deceleration are clearly dissimilar inward of the principal AM isoline. Moreover, the MUM affected by weaker diabatic forcing (e v ¼ 0:5) induces greater positive and negative azimuthal accelerations of the mean flow in the upper-outer part of the eyewall updraft. The perturbation of the symmetric secondary circulation ðu 0 ; w 0 ) that emerges during the growth of either instability mode notably includes a band of eddies along the eyewall updraft. The bands associated with the two MUMs are distinguishable in part by having opposite rotational tendencies at various locations.
Figures 11b,e show the perturbation of kinetic energy density dKE associated with the growth of each MUM. To second-order in the asymmetric mode amplitude, in which the overline denotes an azimuthal average. It has been verified that the bottom line in the second equality involving the density perturbation is negligible (not shown). Moreover, it is seen that the distribution of dKEhere divided by -is similar to that of v 0 regardless of whether e v is 0.5 or 1. The contribution to dKE from the asymmetric fields is well approximated by the following positive definite measure of local wave intensity: KE n q b ðju n j 2 þ jv n j 2 þ jw n j 2 Þ. Figures 11c,f show the spatial distributions of KE n for the two MUMs under present consideration. Both MUMs have their greatest values of KE n near the surface, inward of the radius of maximum wind, in the vicinity of where the vertical vorticity of the basic state (f zb ) has a pronounced maximum. The middle-to-upper tropospheric peaks of KE n are found in distinct locations. Above the surface perturbation, the distribution of KE n corresponding to e v ¼ 0:5 has relatively strong peaks outward of the central part of the eyewall. The instability mode that results from allowing greater diabatic forcing (e v ¼ 1) has its principal middle tropospheric maximum of KE n well within the eyewall updraft.
Differences between the MUMs are also evident in various terms that formally contribute to the growth rate of KE n . Equations (3a)-(3c) imply that The term labelled PC combines tendencies proportional to the radial and vertical shear of the primary circulation of the basic state. SC combines tendencies proportional to u b and the spatial derivatives of the velocity fields of the secondary circulation. BNC is linked to the vertical and radial buoyancy accelerations. AFX primarily represents the convergence of the advective flux of KE n . The included correction is attributable to the small but nonzero divergence of the momentum density of the basic state. PFX primarily represents the convergence of the flux vector associated with forcing by the perturbation of the pressure-gradient. The included correction is attributable to the small but nonzero divergence of the approximate momentum perturbation weighted by h qb . TRB is associated with turbulent momentum transport and (to a lesser extent) sponge-damping near the upper and outer edges of the computational domain. It is worth pointing out that substantial cancellations of the tendency terms often result in a local value of o t KE n ¼ 2k R KE n that is much smaller than its individual parts. Figure 12 illustrates how the value of the diabatic forcing parameter (e v ) affects the volume-integrals of the KE n tendency terms pertaining to the n ¼ 2 MUM of the tropical cyclone. The volume integrals are over the entire domain of the linear model. The results are similar for all e v 0:75. The integral of PC provides the greatest positive contribution to the sum. The component of PC associated with the radial shear of the basic state is dominant. The integral of SC is smaller than that of PC, but often greater than the integral of all terms combined. The integrals of both BNC and TRB are negative and substantial. On the other hand, the integrals of AFX, PFX and their displayed sum are negligible. The budget corresponding to e v ¼ 1 has several distinctive features. The difference between the PC and SC integrals is appreciably reduced. Moreover, the BNC integral is positive. Of lesser significance, the combined integral of AFX and PFX is discernibly negative owing mostly to the corrective component of PFX. Increasing e v to 1.1 moves the strongest peaks of the asymmetric kinetic energy density from the surface to the middle troposphere (not shown). The attendant structural change coincides with notable modifications to the global KE n -budget. For example, the vertical shear component of the PC integral becomes dominant. Moreover, the integral of BNC becomes nearly equal to that of PC.

Sensitivity to the parameterisation of turbulent transport
The MUM associated with arbitrary n generally varies with the parameterisation of small-scale turbulence. Sensitivity to the intensity of turbulent transport is illustrated herein by reducing the value of e k defined in Section 2.3. The minimum value of e k to be considered will be 0.0625, which is slightly below the limit of 0.07 (0.08) that guarantees K m h (K m v ) will uniformly equal the value specified for K m h;min (K m v;min ) in Section 4. Figure 13 shows how reducing e k affects the complex eigenfrequencies of the MUM and the second most unstable eigenmode (SMUM) of linear systems with n ¼ 2 and e v 2 f0:5; 1g. Results are shown for e k ¼ 1, 0.25 and 0.0625. As before, the MUM is provisionally equated to the prevailing instability mode that emerges during a time integration of the linear model initialised with a random distribution of h 0 q on G2. The SMUM is provisionally equated to the prevailing instability mode of a continued integration that filters out the MUM [see Equation (26)]. Both modes are verified to solve the eigenproblem on G4. The displayed data are obtained from the G4 eigensolutions.
Consider first the group of linear systems that allow a medium degree of diabatic forcing (e v ¼ 0:5). Section 5.1 thoroughly described the dominant MUM when e k ¼ 1.
The corresponding SMUM has a lower oscillation frequency and is structurally distinct in having KE n concentrated in the middle troposphere (not shown). Reducing  v, light red). The TRB contribution is decomposed into the primary part attributable to turbulent dissipation (dark cyan) and the much smaller part attributable to sponge damping (light cyan cap). e k introduces a faster instability that overtakes both of the aforementioned eigenmodes. The greater growth rate (k R ) of the new MUM coincides with a greater oscillation frequency (jk I j). The new MUM is also structurally distinct in having KE n largely confined to a shallow layer near the surface (Fig. 14a). Moreover, the global KE n budget is distinguished from that of the original MUM by having a greater vertical shear component of PC, and a minimal contribution from SC (Fig. 14b).
Consider next the set of linear systems that allow relatively strong diabatic forcing (e v ¼ 1). As before, the reader may consult Section 5.1 for a thorough description of the dominant MUM when e k ¼ 1. The corresponding SMUM is similar to that of the equally diffusive system with e v ¼ 0:5. Reducing e k to 0.25 modestly accelerates the instability associated with the original MUM and leads to the appearance of a new SMUM with nearly the same growth rate. Reducing e k to 0.0625 switches the ordering of the preceding instability modes without changing their top-tier status. The new mode is distinguished by having a greater oscillation frequency and a dissimilar distribution of KE n above the boundary layer (Fig. 14c). Moreover, the global KE n budget of the new mode is distinguished by having a greater vertical shear component of PC, and a negative contribution from BNC (Fig. 14d).
It is worth remarking that decreasing the eddy diffusivity often magnifies the importance of higher wavenumber MUMs. For example, reducing e k to 0.25 in a system with e v ¼ 1 allows an n ¼ 3 MUM (Figs. 14e,f) to challenge its n ¼ 2 counterpart for dominance among instability modes with substantial KE n near the surface. While the former oscillates approximately 1.6 times faster than the latter, both MUMs have growth rates of 1:1 Â 10 À4 s À1 .

Relationship to 2D instability theory
It is common practice to explain the instability of the primary circulation of a tropical cyclone in the context of a two-dimensional nondivergent barotropic model (see Appendix D). The foregoing analysis casts doubt on the general adequacy of such an approach. That is to say, the preceding results suggest that the three-dimensionality of the tropical cyclone under present consideration has a major impact on the prevailing mode of instability. The evidence includes MUMs with substantial horizontal vorticity and divergence. The evidence also includes major contributions from SC and/or the vertical shear component of PC to the volume integrated time-derivative of asymmetric kinetic energy (KE n ).
Further insight is gained by directly comparing 2D and 3D instability theory. The 2D analysis requires reduction of the basic state to a circular shear-flow characterised by a 1 D vertical vorticity profile f b ðrÞ. Because the asymmetric kinetic energy density of the instability usually has greatest amplitude in the lower troposphere, f b will be extracted from the q b -weighted vertical average of f zb ðr; zÞ (Fig. 1d) between the sea-surface and z ¼ 2 km. The kinematic viscosity K 2d will be varied between 0 and 4000 m 2 s -1 . The upper limit is roughly 1.4 times the peak value of K m h in the 3D model when e k ¼ 1 (Fig. 5). The nonmonotonic radial variation of f b facilitates a variety of algebraic and exponential instabilities. An algebraic instability is expected to dominate the n ¼ 1 component of an arbitrary disturbance (Smith and Rosenbluth, 1990). The exponentially growing eigenmodes associated with greater azimuthal wavenumbers are readily obtained from a complete numerical solution to the eigenproblem on a stretched radial grid comparable to that of G2. For and imaginary (red) parts of each eigenfrequency are normalised to their respective values (k RÃ ¼ 7:89 Â 10 À5 s À1 and k IÃ ¼ À1:30 Â 10 À3 s À1 ) obtained for the MUM when e v ¼ e k ¼ 1. K 2d ¼ 0, the AMUM corresponds to n ¼ 3, but all MUMs with 2 n 5 have growth rates within 8% of the maximum (Fig. 15a). Increasing K 2d toward 4000 m 2 s À1 diminishes the growth rate of each MUM with greater effect at larger n; ultimately, the exponential instabilities are confined to azimuthal wavenumbers 2 and 3. The preservation of n ¼ 3 dominance (or shared dominance) with increasing viscosity appears to be at odds with the 3D model. For e k ¼ 1 and e v 2 f0; 0:5; 1g, the wavenumber-3 instability modes of the 3D model were found to be subdominant. Figures 15b,c show how the complex eigenfrequencies of the two most unstable n ¼ 2 eigenmodes vary with K 2d .
The two modes are distinguished by their virtually invariant oscillation frequencies that differ roughly by a factor of 2. Decreasing the viscosity from its maximal value is seen to unleash the instability of the high-frequency mode, such that it transitions from SMUM to MUM status as K 2d drops below 2500 m 2 s -1 . Despite the reordering of growth rates, neither the low-frequency mode (Fig. 16a) nor the high-frequency mode (Fig. 16b) radically changes structure with variation of K 2d over the interval under consideration. Except for moderate radial smoothing of the vorticity wavefunction, the unshown modifications linked to greater viscosity are difficult to discern with a casual glance.
To some extent, the low-frequency mode of the 2D system that prevails under conditions of high viscosity resembles the lower tropospheric section of a typical 3D MUM that dominates under moderate diabatic forcing when turbulent transport is parameterized with e k ¼ 1. Figure 16c (16e) depicts the lower tropospheric structure of the 3D MUM corresponding to e v ¼ 0:5 (1). As in the low-frequency mode of the 2D system, the strongest perturbation eddies are centred on the outer edge of the main vorticity annulus. In similar agreement, the prominent inner and outer waves of f 0 z are close to being diametrically out of phase at azimuths where the amplitudes are peaked. On the other hand, seemingly subtle differences cannot be ignored. To begin with, the radii at which the 2D and 3D modes corotate with the circular shear flow (shown by the dashed green circles) do not coincide. In principle, even a slight displacement of a corotation radius can substantially affect the impact of locally enhanced (potential) vorticity stirring on the growth of an instability mode. The nature of any delicate imbalance of various growth and decay mechanisms may also be sensitive to small variations in the relative amplitudes and phases of the primary inner and outer vorticity waves. Moreover, the horizontal velocity perturbations associated with the depicted 3D instability modes have nonnegligible divergence. Figures 16d,f illustrate the divergent (irrotational) components of the modal flow fields obtained from a standard Helmholtz decomposition as explained in Appendix E. The maximum divergent wind speed for e v ¼ 0:5 (1) is an appreciable 16% (27%) of the maximum nondivergent wind speed.
It is notable that (for n ¼ 2) the low-frequency instability modes of both the 2D system and the 3D systems studied in Section 5.2 are superceded by higher frequency modes as viscosity tends toward zero. The low-viscosity 3D MUM corresponding to e v ¼ 0:5 and e k ¼ 0:0625 (Fig. 14a) is fairly similar to its 2D counterpart (Fig.  16b). To begin with, the 3D MUM is confined to a shallow layer near the surface. Moreover, unshown analysis of the horizontal flow in the lower troposphere Computed MUMs with growth rates of order 10 À6 s À1 or less (at high n and appreciable K 2d ) are excluded from the plot, because they are considered virtually neutral and have questionable accuracy. (b) K 2d dependencies of the growth rates of the two most unstable 2D eigenmodes associated with an n ¼ 2 perturbation. (c) As in (b) but for the oscillation frequencies. The extended blue ticks on the vertical axes of (b) and (c) mark the growth rate and oscillation frequency of the n ¼ 2 MUM of the 3D system with e v ¼ 0:5 and e k ¼ 1; the red ticks are the same but for e v ¼ e k ¼ 1. 2D low-freq mode 2D high-freq mode Fig. 16. (a) Vertical vorticity (red and blue), streamlines (black) and corotation circles (dashed green) of the low-frequency mode of the 2D system with K 2d ¼ 10 3 m 2 s À1 . The streamline thickness is directly proportional to the local magnitude of the horizontal velocity perturbation u 0 . (b) As in (a) but for the high-frequency mode. (c) As in (a) but for vertically averaged fields associated with the MUM of the 3D system with e v ¼ 0:5 and e k ¼ 1; the averaging is over a 2 km layer adjacent to the sea-surface. (d) As in (c) but with the streamlines corresponding to the irrotational component of u 0 . (e, f) As in (c, d) but for e v ¼ e k ¼ 1; note that segments of the corotation circle can be found at the corners of both plots. In all subfigures, the axis labels x and y denote horizontal Cartesian coordinates measured from the center of the vortex.
demonstrates that the strongest perturbation eddies are centred on the inner edge of the main vorticity annulus, and that a corotation radius lies in between the primary inner and outer vorticity waves. In good agreement with a key assumption of the 2D model, the maximum magnitude of the divergent component of the lower tropospheric velocity perturbation (averaged over a 2-km layer adjacent to the sea-surface) is merely 6% of the maximum nondivergent wind speed. On the other hand, the 2D model does not provide an entirely accurate picture of the low viscosity perturbation dynamics. The oscillation frequency of the 3D MUM is 0.8 times that of the 2D MUM, and the growth rate is 0.4 (0.6) times that predicted by the 2D model with K 2d ¼ 0 (1000 m 2 s -1 ). Greater facilitation of diabatic forcing (e v ¼ 1) at low viscosity (e k ¼ 0:0625) leads to much greater disparity between the 3D and 2D MUMs. The 3D MUM (Fig.  14c) exhibits complex vertical structure deep into the free troposphere, and the perturbation fields near the surface have more features in common with the low-frequency SMUM of the 2D model (Fig. 16a). Consistently, the oscillation frequency of the 3D MUM is 0.5 times that of the 2D MUM.

Comparison of linear instability theory to NS14
Reducing the general uncertainty of linear instability theory will require refinement of the physics parameterisations. Such refinement will require a comprehensive comparison of theory to state-of-the-art cloud resolving numerical simulations. While a comprehensive refinement effort is beyond the scope of this paper, a comparison of our linear model to the results of one of our earlier simulations is easy and worth reporting. The simulation considered for illustrative purposes corresponds to the three-dimensional moist experiment of NS14 distinguished from others by the following ratio of surface-exchange coefficients: C e =C d % 0:3. The experiment examined the evolution of a random perturbation of an initially axisymmetric category-2 hurricane in CM1. The disturbance followed an initial pattern of development similar to that found in all simulations of the study, including those with larger values of C e =C d and stronger vortices. Specifically, the perturbation spurred asymmetric wave growth energetically concentrated near the surface, and the wave growth engendered a ring of five welldefined mesovortices (Figs. 3 and 6 of NS14). The physics parameterisations utilised in the experiment differed from those described in Section 4 in several notable ways. To begin with, the microphysics parameterisation excluded ice. So as to keep the ratio of surface-exchange coefficients constant over the entire ocean, the drag coefficient was held fixed (along with C e ) at C d ¼ 0:005. Perhaps of greatest significance, K m=h h was an order of magnitude smaller in the vicinity of maximum wind speed. The reader may consult NS14 for further details.
For better compatibility with the model configuration of NS14, the physics parameterisations used presently in computing the linear instability modes differ somewhat from those used previously. Diabatic forcing is given by Equation (9b), butṽ b is calculated under the assumption of liquid-only condensate. The drag coefficient is simplified to C d ¼ 0:005e k . The variables K m h;sm and K m v;sm in Equation (16) are obtained as before, but from the simulation used to generate the basic state of the NS14 vortex. Typical values of both K m h;sm and K m v;sm are between 100 and 200 m 2 s À1 where the pertinent instability modes are concentrated. The lower limits of the eddy diffusivities are given by K m h;min ¼ 100 m 2 s À1 and K m v;min ¼ 40 m 2 s À1 . Figure 17 compares wave growth in the CM1 simulation to that predicted by linear instability theory. The squares show (a) the growth rates and (b) the oscillation frequencies of the primary Fourier components of the asymmetric radial velocity field (du uÀ u) in the simulation. The measurements are made by a straightforward procedure. To begin with, du is vertically averaged over the interval 0 < z < 1:0 km and expanded into a discrete Fourier series with respect to the azimuthal coordinate u. The wavenumber-n Fourier coefficient of the vertically averaged field is denoted du n ðr; tÞ. Time series of the amplitude and phase of du n (for all n between 1 and 8) are obtained from three probes placed 10-km apart on a radial line segment that is centred roughly at the radius of maximum wind. Each time series is taken over the interval 0 t 90 min. Data during the initial adjustment period and near the end of the interval (when incipient mesovortices take form) are generally discarded. The growth rate is obtained from an exponential curve fit to the amplitude data, whereas the oscillation frequency is obtained from a linear regression of the phase data (over 1 oscillation period). The plotted growth rates and oscillation frequencies correspond to their respective means among the three probe measurements; each error bar covers the full range of probe values. The preceding measurements are sensibly associated with the complex eigenfrequencies of MUMs provided that a single growing wave dominates du n . Such a condition appears to be satisfied quite well for 4 n 7. The greater error bars shown for n ¼ 3 and n ¼ 8 indicate that the probe signals are not as clean. Because the time series for n ¼ 1 and n ¼ 2 do not closely resemble those of a single growing wave, they are excluded from the plots.
The diamonds in Fig. 17 show the growth rates and oscillation frequencies of the 3D MUMs predicted by the linear model for the tropical cyclone simulated in NS14. The MUMs were first identified as perturbations dominating solutions of the initial value problem with e v ¼ e k ¼ 1 on a gridG2 with double the resolution of that used in NS14. The diamonds are centred on values of k R and k I obtained by recomputing the eigenmodes on a grid G4 with double the resolution ofG2. Further sensitivity to grid spacing was examined by repeating the computations on the original NS14 grid. Sensitivity to the parameterisations of diabatic forcing and small-scale turbulence were separately examined by reducing e v to 0 and e k to 0.0625 onG2 andG4. The error bar on each diamond covers the full range of values for k R or k I obtained from all configurations of the linear model; the smallness of the error bars indicates robust results. Note that the plotted eigenfrequencies are confined to n between 2 and 8, which correspond to eigenmodes energetically concentrated near the surface (not shown); a slower growing middle-tropospheric MUM associated with n ¼ 1 is excluded from present consideration. Figure 17 clearly demonstrates that the eigenfrequencies of the theoretical and simulated 3D MUMs are in good agreement where the latter are inferred from the cleanest monochromatic signals (4 n 7). On the other hand, both the growth rates and oscillation frequencies are smaller than those of the MUMs associated with an analogous 2D vortex (circles). The 2D vortex under consideration is modelled after the primary circulation of the NS14 tropical cyclone averaged in z over a layer of thickness d adjacent to the surface. The plotted 2D data points correspond to the means taken from 6 configurations in which d 2 f1; 2; 3g km and K 2d 2 f0; 10 3 g m 2 s -1 . As usual, the error bars extend from the minimum to maximum values of the data set for each n.
In this particular case study, the insensitivity of 3D linear instability theory to the degree of diabatic forcing allowed in the model is consistent with the concentration of modal wave activity inward of the eyewall cloud (NS14). Insensitivity to the reduction of e k seems reasonable given the short e-folding times of the MUMs (15-20 min) relative to the minimum applicable time scale for turbulent diffusion, s k minðl 2 h =K h ; l 2 v =K v Þ, in which l h=v is the horizontal/vertical lengthscale relevant to the mode and K h=v is the horizontal/vertical eddy diffusivity. Taking K h=v 200 m 2 s -1 and l h=v ! 10 3 m yields s k ! 83 min.

Conclusion
This paper has proposed a method to account for diabatic forcing and inhomogeneous eddy diffusivities in predicting and analysing the dominant instability modes of numerically simulated tropical cyclones. Excluding explicit moisture equations from the linearised model necessitated a partly intuitive parameterisation of the diabatic forcing S 0 h . The parameterisation considered herein set S 0 h proportional to w 0 with the modulating coefficient e vṽb o z h qb dependent on the local moist thermodynamic conditions of the basic state. A more general parameterisation scheme [Equation (10)] was presented for future consideration. The instability analysis was illustrated for a mature tropical cyclone representative of a category 4 hurricane.
The basic state was generated by an axisymmetric numerical simulation with two-moment cloud microphysics and typical settings for the parameterisation of subgrid turbulence. Initial consideration was given to linear systems having vertical and horizontal eddy diffusivities comparable to those regulating the basic state. With the diabatic forcing parameter e v set to a value between 0 and 1, perturbation growth was commonly dominated by a slowly growing n ¼ 2 eigenmode with deep structure but maximal intensity (KE n ) in the lower tropospheric region of the inner core. The complex eigenfrequency, spatial structure and energetics of the n ¼ 2 MUM were sensitive to variation of e v . Increasing e v from 0 to 0.9 gradually stabilised the mode. Further amplification of e v to 1 introduced a new MUM distinguished in part by having a larger growth rate than any of its predecessors, and by having a slightly positive buoyancy-related contribution to the production of integrated KE n .
Reducing the eddy diffusivities with e v fixed at either 0.5 or 1 generally changed the nature of the n ¼ 2 instability. For e v ¼ 0:5, the original MUM was ultimately replaced by a faster surface-concentrated instability mode whose growth of KE n involved a much smaller fractional contribution from the term directly linked to the secondary circulation. For e v ¼ 1, the original MUM was replaced by a faster instability mode whose growth of KE n distinctly involved a negative contribution from the buoyancy term and a relatively large positive contribution from the tendency associated with the vertical shear of the primary circulation.
Sensitivity of the foregoing analysis to the parameterisations of diabatic forcing and turbulent transport attests to the importance of details in predicting and understanding tropical cyclone instabilities. Improving the predictive skill of the linear model will require reducing the present degree of uncertainty in the aforementioned parameterisations. Refinements of S 0 h and D 0 a will come through a combination of theoretical advancements and testing of the linear model against perturbation growth found in state-of-the-art cloud resolving models.
An initial test of our linear model produced encouraging results. The instability analysis showed very good quantitative agreement with the perturbation growth that leads to mesovortex formation slightly inward of the eyewall cloud in a previously conducted CM1 simulation with relatively low diffusivity (NS14). Such agreement helped validate the dynamical core of the linear model. On the other hand, questions regarding the parameterisation schemes were left unresolved. The instability was theoretically too fast for reasonable variants of turbulent transport to have an appreciable effect on its early development. Moreover, the perturbation seemed largely detached from moist processes (NS14). Accordingly, the instability predicted by the linear model showed little sensitivity to switching e v between 0 and 1. the complex eigenfrequencies of the filtered and unfiltered instability modes hardly differ, regardless of the selected parameters regulating the strengths of diabatic forcing and turbulent transport.
Appendix C. Sensitivity of MUMs to the computational grid In Section 5, a MUM was provisionally equated to the eigenmode that dominates a perturbation within 1 day of initialising the linear model [Equation (20)] with random noise on a grid (G2) having double the resolution-in both r and z -of the CM1 grid used to generate the basic state (G1). The MUMs of G1 are readily found by a similar procedure. Figure C1 displays the complex eigenfrequencies associated with the MUMs of both G1 (blue) and G2 (green) for 0 < n < 8 when the linear model is parameterized with e v ¼ e k ¼ 1. Note that no discernible instabilities could be seen for n ¼ 0 or n ¼ 8. Also shown are solutions to the eigenproblem on a grid (G4, red) with quadruple the resolution of G1. Each G4 data point is obtained from an algorithm that seeks a solution for the complex eigenfrequency (k) closest to that of the G2 MUM.
The system under present consideration exhibits a somewhat complicated sensitivity to grid spacing. The n ¼ 1 mode is virtually invariant with increasing resolution. The oscillation frequency (k I ) of the n ¼ 2 mode is also robust, but the growth rate on G1 exceeds that on G4 by 29%. The values of k R (for n ¼ 2) on G2 and G4 are deemed closer to the continuum limit based on their modest 2% difference. All of the modes on G1 that are shown for 1 n 4 have maximal KE n near the surface. Whereas the properties of the modes with n ¼ 1 and n ¼ 2 are fairly insensitive to increasing resolution, the n ¼ 3 and 4 modes on G1 are fragile and superceded by middle tropospheric instabilities on G2 and G4. All of the dominant instabilities at higher wavenumbers have maximal KE n in the middle troposphere. It is notable that increasing the resolution for cases in which n exceeds 5 markedly accelerates the instabilities. Moreover, the resolution required to establish less than 10% uncertainty in k R at high-n appears to be greater than that of G4. f b , r 2 n W ¼ Z; r 2 n o rr þ r À1 o r Àn 2 =r 2 , and o rr o r o r . A formal solution for the wavefunction of w 0 consistent with regularity at the origin and u 0 ¼ 0 at r ¼ r B is W r ð Þ ¼ À 1 2n ð rB 0 drr r < r > n 1 À r > r B 2n " # Zr ð Þ; (D2b) in which r < (r > ) is the lesser (greater) ofr and r (Schecter et al., 2000). The second outer boundary condition o r ðv 0 =rÞ ¼ 0 combined with u 0 ¼ 0 amounts to Z ¼ 2r À1 dW=dr at r ¼ r B . Substituting (D2b) into both (D2a) and the preceding outer boundary condition eliminates W from the eigenproblem. Subsequent discretization of the radial coordinate yields a standard matrix eigenproblem of the form Mx ¼ kx, in which x is a vector containing the values of Z on each grid point. The 2D results of Section 5.3 and 6 correspond to solutions of the preceding eigenproblem; the outer boundary condition on v 0 is obviated for computations in which K 2d ¼ 0. Selected results were successfully crosschecked against independent solutions of the 3D model set up with a thin barotropic vortex sandwiched in between rigid free-slip walls at z ¼ 0 and 0.5 km.