## Abstract

Quantification of inertia–gravity waves (IGWs) generated by upper-level jet–surface front systems and their parameterization in global models of the atmosphere relies on suitable methods to estimate the strength of IGWs. A harmonic divergence analysis (HDA) that has been previously employed for quantification of IGWs combines wave properties from linear dynamics with a sophisticated statistical analysis to provide such estimates. A question of fundamental importance that arises is how the measures of IGW activity provided by the HDA are related to the measures coming from the wave–vortex decomposition (WVD) methods. The question is addressed by employing the nonlinear balance relations of the first-order *δ*–*γ*, the Bolin–Charney, and the first- to third-order Rossby number expansion to carry out WVD. The global kinetic energy of IGWs given by the HDA and WVD are compared in numerical simulations of moist baroclinic waves by the Weather Research and Forecasting (WRF) Model in a channel on the *f* plane. The estimates of the HDA are found to be 2–3 times smaller than those of the optimal WVD. This is in part due to the absence of a well-defined scale separation between the waves and vortical flows, the IGW estimates by the HDA capturing only the dominant wave packets and with limited scales. It is also shown that the difference between the HDA and WVD estimates is related to the width of the IGW spectrum.

## 1. Introduction

More than 20 yr since the pioneering work by O’Sullivan and Dunkerton (1995), the problem of inertia–gravity wave (IGW) generation by vortical flows has remained an active area of research [see Plougonven and Zhang (2014) for a recent review]. It is related to the fundamental problems of balance (Vanneste 2013; Bühler et al. 2014; Callies et al. 2014; McIntyre 2015) and wave–vortex decomposition (WVD; Bühler et al. 2014; Callies et al. 2014), as well as the more practical problems like parameterization of IGWs (Mirzaei et al. 2014; de la Cámara and Lott 2015) in global models of the atmosphere.

Two streams of research have approached the problem from two different views. In the first stream, the focus has been on the use of balance to carry out WVD in low-order models like the Lorenz–Krishnamurthy equations (Vanneste 2013), in single- and two-layer shallow-water flows (Mirzaei et al. 2012), and in dipolar flows of stratified fluids (Viúdez 2007; Snyder et al. 2007). In the second stream, the focus has been on the main characteristics of the IGWs including the frequency and horizontal and vertical wavelengths (Plougonven and Snyder 2005, 2007; Zülicke and Peters 2006). However, no attempt has been made to connect the two streams by comparing their measures of IGW activity. The importance of the problem has been recognized by Plougonven and Zhang (2007), who used scaling arguments to provide a wave equation describing the forcing of IGWs by the large-scale flow. To quantify the magnitude of IGWs generated by vortical flows, Zülicke and Peters (2006) have presented a harmonic divergence analysis (HDA) and applied it to the horizontal divergence fields of numerical simulations by the MM5. The HDA assumes a zeroth-order balance with vanishing horizontal divergence (*δ*) and thus takes the *δ* as representing the unbalanced field. Further, the HDA captures the wave structure and quantifies the strength of IGWs using the linear wave polarization relations. Previously, the HDA has been applied by Zülicke and Peters (2008) and Mirzaei et al. (2014) for the validation of a bulk parameterization of IGWs generated by jets, fronts, and convection. As its name suggests, the HDA’s working rests on certain assumptions on the wave field like the presence of a locally dominant wavenumber and sufficient separation with the large-scale balanced flow. By construction, the HDA performs well in regions of space filled by the coherent wave packet–like structures [see the analysis in Zülicke and Peters (2006)]. As generally, wave packets of different characteristics may exist simultaneously in various parts of the flow, the HDA is usually applied by partitioning the flow domain to several subdomains based on a judicious choice that can ideally take into account our knowledge of the variations in the background vortical flow and the characteristics of the wave packets present. In practice, however, this may not be possible because of large variations in spatiotemporal complexities of both waves and vortical flows. Therefore, a theoretical underpinning is required for applications of the HDA to estimate the energies of IGWs. The WVD methods based on inversion of a master variable representing vortical flow using various kinds–orders of balance (Warn et al. 1995; Vallis 1996; McIntyre and Norton 2000; Ford et al. 2000; Mohebalhojeh and Dritschel 2000; Bühler et al. 2014; Herbert et al. 2016) can provide the underpinning required. The connection to the notion of balance distinguishes these methods from the spatial filtering procedures used previously in Wang and Zhang (2007), Wei and Zhang (2014), and Wei et al. (2016) in the context of moist baroclinic simulations. A profoundly important aspect of balance to note here is its inherent fuzziness in the sense that, no matter how balance is defined, there can be violations of balance for the flows we are concerned with (Vanneste 2013). Associated with this fuzziness, there is no exact balance and no exact wave–vortex decomposition. Given the constraints set by this fundamental limitation, the waves and vortical flows can only be decomposed in an approximate sense, which can be sufficient for practical purposes.

The current work aims to compare the measures of IGW activity coming from the HDA with those of the WVD methods in the idealized numerical simulations of the dry and moist baroclinic instability by the Weather Research and Forecasting (WRF) Model. This idealized problem has attracted considerable attention in recent years (Waite and Snyder 2013; Wei and Zhang 2014; Wei et al. 2016). To this end, the space of balance relations is explored for the best WVD possible by examining and comparing three different sets of balance relations: first-order *δ*–*γ*, Bolin–Charney, and first- to third-order Rossby number expansion. The global measures of kinetic energy are then extensively compared for the estimates of IGWs by the HDA and WVD methods as determined by the instantaneous inversion of the vertical vorticity field.

The layout of the paper is as follows. Based on the hydrostatic Boussinesq primitive equations in their frictionless form, the balance relations are introduced in section 2 together with the resulting diagnostic equations for the balanced and unbalanced fields. Details of the WVD methods are given in section 3 followed by the presentation of the HDA in section 4. The numerical simulations carried out by the WRF are presented in section 5, with a detailed comparison of IGW activity with unbalanced fields in section 6. The sensitivities to model horizontal resolution and moisture are also examined. Finally, concluding remarks are given in section 7.

## 2. Balance relations

To fix ideas, a balance relation can be thought of as a diagnostic relation between the velocity and mass fields [see McIntyre (2015) for basic concepts and ideas]. The well-known geostrophic balance is the leading-order balance that is obtained in a conventional Rossby number expansion of the primitive equations. Knowing the limitations of the geostrophic balance, the search for exact balance was a main theme of research in the 1980s for both applications in initialization of NWP models and theoretical understanding of atmospheric flows. For exact balance to exist, one needs a balance relation that is exactly held during the time evolution for arbitrary flows of interest. Notwithstanding the nonexistence of exact balance (McIntyre 2015; Vanneste 2013), there may be balance relations that provide highly accurate approximations to a wide range of atmospheric and oceanic flows. The balance represented by balance relations that are minimally violated is called optimal balance (Mohebalhojeh and Dritschel 2000). Our aim here is to look for optimal balance and use it for wave–vortex decomposition. Optimality is judged here based on the extent to which the actual divergent flow is captured by the balanced flow in a global sense to be defined.

### a. Basic formulation

The Boussinesq horizontal momentum equation on an *f* plane is written as

with being the material derivative, **V** = (*u*, *υ*) is the horizontal velocity vector, ∇ is the horizontal gradient operator, *w* is the vertical velocity, *f* is the Coriolis parameter taken to be constant, is the unit vector in the local vertical direction, *ρ*_{0}(*z*) is the reference density, and *p* is the pressure. Taking (**⋅** ∇ ×) and horizontal divergence **(**∇ **⋅)** of (2.1), respectively, the vorticity *ζ* and divergence *δ* equations are obtained as

where the quantity

is the horizontal acceleration divergence, *γ* ≡ ∇ **⋅ ***D***V**/*Dt*, and *J* is the Jacobian operator. It is important to note that ∇ **⋅** and ∇^{2} are both horizontal. Now, taking the partial derivative of (2.4) with respect to time,

substituting (2.2) into (2.5), taking its partial derivative with respect to *z*, dropping the small term involving vertical derivative of *ρ*_{0} based on scaling arguments, we arrive at the following equation:

Using the hydrostatic balance relation, ∂*p*/∂*z* = −*ρg* with *ρ* the density and *g* the acceleration due to gravity and the density tendency as ∂*ρ*/∂*t* = −**V ⋅** ∇*ρ* − *w*∂*ρ*/∂*z* + *Dρ*/*Dt*, one can write the relation

for the last term in the right-hand side of (2.6) in which *C* = −*gDρ*/*Dt* is the term due to compressibility. For the mass continuity equation, the anelastic approximation [see (2.34) in Holton and Hakim (2013)]

is used. Taking the partial derivative of (2.8) with respect to *z*, the following equation is then obtained:

### b. The first-order δ–γ balance

The first-order *δ*–*γ* balance relation is obtained by setting the first time derivative of velocity divergence and acceleration divergence to zero (Mohebalhojeh and Dritschel 2001). Substituting (2.7) and (2.9) into (2.6), applying ∂*γ*/∂*t* = 0 and rearranging its terms, the omega equation is given by

Now, we set , with and *p*′ being the reference and perturbation pressure fields, respectively, and rearrange the linear terms in *w* in the left-hand side of (2.10) to obtain

where with , , , , and . Together with a balance relation for the pressure field, (2.11) will be used for wave–vortex decomposition in section 3a. To estimate the compressibility term *C*, the thermodynamic energy equation

is invoked in which *T* is the temperature, *c*_{p} is the heat capacity per constant pressure, and *q* the heating rate per unit mass. Using the equation of state, hydrostatic balance, and (2.12), the equation

can be obtained for pressure tendency in which *Q* = *q*/*c*_{p}, *κ* = *R*/*c*_{p}, *θ* is the potential temperature, and *R* denotes the specific gas constant. To obtain an estimate of *Q* suitable for slantwise convection, we have used the form for latent heating given by Emanuel et al. (1987) written here in *z* coordinate as

in which *θ*_{e} is the equivalent potential temperature and Γ_{d} and Γ_{m} denote, respectively, the dry and moist adiabatic lapse rates. Solving the first-order differential (2.13) together with (2.7) can give us an estimate of *C*. It should be noted here that for vertical scales smaller than a scale height based on the vertical average of *RT*/*g*, as an approximation, one can drop the second term in the left-hand side of (2.13) and arrive at the following explicit expression for *C*:

When exerted in (2.11), (2.15) will give us the familiar form of the Laplacian of diabatic heating in the omega equation. As a further approximation whose validity can be assessed a posteriori in the sense defined later, the compressibility term *C* is set to zero in the main body of results presented. For a selected case, however, results are also presented using nonzero values of *C* computed based on the thermodynamic information and a perturbative solution of (2.13).

### c. Bolin–Charney balance

The Bolin–Charney balance is constructed by first considering the balance relation

where *u*_{ψ} and *υ*_{ψ} are rotational components of the velocity field when it is decomposed into rotational and divergent parts (**V** = **V**_{ψ} + **V**_{χ}) as

in which *ψ* and *χ* are the streamfunction and velocity potential, respectively, which satisfy

The velocity field is replaced by its rotational part in the third term in the right-hand side of (2.2) and thus (2.6). Now, we can write the time derivative of the acceleration divergence *γ* as

### d. The Rossby number expansion

A well-known method to construct balance relations for the primitive equations is to employ an asymptotic expansion using a small parameter, the most natural candidates for which are the Rossby and Froude numbers. To resolve the issues with secular growth encountered in a conventional asymptotic expansion, Warn et al. (1995) developed the modified expansion where a master variable of balanced dynamics is left unexpanded [for further developments, see Vallis (1996) and Mohebalhojeh (2002)]. Without embarking on a formal Rossby number expansion procedure, the analysis in Mohebalhojeh (2002) for the shallow-water equations is extended to the continuously stratified primitive equations that we are dealing with here. To this end, first the vorticity and divergence equations are rewritten as

where, in (2.21) and (2.22), the leading-order terms have been moved to the left-hand side and the rest of the terms of order Rossby number and higher have been grouped together in the right-hand side. The {Ro} is simply used to denote the order in Rossby number. Next, we take the time derivative of (2.22) to arrive at the counterpart of (2.19) for the Bolin–Charney balance equation,

Now, a procedure similar to that leading to (2.11) but using the exact relation (2.23) instead of the balance relation ∂*γ*/∂*t* = 0 gives us the following form of the omega equation:

which is ready for Ro expansion. Note that, technically, we have not yet defined the balance relations of various orders in Rossby number. To make the exposition compact, the balance relations are introduced once details of their application are given in section 3c.

### e. Vertical mode space

The operator acting on the vertical velocity in the left-hand sides of (2.11), (2.20), and (2.24) is three-dimensional, making the inversion of the corresponding equations complicated. To reduce the dimensionality of the problem, the omega equations can be projected onto the vertical mode space defined by the eigenvectors of the matrix

which is obtained by a standard centered vertical discretization of the operator at *l* = 1, 2,…, *L* uniformly spaced layers numbered from the bottom (*z* = 0) to the top (*z* = *z*_{top}), where zero boundary conditions for *w* are used at *z* = 0 and *z*_{top}, with *L* the number of layers in the vertical direction and Δ*z* the vertical grid spacing. In this way, each of the three omega equations obtained can be compactly written in terms of *L* modified Helmholtz equations by projection onto the vertical modes as (Mohebalhojeh and Dritschel 2004)

## 3. The WVD methods

The three sets of balance relations briefly introduced in section 2 are applied to partition the horizontal velocity field into a balanced part **V**_{b} = (*u*_{b}, *υ*_{b}) and an unbalanced part **V**_{unb} = (*u*_{unb}, *υ*_{unb}) = (*u*, *υ*) − (*u*_{b}, *υ*_{b}). The same partitioning can be done for the other fields like *w* and *δ*. The balanced and unbalanced parts thus defined are assumed to be associated with, respectively, the vortical structures and freely propagating IGWs largely present in geophysical flows (Mohebalhojeh and McIntyre 2007). To carry out the partitioning, one can take the instantaneous field of potential vorticity as the master variable (Warn et al. 1995) and invert it using the balance relations. While it is possible to extend the procedure below to potential vorticity inversion, a much simpler procedure of vorticity inversion is followed below. This simplification is justified for horizontal scales smaller than the Rossby length or Rossby deformation radius characterizing the baroclinic jet. It can be violated for imbalances of large horizontal scale, that is, near inertial frequency IGWs, which are excluded from our analysis on the basis that the IGWs generated in the idealized simulations are dominantly mesoscale features.

### a. The δ–γ

To complete the description in section 2b for the omega equation, one needs to substitute (2.4) in (2.3), set ∂*δ*/∂*t* = 0, and rearrange the result with respect to *p*. In this way, the balance relation for the pressure field is obtained:

The procedure to obtain the balanced part of the velocity field is carried out in an iterative manner starting from the actual vorticity field given by *ζ*^{[0]} = *ζ* = ∂*υ*/∂*x* − ∂*u*/∂*y* with *u* and *υ* the components of the actual velocity field and the condition of zero divergence (i.e., *δ*^{[0]} = 0). The components of the initial, or the first guess, velocity field (*u*^{[0]}, *υ*^{[0]}) are computed by inverting (2.18) while making use of (2.17). Note that the superscript [0] is used for the first guess fields of the iterative procedure. Then (3.1) is inverted to obtain the initial pressure field (*p*^{[0]}).

In the first iteration, the forcing term *F*_{1} of the omega equation [(2.11)] and subsequently the vertical velocity (*w*^{[1]}) is computed by solving the modified Helmholtz equation [(2.26)]. Now, the new divergence field (*δ *^{[1]}) can be determined from (2.8). Then the velocity field is updated to (*u*^{[1]}, *υ*^{[1]}) based on the new *δ *^{[1]}, noting that *ζ* is fixed in this vorticity inversion procedure. For the second iteration, the current solution for the horizontal velocity field (*u*^{[1]}, *υ*^{[1]}) is used to compute *p*^{[1]} from (3.1) and the rest of the steps outlined above are repeated. The horizontal velocity field thus determined, (*u*^{[2]}, *υ*^{[2]}), constitutes the current approximation to the balanced velocity field. In principle, the foregoing iterative process can be carried out arbitrarily to seek for optimal balance (Mohebalhojeh and Dritschel 2001) in the sense defined early in section 2. In this way, the iterative process is treated as an iterative map and each iteration is considered as defining a balance of its own (Leith 1980).

### b. Bolin–Charney

By combining (2.4) and (2.16), the Bolin–Charney balance replaces (3.1) of the first-order *δ*–*γ* balance with

The iterative procedure designed for the Bolin–Charney is the same as that described for the *δ*–*γ*, the differences being that in the vorticity inversion, the pressure field is computed noniteratively from (3.2), and that in the omega equation [(2.26)], *F*_{2} is used instead of *F*_{1}. The latter change involves the inclusion of the first time derivatives of the rotational velocity field in the iterations of the Bolin–Charney omega equation. The required fields (), the superscript [*n*] being *n*th iteration, can be estimated by treating the vorticity (2.2) and the divergence (2.3) equations in an iterative manner starting from ∂*ζ *^{[0]}/∂*t* = −∇⋅ (**V**^{[0]}*ζ *^{[0]}) and ∂*δ*^{[0]}/∂*t* = 0.

### c. Rossby number expansion

The third method to obtain the balanced part of flow is carried out based on the Rossby number expansion of the vertical velocity (*w*) up to as

where *w*^{(0)}, *w*^{(1)}, *w*^{(2)}, and *w*^{(3)} indicate, respectively, the zeroth-, first-, second-, and third-order contributions to the vertical velocity. From (2.24), *w*^{(1)}, *w*^{(2)}, and *w*^{(3)} are computed from the following relations:

where

with *T*^{(0)} = *p*^{(0)}/[*RT*^{(0)}], 1/ *T*^{(1)} = [*ρ*^{(1)}/*ρ*^{(0)} − *p*^{(1)}/*p*^{(0)}]/*T*^{(0)}, and *Y* denoting the right-hand side of (2.13). For brevity, the relations for *C*^{(2)}, *Y*^{(2)}, and ∂*p*^{(2)}/∂*t* have not been presented. Noting that *w*^{(0)} = 0, *w*^{(1)} can be obtained from (3.4). To this end, first *p*^{(0)} is computed from ∇^{2}*p*^{(0)} = *fρ*_{0}*ζ*^{(0)}, that is, the zeroth-order balance *γ*^{(0)} = 0, in (2.4) and using *ζ*^{(0)} = *ζ* as determined by the actual velocity field. Then the zeroth-order velocity field [*u*^{(0)}, *υ *^{(0)}] is obtained by inverting *ζ*^{(0)} and *δ*^{(0)} = 0 using (2.17) and (2.18). Now, *w*^{(1)} is readily determined by inverting (3.4), from which the first-order divergence field *δ*^{(1)} is determined from (2.8) and, consequently, [*u*^{(1)}, *υ*^{(1)}] field from [*δ*^{(1)}, *ζ*^{(1)} = 0]. It is then possible to determine [] field from ∂*ζ*^{(0)}/∂*t* = −∇⋅ [**V**^{(0)}*ζ*^{(0)}] and ∂*δ*^{(0)}/∂*t* = 0. Substituting (2.4) into (2.22) leads to ∇^{2}*p*^{(1)} = 2*ρ*_{0 }*J*[*u*^{(0)}, *υ*^{(0)}] which is applied to obtain the first-order pressure field [*p*^{(1)}]. This is used to determine *w*^{(2)} by solving (3.5) and thus *δ*^{(2)} from (2.8). Therefore, we would be able to compute [*u*^{(2)}, *υ*^{(2)}] field using [*δ*^{(2)}, *ζ*^{(2)} = 0].

To reach the third-order formal accuracy, we carry on the procedure above to obtain by taking the time derivative of (3.4), from which ∂*δ*^{(1)}/∂*t* and then [] are followed. Now *p*^{(2)} can be determined from the following relation:

The computation of the other fields like ∂^{2}*δ*^{(1)}/∂*t*^{2} required for the right-hand side of (3.6) follows a similar procedure as that outlined above (3.8). Then *w*^{(3)} would be obtained by inverting (3.6), from which *δ*^{(3)}, and then using *ζ*^{(3)} = 0, the third-order velocity field [*u*^{(3)}, *υ*^{(3)}] can be determined.

In the iterative procedures adopted for the *δ*–*γ* and Bolin–Charney, the velocity field thus obtained up to *n*th iteration is regarded as the balanced field and denoted by **V**_{b} ≡ . In the Rossby number expansion, the velocity field asymptotically accurate to *n*th order is regarded as the balanced field and denoted by given by

A similar definition is used for the balanced divergence *δ*_{b}. The computations required for each of the wave–vortex decomposition methods are carried out using the Fourier spectral transform and the second-order central differencing in, respectively, the zonal and meridional directions. Details of the solution procedure for the elliptic equations and the computation of the velocity field can be found for a channel model setting in Mohebalhojeh and Theiss (2011). What is helpful to note here is the use of zero boundary conditions for *w* in solving (2.26) and a spectral filter for dealiasing whenever a zonal derivative is computed.

## 4. The HDA

As introduced by Zülicke and Peters (2006) and used in Zülicke and Peters (2008) and Mirzaei et al. (2014), the HDA is applied to determine the characteristics of the IGWs. In application of the HDA, the domain is divided into a number of nonoverlapping sample boxes of *L*_{x,box}, *L*_{y,box}, and *L*_{z,box} lengths in, respectively, *x*, *y*, and *z* directions, covering the analysis domain from *z* = 0 to *z* = 22 km, separately for the troposphere (*z* = 0 to *z* = 11 km) and the stratosphere (*z* = 11 to *z* = 22 km). Unless stated otherwise, the HDA is applied with a medium box size of 2000 × 2000 × 11 km^{3}. As estimated in the appendix of Zülicke and Peters (2006), the maximum detectable wavelength is 0.4 times the box length, which implies typical subsynoptic wavelengths below 800 km in the horizontal direction.

This choice for the box size is based on a meticulous visual inspection of the inertia–gravity wave packets generated whose characteristics need to be captured by the HDA as accurately as possible. Also examined are the large and small box sizes of, respectively, 4000 × 5000 × 11 km^{3} and 1000 × 1000 × 11 km^{3}. The use of the large box size will lead to an approximation for a wavelength of the baroclinic wave and thus gross underestimation of the impacts of variations in the background flow on the wave field. The use of the small box size turns out to underrepresent larger scale and more energetic parts of the IGW spectrum, a point that we will return to in the analysis of the results. One should note that in application of the HDA, only the dominant wavelength within a box is selected. This is optimal for a localized harmonic wave packet in a slowly varying environment but might fail for broader wave spectra or large-scale unbalanced flow components. The dominant wavelengths *λ*_{x}, *λ*_{y}, and *λ*_{z} and variances , , and in the *x*, *y*, and *z* directions are computed by analyzing the empirical autocovariance function of the divergence field in each box. In this way, the intrinsic frequency *ω* can be obtained from the dispersion relation for hydrostatic IGWs (Mirzaei et al. 2014):

in which *f* ≈ 10^{−4} s^{−1} is the Coriolis parameter; *N* ≈ 10^{−2} s^{−1} is the Brunt–Väisälä frequency; *k*_{h} = ()^{1/2} is the total horizontal wavenumber with *k*_{x} and *k*_{y}, respectively, the zonal and meridional wavenumbers; and *k*_{z} is the vertical wavenumber. For a plane wave solution, a straightforward manipulation in spectral space using the relations *ζ* = ∂*υ*/∂*x* − ∂*u*/∂*y* and *δ* = ∂*u*/∂*x* + ∂*υ*/∂*y* gives = ()/, in which denotes amplitude in the Fourier space for the given quantity *X*. This relation holds irrespective of the amplitude of the wave. To estimate energy from the divergence field alone, the linearized form of (2.2), ∂*ζ*/∂*t* = −*fδ*, is then invoked to give = (*i**f* /*ω*), which leads us to the following estimate for IGW kinetic energy of each box:

In HDA, the total divergence variance *s*^{2} = (1/3)() is used to estimate of the IGWs.

## 5. WRF simulations

The numerical setup used is the dynamical core of the Advanced Research version of the WRF (ARW) Model (Skamarock et al. 2008), which integrates the fully compressible, nonhydrostatic form of the primitive equations. The simulations were performed in a channel of *L*_{x} = 4000-km length, *L*_{y} = 10 000-km width, and *L*_{z} = 30-km height on an *f* plane. The fields of interest are much narrower, but the large meridional width prevents unphysical boundary effects. The boundary conditions are periodic in the *x* direction and symmetric (free-slip wall) in the *y* direction. It should be mentioned that preliminary experiments with the open lateral boundaries in the *y* direction led to a significant nonconservation of energy during the lifetime of the baroclinic wave. Given this observation and the fact that the lateral boundaries are sufficiently far away from the main area of baroclinic wave (BCW) activity, the symmetric boundary conditions with their better energy conservation provide a sensible choice for our purpose as in Mirzaei et al. (2014). The bottom boundary condition is specified as a free-slip condition. An absorbing layer with *w*-Rayleigh damping (Klemp et al. 2008) including a damping coefficient of 0.2 s^{−1} is used in the upper 8 km of the model to prevent the reflection of vertically propagating gravity waves.

To construct initial conditions, an ideal two-dimensional uniform distribution of potential vorticity with values of 0.4 and 4 PVU (1 PVU = 10^{−6} K kg^{−1} m^{2} s^{−1}) in the troposphere and stratosphere, respectively, is first inverted, which provides a two-dimensional baroclinic jet for the cyclonic BCW or the cyclonic life cycle LC2 in the classification of Thorncroft et al. (1993). The most unstable normal mode is then superposed on the two-dimensional jet [see Rotunno et al. (1994) and Plougonven and Snyder (2007) for details]. The initial conditions employed here differ from those in Mirzaei et al. (2014) where for the troposphere and stratosphere, respectively, the uniform distributions of 0.7 and 4.8 PVU were used. For both the jet constructed here and that in Mirzaei et al. (2014), the low-level temperature is much higher than what we expect in normal midlatitude conditions. For a more realistic temperature field, the potential temperature is uniformly reduced and then the other fields, including velocity, pressure, and density, are adjusted using the geostrophic and hydrostatic balance relations. Figure 1a shows the potential temperature, zonal velocity, and position of the tropopause in the initial jet.

The model was run with a horizontal resolution of 25 km and a vertical resolution of 250 m in the reference case. Additionally, the sensitivity to resolution was explored by running the model with horizontal resolution of 50 (12.5) km and vertical resolution of 500 (125) m to give a lower (higher)-resolution simulation called LOW (HIGH) (see Table 1). The model also uses a tri-harmonic horizontal hyperdiffusion, *ν*_{h}∇^{6}, described by Knievel et al. (2007) to filter nonphysical structures at the smallest scales of the flow. This explicit hyperdiffusion acts on the wind components, potential temperature, moisture fields, and subgrid turbulence kinetic energy. For the reference run, a hyperdiffusion coefficient of *ν*_{h} = 2.5 × 10^{21} m^{6} s^{−1} was used as in Mirzaei et al. (2014).

To initialize the model, a slightly modified form of the ideal profile of relative humidity due to Tan et al. (2004) was added to the jet where the relative humidity varies in both *y* and *z* directions:

In (5.1), RH is relative humidity, *z*_{RH} = 9000 m, and *z*_{0} = 1000 m; RH_{0} is set to 0.4, 0.55, 0.7, 0.85, and 1.0 indicating M40, M55, M70, M85, and M100 moist cases, respectively, and to 0.85 for the HIGH and LOW resolution cases; and *R*(*y*) is given by

in which *R*_{1} = 0.5, *R*_{2} = 1.0, *y*_{1} = 0.6*y*_{c}, *y*_{2} = 1.0*y*_{c}, and *y*_{c} = 5050 km. The above five moist experiments (M40, …, M100) as well as LOW and HIGH cases (see Table 1) were performed using the Kain–Fritsch cumulus and the Kessler microphysics parameterization schemes. For the dry case (DRY) where water vapor was retained as a tracer in the model (Mirzaei et al. 2014), the cumulus and microphysics parameterization were turned off. The vertical profile of relative humidity after the initialization carried out by the WRF for the M85 has been shown in Fig. 1b. The high values of relative humidity combined with the cold part of the flow result in high values of the convective available potential energy (CAPE) in the south of the domain (Fig. 1c). Comparison of initial conditions in different moist cases shows that the highest values of the CAPE and convective inhibition (CIN) belong to the M100 and M85 cases, respectively (as in Tan et al. 2004), while the initial CAPE and CIN are zero in the rest of the moist cases (Figs. 1c,d).

## 6. Results

### a. Balanced and unbalanced flow pattern

The baroclinic instability of the jet with core speeds in excess of 50 m s^{−1} leads to flows with a Lagrangian Rossby number (Ro) of order unity. This can be seen by considering either Ro = −|(*D***V**/*Dt*)|/*f*|**V**| as in Zülicke and Peters (2006) or Ro = 1/*τf* as in McIntyre (2009), *f* = 10^{−4} s^{−1}, and a jet exit region of about 250-km length and 20 m s^{−1} horizontal velocities appearing during the time evolution (Mirzaei et al. 2014). To give an impression of the structure of BCW and the IGWs generated, related cross sections are shown in Fig. 2 for the growth stage of BCW at day 4 of the M85 case. The growth rate in this case is faster than the MOIST case in Mirzaei et al. (2014) and the EXP100 in Wei and Zhang (2014) such that the maximum amplitude of the BCW occurs 36 and 20 h earlier, respectively. Four wave packets named WP1, WP2, WP3, and WP4 appear at day 4 (Fig. 2a). Generally, for each of the four wave packets, we can identify a counterpart in the previous study by Mirzaei et al. (2014), the only difference of note being that the WP4 is weaker here. In the same way, it can be shown that the WP1 (Fig. 2c), WP2, and WP3 (Fig. 2b) are counterparts of WP4, WP5, and WP1 in Wei and Zhang (2014) and Wei et al. (2016), respectively. For the cross sections shown in Figs. 2b and 2c, the balanced and unbalanced vertical velocity fields, *w*_{b} and *w*_{unb}, respectively, as determined by the second-order Rossby number expansion are shown to demonstrate their close relation with the upper-level jet–front system as well as the surface front (Fig. 3).

To further compare the unbalanced and balanced fields, Fig. 4 shows the horizontal distribution of *δ*_{b}, *δ*_{unb}, *w*_{b}, and *w*_{unb} obtained by the second-order Rossby number expansion related to the baroclinic wave at *z* = 11 km near the tropopause level. The association of the balanced fields with the jet–front system is evident in the quadrupole structure of *w*_{b} resembling that in Fig. 11 of Snyder et al. (2007). The IGW packets exhibit patterns similar to those in Mirzaei et al. (2014). Quantitatively, the most important property to note is that, on average, the extrema of balanced divergence are an order of magnitude smaller than those of the unbalanced divergence. Based on this, one may argue that, as is the case for the zeroth-order Rossby number expansion, the balanced divergence is redundant for the analysis of the wave fields. Since the whole point of the current analysis is the relation between the energy estimates for IGWs captured by the HDA and that of WVD, ever higher accurate estimates for balance are however required. This notion of higher accuracy should be understood within the inherent limitations of the notion of balance arising from the nonexistence of a slow manifold in an exact sense.

### b. BCW energy

By splitting the actual velocity field into zonal mean and the eddy components defined by deviation from the zonal mean, **V** = **V**_{Z} + **V**_{E}, one can write a corresponding partitioning of the kinetic energy to the zonal and eddy components. The time evolution of the resulting eddy component, called eddy kinetic energy and denoted by EKE, for the M85 case illustrates the nonlinear growth stage of BCW involving a triple-peak evolution of EKE (see Fig. S4a in the supplementary material). The initial growth stage of BCW leading to the primary peak at around day 6 is followed by the secondary peak at around day 9 and a smaller, tertiary peak at around day 12.

### c. Energies of vortical flows and IGWs

Given the balanced **V**_{b} and unbalanced **V**_{unb} parts of the velocity field obtained by applying the three WVD methods described in section 3, we are able to obtain estimates for the energies of the vortical parts of the flow and the IGWs. To this end, each of the balanced and unbalanced velocity fields is partitioned into zonal mean and eddy deviation parts like **V**_{b} = **V**_{bZ} + **V**_{bE} and **V**_{unb}= **V**_{unbZ} + **V**_{unbE}. The eddy kinetic energies of the vortical and wave parts of the flow averaged over the entire domain, denoted respectively by EKE_{V} and *K*_{WVD}, come from

where **V**_{bE} = (*u*_{bE}, *υ*_{bE}) and **V**_{unb} = **V**_{unbE} = (*u*_{unbE}, *υ*_{unbE}) because zonal mean flow is exactly balanced. We have used *K*_{WVD} specifically to denote that the estimate for IGW is due to the wave–vortex decomposition to distinguish it from the corresponding estimate by the HDA. As the estimates for each of EKE_{V} and *K*_{WVD} depend on the balance relations used, a name convention is adopted here. That is, we use *X*_{b}(dg^{[n]}), *X*_{b}(bc^{[n]}), and *X*_{b}[re^{(n)}] for the vortical estimates of any given variable *X* coming from, respectively, the *δ*–*γ* and the Bolin–Charney in their *n*th iterations and the *n*th-order Rossby number expansion. For the inclusion of compressibility by consideration of the thermodynamic energy equation in WVD as presented in section 2b, *X*_{b}[re^{(n,c)}] is used to denote the *n*th-order Rossby number expansion. The optimal WVD is taken to be the one with the smallest values of *K*_{WVD} during the main wave generation events at the peak of baroclinic activity.

#### 1) Energy of vortical flow

Given the similarity of behavior in EKE_{V} among the three methods, it suffices to give here only an account of the time evolution of the energy of vortical flow coming from the Rossby number expansion. For the M85 case, although EKE_{V}[re^{(2)}] follows an evolution nearly the same as EKE, [re^{(2)}] =1.22 × 10^{6} J m^{−2} shows a decrease of 1.6% relative to = 1.24 × 10^{6} J m^{−2} . Here, and in what follows, the hat on a variable is used to refer to its peak value during the time evolution of the BCW. The difference between the two energies is related to the energy of unbalanced flow described below.

#### 2) Energy of IGWs by WVD: The *δ* –*γ*

The energy of IGWs as estimated by the first iteration of the *δ*–*γ* in the WVD, that is *K*_{WVD}(dg^{[1]}), follows a similar evolution to EKE and EKE_{V}, but (dg^{[1]}) = 1.79 × 10^{4} J m^{−2} is less than by two orders of magnitude, as shown in Fig. 5 and Table 2. The ratio of (dg^{[1]}) to is estimated to be about 1/69. The results also demonstrate the higher value of (dg^{[1]}) relative to (bc^{[1]}) and [re^{(2)}].

The effect of using a higher number of iterations in the inversion procedure on the energy estimates of IGWs (section 3a) is explored next. While the third and fourth terms of *F*_{1} in the right-hand side of (2.11) are eliminated in the first iteration, making it equivalent to the omega equation of the quasi-geostrophic model, these terms are retained in the next iterations. Results show a marked increase in the energy estimates of IGWs in the second iteration (Fig. 6a), such that (dg^{[2]}) becomes 8 times greater than (dg^{[1]}). The third iteration leads to even more dramatic increase in (not shown here), which points to the divergent nature of the inversion process due to the generation of locally large values of vertical velocities and thus large values of the *F*_{1} term in the right-hand side of (2.26). Beyond the first iteration, the first-order *δ*–*γ* thus leads to improper estimates of IGW activity.

#### 3) Energy of IGWs by WVD: Bolin–Charney

The IGW energy estimates of WVD by the Bolin–Charney balance relation *K*_{WVD} (bc^{[1]}) (Fig. 5) indicate the same evolution as that of *K*_{WVD} (dg^{[1]}) but with lower values, especially at day 8 when *K*_{WVD} peaks. It should be mentioned that while both the Bolin–Charney and *δ*–*γ* start from the same first guess field, the right-hand side of (2.20) differs from that of (2.11) in the first iteration because of the presence of the Jacobian term. For the Bolin–Charney, in the second iteration, two additional terms are added in the computation of the right-hand side of (2.20). The effect of these two terms makes (bc^{[2]}) about 1.2 times greater than (bc^{[1]}), while the (dg^{[2]})/(dg^{[1]}) ≈ 8 (cf. Figs. 6a and 6b). That is, compared with the *δ*–*γ* balance, the results of the Bolin–Charney converge more rapidly with iteration.

#### 4) Energy of IGWs by WVD: Rossby number expansion

The time evolution of the energy estimates of WVD as carried out by the second-order Rossby number expansion follows a trend similar to that of the *δ*–*γ*. Taking a value of about 1.44 × 10^{4} J m^{−2}, [re^{(2)}] is less than both (dg^{[1]}) = 1.79 × 10^{4} J m^{−2} and (bc^{[1]}) = 1.48 × 10^{4} J m^{−2} (Fig. 5 and Table 2).

Comparing the first-, second-, and third-order Rossby number expansion results, one can see a clear change in behavior between different orders of the asymptotic balance relations (Fig. 6c). Whereas (re) undergoes a 4.6% reduction by going from the first to second order, it increases markedly, by 27.1%, when going from the second to the third order. This observation points to the asymptotic nature of the Rossby number expansion in such a way that the minimal value of *K*_{WVD} and thus optimal result comes from the second-order Rossby number. That is, no gain is obtained by going to the third order and higher. As for the *δ*–*γ* with one iteration, the first-order asymptotic expansion is equivalent to the quasi-geostrophic model. The ratio (dg^{[1]})/[re^{(1)}] is about 1.18. This apparent discrepancy between the latter two IGW energy estimates comes from our choice of the starting solution for the pressure field from (3.2) in the *δ*–*γ*, which makes it different from the zeroth-order *γ* = 0 solution used in the first-order vertical velocity field in the asymptotic expansion.

The inclusion of compressibility by using the thermodynamic energy equation in WVD as presented in section 2b (the red dashed–dotted line in Fig. 5 and dashed black line in Fig. 6) leads to an increase in *K*_{WVD} with respect to both the first- and second-order Rossby number expansion that set *C* to zero. This behavior can be understood by noting that in our vertical vorticity inversion WVD, no thermodynamic information is directly available from the actual flow. The approximation of the thermodynamic state in the WVD may result in significant error in the distribution of localized centers of diabatic heating. The adverse effects of such errors on the resulting balanced flow is thought to be responsible for the failure of compressibility in related energy estimates. The extent to which a potential-vorticity-based WVD may overcome this issue remains to be studied in the future. Considering *K*_{WVD} values obtained for the WVD methods examined, the second-order Rossby number expansion with *C* = 0 can be regarded as being optimal.

To further appreciate the impact of balance beyond quasi geostrophy, presented in Fig. 7 are the percentage relative differences in *K*_{WVD} and *l*_{2} norm of *δ*_{unb} as defined by, respectively, 〈{*K*_{WVD}[re^{(n)}] − *K*_{WVD}[re^{(n-1)}]}/ *K*_{WVD}[re^{(n)}]〉 × 100 and 〈*l*_{2}{*δ*_{unb}[re^{(n)}] − *δ*_{unb} [re^{(n-1)}]}/*l*_{2} {*δ*_{unb}[re^{(n)}]}〉 × 100 for *n* = 1, 2, 3. Here, the *l*_{2} norm of a general field *X* is defined by (Σ_{i.j,k}*X *^{2})^{1/2}, where summation is taken over the grid points in the analysis domain. During the main life span of the baroclinic wave over the first 15 days, generally the smallest relative differences are observed between the first- and second-order estimates. From day 15 onward, the flow relaxes to a state of balance of even higher order, and thus, the smallest relative differences become those between the second- and the third-order estimates. Judged based on the maxima of relative differences, the inclusion of the second-order balance has impacted *K*_{WVD}, when comparison is made with that of quasi geostrophy, by near to 20% and 40% for, respectively, *t* ∈ [0, 15] days and *t* ∈ [15, 50] days. The impact may seem small for the flows of order unity Rossby number involved. The smallness should be understood considering that the vorticity-based WVD regards the whole vorticity field as being balanced, which may lose validity at large scales. Greater impacts are expected if a potential-vorticity-based WVD is used. The same impacts are nearly 7% and 8% for *l*_{2}(*δ*_{unb}). In the absence of an exact solution, the 20% change in the quasi-geostrophic result can also serve as an estimate for uncertainty on *K*_{WVD} in the time interval *t* ∈ [0, 15] days. The corresponding uncertainty on the second-order result for *K*_{WVD} in *t* ∈ [15, 50] days is estimated at around 5% when comparison is made by the results of the third-order Rossby number expansion.

### d. Energy of IGWs by HDA

Here, the IGW kinetic energy averaged over the horizontal domain is estimated by

in which is the box-averaged density, is the kinetic energy in the box *i* as described in section 4, and the summation is over all the boxes.

We first explore the possibility of using unbalanced divergence instead of actual divergence in the HDA. This change in the HDA results in 1.6% reduction in but leaves the overall pattern unaffected. The time evolution of *K*_{HDA} for the M85 case estimated by actual divergence has been shown in Fig. 5. Although *K*_{HDA} involves three peaks as *K*_{WVD} of the *δ*–*γ*, Bolin–Charney, and Rossby number expansion, the major peak of the IGW energy estimate of the HDA = 0.64 × 10^{4} J m^{−2} occurs at day 6, that is, two days earlier than the time when the peaks of (dg^{[1]}), (bc^{[1]}), and [re^{(2)}] occur.

Moreover, the IGW energy estimates obtained by the HDA are clearly lower than those of the three WVD methods. The lowest value for WVD, [re^{(2)}], is nearly twice the value of (Fig. 5 and Table 2). On the other hand, the ratio of to for the M85 is assessed to be about 1/161, which is comparable to the corresponding ratio of 1/200 reported in Mirzaei et al. (2014) despite the differences in the setup. Two sensitivity tests have also been carried out. First, the results for increasing *N* to 2 × 10^{−2} s^{−1} in the stratospheric sample boxes of the HDA (see Fig. S1) showed little difference from that in Fig. 5. Second, the HDA was also applied using the large and small box sizes defined in section 4. As shown in Fig. 5, the IGW energy estimates by HDA exhibit a large degree of sensitivity to the choice of the box size. The effect is such that halving or doubling the box size in the horizontal direction, relative to the medium box size selected here, will lead to a substantial change in *K*_{HDA} values. This change results from the limitation of the dominant wave captured by the HDA to less energetic, shorter scales as the box size decreases. There is also a change in variance: the integral over the spectrum extends down to the smallest wavelength, which is inversely proportional to the box length. The sensitivity to box size further justifies our quest for a detailed comparison with the results of WVD methods.

### e. Impact of initial humidity on energy

A quantitative assessment of the impact of initial humidity on the growth rate of the BCW shows that the peak value of eddy kinetic energy increases by about 17.6% from the DRY to the M85, whereas the corresponding increase reported by Mirzaei et al. (2014) was 15%. Moreover, the (DRY) is 1.05 × 10^{6} J m^{−2}, nearly twice the 0.58 × 10^{6} J m^{−2} value obtained in Mirzaei et al. (2014). These differences between the current study and Mirzaei et al. (2014) could come from differences in initial conditions and initial humidity profiles. Note that in the current initial jet, overall, we have lower temperatures in the troposphere than the initial jet in Mirzaei et al. (2014). Despite the fact that the initial relative humidity varies with the same 15% interval between the moist cases, the relative enhancement in from M40 to M55, M55 to M70, M70 to M85, and M85 to M100 cases are about 1.3%, 1.3%, 4.4%, and 4.5%, respectively (Table 2). Therefore, the impact of moisture does not obey a simple linear relation. A similar behavior is observed in the growth rates reported in Wei and Zhang (2014) (see their Fig. 2a). When compared with the DRY case, the largest relative increase in seen for the M100 amounts to 23%.

Results show that the increase in initial humidity markedly affects the *K*_{WVD} values obtained by the three WVD methods. From the DRY to M100, while there is a 23% increase in (Table 2), (dg^{[1]}), (bc^{[1]}), and [re^{(2)}] increase by about 76.8%, 69.4%, and 67.4%, respectively. Our results also show that (dg^{[1]}) is greater than both (bc^{[1]}) and [re^{(2)}] in all cases. The smaller values of sensitivity to humidity for *K*_{WVD} [re^{(2)}] and *K*_{WVD} (bc^{[1]}) when compared to that for *K*_{WVD} (dg^{[1]}) is consistent with the measures of the global kinetic energy of IGWs by the WVD in Fig. 5.

Investigating the impact of increasing initial humidity on *K*_{HDA} as obtained by applying HDA to the actual divergence shows that is enhanced by about 129.4% from the DRY to M100. The time evolution of the IGW kinetic energy has also been computed by applying HDA to unbalanced divergence fields of the first iterations of the *δ*–*γ* and the Bolin–Charney, as well as the second-order Rossby number expansion, denoted respectively by *K*_{HDA} (dg^{[1]}), *K*_{HDA} (bc^{[1]}), and *K*_{HDA} [re^{(2)}]. Results show that increasing the initial humidity from the DRY to the M100 changes (dg^{[1]}), (bc^{[1]}), and [re^{(2)}] by about 192.1%, 148.4%, and 132.2%, respectively. The latter values are overall higher than the corresponding values for as well as [i.e., that obtained by applying HDA to the actual divergence (Table 2)].

The sensitivity to moisture content is summarized in Fig. 8, which presents the scatterplots of (DRY), (DRY) for the second-order Rossby number expansion, and (DRY). Here, the denotes time average taken from day 1 to day 14, during which the cases exhibit meaningful differences in the energy estimates. Being in overall agreement with the foregoing results based on the peak values, the scatterplots for the time averages point to transition to steeper slopes and thus higher moisture sensitivities around the M55 case for both and .

### f. Sensitivity of energy to horizontal resolution

The time evolution of EKE at different horizontal resolutions including the LOW (50), M85 (25), and HIGH (12.5 km) cases has been shown in Fig. 9a. While increases by about 1.6% from the LOW to the M85, it increases by about 6.4% from the M85 to the HIGH. Therefore, increasing resolution sharply intensifies , especially at the highest resolution examined where finescale processes are captured and affect BCW energy.

The sensitivity to horizontal resolution of *K*_{WVD} by the second-order Rossby number expansion has been shown in Fig. 9b. By going from LOW to HIGH, (dg^{[1]}), (bc^{[1]}), and [re^{(2)}] increase by about 76%, 60.5%, and 59.8%, respectively (Table 2). The corresponding increase in is only 8.2%. Therefore, the sensitivity of *K*_{WVD} estimates to resolution is greater than that of EKE, irrespective of the method used to carry out WVD. This means that IGW generation does not obey a linear relation with the strength of the vortical flow. It is also worth mentioning that *K*_{WVD}(dg^{[1]}) and *K*_{WVD} [re^{(2)}] exhibit, respectively, the highest and lowest sensitivity to resolution. As sensitivity to initial humidity, this is consistent with the behavior seen in Fig. 5 in the sense that closeness to optimality decreases sensitivity to both moisture content and resolution.

In the same way, the sensitivity to horizontal resolution of *K*_{HDA} as determined by applying HDA to the actual divergence is explored in Fig. 9c. The increase in from LOW to HIGH is estimated to be about 128.2%. In comparison, (dg^{[1]}), (bc^{[1]}), and [re^{(2)}] increase from LOW to HIGH by 195.2%, 139.5%, and 138.9%, respectively. The quantitative details are presented in Table 2. The IGW energy estimates of HDA as applied to each of the unbalanced divergence fields exhibit higher sensitivity to resolution than the estimate obtained when the actual divergence is used. At the same time, in each measure, the IGW energy estimates of HDA exhibit sensitivities to horizontal resolution that are greater than twice the corresponding sensitivities of *K*_{WVD}.

### g. Kinetic energy and divergence spectra

More insight can be gained on the nature of the relation between the measures of IGWs by the WVD and HDA methods by looking into the kinetic energy and divergence spectra. The two-dimensional divergent kinetic energy per unit mass in the spectral space spanned by the discrete Fourier basis functions is given by (Waite and Snyder 2013)

where *S*(*k*_{x}, *k*_{y}) = *δ*(*k*_{x}, *k*_{y})*δ**(*k*_{x}, *k*_{y}), *δ*(*k*_{x}, *k*_{y}) is the discrete Fourier transform of the horizontal divergence, and the superscript * denotes complex conjugate. A one-dimensional *k*-dependent spectrum can be obtained by summation over circular bands of constant *k* as in (Errico 1985)

for *k* = *l*Δ*k*, *l* = 0, 1, 2, …, with *N*_{t} = ()^{1/2}, *N*_{x} and *N*_{y} being the number of grid points in the zonal and meridional directions, and Δ*k* = 2π/(Δ*yN*_{y}). To make the main points, here, the vertically averaged over *z* ∈ [0, 22]-km horizontal kinetic energy spectra [*K*_{div}(*k*)] and the power spectra [*S*(*k*)] of the actual, balanced, and unbalanced divergence as given by the second-order Rossby number expansion are presented at day 8 when *K*_{WVD} [re^{(2)}] peaks for the M85 (Figs. 10a,c). The corresponding spectra at *z* = 11 km (Figs. 10b,d) are also presented to make it possible to visually compare the shape of the spectra with the horizontal divergence features seen at the same level in Figs. 4a and 4b.

The results shown in Fig. 10a indicate that for *k*/Δ*k* < 2 or in planetary scales (*λ* > 5000 km), while the kinetic energy spectra of the actual and unbalanced divergence coincide, the energy spectrum of the balanced divergence is about one order of magnitude smaller than that of the unbalanced component. There is hardly any sign of scale separation between the balanced and unbalanced flows. At the synoptic scales around *k*/Δ*k* = 3 (*λ* ~ 3000 km), which is associated with the scale of the BCW, the power spectra of balanced and unbalanced divergence reach the same peak value and both are smaller than that of the actual divergence. For 4 < *k*/Δ*k* < 15, the actual and unbalanced divergence spectra exhibit the same slow rate of decrease with increase in horizontal wavenumber. This shallow spectrum is then followed by a steep fall for *k*/Δ*k* > 40, which manifests the scales affected by numerical diffusion.

To make a better comparison between the wave–vortex decomposition and the HDA (see section 4) results, a high-pass spectral filter that retains wavelengths of smaller than 1000 km has been applied to the wave fields of the WVD by the second-order Rossby number expansion for the M85 case. The peak energy of the wave fields *K*_{WVD}[re^{(2)}], occurring at day 8, is estimated to be about 1.06 × 10^{4} J m^{−2} for the high-pass-filtered wave field. This amounts to nearly 26.4% reduction of the peak energy compared to the unfiltered wave field with 1.44 × 10^{4} J m^{−2} and brings the energy estimate of the WVD closer to the = 0.64 × 10^{4} J m^{−2} so that [re^{(2)}]/≈ 1.66. Therefore, near to 50% of the difference between [re^{(2)}] and is related to large scales for which better estimates are expected if the WVD procedures are carried out by inverting potential vorticity.

With regard to the broad maxima seen in both the vertically integrated and *z* = 11-km power spectra of the *δ* and *δ*_{unb}, as a rough estimate, one can cut the mesoscale power between *λ*_{min} = 200 km (*k*_{max} = π × 10^{−5} m^{−1}) and *λ*_{max} = 800 km (*k*_{min} = π/4 × 10^{−5} m^{−1}) and turn it into a boxcar spectrum for which one can write the density *S*_{k} = *s*^{2}/(*k*_{max} − *k*_{min}) with *s*^{2} = *S*_{k }*dk* the variance of the *δ*. In this way, the mesoscale part of the divergent kinetic energy is found to be equal to *K*_{div}(mso) = *s*^{2}/(2*k*_{min}*k*_{max}), which is almost the same for each of the *δ* and *δ*_{unb} fields. At the same time, one can show that the HDA sets the dominant wavenumber as the mean wavenumber [*k*_{mean} = (*k*_{max} + *k*_{min})/2 = 5π/8 × 10^{−5} m^{−1}] of the boxcar spectrum, from which we will have *K*_{div} (HDA) = *s*^{2}/(2) as the HDA estimate for the divergent kinetic energy and thus *K*_{div} (meso)/*K*_{div} (HDA) =/(*k*_{max}*k*_{min}) ≈ 1.56. This ratio is indeed close to the [re^{(2)}]/ ≈ 1.66 estimate obtained above by taking into consideration only scales smaller than 1000 km. Therefore, if one multiplies the *K*_{HDA} estimates by a factor because of the broad bandwidth of the mesoscale spectral maxima as opposed to a vanishing bandwidth originally assumed in the HDA, one can reach closer estimates from the HDA to that of the optimal WVD. The broad bandwidth is indeed consistent with the fine transverse scale of the wave packets seen in Fig. 4b.

From a comparison of Fig. 10a with Fig. 10b and Fig. 10c with Fig. 10d, it is also apparent that the spatial structure of IGWs may vary significantly with height. To examine that, the *K*_{WVD}[re^{(2)}] and *K*_{HDA} have been computed separately for the [0, 4]-, [4, 8]-, [8, 12]-, and [12, 22]-km layers (Fig. 11 and Table 3). In terms of the peaks ratio, [re^{(2)}]/, the maximum discrepancy between the two measures is observed in the layer containing the tropopause (i.e., in the source region where the sharpest horizontal and vertical gradients of potential vorticity form) leading to the broadband mesoscale spectra. In terms of the instantaneous ratio, the discrepancy decreases from the lower troposphere upward to the stratosphere where the wavelike character of the divergence field is well established. Further, the time lag between the stratospheric and tropospheric values of *K*_{WVD} [re^{(2)}] and *K*_{HDA} as well as the substantial reduction in the instantaneous values of the ratio {*K*_{WVD} [re^{(2)}]/*K*_{HDA}} (Fig. 11c) seen above the tropopause layer are clear signs of the dominance of vertical propagation over local generation in the stratosphere.

## 7. Concluding remarks

The main objective of the work presented was to determine the quantitative relation between estimates of the wave–vortex decomposition (WVD) and harmonic divergence analysis (HDA) methods for IGWs generated during the evolution of complex vortical flows involving jets, fronts, and convection. This is a problem of fundamental theoretical importance that may have practical implications for the parameterization of IGWs in general circulation models of the atmosphere, for example. The problem was addressed using the cyclonic life cycle of moist baroclinic waves by the WRF Model. Several experiments designed based on various degrees of moisture available in the initial state, from the fully dry to partially saturated, were performed to examine the sensitivity to moisture and moist convection as well. For WVD, instantaneous fields of vorticity were inverted using the balance relations of the first-order *δ*–*γ*, the Bolin–Charney, and a formal Rossby number expansion of first to third orders. In every case examined, the minimal measure of IGW activity came from the second-order Rossby number expansion. The minimal measure was regarded as being the optimal. Given the closeness of the best estimates from the three sets of balance relations as well as between the second- and the third-order Rossby number expansion, one can be confident that the space of balance relations has been sufficiently explored. For HDA’s application, results of a reference partitioning of the domain to boxes of 2000 × 2000 × 11 km^{3} volume were based for comparison, considering our desire to capture the characteristics of IGWs present in various parts of the domain. For the WVDs that can be retained as reasonable (excluding dg^{[2]}), we find that the different WVD estimates are, when integrated over the domain, within 30% of each other (as estimated from Fig. 6). This can then be contrasted with the factor-2–3 difference between HDA and WVD. Then again, as explained, the main cause of the difference between WVD and HDA estimates comes from the limitations of HDA (single wave packet) and the assumed scale separation. In a more general view, we cannot evaluate the quality or correctness of the different approaches to estimate imbalances and IGWs. For such an exercise, it would be necessary to define synthetic data fields where the “truth” is known by construction and to which the different methods are applied. Such a study is more complicated than can be included in the work presented. Here, instead, we aimed to interpret and understand the results for data coming from a particular setting.

Based on the peak values of the estimates for global kinetic energies of IGWs, , and , generally, the resulting / ratio varied between two and three with the lowest and highest values seen for the cases with the maximum initial relative humidities of 100% and 55% referred to by M100 and M55, respectively. While both the WVD and HDA estimates increase monotonically with increase in moisture content, the / ratio increases from the DRY to the M55 and then decreases more dramatically from the M55 to the M100. This nonlinear moisture dependence of the ratio reflects the fact that the validity of basic assumptions under which the HDA works (i.e., scale separation and dominance of wave packet–like structures in IGWs) is itself moisture dependent. This is supported by the observation that the M55 (M100) case exhibits the highest power in the unbalanced divergent kinetic energy at scales larger (smaller) than about 10^{3} km (see Fig. S8). Limiting the IGWs to subsynoptic scales (i.e., scales smaller than 10^{3} km), markedly smaller values for the / ratio can be obtained. For the M85 case, for example, the ratio decreases from 2.25, if all scales are considered in the WVD, to 1.66 if only scales smaller than 10^{3} km are treated by the WVD. Further, the broad bandwidth of the mesoscale maxima present in the actual and unbalanced divergence spectra suggests that one may include a factor of spectral broadening in the HDA. In this way, for the M85 case, by multiplying the energy estimates of the original HDA for IGWs by a factor of 1.7, one may bridge the gap between the WVD and HDA estimates at subsynoptic scales in these idealized experiments. In more practical situations, in different configurations and regimes of flow, the factor may be different. A main finding was that the WVD by the first-order Rossby number expansion leads to results near to optimal. For those global models that are equipped with fast elliptic equation solvers, to determine a flow-dependent spectral broadening factor, it should be feasible to carry out a first-order WVD every few time steps. It remains to see the effect of the resulting improvement of the HDA in the bulk parameterization schemes of the kind undertaken in Mirzaei et al. (2014), which we are currently working on.

## Acknowledgments

We thank the University of Tehran and IAP for providing support during this research. The studies of C. Z. were partly funded by the Deutsche Forschungsgemeinschaft through Grant ZU 120/2-1 for the research unit FOR 1898 (multiscale dynamics of gravity waves/spontaneous imbalance).

## REFERENCES

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

*Encyclopedia of Atmospheric Sciences*, 2nd ed. J. Pyle and F. Zhang, Eds., Vol. 2, Elsevier, 298–303.

## Footnotes

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175_JAS-D-16-0366.s1.

This article is included in the Multi-Scale Dynamics of Gravity Waves (MS-GWaves) Special Collection.

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