Intensification of tilted atmospheric vortices by asymmetric diabatic heating
Tom Dörffel Ariane Papke Rupert Klein Natalia Ernst Piotr Smolarkiewicz












T. Dörffel and A. Papke and R. Klein FB für Mathematik und Informatik, Freie Universität Berlin, Berlin, Germany

              rupert.klein@fu-berlin.de
          N. Ernst Zuse Insitute Berlin, Berlin, Germany 

          P. Smolarkiewicz National Center for Atmospheric Research, Boulder, CO 80307, USA





Intensification of tilted atmospheric vortices by asymmetric diabatic heating
    Tom Dörffel 
        Ariane Papke 
        Rupert Klein 
        Natalia Ernst 
        Piotr K. Smolarkiewicz
    Received: date / Accepted: date
==================================================================================================================









<cit.> studied the nonlinear dynamics
of strongly tilted vortices subject to asymmetric diabatic heating by asymptotic methods.
They found, , that an azimuthal Fourier mode 1 heating pattern can intensify or
attenuate such a vortex depending on the relative orientation of tilt and heating
asymmetries. The theory originally addressed the gradient wind regime which,
asymptotically speaking, corresponds to vortex Rossby numbers of order unity
in the limit. Formally, this restricts the applicability of the theory to rather weak
vortices in the near equatorial region. It is shown below that said theory is, in
contrast, uniformly valid for vanishing Coriolis parameter and thus applicable to
vortices up to low hurricane strengths. In addition, the paper presents an
extended discussion of the asymptotics as regards their physical interpretation and
their implications for the overall vortex dynamics. The paper's second
contribution is a series of three-dimensional numerical simulations examining the effect
of different orientations of dipolar heat release on idealized tropical cyclones.
Comparisons with numerical solutions of the asymptotic equations yield evidence that supports the original predictions. In addition, the influence of asymmetric diabatic heat release on the time evolution of centerline tilt is analysed further, and a steering mechanism based on the orientation of the heating dipole is revealed.








§ INTRODUCTION



Atmospheric vortex intensification and the associated evolution of vortex structure remain
a topic of intense investigations. As <cit.> point out in their
review article, intricate interactions of boundary layer processes, moist thermodynamics,
multiscale stochastic deep convection, and the vortex-scale fluid dynamics produce the
observed, sometimes extremely rapid intensification of incipient hurricanes. They also
emphasize that, despite the valuable insights that have been gained in many studies of
idealized axisymmetric flow models, asymmetries of vortex structure, convection patterns,
and boundary layer structure have been observed to be important for vortex intensification
in real-life situations.

This study focuses on the principal response mechanisms of Tropical cyclone-like (TC) atmospheric vortices
to asymmetric diabatic heat release following the theory of <cit.>.
Therefore, we analyze both the structure and intensity of the bulk vortex above the
boundary layer under the influence of different configurations of asymmetric heating.
From <cit.>, among others, we adopt
the point of view that latent heat release from condensation can be modeled, with
limitations, by external diabatic heat sources in dry air. In the cited studies,
non-axisymmetric heating patterns were shown to have at most a small effect on vortex
strength within the framework of linearizations about an axisymmetric upright vortex.
These results of linear theory were corroborated in <cit.>
by comparison with fully nonlinear three-dimensional simulations.
By both, analytical and numerical examination, we will see that the particular flow
structure of a strongly tilted vortex allows for a leading order intensification mechanism
based on asymmetric heating that cannot be captured for linearized (weak) vortex tilts.

Investigating incipient hurricanes that develop in the tropical
Atlantic, <cit.> revealed,
that such vortices can exhibit
very strong tilt. Thus, for instance, the locations of the vortex center at heights
equivalent to the 925 and 200 pressure
levels are located about 200 apart, e.g., in
<cit.> and <cit.>.
This amounts to an overall vortex tilt at a scale comparable to the vortex diameter,
, to a situation that clearly does not allow for linearizations about an upright
vortex. In fact, in characterizing the wind field of Hurricane Norbert,
<cit.> already utilized the concept of a height-dependent
vortex center, i.e., of time dependent centerline, and this is one of the key
structural aspects in the analysis of <cit.>, which we revisit in
this paper.

<cit.> analyzed the
dynamics of such strongly tilted atmospheric vortices in the gradient wind regime by
matched asymptotic expansions. They obtained a closed coupled set of evolution
equations for the primary circulation structure and the vortex centerline, and
demonstrated that in a strongly tilted vortex symmetric and asymmetric heating
patterns can have a comparable impact on vortex intensity. As by its very definition
the gradient wind regime is restricted to vortex Rossby numbers of order unity,
this theory has thus far been considered applicable only to rather weak vortices
with intensities relatively far from the interesting stage of the tropical
storm/hurricane transition <cit.>.

To allow for vortices in this transition regime, we consider here in the
first part of the paper the dynamics
of meso-scale atmospheric vortices ∼100 that extend
vertically across the depth of the troposphere ∼10
but feature large vortex Rossby number ≫ 1. We use the asymptotic
techniques introduced by <cit.> and recycle many of their technical
steps. As indicated in fig. <ref>, we assume vortices with nearly
axisymmetric core structure at each horizontal level, and we allow for strong vortex
tilt such that the vortex centers observed at different heights may be displaced
horizontally relative to each other by distances comparable to the vortex core
size .




One of the main findings of <cit.> was the following evolution
equation for the primary circulation described by the axisymmetric leading-order
circumferential velocity, , valid for time scales large compared
to the vortex turnover time scale,


    u_θt
    + w_0u_θz
    + u_r,00(u_θr + u_θ/r + f_0
      )
    = - u_r,*(u_θ/r + f_0
      )  .


Here (t,r,z) are the appropriately rescaled time, radial, and vertical coordinates,
f_0 is the Coriolis parameter, and w_0 and u_r,00 are the axisymmetric components of
the vertical and radial velocities induced by the axisymmetric mean heating patterns,
also properly rescaled. The apparent radial velocity
u_r,* results from an interaction of the vortex tilt with the asymmetric first
circumferential Fourier mode of the vertical velocity. In particular,


    u_r,*
    = 1/2π∫_-π^π w  ·z  dθ  ,


where (t,z) is the time dependent vortex centerline position at height z
(see fig. <ref>), w is the full vertical velocity, and
= cos(θ) + sin(θ) is the radial unit vector of a
horizontal polar coordinate system attached to the centerline.

X(t,z) itself is governed by the centerline equation:


    Xt =  u_s + ( X ·∇)  u_s + ln1/δ k × M_1 +  k ×


u_s(t,z) expresses the background wind profile, δ is a small number according to the asymptotic scaling, k the vertical unit vector, M_1 is a weighted curvature measure of X and  evaluates Fourier-1 modes of vertical velocity resulting from both, diabatic heating and adiabatic balances within a tilted vortex.
In the adiabatic case and without vertical wind shear equation (<ref>) simplifies to a linear Schrödinger-like equation exhibiting undamped precession of eigenmodes.
More details on the expressions for M_1 and  follow in the further course of this article.

The main findings of present work are:




  * 
The evolution equation from eq:TempevoluthetaIntro is uniformly valid
as f_0 → 0 so that it holds, in particular, also for ≫ 1, , for
vortices of hurricane strength.
In fact, we argue the for ≥1 the structure of the leading-order equations does not change.


  * 
The mechanism of vortex spin-up by asymmetric heating of a tilted vortex is
traced back analytically to an effective circumferential mean vertical mass
flux divergence that arises when the first Fourier mode diabatic heating
and the vortex tilt correlate negatively.


  * Asymmetric heating introduces a forcing of the vortex motion which can accelerate/decelerate the centerline precession and/or increase/decrease its tilt depending on the relative orientation of tilt and heating dipole.


  * 
Equation eq:TempevoluthetaIntro can be recast into a balance equation
for kinetic energy, = ρ_0^2/2,


    (r 0pt9pt)_t
    + (r u_r,00[ + p]0pt9pt)_r
    + (r _0 [ + p]0pt9pt)_z
    =  rρ_0/N^2^2(Θ· Q_Θ)_0


in line with the theory by <cit.> for available potential energy (APE) generation.

Here p is the relevant pressure perturbation,
Θ, Q_Θ are the potential temperature perturbations and the diabatic heating, respectively, and (·)_0 corresponds to the axisymmetric mean.
In the current case, it encodes the correlation of potential temperature perturbation and diabatic heating.
N and  are the Brunt-Väisälä frequency
and the background potential temperature stratification, respectively.
Equation eq:KineticEnergyBudget states that, except for a conservative
redistribution of kinetic energy due to advection and the work of the pressure
perturbation, p̃, positive correlations of diabatic sources and potential temperature
perturbations generate the potential energy available for increasing the kinetic
energy of the vortex.

<cit.> study the effects of asymmetric diabatic heating on vortex
strength in a linearized model. One of their conclusions is that “... purely
asymmetric heating generally leads to vortex weakening, usually in terms of the
symmetric energy, and always in terms of the low-level wind.” The present theory
shows that this conclusion does not hold up in case of a strongly tilted vortex,
but that in this case symmetric and suitably arranged asymmetric heating have
vortex intensification efficiencies of the same order of magnitude.


  * 
The theory compares favorably with three-dimensional numerical simulations based on
the compressible Euler equations.




To arrive at these results, we first recount the governing equations and the principles
of our analytical approach in section <ref>, and then revisit the
derivations by <cit.>.
A discussion of the scaling regime is given in section <ref> to
investigate the influence of the Coriolis effect (item (<ref>)), and the asymptotic vortex core expansion is carried out in section <ref> analytically supporting the physical interpretation of the asymmetric intensification mechanism given in item (<ref>).
In section <ref> we establish the kinetic energy balance of item (<ref>).
Section <ref> presents results of the theory in comparison with three-dimensional computational simulations to corroborate item (<ref>). Conclusions and an outlook are provided in
section <ref>.







§ DIMENSIONLESS GOVERNING EQUATIONS AND DISTINGUISHED LIMITS








 §.§ Governing equations



The dimensionless inviscid rotating compressible flow equations for an ideal gas with
constant specific heat capacities in the beta plane approximation form the basis
for the subsequent asymptotic analysis:




   t
     
        + ·
         
            + w   z
             
                + 1/^2 1/ρ   p
                 
                    + 1/(1 + βy)  k ×
                     
                        = 0  ,
                      

   wt
     
        + ·w
         
            + w  wz
             
                + 1/^21/ρ  pz
                 
                   
                     
                        = -1/^2 ,
                      


   ρt
     
        + ·ρ
         
            + w   ρz
             
                + ρ·   
                    + ρwz
                     
                        = 0  ,
                      

   Θt
     
        + ·Θ
         
            + w  Θz
             
               
                 
                   
                     
                        = Q_Θ   ,
                      

   
     
       
         
           
             
               
                 
                   
                      ρΘ    = p^1/γ



Here p,ρ,Θ,,w are rescaled pressure, density, potential temperature, and the
horizontal and vertical velocities, and γ is the specific heat ratio.

The three-dimensional
gradient is =  +  ∂/∂ z with the horizontal gradient
=  ∂/∂ x +  ∂/∂ y, the
zonal, meridional, and vertical coordinates (x,y,z), and the related unit vectors
(,,). Finally, t is the time variable and Q_Θ is a
diabatic source term.

Table <ref> lists general characteristics of the near-tropical atmosphere.
Together with the combined values in Table <ref> they form reference values for non-dimensionalization.
Let an asterisk denote dimensional quantities, then the unknowns and coordinates in eq:CompressibleFlowEquationsDimless are


    p = p^*/p  ,
      ρ = ρ^*/ρ  ,
      
    (, w) = (^* , w^*)/u  ,
      
    (, z) = (^*,z^*)/  ,
      
    t = t^*u/ .


Note that u/ is an estimate of the large-scale thermal wind shear, and =  x +  y is the horizontal coordinate vector.





In deriving the dimensionless equations eq:CompressibleFlowEquationsDimless
using the quantities from tables <ref> and
<ref> the Mach, internal wave Froude, and Rossby numbers, and the β-parameter


    [                     =  u/√( R T)          ≈ 3.4· 10^-2;                     =        u/N          ≈ 1.1· 10^-1 ] ,
        [                   =       u/f         ≈      13.3;                   =       β/f         ≈ 2.7·10^-3 ]


emerge naturally.
These are replaced with functions of a single small expansion parameter
≪ 1 through the distinguished limits


    = ^3/2  ,
       = /N  ,
       = 1/ f  ,
       = ^3 β  ,


in line with the multiscale asymptotic modelling framework of <cit.>.
Here (N, f, β) = 1 as → 0, with
concrete values


    N = 0.91 ,
      
    f = 0.75 ,
      β = 2.7


derived from eq:GeophysicalParameters for = ^2/3 = 0.105.
Replacing the characteristic numbers according to (<ref>) we get




   t
     
        + ·
         
            + w   z
             
                + 1/^3 1/ρ   p
                 
                    + (f + ^3βy)  k ×
                     
                        = 0  ,
                      

   wt
     
        + ·w
         
            + w  wz
             
                + 1/^31/ρ  pz
                 
                   
                     
                        = - 1/^3 ,
                      


   ρt
     
        + ·ρ
         
            + w   ρz
             
                + ρ·   
                    + ρwz
                     
                        = 0  ,
                      

   Θt
     
        + ·Θ
         
            + w  Θz
             
               
                 
                   
                     
                        = Q_Θ   ,
                      

   
     
       
         
           
             
               
                 
                   
                      ρΘ    = p^1/γ   .



Whereas
f and β appear explicitly in (<ref>), N characterizes the background stratification of potential temperature and will be invoked below where we define the initial conditions for the vortex flow.

Equations eq:AsymptoticFlowEquations will form the basis for the
subsequent asymptotic analysis for ≪ 1, although much of the expansions
will proceed in terms of the small parameter


    = √() .









§ SCALING REGIME FOR LARGE VORTEX ROSSBY NUMBER AND STRONG TILT








 §.§ Vortex core size, intensity, and evolution time scale



Vortex core sizes of 50 to 200 are typical
for tropical storms and hurricanes, and the storm/hurricane threshold lies at wind
speeds of 30 <cit.>.
With ^2 ≡∼ 1/10, ∼10, and
u∼10, these data correspond well with


    ∼/^2 ≈100 ,
      ∼u/≈33  ,
       p_ v∼^4p ,











for a characteristic vortex core size , a typical wind speed, and the associated
depression in the vortex core, respectively. Note that these scalings deviate from
those adopted by <cit.>, who considered systematically larger
radii of the order ∼/^3 needed for direct matching to a
quasi-geostrophic large scale outer flow. From their work we recall, however, that the
vortex core structure and tilt develop on a time scale t_ v that is by
1/^2 longer than the vortex core turnover time scale
t_ to =  /. Thus, in view of eq:CoreScalings, we will follow
the vortex core evolution on the time scale


    t_ v
    = t_ to/^2
    = 1/^2/^2/u
    = t/^3∼10 .


The scalings in eq:CoreScalings and eq:TimeScales include the regime of
“rapid intensification”, defined by NOAA's National Hurricane Center
[http://www.nhc.noaa.gov/aboutgloss.shtml] to denote maximum wind accelerations
of 30  kt∼15 in 24.

Also, the adopted scalings describe a vortex in the cyclostrophic
regime since


    /u^2^2/r = 1    whereas    /u_θ^2f = 1/u/u =  ,


, the Coriolis term is subordinate to the centripetal acceleration in the horizontal
momentum balance in this regime. Accordingly, the vortex Rossby number is large,


    _ v
    = u_ max/f_0
    = u_ max/u/
    = ^-2-1+2
    = 1/ .








 §.§ Co-moving coordinates for a strongly tilted vortex



Following <cit.>, we resolve the flow dynamics on the vortex
precession and core evolution time scale t_ v from eq:TimeScales. The
appropriate time coordinate is


    = ^3 t  .


For the core structure analysis we introduce vortex centered horizontal coordinates


    =
    1/^2((, z) + )


where (, z) is the horizontal position of the vortex centerline at height z
and  is the relative horizontal offset. With this scaling  resolves
the core scale  from eq:CoreScalings and the centerline covers
comparable distances. This justifies the notion of “strong tilt”.

In the sequel we use polar coordinates in horizontal planes, ,


    = x  + y     where    {[  x = cosθ  ;; y =  sinθ  ; ][ = cosθ - sinθ; = sinθ + cosθ ].


with  and  the radial and circumferential unit vectors, respectively.
The transformation rules for derivatives in these coordinates read



    =   ^2(
              + 1/θ)   ≡  ^2  ,
          
    z|_t,x,y    =   z|_,,θ
               -  · ,
          
    t|_x,y,z   =  δ^3 ( |_,θ,z        -  ·)  .



The horizontal velocity is decomposed into the vortex' motion plus the relative velocity,


    =  + (u_r  + u_θ )  .


For later reference, here are the centerline represented in the (, ) basis,


    =
    (Xcosθ   + Ysinθ)  
       + (- X sinθ + Y cosθ)    ,


and the Fourier expansion of functions of the circumferential angle,
θ,


    F(θ) = F_0 + ∑_n (F_n1cos(nθ) + F_n2sin(nθ))  .


Note that we have exchanged the roles of F_n1 and F_n2 relative to their use
in <cit.> as this will streamline the analysis of the orientation of a dipolar field characterized by F_1 = (F_11, F_12)^T relative to the tilt vector ∂ X /∂ z.






 §.§ Vortex core expansion scheme



The circumferential velocity is expanded as



  
     (t,,z;ε)           =   0(t, , z)
         + 1(t, , z)
         + 2(t, ,θ, z)
         +   ,
  
 
  
     (t,,z;ε)           =          
         + .    2(t, ,θ, z)
         +   .



, non-axisymmetry relative to the centerline is allowed for scaling orders from
u upwards. Across the core size length scale, ,
such asymmetries induce horizontal divergences of order
/ ∼u / ( / ^2) = ^3 u/,
see eq:CoreScalings. Since the flow field is anelastic to leading order
as derived below, this implies the vertical velocity scaling,


  
  w(t,,z;ε)
  = ^3 0(t, ,θ, z) + δ^3  .


Expansions for the thermodynamic variables are anticipated as follows,



  
            =        + ^2       + ^4  (4 + _4)
         + ^5  (5 + _5)
         + ^5  ,
  

  
      ρ       =        + ^2       + ^4  ( 4 + _4)
         + ^5  ( 5 + _5)
         + ^5  ,
  

  
      Θ       = Θ_0
         + ^2        + ^4  ( 4+ _4)
         + ^5  ( 5+ _5)
         + ^5  ,



<cit.>.
In eq:MesoExpansion, the variables
(p_0, p_2, ρ_0, ρ_2, Θ_2)(z) describe the stationary
background (Θ_0 is a constant), (_i, _i, _i)(, z), are higher-order
horizontal means, and
(i, i, i)(, , θ, z)
are the quantities of prime interest.

Note that, owing to the Fourier representation defined in eq:Fourier this
notational convention “overloads” the subscript (·)_0 with a double-meaning,
but the distinction should always be clear from the context.

The vortex centerline position is expanded as


    (, z) =  0(, z) + ^1 .









§ ASYMPTOTIC ANALYSIS OF THE CORE STRUCTURE EVOLUTION



This section revisits the analysis of <cit.> for large vortex Rossby
numbers focusing on the evolution equation for the primary circulation.






 §.§ Asymptotic equation hierarchy for the vortex core







The governing equations transformed to the co-moving coordinates are provided in appendix <ref>. Inserting the expansion scheme from the previous section
we obtain




- (u_θ)^2/
+ 1/ρ_0 4
      = 0 ,
            4θ
          = 0
        


- 2 u_θ1/
+1/ρ_0 5
- f_0 u_θ
      = 0 ,
            5θ
          = 0



from the horizontal momentum balance at leading and first order, respectively.
Each line in eq:HorMomLeadingFirst displays the respective radial balance
first and the circumferential balance as the second equation.
We observe from the radial component in eq:HormomLeading that the vortex
is in cyclostrophic balance to leading order which implies large vortex
Rossby number. The Coriolis effect enters only as a first-order perturbation
in the present regime as seen in the radial component of eq:HormomFirst.
The pressure perturbations p4, p5 inherit the assumed axisymmetry
of 0, 1 thanks to the leading and first order
circumferential momentum balances in eq:HormomLeading and eq:HormomFirst,
respectively.

The full second order horizontal momentum
equations are listed in appendix <ref>,
equations eq:HormomSecondApp, but for the rest of the paper we only need the
circumferential average of the circumferential component mom_theta_32.
Letting ψ_0 ≡1/2π∫_-π^πψ(θ)   dθ
denote the circumferential average of some θ-dependent variable ψ in
line with eq:Fourier, we have


    t
    + 0_0 z
    + u_r,02( 
    + /)
    - u_r,*2
    = 0  ,


where


    u_r,*2 = (0·0z)_0 .


The flow is hydrostatic up to third order, ,
p_iz = -ρ_i (i = 1,...,4), whereas


    ∂4/∂ z - ∂0/∂ z·4 = - 4  .


The leading and first order velocities are horizontal and axisymmetric according to
Meso:HorizontalVelocity, Meso:VerticalVelocity and thus divergence free.
The second order velocity is subject to an anelastic divergence constraint obtained
from the mass balance,


    ρ_0/(
      (2)
      + 2θ)
    + ∂/∂ z( ρ_0 0)
    - ∂0/∂ z·∇̂_ (ρ_0 0) = 0  .


Similarly, the first non-trivial potential temperature transport equation reads


    0/4θ
    + 0d Θ_2/d z
    = Q_Θ0  ,


and the equation of state relates the thermodynamic perturbation variables
through

    4
    =
    ρ_0
    (
      4/γ_0 - 4/Θ_0)     .








 §.§ Temporal evolution of the vortex structure



<cit.> observed that with the aid of eq:HorMomLeadingFirst and
eq:UrStarDef–state42, and given the vortex tilt, 0z,
as well as the diabatic source term, 0, one may interpret
eq:MomThetaTwoAverage as a closed evolution equation for the leading
circumferential velocity, 0.

To corroborate this, we use the Fourier decomposition, eq:Fourier, for 0
and  the representation of the centerline representation in polar coordinates from
eq:CenterlinePolar to obtain


    u_r,*2
    =
    (0·0z)_0
    =
    1/2[
      0_110z + 0_120z]  .


Expressions for 0_0 and 0_1k for k = 1,2
follow from the Θ–transport equation in potentialtemp2,


    0_0 Θ_2z =  Q_Θ,00  ,
        0_1kΘ_2z = Q_Θ,1k0 -
          (-1)^k u_θ/4_1[3-k] .


Since 4 is axisymmetric (see eq:HormomLeading),
4_1k≡ 0 and the equation of state, state42,
yields 4_1k/Θ_0 = - 4_1k/ρ_0.
With this information, the vertical momentum balance vertmom2 yields


    _114/Θ_0
      = - _114/ρ_0
      =  - 1/ρ_0∂0/∂ z∂4/∂  ,
        _124/Θ_0
      = - _124/ρ_0
      = - 1/ρ_0∂0/∂ z∂4/∂  .


Using the cyclostrophic balance in eq:HormomLeading to eliminate
4, and going back to w52 we obtain
expressions for the 0_1k in terms of
0, 0z, and 0,


    0_1kΘ_2z
      = Q_Θ,1k0 -
             (-1)^k Θ_0 X_[3-k]0z (u_θ^(0))^3/^2    (0pt10ptk = 1,2
    )  ,


where 0_1 ≡0 and 0_2 ≡0.
Upon insertion of this result in eq:UrStar, the second term on the right cancels,
so that


    u_r,*2
    =
    1/2  dΘ_2/dz[
      Q_Θ,1100z + Q_Θ,1200z]
    ≡1/2  dΘ_2/dzQ0_Θ,1·0z  .


Here we have interpreted the cosine and sine Fourier-1 components of Q_Θ0
as the components of a heating dipole vector, Q_Θ, in the horizontal plane.

To find a corresponding expression for u_r,02 (see the third term in eq:MomThetaTwoAverage), consider the circumferential average of mass continuity, mass2.
A brief calculation yields


    (ρ_0 u_r,02)
    + ∂(ρ_0 0_0 )/∂ z
    - 1/2[
    0z(ρ_0 0_11) + 0z(ρ_0 0_12)]
    = 0


or, equivalently,


    (ρ_0 [u_r,02 - u_r,*2])
    + ∂(ρ_0  0_0 )/∂ z
    = 0


with u_r,*2 defined in eq:UrStar. Exploiting w5211 in that definition
and integrating in  requiring that u_r,02 be finite at
=0 we find


    u_r,02 =
    
    
    
    
    u2_r,00 + u_r,*2  ,


where


    u2_r,00
    =
        - 1/∫_0^r/ρ_0z(ρ_0 Q_Θ,00/ dΘ_1/ d z)
          d r  .


With w52 (first equation), eq:UrStarFinal, eq:urzeroone, and eq:UradZeroZero
we have now indeed expressed w0_0, u_r,02, and u_r,*2 in
terms of 0, 0z, and Q_Θ0 as announced. In the
sequel, we may thus derive from eq:MomThetaTwoAverage
how vortex tilt and diabatic heating affect the evolution of the primary circulation.


The results in this section match the corresponding result by <cit.> with
the Coriolis parameter f_0 set to zero. This corroborates our statement (<ref>)
in the introduction that the vortex amplification/attenuation mechanism described
in their work does not depend on the vortex Rossby number being at most of order unity.







§ DISCUSSION OF THE ASYMMETRIC INTENSIFICATION/ATTENUATION MECHANISM








 §.§ The influence of asymmetric heating on the primary circulation



As elaborated in the previous section, eq:MomThetaTwoAverage describes the
evolution of the primary circulation in response to external diabatic heating in the
present vortex flow regime. Aiming to separate the influence of heating asymmetries from
those of axisymmetric effects, we recall from eq:urzeroone that the net
circumferentially averaged radial velocity is entirely a response to diabatic effects,
and that it consists of one part, u2_r,00, which, according to
eq:UradZeroZero is induced by axisymmetric heating, and a second part,
u_r,*2, which, according to eq:UrStarFinal, arises from first Fourier
mode asymmetric heating patterns. Using this decomposition in eq:MomThetaTwoAverage,
we rewrite the equation as


    t
    + 0_0 z
    + u_r,002( 
    + /)
    = - u_r,*2/ ,


which is the large-Rossby version of equation eq:TempevoluthetaIntro announced in the
introduction. In this equation, the left hand side captures the influence of the
axisymmetric dynamics and diabatic heating, whereas the right hand side covers all
effects due to the interaction of asymmetric heating and vortex tilt.






 §.§ Mechanics of vortex intensification by asymmetric heating of a tilted vortex





In the following section we analyze the leading-oder mass balance relations given in (<ref>) and (<ref>). We furthermore argue that u_r,*2 given in (<ref>) plays a crucial role in explaining the spin-up mechanism based on asymmetric diabatic heating. In this context we note that, according to (<ref>), the first order Fourier modes of the vertical velocity involve a contribution from diabatic heating (first term) and one due to the adiabatic dynamics (second term). It is only the contribution by diabatic heating that has an impact on u_r,*2 as seen in (<ref>).

Depicting the situation of asymmetric heating anti-parallel to the tilt in figure <ref> we observe that the suitably arranged vertical motions can generate a mass flux through the boundary of the tilted disc control volume with the coordinate interval (z, z+ Δ z). Considering mass continuity in centerline-attached coordinates in (<ref>), we can identify the term in brackets as the axisymmetric mean of the vertical mass flux. Equation (<ref>) reveals that this expression is equal to a horizontal mass flux governed by u_r,*2. We therefore conclude, that the net vertical outflow in figure <ref> is compensated in the present balanced vortex situation by a net horizontal inflow to preserve continuity.

This gives an additional spin-up mechanism which exploits the vertical (tilted) structure of the vortex to gain angular momentum by moving air masses from larger radii to the center of the vortex. In contrast, the opposite orientation of diabatic Fourier-1 modes leads to an attenuation of the vortex by transporting angular momentum away from the center.
We therefore claim that by this mechanism it is possible to influence the overall strength of an atmospheric vortex as will be demonstrated in section <ref>.

This should settle the announcement of (item <ref>) in the introduction.






 §.§ Energy budget for the externally heated vortex



Here we elaborate on how the asymmetric diabatic heating is transferred to kinetic
energy of the primary circulation in a tilted vortex.
This will be particularly useful in assessing the derived equations within the framework of Available Potential Energy <cit.>.

To this end we multiply
eq:MomThetaTwoAverageAsymmetriesSplitOff by ρ_0  0,
use the θ-averaged leading-order mass balance from anelastic and recast
the advective terms in conservation form to obtain,


    t(ρ_0^2/2)
    + (ρ_0 u_r,002^2/2)
    + z(ρ_0 _0 ^2/2)
    = -   u_r,02p4 .


Here we have dropped the 0 superscript on 0 and
w0 to simplify
the notation, and we have used the cyclostrophic radial momentum balance from
eq:HorMomLeadingFirst to introduce the pressure gradient on the right.


This reveals the change of kinetic energy (left hand side) to result from the work of
the pressure force due to the mean radial motion (right hand side). Some straightforward
but lengthy calculations, the details of which are given in
appendix <ref>, yield a direct relation of the kinetic
energy balance in eq:KineticEnergyBalanceI to the Lorenz' theory of generation
of available potential energy (APE) by diabatic heating,


    ()_t
    + ( u_r,002)_
    + ( w_00)_z
       = ρ_0/dΘ_2/dz1/Θ_0[Θ_04Q_Θ,00
         +  1/2Θ_14·Q_Θ,10] , 
       = rρ_0/N^2Θ_0^2(Θ4· Q_Θ0)_0


where =  + p4, and
(Θ,Q_Θ)_1
= (Θ, Q_Θ)_12 + (Θ, Q_Θ)_11
are the dipole vectors spanned by the first circumferential Fourier components of the
fourth order potential temperature perturbation, Θ4, and of the diabatic
heating function, Q_Θ0, respectively.

Equation (<ref>) poses the differential form of kinetic energy balance.
To end up with an integral form as presented in <cit.> we make use of the Gauss's theorem which allows us to drop the radial and vertical derivative assuming u_r,002 and w_02 vanish for sufficiently large  and z respectively.
To achieve this condition for u_r,002 (<ref>) shows that we do not only need Q_Θ,0 such that the integral converges for large radii but we need the integral to converge to zero.
When assuming a concentrated pattern of heat release with amplitude 10^-4 close to the vortex center (in the eyewall) over a surface of ∼ 10^4 km^2 it would need to be counteracted by contributions with opposite sign, , cooling, but over a much larger surface of ∼ 10^6 km^2.
This simple scale approximation reveals cooling rates of 0.1K/d which is by an order of magnitude smaller than what is observed by radiative cooling <cit.>.

For the total (integrated) kinetic energy E_k we find


    dE_k/dt = 2π∫_0^∞∫_0^∞ρ_0/N^2Θ_0^2( Θ4 Q_Θ0)_0   dr dz ,


On the one hand, <cit.> balanced the kinetic energy with the conversion rate from APE to kinetic energy (C) and the dissipation rate (D) where the latter is neglected here.
On the other hand, the expression on the right-hand side coincides with the generation rate (G) of APE (see appendix <ref> for details).
Therefore, no APE accumulates in the present flow regime as it is directly converted to kinetic energy (at leading order).
This is the result of the timescale used in the asymptotic analysis as conversion between APE and kinetic energy is accomplished by the advective and pressure-velocity fluxes on faster timescales.

In line with <cit.> and announced in the introduction in item <ref> this result shows that positive correlations of temperature perturbation and diabatic heat release lead to the increase of kinetic energy.

The precise form of the right hand
side of eq:KineticEnergyBudget as announced in the introduction
(item <ref>) is obtained from eq:EKinBalance by realizing that
(1/Θ_0) dΘ_2/dz is the dimensionless representation of N^2, the square
of the Brunt-Väisälä frequency, and that the constant Θ_0 is the
leading-order dimensionless background potential temperature
= T (Θ_0 + 1).


<cit.>, extending prior similar studies, investigate the influence
of asymmetric diabatic heating on vortex intensification on the basis of a linearized
anelastic model that includes a radially varying base state and baroclinic primary
circulation. Their central conclusions are that (i) asymmetric heating patterns
quite generally tend to attenuate a vortex, that (ii) there are situations in which
they can induce amplification, but in these cases their influence is (iii) generally
rather weak. In fact, they state in their section e: “... purely asymmetric heating
generally leads to vortex weakening, usually in terms of the symmetric energy, and
always in terms of the low-level wind.”. Equation eq:EKinBalance shows, in contrast, that purely
asymmetric heating in a tilted vortex can intensify or attenuate a vortex depending
on the arrangement of the heating pattern relative to the tilt, and that the efficiencies
of symmetric and asymmetric heating in generating kinetic energy are of the same order
in the asymptotics as claimed in (item <ref>) of the introduction.






 §.§ Diabatic forcing of the centerline motion



Together with the previous discussion we want to highlight some aspects of the effects of asymmetric heating on the vortex centerline motion.
Examinations in appendix <ref> of the constituents of the centerline equation reveal that the term k × splits into an adiabatic and a diabatic contribution due to the linear dependency of  on the vertical velocity dipole w_1 and w_1 being composed of an adiabatic and a diabatic contribution (see (<ref>)).
In particular, we find that the adiabatic expression is of the same form as the diabatic one but evaluated with the adiabatic vertical velocity:


    w_1,ad = - [  0  1; -1  0 ]1/Θ'_2 W ∂_z  X
        = R̂_-π/21/Θ_2' W ∂_z  X ,


with


    W = u_θ/r( u_θ^2/r + f u_θ) ,


and the rotation matrix R̂_θ_0.

Inserting (<ref>) into eq. (<ref>)  results in a linear differential operation on X which, interpreted as Hamiltonian of a (complex-valued) Schrödinger equation (<ref>), leads to a purely real spectrum, , to precession of the centerline in the complex (x-y) plane.

The diabatic motion of the centerline on the other hand results from inserting some non-trivial w_1,dia into  which in our case shall be a rotated version of (<ref>):


    w_1,dia = R̂_θ_01/Θ_2' W ∂_z  X .


Here θ_0 is the relative orientation of the diabatic vertical velocity dipole relative to the tilt.
Clearly, w_1,dia coincides with w_1,ad for θ_0=-π/2.
Therefore, such a diabatic heating vertical velocity pattern will result in an additional contribution to the centerline motion of the same orientation and magnitude as the adiabatic contribution.
By numerical experiments we observed that the adiabatic vertical velocity dipole, oriented -π/2 relative to the tilt leads to a centerline motion in the direction +π/2.

In contrast by varying θ_0 we expect a diabatic centerline forcing which is oriented -θ_0 relative to the tilt.
Therefore, for our formulation of diabatic heating the rotation angle θ_0 determines whether the diabatic forcing leads to acceleration (θ_0=-π/2) or deceleration (θ_0=π/2) of the centerline precession or to increasing (θ_0=π) or decreasing the tilt of the centerline.
[Note that the resulting linear operator governing the diabatic centerline motion effectively is the rotated version of the adiabatic operator (multiplied with complex number on the unit circle), , its spectrum is rotated in complex plane as w_1 is rotated.
Real parts of the spectrum lead to precession while imaginary parts lead to growth/damping of the centerline amplitude.]

In the further course of diabatic experiments we make use of our findings and specify heating patterns of the form


    Q_1,θ_0 := Θ_2'  w_1,dia = _θ_0 W ∂_z  X.








§ COMPARISON BETWEEN ASYMPTOTIC MODEL AND THREE-DIMENSIONAL SIMULATIONS



In the course of this paper we presented the derivation of a system of PDEs governing the leading-order dynamics of a tilted tropical cyclone in weak cyclostrophic regime (eqs. (<ref>) and (<ref>)).

The aim of this section is twofold:
First, we want to validate the reduced model equations against the three-dimensional Euler equations, , first principle equations (<ref>), by solving both sets of equations in a suitable numerical framework.
To this end, we follow the work of <cit.> who presented the first results on this scenario.
The second goal is to highlight principal mechanisms that are activated by purely asymmetric heating.
Therefore, after analyzing the adiabatic dynamics of an initially tilted vortex, we continue by constructing a prototypical asymmetric diabatic heating pattern which will be imposed on the vortex under different angles relative to the tilt.
We will refer to these experiments as adiabatic (reference simulation), stagnation, intensification and attenuation according to their influence on tilt and centerline.

The quasi-two-dimensional equations (<ref>) and (<ref>) are solved by a combination of appropriate numerical methods details of which are presented in section <ref>.
For the three-dimensional simulations the general purpose atmospheric flow model EULAG <cit.> provides efficient integration strategies for equations (<ref>), and its compressible model was used. These two alternative representations of the tilted vortex flows will be referred to as asymptotic and three-dimensional simulations.

Note that, although we worked out the scaling in the cyclostrophic regime (= 1/δ), initial data in the following sections will be in the gradient-wind regime (= 1).[Experiments have shown substantial damping of the centerline tilt amplitude when initializing the vortex in the cyclostrophic regime.
Furthermore, the special formulation of asymmetric diabatic heat release sensitively depends on the magnitude of tangential velocity causing numerical instabilities in the cyclostrophic regime.]
Therefore, in the numerical treatment terms involving the Coriolis parameter f, even if reasonably small, will not be neglected.



 §.§ Numerical setting and initial data



With the asymptotic analysis above we demonstrated that a tilted vortex evolves at leading-order on a time scale comparable to the synoptic time scale.
Higher-order dynamics occurs in the presence of initial perturbations (excitation of higher-order asymptotic modes) on faster time scales.
However, we are interested in the leading-order effects only, hence we construct initial data to closely reproduce the leading-order symmetries imposed by the asymptotic analysis of section <ref> allowing to solve for solutions on the slowest varying manifold.

In the case of an adiabatic vortex the tangential velocity equation (<ref>) is trivially stationary,


    t u_θ0=0 ,


and the centerline equation (<ref>) becomes a Schrödinger-type equation,


    i ∂_t X = Ĥ X ,


with a Hamiltonian Ĥ depending on the tangential velocity u_θ(r,z) and the z-dependent background profiles ρ̅, p̅, Θ̅.
Note, that we have introduced X := ( X · i) +  i ( X · j) as the complex representation of X
(see appendix <ref>).

With suitable boundary conditions the governing Hamiltonian takes the form of a Sturm-Liouville operator and therefore exhibits a real, discrete spectrum <cit.>, which sets the precession frequency of each eigenmode.
The first non-trivial eigenvalue corresponds to the slowest varying solution and posses a cosine-like eigenfunction, , the simplest tilted solution.[Details on that are skipped but can easily be validated by numerical means.]

As tilt is crucial for coupling asymmetric diabatic heating modes to the leading-order vortex dynamics, we prescribe an initially barotropic tangential velocity profile u_θ(r), corresponding to a Gaussian vorticity profile, q(r),


    = q_m 1 - e^-σ^2 r^2/2 σ^2 r ,


where the radial vorticity profile reads


    q(r) = q_m e^-σ^2 r^2,


and choose the first non-trivial eigenfunction of the corresponding Hamiltonian to define the initial centerline geometry, scaled to a displacement of 160 km, see figure <ref>.



For the sake of complying with horizontal boundary conditions imposed by EULAG and to avoid reflective features near the boundary of the computational domain, the initial data for the three-dimensional simulations are smoothly transitioned to zero at some finite radius by applying a mollifier:


    m(r) =
      sin^2( π/2r/r_0)    ,    r < r_0
    
      1     ,   r_0 < r < r_1 
    cos^2( π/2r-r_1/_̊1-r_∞)    
        ,   r_1 < r < r_∞
    
      0     ,   r > r_∞


r_1 = 1 250  and r_∞ = 1 750  are the radii where the mollifier starts and where it reaches full-suppression. In addition to that the profile within r_0 = 100  is adjusted to preserve differentiability at the origin.

As mentioned earlier, the solution depends on the background state of the atmosphere which is determined by the potential temperature profile


    Θ̅(z) = T exp(N^2/g  z) .


In the three-dimensional case the vortex is embedded into a domain of 4000  extent in
both horizontal directions, 10  in the vertical, and a damping layer surrounds the domain
near the horizontal boundaries to suppress gravity waves emerging from the inner core and to keep them from reflecting.

Along the structural properties of the tropical cyclone, i.e., an inner core and smooth transition to the quasi-geostrophic far field, we use the static mesh refinement capability of EULAG <cit.> and map equidistant coordinates onto a grid focused at the inner core.
The actual mapping is accomplished by


    [ x_p; y_p ]
     = c_1
      [ x_c; y_c ]
      +
      c_2
      [ x_c^α; y_c^α ] ,


where α = 5, c_1,2 = 1/2.
(x,y)_p,c are normalized coordinates on the domain [-1,1]^2.
Figure <ref> demonstrates how the horizontal grid is focused towards the center of the computational domain.




The asymptotic equations are solved on a regular equidistant tilted polar grid on the domain (r,z) ∈ [0, 1.12] × [0, 12.5] in dimensionless units covering roughly a tilted cylinder of 1000   around the centerline and the full vertical extent.

In the further course, we will compute the diabatic heating from eq. (<ref>) which involves the reconstruction of both, the current centerline position and circumferentially averaged tangential velocity.
For the first, we compute the center of mass of the vorticity field at each horizontal level:


    X = ∫ k · (∇× u_||)  x   dx dy


We follow the ideas of <cit.> and the implementation outlined in <cit.>.

The circumferential mean of u_θ then is computed by the Biot-Savart integral


    u_θ,0 = Γ(r,z)/2π r = 1/2 π r∫_B_r( X) k · (∇× u_||)   dx dy ,


where B_r( X) = {( x_||∈ℝ^2 | ( x_|| -  X)^2 < r^2 } denotes the circular domain centered at X with radius r.



 §.§ Results and discussion



In the following subsection we will present results of numerical simulations solving either the full three-dimensional Euler equations (<ref>) or the reduced asymptotic equation (<ref>) for the primary circulation velocity u_θ0, and the centerline evolution (<ref>) explained in detail in App. <ref> and following.



  §.§.§ Adiabatic vortex



As a reference for the following experiments we first investigate the dynamics of a tilted, adiabatic vortex.
In subsection <ref> we constructed initial data to follow the first non-trivial eigenmode of the governing (adiabatic) Hamiltonian and we found stationarity (eq. (<ref>)) of the mean tangential wind.
Hence, from the structure of eq. (<ref>) we expect undamped precession of the eigenfunction.

Figure <ref> compares the results of both adiabatic simulations, three-dimensional and asymptotic, with initial data as discussed.
Though exhibiting small-scale oscillation and damping, the three-dimensional simulation (fig. <ref>, right panel) compares well with its asymptotic analog.
Time scales of one precession are 5.5 days for the asymptotic and 6.5 days for the three-dimensional simulation.
This difference leads to a deviation of the final positions of the centerline but given that the effective expansion parameter δ = √()∼ 1/3, this is well within the error bounds of the leading order solution.















The asymptotic analysis revealed non-trivial leading-order balances for w and the thermodynamic quantities according to tilt, gradient-wind (cyclostrophic) and hydrostatic balance.
In figure <ref> both, w and Θ̃ are visualized by representative horizontal slice at 5000  at t = 6.5 comparing asymptotic values (black contours) with the 3D numerical simulation results (color-coded).
The tilt vector ∂_z  X is indicated by an arrow.
Qualitative similarities are rather apparent, while deviations are again well within the asymptotic truncation order δ∼ 1/3.
Both figures demonstrate the alignment of dipolar perturbations relative to the tilt as w_1 and Θ'_1 are rotated by -90  and 180, respectively.
This is in agreement with the findings of <cit.>.






  §.§.§ Stagnation



The stagnation test follows the idea that the choice of θ_0 = π/2 in (<ref>) leads to deceleration of the centerline precession by canceling the term  in (<ref>).

Furthermore, it has no impact on the leading-order tangential velocity as we immediately see when neglecting symmetric vertical motions due to symmetric diabatic heat release resulting in the tangential velocity time evolution


    u_θt = -u_r,*(u_θ/r + f )


and realizing that u_r.* is the projection of the diabatic heat release onto the tilt vector,


    u_r,* = 1/2( ∂_z  X · Q_1,π/2)
              = 1/2( ∂_z  X · W R_π/2∂_z  X ) ,


which vanishes due to orthogonality.

By construction, inserting the heating Q_1, π/2 into (<ref>) satisfies w_1 ≡ 0 up to leading order.
Figure <ref> shows the expected behavior for the three-dimensional simulations.
Q_1, π/2 attenuates vertical velocity by a factor of ∼δ, , canceling the leading-order expansion mode of w.
The residual, depicted in the right panel may refer to higher-order expansion modes.



Furthermore, the leading-order tangential velocity is not affected by asymmetric heating of that orientation.
Figure <ref> presents both the time series of maximum tangential wind (blue) and of maximum heating (red).
Only small variations of the tangential velocity are apparent and, as we will see in the next subsection, this changes substantially when we alter the orientation of Q_1, θ_0.



Experiments (not shown here) revealed instabilities caused by small perturbations due to discretization errors:
The diabatic heating will exhibit large amplitudes where the local tilt ∂_z  X (being result of the reconstruction (<ref>)) is large.
This affects the local velocity and as a consequence the spectral properties of the Hamiltonian of eq. (<ref>) (projecting onto higher frequency modes) increasing small-scale oscillatory features of the centerline X.
Hence, we need make sure to maintain a certain regularity of X to avoid obscuring the effect under consideration by triggering this feedback loop.
We achieve this by restricting the heating to a concentrated pulse by applying a time-dependent amplitude factor of the shape


    f(t) = a exp( (t-b)^2/2c^2)


For the current setting a=1, b=5.5 and c=1 /2√(2 ln 300).

The heating distribution Q_1,π/2 constructed in such a way satisfies u_r,* to vanish (cf. (<ref>)) and additionally, by canceling w_1,  is canceled in the centerline equation of motion (<ref>). As discussed in the introduction of section <ref>, we are initially in the regime of gradient-wind balance.
Hence M_1 remains the only contribution to the centerline equation of motion.
Numerical evaluations show that the amplitude of M_1 is about 1/6 of the amplitude of the  in the adiabatic case which is why we observe a significant slow down of the centerline precession in the asymptotic case in figure <ref>, left panel, during the heating time interval (cf. fig <ref> for reference).
Although not as prominently as with the asymptotic simulation due to the aforementioned discretization errors, the centerline also slows down in the three-dimensional simulation.
I becomes also present that due to heating the shape of the centerline becomes superimposed by higher-frequent features that are excited due to the slight misalignment of the diabatic heating dipole.





  §.§.§ Intensification



In subsection <ref> we constructed a dipolar heating pattern and by aligning its orientation π/2 relative to the centerline tilt we found stagnation of the centerline precession as well as suppression of the vortex-induced vertical velocity dipole.
However, what is the influence of this heating dipole when oriented with θ_0 = π relative to the tilt?
θ_0=π determines the sign of (<ref>) negative and as a consequence the right-hand side of (<ref>) positive, , leading to intensification.
To examine the effect of the orientation of a dipolar heating pattern without altering its amplitude by the dependency of Q_1,θ_0 on u_θ and ∂_z  X we fix u_θ and | Xz| in (<ref>) to values at t=5 but still keep track of the orientation of Xz.

In contrast to figure <ref>, where we did not see any sizeable impact of diabatic heating on the mean circumferential velocity, figure <ref> displays a clear increase for both the asymptotic and three-dimensional equations.
We notice that the intensification is more effective in the asymptotic simulation (∼ 3 m/s) compared to the three-dimensional simulation (∼ 1.5 m/s), even though the heating amplitude of the asymptotic simulation is tuned to meet the three-dimensional simulation.
This may be caused by the fact that the centerline in the three-dimensional simulations is a reconstruction from the flow field and hence affected by inaccuracies leading to non-optimal alignment of the heating dipole.
Furthermore, even with rather high resolution the three-dimensional simulation is affected by numerical damping leading to a decrease of the centerline tilt over time which further reduces the impact (cf. eqs. (<ref>) and (<ref>)) even with Q_1,θ_0 computed by (<ref>) with values of ∂_z  X just before the heating.
Figure <ref> reveals the tilt dynamics without and with diabatic heating for both simulation approaches. In the three-dimensional simulations the tilt is more distorted than in the asymptotic simulations and the effect on the centerline tilt is less pronounced.
However, the overall behavior is comparable between the two simulations.




From the discussion in subsection <ref> we concluded that a heating pattern with θ_0 = π would lead to an increase of the centerline tilt.
This behavior could be validated in the current situation as seen in figure <ref> for both simulations, asymptotic and three-dimensional.
As the increase in tangential velocity is not as efficient for the three-dimensional simulation as it is for the asymptotic one, after the period of heating, the centerline precesses with higher angular frequency in the asymptotic case.
Again, this may be due to inaccuracies in the three-dimensional simulation when computing and aligning the heating dipole.




In section <ref> we discussed the intensification mechanism as result of the effective vertical velocity dipole resulting from the superposition of both adiabatic and diabatic contributions.
In the current situation the diabatic vertical velocity dipole is oriented by θ_0=π while the adiabatic dipole is oriented by -π/2.
As both have comparably the same amplitude we would expect -3π/4 rotation between tilt and resulting vertical velocity which is verified by results shown in the right panel of figure <ref>.




Although the intensification is rather weak, there is evidence that the orientation of the asymmetry of the diabatic heat release matters for the evolution of the circumferential velocity.
In fact, by allowing for stronger heating by increasing the duration of the heat pulse the effect on the tangential velocity is stronger (see figure <ref>), but it also affects the structure of the centerline in a more profound manner.




Figure <ref> also demonstrates that the system's kinetic energy increases as a consequence of the imposed asymmetric heating in line with our findings of section <ref>.






 §.§ Attenuation



The final experiment of this work consists in switching the heating dipole pattern to a configuration where we expect attenuation of the vortex and vertical alignment of the centerline.
Following equation (<ref>) attenuation corresponds to θ_0=0, , positive heating in the direction of the centerline tilt.
Again, tilt amplitude and tangential velocity of eq. (<ref>) are set to initial values to avoid non-linear feedback and restricted by an amplitude factor to act over a short interval only.

In figure <ref> for both, asymptotic and three-dimensional simulations, the centerline aligns when forced by heating.



In addition, figure <ref> demonstrates the reduction of integrated kinetic energy due to the attenuating heating dipole.






 §.§ Summary of tilt dynamics


We want to emphasize the effect of asymmetric diabatic heating on the centerline tilt by analyzing the time series of tilt amplitude as measured by an L_2-norm:


    ‖ Xz‖ := ∫_0^z_top( √( Xz^2))   dz,      z_top = 10 ,


With figure <ref> this quantity is plotted for all types of experiments performed in the course of this work.
It confirms on the one hand, that the orientation of the diabatic heating dipole correlates with changes in the tilt amplitude.
Aligning the heating dipole with the tilt (attenuation case) leads to decreasing the tilt, , the vortex aligns.
This situation turns around if the heating dipole is rotated 180 (intensification case) in that the tilt further increases.
The stagnation configuration (tilt and heating dipole are perpendicular) leads to no significant alterations from the adiabatic behavior for the asymptotic simulations.

On the other hand, we see that the three-dimensional representation of the initial data does not project as well onto the first eigenmode of the centerline as in the asymptotic cases.
Besides features which oscillate on a time scale of less than one day all three-dimensional simulations exhibit an additional oscillation on the scale of roughly six days.
However, in every case, asymptotic and three-dimensional, the effects of activating an asymmetric diabatic heating dipole are superimposed onto the adiabatic reference.
The stagnation simulation follows the adiabatic reference on both cases but exhibits slight distortions in the three-dimensional case.
As discussed before this is the reason for restricting the diabatic heating to a short time period as this effect would increase for longer periods.
However, for both experiments, intensification and attenuation, the tilt amplitude follows approximately the adiabatic reference curve on an increased and lowered level, respectively.

We argue that the response of asymmetric diabatic heating for the three-dimensional simulations is not as direct as it is for the asymptotic analog simulations.
This is probably due to imbalances excited through the diabatic heating, , by misalignment of the heating dipole, and limits the efficacy of the heating in (<ref>) in influencing tilt and circumferential velocity.
Nonetheless, in all cases the intensification configuration increases the centerline tilt while the attenuation configuration decreases it.




The analysis above corroborates our initial statement that the orientation of a purely dipolar diabatic heating pattern intensifies and shears a atmospheric vortex apart in the anti-parallel orientation of heating dipole and tilt while it attenuates and aligns the vortex vertically in the parallel orientation.











§ CONCLUSIONS AND OUTLOOK



With the present work we have extended the results of <cit.> to vortex Rossby numbers larger than unity corresponding to wind speeds of 𝒪(30)  and showed that the principal structure of the reduced asymptotic model equations is retained in this limit. We found that the validity of the equations holds from the Gradient-wind regime up to (modest) cyclostrophic vortices. This corresponds to hurricanes of strength H1 on the Saffir-Simpson scale.

Current models of tropical cyclone intensification rely on organized symmetric heating in an upright vortex <cit.>.
Observations show that in incipient tropical storms the level of organization of convection is weak compared to that of mature storms.
Our findings of leading-order effects of asymmetric diabatic heating both on the strength of the primary circulation and on the vortex tilt may add a new route of acceleration from tropical storms to mature hurricanes.

Here, we have focused on highlighting the potential effect of asymmetric diabatic heating release on the mean tangential velocity, which we argued, can be of the same order of magnitude as symmetric diabatic heat release.
Although the gain of horizontal wind was limited as we restricted to heating configurations which maintained the overall flow structure which the asymptotic analysis was based on, we expect to see potentially stronger efficiencies in nature due to the self-regulation of moist convection in a sheared environment.
Moist thermodynamics has been replaced here, however, by artificial diabatic source and sink terms neglecting effects of water phase transitions. As we argued, the resulting asymmetric pattern of vertical velocity is the driver for intensification/attenuation.
In ongoing research, to be published elsewhere in the future, we find that mean vertical mass fluxes of an ensemble of convective towers can provide for similar effects.
Thus, despite its somewhat artificial setup, the present study does reveal an interesting physical mechanism.

In this work we did not discuss the interaction of environmental shear with the TC. Observations of TCs just before rapid intensification often show a phase of relatively stationary configuration of background wind shear <cit.>, vortex tilt, and asymmetries of convection <cit.>.
Further, we explicitly restricted to initially barotropic velocity distribution to avoid interaction due to baroclinicity.
Still, these interactions may find an explanation in the framework of the present asymptotic model and are subject to current investigation.










T.D., A.P., and R.K.'s work has been funded by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114 “Scaling Cascades in Complex Systems”, Project Number 235221301, Project (C06) “Multi-scale structure of atmospheric vortices”, and by the Helmholtz Society of Research Institutions for funding through the “GeoSim” Graduate College.
The authors also thank the European Centre for Medium Range Weather Forecast for supporting this work by R.K's ECMWF research fellowship as well as hosting T.D., A.P., and R.K. for research stays.
The authors acknowledge the North-German Supercomputing Alliance (HLRN) as well as the German Climate Computing Center (DKRZ) for providing HPC resources that have contributed to the research results reported in this paper.
The authors also gratefully thank Olivier Pauluis, Remi Taullieux, and Mike Montgomery for many fruitful discussions which have helped strengthen our interpretation of the asymptotic results.
T.D. and R.K. further thank Sundararaman Gopalakrishnan, Frank Marks, Paul Reasor, and Dave Nolan for their hospitality and insightful discussions during a research stay at NOAA/HRD and University of Miami.
NCAR is sponsored by the National Science Foundation.








§ GOVERNING EQUATIONS IN THE CO-MOVING COORDINATES



Transforming (<ref>) to the vortex-centered
coordinates from section <ref>
using eq:DerivativeTransformations and defining
≡ t and = u_r  + u_θ
we find


  
  
    (+)
        + 1/·
        + /^3[z
                  - ·] (+)
           
    
        + 1/^71/ρ p
        + 1/^2 f × (+)    =
        0   ,
        
    w
        + 1/·w
        + /^3[z
                - ·]    
    
        + 1/^10 (1/ρ [z
                  - ·]
        + 1
        )
            
        = 0  ,
        
    ρ
        + 1/·(ρ)
        + 1/^3[z
            - ·] (ρ)
            =
        0  ,
        
    Θ
        + 1/·Θ
        + /^3[z
                   - ·] Θ   =
        Q_Θ  .










§ FULL SECOND ORDER HORIZONTAL MOMENTUM BALANCES






    /2θ
    - 2 2/
    - (1)^2/
    +  ·0z0/
    
    + 1/ρ_06
    -  ρ_1/ρ_0^24
    
    
    - f_0 1
    = 0




    t
    + 0z
    + 2(  + /)
    + /2θ
    
    - 0·0z
    +  1/ρ_0 6θ
    
    
    = 0










§ DERIVATION OF THE KINETIC ENERGY BUDGET EQ:EKINBALANCE



We start from eq:MomThetaTwoAverageAsymmetriesSplitOff, which is equivalent to
(4.21) in <cit.> for f_0 = 0. This is verified straightforwardly
by using u_r,02 = u_r,002 + u_r,*2. The equation is
multiplied by ρ_0 0, and we use the mass conservation law
in the form of eq. anelastic, ,


    (ρ_0 u_r,002)_
    + (ρ_0 w_00)_z = 0 ,


to generate the advective transport terms of kinetic energy in conservation form.
We let


    = ρ_0 0^2 / 2


and obtain


    ( e_k)_t
    + ( u_r,002)_
    + ( w_00)_z
    = - (u_r,002 + u_r,*) p4 .


Focusing on the right-hand side of this equation, we rewrite the first term as


    [                                                    u_r,002p4                                                            = ( u_r,002p4)_
          + ( w_00p4)_z
          - ρ_0 w_00(p4/ρ_0)_z;                                                                                                                                     - p4/ρ_0[(ρ_0 u_r,002)_
              + (ρ_0 w_00)_] . ]


The square bracket vanishes according to eq:MassConservationApp, while we observe
that by combining the axisymmetric part of the hydrostatic balance in
vertmom2 with the equation of state in state42 to replace
(p4/ρ_0)_z, and w52 to replace w_00 one finds


    ρ_0 w_00(p4/ρ_0)_z
    = ρ_0 Q_Θ,00/dΘ_2/dz4_0/Θ_0 .


Next we rewrite the second term on the right of eq:EKinBalanceAppI
using the definition of u_r,*2 in eq:UrStar, and the
first Fourier modes of the vertical momentum balance in vert11newtest2
to find


    u_r,*2p4
    = 1/2ρ_0 w_1·Θ_14/Θ_0
    = 1/2ρ_0 Q_Θ,10/dΘ_2/dz·Θ_14/Θ_0 .


To obtain the second equality we have used the asymmetric WTG-law from
w5211 and the fact that the second term in that equation contributes
a component to w_1 that is orthogonal ∂X0/∂ z,
and thus also orthogonal to Θ_14.

Insertion of eq:PressureWorkTerms–eq:AsymmetricHeatingPart generates
the desired equation eq:EKinBalance.








§ RELATION TO LORENZ' THEORY OF AVAILABLE POTENTIAL ENERGY



We want to strengthen the fact that our outlined theory is compliant with Lorenz' concept of available potential energy (APE) <cit.>.
We will see, under certain assumptions, that the kinetic energy generation of equation (<ref>) is identical to Lorenz' expressions of APE generation and conversion to kinetic energy.

Therefore, we start with his equations (16), (20), (17), (18):


    At   = -C + G  , 
     Kt   = C       , 
    
     C                 = -R/g∫_0^1/pTω  dp 
    
     G                 =  1/g∫_0^Γ_d/Γ_d - ΓT'Q'/ T  dp  .


Here, A is the average available potential energy, K the mean kinetic energy, C the conversion of APE to kinetic energy, G the generation of APE (all per unit surface),  the surface pressure ω = pt = pz w = -g ρ w the vertical velocity in pressure coordinates, Γ=  Tz the  lapse rate of temperature, Γ_d=g/ the dry-adiabatic lapse rate, Q = Θ Q_Θ the heat rate, (·) the horizontal mean, and (·) the deviation from the mean.
In our case we neglect friction, hence equation (21) of <cit.> is trivially zero.

We can make use of the results from above, as well as  those from <cit.>, that the leading-order expansion modes of the thermodynamical quantities are horizontally homogeneous. For example, for potential temperature we have


    Θ   = T( Θ_0 + δ^2Θ_2 + δ^4) 
    Θ   = T( δ^4Θ4 + δ^5)


Furthermore, we can find derived expression:


    T = Θπ
     T/ T = Θ/Θ + π/π


Additionally, we make use of the circumstance, that the mean of an asymmetric field is zero.

With that we get for equation (<ref>)


    G = g∫_0^H ρ/ΘΘz(Q_ΘΘ)  dz ,


which is the the right-hand side of equation (<ref>) in physical dimensions, integrated over the whole domain, where Q_Θ is decaying sufficiently fast.

With (<ref>) we get for the resulting vertical velocity


    w = 1/Θz( Q_Θ - u_θ/rΘ'_1 · e_r )


in physical units.
With that we get for the horizontal mean ( T ω)


    ( T ω) = -g/ R p/ΘΘz( Θ Q_Θ) .


Inserted into (<ref>) we get for C the same expression as for G in (<ref>).

Finally, we conclude, that the generation of kinetic energy is equivalent to the generation of APE in Lorenz' theory and furthermore, that the expressions of generation APE and conversion into kinetic energy are identical for the setup of an asymmetrically heated vortex.
All the APE, which is generated by heating is converted immediately into kinetic energy.







§ THE CENTERLINE EQUATION OF MOTION



As preparation for appendix <ref> give more details on the centerline equation of motion (<ref>) and provide all the necessary information to close the system of equations (in conjunction with (<ref>)).

We will present the missing terms necessary to close equation (<ref>) in section <ref>, provide a split into adiabatic and diabatic contributions in section <ref> and further deepen the analysis of this equation in section <ref> for being able to construct a stable and efficient numerical scheme to solve the asymptotic equations.



 §.§ Formulation of Paeschke et al. (2012)



In the course of this work we argued that the equations derived by <cit.> are valid in the vanishing-Coriolis case, , f_0 → 0.
The structure of the equations stays essentially the same, only terms proportional to f_0 are to be dropped.

For completeness below we present the remaining expressions for M_1 and  necessary to solve the system (<ref>) and (<ref>):


    M_1 = f^2/4πρ_0 Γ∂_z
        ( ρ_0 Γ^2/Θ_2'∂_z  X )


= L[] = _π/2 L[_1] depends the following expressions:


  
    _1      = _1 + _1 + _1 + _1 , 
    _1      = ∂_r ( r  w_1 ∂_z u_θ) , 
    _1     = _1 + H_s(r-1) 1/r^2 I_1 , 
    _1      = r( ζ + f )  W_1 , 
     I_1    = Γ/2π_-π/2 M_1 , 
    _1      = (∂_r ϕ_1) (r∂_r ζ) , 
    _1      = (w_0 u/r - ∂_r(r w_0 ∂_r u_θ)) ∂_z  X ,

  


where L[·] is an integral operator,


    L[] = π/Γ∫_0^∞ r (r)  dr ,


and _θ_0 denotes the matrix of two-dimensional rotation by an angle θ_0.

Equations (<ref>) – (<ref>) involve terms which are resolved in terms of u_θ, X, Q_0 and Q_1:


    W_1       = -1/ρ_0∂_z (ρ_0  w_1) 
    
      w_0             = Q_0/Θ'_2
     w_1       = 1/Θ_2'(  Q_1 + W _-π/2∂_z  X ) 
    
      W               = u_θ/r( u_θ^2/r + fu_θ) 
    ζ   = 1/r(r u_θ)r
    ϕ_1    = -r ∫_r^∞1/r̅^3∫_0^r̅r̅̅̅^2  R_1   dr̅̅̅  dr̅
     R_1       =  W_1 + 1/2(∂_r w_0) ∂_z  X


The expressions above give rise to the reformulated equation of centerline motion:


    ∂_t  X =  u_s +
        1/2ln (δ)  k × M_1 - L[_1]




 §.§ Split into diabatic and adiabatic contributions



To reveal more clearly the structure of the centerline evolution equation,
we recall the formula w5211 for the vertical velocity Fourier modes,
which separates diabatic from adiabatic effects. Rewriting this formula
in the dipole vector notation we have


    w_1           =  w_1, + w_1, , 
     w_1,   = 1/Θ_2' Q_1 , 
     w_1,   = 1/Θ_2' W _-π/2∂_z  X
        =: _-π/2ŵ X ,


and realize that the adiabatic part is a linear (differential) operation on X:


    ŵ = W/Θ_1'∂_z


From the equations in the previous subsection we see that by linearity of the expressions in w_1 we can assemble  (and ultimately ) by linear superpositions of linear operations on X (operators are symbolically denoted by ) and in general nonlinear diabatic expressions.

_1 then becomes:


    _1
           = _1, + _-π/2ℋ̂ X 
    _1,   = ∂_r ( r  w_1,) 
    ℋ̂ X
           = ∂_r ( (ŵ X) ∂_z u_θ)


We first need to evaluate the expression for W_1,


    W_1
           =  W_1, + _-π/2Ŵ X
     W_1,   = -1/ρ_0∂_z (ρ_0  w_1,) 
     W_1,   = -_-π/21/ρ_0∂_z (ρ_0 ŵ X) := _-π/2Ŵ_1, X


to split _1 accordingly:


    _1    = _1, + _-π/2ℐ̂ X 
    _1,   = r( ζ + f )  W_1,
    ℐ̂ X    = r( ζ + f ) Ŵ_1, X


Together with M_1


    M_1    = f^2/4πρ_0 Γ∂_z
        ( ρ_0 Γ^2/Θ_1'∂_z  X ) 
       =:M̂ X


we get the following split for _1:


    _1    = _1, +
        _-π/2ℐ̂ X +
          H_s(r-1) 1/r^2Γ/2π_-π/2M̂ X 
       =: _1, + _-π/2ℐ̂̃̂ X


Performing the split of R_1,


    R_1    =
         W_1, + _-π/2Ŵ X + R_Q,0∂_z  X  , 
    
      R_Q,0   = 1/2(∂_r w_0)


ϕ_1 divides into


    ϕ_1    =
        ϕ_1, + _-π/2ϕ̂ X + ϕ_Q,0∂_z X 
    ϕ_1,   =
        -r ∫_r^∞1/r̅^3∫_0^r̅r̅̅̅^2
           W_1,  dr̅̅̅  dr̅
    ϕ̂ X    =
        -r ∫_r^∞1/r̅^3∫_0^r̅r̅̅̅^2
          Ŵ X   dr̅̅̅  dr̅
    ϕ_Q,0   =
        -r ∫_r^∞1/r̅^3∫_0^r̅r̅̅̅^2
          R_Q,0  dr̅̅̅  dr̅


and therefore _1 into


    _1    =
        _1, + _-π/2𝒥̂ X + 𝒥_Q,0∂_z  X . 
    _1,   = ∂_r (ϕ_1,)(r∂_r ζ) 
    𝒥̂ X    = ∂_r ( ϕ̂ X)(r∂_r ζ) 
    𝒥_Q,0   = (∂_r ϕ_Q,0) (r∂_r ζ)


For Q_1 we find the simple shorthand


    _1    = ( w_0 u/r - ∂_r(r w_0 ∂_r u) ) ∂_z  X 
       =: 𝒬_0 ∂_z  X .











Our final result is identifying three different contributions to , a (generally nonlinear) diabatic term, an advective term, and a linear Sturm-Liouville-type operator acting on X:


    = L[_1,] + L[𝒥_Q,0∂_z  X+  𝒬_0 ∂_z  X] + R_-π/2 L[ℋ̂ + ℐ̂̃̂ + 𝒥̂]  X





 §.§ Characteristic structure of the vortex centerline equation



We continue rephrasing the original centerline tendency equation to further emphasize its structure.
By trivially identifying ℝ^2 with ℂ we symbolically transform two-dimensional vectors a = (a_x, a_y) ∈ℝ to a = a_x + i a_y ∈ℂ.
Operations such as (k ×·) and _π/2 become multiplications with i.
By that we can identify a substructure of the equation to be of Schrödinger-type, and advective contribution and sources, generally dependent on X, u_θ, and the coordinates (r,z,t) but without any further specification:


    i(∂_t X + L[(𝒥_Q,0+  𝒬_0) ∂_z X])    = - 1/2lnδM̂ X - L[ℋ̂ + ℐ̂̃̂ + 𝒥̂] X  + 
        +  i u_s - i L[ℋ_Q,1 +ℐ_Q,1 + 𝒥_Q,1]


By identifying structural components the centerline equation takes the form


    i(∂_t X + A ∂_z X)    = Ĥ X + i Q + i u_s


Note that the left-hand sides of eq. (<ref>) takes the form of an advection operation while the right-hand side involve linear and non-linear source terms.

<cit.> pointed out that the adiabatic time evolution of the vortex centerline poses an eigenproblem.
For that case the (homogeneous) centerline equation in the complex plane is written as


    i X_ht = Ĥ X_h


As the spectrum ω_k of Ĥ is real, (<ref>) can be interpreted as a Schrödinger-type equation, hence eigenmodes X_k precess with the angular frequency ω_k in the complex plane.
For the adiabatic problem we therefore find, that


    _k(t) = _ω_k t_k(t=0) .


For the numerical experiments presented in section <ref> we use the first non-trivial eigenmode for initialization corresponding to the non-zero eigenvalue with the smallest magnitude.

















































































































































§ NUMERICAL SCHEME FOR ASYMPTOTIC EQUATIONS



By providing a closure (Q_0,  Q_1) = F( X, u_θ, t) equations (<ref>) and (<ref>) form a closed set of partial differential equations which in general cannot be solved analytically.
Thus, we further analyze the structure of these equations seeking for an adapted numerical method to allow for efficient and stable time integration.

<cit.> first presented a numerical scheme for solving the coupled system (<ref>) and (<ref>).
While he followed a method-of-lines approach discretizing the spatial derivatives by fourth-order approximation to solve the resulting system of ordinary differential equations by generic integrators we try to make use of the structure revealed in appendix <ref>.
We revisited all the equations presented by <cit.> needed for closure and further performed the split of (<ref>) into linear and non-linear contributions which lead to a quasi-Hamiltonian substructure giving rise to the numerical scheme to be presented subsection <ref>.
In addition to the centerline equation (<ref>) time evolution of the tangential velocity is re-written as


    ∂_t u_θ + u_r,001/r∂_r(ru_θ) + w_0 ∂_z u_θ   = - u_r,*( u_θ/r + f ) - u_r,00f


to identify an advection term (in polar coordinates) on the left-hand side and source terms on the right-hand side.
The integration scheme of this equation is presented in subsection <ref>

By our analysis we learned that the system of equations (<ref>) and (<ref>) is assembled by prototypes of partial differential equations that are i) advection equation, ii) Schrödinger equation, and iii) non-linearly coupled ordinary differential equations (source terms).

Aiming at the scope of this work we restrict to asymmetric diabatic heat release which allows to drop all terms referring to symmetric vertical and radial motions.
We further neglect background wind shear.
Hence we drop all the advective terms from eqs. (<ref>) and (<ref>) and set u_s=0.



 §.§ Integration of centerline



For the integration of the centerline position the general solution strategy is to integrate nonlinear source terms by the trapezoidal rule and the Sturm-Liouville operator (after appropriate spatial discretization) by the implicit midpoint rule.
The choice of the latter is based on the idea of preserving unitarity during integration of the linear part of the equation.
We also dropped contributions according to shear as they are not discussed in section <ref>.
The composition of the sub-steps reads:


    X^*         = 1/2Δ t Q(u^n, X^n, t^n) + X^n 
    
      X^**   = ( 1 - 1/2i Δ t Ĥ(u^n+1/2)) X^* 
    
      X^***   = ( 1 + 1/2i Δ t Ĥ(u^n+1/2))^-1 X^**
    
      X^n+1   = 1/2Δ t + Q(u^n+1, X^n+1, t^n+1) + X^***


1 represents the identity operator.
For the evaluation of Ĥ we need u^n+1/2 which we obtain by the first-order predictor u^n+1/2 = u^*.



 §.§ Integration of tangential velocity



The integration of the tangential velocity equation is accomplished by applying the trapezoidal rule to the source term.
Note that after dropping symmetric diabatic contributions all differential expression in equation (<ref>) vanished and we only have to integrate the non-differential source term proportional to u_r,*.
The integration scheme for one timestep Δ t reads


    u^*        = -1/2Δ t u_r,*^n ( u^n/r + f )  + u^n  ,
    
      u^n+1   = -1/2Δ t u_r,*^n+1( u^n+1/r + f )  + u^* ,


For clarity we dropped the _θ subindex.
It becomes obvious that the final step of the integration involves an implicit solution strategy as the term u^n+1_r,* depends on both X and u at time level n+1.



 §.§ Coupled Integration


The above stated integration scheme involves information from previous sub-steps for the coupled integration.
In both, explicit and implicit forcing, the equations are couples to each other.
For the implicit (finalizing) step it is necessary to iterate solving X^n+1 and u^n+1 with second order:


    X^n+1,0   = X^***
    
      u^n+1,0   = u^**
    
      X^n+1, ν   = 1/2Δ t  + Q(u^n+1,ν-1, X^n+1,ν-1, t^n+1) + X^****
    
      u^n+1,ν   = -1/2Δ t u_r,*^n+1, ν( u^n+1,ν/r + f )  + u^*


This integration strategy is adopted from literature on implicit methods for fluid dynamics <cit.>.



 §.§ Details on the spatial discretization


The equations are discretized on an equidistant grid allowing for straightforward finite-difference approximations of the derivate operators.
Boundary conditions are accommodated by a extending the grid covering the physical domain plus a ghost layer of two cells.
Solution values are stored in cell-centers while first derivate are computed typically on the corresponding faces.

Prototypical differential expressions such as α∂_z(β∂_zψ) are discretized as:


    .α∂_z(β∂_zψ)|_z=z_i = 1/Δ z^2α_i (  β_i+1/2( ψ_i+1 - ψ_i)- β_i-1/2( ψ_i - ψ_i-1) )


Integrals are computed via the trapezoidal rule where ghost cells are obtained by quadratic extrapolation.
We further include an option to apply hyper-viscosity to the centerline stabilizing time integration with activated diabatic heating.
Further details may be taken from the source code available on demand by the corresponding author.








§ DETAILS ON THE NUMERICAL IMPLEMENTATION





 §.§ Dimensional variables


Though the derivation outlined above is carried out in terms of non-dimensional variables for the actual implementation into EULAG we used dimensional quantities of which some details will be presented in this section.
In the spirit of asymptotic analysis for reconstructing dimensional variables and formulate leading-order relations by using leading-order or next-to-leading order modes.

Before presenting specific relation which arise from the asymptotic analysis, we want to relate the expansion modes with mean background values, denoted by bars (·) and perturbations, denoted by primes (·).


 ρ   = ρ( ρ_0 + δ^2 ρ_2 + δ^4 ρ̂_4 + δ^5 )
         ρ   = ρ ( δ^4 ρ̂4 + δ^5 ) 

 
 p       = p ( p_0 + δ^2 p_2 + δ^4 p̂_4 + δ^5 )
         p    = p ( δ^4 p̂4 + δ^5 ) 

 
  Θ   = T ( Θ_0 + δ^2 Θ_2 + δ^4 Θ̂_4 + δ^5 )
         Θ   = T ( δ^4 Θ̂4 + δ^5 )


Furthermore, we have


    u_θ   =  u · e_θ = u( 1/δ u_θ0 + 1) 
    
     w    =  u( δŵ1 + δ^2)


A trivial, but useful observation is (·) = 0.



  §.§.§ Pressure


Equation (<ref>) balances pressure gradient with radial forces.
We find the dimensional version as


    1/ρpr = u_θ^2/r .




  §.§.§ Potential temperature


For the deviation of potential temperature from its background mean value Θ we have for the first Fourier modes


    Θ_1 = -Θ/g1/ρpr Xz .





  §.§.§ Vertical velocity


In general, , for arbitrary heating, the vertical velocity takes the following form in physical dimensions:


    w = 1/Θz( Q_Θ + Θ/gu_θ^3/r^2( _-π/2 Xz) · e_r )





  §.§.§ Diabatic heating


Finally, we defined a heating dipole aligned with an angle θ_0 relative to the tilt direction:


    Q_Θ^θ_0 =
              Θ/gu_θ^3/r^2( _θ_0 Xz) · e_r





 §.§ Convergence results



In addition to comparing both simulation results of three-dimensional and asymptotic equations we check for convergence of each numerical schemes.
For both simulations we check for self-consistency by comparing results at increasing resultion with high-resolved reference.
Convergence of the EULAG simulations is displayed in figure <ref> showing second-order convergence.




The numerical scheme solving the asymptotic equation on the other hand is tested by evolving


    X(t=0) = (cos(zπ/z_max), 0)^T


for T=0.1 (in asymptotic units) and comparing solutions of different resolutions against a reference solution with 1280 grid points.
The results are plotted in figure <ref> indicating second-order convergence as expected.












spbasic   
   

