## Abstract

This study proposes that secondary eyewall formation (SEF) of tropical cyclones (TCs) can be attributed to an instability of flow in the free atmosphere coupled with Ekman pumping. Unstable solutions of a 1.5-layer shallow-water system are obtained under fast–wind speed conditions in the free atmosphere. The instability condition derived in the linear model indicates the importance of the ratio of angular velocity to vorticity, and the condition is more likely to be satisfied when the ratio is large and its radial gradient is positive. Thus, fast angular velocity, low absolute vertical vorticity, small negative radial gradient of angular velocity, and large negative gradient of vertical vorticity are favorable. Eigenvalue analyses are performed over a wide range of parameters using a vorticity profile with an infinitesimal secondary maximum. The growth rate increases with vorticity outside the radius of maximum wind (RMW), the radius of the secondary vorticity maximum, its magnitude, and the Rossby number defined by maximum tangential velocity, the RMW, and the Coriolis parameter. Furthermore, the growth rate is positive only between 2 and 7 times the RMW, and it is negative close to or far outside the RMW. These features are consistent with previous studies on SEF. A dimensionless quantity obtained from the unstable condition in the linear theory is applied to SEF events simulated by two different full-physics numerical models; increases several hours before a secondary peak of tangential velocity forms, suggesting that the initial process of SEF can be attributed to the proposed mechanism.

## 1. Introduction

Strong tropical cyclones (TCs) often experience secondary eyewall formation (SEF; Willoughby et al. 1982), in which a new eyewall forms outside the original eyewall. It is usually followed by an eyewall replacement cycle, and hence, predicting SEF is important for accurate TC forecasts. However, there is no widely accepted theory for SEF despite the significant effort by many researches, especially in the last decade. We here present a new theory focusing on flow interaction between the boundary layer and the free atmosphere.

SEF often occurs at a radius of 2–4 times the radius of maximum wind (RMW) in recently intensified TCs (Willoughby et al. 1982; Black and Willoughby 1992; Dodge et al. 1999; Houze et al. 2007; Kossin and Sitkowski 2009; Kuo et al. 2009; Wu et al. 2016). The key parameters and processes for SEF proposed by previous studies are inertial instability (Willoughby et al. 1982), terrain effect (Hawkins 1983), ice microphysics (Willoughby et al. 1984), inflow surge (Molinari and Skubis 1985; Molinari and Vollaro 1989), vortex Rossby waves (Montgomery and Kallenbach 1997; Terwey and Montgomery 2003), potential-vorticity intrusion from upper layer (Nong and Emanuel 2003), axisymmetrization of vorticity disturbance (Kuo et al. 2004, 2008), the rectification of eddies in the *β* skirt outside the primary eyewall (Terwey and Montgomery 2008), supergradient wind (Huang et al. 2012; Wu et al. 2012), wind field expansion (Rozoff et al. 2012; Sun et al. 2013), radial vorticity gradient outside the RMW where relative vorticity is low (Kepert 2013), top-down process associated with interaction with upper-level trough (Leroux et al. 2013), rainband heating (Zhu and Zhu 2014; Zhang et al. 2016), and upper-level outflow interactions (Dai et al. 2017).

Kepert (2013) showed the importance of the radial gradient of relative vorticity in a low-vorticity environment using both linear and nonlinear boundary layer models. According to the Ekman theory for a rotating fluid, Ekman pumping velocity is proportional to the radial vorticity gradient and inversely proportional to the square of the absolute vorticity itself. Kepert (2013) and Kepert and Nolan (2014) proposed a positive feedback that can lead to SEF: when the radial gradient of vertical relative vorticity is steep, Ekman pumping velocity is large, the upward velocity is enhanced by convection and causes further convergence, and the vorticity gradient steepens. However, the previous studies consider vorticity profiles *after* SEF and secondary wind maximum, in which a secondary peak of tangential velocity is already present outside the RMW. From ensembles of numerical simulations, Zhang et al. (2016) showed that the diabatic heating outside the RMW in the upper layer proceeds downward, changing the radial vorticity profile, and the secondary eyewall develops through the interaction between convection and boundary layer processes. Although the interaction might work once the vorticity gradient is steep, the proposed mechanism is incomplete because it offers no explanation of conditions under which outer rainbands produce such a strong vorticity gradient (otherwise, rainbands would always cause SEF). Therefore, the underlying mechanism of SEF, especially for the initial process (i.e., when and where SEF occurs), is unclear.

This study proposes a linear instability via an interaction between Ekman pumping and the flow in the free atmosphere, which is not associated with moist convection, unlike the previous theories, and hypothesizes that very early processes of SEF are attributed to this instability. Once the initial processes of SEF proposed in this study create a local maximum of vertical velocity or tangential velocity outside the primary eyewall, the heating–convergence feedback might work to develop the secondary eyewall, followed by an eyewall replacement cycle (Kepert 2013; Kepert and Nolan 2014).

The rest of the paper is organized as follows. We develop a linear theory for a 1.5-layer shallow-water system in Cartesian coordinates and eigenvalue analyses and numerical integrations for a nonlinear equation set are conducted in section 2. The theory is extended to cylindrical coordinates and a condition for the instability to occur is derived and linear stability analyses of an axisymmetric 1.5-layer shallow-water system are performed in section 3. The unstable condition is tested by applying it to simulation results by two nonhydrostatic full-physics models in section 4. The results are discussed by comparing them with previous studies in section 5. This study is concluded in section 6.

## 2. Extension of the Ekman theory in Cartesian coordinates to high–wind speed conditions

### a. Classic Ekman theory

Ekman pumping is a vertical motion at the top of the planetary boundary layer induced by surface friction (Ekman 1905; Holton 1992). A cyclonic flow with Ekman pumping transports angular momentum in the interior layer outward, and the flow decays much more quickly than would occur from viscous diffusion (vortex spindown; Holton 1992; Eliassen and Lystad 1977).^{1} This can be explained by considering a one-dimensional 1.5-layer shallow-water system on an *f* plane with a semislip lower-boundary condition: a free atmosphere plus a boundary layer (cf. Fig. 1, but velocity in the free atmosphere is weak and its horizontal variations are slow). We can obtain an analytical solution for the evolution of perturbations of relative vorticity, and the solution indicates that the vorticity perturbation decays exponentially. A key relationship in the classic Ekman theory is that perturbation vertical velocity at the top of boundary layer is proportional to the horizontal gradient of the *y*-velocity perturbation :

Details to obtain this relationship in the 1.5-layer shallow-water system can be found in appendix A.

The Ekman damping process for velocity perturbations can be interpreted as follows: once a positive is added at the top of the boundary layer (left-hand side of Fig. 2), it produces positive *x*-velocity perturbation in the free atmosphere for mass conservation at larger *x* and negative at smaller *x*, which decelerates and accelerates, respectively, *y*-velocity perturbation by transporting basic-state absolute vorticity. Thus, the horizontal gradient of perturbation *y* velocity is negative around the peak. Assuming that in the boundary layer is changed in a similar way, by mass divergence, tends to be reduced, and hence, the initially added is damped. This is because of (1).

The classic Ekman theory is widely accepted and has provided many insights into geophysical fluid dynamics. For example, the effects of Ekman pumping on Rossby waves are modeled by adding vertical velocity at the top of the boundary layer that is proportional to vertical relative vorticity (e.g., James 1995; Vallis 2017). Although the theory can capture fundamental processes for various phenomena, the assumptions may limit its application. In particular, it is assumed that the background flow is not fast and its spatial gradient (relative vorticity) is not large. Hence, the theory may not be applicable to some circumstances such as jets or tropical cyclones in which the horizontal velocity is fast in the lower free atmosphere.

### b. Extension of Ekman theory to fast velocity conditions

We here hypothesize that the relationship in (1) can possibly be reversed when we include the advection in the boundary layer. If this were to occur, would be negatively proportional to the gradient of ; that is, (right-hand side of Fig. 2). According to the mechanism described above, a negative gradient of is produced around a positive peak of because of mass conservation and horizontal transport of absolute vorticity, which is invariant even in fast flow conditions. In this case, however, the negative velocity gradient tends to produce *positive * around the initial peak. It means that is amplified and thus perturbations will grow exponentially.

To demonstrate the hypothesized instability, equations for the 1.5-layer shallow-water system on an *f* plane [(A7)–(A11)] are linearized. The basic state satisfies the geostrophic balance in the free atmosphere. In addition, we here assume that the basic-state *y* velocity is strong so that the advection of by basic state *x* velocity in the boundary layer is nonnegligible. A diagnostic equation for is obtained from the momentum equation for the *y* direction [(A11)] as

The linearized equations for mass and momentum conservation in the free atmosphere are

and those in the boundary layer are

We have neglected correlation terms between perturbations (e.g., ). The fundamental differences from the classic Ekman theory in the present 1.5-layer system are the advection of by [the second term in the lhs of (7)] and the nonuniform basic-state *y*-velocity [, and hence, ], the first term in the lhs of (7) and the rhs of (5).

Substituting (7) into (6) yields an equation for perturbation vertical velocity at the top of boundary layer:

where the coefficients consisting of basic-state quantities are

In the classical theory, and , because the variation of is small compared to *f*, which results in as . The terms appearing in the rhs of (8) in addition to that in the classic theory are due to the presence of basic-state flow in the boundary layer and the spatial variability of . As a result, the coefficient defined by (10) is no longer positive definite and can possibly be negative. The hypothesis indicates that the instability occurs when .

A series of numerical calculations for the equations for the free atmosphere [(3)–(5)] and vertical velocity [(8)] have been performed to test the hypothesized instability by artificially setting the coefficients , , and ( appendix C). It was confirmed that the perturbation of vertical velocity grows only when , which is consistent with the proposed theory, and all the other cases such as or show decaying solutions.

### c. The condition for instability

We here analytically derive the condition for the instability (i.e., ) in the 1.5-layer system. With some algebra and using the fact that results in an alternative form of the condition

If we consider cases with , the most favorable case to satisfy the inequality (12) is when both the first- and second-order derivatives of velocity and are negative (see appendix B). In other words, the conditions correspond to when the relative vorticity is negative and the spatial gradient of vorticity is also negative. This flow field is expected to appear outside the peak of a jet (at larger *x* than the peak velocity) where the velocity gradient is negative and the second-order derivative is also negative up to the inflection point. The unstable condition (12) can be rewritten as

where

Thus, is a function of the spatial distribution of and the lhs of (13), or equivalently, the sign of depends on the velocity distribution as well as Ro. To satisfy the unstable condition (13), must be positive, and both Ro and should be large. Specifically, a larger magnitude of and and smaller magnitude of are more favorable. As expected above and will be shown later, these conditions are satisfied in a jetlike flow.

### d. Eigenvalue analysis

We conducted numerical eigenvalue analyses for the 1.5-layer system to seek unstable modes produced by the proposed mechanism. The linearized equations for the free atmosphere [(3)–(5)] and those for the boundary layer [(6) and (7)] can be combined into a single equation for :

where

In the classic theory in which the basic-state velocity is not strong and varies slowly in space, , and . The velocity perturbation that is a function of space *x* and time *t* is decomposed as , where *σ* is the growth rate. When Re, the disturbance grows exponentially. Substituting this into (15) yields the eigenvalue equation for the growth rate:

where and are matrices for the lhs and rhs of (15), respectively, and is a vector form of .

The eigenvalue equation [(16)] was discretized on a staggered grid system and eigenvalues and eigenvectors were computed with MATLAB. As in the test calculations in appendix C, the spatial derivatives in (15) were discretized by the fourth-order central difference scheme, and periodic boundary conditions were used. We have tested the second-order difference scheme and found that the results are insensitive to the discretization scheme. The nondimensional domain width was 8 with 201 grid points and uniform grid spacing of 0.04. The initial velocity field was a jetlike profile with a peak at and a spatial-scale , where *L* is a length scale. The velocity was computed from *x* = 0 to 4, and then the value at {i.e., } was subtracted to make the velocity zero at . To satisfy the periodic boundary condition for , one other velocity profile that is symmetric to that from *x* = 0 to 4 with the same peak magnitude and width but opposite sign was added from *x* = to 8; that is, (Sugimoto et al. 2007). The basic-state depth was in geostrophic balance with .

Figures 3a and 3b show spatial profiles of , vertical absolute vorticity , and in (13) for a representative case, in which Ro = 0.70, Fr = 0.10, and , corresponding approximately to *V* = 20 m s^{−1}, , and . The velocity has a peak value at and gradually decays in each side of peak, and the vorticity has positive and negative peaks at *x* = 1.52 and 2.48. increases with *x* up to the maximum at and sharply decreases from there to . The unstable condition, (13), is met in this case between *x* = 1.88 and 2.36. As discussed in the previous section, the condition is more likely satisfied to the right of the peak of the jet where both the first-order and second-order derivatives of the basic-state velocity are negative. The range satisfying the unstable condition is highlighted in Fig. 3.

The eigenmodes with the maximum eigenvalue are shown in Figs. 3c and 3d. In this case, the real part of the maximum eigenvalue (i.e., growth rate) is 18.73 and corresponds to an *e*-folding time of 1.33 × 10^{4} s (≈3.7 h). The vertical velocity was calculated using (8) from the eigenmode of , and then was diagnosed from (3) to (5). The plot of has a sharp positive peak at . The eigenvector of also has sharp peaks around the peak. Finally, is largely negative at the positive peak of and positive at the negative peak. These profiles of velocities are basically consistent with the proposed theory.

We conducted a series of eigenvalue analyses for wide ranges of Ro and the spatial parameter that controls the width of the velocity profile and hence changes the spatial derivatives of . Specifically, monotonically increases with . According to the theoretical condition for the instability in (13), Ro and are the key parameters to determine the stability. The parameter-sweep ranges are listed in Table 1. Meanwhile, Fr and were fixed as 0.1 and 1.0, respectively, in all the calculations, because they do not affect the stability criteria.

Figure 4 displays a phase diagram in space of two parameters: Ro and the domain maximum of . The neutral curve, , is shown by the black line. The circles indicate experiments in which the real part of the maximum eigenvalue (growth rate) is positive and the domain-maximum is greater than 1.0. The × marks indicate experiments that do not satisfy both the criteria. Note that the gray circles represent experiments in which the real part of the eigenvalue is positive but the minimum value of basic-state absolute vorticity in the domain is negative, which is generally unphysical. The eigenvalue analyses show that the vertical velocity amplifies when the theoretical condition (13) is met. The stable and unstable regimes are clearly separated by the neutral curve . Beyond the neutral curve, the solutions are mostly unstable, and growing modes are obtained. It should be noted that the negative absolute vorticity regime also appears beyond the neutral curve. The unstable regime with the minimum vorticity greater than 0 (black circles) is sandwiched in a narrow range between the neutral curve and the negative vorticity regime. It is indicated that the instability will occur when the unstable condition is satisfied but under limited environmental conditions.

### e. Numerical integrations of nonlinear equation set

A series of numerical integrations for the nondimensional form of the nonlinear equation set (A7)–(A11) were performed with variations of the two parameters, Ro and . We added a tendency term to the momentum equation for the *x*-direction in the free atmosphere (A8) to allow ageostrophic components. The dimensionless form of *u* equation used here is

where and is determined by using the parameter values shown in Table 2. The mass conservation and *y*-momentum equations in the boundary layer [(A10) and (A11)] were combined into an equation for *w*.

The equations were solved numerically using a finite difference method. As in the eigenvalue analysis, the quantities were allocated on a staggered grid system. The time and space derivatives in the equations were discretized by the third-order Runge–Kutta method and by the fourth-order central difference scheme, respectively. At the first time step, *w* was estimated from *υ*, and then the equations for the free atmosphere were stepped forward. The integration was conducted up to , dimensionally equivalent to 2.3 days, with a time step of . The domain size, grid spacing, and boundary condition were the same as in the eigenvalue analysis. The profile of *υ* used in the eigenvalue analysis (Fig. 3a) was also applied to the initial field, which was balanced with the fluid depth. The experimental parameters and their range were also identical to the eigenvalue analysis.

Figure 5 shows spatial distributions of velocities at *t* = 0.02 in the control case, in which Ro = 0.80 and . The vertical velocity has a positive peak at *x* = 2.40, which is the same location as the eigenvalue analysis (cf. Fig. 3c) and a negative peak near the positive one. There is a sharp peak in the *u* distribution around the positive peak of *w*, whereas *u* is negative at the both sides of the peak. The negative *u* transports absolute vertical vorticity around the peak of *υ* where the vorticity is not weak compared with the negative *u* region outside of the peak. As a result, *υ* is amplified around the peak of initial distribution. The velocity distributions in the nonlinear integration are qualitatively consistent with the eigenvalue analysis and the present hypothesis.

Figure 6a displays a space–time cross section of *υ* in the control case and time series of domain-maximum *w*. The maximum *w* amplifies with time. The maximum *w* keeps growing until the end of the simulation, and the locations of velocity peaks seen in Fig. 5 are nearly fixed in space. The spatial profiles of velocity during the integration (Fig. 5) strongly suggest that the proposed feedback works to amplify *w*. Since no explicit diffusion is incorporated into the discretized equation set, the linearly unstable modes continue to grow.

Figure 7 depicts a phase diagram in Ro and the domain-maximum . The experiments in which the domain-maximum *w* amplifies and exceeds 5 m s^{−1} are observed in the unstable regime derived from the linear theory, Ro. The amplifying regime is clearly separated by the neutral curve from the regime in which the maximum *υ* decays after the initiation of integration (× marks in the figure). The theoretically derived unstable regime also includes a regime in which the domain-minimum absolute vertical vorticity is negative at the initial time. These features are quite similar to those obtained in the eigenvalue analyses (cf. Fig. 4). It should be noted that unlike the eigenvalue analyses, *w* does not amplify in several cases close to the neutral curve. In these cases, the maximum *w* decays rather than amplifies. The failure of unstable modes to grow when the parameters are close to the neutral curve may be due to nonlinear terms.

## 3. Theory in cylindrical coordinates

We here extend the theory developed in section 2 to cylindrical coordinates. The dimensionless governing equations are obtained in appendix D. The linearized mass conservation and radial and tangential momentum equations for free atmosphere are given as

and linearized equations for boundary layer are

where and is basic-state absolute vertical vorticity. The basic-state radial velocity in the boundary layer is given by the azimuthal momentum equation [(D12)]:

The perturbation vertical velocity at the top of boundary layer is obtained by combining (21) and (22) with the use of (23):

where

and

are the ratio of angular velocity to absolute vorticity. When the wind speed is not so fast in the free atmosphere that the advection term by the basic-state velocity is negligible and *υ* does not vary in space, , and . This converges to the classic theory where is proportional to the perturbation vorticity.

The necessary condition for the proposed instability is . Using (23), this condition can be written as

The condition is more likely to be satisfied as the factor in parentheses increases. Since we only consider cases in which and , and hence , the second term in the parentheses is always positive, whereas the first term can be either positive or negative. Thus, the terms in the parentheses are positive once is larger than 0, which can be written as

Since the angular velocity is larger than 0, it is rewritten using (28) as

It indicates that, for the instability, the radial gradient of vorticity should be more negative than the right-hand side, because rarely takes positive values outside the RMW. Furthermore, it is more favorable if the absolute vorticity is small, angular velocity is large, magnitude of radial gradient of angular velocity is small, and magnitude of radial gradient of absolute vorticity is large. The importance of small absolute vorticity and large negative gradient of vorticity is consistent with Kepert (2013).

Equation (24) for is a function of , and hence, the rhs of condition (30) is determined by the radial profile of and Ro. To separate Ro and terms associated with, the condition, , is written as

where, in a vortex,

The value of tends to increase with the larger angular velocity and absolute value of second-order derivative of velocity and smaller first-order derivative of velocity , which is similar to the result obtained above.

### Eigenvalue analysis in a cylindrical coordinate

To investigate the stability of the linear system developed in the previous section, linear stability analyses are conducted by solving an eigenmode equation for . The equations for the free atmosphere [(18)–(20)] and the equation for vertical velocity [(24)] are combined into a single equation for :

where

Substituting a waveform solution, , where *σ* is the growth rate, into (35) yields the following eigenvalue equation:

where and are the matrices in the lhs and rhs of (35) and is the vector form of .

The basic-state quantities were produced from a radial profile of vorticity in Kepert (2013) that used a smoothing function from Willoughby et al. (2006):

where is a weight function (Willoughby et al. 2006) and is the same as in Kepert (2013). The filter width was 15 km. The vorticities in the regions A, B, and C are

In the analyses, was given by Ro, where was fixed as 5 × 10^{−5} s^{−1}. A radial profile of dimensional tangential velocity was obtained by radial integration of (37). Then was normalized by the maximum value. Fluid depth and radial velocity in the boundary layer were obtained by gradient wind balance and (23).

The radial dimension was discretized on a staggered grid system with uniform grid spacing of 0.04. The domain size was 80.0 (i.e., 80 times the RMW) with 2001 grid points. The radial derivatives were discretized by a fourth-order central difference scheme. At the boundaries (*r* = 0 and 80), .

Figures 8a and 8b show radial profiles of basic-state and in the control case (Ro = 50, = 4.0, , , ); has a peak at and decays monotonically outside the RMW. The vorticity is constant inside the RMW and rapidly decreases around the RMW. There is a small local decrease of vorticity at *r* = 4.0. However, this jump is more than 10 times smaller than that in Kepert (2013). Hence, unlike most of the cases in Kepert (2013), has no clear secondary peak around *r* = 4.0. The ratio of angular velocity to absolute vorticity has a peak outside the RMW and a small peak around *r* = 4.5. also has a sharp peak immediately outside the RMW and a peak around *r* = 4.0, where the radial gradient of is large. The unstable condition (33) is satisfied in in this case, which is highlighted in the panels.

The eigenmodes of , , and are shown in Figs. 8c and 8d. In the cylindrical coordinate model with the hurricane-like wind profile, Fig. 8b also shows a large value of around the RMW. Not surprisingly, a very unstable mode is also excited in that region. To isolate the instability in the SEF region, we select the unstable mode that has the fastest growth rate of but also has significant amplitude of at *r* > 1.5 × RMW. The , , and profiles of this mode are shown in Figs. 8c and 8d; the instability around the RMW is also evident. The nondimensional growth rate is 52.4, and the corresponding dimensional *e*-folding time is . The value of is large and positive around *r* = 4.0 where has a peak; is largely negative (positive) outside (inside) the radius of *r* = 4.0, which is opposite to in phase.

We examine the sensitivity of the solution to Ro and vortex parameters , and . Each parameter out of the five is changed to examine the sensitivity, while the other four are fixed as those in the control case (shown above). Table 3 shows the range of parameters. The minimum and maximum values of Ro are 20.0 and 87.5, which approximately correspond to *υ*_{m} = 25.0 and 109.4 m s^{−1}, provided that *r*_{m} = 25 km and . The ratio of the radius of the second “vorticity jump” between the region B and C to the RMW varies from 2 to 11. The coefficient controlling vorticity outside the RMW is changed from 0.005 to 0.05. The magnitude of secondary vorticity jump is between 0.1 and 1.0. The radial decay parameter *α* is changed from 0.1 to 1.0.

Figure 9 shows the growth rate as a function of experimental parameters. The growth rate increases with Ro, , , and and decreases with . Thus, the growth rate is larger with stronger TC intensity, smaller RMW, and lower Coriolis parameter. The secondary peak of velocity tends to form more outside. Interestingly, the growth rate is positive only between *r* = 2.5 and 7.0: the disturbance cannot grow if it is too close to the RMW or too far from the RMW. The growth rate is large with large , corresponding to large vorticity and hence fast velocity outside the RMW. When the vorticity outside the RMW is not strong enough , unstable solutions are not obtained. Both the dependence on and indicate that the growth rate is larger when the radial vorticity gradient at the secondary vorticity jump is steeper. The large vorticity and fast velocity is consistent with the wind field discussion that has been discussed in the previous studies.

To summarize the series of parameter sensitivity experiments, the relationship between the growth rate and divided by is shown in Fig. 10. The growth rate monotonically increases with . Whereas no analytical formula is obtained for *σ*, the figure indicates is a good scaling factor for the growth rate in the system; , and the first term is written as

Since and are positive and is never positive, a large negative radial gradient of is required for positive . Therefore, larger , smaller , smaller , and larger is more favorable. These are consistent with the analytically derived unstable condition in the previous section.

## 4. SEF simulated in full-physics models

We examine SEF simulated by two nonhydrostatic, full-physics numerical models. The unstable condition (33) was applied to the simulation results to test the hypothesis that the proposed theory accounts for the initial process of SEF. To calculate the unstable condition, (33) is dimensionalized as

Unlike the unstable condition for the linear model (33), Ro does not appear in the inequality as it is eliminated in the dimensionalization.

### a. Axisymmetric model

A long-term simulation was performed using an axisymmetric version of Cloud Model 1 (CM1; version 17; Bryan and Fritsch 2002), which solves fully compressible, nonhydrostatic equations. The experimental setting was basically the same as Hakim (2013) in which a number of eyewall replacement cycle events were simulated during his long-term simulation for 500 days. The differences from Hakim’s (2013) simulation are listed below. The initial mass, temperature, and water vapor mixing ratio was obtained from Jordan’s (1958) tropical mean sounding for hurricane season. The Coriolis parameter was 5 × 10^{−5} s^{−1}, and the sea surface temperature was 301.15 K, both of which were fixed during integration. The radiative effects were incorporated by the NASA Goddard scheme, which was called every 5 min. The integration period was 60 days (=1440 h) with an output interval of 3 h for prognostic quantities. The grid spacing for the radial direction was 2 km, whereas the spacing for the vertical direction was stretched from the bottom to top. The minimum grid spacing for the vertical direction was 0.25 km at the lowest model layer, stretching to 1.0 km near the top.

Figurea 11a and 11b show time series of TC intensity defined as the maximum tangential velocity at *z* = 2 km and its radius. The TC achieves maximum intensity around *t* = 340 h after an intensification phase, decreases to about 60 m s^{−1}, and then is approximately steady afterward, which is qualitatively consistent with the long-term simulations in previous studies (Hakim 2011, 2013). We shall call the period after *t* = 500 h a quasi-steady state. During the quasi-steady state, rapid changes in intensity and RMW are observed several times. As pointed out by Hakim (2011), the intensity oscillation appears to be due to the eyewall replacement cycle associated with SEF.

Figure 11c depicts a radius–height cross section of tangential and radial velocities and mixing ratio of hydrometeors averaged during the quasi-steady state. Both the velocity and cloud fields are consistent with the observational studies (e.g., Hawkins and Imbembo 1976): the tangential velocity has a peak at the top of the boundary layer (~1.5 km) and decays in both the radial and vertical directions, the inflow and outflow regions appear at the surface and around the 15-km level, and the hydrometeor mixing ratio is large around the RMW from the 1.5- to 15-km altitudes.

Figure 11d depicts a radius–time cross section of tangential velocity at *z* = 2 km. The RMW is about 60 km during the quasi-steady state as seen in Fig. 11b. Interestingly, a strong velocity region periodically forms about every 80 h, and the region proceeds inward after formation. The mixing ratio of hydrometeors and vorticity fields shows similar behavior to the tangential velocity (not shown). It is suggested that this cyclic process is the eyewall replacement cycle and formation of strong tangential velocity outside the RMW is SEF.

The SEF events were extracted from the simulation by applying the following conditions: the RMW is greater than that at the previous and next output time, the increase from the previous step is greater than 10 km, and the new RMW is more than 20 km beyond its average value during the quasi-steady state (~60 km in this case). Applying the conditions, six SEF events were detected in the quasi-steady period, and the six events were combined to make composites.

Figure 12a displays a radius–time cross section of tangential velocity at *z* = 2 km subtracted from the temporal average from to +40 h and , which are composited from the detected SEF events. The value of is calculated from at *z* = 2 km at each output time after taking 1–2–1 average for for 20 times in the radial direction. A secondary peak of tangential velocity forms at *r* = 163 km, which is about 100 km outside the RMW, 20 h before the time of SEF. The peak amplifies as it propagates inward. The becomes significant 3 h before the secondary peak forms, and the unstable condition is satisfied; is large around the RMW for the entire time frame, which is obviously due to the large radial gradient of vertical vorticity.

Figure 12b shows the time series of , , and , which are radially averaged from *r* = 135 to 160 km, where the secondary peak of tangential velocity forms; clearly has a large value at , which is immediately before the formation of the secondary peak. Since is large and negative at this time, the increase in results from the enhanced radial vorticity gradient. It suggests that the proposed linear instability works to amplify the tangential velocity to form a secondary peak outside the RMW.

### b. Three-dimensional model

We conducted a numerical simulation using a three-dimensional full-physics numerical model, Weather Research and Forecasting (WRF) Model, version 3.4.1 (Skamarock et al. 2008). The experimental setting is described in appendix E.

Figures 13a and 13b show time series of TC intensity and RMW of the simulated TC. The TC experiences two different intensification phases and becomes quasi steady after *t* = 45 h. The intensity weakens from *t* = 65 to 76 h and reintensifies afterward. The RMW gradually decreases, jumps outward at *t* = 76 h, and then decreases again. Structural evolution of the simulated TC (shown below) indicates that the sudden changes in TC intensity and RMW are associated with SEF and the eyewall replacement process. We define *t* = 76 h as the time of SEF.

Figure 13c depicts a radius–height cross section of the azimuthally averaged hydrometeor mixing ratio, tangential velocity, and radial velocity, which are temporally averaged from *t* = 48 to 54 h. Figure 13d depicts a radius–time cross section of azimuthally averaged tangential velocity at *z* = 2 km and diabatic heating rate that is averaged vertically from *z* = 2 to 10 km. The tangential velocity intensifies with decreasing RMW to *t* = 45 h, and then the RMW is approximately constant. A notable point is that the tangential velocity field expands with time after the TC achieves the maximum intensity. The temporal changes in structure and intensify during the SEF event are consistent with observations.

Figure 14a displays a radius–time cross section of the deviation of tangential velocity from the average between *t* = 48 and 96 h, diabatic heating averaged vertically from *z* = 1 to 10 km, and . The velocity field is smoothed by the 1–2–1 filter for 20 times before is calculated. Only the contours of 1.0 and 2.0 are shown for . The velocity deviation increases (i.e., a secondary peak of tangential velocity forms) around *r* = 125 km from *t* = −17 h (*t* = 59 h of simulation time). The secondary peak intensifies while proceeding inward, which is accompanied with strong diabatic heating. The rapid inward motion of the secondary peak is consistent with Kepert (2017). It suggests that the convection–convergence feedback plays a role in intensifying the secondary peak of tangential velocity. The time of secondary-peak formation is approximately the same as the increase in diabatic heating in the outer region. The point is that increases from , which is 1 h before the formation of a secondary peak.

Figure 14b shows a times series of , , and radial gradient of vorticity, which are averaged between *r* = 115 and 135 km where the secondary peak forms; is large enough to satisfy the unstable condition from 19 to 15 h before the time of SEF and is otherwise smaller than that satisfying the condition. The increase in results from enhanced and the radial gradient of vertical vorticity. The velocity field is expanding in this period (cf. Fig. 13d), indicating that the angular velocity is large, which is favorable to increase . In conclusion, a secondary peak of tangential velocity forms 17 h before the time of SEF, and is above the criteria of instability 18 h before the time of SEF.

## 5. Discussion

The theory indicates that a small gradient of vorticity can produce large such that the unstable condition for the linear instability is satisfied. Both Kepert (2013) and the present study show that vertical velocity can be strengthened at a radius where the radial gradient of vorticity is strong. However, the initial radial gradient of vorticity considered in Kepert (2013) is more than 10 times stronger than that in the present study. This study suggests that the vertical velocity can be amplified by the proposed instability under a much smaller gradient of vorticity than used by Kepert (2013) if the vorticity gradient increases at a radius where is large. Once the vertical velocity and radial vorticity gradient increase up to the magnitude of Kepert (2013), the outer eyewall can intensify via the interaction between the convection and low-level convergence, as Kepert proposed.

The initial amplification of the radial gradient of vorticity to satisfy the unstable condition might be caused by some processes such as diabatic heating in rainbands (Rozoff et al. 2012; Zhang et al. 2016). The present 3D simulation shows little convective organization before the secondary peak of tangential velocity forms (Fig. 14a). Rather, many isolated convective cells are present, but the accumulated effect of these can produce the needed vorticity gradient. Once a secondary peak of tangential velocity forms, a symmetric convective region is organized around that radius. Subsequently, the secondary peak and organized convection may amplify together through the feedback.

Figure 15 illustrates how the proposed feedback mechanism works, with a comparison to the classic Ekman feedback (and damping). Starting from perturbation of vertical velocity with a positive peak, the radial and tangential velocities in the free atmosphere will be produced (cf. Fig. 2). This change will result in faster velocity in the *x*-direction in the boundary layer to the left of the peak of the initial perturbation of vertical velocity and slower velocity to the right. Thus, the responding flow in the boundary layer tends to cause divergence and decrease *w*. On the other hand, when the unstable condition is satisfied, the change in flow fields in the free atmosphere rather produces slower velocity to the left and faster to the right, which increases in *w*. This would happen if the decrease in is significant because . Furthermore, it is more likely to happen when , and if , it will not occur.

## 6. Conclusions

We demonstrate an instability via interaction between the free atmosphere and the Ekman layer based on the classic Ekman pumping theory and hypothesize that the very early process of SEF in TCs can be attributed to this instability. The classic Ekman theory indicates that vertical motion at the top of the boundary layer is proportional to the spatial gradient of velocity (relative vorticity), and this results in a decay of velocity perturbations in the free atmosphere. Here, it was hypothesized that once the Ekman pumping velocity is *negatively* proportional to the velocity gradient under some special conditions, disturbances can grow exponentially.

The instability happens when velocity is fast and changes rapidly in the horizontal direction. Using a 1.5-layer shallow-water model on the *f* plane [(3)–(7)], we analytically derived an unstable condition (13), in which the proportionality coefficient of relative vorticity perturbation in the equation is negative; that is, . In particular, the unstable condition is more likely to be satisfied for larger Ro, faster basic-state velocity , smaller magnitude of negative first-order derivative of , and larger magnitude of second-order derivative of provided that is greater than 0.

Eigenvalue analyses were performed for the linearized equation set. A positive eigenvalue (growth rate) was obtained, and the eigenvector for the maximum eigenvalue showed a velocity field consistent with the hypothesis. In the parameter space of Ro and that controls the unstable condition, both stable and unstable regimes were obtained. The theoretical criteria of instability explains well the border of the unstable regime. Numerical integrations of a nonlinear equation set showed consistent results with the eigenvalue analyses.

We extended the theory to cylindrical coordinates. By linearizing the nondimensionalized equations, it was shown that the instability occurs where Ro is greater than 1. The unstable condition is more likely satisfied for larger Ro, , and . Eigenvalue analyses of the linearized 1.5-layer system were performed for varying parameters for the basic-state vortex. The eigenmodes of velocity were consistent with the proposed theory.

The hypothesis was tested by computing the dimensional values of in two different full-physics model simulations. We performed a long-term simulation using the axisymmetric version of CM1, and several cases of SEF and eyewall replacement cycle were successfully simulated. By extracting SEF events and compositing them, it was shown that a secondary peak of tangential velocity formed around 2.5–3.0 times the RMW 22.5 h before SEF, which was defined as the time when the tangential velocity at the secondary peak is stronger than the primary peak. The value of increased around this radius immediately before the secondary peak formed. The enhancement of resulted from an increasing radial gradient of vorticity with large .

Another SEF was simulated using WRF in an idealized environment with nonzero shear. A secondary peak of tangential velocity formed around 5.0–5.5 times the RMW and intensified while proceeding inward, while accompanied by strong convective heating. The local value of increased several hours before the formation of the secondary wind maximum.

The two modeling results suggest that the initial process of SEF, the formation of a secondary peak of a tangential velocity, may be attributed to the instability, and thus derived from the linear theory can potentially be used as a predictor of SEF. Once the secondary peak is amplified by the mechanism, the convection–convergence feedback proposed by Kepert (2013) works to develop the secondary eyewall. Whereas the simulations were conducted by two different models, further verification would be needed using especially observational data as well as other real-case simulations.

## Acknowledgments

Y. Miyamoto was supported by JSPS Scientific Research 26–358 for the JSPS fellowship program for overseas researchers. D. S. Nolan was supported by the NSF through Grant AGS-1654831. This work was supported by Keio University Academic Development Funds for Joint Research. The authors thank Dr. Jeff Kepert and two anonymous reviewers for their careful reviews and Drs. Hiroshi Taniguchi and Shigenori Otsuka for fruitful discussions.

### APPENDIX A

#### Equations for the One-Dimensional 1.5-Layer Shallow-Water System in Cartesian Coordinates

This appendix describes the governing equations for the one-dimensional, 1.5-layer shallow-water system on an *f* plane (cf. Fig. 1) considered in the main text. For the *x* direction that is orthogonal to the direction with fast flow speed, fluid depth and velocity are in geostrophic wind balance in both the free atmosphere and the boundary layer. A bulk drag formula is used to represent the surface friction. No internal diffusion is included in either layers. The mass and momentum conservation equations for the free atmosphere are

and those for the boundary layer are

where the asterisk indicates dimensional quantities, is the time, is the direction in which the fluid depth varies, is orthogonal to , is the fluid depth, and are the velocities in *x* and *y* directions, is the vertical velocity at the top of boundary layer, is the gravitational acceleration, is the Coriolis parameter, is the absolute vertical vorticity, is the drag coefficient, and the subscript *b* represents the boundary layer. We have assumed that the *y*-velocity and horizontal gradient of fluid depth in the boundary layer are equal to those in the free atmosphere.

The variables are normalized as

where *L*, *H*, *T*, and *V* are the representative scales of length, depth, time, and velocity, respectively. The dimensionless form of the equations is

and those for the boundary layer are

where is the dimensionless absolute vertical vorticity, is the Rossby number, is the Froude number, and .

To obtain the solution for Ekman pumping in the classic theory (e.g., Holton 1992), we linearize the equations by assuming that the basic-state velocity is negligible:

and those in the boundary layer are

where the bar and prime indicate the basic state and perturbation, respectively. Combining (A12) and (A14) yields the tendency equation for :

As a result, we obtain the tendency equation for the perturbation vorticity :

The equation can be integrated, and the solution is

where is initial value of relative vorticity perturbation. The decaying time scale of vorticity perturbation is short as the effect of surface friction is large. Since the semislip boundary condition is applied, faster velocity results in stronger friction.

### APPENDIX B

#### Condition for the Instability in Cartesian Coordinates

We here examine the condition for the instability in Cartesian coordinates, which is , where is defined in (10). The condition can be written as

which is identical to (12). The points of contact on the Ro axis of the curve defined by the lhs can be obtained as

We need to consider several different cases depending on the sign of coefficients of the first and second terms, whereas the number of cases can be reduced since Ro > 0. First, when the coefficient of the first term in (B1), , is positive, the inequality (B1) is satisfied only when the coefficient of the second term is negative. In this case, the unstable range is limited in , because the lhs of (B1) is concave upward. When the second term is positive, both the points of contact, and , are negative, and hence, Ro must be negative to satisfy (B1), which is not physical. In contrast, when the coefficient of the first term is negative, the lhs is concave downward and . Hence, the unstable range is much wider . More specifically, the unstable range is wider when the sign of the coefficient of the second term is negative compared with when the sign is positive. The points of contact on the Ro axis when the sign is positive and negative correspond to and of the first case with positive coefficient of , respectively. In conclusion, the inequality is more likely satisfied when both the coefficients, and , are negative. In fact, the eigenvalue analysis shows that the eigenmode of has a maximum where and (cf. Fig. 3).

### APPENDIX C

#### Test Calculations for the Instability

To verify the hypothesized sensitivity of solutions to the coefficient , some test calculations were performed by numerically integrating the equations for the free atmosphere [(3)–(5)] and vertical velocity [(8)] by artificially setting the coefficients , , and . For the calculations, the equations were combined into a single diagnostic equation for . Then the equations were discretized in time and space on a staggered grid system. The third-order Runge–Kutta method was used for time integration. The spatial derivatives were discretized by the fourth-order central difference scheme. The periodic boundary condition was applied to the lateral boundary. The domain width was 8 with 201 grid points and uniform grid spacing of 0.04. The integration time was 150 with the time step of 0.001.

The initial perturbation horizontal velocities were zero everywhere at the initial time. White noise perturbations with amplitude of 1.0 were added to the vertical velocity. Given , was diagnosed from the combined equation using the Gauss–Seidel method. The obtained was used to calculate temporal change in in (5). Then was obtained by using (8). The series of calculations was repeated until the end of integration.

The present hypothesis indicates that an instability occurs when is negative. Hence, we artificially set one of the three coefficients in (9)–(11) to ±1 and the other two to 0 rather than the values determined explicitly from velocity distribution. All combinations for the coefficients were tested. Figure C1 shows a time series of maximum vertical velocity in the test calculations. The experiment with is the only one that has an exponentially growing solution. The flow field (not shown) implies that the vertical velocity is amplified through the mechanism described above. In the other calculations, perturbations decay after the integration is initiated. The results show the presence of the linear instability when is negative.

### APPENDIX D

#### Equations for the Rapidly Rotating, Axisymmetric, 1.5-Layer Shallow-Water System

This appendix introduces the governing equations of a rapidly rotating, axisymmetric, 1.5-layer shallow-water system considered in this study. The mass and momentum conservation equations for the free atmosphere in cylindrical coordinates are given by

In the boundary layer, we assume that the tangential velocity is the same as that in the free atmosphere (i.e., ), and a bulk aerodynamic formula is used as the surface stress. The equations are

where is the time, is the radial direction, is the fluid depth, and are the radial and tangential velocities, is the vertical velocity at the top of boundary layer, is the gravitational acceleration, is the Coriolis parameter, is the absolute vertical vorticity, is the drag coefficient, the subscript *b* represents the boundary layer, and the asterisk indicates dimensional quantities. Equation (D5) indicates the balance between frictional destruction of absolute angular momentum and radial advection. This assumption has been applied by previous studies (Ooyama 1969). The equations are normalized as

where , and *T* are the scales for radius, the depth of free atmosphere, the boundary layer depth, and time, respectively.

The dimensionless equations for the free atmosphere are

and those for the boundary layer are

where , is the dimensionless absolute vertical vorticity, is the Rossby number, is the Froude number, and .

### APPENDIX E

#### Details of the WRF Simulation

The three-dimensional, full-physics simulation used in section 4b was produced using the WRF Model, version 3.4.1. The simulation depicts the early development and then rapid intensification of a tropical cyclone in a large, zonally periodic channel using the idealized modeling framework of Nolan (2011), with wind shear balanced by a meridional temperature gradient. However, in many other aspects, the simulation is modeled after the hurricane nature run simulations described by Nolan et al. (2013) and Nolan and Mattocks (2014). The simulation was originally intended to be one of a larger set of idealized simulations used by Koltz and Nolan (2018, manuscript submitted to *Wea. Forecasting*) for observing system tests, but its robust SEF and ERC made it useful for this paper.

The outer domain has 240 × 180 grid points in the zonal and meridional directions with 27-km grid spacing. Three nested, vortex-following grids are used with 9-, 3-, and 1-km grid spacings, with 180 × 180, 360 × 360, and 480 × 480 grid points, respectively. There are 60 vertical levels between the surface and 20-km altitude. The grid spacings in the WRF vertical coordinate are the same as used in the nature runs (see Fig. 2 of Nolan et al. 2013), and the physical parameterizations are also all identical, including the use of the one-dimensional mixed-layer cooling model of Pollard et al. (1972) available in WRF, version 3.4.1.

The initial vortex and background flow are very similar to the simulations described in Nolan (2011) and Nolan and McGauley (2012). The outer domain is initialized with mean easterly flow of 5 m s^{−1} at the surface, which smoothly increases to 5 m s^{−1} of westerly flow between 850 and 20 hPa. An axisymmetric vortex, balanced by temperature and pressure anomalies, is introduced into the eastern end of the domain. The symmetric initial vortex has a maximum tangential wind of 15 m s^{−1} at *r* = 135 km and *z* = 1.5 km. The radial profile of the tangential wind varies as a modified Rankine vortex with decay parameter *a* = 1/3 and decreases with height using the same analytical function as in (4.2) of Stern and Nolan (2011).

The Coriolis parameter is constant across the channel domain with value . The SST is set to 28°C along the centerline of the domain but varies meridionally. This meridional variation matches the temperature variation that exists at 5-km altitude because of the thermal wind balance. Nolan (2011) used this modification to suppress excessive convection in the northern part of the channel, where otherwise the temperature difference between the SST and the midtroposphere would be substantially increased. The SST variation is approximately linear, ranging from 31.5° to 24.5°C from the south to north edges of the channel.

Finally, as the simulation proceeds, the wind, temperature, and moisture fields on the outer domain are relaxed back to their initial values with a 24-h relaxation time scale. This relaxation does not apply on the nested grids. The relaxation keeps the environmental sounding and wind shear profile around the tropical cyclone roughly constant as the storm develops.

## REFERENCES

*An Introduction to Dynamic Meteorology.*3rd ed. Academic Press, 511 pp.

*Introduction to Circulating Atmospheres.*Cambridge University Press, 448 pp.

*Cyclones: Formation, Triggers and Control*, K. Oouchi and H. Fudeyasu, Eds., Nova Science Publishers, 1–36.

*31st Conf. on Hurricane and Tropical Meteorology*, San Diego, CA, Amer. Meteor. Soc., 91, https://ams.confex.com/ams/31Hurr/webprogram/Paper244751.html.

*Atmospheric and Oceanic Fluid Dynamics.*Cambridge University Press, 964 pp.

*Dynamics and Predictability of Large-Scale, High-Impact Weather and Climate Events*, Cambridge University Press, 168–175.

## Footnotes

^{a}

Current affiliation: Faculty of Environment and Information Studies, Keio University, Kanagawa, Japan.

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

We here consider the atmosphere in which the interior layer is located above the boundary layer. The basic concept of the present theory can be applied to the ocean and other flows.