Consequences of regularizing the Sawyer–Eliassen equation in balance models for tropical cyclone behaviour

The consequences of regularizing the Sawyer–Eliassen equation to calculate the stream function for the axisymmetric secondary circulation of a tropical cyclone are explored. Regularization is an ad hoc procedure in which the coefficients of the equation are suitably modified to replace negative values of the discriminant by small positive values, thereby ensuring that the equation is globally elliptic. The consequences of the procedure may be understood in terms of the analogue behaviour of a stretched membrane subject to a particular force distribution. Several regularization procedures are assessed by comparing the azimuthally averaged radial flow from a three‐dimensional numerical simulation of a tropical cyclone with that from an axisymmetric balance calculation of the Sawyer–Eliassen equation, forced by diabatic and frictional terms diagnosed from the simulation. The comparison shows that the largest challenge for regularization occurs in regions of inertial instability, especially when the diagnosed forcing overlaps with such regions. In the example shown, the diagnosed balanced flow is sensitive to the particular regularization procedure and none of the procedures examined give a flow that is structurally and quantitatively close to that obtained from the numerical solution in and near the region of regularization. The flow in regions of large vertical shear that are common in the lower part of the boundary layer is less sensitive to the regularization procedure, even though such regions are ones in which there is (frictional) forcing. Nevertheless, there are comparatively large differences between the low‐level inflow in the azimuthally averaged numerical solution and the axisymmetric balance solution. These differences can be attributed to the intrinsic lack of balance in the boundary layer. This finding, together with the issues associated with regularization, is further confirmation that balance dynamics is unable to adequately capture the flow in the boundary layer, contrary to recent claims.

wind balance (Willoughby, 1979). Regions where thermal wind balance does not hold include the frictional boundary layer and the upper tropospheric outflow layer. Assuming that thermal wind balance holds everywhere enables one to derive an equation for the stream function of the overturning circulation driven by diabatic heating and near-surface friction, processes that, in the absence of such a circulation, would drive the vortex away from thermal wind balance. This stream-function equation is generally referred to as the Sawyer-Eliassen (henceforth SE-) equation and has a variety of forms. For example, Sundqvist (1970a) derived a form of the equation in pressure coordinates, Shapiro and Willoughby (1982) used log-pressure as vertical coordinate, while Smith et al. (2005) derived a very general form in radius-height coordinates. A rather simple form can be derived by making the Boussinesq approximation (e.g. Montgomery and Smith, 2017) and a mathematically elegant form is obtained by using potential radius, R, instead of physical radius r (e.g. Schubert and Hack, 1982;Schubert and Alworth, 1987). The potential radius is defined by 1 2 2 = + 1 2 2 , where v the tangential velocity and f is the Coriolis parameter. Since R 2 = 2 M/f , where M is the absolute angular momentum, R-surfaces are surfaces of constant M.
The SE-equation is a key equation in the formulation of a prognostic axisymmetric balance theory for tropical cyclone evolution (Sundqvist, 1970a;1970b;Schubert and Alworth, 1987;Möller and Smith, 1994;Smith and Wang, 2018) and it has formed a basis for many diagnostic studies of tropical cyclone structure. In the latter studies, the SE-equation is solved diagnostically for the secondary circulation in the presence of a prescribed forcing mechanism (or mechanisms), possibly with an examination of the instantaneous tangential wind tendency accompanying the calculated overturning circulation (e.g. Smith, 1981;Schubert and Hack, 1982;Shapiro and Willoughby, 1982;Hack and Schubert, 1986;Rozoff et al., 2008;Bui et al., 2009;Pendergrass and Willoughby, 2009;Wang and Wang, 2013;Abarca and Montgomery, 2014;Smith et al., 2014;Ohno and Satoh, 2015;Heng and Wang, 2016;Heng et al., 2017). In some of these diagnostic studies, the axisymmetric vortex structure and the forcing functions (e.g. diabatic heating rate, frictional forcing) are obtained by a suitable azimuthal average of numerical model output (Bui et al., 2009;Wang and Wang, 2013;Abarca and Montgomery, 2014;Smith et al., 2014;Ohno and Satoh, 2015;Heng and Wang, 2016;Heng et al., 2017).
The solution of the SE-equation requires that the equation be globally elliptic, a condition that is usually satisfied by the choice of the vortex in idealized diagnostic studies, but is frequently not satisfied when the axisymmetric vortex structure is determined as an azimuthal average from the numerical model output of a tropical cyclone simulation. Even in prognostic balance theories that start from a state in which the SE-equation is globally elliptic, localized regions ultimately develop in which the ellipticity is violated Smith and Wang, 2018). When this happens, the solution can be carried forwards in time only by adjusting the coefficients in the SE-equation in the unstable regions to keep the equation elliptic globally. Nevertheless, regularization does not suppress the development of instabilities and the extended solutions ultimately break down, thereby limiting the time over which the balance model can be integrated.
In essence, regularization is an ad hoc procedure and various methods have been used. One method was devised by Möller and Shapiro (2002) in a case-study of Hurricane Opal (1995) and modifications thereof were used by Bui et al. (2009), Smith et al. (2014 and Smith and Wang (2018). An alternative method was suggested by Wirth and Dunkerton (2006), who effectively flattened out the M-surfaces in regions where the flow becomes inertially unstable, that is, where M/ r < 0. Despite the fact that the flattening out was accomplished using a scheme that globally conserves angular momentum, it makes the SE-equation parabolic in these regions, but not elliptic as required by the code they used to solve the equation. 1 Some of the authors referenced above have not checked whether their SE-equation is globally elliptic (e.g. Sundqvist, 1970a;1970b;Ohno and Satoh, 2015), raising questions about the convergence of their solutions.
The purpose of the present article is to develop a framework for exploring and understanding some of the local and global consequences of regularization and to investigate improved ways to carry out the regularization. In doing so, we highlight some fundamental limitations of regularization.

THE SAWYER-ELIASSEN EQUATION
The most general form of the SE-equation in cylindrical coordinates (r,z) may be written as where r is the radius, z is the height, is the density, is the stream function for the secondary circulation, = 1/ is the inverse of potential temperature , C = v 2 /r + fv is the 1 The solution code is described in the appendix of Wirth (1995), who states that "The resulting finite-difference equation is solved with the help of a multigrid algorithm (routine D03EDF) from the NAG Fortran library, which readily returns the desired solution so long as the equation is elliptic everywhere in the domain".
sum of centrifugal and Coriolis forces per unit mass, f is the Coriolis parameter (assumed constant), v is the tangential velocity component, = f + 2v/r is twice the absolute angular velocity, g is the acceleration due to gravity, = (1/r) (rv)/ r is the vertical component of relative vorticity,̇= d /dt is the material derivative of the diabatic heating rate anḋis the tangential momentum sink associated with the near-surface frictional stress.
where = 2 , = , = 2 g , = /( r), N 2 is the static stability, 2 g is the generalized inertial stability, and B is the baroclinicity. The last three quantities are given by the expressions: where I 2 = ( + f ) is the inertial stability squared. When has been determined, the radial and vertical velocity components, u and v, may be obtained using the formulae: u = −[1/( r)]( / z) and w = [1/( r)]( / r), respectively, which ensure that the continuity equation is satisfied.
The discriminant of the SE-equation, Δ, is given by The equation is locally elliptic if Δ > 0, locally hyperbolic if Δ < 0 and locally parabolic if Δ = 0. It can be shown that Δ is proportional to the potential vorticity, PV: that is, so that regions where the SE-equation is hyperbolic correspond with regions of negative PV, equivalent to the flow being symmetrically unstable. Regions where Δ < 0 are where the flow is inertially unstable ( 2 g < 0), statically unstable (N 2 < 0) or where the baroclinicity, a measure of the vertical shear, is sufficiently large ( 2 > ). In general, for tropical-cyclone-scale vortices, the coefficients of the highest derivatives in the SE-equation are functions of r and z, and numerical methods are called for to obtain solutions. Moreover, the complex nature of the coefficients makes it difficult to determine the consequences of any regularization method. For that reason, it is helpful to step back and investigate an analogous problem with a simpler partial differential equation.

THE MEMBRANE ANALOGY
One of the simplest physical problems for understanding the behaviour of elliptic second-order partial differential equations is the equilibrium displacement of a stretched membrane subject to a distribution of forces normal to the membrane. Two examples, those of a point force and point force dipole are sketched in Figure 1. In a rectangular coordinate system (x,y), the membrane displacement Z(x,y) satisfies the Poisson equation: where F(x, y) is the imposed force. Here, positive F corresponds to an upward force acting on the membrane. In the case of a square domain (0 ≤ x ≤ 1, 0 ≤ y ≤ 1) with zero displacement along the boundary and a point force at the centre [ ( , ) = ( − 1 2 ) ( − 1 2 )], one can use one's intuition to see that the solution for the membrane displacement has to be a maximum at the point of forcing with closed contours that are near-circular in the vicinity of the forcing and approach a square with smoothed corners near the boundaries. This intuition is confirmed by the numerical solution for a concentrated forcing 2 shown in Figure 2a. This and other solutions that follow are obtained using the same over-relaxation procedure described by Bui et al. (2009).
For a membrane with the property that it deforms more easily in the x-direction than in the y-direction, the membrane displacement satisfies an equation of the type 2 2 + 2 2 2 = − ( , ), where is a constant smaller than unity. In the case where = 0.1, the solution with the forcing function in Figure 2a   shown in Figure 2b. In this case, the maximum displacement amplitude has increased and the membrane displacement has become confined in the y-direction, barely feeling the boundaries in that direction. Note that a transformation of the y-coordinate in Equation 7 to Y = y/ would lead to the same equation as Equation 6, but with Y replacing y and in (x,Y) space the solution would be similar to that in Figure 2a, but the domain would be larger in the Y-direction. It follows that the solution of Equation 7 is simply a stretched version of Equation 6 in the y-direction if > 1 and a shrunken version of Equation 6 if < 1. These solutions may be used to understand the consequences of regularization, which would be equivalent to solving Equation 6 over much of the domain, but solving Equation 7 over a limited region with some small value of . We shall refer to this region as the" region of regularization" and take it to be a square that includes or excludes the small region of forcing shown in Figure 2a. Figure 3a shows the solution when the region of regularization is confined to the dot-dashed square shown. In that region the displacement contours are flattened as in Figure 2b, but as the boundaries are approached the solution is similar to that in Figure 2a. Nevertheless, as in Figure 2b, the amplitude of the maximum displacement is larger than in Figure 2a. The effect of flattening of the displacement contours is seen in the difference field shown in Figure 3b.
When the region of regularization is situated away from and to the left of the forcing (Figure 3c), the maximum displacement is still larger than in Figure 2a, and in fact, the displacement is larger everywhere with the maximum difference located inside the region of regularization ( Figure 3d). When the region of regularization is to the right of the forcing, one may expect a similar pattern of displacement, but with the enhanced values to the right instead of the left of the forcing.
When the region of regularization is situated away from and below the forcing (Figure 3e), the maximum displacement is smaller than in Figure 2a and the displacement contours are again flattened out inside the region of regularization. The effect in this location is to reduce the amplitude of the displacement everywhere (Figure 3f), the maximum difference between the control case being on the border of the domain of regularization closest to the location of forcing. When the region of regularization is above the forcing, one may expect a similar pattern of displacement, but mirror imaged in the x-axis.
When the region of regularization is located on the diagonal to the right of and above the forcing (Figure 4a), the flattening of the contours there leads to a dipole pattern of deviation displacement from the control calculation ( Figure 4b) with a negative deviation in the upper portion of the regularization region and a positive deviation in the lower half, the maximum being located on the lower boundary of that region.

A MORE REALISTIC CONFIGURATION VIS-Á-VIS THE ATMOSPHERE
In the atmosphere, diabatic heating appears in the SE-equation as a dipole of forcing oriented principally in the radial direction. Moreover, a domain that has a large aspect ratio (length to depth) and is open at its lateral boundary is more appropriate in the atmospheric context. Assuming a two-dimensional flow configuration in rectangular coordinates (x,z), the SE-equation for a resting atmosphere would have the form The dipole forcing corresponding to that produced by an idealized line of diabatic heating from deep convection would look something like that in Figure 5a, the dipole being related primarily to the radial gradient of the heating (see Equation 1). The forcing is located relatively close to the z-axis which is chosen to be a closed boundary ( = 0) analogous in rectangular geometry to the axis of rotation of an axisymmetric vortex in cylindrical coordinates. The solution for the stream function induced by this forcing for 3 = 0.1 is shown in Figure 5b, assuming that the right boundary of the domain is open and that / x = 0 along it. The upper and lower boundaries are taken to be closed with = 0 there. The pattern of "lateral velocity", u = − / z, corresponding with the stream function is shown in Figure 5c. As expected, and in analogy to the situation in axisymmetric geometry 3 In the atmosphere in middle latitudes, a more typical value for would be Figure 3, but with the region with = 0.1 moved along the diagonal to the right of and above the forcing (e.g. Shapiro and Willoughby, 1982), the stream function shows two cells of circulation with ascent along the axis of the forcing and within the forcing region itself, and descent elsewhere. Beyond the forcing there is inflow in the lower troposphere and outflow in the upper troposphere and this inflow and outflow pattern extends to the right boundary. Decreasing the value of would increase the lateral scale of the outer circulation cell, leading to a larger flow through the right boundary and less recirculation within the domain (not shown). Figure 5d shows the analogous solution when the coefficient in Equation 8 is reduced to a constant value 0.01 in the rectangle in the "upper troposphere" shown. This configuration is analogous to the procedure of regularizing the SE-equation in regions in the upper troposphere where the flow becomes inertially unstable, equivalent in Equation 8 to 2 becoming locally negative. As described in Möller and Shapiro (2002), regularization involves effectively setting 2 equal to some small positive value in such a region. From the understanding gained in section 3, we know that regularization in the rectangular region shown in Figure 5d will have the effect of flattening out the streamlines in the rectangle. Moreover, the effect is global, but diminishes in magnitude with increasing distance from the region of regularization. As can be seen in the figure, this is precisely what happens. Figure 5e shows the lateral velocity component, u = − / z, derived from the stream function shown in Figure 5d. 4 The effect of "regularization" is to destroy the symmetry of the inflow and outflow regions beyond the forcing, leading to stronger outflow in the upper troposphere, albeit concentrated in a shallower layer than the inflow. The difference between the pattern of inflow and outflow between the regularized solution in Figure 5e and the unregularized solution in Figure 5c is shown in Figure 5f. Significantly, the outflow is strengthened throughout much of the upper troposphere, not only in the region of regularization, with the maximum increase near the lower boundary of the region of regularization. Elsewhere, the radial flow is decreased, that is the inflow has increased, but the maximum decrease occurs a little below the region of regularization. Figure 5g-j show similar fields to those in Figure 5d,e, but where the region of regularization is moved radially outwards (g,h) or inwards (i,j). When the region of regularization is moved outwards, the radial flow is able to rise higher before the flattening occurs (compare (g) with (d)), but when the region is moved inwards with the inner boundary at the axis of forcing, the flattening occurs almost immediately as the flow exits the updraught produced by the forcing (compare (i) with (d)), The consequences for the radial flow are shown in (h,j), respectively. In the former case, the outflow layer extends over a deeper layer than in (e), but the outflow is weaker, whereas, in the latter case, the radial flow is more confined in the vertical, but much stronger than in (e).
As will be discussed in section 6, the structural changes brought about by regularization shown in Figures 5 provide an understanding of possible consequences of regularization in solving the SE-equation itself.

METHODS OF REGULARIZATION
As discussed in section 1, the main purpose of regularization in solving the SE-equation is to remove any regions of symmetric instability (Δ < 0). This removal can be achieved by replacing the corresponding negative coefficients (N 2 or 2 g ) with small positive values and/or by suitably decreasing the coefficient B. Alternatively, the removal can be achieved by sufficiently increasing the magnitude of the inertial stability, static stability, or both (a procedure adopted by Möller and Shapiro (2002)).   and Smith and Wang (2018) calculate the minimum value of 2 g in the region where Δ < 0, say 2 g min , then remove the negative values of 2 g by adding | 1.001 2 g min |. Further, at points where N 2 < 0, which typically do not coincide with those where 2 g < 0, N 2 is set equal to 10 −8 s −2 . Finally, if Δ is still less than or equal to zero, which, when | 2 g | is made small and positive is frequently the case, 2 is replaced with 1 2 at the grid point in question. Heng and Wang (2016)  A different procedure is adopted by Möller and Shapiro (2002). In regions where Δ < 0 they increased the value of (and thereby a ) so that, effectively, 5 Δ has some small positive threshold value. No other quantities appearing in the SE-equation coefficients are altered so that, in particular, the new value of is not consistent with the local structure of v. Of course, this is a property of the other schemes as well. Apparently, in the vortex examined by Möller and Shapiro, regions of negative Δ were due largely to the occurrence of inertial stability, 2 g < 0. The hope in all these studies has been that, whatever procedure is used, regularization will lead to a useful balanced solution, at least in regions remote from those where regularization is needed, but there are some subtle differences between the procedures that have consequences for the diagnosed structure of the balanced solution.
As noted above, the Möller and Shapiro procedure differs from the others in that, at points where Δ < 0, 2 g is increased in magnitude, even if the point with Δ < 0 is a consequence of large vertical shear and not necessarily because the flow is inertially unstable ( 2 g < 0). In contrast, in the other schemes, points with inertial instability are removed first by setting 2 g to be a small positive number, typically smaller in magnitude than its original magnitude. The analysis in section 3 points to a different local response to forcing depending on which regularization procedure is adopted and therefore to a different structure of the balanced solution in and near the region where Δ < 0.

AN IDEALIZED THREE-DIMENSIONAL NUMERICAL SIMULATION OF A TROPICAL CYCLONE
We apply now the insights gained above to assess the applicability of balance theory in analysing the secondary circulation of an idealized three-dimensional numerical simulation of tropical cyclone evolution on an f -plane. The simulation is similar to the one described by Kilroy et al. (2016), but uses the CM1 model (Bryan and Fritsch, 2002) with a horizontal grid spacing of 1 km and a vertical grid spacing of 100 m. These data were kindly provided by Gerard Kilroy. Figure 6a shows the azimuthally averaged and 3 h time-averaged radial and tangential velocity components from the foregoing simulation at 32 h. At this time the vortex was undergoing a period of rapid intensification. The main features of the simulation are similar to those described in many previous studies (see Montgomery and Smith, 2017 and references). There is a shallow layer of strong inflow near the surface and one of strong outflow centred at about 12 km. There is a shallow region of marked outflow just above the boundary layer, where the inflow terminates and ascends into the developing eyewall updraught. This updraught is indicated by the contour of vertical velocity equal to 0.25 m⋅s −1 . There is a region of weaker inflow in the lower troposphere, with a shallow layer of weak inflow just below the main outflow layer. The maximum tangential wind speed occurs at a height 600 m and radius of 38 km, within the layer of strong inflow.
The mean tangential wind field in Figure 6a is used to obtain balanced density and potential temperature fields using the method described by Smith (2006). In turn, these fields are used to evaluate the coefficients on the left-hand side of the SE-equation. The forcing terms on the right-hand side of the SE-equation arising from diabatic heating and friction are diagnosed also from the time-and azimuthally averaged model output. The structure of the combined forcing is shown in Figure 6b. The main region of positive forcing is near the inner edge of the main region of ascent, the region within the yellow contour in Figure 6a. At larger radii, mostly beyond a radius of 28 km and inside a radius of about 90 km, there are narrow strips of negative forcing, punctuated by even narrower strips of positive forcing. There is a shallow region of negative forcing below a height of about 1.5 km and inside a radius of 50 km. This feature is associated with the inner-core boundary layer. Figure 6b shows also the regions where the discriminant of the SE-equation, Δ, is negative. The main area of negative Δ is located in the mid-to upper troposphere between radii of approximately 70 to 180 km, much of it overlapping with the main outflow layer. This region, together with a much smaller region near the outer boundary between about 9 and 10 km in height, is associated with the generalized inertial stability, 2 g , being negative. A shallow finger of negative Δ located just above 2 km height and extending to nearly 30 km in radius is associated with static instability N 2 < 0 and a shallow (less than 400 m deep) surface-based layer of negative Δ is associated with large vertical shear where 2 > . All of these regions require regularization in order to solve the SE-equation.
Because the region of static instability is so small, the flow therein appears to be little influenced by the regularization. For this reason we do not examine other methods to regularize the SE-equation in such regions. More details about the regularization of the equation in regions where 2 > are discussed in section 7.  Figure 6b using two regularization schemes. They show also fields of the ratio 2 g /N 2 , which is the same as the ratio ∕ in Equation 2 and is analogous to the quantity 2 in Equation 8. In the scheme in Figure 6c, which we refer to as Scheme A, regions of negative 2 g ( , ) are removed by adding | 1.001 2 g min |. In Figure 6d, a new procedure is adopted in which negative values of 2 g ( , ) are removed by adding the local value | 1.001 2 g ( , ) |. This procedure, which we refer to as Scheme B, has the advantage of avoiding artificially For v, every 5 m⋅s −1 from 0 to 30 m⋅s −1 ; for u, every 0.5 m⋅s −1 from ±0.5 to ±2 m⋅s −1 and every 3 m⋅s −1 from ±2 to ±20 m⋅s −1 . For 2 g ∕ 2 , every 0.2 units from 0.2 to 1 unit, 1 unit = 1 × 10 −3 . Shading as shown in the colour bar sharp gradients of 2 g at the boundary of the region where Δ < 0. Even so, the reduction in the magnitude of 2 g generally requires a reduction in the magnitude of B to keep Δ > 0.

Two regularization schemes
In this section, 2 is replaced with 0.99 at each grid point where Δ remains negative after modifying 2 g . Other possibilities are explored in section 7.
While both schemes capture the broad features of the secondary overturning circulation, the flow structure in the upper troposphere shows considerable differences, principally in the region of regularization and regions adjacent to it. As shown in section 3, a small inertial stability as in the modified scheme provides for an enhanced response of the radial velocity component to the forcing, while the larger inertial stability in the original scheme acts to inhibit the radial flow. Figure 5f provides a clue to understanding this behaviour. It shows that there is enhanced radial outflow flow just inside the region of reduced inertial instability and enhanced radial inflow or reduced outflow below that region. Conversely, if the inertial stability in the region of regularization is increased in magnitude beyond that of the surrounding values, there is enhanced radial inflow or reduced outflow above the boundary of regularization (not shown).
The new regularization scheme appears to reproduce the flow structure in the numerical model somewhat better than the original scheme, although the layer of inflow just below the main outflow layer is much too strong. The maximum outflow in the upper troposphere in the balance solution is 20.2 m⋅s −1 compared with 11.2 m⋅s −1 in the numerical model, while the maximum upper-level inflow is 10.7 m⋅s −1 compared with only 1.8 m⋅s −1 in the numerical model. Thus, even with the new procedure for replacing negative values of 2 g ( , ), the regularized balance solution does a relatively poor job in capturing the outflow and inflow strengths in the numerical model.

The Möller and Shapiro scheme
An alternative regularization scheme, which we refer to as Scheme C, is to set = 2 ∕0.99 , whereupon it is not necessary to change . This scheme was suggested by Möller and Shapiro (2002). The results of this scheme are shown in Figure 6e. While the radial flow structure in this figure is closer to that in Figure 6d than that in Figure 6c, it is no improvement in relation to the numerical solution in Figure 6a. In this case, the maximum outflow in the upper troposphere is 20.2 m⋅s −1 , the same as before, but the maximum upper-level inflow is slightly larger, 11.2 m⋅s −1 instead of 10.7 m⋅s −1 , making the agreement with the numerical solution slightly worse.

6.3
The issue of forcing overlapping with regions in which < 0 As indicated in Figure 6b, there is considerable overlap between the total forcing distribution due to heating and friction and the primary region where the flow is inertially unstable. Figure 3a shows that this is a situation where the response to the forcing is particularly large in amplitude. This finding may explain why the magnitude of upper-tropospheric inflow and outflow shown in Figure 6e is overestimated. To examine this possibility, we show in Figure 6f the solution to the SE-equation analogous to that in Figure 6e, but with the forcing function set equal to zero in the upper-level region of non-positive discriminant. While the maximum outflow and maximum inflow are indeed reduced in comparison with those in Figure 6e, the second layer of outflow centred at a level of about 8 km has strengthened considerably and this layer is not even present in the numerical calculation in Figure 6a. This feature is presumably a result of the artificially large vertical gradient of the forcing on the boundary of the main regularization region, which is introduced by setting the forcing abruptly to zero inside the region of regularization.
The inability of the SE-calculation to capture quantitatively the upper-level structure seen in the numerical calculation could be the fact that the flow in the numerical model is nowhere near axisymmetric at the time shown (32 h). An alternative, but not necessarily mutually exclusive explanation would be that the inability is simply a consequence of regularizing the SE-equation. Based on the understanding gained in section 4, this would seem to be the most likely scenario, since the flow in the lower half of the troposphere is somewhat better captured by the balance calculation, except in a shallow layer near the surface. The flow in the near-surface layer, which is one that overlaps also with negative discriminant of the SE-equation, is examined in the next section.
Despite the large differences in the structure of inflow and outflow in the middle and upper troposphere in Figure 6c-e as a result of the different regularization schemes, there is little difference in the lower troposphere and there are only small differences in Figure 6f, in which the forcing is suppressed in the region requiring regularization. We conclude that the boundary-layer inflow is at most weakly influenced by the regularization of regions of inertial instability in the upper troposphere. This is counter to the claim by Heng et al. (2017) that " … the boundary layer inflow in the balanced response is very sensitive to the adjustment to inertial stability in the upper troposphere … ". A more detailed examination of boundary-layer structure is shown in the next section.

REGULARIZATION IN REGIONS OF LARGE VERTICAL SHEAR
In regions of large vertical shear, Δ may become negative on account of 2 exceeding . Typically, such regions occur in a shallow surface-based layer within the friction layer itself. As explained in section 5, one method for removing the negative discriminant is to set 2 = 1 2 or, perhaps preferably, 2 = 0.99 to make Δ positive, but closer to zero.

Exploitation of the membrane analogy
The consequences of redefining are illustrated in three idealized calculations shown in Figure 7. These calculations involve solutions of Equation 8 with 2 = 0.01, as in Figure 5, and with an idealized surface-based layer of forcing, −F(x, z) shown in Figure 7a. This forcing distribution is analogous to that involving the vertical gradient oḟin Equation 1. Figure 7b shows the stream function (x, z) induced by the forcing distribution shown in (a). Figure 7c shows the radial velocity derived from the stream function shown in (b), highlighting the fact that there is inflow in the region of forcing and outflow above it, but because of the implied strong vertical stability in the value chosen for , the maximum outflow occurs at low levels, just above the layer of forcing. Figure 7d shows the stream function for a similar calculation to that in (b), but when a term 2 2 ∕ is added to the left-hand side of Equation 7 in a layer that has half the depth of the layer of forcing. The constant is chosen to be equivalent to setting 2 = 0.99 in the SE-equation in cases where the vertical shear is large and would otherwise make the discriminant Δ negative. Moreover, when taking the square root, the sign of should be preserved as it was before regularization.
Comparison of Figure 7d with Figure 7b shows that the effect on the stream function from the inclusion of the term involving in Equation 7 is minimal, producing a slight clockwise rotation of the streamlines in the region of non-zero . Such rotation was explained in the classic paper by Shapiro  Willoughby (1982): see especially their Figure 1 and related discussion. The effect is mainly discernible in the slight elevation of the streamlines in the inner region (x < 5) and in the slight depression of the streamlines in the outer region (x > 5). Figure 7e shows the lateral component of flow in this case, which should be compared with Figure 7c. In essence, the "regularization" has reduced both the surface-based inflow and the outflow above it on the inner side of the forcing and has enhanced both the inflow and outflow on the outer side of the forcing. These effects are highlighted in (f), which shows the difference between the lateral flow in (e,c). A comparison of Figures 7g,h indicate that the "inflow" on the right side of the regularization region has been strengthened, which is the situation in our simulation as discussed in the next subsection.

7.2
Low-level comparison between the numerical simulation and the balance calculation Figure 8 shows similar fields to those in Figure 6, but focusing on the low-level flow structure in the numerical model simulation and in the calculation of the balanced response to the total forcing due to heating and friction. Figure 8a shows the flow structure of Figure 6a in the lowest 3 km, while Figure 8b shows the structure of the forcing (Figure 6b) in this region together with the regions where the SE-equation requires regularization. The region where large vertical shear leads to a need for regularization is rather shallow, less than 400 m deep, extending from a radius near 30 km. Based on the idealized calculations in Figure 7, the effect of the regularization required in this layer would be expected to be minimal and unlikely to account for the difference in low-level structure between the numerical solution in Figure 8a and the balance solution shown in Figure 8c. (The latter figure shows just the lower 3 km of Figure 6d.) While the maximum inflow in the numerical solution is 11 m⋅s −1 , that in the balance calculation is only about 8.3 m⋅s −1 . Moreover, the radial location of the maximum inflow occurs at a much smaller radius (41 km) compared with the radius in the balance calculation (174 km). One possible reason is that the inflow inside the regularization area to the right side of the momentum forcing maximum has been for v, every 5 m⋅s −1 from 0 to 30 m⋅s −1 ; for u, every 0.2 m⋅s −1 from ±0.2 to ±1 m⋅s −1 and every 1 m⋅s −1 from ±1 to ±10 m⋅s −1 . Shading as shown in the colour bar enhanced, as in Figure 7e. However, the analysis of the previous subsection suggests that this effect would not be large enough to explain the large difference between Figure 8c and Figure 8a. Figure 8d shows a similar regularization scheme to Scheme A, but in regions of large baroclinicity, 2 is replaced with 0.5 . We refer to this as Scheme D. This is the scheme used by Bui et al. (2009 and Smith and Wang (2018). There are only small differences from the fields shown in (c) and these are confined to the vicinity of the regularization region. From this result it would appear that the flow in regions of large vertical shear that are common in the lower part of the boundary layer is less sensitive to the regularization procedure than that in regions of inertial instability.
In support of our conclusion at the end of subsection 6.3 that the boundary-layer inflow is at most weakly influenced by the regularization of regions of inertial instability in the upper troposphere, Figure 8e,f show just the lower 3 km of Figure 6c,e. The boundary-layer structures in Figure 6c-e are almost identical and they even have the same magnitude of maximum inflow (8.3 m⋅s −1 ) at the same radius (174 km).
Clearly, the balance solution poorly captures the boundary-layer inflow in the numerical calculation, a finding consistent with the study of Bui et al. (2009) and the more recent calculations of Montgomery and Persing (personal communication, 2019). The finding is clearly at odds with one of Heng et al. (2017) who claim that "balanced dynamics can well capture the secondary circulation in the full-physics model simulation even in the inner-core region in the boundary layer", but is supported by a scaling analysis of the boundary-layer equations, which shows that the unbalanced (nonlinear) terms are important in the inner-core region of a tropical cyclone (Vogl and Smith, 2009). It is supported also by the finding of Vogl and Smith (2009) that even a linear (but unbalanced) approximation to the boundary-layer equations is a poor representation of the inner-core boundary layer of a tropical cyclone.

DISCUSSION AND CONCLUSIONS
We have developed a framework for exploring the consequences of regularizing the Sawyer-Eliassen equation to diagnose the stream function for the axisymmetric secondary circulation of a tropical cyclone subject to a given distribution of diabatic forcing and tangential frictional stress. Regularization amounts to adjusting the coefficients of the equation in regions where the discriminant is negative to ensure that the equation is globally elliptic. The possible consequences of regularization have been explored using the analogue behaviour of a stretched membrane subject to a particular force distribution.
Regularization is required in three regions: (a) regions where the flow is inertially unstable, (b) regions where it is statically unstable, and (c) regions where the baroclinicity is large. Regions of large baroclinicity are typically ones of large vertical shear. In numerical models of tropical cyclones, regions of azimuthally averaged inertial instability are generally the most extensive, while regions of static instability are typically small in areal extent. Regions where the azimuthally averaged baroclinicity is large are typically confined to the lower part of the frictional boundary layer, where the vertical shear is large. However, setting the inertial stability to be small and positive in regions of inertial instability generally requires the baroclinicity to be reduced in magnitude as well to keep the discriminant of the Sawyer-Eliassen equation positive. Possible improvements in the procedure for regularizing in cases (1) and (3) were suggested.
A comparison of the azimuthally averaged radial flow from a three-dimensional numerical simulation of a tropical cyclone with those from an axisymmetric balance calculation of the Sawyer-Eliassen equation forced by diabatic and frictional terms from the numerical simulation was presented. Important findings from this comparison are: 1. The largest uncertainty in the integrity of the balance solutions results from the regularization in regions of inertial instability, especially when the diagnosed forcing overlaps with such regions. In the example shown, where there is some overlap of this type, the diagnosed balanced flow is sensitive to the particular procedure for regularization and none of the schemes produced a flow that was structurally and quantitatively close to that obtained from the numerical solution. 2. Regularization in regions of large vertical shear that typically occur in the lower part of the boundary layer is less problematic, even though such regions are ones in which there is forcing. The reason is that a modification of the coefficient B in the SE-equation leads to a rotation of the stream function response, but the degree of rotation is constrained by the proximity of the lower boundary. 3. On account of (2), the large difference found between the low-level inflow in the azimuthally averaged numerical solution and that in the axisymmetric balance solution is further indication that balance dynamics is unable to adequately capture the flow in the boundary layer, contrary to recent claims.
While balance ideas have played a central role in the development of a theoretical framework for understanding tropical cyclone dynamics, the application of such ideas to diagnose the results of numerical simulations almost always requires that the Sawyer-Eliassen equation be regularized. Regularization is intrinsically an ad hoc procedure and some methods may be better than others. Exploitation of the membrane analogy as outlined herein would seem to offer a useful framework for assessing the integrity of such procedures and their possible limitations. Our analysis suggests, however, that regularization introduces uncertainties in the integrity of balance solutions to a degree that much caution is called for in the use of such solutions for "explaining" tropical cyclone structure, especially within and near the regions which have been regularized.