Path properties of atmospheric transitions: illustration with a low-order sudden stratospheric warming model

The polar winter stratosphere typically supports a strong, cyclonic polar vortex, maintained by the thermal wind relation and meridional temperature gradient. A sudden stratospheric warming (SSW) event is a large excursion from this normal state, which can take many different forms. In split-type SSWs, the vortex splits completely in two. In displacement-type SSWs the vortex displaces far away from the pole (these can be considered wavenumber-2 and wavenumber-1 disturbances, respectively). Both types are “major” warmings, as the mean zonal wind reverses. In a “minor” warming, the zonal wind slows down significantly without completely reversing (Butler et al., 2015).

SSW is a rare event occurring about twice every three years, depending on the definition used (Butler et al., 2015). Its effects can propagate downward into the troposphere, altering the tropospheric jetstream and inducing extreme midlatitude surface weather events, including cold spells and precipitation (Baldwin and Dunkerton, 2001; Thompson et al., 2002). Abrupt cold spells severely stress infrastructures, economies and human lives, and every bit of extra prediction lead time is helpful for adaptation. Unfortunately, numerical weather prediction struggles to forecast SSW at any lead time longer than about two weeks (Tripathi et al., 2016). Understanding SSW is therefore important for practical forecasting as well as science, but this task remains difficult. Several different geophysical fields are often used as indices of SSW onset. One simple indicator is zonal-mean zonal wind at $60^{\circ}$N, which defines thresholds for minor and major warming (Charlton and Polvani, 2007; Butler et al., 2015). Another common indicator is the 10hPa geopotential height field, which was used by Inatsu et al. (2015) to estimate a fluctuation-dissipation relation in its leading empirical orthogonal functions (EOFs). Many studies have examined SSW precursors and dominant pathways through simulation and observation. Limpasuvan et al. (2004), for instance, catalogued the various wavenumber forcings, heat fluxes and zonal wind anomalies that accompanied each stage of SSW events from reanalysis data. While planetary wave forcing from the troposphere is an accepted proximal cause of SSW, the polar vortex’s susceptibility to such forcing, or “preconditioning”, is a nontrivial and debated function of its geometry (Albers and Birner, 2014; Bancalá et al., 2012). Tropospheric blocking is also thought to be linked to SSW; Martius et al. (2009) and Bao et al. (2017) found blocking to precede many major SSW events of the past half century. The diversity and complex life cycle of SSWs makes it difficult to build a unified picture of their onset.

In this article, we are interested in developing a detailed understanding of transition events between two states, at least one of which is typically long-lived. Consider, for example, a particle with position $x(t)$ moving in the double well potential energy landscape $V(x)=\frac{x^{4}}{4}-\frac{x^{2}}{2}$ (illustrated in Figure 1) and forced by stochastic white noise $\dot{W}$: $\dot{x}=-V^{\prime}(x)+\sigma\dot{W}$. If the system starts in the left well, it will tend to remain there a while, but occasionally the stochastic forcing will push it over the barrier into the right well. The natural predictor for this event is the committor: the probability of reaching the right well before the left well.

We denote this function by $q(x)$, which solves the Kolmogorov backward equation (to be introduced later). For this simple system the equation takes a form which can be solved exactly:

Note that the boundary conditions are implied by the probabilistic interpratation.

The committor is plotted in the right panel of Figure 1 for various noise levels. N.B., the potential landscape picture is not fully general, but is a useful mental model. The equations that determine $q(x)$ will be presented in section 2.

In the case of SSW, the long-lived states are the steady and disturbed circulation regimes of the stratospheric polar vortex. Recent work by Yasuda et al. (2017) has studied SSW in an equilibrium statistical mechanics framework, with these two stable states as saddle points of energy functionals. TPT takes a complementary non-equilibrium view, describing the long-time (steady-state) statistics of trajectories between the two states. For example, TPT introduces a probability density of reactive trajectories (or “reactive density”) indicating the regions where trajectories tend to spend their time en route from A to B. The system is said to be reactive at a point in time if it has most recently visited $A$ and will next visit $B$. The associated probability current of reactive trajectories (or “reactive current”) indicates the preferred direction and speed of transition paths. These detailed descriptors of the mechanism underlying a rare event can be expressed in terms of probabilistic forecasts like the committor $q(x)$, the probability of entering state A before reaching state B from a given initial condition $x$ (not in either A or B). The committor is the ideal probabilistic forecast in the usual variance-minimizing sense of conditional expectations (Durrett, 2013). Any other predictor of a transition derived through experiments and observations, such as vortex preconditioning and forcing at different wavenumbers (Albers and Birner, 2014; Bancalá et al., 2012; Martius et al., 2009; Bao et al., 2017) necessarily corresponds to an approximation of the committor.

Ensemble simulation is a commonly used method to estimate the committor at a single given initial condition by measuring the fraction of ensemble trajectories that achieve the rare event. This is challenging because many simulations are needed to generate enough rare events for significant statistical power. Recent and ongoing work aims to channel this computing power more efficiently in weather simulation using importance sampling and large deviation theory (Hoffman et al., 2006; Weare, 2009; Vanden-Eijnden and Weare, 2013; Ragone et al., 2018; Dematteis et al., 2018; Plotkin et al., 2019; Webber et al., 2019). The committor, and other quantities of interest, are fundamentally averages over sample trajectories. We can also express them as solutions to a concrete set of partial differential equations (PDEs) using basic stochastic calculus. TPT provides a framework to exploit these quantities and enhance our understanding of rare events, both from simulation data and from the fundamental equations of motion.

TPT has been applied primarily to molecular dynamics simulation to determine reaction rates and pathways of complex conformational transitions (Vanden-Eijnden and E, 2010; Metzner et al., 2006; E and Vanden-Eijnden, 2006); however, the framework does not depend on the details of the underlying dynamical system. TPT deals particularly with stochastically forced systems, such as Brownian dynamics of particles. Stochastic forcing applies quite generally; while the climate system is deterministic in principle, nonlinear interactions between resolved and unresolved scales inevitably leads to resolution-dependent model errors that can be approximated as stochastic. Hasselmann (1976) originally formulated stochastic climate models to capture the influence of quickly evolving “weather” variables on the slowly evolving “climate” variables. Stochastic parameterization remains an active area of research. For example, Franzke and Majda (2006) had success in capturing energy fluxes of a 3-layer quasigeostrophic model by projecting onto ten EOF modes and treating the remainder as stochastic forcing. Kitsios and Frederiksen (2019) addressed the challenge of designing consistent numerical schemes for subgrid-scale parametrization. Deep convection in the atmosphere and turbulence in the ocean boundary layer are two examples of multiscale processes that are especially challenging to resolve.

The aim of this article is to introduce the key quantities and relations describing the path properties of rare atmospheric events. Computing those quantities for more complicated systems than the low order model studied here is a significant and, we argue, worthwhile challenge. We do not address that computational challenge here. Instead we note that development of approximation techniques for TPT and related quantities in high dimensional settings is an active area of research Thiede et al. (2019); Bowman et al. (2009); Chodera and Noe (2014). These methods incorporate both simulated and observational data and we outline them briefly in the conclusion.

The paper is organized as follows. Section 2 describes the dynamical model we use, building on work by Ruzmaikin et al. (2003) and Birner and Williams (2008). Section 3 describes the mathematical framework of TPT, with detailed, but informal, derivations mainly put in the supplement. Section 4 explains the methodology and results particular to this model.

Holton and Mass (1976) studied “stratospheric vacillation cycles,” a certain kind of minor warming in which zonal wind oscillates on a roughly seasonal timescale. They posited a mechanism of wave-mean flow interaction, which continues to be an important modeling paradigm. The quasi-geostrophic equations are confined to a $\beta$-plane channel from $60^{\circ}$ to the north pole, and the streamfunction is perturbed from below by orographically induced planetary waves, specified through the lower boundary condition. The Holton-Mass model combines the zonal-mean flow equations,

and the linearized quasigeostrophic potential vorticity equation,

These equations represent (2) zonal and (3) meridional momentum balance, (4) conservation of energy, (5) conservation of mass, and (6) a combination of all these for the perturbations. Overbars and primes represent zonal averages and perturbations. $\Phi$ is the geopotential height; $H=7\,km$ is the atmospheric scale height; $z=-H\ln(\frac{p}{p_{0}})$ is log-pressure; $\rho_{s}=\rho_{0}e^{-z/H}$ is a standard density profile; $\alpha=\alpha(z)$ is an altitude-dependent damping coefficient; and $\overline{T}^{*}=\overline{T}^{*}(y,z)$ is the radiative equilibrium temperature. Holton and Mass (1976) studied the solution consisting of single Fourier modes in the zonal and meridional directions:

where $\psi^{\prime}$ is the zonal perturbation of $\psi=(g/f_{0})\Phi$. $k=2/(a\cos 60^{\circ})$ and $\ell=3/a$ are the zonal and meridional wavenenumbers, where $a$ is Earth’s radius. These wavenumbers are commonly observed in real SSWs and used in theoretical studies (Birner and Williams, 2008; Ruzmaikin et al., 2003; Yoden, 1987; Holton and Mass, 1976); a split-type SSW is an extreme wavenumber-two perturbation. The lower boundary condition at $z=z_{B}$ (the tropopause) is

where $h$ is a topographically induced perturbation to geopotential height at the tropopause. Holton and Mass found that for a certain range of $h$, this system has qualitatively different regimes: a steady eastward zonal flow close to radiative equilibrium, and a weaker zonal flow with quasi-periodic “vacillations” from eastward to westward, even under constant forcing. Each vacillation cycle consists of a sudden warming and cooling over the timescale of weeks. Although these individual cycles are interesting weather events unto themselves, in this paper we think of the vacillations as occurring within a general climate regime that is conducive to sudden warming, as opposed to the steady flow state, which is not. Transitions between these two regimes, which we focus on here, are more accurately described as climatological shifts than weather events. The study by Ruzmaikin et al. (2003) varies $h$ on an interannual timescale, with each single winter season occupying one of the two stable states and generating its daily weather accordingly. Hence, for this paper we will use the term “climate transitions.”

The original Holton-Mass model discretizes the above PDE with finite differences across 27 vertical levels, which is assumed to be close to a continuum limit. Following several studies at this resolution (Holton and Mass, 1976; Yoden, 1987; Christiansen, 2000), Ruzmaikin et al. (2003) did the most severe truncation possible, resolving only three vertical levels (including fixed boundaries) for easy analysis and exploration of parameter space. This reduces phase space to only three degrees of freedom: $U(t)$, which modulates $\overline{u}$ as a sine jet; $X(t)=\text{Re}\{\Psi(t)\}$; and $Y(t)=\text{Im}\{\Psi(t)\}$. $X$ and $Y$ modulate the amplitude and phase of the perturbation streamfunction:

Carrying the Ansatz through the QG equations, Ruzmaikin et al. (2003) derived the following system:

The primary control parameter, $h$, represents topographic forcing and other sources of planetary waves, such as land-sea ice contrast. While Ruzmaikin et al. (2003) also varies $\Lambda$, representing vertical wind shear, we will only vary $h$ and set $\Lambda$ constant. Time derivatives $\dot{h}$ and $\dot{\Lambda}$ are zeroed, removing transient forcings such as seasonality effects. Appendix A and Ruzmaikin et al. (2003) have a detailed list of parameters.

Remarkably, this hugely simplified model retains the qualitative structure of the Holton-Mass model as a bistable system for a certain range of $h$ between the critical values $h_{1}\approx 20\,m$ and $h_{2}\approx 160\,m$, as shown in the bifurcation diagram of Figure 2. Blue points represent the normal state of the vortex, in approximate thermal wind balance with the radiative equilibrium temperature field (henceforth called the “radiative solution”). Red points represent a disturbed vortex, with weaker zonal wind and vacillations. This climatological regime supports more SSW events, and is henceforth called the “vacillating solution.” We use the same blue-red color scheme consistently here to represent these two states. Transitions between them happen on interannual time scales, affecting each year’s likelihood of SSW events. The structure of transitions is illustrated in Figure 3: as $h$ increases slowly past the bifurcation threshold $h_{2}$, the system enters a series of rapid, large-amplitude oscillations that spiral into the weaker-circulation state.

In Figures 2 and 3, transitions require crossing the bifurcation threshold $h_{2}$, where the radiative solution ceases to exist. Birner and Williams (2008) introduced additive white-noise forcing in the $U$ variable to model unresolved gravity waves and found that these perturbations were sufficient to excite the system out of its normal state and into a vacillating regime. In Figure 4 we illustrate stochastic trajectories of the system for three different (fixed) values of $h$. (For numerical reasons we also add a small amount of independent white noise to $X$ and $Y$ variables). Even when $h$ is far below $h_{2}$, transitions still occur, and in fact the preference for the vacillating solution branch increases quickly with $h$.

Birner and Williams (2008) used direct numerical simulation and the Fokker-Planck equation to calculate long-term occupation statistics, i.e., how much time on average was spent in each regime and the mean first passage time before a transition to the vacillating regime, all for a range of forcing and noise levels. Our approach differs in both target and methodology. We aim to characterize the transition process between the two states, to monitor its progress in real time, as well as to describe statistics of the transition over many realizations. Methodologically, transition path theory phrases these questions in terms of the generator of the stochastic process, a differential operator that encodes all information about the behavior of the process.

TPT characterizes the steady state statistics of transitions between states. In this section, we introduce key quantities needed to introduce TPT as applied to the Ruzmaikin model to obtain a more complete picture than we get from the sample paths shown above. For mathematical details, see the supplement and background literature (E and Vanden-Eijnden, 2006; Metzner et al., 2006; Vanden-Eijnden and E, 2010).

We begin this section by describing the kinematic path characteristics of the process in its three-dimensional phase space, according to the quantities described above. Following this purely geometrical description, we will suggest some dynamical interpretations and compare with previous studies. Finally, we will map statistical features as functions of background parameters.

The Ruzmaikin model is attractive for demonstrating use of the tools introduced in Section 3 due to its low-dimensional state space, in which PDEs can be solved numerically using standard methods such as our finite volume scheme. We tested the committor’s accuracy empirically by randomly selecting 50 cells in our grid (this is $0.04\%$ of the grid) and evolving $n=60$ stochastic trajectories forward in time from each, stopping when they reach either set $A$ or set $B$. The fraction of trajectories starting from $z$ that first reach $B$ is taken as the empirical committor at point $z$ and is denoted $\widetilde{q}^{+}$, which we compare with the predicted committor $q^{+}$ from the finite volume scheme. Figure 5 clearly demonstrates the usefulness of the committor for probabilistic forecasting. The left column displays the committor calculated from finite volumes, averaged in the $Y$ direction, for two different forcing levels $h$. The right column shows a scatterplot of $\widetilde{q}^{+}$ vs. $q^{+}$ at the 50 randomly selected grid cells. We expect the points to fall along the line $\widetilde{q}^{+}=q^{+}$ with some spread proportional to the standard deviation of the binomial distribution, $\sqrt{q^{+}(1-q^{+})/n}$. Approximating this as a Gaussian, we have plotted red curves enclosing the 95% confidence interval, which indeed contains approximately 95% of the data points. Statistical sampling error can explain the observed level of deviation, although discretization error (from solving the PDE and from time-stepping) may also contribute.

Figure 5 also shows how $q^{+}$ responds to increasing $h$, even far below the bifurcation threshold: the committor values throughout state space become rapidly skewed toward unity (meaning more red in the picture). This means that even slight perturbations can kick the system out of state $A$ toward state $B$. Another indicator is the “isocommittor surface”, the set of points $z$ such that $q^{+}(z)=1/2$; that is, the system has equal probability of next entering set $A$ or $B$. In the left-hand column this is the set of gray points (averaging out the variable Y). For low forcing, this surface tightly encloses set $B$, meaning the system must wander very close before a transition is imminent. For high forcing values, the isocommittor hugs set $A$ more closely, meaning that small perturbations from this normal state can easily push the system into dangerous territory. In Figure 6, the isocommittor is shown as a set of gray points in a 3D plot viewed from various vantage points. In the low-$U$ region, the isocommittor resembles a spiral staircase, reflecting the spiral-shaped stable manifold of the fixed point in set $B$. Different initial positions with the same streamfunction phase, differing only slightly in the $U$ direction, can have drastically different final destinations. These spiral surfaces are responsible for the blue lobes in the lower part of Figure 5, but they disappear at higher noise.

Figures 7 and 8 display numerical solutions of the equilibrium density $\pi$ and reactive density $\rho_{R}\propto\pi q^{-}q^{+}$ for two forcing levels. While $\pi$ indicates where $Z_{t}$ tends to reside, $\rho_{R}$ indicates where $Z_{t}$ resides given a transition from $A$ to $B$ is underway. As $h$ increases, even far below the bifurcation threshold, $\pi$ responds strongly, shifting weight toward state $B$. On the other hand, the reactive density displays similar characteristics for all $h$ values. In the $XU$ plane, the two lobes of high reactive density surrounding $A$ indicate that zonal wind tends to remain strong for a while before dipping into the weaker regime. Viewing the same field in the $XY$ plane (Figure 8) reveals a halo of intermediate density about set $A$. While many different motions would be consistent with this pattern, the coming figures verify that the early stages of transition have circular loops in the $XY$ plane, meaning zonal movement of the streamfunction’s peaks and troughs. The exact streamfunction phase corresponding to the $(X,Y)$ position is calculated as follows. Recall the streamfunction is $\psi^{\prime}=\text{Re}\{\Psi(t)e^{ikx}\cos(\ell y)\}$, where $\Psi=X+iY$. In polar coordinates, $\Psi=\sqrt{X^{2}+Y^{2}}e^{i\phi}$, where $\phi=\tan^{-1}(Y/X)$. The full streamfunction is

where $\lambda$ is longitude.

The angle from the origin in the $XY$ plane indicates the zonal streamfunction phase, and circular motion indicates zonal movement. (This “looping” motion is indeed shared by the transition path samples shown in Figures 9 and 10, to be described later.)

The darkest (most-trafficked) region of this loop is the sector $\frac{\pi}{4}\lesssim\phi\lesssim\frac{\pi}{2}$. The relationship between $(X,Y)$ and $\lambda$ indicates $\psi^{\prime}$ is maximized at longitudes $\lambda=\{-\frac{\phi}{2},\pi-\frac{\phi}{2}\}$. As the maximum reactive density occurs around $\phi=\frac{3\pi}{8}$, the streamfunction peaks are at $\{-\frac{3\pi}{16},\frac{13\pi}{16}\}\approx\{326^{\circ},146^{\circ}\}$. What is the significance of this phase relative to the lower boundary forcing? Recalling the forcing form $\text{Re}\{\Psi(z_{B},t)e^{ikx}\}=h\,\text{Re}\{e^{i2\lambda}\}\propto\cos(2\lambda)$, the bottom peaks are located at $\lambda=\{0,\pi\}$. Hence, the bulk of the transition process happens when the perturbation streamfunction at the mid-stratosphere lags the lower boundary condition by $\frac{3}{16}\pm\frac{1}{16}$ of a wavelength. Meanwhile, the $XU$ plane reveals what happens to the zonal wind speed during the SSW transition. The high-reactive density region discussed above coincides with the crescent-shaped bridge of high density between the sets in Figure 7. This suggests that in an SSW, the zonal wind weakens while the streamfunction stays in that particular phase window.

The pictures of reactive density suggest that reactive trajectories tend to loop around set $A$, physically meaning the streamfunction tends to travel in one direction before slowing down, but they technically convey no directional information to explicitly support this claim. For this, we turn to the reactive current. We computed the discrete-space effective current matrix $f_{ij}^{+}$, directly from the finite volume discretization and numerical solutions of $\pi$, $q^{+}$ and $q^{-}$. Physically, this matrix represents the flux of a vector field from grid cell $i$ to cell $j$. From this we calculated the maximum-current paths as described in Metzner et al. (2009) and displayed the results in Figure 9 for a forcing level of $h=30\,m$ (other levels are qualitatively similar). Both the $XU$ and $XY$ views are shown. Superimposed on these paths are seven actual reactive trajectories that occurred during a long stochastic simulation, to demonstrate features that are captured by the dominant pathways. The dominant path from $A$ to $B$ indeed contains a half-loop in the $XY$ plane in the clockwise direction, which means an eastward phase velocity. With a smaller set $A$, this dominant path would contain more of these loops. However, during the next transition stage, the streamfunction slows to a halt at the phase angle $\phi=\frac{\pi}{2}$, doubles back and travels westward as zonal wind loses strength. The smear of high density in the neighborhood $\phi\sim\frac{3\pi}{16}$ therefore comprises not only a precipitous drop in zonal wind (which happens at the edge of that region) but also a backtrack, this time with weaker background zonal wind. This behavior is borne out by the trajectory samples, which vacillate in the upper-middle section of the $XY$ plot. These paths are displayed as space-time diagrams of the streamfunction in Figure 10. In panel (a), the dominant path’s two loops correspond to two troughs moving east past a fixed longitude before the slowdown. The random streamfunction trajectories shown in (b)-(g) do not follow this representative history exactly, but they do combine elements of it: steady eastward wave propagation followed by slowdown and reversal. Each stage can have multiple false starts. Notably, the slowdown consistently happens at the same phase, with peaks at $\sim 120^{\circ}E$ and $300^{\circ}E$, at roughly the same phase as found from the density plots. In fact, the figures show a brief slowdown every time the streamfunction passes this phase. This can be thought of as a representation of blocking events that often accompany sudden stratospheric warmings. The third transition path shown is an exception to the general pattern, making a final turn toward the east instead of to the west. This outlier of a reactive trajectory can also be seen in Figure 9, as the single green trajecory that decreases in $X$ before decreasing in $U$ instead of the other way around.

This kinematic sequence of events has a dynamical interpretation with precedent in prior literature. A critical ingredient of SSW is meridional eddy heat flux, which in this model takes the form $\overline{v^{\prime}\Phi_{z}^{\prime}}\propto hY$ (see the supplement for a detailed derivation.) This term shows up explicitly as a negative forcing on $U$ in equation (15), showing that a reduced equator-to-pole temperature gradient in turn weakens the vortex via the thermal wind relation. The association of heat flux with SSW has been demonstrated in reanalysis (Sjoberg and Birner, 2012) and in detailed numerical simulations of internal stratospheric dynamics, even with time-independent lower boundary forcing (Scott and Polvani, 2006). This relationship favors the phase $\phi=\frac{\pi}{2}$ as the most susceptible state for SSW onset, which is exactly picked out by the dominant transition pathway in Figure 9.

However, immediately after the wind starts weakening at $\phi=\frac{\pi}{2}$, where the streamfunction lines up with its lower boundary condition, the phase velocity reverses, giving rise to the westward lag of $\frac{3\pi}{16}$ we observed in the reactive density. A similar phase lag has also been observed in more detailed numerical studies. For example, Scott and Polvani (2006) observed a lag of $\frac{\pi}{2}$ across the whole stratosphere ($\frac{\pi}{4}$ at the mid-level), quite similar to our result. They found that vortex breakup was preceded by a long, slow build-up phase in which the vortex became increasingly vertically coherent, only to be ripped apart by an upward- and west-propagating wave. In an experiment with slowly-increasing lower boundary forcing, Dunkerton et al. (1981) saw a phase lag across the whole stratosphere that increased from $\sim 100^{\circ}$ to $\sim 180^{\circ}$ ($50^{\circ}$ to $90^{\circ}$ between the lower boundary and the mid-stratosphere) over the course of the warming event. They attribute this phase tilt to the zonal wind rapidly reversing and carrying the streamfunction along. The weakening zonal wind simply removes the Doppler shift from the Rossby wave dispersion relation, $\omega=Uk-\beta k/(k^{2}+\ell^{2})$, allowing the waves to revert to their preferred westward phase velocity. This balance is also clear in equations (13-14): ignoring the damping and forcing terms, the dynamics read $[\dot{X},\dot{Y}]=[(sU-r)Y,-(sU-r)X]$. For weak $U$, rotation in the $XY$ plane is anti-clockwise and phase speed is westward.

Let us re-emphasize the probabilistic interpretation of reactive density. We have found that transitions from $A$ to $B$ are accompanied by anomalous increases in meridional heat flux. In other words, $\mathbb{P}\{\text{High heat flux$|$SSW}\}$ is high. As noted earlier, this does not imply that $\mathbb{P}\{\text{SSW$|$high heat flux}\}$ is also high. The identified streamfunction phase is not a “danger zone” in the sense that a trajectory entering this region is at higher risk of falling into state $B$; the committor alone conveys that information. Rather, a trajectory is highly likely to pass through that region given that it is reactive. Notably, reactive trajectories are unlikely to take a straight-line path from $A$ to $B$ with $U$, $X$ and $Y$ changing linearly. This unrealistic path would represent a zonally stationary streamfunction growing steadily in magnitude, while zonal wind falls off gradually. At higher noise levels, however, the system would be increasingly dominated by pure Brownian motion, and such a path would become more plausible.

We now turn to a quantitative comparison of committors and transition rates for different forcing and noise levels. These trends illustrate the effects of modeling choices and global change on the climatology of SSW. Planetary wave forcing, $h$, varies across days and seasons as well as different planets. The strength of additive noise, $\sigma$, is a modeling choice intended to represent gravity wave drag. Different stochastic parametrizations will vary in their effective $\sigma$ value, and it is important to understand the sensitivity of SSW to model choices (Sigmond and Scinocca, 2010). Furthermore, long-term climate change may cause both parameters to drift, altering the occurrence of SSW-induced severe weather events.

The measure of the relative “imminence” of a vacillating solution vs. a radiative solution, as described in the background section on committors, is the equilibrium density-weighted average committor, denoted $\langle q^{+}\rangle_{\pi}$. Figure 11(a) shows this quantity for $25\leq h\leq 45$ ($m$) and $0.4\leq\sigma_{3}\leq 2.0$ ($m/s/day^{1/2})$. Two trends are clearly expected from the basic physics of the model. First, as seen in Figures 7 and 8, $\langle q^{+}\rangle_{\pi}$ should increase with $h$. Second, in the limit of large noise and infinite domain size, the dominance of Brownian motion will smooth out the committor function and make $\langle q^{+}\rangle_{\pi}$ tend to an intermediate value between zero and one. On the other hand, as noise approaches zero, the dynamical system becomes increasingly deterministic, and the ultimate destination of a trajectory will depend entirely on which basin of attraction it starts in. The boundary, or separatrix, between these two basins is the stable manifold of the third (unstable) fixed point. In the case of a potential system, of the form $\dot{x}=-\nabla V(x)$, this boundary would be a literal ridge of the function $V$. Our system admits no such potential function, but this is a useful visual analogy. The committor function becomes a step function in the deterministic limit, with the discontinuity located exactly on this boundary. The addition of low noise moves the committor-$\frac{1}{2}$ surface away from the separatrix, possibly asymmetrically: one basin will shrink, becoming more precarious with respect to random perturbations, while the other will expand, becoming a stronger global attractor. Which basin will shrink is not evident a priori, so we compute the averaged committor, $\langle q^{+}\rangle_{\pi}$, as a summary statistic which will increase when the basin of $B$ expands.

Figure 11(a) plots the trends in $\langle q^{+}\rangle_{\pi}$ as a function of $h$ (along the horizontal axis) and $\sigma_{3}$ (along the vertical axis). The two basic hypotheses are verified: $\langle q^{+}\rangle_{\pi}$ increases monotonically as $h$ increases, and $\langle q^{+}\rangle_{\pi}\sim\frac{1}{2}$ as $\sigma$ increases, no matter the value of $h$. The less predictable behavior is in the range $h=35\,m$, $0.75\leq\sigma_{3}\leq 1.0$, where $\langle q^{+}\rangle_{\pi}$ displays non-monotonicity with respect to noise, at low noise levels. As $\sigma_{3}$ increases, the average committor increases from $\sim 0.4$ to $\sim 0.6$, and then decreases again. The four committor plots at the bottom of 11 illustrate the trend graphically. At low noise, the $A$ basin includes winding passageways leading from the small-$U$ region back to $A$. Small additive noise closes them off, effectively expanding the $B$ basin. As noise increases and Brownian motion dominates the dynamics, committor values everywhere relax back to less-extreme values, reflecting the unbiased nature of Brownian motion.

Despite the coarse grid resolution, the first-order effect of noise is clear. At the low and high margins of $h$, where falling into state $A$ and $B$ respectively is virtually certain, an increase in noise decreases this virtual certainty, and the trend continues at larger noise to attenuate $\langle q^{+}\rangle_{\pi}$ to its limit of 1/2. The middle $h$ range, however, behaves differently. Whereas $h=35\,m$ appears to balance out the basin sizes at low noise, a slight noise increase tends to kick the system out of the $A$ basin and toward $B$, more so than the other way around. At higher noise, the committor relaxes back to 1/2. Examining the committor fields, it’s clear that the isocommittor surface location does not move back and forth; rather, it moves toward $A$, and then the rest of the field flattens out.

Figure 12 shows trends in the return times of SSW with varying $h$ and $\sigma_{3}$. There are several different return times of interest. The first, shown in panel (a), is the total expected time between one transition event and the next, whose reciprocal is the rate $R_{AB}$, the number of transitions per unit time. The return time is a symmetric quantity between $A$ and $B$, since every forward transition is accompanied by a backward one. Among the parameter combinations, $(h,\sigma_{3})\approx(35\,m,0.75\,m/s/day^{1/2})$ is the one which minimizes return time, or equivalently maximizes the transitions per unit time. $h=35\,m$ is a forcing level which approximately balances out the time spent between the two sets, making transitions relatively common. At lower noise, transitions are exceedingly rare, and at higher noise the two states cease to be long-lived. However, this symmetric quantity does not capture information about the relative speed of transition from $A$ to $B$ vs. from $B$ back to $A$. Panel (b) shows a different passage time, which is the average time between the end of a backward $(B\to A)$ transition and the end of the next forward ($A\to B$) transition, which we call $T_{AB}$. This is computed as the reciprocal of the rate constant $k_{AB}$, as described in the previous section. In other words, the stopwatch begins when the system returns to $A$ after having last visited $B$, and ends when the system next hits $B$. This metric is asymmetric: a smaller $A\to B$ return time indicates that the forward transition is faster than the backward transition. Panel (b) shows the complementary $B\to A$ return time, $T_{BA}$. Unsurprisingly, an increase in $h$ causes a decrease in $T_{AB}$ and an increase in $T_{BA}$ regardless of the noise level. The noise level has a less obvious effect. Whereas $T_{BA}$ decreases monotonically with increasing noise, regardless of $h$, the forward time $T_{AB}$ is minimized by a mid-range noise level of $\sigma_{3}\approx 0.75\,m/s/day^{1/2}$. This is another reflection of the bias toward state $B$ that is effected by adding noise to a very low baseline.

Transition path theory is a framework for describing rare transitions between states. We have described TPT along with a number its key ingredients like the forward and backward committor functions. While TPT has been applied primarily to molecular systems, we believe it offers valuable insight into climate and weather phenomena such as sudden stratospheric warming, primarily through committors, reactive densities, and reactive currents. Of interest apart from its role in TPT, the forward committor defines an optimal probabilistic forecast, borne out by direct numerical simulation experiments. The reactive densities and currents describe the geometric properties of dominant transition mechanisms at low noise. In applying TPT to a noisy, truncated Holton-Mass model, we find that transitions tend to begin with a drop in mean zonal wind and a reversal of the streamfunction’s phase velocity at a particular streamfunction phase. This is consistent with the significance of blocking precursors to SSW as found in Martius et al. (2009), insofar as this idealized model can represent them. We also find that noise has a non-monotonic effect on the overall preference for a vacillating state, measured by the average committor. At a forcing of $h\approx 35\,m$, where the isocommittor surface (essentially the basin boundary of the deterministic dynamics) divides the space approximately in half, we find that raising the noise tilts the balance decisively toward the vacillating solution. Still larger noise evens the whole field out. The transition rate constant shows a similar dependence on $h$ and $\sigma_{3}$.

In future work, we plan to scale these methods to more realistic and complex models as well as observational data, where predicting SSW remains an active area of research. In this high-dimensional setting, the generator may be unknown or computationally intractable. Any state space with more than $\sim 5$ degrees of freedom is beyond the reach of a finite-volume discretization, because the number of grid cells increases exponentially with dimension. However, there is a generic insight that physical models evolve on very low-dimensional manifolds within the full available state space. A growing body of research in molecular dynamics (Thiede et al., 2019), fluid dynamics (Giannakis et al., 2018; Froyland and Junge, 2018), climate dynamics (Giannakis and Majda, 2012; Sabeerali et al., 2017), and general multiscale systems (Harlim and Yang, 2018; Berry et al., 2015; Giannakis, 2015, 2019) exploits the intrinsic low-dimensionality to represent the infinitesimal generator more efficiently. While here we represented the generator as a finite-volume or finite-difference operator on a grid, one can also write it in a basis of globally coherent functions, such as Fourier modes, or more generally harmonic functions on a manifold. Given only data, without an explicit form of the dynamics, this manifold and the basis functions can be estimated from (for example) the diffusion maps algorithm, and the generator’s action on this basis can be approximated from short trajectories. These ideas can be applied to computing the dynamical statistics that have been the focus of this paper (Thiede et al., 2019). We hope that these techniques will enable more efficient observation strategies for targeted data assimilation procedures with the goal of tracking the progression of specific extreme events, including hurricanes and heat waves as well as sudden stratospheric warming.

Data from this study will be made available upon request.