Waveguidability of idealized midlatitude jets and the limitations of ray tracing theory

Ray paths of stationary Rossby waves emanating from a local midlatitude source are usually refracted equatorward. However, this general tendency for equatorward propagation is mitigated by the presence of a midlatitude jet that acts as a zonal waveguide. This opens up the possibility of circum-global teleconnections and quasi-resonance, which suggests that the ability to guide a wave in the zonal direction is an important jet property. This paper investigates waveguidability of idealized midlatitude jets in a barotropic model on the sphere. A forced-dissipative model configuration with a local source for Rossby waves is used in order to quantify waveguidability by diagnosing the latitudinal distribution of waviness in a longitudinal sector far downstream of the forcing. Systematic sensitivity experiments show that waveguidability increases smoothly with increasing jet amplitude and with decreasing jet width. This result is contrasted with the predictions from two idealized theoretical concepts based on (1) ray tracing as derived from Wentzel–Kramers–Brillouin (WKB) theory and (2) a sharp jet with a zonally oriented front of potential vorticity. The existence of two so-called turning latitudes, which is the key diagnostic for a zonal waveguide according to ray tracing theory, turns out to be a poor predictor for the dependence of waveguidability on jet amplitude and jet width obtained in the numerical simulations. By contrast, the meridional gradient of potential vorticity correlates fairly well with the diagnosed waveguidability. The poor predictions from ray tracing are not surprising, because the underlying WKB assumptions are not satisfied in the current context. The failure of WKB is traced back to the properties of the underlying equations, and a heuristic argument is presented to elucidate the potential of the potential vorticity (PV) gradient to act as a proxy for waveguidability.

1 Introduction

Rossby waves are a ubiquitous feature of the atmospheric flow in the upper troposphere (Rossby, 1940; Rhines, 2002). They can transfer energy and momentum across large distances and sometimes give rise to teleconnections (Wallace and Gutzler, 1981; Hoskins and Karoly, 1981; Branstator, 2002). The present paper focuses on Rossby waves in midlatitudes, where they may occur in the form of Rossby wave packets (Wirth et al., 2018).

An important aspect in connection with midlatitude Rossby waves is the extent to which they are ducted in the zonal direction. As is well known, there is a general tendency for Rossby waves to be refracted equatorward owing to the sphericity of the Earth (Hoskins and Karoly, 1981; Hoskins and Ambrizzi, 1993); in practice, this holds true to the extent that the background flow varies smoothly. However, the existence of a zonal jet may change the situation and lead to preferential propagation in the zonal direction. The latter situation is often referred to as a zonal waveguide (Branstator, 2002; Schwierz et al., 2004; Martius et al., 2010; Branstator and Teng, 2017). An interesting question is as follows: which properties of the background atmosphere constitute a strong zonal waveguide? Earlier work of Manola et al. (2013) suggests that strong and narrow jets are associated with high waveguidability, although their work does not suggest a mechanistic explanation.

To the extent that the midlatitude background flow represents an efficient waveguide, this may lead to circumglobal Rossby waves, which in turn can result in circumglobal teleconnections (Branstator, 2002; Ding and Wang, 2005; Feldstein and Dayan, 2008; O'Reilly et al., 2018). If, in addition, such circumglobal Rossby waves turn quasi-stationary, this has two implications: first, weather is synchronized across long distances (Kornhuber et al., 2019), and, second, there is the possibility of constructive interference of the Rossby wave with itself, which may lead to quasi-resonance in the case of stationary forcing (Petoukhov et al., 2013; Kornhuber et al., 2017b). In fact, a number of previous authors argued that an increased tendency for quasi-resonance may be responsible for the increased occurrence of recent Northern Hemisphere weather extremes (Petoukhov et al., 2013, 2016; Coumou et al., 2014; Stadtherr et al., 2016; Kornhuber et al., 2017a; Mann et al., 2017). The issue is interesting because quasi-resonance would imply that small changes in the properties of the jet and its associated waveguidability can cause large changes in Rossby wave amplitude. This scenario provides additional motivation for a comprehensive understanding of waveguidability.

A key argument in the work of previous authors to support the theory of quasi-resonance draws on the ideas of ray tracing, which in turn are based on the Wentzel–Kramers–Brillouin (WKB) theory (Lighthill, 1967). Ray tracing allows one to diagnose the refraction of ray path along which the waves propagate. In this theory, a zonal jet turns into an efficient waveguide as soon as a region of wavelike propagation in the neighborhood of the jet is straddled by two so-called turning latitudes on either side of the jet, where a suitably defined “refractive index” turns zero. Not surprisingly, the existence of two turning latitudes plays a crucial role in the analysis of the above-quoted papers in their attempt to connect the observations with the theory of quasi-resonance (e.g., Fig. 4 in Petoukhov et al., 2013, Fig. 2 in Petoukhov et al., 2016, Fig. 1 in Kornhuber et al., 2017a).

The basic tenet of WKB theory is that the scale Δ_w of the wave must be much smaller than the scale Δ_bg on which the background flow varies, i.e.,

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{{\mathrm{\Delta }}_{\mathrm{w}}}{{\mathrm{\Delta }}_{\mathrm{bg}}}}\ll \mathrm{1}.\end{array}$$

Unfortunately, this assumption is often violated in connection with Rossby waves. For instance, in the analysis of Kornhuber et al. (2017a), the two turning latitudes are separated by about 10^∘ ≈ 1000 km (see their Fig. 1). Assuming that this distance corresponds to one-half wavelength in the meridional direction and that the associated scale is given by the wavelength divided by 2π, one obtains Δ_bg≈320 km; but at the same time the zonal wavelength of the waves in question (wavenumber 7 at midlatitudes) is on the order of 4000 km, which means that Δ_w≈640 km. In this situation, the relation (Eq. 1) is badly violated and one would not expect WKB theory to be applicable at all.

This state of affairs motivates the current work, in which we consider the waveguidability of midlatitude jets in an idealized modeling framework and investigate the validity of ray tracing arguments. Our approach is partly based on the work of Manola et al. (2013), but with modifications and extensions in a number of ways. Key to our analysis is a forced-dissipative model configuration, which allows us to explicitly simulate the propagation of stationary Rossby waves emanating from a local source. To the extent that the path of propagation is confined in the zonal direction, the scenario will be associated with a high degree of waveguidability. These results will be compared with the predictions from ray tracing theory. We will see that the existence of two turning latitudes is a poor predictor for waveguidability. Instead, the strength of the meridional gradient of potential vorticity turns out to be a better proxy for waveguidability in our framework. In addition, we will corroborate our results by analyzing the underlying equations, which is facilitated by the idealized nature of our model.

The paper is organized as follows. First in Sect. 2 we present our model, its configuration, and our method of numerical solution. Section 3 reviews the relevant theoretical concepts before our results are presented and compared with the theoretical concepts in Sect. 4. Section 5 seeks a deeper understanding of the numerical results through reference to the linearized equations, and finally Sect. 6 gives a short summary and our conclusions.

2 Model configuration and numerical solution

We consider non-divergent barotropic flow on a sphere in a forced-dissipative configuration. In our analysis we focus on the stationary part of the solution, which is obtained through temporal averaging. The key model variable is absolute vorticity q, which plays the role of potential vorticity (PV) in the barotropic model. It is given by

$$\begin{array}{}\text{(2)}& q(\mathit{\lambda },\mathit{\varphi },t)=\mathrm{2}\mathrm{\Omega }\sin \mathit{\varphi }+\mathit{\zeta },\end{array}$$

where

$$\begin{array}{}\text{(3)}& \mathit{\zeta }(\mathit{\lambda },\mathit{\varphi },t)={\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}{\displaystyle \frac{\partial v}{\partial \mathit{\lambda }}}-{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(u\cos \mathit{\varphi }\right)\end{array}$$

is relative vorticity, v=(u, v) is the horizontal wind, λ is longitude, ϕ is latitude, t is time, a is the Earth's radius, and Ω is the angular velocity of the Earth's rotation. The dynamics are determined by the following model equation

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{D}q}{\mathrm{D}t}}=-{\mathit{\lambda }}_{\mathrm{r}}\left(q-q_{\mathrm{0}}\right)+F,\end{array}$$

where $\mathrm{D}/\mathrm{D}t=\partial /\partial t+\mathit{v}\cdot \mathrm{\nabla }$ denotes the material derivative, λ_r is a damping rate,

$$\begin{array}{}\text{(5)}& q_{\mathrm{0}}\left(\mathit{\varphi }\right)=\mathrm{2}\mathrm{\Omega }\sin \mathit{\varphi }-{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}\mathit{\varphi }}}\left(u_{\mathrm{0}}\cos \mathit{\varphi }\right)\end{array}$$

is a zonally symmetric background PV (with u₀(ϕ) denoting the corresponding background zonal wind), and F represents the forcing. For later reference we note that

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\mathrm{d}{q}_{\mathrm{0}}}{\mathrm{d}\mathit{\varphi }}}=\mathit{\beta }\left(\mathit{\varphi }\right)-{\displaystyle \frac{\mathrm{d}}{a\mathrm{d}\mathit{\varphi }}}\left[{\displaystyle \frac{\mathrm{1}}{\cos \mathit{\varphi }}}{\displaystyle \frac{\mathrm{d}}{a\mathrm{d}\mathit{\varphi }}}\left(u_{\mathrm{0}}\cos \mathit{\varphi }\right)\right]\end{array}$$

with

$$\begin{array}{}\text{(7)}& \mathit{\beta }\left(\mathit{\varphi }\right)={\displaystyle \frac{\mathrm{2}\mathrm{\Omega }}{a}}\cos \mathit{\varphi }.\end{array}$$

The forcing is implemented as pseudo-orographic forcing with a rather local orography, i.e.,

$$\begin{array}{}\text{(8)}& F=\mathrm{div}\left(f_{\mathrm{0}}{\mathit{v}}_{\mathrm{f}}h\right)\end{array}$$

with $f_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$ s⁻¹, v_f=(u_f, 0),

$$\begin{array}{}\text{(9)}& u_{\mathrm{f}}\left(\mathit{\varphi }\right)=U_{\mathrm{f}}\cos \mathit{\varphi }\end{array}$$

(using U_f=15 m s⁻¹), and with a Gaussian-shaped orography

$$\begin{array}{}\text{(10)}& h(\mathit{\lambda },\mathit{\varphi })=h_{\mathrm{0}}\exp \left(-{\displaystyle \frac{{\left(\mathit{\lambda }-{\mathit{\lambda }}_{\mathrm{0}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{\lambda }}^{\mathrm{2}}}}-{\displaystyle \frac{{\left(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{\varphi }}^{\mathrm{2}}}}\right),\end{array}$$

where h₀=0.3, σ_λ=10^∘, σ_ϕ=10^∘, λ₀=30^∘ E, and ϕ₀=45^∘ N. Our formulation (Eq. 8) for the forcing has the advantage that it integrates to zero upon global integration, which means that it creates equal amounts of positive and negative local PV anomalies (in contrast with the forcing used by Manola et al., 2013). Note also that our forcing F contains a fixed and specified flow field v_f and is, therefore, independent of the actual flow v. The forcing is localized in the sense that the spatial scale of the pseudo-orography h is much smaller than the planetary scale.

The dissipative term in Eq. (4) is designed such that in the unforced case (F=0) the background state q₀(ϕ) is a stationary solution of the equation. Waves, i.e., deviations from zonal symmetry, arise to the extent that F is nonzero. But even for nonzero forcing the solution is gradually relaxed back towards the zonally symmetric background state. As relaxation parameter we use λ_r=(7 d)⁻¹, which was chosen such that Rossby wave packets emanating from the pseudo-orography are just about able to propagate once around the Earth before they get completely dissipated. This facilitates the interpretation of our solutions because it deliberately avoids constructive interference of a Rossby wave train with itself and, hence, the phenomenon of quasi-resonance. It contrasts our approach from that of Manola et al. (2013), who used a slightly smaller value for λ_r and obtained weak resonance phenomena.

Any deviation from purely zonal flow will be considered to be an eddy or, more specifically, a wave. Correspondingly, we will use maps of the meridional wind v in order to visualize the waves.

The flow is initialized with the zonally symmetric background state. The latter is defined by specifying a zonally symmetric zonal flow field u₀(ϕ) and computing the corresponding q₀(ϕ) via Eq. (2). Our main strategy includes a set of sensitivity experiments in which we impose different background flows u₀(ϕ) and investigate their impact on the stationary part of the solution. More specifically, the background flow is defined as a superposition of a solid-body rotation u_SB(ϕ), a midlatitude jet U_J(ϕ), and a small linear component L(ϕ), i.e.,

$$\begin{array}{}\text{(11)}& u_{\mathrm{0}}\left(\mathit{\varphi }\right)=u_{\mathrm{SB}}\left(\mathit{\varphi }\right)+u_{\mathrm{J}}\left(\mathit{\varphi }\right)+L\left(\mathit{\varphi }\right).\end{array}$$

The solid-body rotation part is given by

$$\begin{array}{}\text{(12)}& u_{\mathrm{SB}}\left(\mathit{\varphi }\right)=U_{\mathrm{SB}}\cos \mathit{\varphi }\end{array}$$

with U_SB=15 m s⁻¹, and the jet part is given by a Gaussian profile

$$\begin{array}{}\text{(13)}& u_{\mathrm{J}}\left(\mathit{\varphi }\right)=U_{\mathrm{J}}{e}^{-\frac{{\left(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{J}}^{\mathrm{2}}}},\end{array}$$

with ϕ₀=45^∘ N. The latter choice is motivated by the work of Manola et al. (2013), who showed that a Gaussian profile fits the shape of observed jets fairly well. The two parameters U_J and σ_J characterize the amplitude and width of the jet, respectively. The value of U_J will be varied between 0 and 40 m s⁻¹. For the width parameter we choose σ=5^∘ unless stated otherwise, and we will refer to this choice as a “narrow” jet. Generally, the sum of u_SB(ϕ) and u_J(ϕ) would be nonzero at the poles, which is undesirable. We, therefore, included the linear function $L\left(\mathit{\varphi }\right)=a+b\mathit{\varphi }$ on the right-hand side of Eq. (11) with the parameters a and b chosen such that the resulting total profile u₀(ϕ) turns zero at both poles. Quite deliberately we choose our background wind to be positive for all latitudes in order to allow stationary wave propagation. As explained below, this enables us to diagnose the waveguidability of the jet in a meaningful manner. For later reference we note that

$$\begin{array}{}\text{(14)}& {\left.{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{u}_{\mathrm{J}}}{\mathrm{d}{\mathit{\varphi }}^{\mathrm{2}}}}\right|}_{\mathit{\varphi }={\mathit{\varphi }}_{\mathrm{J}}}=-{\displaystyle \frac{U_{\mathrm{J}}}{{\mathit{\sigma }}_{\mathrm{J}}^{\mathrm{2}}}}.\end{array}$$

Numerical solutions are obtained using a standard pseudo-spectral scheme. To ensure numerical stability, we add a hyperdiffusion term

$$\begin{array}{}\text{(15)}& T_{\mathrm{h}}=D{\mathrm{\nabla }}^{\mathrm{4}}\left(q-q_{\mathrm{0}}\right)\end{array}$$

to the right-hand side of Eq. (4) with D=10¹⁵ m⁴ s⁻². The term T_h is generally small in comparison with the other terms in Eq. (4) and will be neglected in our theoretical analysis. The code uses triangular truncation at wavenumber N=72. Time-stepping is done using the leapfrog scheme with a time step of Δt=20 min, combined with a Robert–Asselin filter in order to damp the computational mode.

For some choices of the background wind our numerical solutions turn out to be highly transient. This transience is due to barotropic instability in the case of a strong narrow jet. In all our experiments these unstable modes are characterized by a large eastward phase velocity. This is illustrated in Fig. 1 showing a Hovmöller diagram of the transient part of the meridional wind, $v^{\prime }=v-\stackrel{\mathrm{‾}}{v}$, where the overbar denotes the time average. The Hovmöller diagram indicates how barotropic instability is triggered downstream of the orography (which is located at λ₀=30^∘). During the first 10 d the instability materializes in the form of an eastward-propagating Rossby wave packet with nonzero (eastward) phase speed. Later this wave packet disperses in the zonal direction, and after about day 30 it has turned into an almost pure sine wave with eastward phase velocity. This behavior allows us to eliminate the transient part of the solution by simple time averaging. More specifically, for each model run we performed a 100 d long model integration and extracted the stationary part $\stackrel{\mathrm{‾}}{v}$ of the flow by averaging over the last 90 d (the initial 10 d were discarded in order to eliminate initial transients).

Figure 1 Development of a barotropically unstable initial state for a model setup with a strong narrow jet at 45^∘ N. The plot shows a Hovmöller diagram of the meridional wind component of the transient part of the flow, $v^{\prime }=v-\stackrel{\mathrm{‾}}{v}$, averaged over 30–60^∘ N.

For illustration we consider two very different background flows, one with pure solid-body rotation (corresponding to U_J=0 m s⁻¹) and one with a strong narrow jet superimposed (corresponding to U_J=40 m s⁻¹). The two latitudinal profiles U_J(ϕ) are depicted in Fig. 2. The corresponding numerical solutions for $\stackrel{\mathrm{‾}}{v}$ are presented in Fig. 3. In the case of the solid-body rotation, the transients have died out after about 10 d such that the solution for v at later times (not shown) is practically indistinguishable from the time average $\stackrel{\mathrm{‾}}{v}$ displayed in Fig. 3a. By contrast, in the strong jet case there is strong transience throughout the integration, but these transients are effectively eliminated by the time averaging, leaving us with the forced stationary part $\stackrel{\mathrm{‾}}{v}$ displayed in Fig. 3b.

Figure 2 Latitudinal profiles of the background zonal wind u₀(ϕ) for the strong narrow jet case (solid line) and for the pure solid-body rotation case (dashed line).

Figure 3 Illustration of the numerical solution for an initial state with (a) pure solid-body rotation and (b) a strong narrow jet superimposed on solid-body rotation. Displayed in both cases is the time-averaged meridional wind $\stackrel{\mathrm{‾}}{v}$. The oval-shaped yellow contour is the contour of the orography h(λ, ϕ) at 0.3 h₀.

Apparently, for pure solid-body rotation (Fig. 3a) there is a wave train emanating from the pseudo-orography; downstream of the source the individual troughs and ridges develop a strong NE–SW tilt and the wave train crosses the Equator and disperses into the Southern Hemisphere (Hoskins et al., 1977; Hoskins and Karoly, 1981). At the same time, the wave signal is damped with increasing distance from the wave source, which is due to our relaxation term in Eq. (4). The cross-equatorial propagation implies that some 180^∘ downstream of the Rossby wave source the wave train is found in the Southern Hemisphere, and the Northern Hemisphere is practically devoid of any wave signal at these longitudes.

By contrast, the solution for the strong jet case (Fig. 3b) indicates that the majority of the wave signal remains in the Northern Hemisphere midlatitudes, and only a rather small fraction of the wave signal follows a great circle into the Southern Hemisphere. This behavior is consistent with the notion that a strong narrow jet acts as a waveguide (Manola et al., 2013). As a consequence, in this case the majority of the wave signal 180^∘ downstream of the Rossby wave source has remained in the Northern Hemisphere midlatitudes. Later in Sect. 4.1 we will use the striking difference in behavior between these two contrasting scenarios to define a quantitative measure for waveguidability.

3 Theoretical concepts

Before we do so, however, in this section we review some well-known concepts for the analysis of stationary Rossby waves and their propagation in a spherical domain. This will help us to interpret our results in the later parts of the paper.

First, we linearize the model equation (Eq. 4) about the background state q₀(ϕ) and obtain the following equation for the (small) perturbation $q^{\prime }=q-q_{\mathrm{0}}$,

$$\begin{array}{}\text{(16)}& \left({\displaystyle \frac{\partial }{\partial t}}+u_{\mathrm{0}}{\displaystyle \frac{\partial }{a\cos \mathit{\varphi }\partial \mathit{\lambda }}}\right)q^{\prime }+v^{\prime }{\displaystyle \frac{\mathrm{d}{q}_{\mathrm{0}}}{a\mathrm{d}\mathit{\varphi }}}=S^{\prime },\end{array}$$

where

$$\begin{array}{}\text{(17)}& S^{\prime }=-{\mathit{\lambda }}_{\mathrm{r}}{q}^{\prime }+F.\end{array}$$

Note that the forcing F is now assumed to be small and, hence, a perturbation term. Expressing, as usual (see, e.g., Vallis, 2006), the perturbation variables q^′ and v^′ in terms of a perturbation streamfunction ψ^′, this becomes

$$\begin{array}{}\text{(18)}& \left({\displaystyle \frac{\partial }{\partial t}}+u_{\mathrm{0}}{\displaystyle \frac{\partial }{a\cos \mathit{\varphi }\partial \mathit{\lambda }}}\right){\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\psi }}^{\prime }+{\displaystyle \frac{\partial {\mathit{\psi }}^{\prime }}{a\cos \mathit{\varphi }\partial \mathit{\lambda }}}{\displaystyle \frac{\mathrm{d}{q}_{\mathrm{0}}}{a\mathrm{d}\mathit{\varphi }}}=S^{\prime },\end{array}$$

where ∇² is the two-dimensional Laplace operator in spherical coordinates.

For further progress, it turns out it is convenient to perform a coordinate transformation corresponding to a Mercator projection of the sphere onto a plane (see Hoskins and Karoly, 1981; Hoskins and Ambrizzi, 1993). This is achieved by defining

$$\begin{array}{}\text{(19)}& x=a\mathit{\lambda },\end{array}$$

$$\begin{array}{}\text{(20)}& y=a\ln \left({\displaystyle \frac{\mathrm{1}+\sin \mathit{\varphi }}{\cos \mathit{\varphi }}}\right).\end{array}$$

This transformation is valid everywhere except at the poles; it transforms Eq. (18) into

$$\begin{array}{}\text{(21)}& \left({\displaystyle \frac{\partial }{\partial t}}+u_{\mathrm{M}}{\displaystyle \frac{\partial }{\partial x}}\right)\left({\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\psi }}^{\prime }}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\psi }}^{\prime }}{\partial y^{\mathrm{2}}}}\right)+{\mathit{\beta }}_{\mathrm{M}}{\displaystyle \frac{\partial {\mathit{\psi }}^{\prime }}{\partial x}}=S^{\prime },\end{array}$$

where

$$\begin{array}{}\text{(22)}& u_{\mathrm{M}}={\displaystyle \frac{u_{\mathrm{0}}}{\cos \mathit{\varphi }}}\end{array}$$

and

$$\begin{array}{}\text{(23)}& {\mathit{\beta }}_{\mathrm{M}}=\cos \mathit{\varphi }{\displaystyle \frac{\mathrm{d}{q}_{\mathrm{0}}}{a\mathrm{d}\mathit{\varphi }}}\equiv {\displaystyle \frac{\mathrm{d}{q}_{\mathrm{0}}}{\mathrm{d}y}}.\end{array}$$

Assuming u_M and β_M to be constants, one can look for plane wave solutions of the form

$$\begin{array}{}\text{(24)}& {\mathit{\psi }}^{\prime }(x,y,t)=\widehat{\mathit{\psi }}{e}^{i(kx+ly-\mathit{\omega }t)},\end{array}$$

where k and l are the wavenumbers in the zonal and meridional directions on the Mercator projection, respectively, and ω is the angular velocity. Note that k relates to the dimensionless zonal wavenumber s through

$$\begin{array}{}\text{(25)}& s=ka.\end{array}$$

Perturbations of the form in Eq. (24) are solutions of Eq. (21) if the following dispersion relation is satisfied:

$$\begin{array}{}\text{(26)}& \mathit{\omega }=u_{\mathrm{M}}k-{\displaystyle \frac{k{\mathit{\beta }}_{\mathrm{M}}}{k^{\mathrm{2}}+l^{\mathrm{2}}}}.\end{array}$$

Obviously, the advantage of using the new coordinates (x, y) is that both the Eq. (21) and the associated dispersion relation (Eq. 26) are formally identical to their Cartesian β-plane version, if the zonal basic flow u₀ is replaced by u_M and the planetary vorticity gradient β is replaced by β_M.

For later reference we note that for a background flow with pure solid-body rotation as in Eq. (12), one obtains $u_{\mathrm{M}}=U_{\mathrm{SB}}=\mathrm{const}$ and

$$\begin{array}{}\text{(27)}& {\mathit{\beta }}_{\mathrm{M}}=\mathrm{2}\left(\mathrm{\Omega }+{\displaystyle \frac{U_{\mathrm{SB}}}{a}}\right){\displaystyle \frac{{\cos }^{\mathrm{2}}\mathit{\varphi }}{a}}.\end{array}$$

The latter is obviously not a constant, which, hence, violates the assumption made above when deriving the dispersion relation.

Figure 4 Diagnosing the zonally symmetric background state for three different cases: solid-body rotation (a–c), a weak narrow jet with U_J=10 m s⁻¹ (d–f), and a strong narrow jet with U_J=40 m s⁻¹ (g–i). Panels (a, d, g) show the meridional wind profile u₀(ϕ), (b, e, h) show β_M(ϕ), and (c, f, i) show the dimensionless stationary wavenumber ${\widehat{K}}_{\mathrm{s}}\left(\mathit{\varphi }\right)$, where the negative values (shading) represent minus the imaginary part of ${\widehat{K}}_{\mathrm{s}}$.

3.1 WKB approximation and ray tracing

If the factors u_M(ϕ) and β_M(ϕ) in Eq. (21) depend on latitude instead of being constants, further analytical progress is possible by assuming that the spatial variation in u_M and β_M is small in the sense that their derivatives can be neglected in comparison with the rate of change resulting from the wave's phase. This method is often referred to as WKB approximation (see, e.g., Lighthill, 1967). In this framework, a plane wave of the form in Eq. (24) can still be considered to be an approximate solution of Eq. (21). For the special case of stationary waves, this immediately leads to

$$\begin{array}{}\text{(28)}& k^{\mathrm{2}}+l^{\mathrm{2}}={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{M}}}{u_{\mathrm{M}}}}\end{array}$$

or, using Eq. (25),

$$\begin{array}{}\text{(29)}& a^{\mathrm{2}}{l}^{\mathrm{2}}={\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}-s^{\mathrm{2}},\end{array}$$

where the square of the dimensionless stationary wavenumber is defined as

$$\begin{array}{}\text{(30)}& {\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}:=a^{\mathrm{2}}{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{M}}}{u_{\mathrm{M}}}}.\end{array}$$

Note that ${\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}$ inherits the property of “weak variation” from u_M and β_M. Apparently, for a given zonal wavenumber s, the solution (Eq. 24) is wavelike in those regions where the meridional wavenumber l²>0. In these regions the Rossby wave packets are free to propagate. For illustration, Fig. 4 shows profiles of u₀, β_M, and ${\widehat{K}}_{\mathrm{s}}=({\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{2}}$ for pure solid-body rotation (panels a, b, c), for the strong narrow jet (panels g, h, i), and for a weaker jet (panels d, e, f) representing a scenario intermediate between the first two. The profiles of ${\widehat{K}}_{\mathrm{s}}$ in Fig. 4 indicate that for solid-body rotation wave propagation in the WKB sense is possible over most of the sphere except in the regions close to the poles. In contrast, for background flows with a jet, both flanks of the jet are characterized by a strip of relatively low (Fig. 4f) or even imaginary (Fig. 4i) values of ${\widehat{K}}_{\mathrm{s}}$; depending on the zonal wavenumber, this potentially prohibits wave propagation in these strips.

Ray tracing arguments consider the propagation of Rossby waves in those regions which allow wave propagation (e.g., Hoskins and Karoly, 1981; Hoskins and Ambrizzi, 1993). It can be shown that the vector of the group velocity is refracted towards larger values of ${\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}$ as one proceeds in the direction of the propagation; in addition, latitudes where l=0 are associated with a purely zonal group velocity, which implies that ray paths return back to latitudes where the wave packet came from. The latter is only valid as long as the background state does not vanish at those latitudes, which is satisfied for all background states used in our study. Latitudes where l(ϕ) turns zero are aptly called “turning latitudes”.

Turning latitudes can easily be diagnosed from profiles of ${\widehat{K}}_{\mathrm{s}}$ as those latitudes where the function ${\widehat{K}}_{\mathrm{s}}\left(\mathit{\varphi }\right)$ intersects the straight line s=const. Let us, for illustration, consider a wave with s=4. Figure 4c then indicates that for solid-body rotation there are two turning latitudes, one at 60^∘ N and one at 60^∘ S. This is consistent with our numerical solution from Fig. 3a which shows wavelike behavior throughout most of the domain. By contrast, both cases with a jet superimposed (Fig. 4f and i) indicate that waves emanating from a source at 45^∘ N encounter two turning latitudes, one at 50^∘ N and another one at 35^∘ N. We note in passing that the strong jet case displayed in Fig. 4g, h, i features regions with a negative meridional gradient of absolute vorticity and, hence, of β_M (Fig. 4h). This opens up the possibility of barotropic instability (Charney and Stern, 1962), which we do observe in our numerical simulations.

The jet scenario is particularly interesting in the present context. Whenever there are two turning latitudes (l=0) separated by a region of wavelike propagation (l>0), ray tracing arguments predict that Rossby wave packets oscillate between the two turning latitudes as they propagate eastward. This means that they are effectively trapped between the two turning latitudes and, thus, ducted in the zonal direction. In other words, ray tracing theory suggests that the existence of two turning latitudes is tantamount to a perfect zonal waveguide. In fact, this seems to be broadly consistent with our numerical solution shown in Fig. 3b, where the majority of the wave signal is ducted in the zonal direction: the corresponding function ${\widehat{K}}_{\mathrm{s}}\left(\mathit{\varphi }\right)$ shows the existence of two turning latitudes for any wavenumber $\mathrm{1}\le s\le \mathrm{7}$, straddling a wavelike region in the core of the jet.

Let us return to the case of a background flow with pure solid-body rotation according to Eq. (12), yielding

$$\begin{array}{}\text{(31)}& {\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}\left(\mathit{\varphi }\right)=\mathrm{2}\left({\displaystyle \frac{\mathrm{\Omega }a}{U_{\mathrm{SB}}}}+\mathrm{1}\right){\cos }^{\mathrm{2}}\mathit{\varphi }\end{array}$$

(see Fig. 4c). For this case, ray tracing theory predicts that ray paths are generally refracted equatorward, because the function ${\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}\left(\mathit{\varphi }\right)$ maximizes right at the Equator. In fact, Hoskins and Karoly (1981) showed explicitly that ray paths on a solid-body rotation background flow are identical to great circles. In addition, ${\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}$ goes to zero at both poles such that at some point the expression on the right-hand side of Eq. (29) turns negative, leading to a turning latitude close to the pole. The crucial question in our context is whether the addition of a jet leads to a second (more equatorward) turning latitude which prevents the general equatorward refraction and forces the ray path back into the zonal direction. The latter scenario would be equivalent to a perfect zonal waveguide. It is this prediction from ray tracing theory that led previous authors to systematically search for the occurrence of two turning latitudes straddling a region of wave propagation. Interestingly, this prediction is not fully supported by our numerical simulations, as we will show in Sect. 4.

3.2 PV front

Another school of thought associates a midlatitude zonal waveguide with the existence of a sharp meridional gradient of PV (Schwierz et al., 2004; Martius et al., 2010). A highly idealized model that represents this idea would be a zonally oriented front of PV in the basic state separating two regions with a completely homogeneous PV distribution (Platzmann, 1949). Note that the idealization made in this model is opposite to that in the WKB approximation: a PV front is equivalent to a discontinuous jump in the background PV, while the WKB approximation assumes the background state to vary very gently.

The background flow associated with a PV discontinuity is a westerly jet with a cusp-like peak at the latitude of the discontinuity. Solutions of the linearized equations can be found which are wavelike in the zonal direction and evanescent in the meridional direction like

$$\begin{array}{}\text{(32)}& {\mathit{\psi }}^{\prime }=\widehat{\mathit{\psi }}{e}^{-\mathrm{sgn}\left(y-y_{\mathrm{0}}\right)k\left(y-y_{\mathrm{0}}\right)}{e}^{i(kx-\mathit{\omega }t)},\end{array}$$

where we assumed Cartesian geometry for simplicity, where sgn(…) denotes the sign function and where the other symbols have their usual meaning. The resulting dispersion relation is similar (albeit not identical) to the dispersion relation of Rossby waves on a β plane with constant background flow (see, e.g., Schwierz et al., 2004); in particular, both kinds of waves have westward phase propagation with respect to the basic flow. More importantly, however, these interfacial waves can only propagate in the zonal direction, and their amplitude is exponentially damped away from the PV front. This is consistent with the concept of the key importance of PV gradients for the existence of Rossby waves. It transpires that in this framework a zonal waveguide is present to the extent that the PV distribution has a zonally oriented sharp PV gradient separating two regions with near-homogeneous PV on both sides of the sharp gradient. As we will see, our numerical solutions turn out to be broadly consistent with this idea.

4 Investigating waveguidability

Based on the theoretical concepts sketched out in the previous section, we now define waveguidability in the framework of our numerical model, explore it systematically for various background flows, and compare the results with predictions from ray tracing theory.

4.1 Definition of waveguidability

As we argued before, background solid-body rotation can be considered a reference scenario in which the waves propagate along great circles from one hemisphere to the other. As a consequence, wave activity emanating from a local northern hemispheric Rossby wave source can be expected to be located almost entirely in the Southern Hemisphere some 180^∘ downstream (in longitude) of the source region. We interpret this scenario as the lack of waveguidability, and our diagnostic (to be defined shortly) should reflect this by producing a very small value.

These considerations motivated us to use the following method to quantify waveguidability. Introducing wave enstrophy of the stationary part of the solution as

$$\begin{array}{}\text{(33)}& \mathcal{E}(\mathit{\lambda },\mathit{\varphi })={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\stackrel{\mathrm{‾}}{q}}^{*}\right)}^{\mathrm{2}},\end{array}$$

where $q^{*}=q-\left[q\right]$ is the deviation of absolute vorticity from the zonal mean and the overline denotes the time average, we define a probability density function as

$$\begin{array}{}\text{(34)}& P\left(\mathit{\varphi }\right)={\displaystyle \frac{\cos \mathit{\varphi }}{N}}\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}\mathcal{E}(\mathit{\lambda },\mathit{\varphi })\mathrm{d}\mathit{\lambda },\end{array}$$

where a suitable downstream sector is defined through ${\mathit{\lambda }}_{\mathrm{1}}\le \mathit{\lambda }\le {\mathit{\lambda }}_{\mathrm{2}}$ with λ₁=180^∘ and λ₂=270^∘ and where N represents a normalization factor to guarantee that $\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}P\left(\mathit{\varphi }\right)\mathrm{d}\mathit{\varphi }=\mathrm{1}$. Note that our definition of P(ϕ) includes a factor cos ϕ in order to account for the variation in surface area with latitude. In other words, the function P(ϕ) quantifies the likelihood of encountering wave enstrophy ℰ as a function of latitude in the downstream sector [λ₁, λ₂]. Figure 5 shows this probability density function for our two background profiles from Fig. 2. Apparently, for solid-body rotation most of the wave enstrophy in the downstream sector is found in the southern midlatitudes, consistent with Fig. 3a. On the other hand, for the strong narrow jet, a large fraction of the wave enstrophy remains in the northern midlatitudes, which is consistent with Fig. 3b.

Figure 5 Probability density function P(ϕ) as defined in Eq. (34) for two different background flows: (a) solid-body rotation and (b) a strong narrow jet with U_J=40 m s⁻¹. The grey shading indicates the northern midlatitudes, where the forcing is located. Broadly speaking, the function P(ϕ) quantifies where the “wavines” can be found far downstream of the Rossby wave source.

The probability density function P(ϕ) is now used to define waveguidability 𝒲 as the probability of encountering downstream wave enstrophy in the northern midlatitudes at ${\mathit{\varphi }}_{\mathrm{1}}\le \mathit{\varphi }\le {\mathit{\varphi }}_{\mathrm{2}}$, i.e.,

$$\begin{array}{}\text{(35)}& \mathcal{W}=\underset{{\mathit{\varphi }}_{\mathrm{1}}}{\overset{{\mathit{\varphi }}_{\mathrm{2}}}{\int }}P\left({\mathit{\varphi }}^{\prime }\right)\mathrm{d}{\mathit{\varphi }}^{\prime },\end{array}$$

where we set ${\mathit{\varphi }}_{\mathrm{1}}=\mathrm{30}{}^{\circ }$ N and ϕ₂=60^∘ N. By design, 𝒲 is a number between 0 and 1, and it will be expressed in percent in the following. For the case of solid-body rotation (Fig. 5a), we obtain 𝒲=7.5 %, indicating that this background flow has a very low capacity to duct waviness in the zonal direction. On the other hand, for the strong narrow jet (Fig. 5b) we obtain 𝒲=75.4 %; although this value is still less than 100 %, it is an order of magnitude larger than for pure solid-body rotation.

Our method to quantify waveguidability is based on, but also extends, the ideas from Manola et al. (2013). A key difference is that we diagnose the meridional distribution of waviness only in a downstream longitudinal sector rather than circumglobally, which allows us to associate a percentage (i.e., a number between 0 and 100) with the waveguidability of any jet. Nevertheless, our definition of waveguidability is neither universal nor exact. It is not universal because it only makes sense in our particular model configuration where we consider the stationary response from a local Rossby wave source in a forced-dissipative framework; it is not exact because it depends on model details like, for example, the exact value of the damping parameter λ_r, the choice of the boundaries for the downstream region [λ₁, λ₂], and the definition of northern midlatitudes [ϕ₁, ϕ₂]. Nevertheless, this diagnostic is appropriate for our purpose, because we will only study relative changes of 𝒲 as we change the background flow.

4.2 Waveguidability as a function of jet amplitude and jet width

We now systematically vary the amplitude U_J and the width σ_J of our background jet (Eq. 13) and use our numerical solutions to compute the value of 𝒲 as explained above. The result is shown in Fig. 6 (red squares and thin red line). Apparently, there is a very smooth and gradual variation in 𝒲 with both U_J and σ_J extending from rather low values around 10 % to rather high values close to 80 % (within the ranges considered for U_J and σ_J). Qualitatively the behavior is entirely consistent with the results of Manola et al. (2013), who suggested increasing waveguidability for increasing jet amplitude and decreasing jet width.

Figure 6 Dependence of waveguidability 𝒲 on (a) jet strength U_J (with σ_J=5^∘ kept fixed) and (b) jet width σ_J (with U_J=40 m s⁻¹ kept fixed). The red squares (connected by a thin red line) represent the values of 𝒲 diagnosed from the numerical solutions. The solid blue line represents the prediction from ray tracing theory using zonal wavenumber s=4. The blue shaded area indicates the uncertainty of the ray tracing prediction associated with the fact that the numerical solution does not have a fixed single zonal wavenumber (see text for details).

We tested the sensitivity of our results with respect to the particular metric (Eq. 33) for quantifying waviness of the stationary eddies. To this end we repeated the calculations, but using eddy kinetic energy,

$$\begin{array}{}\text{(36)}& \mathcal{K}(\mathit{\lambda },\mathit{\varphi })={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{{\left({\stackrel{\mathrm{‾}}{u}}^{*}\right)}^{\mathrm{2}}+{\left({\stackrel{\mathrm{‾}}{v}}^{*}\right)}^{\mathrm{2}}\right\},\end{array}$$

as an alternative metric. The corresponding dependence of waveguidability 𝒲 on jet amplitude and jet width is very similar to the behavior shown in Fig. 6, suggesting that our results do not sensitively depend on the choice of the metric.

4.3 Comparison with theoretical expectations

The numerical results from Fig. 6 are now contrasted with the predictions from theoretical concepts that we sketched in the previous section. We start with ray tracing theory, in which the profile of ${\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}$ is the key diagnostic and the existence of two turning latitudes indicates a perfect waveguide. The difficulty in applying this theory is that it is not clear, a priori, which zonal wavenumber to consider as the relevant one. Our numerical solution generally indicates an entire range of wavenumbers (Fig. 7). For the background profiles with solid-body rotation and for weak jets, the Fourier amplitudes maximize at s=4 (Fig. 7a and b), while for the stronger and narrower jets the Fourier amplitudes maximize at s=3 (Fig. 7c). Taking s=4 as the relevant zonal wavenumber, we checked whether there are two turning latitudes bounding a wavelike region within the jet. Whenever this turned out to be the case, we associated this jet with the value 𝒲=100 %, otherwise we set 𝒲=0 %. This literal interpretation of ray tracing theory is depicted as the solid blue line in Fig. 6. It is clear that the resulting discontinuous dependence of 𝒲 on U_J and σ_J (thick blue line) is in stark contrast with the very gradual behavior obtained in the numerical solutions (red squares and thin red line). In addition, waveguidability diagnosed from the numerical solutions is much less than 100 % in a large part of the range of parameters where ray tracing theory predicts a perfect waveguide. For instance, a narrow jet with U_J≈11 m s⁻¹ yields 𝒲≈26 % in our simulations – rather than 100 % as predicted by ray tracing theory.

Figure 7 Zonal Fourier spectrum of $\stackrel{\mathrm{‾}}{v}$ averaged over 30–60^∘ N, for (a) solid-body rotation, (b) a weak narrow jet with U_J=10 m s⁻¹, and (c) a strong narrow jet with U_J=40 m s⁻¹.

One could argue that the above interpretation of ray tracing theory is too brute force because it does not account for the existence of an entire spectrum of zonal wavenumbers (see Fig. 7). We, therefore, considered another extreme scenario by assuming that all zonal wavenumbers have equal amplitude. The value of 𝒲 for this scenario was obtained as the fraction of the number of wavenumbers with two turning latitudes divided by the number of all wavenumbers allowing wavelike behavior within the jet region. The resulting values are represented as the other edge of the blue-shaded area in Fig. 6. Apparently, there is now a more gradual transition from 𝒲=0 % to 𝒲=100 % as jet amplitude and jet width are varied, but the transition is still significantly steeper than for the values obtained from the simulations. In addition, the scenario with equal amplitude for all relevant zonal wavenumbers is not realistic either, because even for very weak jets our numerical solutions indicate a maximum spectral amplitude as some intermediate wavenumber. In the end, the best estimate for a “fair” prediction from ray tracing theory is a line which is located somewhere in the middle of the blue shaded area in Fig. 6. This yields a fairly steep (albeit not discontinuous) transition in both panels of the figure, which – again – is in stark contrast with the very gradual curves obtained from the simulations. But even for jets where the ray tracing interpretation is unambiguous, its prediction is still far off the result of our simulations: for the narrow jet with U_J=20 m s⁻¹, ray tracing theory unambiguously predicts 𝒲=100 %, while our analysis suggests a much smaller value of 𝒲≈50 %.

It is enlightening to focus on the behavior in Fig. 6a for 20 m s⁻¹ ≤ U_J ≤ 40 m s⁻¹. Within this range, the numerical solutions indicate an increase in 𝒲 from 49 % to 75 %. On the other hand, the ray tracing prediction is 100 % for all those jets, and this prediction is unambiguous because there are two turning latitudes for any zonal wavenumber in question. This behavior is due to the fact that within this range of U_J values the profiles of ${\widehat{K}}_{\mathrm{s}}\left(\mathit{\varphi }\right)$ are practically independent of U_J. This, in turn, is related to the fact that for narrow strong jets the meridional gradient of absolute vorticity is dominated by the meridional curvature of the background wind field. In this case the meridional gradient of background PV (Eq. 6) can be approximated as

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{\mathrm{d}{q}_{\mathrm{0}}}{a\mathrm{d}\mathit{\varphi }}}\approx -{\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{u}_{\mathrm{0}}}{\mathrm{d}{\mathit{\varphi }}^{\mathrm{2}}}},\end{array}$$

and, using Eqs. (22), (23), and (30), one obtains

$$\begin{array}{}\text{(38)}& {\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}\approx -{\displaystyle \frac{{\cos }^{\mathrm{2}}\mathit{\varphi }}{u_{\mathrm{0}}}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{u}_{\mathrm{0}}}{\mathrm{d}{\mathit{\varphi }}^{\mathrm{2}}}}.\end{array}$$

For our Gaussian jet with a fixed width σ_J, this together with Eq. (14) yields

$$\begin{array}{}\text{(39)}& {\left.{\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}\right|}_{{\mathit{\varphi }}_{\mathrm{J}}}\approx +{\displaystyle \frac{{\cos }^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{J}}}{{\mathit{\sigma }}_{\mathrm{J}}^{\mathrm{2}}}}\end{array}$$

at the jet latitude ϕ_J, which is apparently independent of jet amplitude U_J. This argument explains why ray tracing theory is unable to predict the substantial increase in 𝒲 when increasing the jet amplitude from 20 to 40 m s⁻¹. We conclude that ray tracing theory does a poor job in predicting the variability of waveguidability as we vary either the jet amplitude or the jet width.

How does the other theoretical concept, in which the meridional gradient of PV plays the key role for waveguidability, fare? To the extent that the PV gradient is a relevant proxy for waveguidability, we expect a close correlation between 𝒲 and β_M at jet latitude. For this reason, we plot in Fig. 8 the values of 𝒲 together with the values of β_M at ϕ_J after we non-dimensionalized β_M by multiplying with $f_{\mathrm{s}}=\mathrm{8.5}\times {\mathrm{10}}^{\mathrm{11}}$ m s⁻¹ (with the same f_s in both panels). Regarding the dependence of 𝒲 on jet amplitude U_J, the PV gradient (viz. β_M) turns out to be a very good predictor, showing a smooth increase from small to large values. In particular, the PV gradient increases for values of U_J in the range 20 m s⁻¹ ≤ U_J ≤ 40 m s⁻¹ (unlike the ray tracing prediction), which is also visible in Fig. 4e and h. Regarding the dependence of 𝒲 on jet width σ_J, the PV gradient still predicts a smooth increase from small to large values, although the blue line is less linear than the red line. Note that, at jet latitude ϕ_J, the dependence of β_M on U_J and σ_J can be obtained by combining Eq. (23) with Eq. (6) and using Eq. (14). This yields

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{\mathit{\beta }}_{\mathrm{M}}\left(U_{\mathrm{J}},{\mathit{\sigma }}_{\mathrm{J}}\right)& \approx {\displaystyle \frac{\mathrm{2}\mathrm{\Omega }}{a}}{\cos }^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{J}}+{\displaystyle \frac{\cos {\mathit{\varphi }}_{\mathrm{J}}}{a^{\mathrm{2}}}}{\displaystyle \frac{U_{\mathrm{J}}}{{\mathit{\sigma }}_{\mathrm{J}}^{\mathrm{2}}}}\\ & =A+B{\displaystyle \frac{U_{\mathrm{J}}}{{\mathit{\sigma }}_{\mathrm{J}}^{\mathrm{2}}}}\end{array}\end{array}$$

with positive constants A and B, which explains the linear increase in β_M with increasing U_J in Fig. 8a and the smooth but nonlinear decrease in β_M with increasing σ_J in Fig. 8b.

Figure 8 Same as Fig. 6, except that here the solid blue line represents the non-dimensionalized β_M (see text for details).

Figure 9 Results from the parameter sweep involving 5×9 model simulations with five different values for σ_J and nine different values for U_J: (a) waveguidability 𝒲 and (b) scaled β_M (see text for details).

Overall, a comparison between Figs. 6 and 8 shows that, in our model experiments, waveguidability correlates much better with the meridional PV gradient than with the predictions from ray tracing theory. This result is synthesized in Fig. 9, which shows the results of a systematic parameter sweep varying σ_J between 5 and 15^∘ in steps of 2.5^∘ and varying U_J between 0 and 40 m s⁻¹ in steps of 5 m s⁻¹. Apparently, the variation in 𝒲 and β_M in this parameter space is qualitatively similar.

5 Analysis of the linearized equations

We now aim to obtain a deeper understanding for why WKB-based ray tracing theory provides partly misleading results regarding waveguidability. For this purpose we linearize the model equation (Eq. 16) about the background state and restrict attention to stationary flow, yielding

$$\begin{array}{}\text{(41)}& {\displaystyle \frac{u_{\mathrm{0}}}{a\cos \mathit{\varphi }}}{\displaystyle \frac{\partial q^{\prime }}{\partial \mathit{\lambda }}}+{\displaystyle \frac{v^{\prime }}{a}}{\displaystyle \frac{\mathrm{d}{q}_{\mathrm{0}}}{\mathrm{d}\mathit{\varphi }}}=-{\mathit{\lambda }}_{\mathrm{r}}{q}^{\prime }+F.\end{array}$$

Note that this equation does not imply any weak-variation assumptions, which means that the coefficients u₀ and q₀ may be arbitrary functions of latitude. This means that Eq. (41) can be applied to a broader class of problems than covered by WKB theory.

We, first, investigate whether or not our numerical solutions are close to linear. For this purpose we repeated the simulations with the forcing amplitude h₀ reduced by a factor of 10⁴ and the resulting perturbation fields increased by a factor of 10⁴. This so-called pseudo-linear solution is shown in Fig. 10b and compared with the full nonlinear solution shown in Fig. 10a. As it turns out, both solutions are very similar, especially regarding their global appearance. There are only small local differences in the northern midlatitudes, where the pseudo-linear solution is somewhat more peaked, suggesting that the effect of the nonlinear terms in the equation is similar to weak diffusion. Yet, the overall character of the solution is well captured by linear dynamics.

Figure 10 Numerical solution of the stationary part $\stackrel{\mathrm{‾}}{{\mathit{\psi }}^{*}}$ of the perturbation streamfunction for the weak narrow jet with u_J=10 m s⁻¹: (a) full nonlinear solution, and (b) pseudo-linear solution. The contours are chosen to be symmetric about zero and the contour interval is 0.4×10⁷ m² s⁻¹ (positive contours solid, negative contours dashed).

We now write the perturbation PV in terms of the perturbation streamfunction,

$$\begin{array}{}\text{(42)}& q^{\prime }={\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\psi }}^{\prime }.\end{array}$$

The zonal periodicity of our domain allows a Fourier series expansion of both F and ψ^′ like

$$\begin{array}{}\text{(43)}& F=\sum _{s=\mathrm{0}}^{\mathrm{\infty }}\widehat{F}(s,\mathit{\varphi })e^{is\mathit{\lambda }},{\mathit{\psi }}^{\prime }=\sum _{s=\mathrm{0}}^{\mathrm{\infty }}\widehat{\mathit{\psi }}(s,\mathit{\varphi })e^{is\mathit{\lambda }}.\end{array}$$

This yields the following equation for the (possibly complex) Fourier coefficients

$$\begin{array}{}\text{(44)}& \begin{array}{rl}{\displaystyle \frac{(\mathrm{1}-i\stackrel{\mathrm{̃}}{\mathit{ϵ}})}{\cos \mathit{\varphi }}}{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(\cos \mathit{\varphi }{\displaystyle \frac{\partial \widehat{\mathit{\psi }}}{\partial \mathit{\varphi }}}\right)& +\left[{\stackrel{\mathrm{̃}}{K}}_{\mathrm{s}}^{\mathrm{2}}-{\displaystyle \frac{(\mathrm{1}-i\stackrel{\mathrm{̃}}{\mathit{ϵ}})s^{\mathrm{2}}}{{\cos }^{\mathrm{2}}\mathit{\varphi }}}\right]\widehat{\mathit{\psi }}\\ & =-i{\displaystyle \frac{a^{\mathrm{3}}\cos \mathit{\varphi }}{s{u}_{\mathrm{0}}}}\widehat{F},\end{array}\end{array}$$

with the two dimensionless parameters

$$\begin{array}{}\text{(45)}& \stackrel{\mathrm{̃}}{\mathit{ϵ}}=a\cos \mathit{\varphi }{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{r}}}{s{u}_{\mathrm{0}}}}\end{array}$$

and

$$\begin{array}{}\text{(46)}& {\stackrel{\mathrm{̃}}{K}}_{\mathrm{s}}^{\mathrm{2}}=a^{\mathrm{2}}{\displaystyle \frac{\frac{\mathrm{1}}{a}\frac{\mathrm{d}{q}_{\mathrm{0}}}{\mathrm{d}\mathit{\varphi }}}{u_{\mathrm{0}}}}.\end{array}$$

Note that ${\stackrel{\mathrm{̃}}{K}}_{\mathrm{s}}^{\mathrm{2}}$ and ${\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}$ are related through ${\stackrel{\mathrm{̃}}{K}}_{\mathrm{s}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\varphi }={\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}$. Equation (44) is an inhomogeneous Helmholtz-like equation, and the appropriate boundary condition is $\widehat{\mathit{\psi }}=\mathrm{0}$ at both poles. The term in angle brackets, called T(ϕ), may be either positive or negative. For positive T, the character of the (homogeneous part of the) solution $\widehat{\mathit{\psi }}\left(\mathit{\varphi }\right)$ is oscillatory, while for negative T the character of the solution $\widehat{\mathit{\psi }}\left(\mathit{\varphi }\right)$ is exponentially evanescent. Our profiles of ${\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}\left(\mathit{\varphi }\right)$ shown in Fig. 4 indicate that the solution has an oscillatory character for solid-body rotation or for a weak jet, as long as the zonal wavenumber s is small enough. However, in the case of a strong jet there are small regions with evanescent behavior on both sides of the jet where ${\widehat{K}}^{\mathrm{2}}$ is negative. In either case, the solution $\widehat{\mathit{\psi }}$ at a given latitude not only depends on the local coefficients $\stackrel{\mathrm{̃}}{\mathit{ϵ}}$ and ${\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}$ at that latitude, but it also “feels” the values of these coefficients at other latitudes. In particular, the solution can “tunnel through” regions where T<0, if that region is small enough and if there is another oscillatory region beyond. This is in stark contrast with ray tracing theory, which is purely local in the sense that the direction of the ray path and the amplitude of the solution are given by the local properties of the basic state. The ray path solution only needs two turning latitudes in order to be completely confined in the meridional direction, resulting in a perfect zonal waveguide. By contrast, the mathematical character of Eq. (44) indicates that there may be situations in which ray tracing gives misleading results. In particular, the possibility of tunneling through finite regions of exponential behavior may explain why our numerical solutions show a weak wave signal in the Southern Hemisphere even for the strong narrow jet (Fig. 3b). Incidentally we note that a similar kind of non-localness is obtained in the framework of PV thinking due to the elliptic nature of the equation for PV inversion (Hoskins et al., 1985).

The other question we want to discuss here is why the PV gradient within the jet is possibly a more appropriate proxy for its waveguidability. Let us, for the moment, consider the unforced, undamped version of Eq. (44), which reads

$$\begin{array}{}\text{(47)}& {\displaystyle \frac{\mathrm{1}}{\cos \mathit{\varphi }}}{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(\cos \mathit{\varphi }{\displaystyle \frac{\partial \widehat{\mathit{\psi }}}{\partial \mathit{\varphi }}}\right)+\left({\stackrel{\mathrm{̃}}{K}}_{\mathrm{s}}^{\mathrm{2}}-{\displaystyle \frac{s^{\mathrm{2}}}{{\cos }^{\mathrm{2}}\mathit{\varphi }}}\right)\widehat{\mathit{\psi }}=\mathrm{0}.\end{array}$$

For weak damping, solutions of the full problem (Eq. 44) should locally be close to solutions of Eq. (47). Interestingly, all the information about the background atmosphere in the latter equation is contained in the dimensionless parameter ${\stackrel{\mathrm{̃}}{K}}_{\mathrm{s}}^{\mathrm{2}}={\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}/{\cos }^{\mathrm{2}}\mathit{\varphi }$. It implies that the solution of this equation can only depend on the combination of background fields contained in ${\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}$, and this should lead to a better correlation of waveguidability with ${\widehat{K}}_{\mathrm{s}}^{\mathrm{2}}$ than with β_M. However, this is in conflict with our numerical results from Sect. 4.2, which show a better correlation of waveguidability with β_M than with ${\widehat{K}}_{\mathrm{s}}$. The only possible conclusion is that our results about waveguidability and their dependence on the background flow are essentially determined by the forced-dissipative nature of our model configuration.

The remaining question in this line of argument is why the global character of the solution should depend on dissipation in an essential manner. At this point we suggest a heuristic argument according to which the global character of the solution is a result of a competition between the zonal speed of Rossby wave propagation (in the group velocity sense) and the dissipation. While the rate of dissipation is constant and fixed (being equal to λ_r), the speed of zonal propagation may depend, inter alia, on the meridional gradient of PV. The last idea would, indeed, be consistent with the concept of downstream development in the framework of PV thinking (see Fig. 8 in Wirth et al., 2018 for details): a given PV anomaly at the leading edge of an Rossby wave packet induces a new PV anomaly on the downstream side, and the speed of downstream propagation should be related to the strength of the meridional PV gradient. In competition against the omnipresent dissipation, those parts of the Rossby wave field that have a larger propagation speed can travel longer distances before they are dissipated. In a sense, regions with large PV gradients are like “express lanes” for Rossby wave packets. On the other hand, those parts of the wave that deviate from the “express lane” soon find themselves in regions with a much smaller meridional PV gradient, implying much slower zonal propagation. Since these wave packets are subject to the same damping rate, they can travel only much smaller distances before they are dissipated. It follows that 180^∘ downstream of our local forcing only those wave packets that used the “express lane” have survived significant damping, and this might explain why waveguidability as defined in Sect. 4.1 is large for jets with large β_M.

6 Summary, discussion, and conclusions

In this paper we investigated the waveguidability of idealized midlatitude jets using a barotropic model for flow on a sphere. Our focus was on stationary waves in a forced-dissipative model configuration with a local source for Rossby waves. This setup allowed us to quantify waveguidability (in percent) by diagnosing the latitudinal distribution of the explicitly simulated wave enstrophy in a longitudinal sector far downstream of the forcing. We carried out a systematic sensitivity study by varying the amplitude and the width of the jet.

Our numerical solutions indicate that waveguidability increases smoothly with increasing jet amplitude and with decreasing jet width. This result is in contrast to the prediction from WKB-based ray tracing theory, where the waveguide is either perfect or nonexistent, depending on the existence or nonexistence of two turning latitudes (one on either side of the jet). Another weakness of ray tracing theory is that the profile of the all-important stationary wavenumber saturates at some point for fairly narrow and strong jets; this theory would, therefore, not predict any further increase in waveguidability as the jet's amplitude is further increased, which is in conflict with the simulated behavior. To be sure, the overall poor predictive power of ray tracing theory in the current application is not surprising because the underlying WKB assumptions are not satisfied. We conclude that a literal application of ray tracing theory can be misleading; in particular, the existence of two turning latitudes is by no means a binary event at which the waveguidability suddenly changes in a fundamental way.

In addition, we found that waveguidability in our simulations correlates much better with the meridional gradient of background PV at jet latitude than with the predictions from ray tracing theory: both the meridional PV gradient and waveguidability vary smoothly as a function of jet amplitude and jet width. This result is consistent with the idealized concept of Rossby waves on a zonally oriented PV front suggesting that sharp meridional gradients of PV are conducive to zonal Rossby wave propagation. The good correlation with the background PV gradient was argued to be due to the forced-dissipative character of our model setup. A heuristic argument was presented based on the idea that the reach of a Rossby wave packet is determined by a balance between the dissipation rate and the speed of propagation of the Rossby wave packet; to the extent that the zonal speed of propagation is proportional to the meridional PV gradient, waveguidability in the zonal direction should, indeed, be related to the PV gradient.

Analysis of the linearized model equation (but without the WKB approximation) indicates that the amplitude of the forced stationary solution is governed by an inhomogeneous one-dimensional Helmholtz equation. This implies that the solution at a given grid point depends on the spatial distribution of the coefficients in a nonlocal fashion. Furthermore, a wave-like solution can tunnel through small regions of evanescent behavior, which means that the existence of two turning latitudes is not necessarily equivalent to a perfect waveguide. By contrast, the ray path solution only depends on the local values of the coefficients which prevents the possibility of tunneling. This argument helps to understand why the literal interpretation of ray tracing theory can be misleading and why the existence of two turning latitudes is often associated with values of waveguidability well below 100 %.

Overall, our analysis is not very affirmative regarding the utility of ray tracing for diagnosing waveguidability of Rossby waves in midlatitude jets. This seems to be in conflict with results from earlier publications such as Hoskins and Ambrizzi (1993), who extended the WKB concept to longitudinally varying background flows and showed that the spatial distribution of the stationary wavenumber has a certain amount of predictive power. Later, it was even argued that this kind of analysis may be useful in explaining the occurrence of persistent regimes and climate extremes (Hoskins and Woollings, 2015). The point is that these applications mostly rely on a broad qualitative interpretation of ray tracing. In fact, a comparison of the global distribution of the fields of ${\widehat{K}}_{\mathrm{s}}$ and β_M in Hoskins and Ambrizzi (1993, their Fig. 3) indicates that both fields suggest broadly the same behavior regarding the horizontal propagation of Rossby waves on the sphere. In other words, as long as the interpretation does not focus on the details of the diagnosed fields, it may well be qualitatively consistent with the observed behavior. What we have shown here is that the literal application of ray tracing and, in particular, the search for two turning latitudes as indication for a perfect waveguide must be considered with care.

Despite these caveats, an analysis that focuses on the existence of two turning latitudes in the background flow may still have some merit. As suggested by our Fig. 6, such a binary definition of events distinguishes background states with smaller values of 𝒲 from background states with larger values of 𝒲. In a statistical analysis which – from a large sample – distinguishes episodes with and without two turning latitudes, the sample with two turning latitudes will, on average, have larger waveguidability than the complementary sample without two turning latitudes. To the extent that the sample with two turning latitudes is associated with specific properties which the complementary sample is lacking, these properties are likely to have a systematic association with waveguidability.

On the more constructive side, our work suggests that the meridional PV gradient might be a suitable proxy for waveguidability. Obviously, this hypothesis is currently lacking any underlying closed theory. In addition, it should be tested in more realistic modeling frameworks, because our barotropic model entirely lacks the partitioning of PV into a contribution due to vorticity and a contribution due to static stability. By contrast, Ertel PV in a primitive equation model would be able to capture the sharp gradients of PV on isentropes at the boundary between the upper troposphere and the lower stratosphere (Martius et al., 2010). It is left for future work to explore suitable diagnostics for waveguidability in such more realistic frameworks.