Emergence of representative signals for sudden stratospheric warmings beyond current predictable lead times

Major sudden stratospheric warmings (SSWs) are extreme wintertime circulation events of the Arctic stratosphere that are accompanied by a breakdown of the polar vortex and are considered an important source of predictability of tropospheric weather on subseasonal to seasonal time scales over the Northern Hemisphere midand highlatitudes. However, SSWs themselves are difficult to predict, with a predictability limit of around one to two weeks. The predictability limit for determining the type of event, i.e., wave-1 or wave-2 events, is even shorter. Here we analyze the dynamics of the vortex break5 down and look for early signs of the vortex deceleration process at lead times beyond the current predictability limit of SSWs. To this end, we employ a mode decomposition analysis to study the potential vorticity (PV) equation on the 850 K isentropic surface by decomposing each term in the PV equation using the empirical orthogonal functions of the PV. The first principal component (PC) is an indicator of the strength of the polar vortex and starts to increase from around 25 days before the onset of SSWs, indicating a deceleration of the polar vortex. A budget analysis based on the mode decomposition is then used to 10 characterize the contribution of the linear and nonlinear PV advection terms to the rate of change (tendency) of the first PC. The linear PV advection term is the main contributor to the PC tendency at 25 to 15 days before the onset of SSW events for both wave-1 and wave-2 events. The nonlinear PV advection term becomes important between 15 to 1 days before the onset of wave-2 events, while the linear PV advection term continues to be the main contributor for wave-1 events. By linking the PV advection to the PV flux, we find that the linear PV flux is important for both types of SSWs from 25 to 15 days prior 15 to the events but with different wave-2 spatial patterns, while the nonlinear PV flux displays a wave-3 wave pattern, which finally leads to a split of the polar vortex. Early signs of SSW events arise before the one to two week prediction limit currently observed in state-of-the-art prediction systems, while signs for the type of event arise at least one week before the event onset.

the breakdown of the polar vortex should be distinct between displacement (wave-1) and split (wave-2) events and should also be distinguishable from normal winter days (without SSWs). Therefore, understanding the dynamics of the vortex disruption and identifying signals that contribute to the vortex deceleration are crucial for improving the predictability of SSWs and of each type of event, and ultimately, of the weather at the Earth's surface.
Since the stratospheric circulation is well described by Ertel's potential vorticity (PV) (McIntyre, 1982), the evolution of the 60 polar vortex during SSWs can also be captured by the changes in the values and structure of the PV in the stratosphere. As discussed above, while the polar vortex undergoes breakdown in each major SSW event, the associated vortex structures are different for the two types of SSW events. Decomposing the PV into an empirical orthogonal function (EOF) basis, we can identify the PV structure that best describes the weakening of the polar vortex and subsequently investigate how its corresponding principal component (PC) changes with time. By projecting the other variables from the PV equation (i.e., the zonal and meridional wind) onto the EOF basis from PV, one can analyze the contribution of each term of the equation to the changes in the PC time series in order to identify the dynamical processes that are the most relevant for the weakening of the polar vortex. This approach was proposed by Aikawa et al. (2019) and called mode decomposition analysis. Aikawa et al. (2019) applied the mode decomposition analysis to diagnose the atmospheric blocking development in the Eastern Pacific and Central Atlantic, and demonstrated that the blocking index can be faithfully reconstructed using only the first 10 EOF modes. The vorticity 70 equation was then decomposed into three terms (i.e., linear advection, nonlinear mode-to-mode interaction, and dissipation), and their contribution to the combined time evolution of the first 10 PC time series was subsequently investigated. Their results showed that the nonlinear interaction terms contribute to the increase in the amplitude of the blocking index in both regions (Eastern Pacific and Central Atlantic), but that their contributions are different. Since each term in the vorticity equation can be linearly reconstructed using the EOF modes (that correspond to specific spatial patterns) and the PCs, this method allowed the 75 identification of the wind and vorticity patterns that are crucial for the development of the blocking. As the results from Aikawa et al. (2019) indicate the effectiveness of mode decomposition analysis in studying dynamically-driven events, we use the same method to study the development of SSWs, which are also driven by wave dynamics (e.g., Matsuno, 1971). Indeed, we find that the vortex weakening can be represented by the evolution of the EOF modes extracted from PV. We therefore employ a budget analysis of the PV equation in the stratosphere to quantify the contribution of each EOF mode to the dynamical processes that 80 lead to the deceleration of the polar vortex and the subsequent onset of SSWs.
The onset of SSW events is associated not only with the anomalously large excitation of wave activity in the troposphere (Matsuno, 1971;Polvani and Waugh, 2004;Lindgren et al., 2018), but also with the stratospheric mean state and stratospheric wave anomalies prior to SSWs (Hitchcock and Haynes, 2016;Jucker, 2016;Birner and Albers, 2017;de la Cámara et al., 2019). Moreover, split and displacement SSW events exhibit distinct pre-warming evolutions (Charlton and Polvani, 2007;85 Matthewman et al., 2009;Bancalá et al., 2012;Albers and Birner, 2014). For example, the zonal wavenumber of the wave flux leading to the breakdown of the polar vortex can be different for the two types of SSW events (e.g., Bancalá et al., 2012).
Some studies suggest that the explosive growth of wave amplitude is triggered by resonant behavior, which is also different between the two types of SSW events (e.g., Esler and Matthewman, 2011;Matthewman and Esler, 2011). Albers and Birner vortex preconditioning that is conducive to developing the respective split and displacement SSWs events. Given the distinct dynamical developments of the two types of SSWs, one should be able to observe different evolutions of the PV terms by the mode decomposition analysis for each type of SSW event. Since displacement (split) events are mainly related to wave-1 (wave-2) wave activity, in this study, SSWs are classified into wave-1 and wave-2 events based on the dominant wavenumber that leads to the breakdown of the vortex. As some studies point out that the dynamical process of the vortex breakdown starts 95 earlier than the current predictable lead time of two weeks (Polvani and Waugh, 2004;Jucker and Reichler, 2018), one would expect to see signals indicative of the vortex breakdown appearing in the mode equation budget before the onset of the SSWs.
Our goal in the current study is therefore to identify signals that are representative of SSWs, i.e., distinct from normal winter days, ahead of the vortex breakdown with lead times longer than two weeks, and to distinguish onset signals for wave-1 and wave-2 events beyond the currently achieved predictable lead times. 100 The paper is organized as follows. Section 2 describes the data used in the analyses and the methodology behind the mode decomposition analysis. Section 3 shows the results of the analysis and their implications. Section 4 further provides the physical interpretation of the signals found in the mode equation budget by linking them with wave-mean flow interactions.
Conclusions are given in Section 5.
2 Data and methodology 105

Data and EOF basis
We use two data sets for analysis in this paper: 1) the ERA-Interim reanalysis (Dee et al., 2011), and 2) simulations from an intermediate complexity configuration of the Isca model (Vallis et al., 2018) (hereafter the Isca model). This version of the Isca model uses the model configuration from Jiménez-Esteve and Domeisen (2019), and it uses a T42 horizontal resolution and 50 vertical levels up to 0.02 hPa with 25 levels above 200 hPa. The model includes moist and radiative processes through 110 evaporation from the surface and fast condensation. Water vapour in the atmosphere interacts with a multiband radiation scheme (Mlawer et al., 1997) and a simple Betts-Miller convection scheme (Betts and Miller, 1986). The CO 2 concentration is fixed at 300ppm and the seasonal cycle of ozone in the stratosphere is prescribed based on the ERA-Interim (1979-2016 climatology. For the lower boundary conditions, the Isca model uses realistic topography and the continental outline from the ECMWF model, and sea surface temperatures (SST) are prescribed. The model does not include a representation of clouds, interactive 115 chemistry, or gravity wave drag. In this paper, we use the experiment that uses prescribed strong El Niño-like SST anomalies as described in Jiménez-Esteve and Domeisen (2019). The motivation to use this model experiment is that it produces a realistic climatology and a frequency of SSWs that is similar to the reanalysis.
For ERA-Interim, we use daily mean fields of potential vorticity (PV) P , zonal wind u, and meridional wind v at the 850 K isentropic level from 1979-2018 with a horizontal resolution of 2.5 • × 2.5 • . Only data north of 30 • N (24 latitude values by 120 144 longitude values leading to a total of D = 3456 grid points) from the winter season (October to April, 8490 days in total) are included in the analysis. For the Isca model simulations, the daily vertical gradient of potential temperature (θ), zonal wind u, and meridional wind v are interpolated to the 850 K isentropic surface from pressure levels. The Isca model data contains a total of 130 years, corresponding to 27300 winter days. The PV is computed as

125
where P is Rossby-Ertel's PV (Hoskins et al., 1985), ζ θ is the relative vorticity on the isentropic surface, θ is the potential temperature, f is the Coriolis force, g is the gravity, and p is the pressure. Note that in this work potential vorticity (PV) refers to Rossby-Ertel's PV.
The PV data (for both the ERA-Interim reanalysis and the Isca model) is further decomposed into: 1) the daily climatology, obtained by computing the daily mean values of PV over all available years, and 2) daily anomalies with respect to the 130 climatological seasonal cycle. Given the limited number of years in the reanalysis, we also applied a 30-day running mean to the daily climatology to remove the high-frequency variability and repeat the analysis performed in the study. The results are almost identical to the ones using climatologies without low-pass filtering, which is mainly due to the fact that applying PCA acts as a filtering and indicates the importance of the low-frequency components in the evolution of the polar vortex. In the following figures, we show the results using climatologies without low-pass filtering. The EOF modes of the PV and their 135 corresponding PC time series (later used in the mode decomposition analysis) are obtained by employing principal component analysis (PCA) on the daily PV anomalies at 850 K. The procedure consists of applying PCA twice, as described in the following. We apply a first PCA only to the PV data around the onset date of all SSWs (from -10 to +5 days around the onset date). The motivation behind this first PCA is to obtain a mode that captures the characteristics of SSWs. The spatial pattern of the resulting first EOF mode (E 1 ∈ R D ) is shown in the first panel in Figure 1, with a wavenumber 0 structure centered at 140 the pole. Next, we project the whole winter data (October to April) onto the subspace orthogonal to E 1 by subtracting from the winter data its projection onto E 1 , and then perform a second PCA on this projection. The resulting data does not contain any information from E 1 since it is in the space orthogonal to E 1 , and yields a total of D = 3456 modes {E 2 , . . . , E D+1 } (equivalent to the number of grid points). Combining the first EOF mode E 1 from the first PCA with the D EOF modes from the second PCA forms an orthogonal basis for the whole winter data, referred to in the following as the "combined set of basis 145 vectors". Figure 1 shows the spatial pattern of the first 10 EOF modes that explain together ≈71% of the variance of the PV anomalies of all winter days in the ERA-Interim data. The spatial patterns of the first 10 EOF modes of the PV anomalies in the Isca model data is shown in Figure C1 in Appendix C. We note that the EOF modes cannot be interpreted by default as physical modes as discussed also in Dommenget and Latif (2002) and Monahan et al. (2009). Here, since E 1 is derived using only days around SSWs (-10 to +5 days around the onset date of SSWs) and displays a clear zonally symmetric structure, the 150 physical process represented by the variation of E 1 can be interpreted as representing the changes in the strength of the polar vortex during SSWs. In ERA-Interim, the variance explained by E 1 (15.8%) for the whole winter data is slightly smaller than that of E 2 (16.8%), as shown in Figure 1. The EOF modes are orthogonal and therefore satisfy the following relationship of the inner product ·, · ,

155
where δ mn is the Kronecker delta function, and Ω is the region poleward of 30 • N. As mentioned above, the main motivation for applying the double PCA approach to obtain the EOF basis rather than directly applying PCA once to all the winter days is that E 1 is computed from the days around the onset day of SSWs, and its variability is therefore closely linked to SSW events, as can be seen from Figure 2. The PC time series of all winter days associated with E 1 (Figure 2) is highly correlated with the polar-cap averaged temperature at 30 hPa, thus being a good representation of SSWs (Blume et al., 2012). The evolution of the 160 first PC (for all winter data) therefore enables us to better understand the vortex breakdown process.

Mode decomposition analysis
Following the methodology from Aikawa et al. (2019), the mode decomposition analysis is applied to the PV conservation equation to study the dynamical development of SSW events. The anomalous PV equation on an isentropic surface can be written by separating the daily anomalies from the daily climatological mean for 1979-2018 as where P is the PV, V V V is the wind vector field, and F is the forcing term (computed here as residual). The subscript c refers to daily climatology, and a refers to daily anomalies. The anomalous PV tendency ∂Pa ∂t equals the sum of linear effects that consist of the advection of daily anomalous PV by climatological wind vector (V V V c · ∇P a ) and the advection of climatological PV by the anomalous wind vector (V V V a · ∇P c ), and the nonlinear effects of the advection of anomalous PV by the anomalous 170 wind vector (V V V a · ∇P a ).
For the mode decomposition analysis of the PV equation, we focus only on the time frame between 50 to 1 days prior to the onset day of SSW events as the goal of this study is to identify signals that precede SSWs. We project the PV field (obtained by concatenating the data from -50 to -1 days before the onset of all SSWs) onto the first d =1000 modes of the combined set of basis vectors {E 1 , E 2 , . . . , E d } and get the corresponding PC time series, denoted as {A 1 , A 2 , . . . , A d }. We found that the 175 first d =1000 modes (which in total explain 99.9% of the variance of PV as shown in Figure A1a in Appendix A) are sufficient to reproduce the actual rate of change of the PCs. The PV daily anomalies (P a ) can be expressed as the linear combination For the wind vector daily anomalies (V V V a ), the temporal evolution of PV is correlated with that of the wind fields. Therefore, it is possible to obtain a set of spatial patterns

185
where C k = A k , A k is the eigenvalue associated with mode E k ; L A kn = E k ,V V V c · ∇E n and L B kn = E k ,U U U n · ∇P c are the inner products between the linear advection terms (V V V c · ∇E n and U U U n · ∇P c ) and E k ; N kmn = E k ,U U U m · ∇E n is the inner product between the nonlinear advection term (U U U m ·∇E n ) and E k ; and F k = E k , F a is the residual term. The sum of the two linear advection terms gives the total linear advection term. Using Eq. (4), we then compute the contribution of each mode to the linear and nonlinear advection terms and thus to the total rate of change of A k to determine which modes (or combinations 190 of modes) play an important role in identifying signals distinguishing SSWs from normal winter days.
From the power spectrum of each PC time series A k (see Figure A1b in Appendix A), the power of A 1 to A 25 is concentrated at frequencies lower than once a week, which is different from the power spectrum of the other PCs beyond A 25 . Based on these power spectra, here we consider the associated modes E 1 to E 25 as low modes, which together explain around 85% of the total variance, and modes beyond E 25 as high modes. To separate the contributions from low and high EOF modes, the 195 summation over all modes from Eq. (4) is divided into the summation of low modes and that of high modes. Thus, Eq. (4) can be written as where l = 25 represents the nth mode that is treated as the last low mode; L low and L high are the linear low and high mode terms, respectively; N low−low , N low−high and N high−high are the contributions from nonlinear low-low, low-high, and highhigh mode interactions, respectively; L total is the total linear PV advection as L total = L low + L high and N total is the total nonlinear PV advection as N total = N low−low +N low−high +N high−high . We performed sensitivity tests for a range of l values, 205 and our results and conclusions are robust for values of l > 25. Moreover, the contribution from modes higher than l =25 is very small (not shown). In this study, we only focus on the rate of change of the first PC time series dA1 dt (k = 1) due to the strong relationship between A 1 and SSWs.

SSW definition
The major SSW definition follows the criterion of Charlton and Polvani (2007), based on the reversal of the daily zonal-mean 210 zonal winds at 10 hPa and 60 • N, and a return to westerlies afterward for at least 10 consecutive days before the final warming.
Two events in the same season are treated as distinct SSWs if they are separated by at least 20 days. The central dates of split and displacement SSW events in ERA-Interim before 2014 are taken from Karpechko et al. (2017) and consist of 11 split events and 12 displacement events. In addition, we added the SSW event in 2018, which is classified as split event based on the vortex geometry (Charlton and Polvani, 2007). The definition of wave-1 and wave-2 SSW events is based on the eddy heat 215 flux at 100 hPa and 60 • N similar to Bancalá et al. (2012). If the wave-2 component of eddy heat flux is larger than the wave-1 component by 15 Kms −1 in the period of -2 to 0 days before the SSW event for at least 1 day, then the SSW event is classified as a wave-2 event, otherwise as a wave-1 event. Note that the time window around the SSW event (day -2 to day 0) that is used to classify the type of event is shorter than that in Bancalá et al. (2012). The reason for using a shorter window is to reduce the overlap between the time interval used to define the type of SSW event and the lead times that emerge as relevant in 220 the predictability of wave-1 vs. wave-2 events. According to this definition, there are 18 wave-1 and 7 wave-2 SSW events in ERA-Interim. Note that not all split events are dominated by wave-2 wave flux (only 6 out of 12 split SSWs are also classified as wave-2 events), while one displacement event is dominated by wave-2 wave flux according to this definition.

Interpretation of results from mode decomposition based on PV flux
Up to this point, our mode decomposition analysis does not employ any explicit approximation. In this section, we will demon-225 strate how L total and N total are related to the poleward PV flux, which can help to illustrate the physical interpretation of the results obtained from the mode decomposition analysis.
In order to provide a physical interpretation for each term in Eq. (5) and understand the physical process that can lead to the disruption of the polar vortex, we introduce the poleward flux of PV on an isentropic surface given by Tung (1986) as

230
where F is the Eliassen-Palm flux (EP flux), P is the potential vorticity, v is the meridional wind, ρ θ is the density in isentropic coordinates, defined as ρ θ = − 1 g ∂p ∂θ , and φ is the latitude. The brackets denote the zonal mean, one asterisk denotes the deviation from the zonal mean, and two asterisks denote the deviation from the density-weighted zonal average, as in The formulation of the PV flux in Eq. (6) is also equivalent to that defined using v * * , as in Eq. (4.5) in Tung (1986). The 235 left-hand term in Eq. (6) is the EP flux pseudo-divergence and the right-hand term in Eq. (6) is the zonal-mean northward flux of PV on the isentropic surface. According to Tung (1986), the PV flux corresponds to the pseudo-divergence of the EP flux along isentropic surfaces and acts as the net eddy forcing term of the mean flow. Note that the pseudo-divergence of the EP flux is used here (instead of the divergence) due to the fact that the density on isentropic surfaces changes with time, which is the main difference from the conventional EP flux (Edmon et al., 1981).

240
Using the concepts of PV flux and EP flux pseudo-divergence, we rewrite Eq. (3) as The second to fourth terms on the left-hand side are the linear and nonlinear terms of the density-weighted PV flux divergence, while the last term is the local density tendency. Since the density-weighted PV is proportional to vorticity, the second to fourth terms can be interpreted as the dynamical contribution to the PV evolution. On the other hand, since the density tendency is 245 inversely proportional to the temperature tendency, the fifth term can be interpreted as the thermodynamic component of the PV evolution. Equation (8) is derived by converting the advection terms in Eq.
(3) into density-weighted flux divergence using the continuity equation. Taking the inner product between Eq. (8) and E 1 and neglecting the longitudinal variation of E 1 given its wavenumber-0 structure, i.e., E 1 ≈ E 1 (φ), we obtain an approximation of the linear and nonlinear terms in equation Eq. (8) (a detailed derivation is provided in Appendix B): where v is the meridional wind, a is the radius of the Earth, and φ is the latitude with φ 1 = 30 • N and φ 2 = 90 • N. Next, we combine the mode equation budget from Eq. (5) with the zonal-mean PV flux from Eq. (6). Based on the relation shown in Eq. (6) and the fact that ∇ · F ∝ ∂[u] ∂t from the zonal momentum equation, one can further obtain the following relation 255 Equation (10) shows that the meridional gradient of zonal- ∂φ is connected to the vorticity tendency ∂ζ θ ∂t , which is the dynamical component of the rate of change of A 1 .

Results of mode decomposition: mode equation budget
As shown in Figure 1, the first EOF spatial pattern E 1 of the PV daily anomalies in ERA-Interim takes the shape of a wavenumber-0 structure with a negative anomaly at the pole. Thus, a positive 1 (negative) value of A 1 (PC time series of 260 E 1 ) indicates a weakening (strengthening) of the polar vortex. For example, Figure 2a shows the corresponding A 1 for all winter days of ERA-Interim (8490 days in total) and the red vertical lines indicate the onset day of the SSW events. Before SSW events occur, A 1 increases significantly and is strongly positive on the onset day, indicating a weakening and breakdown of the polar vortex. Similar EOF spatial patterns are found for the Isca model data ( Figure C1 in Appendix C), together with a similar increase in A 1 when approaching the SSW central day ( Figure 2b). The Isca model data consists of a total of 130 years shown as an example. The PC time series for winter days of the remaining 90 years in the Isca model can be seen in Figure C2 in Appendix C. The red vertical lines indicate the onset dates of SSW events. that are outside the 2.5th to 97.5th percentile range of normal winter days values, which is computed via a bootstrapping procedure described as follows. In each bootstrap sample, we randomly select with replacement 25 sets of 50 consecutive non-SSW winter days (excluding the 50 days before each SSW) across the different years. We then calculate the mean of these 25 sets of non-SSW days to represent the "composite" of normal winter days. We repeat the bootstrap resampling procedure B =1000 times and compare the SSW composite against the 2.5th and 97.5th percentile of the B bootstrap samples. The same 280 procedure is applied also for the two types of SSWs and for the Isca model data, using the number of SSWs in each dataset as the number of sets of consecutive non-SSW winter days in the bootstrapping. If we only select 50 consecutive non-SSW winter days in December, January, February, and March (months when SSWs occur) to reflect the temporal distribution of SSWs, we obtain qualitatively similar results as when using non-SSW days for all winter months. Here we present the results with the bootstrapping using non-SSW days for all winter months. The reconstructed A 1 tendency (black line in Figure 3) is computed 285 from Eq. (5). The left panels show each term in Eq. (5), while the right panels show the combined effect of all linear (red) and all nonlinear (green) terms without separating the contributions from low and high EOF modes. Figure 3a indicates that around 25 days before the central day of SSWs dA1 dt starts to increase, and the increase shows a steeper slope around 10 days before the SSW event, leading to a large positive A 1 on the central day as shown in Figure 2a. Along with the increase of dA1 dt , L low (red) also increases and is well correlated with the A 1 tendency (r = 0.8). In fact, L low starts to increase from 35 days before 290 the events but its effect is offset by other terms and it becomes the only contributor to the increase of dA1 dt at 25-to 15-day leads. The nonlinear term N low−low (green) shows a rapid increase around two weeks before the SSW event and, together with the linear term L low , significantly contributes to the changes in A 1 . The high-frequency components are overall weaker than the low-frequency terms, especially the N high−high (cyan), but the N low−high (orange) has large variations and tends to offset the effect of L low−low at -25 to -15 days before the vortex weakening.

295
The contributions from each term are different between the composites of wave-1 vs. wave-2 events (Figures 3c and 3e (wave-2) events, while the nonlinear (linear) terms are less important or even counteracting the increase of dA1 dt . Therefore, the relative importance of the linear and nonlinear terms emerges as a good indicator of the type of SSW events with a lead time of around 10 days prior to the events.
The relative importance of the linear and nonlinear advection terms for the two types of SSW events is similar to that of the stratospheric wave amplitude of the wavenumber 1 and 2 components for wave-1 and wave-2 SSW events as shown in Bancalá 315 et al. (2012). From their composite analysis of the wave-2 SSW events, the wave-2 component of the geopotential height anomaly in the stratosphere is significantly positive during day -10 to day 0, and the anomalous increase of wave-1 component was found in the period from day -30 to day -10. In Figure 3 we find that the linear and nonlinear advection terms behave similar to the wave-1 and wave-2 components of the upward propagating wave activity. Consistent with the two periods in from day -25 to day -15 and the other from day -10 to day 0. Since the displacement (split) SSW events are mainly attributed to the enhanced upward propagation of wave-1 (wave-2) (Nakagawa and Yamazaki, 2006;Bancalá et al., 2012), we also look into the contributions of linear and nonlinear terms to dA1 dt for displacement and split SSW events (shown in Figure D1 in Appendix D), with very similar results. Comparing Figures D1a,b with Figures 3c,d, Figure 3d and 3f (for displacement and split SSWs, see Figure D1) is comparable to the results in 330 Smith and Kushner (2012), who showed that displacement (split) SSW events are preceded by pronounced linear (nonlinear) vertical wave activity. Our results suggest that the linear and nonlinear contributions are more strongly related to the dominant wavenumber wave forcing than to the vortex geometry.
In order to examine the significance in the differences and the robustness of the relative importance of the linear and nonlinear advection terms for the two types of SSW events, we perform bootstrapping on individual wave-1 and wave-2 events with 335 replacement, respectively. We repeat the resampling B =1000 times and compute the means and the standard deviation for the sum of linear and nonlinear terms. The results of the bootstrapping are shown in Figure 4. There is almost no overlap between the ±1 standard deviation of wave-1 (red) and wave-2 (black) events in neither the total linear advection (Figure 4a) nor the total nonlinear advection (Figure 4b) term at one week before the onset of the events. In particular, the separation of the wave-1 and wave-2 events in the total linear advection term is as early as 10 days before the onset of the events. Figure 4 demonstrates 340 the significance and robustness of the differences in the contribution of the linear and nonlinear advection terms to dA1 dt in the wave-1 and wave-2 SSW events, respectively, at least up to one week before the events. Another point to highlight is that significant anomalies of the linear terms are observed around 20 days before both types of SSW events, which is beyond the current predictability limit of SSWs of one-two weeks.

345
Given the limited number of SSW events in the reanalysis data and to further examine the characteristics and robustness of the linear and nonlinear terms contributions to the vortex breakdown, we now apply the same analysis as for ERA-Interim to the output of the Isca model experiment. We use the methodology from Section 2 to extract the EOF modes (spatial patterns) and apply the mode decomposition analysis to the data concatenating 50 to 1 days prior to the 78 SSWs present in the model data.
The EOF spatial patterns derived from the model output are similar to those derived from ERA-Interim, especially the first 10 350 EOF modes ( Figure C1), indicating that the model is able to capture the PV features as in the reanalysis. standard deviation of a bootstrapping using B =1000 samples. Figure 5 shows the results of the A 1 budget for the SSW composites. Similar to the results in ERA-Interim (Figure 3), the linear term L low starts to increase at around day -25, but the increase in dA1 dt starts at around day -10 (Figure 5a), which is later than that in ERA-Interim (around day -25). In the period of day -10 to day 0 of SSWs, N low−high increases rapidly.
The increasing N low−high and L low lead to a rapid increase of dA1 dt . We note that in the 2 days before the onset date, the 355 magnitude of the wave-1 and wave-2 wave fluxes is similar (with differences < 5Kms −1 ) in some SSW events simulated in the Isca model. In order to make a clearer separation between wave-1 and wave-2 SSWs, we exclude events with the difference in the magnitude of wave-1 and wave-2 heat flux at 100 hPa and 60 • N smaller than 5 Kms −1 as these events cannot be clearly categorized as either wave-1 or wave-2 events. Therefore, in the model we have 36 wave-1, 27 wave-2, and 15 unclassified events. Figure 5c shows that L low for SSWs increases and starts to differentiate from normal winter days starting at a lead of 360 20 days for wave-1 SSWs. Different from the A 1 budget of wave-1 SSWs in ERA-Interim (Figure 3c), N low−high starts to increase from day -7 and becomes an important contributor to dA1 dt . For wave-2 SSWs, Figure 5e shows that N low−high is the main contributor to the increase of dA1 dt from day -10, and N low−low as well as L high are the second largest contributors to dA1 dt from day -5, which is different from the evolution of L high for wave-2 SSWs in ERA-Interim (Figure 3e). In both types of events, L low starts to increase at around day -20, which helps to weaken the polar vortex in the preconditioning stage and is similar to the evolution of L low (with a smaller amplitude) in the same period in ERA-Interim. The effects of L total and N total are shown in the right panels in Figure 5. The evolution of L total and N total , and thus of dA1 dt , for the Isca model ( Figure  5b) is similar to the evolution of these terms for the ERA-Interim data (Figure 3b), indicating that the Isca model successfully reproduces the vortex breakdown in the 10 days preceding the SSWs. The increase of dA1 dt can therefore be used to predict the occurrence of SSWs with one-two weeks lead time. Even though the distinct increase of dA1 dt only shows up at around day -10, L total actually increases as early as day -29. However, this amplification in the linear term is offset by the nonlinear and forcing terms, which leads to a near-zero dA1 dt . Figures 5d and 5f show the evolution of L total and N total of the wave-1 and wave-2 SSW composites, respectively. Different from wave-1 events in ERA-Interim, the linear and nonlinear terms are equally important. However, when comparing wave-1 with wave-2 composites, one can still see the difference in the relative importance of the linear and nonlinear terms for the two types of SSWs. The nonlinear term is stronger in wave-2 SSWs than in 375 wave-1 SSWs and is more than twice as large as the linear term one week before the central day (Figure 5f). Thus, our finding from the ERA-Interim that the linear (nonlinear) term is important for wave-1 (wave-2) events is also true for the Isca model data. The main differences compared with the reanalysis data are that dA1 dt exhibits a substantial increase only from day -10, and the variations of all terms in Eq. (5) are overall small before day -10, which could potentially limit the predictability of SSWs in the Isca model. Different reasons might be able to explain the differences in results between the Isca model and the 380 reanalysis, e.g., different model complexities (e.g., lack of parameterizations for gravity waves breaking and interactive ozone chemistry in the Iscal model) or the coarse model horizontal resolution (T42), both of which might lead to an underestimation of some of the high-frequency variability in comparison with the reanalysis. Another important difference is that in the Isca model the nonlinear term displays larger values for both types of SSWs compared to the reanalysis. This latter behavior might be related to the stronger SSW-sensitivity to wave-2 forcing in the idealized models (e.g., Iscal model) where most of the SSW 385 are likely triggered by wave-2 activity (Gerber and Polvani, 2009), in contrast with wave-1 which seems to be more dominant in more complex general circulation models or reanalysis datasets (see for example Figure 3 in Ayarzagüena et al., 2018). Note that even when classifying an SSW event as a wave-1 event in the model, its wave-2 component, although weaker than the wave-1 component, might still play an important role in the overall evolution of the event.
4 Physical interpretation of the mode decomposition 390 In our analyses above, we found that the persistent positive values of dA1 dt and its contributors that emerge during the vortex breakdown, such as during SSWs, are significantly different from the values observed during normal winter days. Additionally, the signals identified as representative of wave-1 and wave-2 events are also different. We also observed that the signals that are characteristic of SSWs emerge as early as 20-25 days before the onset of SSWs. Given that these results hint that SSWs are potentially predictable at longer lead times, i.e., beyond the current predictability limit of one-two weeks, in this section 395 we provide a physical interpretation of these signals that we identified through the mode decomposition analysis.
As demonstrated in Section 2.4, the total linear and nonlinear advection terms in Eq. (5) are closely linked to the PV flux divergence, which offers a more intuitive interpretation in an Eulerian framework. In this section, we use ERA-Interim data to illustrate the physical interpretation of the increase of the linear and nonlinear advection terms in Eq. (5). The motivation to introduce the PV flux into the PV equation is that the zonal-mean PV flux is connected to the pseudo-divergence of the EP flux 400 along isentropic surfaces, and thus acts as the net eddy forcing term of the mean flow, thus allowing for connections with the theory of wave-mean flow interaction (i.e., Eq. (6) and Eq. (10)). According to McIntyre and Palmer (1983), the wave activity of planetary waves is converted to the angular momentum of the mean flow, which violates the non-acceleration condition (Charney and Drazin, 1961), leading to the reversal of the mean flow. To understand the importance of the zonal-mean PV flux during the development of SSWs, we decompose the zonal-mean PV flux into different components as in Ayarzagüena et al. 405 (2011), where the subscript c represents daily climatology, and a represents daily anomalies. On the right-hand side of Eq. (11), the first term corresponds to the climatological planetary waves, the second and third terms correspond to the interaction between the climatological planetary waves and the daily anomalies, and the fourth (last) term corresponds to the interaction between 410 daily anomalies. Similar to Eq. (3), the second and third right-hand terms can be viewed as linear components, and the fourth term as the nonlinear component. Figure 6 shows the composite of zonal-mean poleward PV flux averaged north of 45 • N and its decomposition in Eq. (11) as a function of lead time ahead of SSW events. The first term on the right-hand side of Eq. (11) can be seen as a constant since its variation with time is very small as shown in Figure 6. When approaching the onset day of SSWs, the zonal-mean PV flux (black) becomes increasingly negative, indicating a weakening of the polar vortex. This 415 further decrease of the negative zonal-mean PV flux is mainly due to the linear interaction between the climatological planetary waves and the daily anomalies (red) and the nonlinear interaction between anomalies (green), which correspond to the total linear and nonlinear PV advection terms (right columns in Figure 3), respectively. Even though the climatological planetary waves (blue) also have negative contribution to the total PV flux, the variations are very small with time. Similar to the distinct contributions of the linear and nonlinear PV advection terms in the wave-1 and wave-2 SSW composites, the negative total 420 zonal-mean PV flux is mainly due to its linear component in wave-1 SSWs (red in Figure 6b) and its nonlinear component in wave-2 SSWs (green in Figure 6c) of the warming (day -25 to day -15), the nonlinear PV flux has a very low magnitude as shown in Figure 7g. Since previous studies suggested that split SSWs have a predominantly barotropic structure (Manney et al., 1994;Matthewman et al., 2009;Albers and Birner, 2014), we investigate the vertical structure of the nonlinear PV flux in the wave-2 SSWs composite. Figure  8b shows the longitude-height cross section of the nonlinear PV flux of the wave-2 SSWs composite at day -5 as an example.
As can be seen in Figure 8b, the wave pattern shown in Figure 7h extends throughout the stratosphere and displays barotropic characteristics for wave-2 events. On the other hand, the longitude-height cross section of the linear PV flux at day -5 of the wave-1 SSWs in Figure 8a displays a more baroclinic structure in the Eurasia and Pacific regions. The spatial pattern of the nonlinear PV flux in the composite of wave-1 events exhibits substantial transient fluctuations without a clear wave pattern 460 before day -9 (Figure 7e), and shows a more stable and organized spatial pattern in the period from day -9 to day 0 ( Figure 7f).
However, the magnitude of the nonlinear PV flux in wave-1 events is smaller than its linear flux counterpart and also smaller than the nonlinear PV flux in the wave-2 events composite.
Even though the spatial patterns of the linear PV flux in the two types of SSWs show a wave-2 pattern in the period of 20 to 10 days preceding the central day of SSWs, the locations of the maximum and minimum PV flux shift around 30 • in longitude in 465 the Pacific and North America regions (Figures 7a and 7c). In the first period from day -10 to day 0, the wave pattern shown in Figure 7h sets in and leads to the final split of the vortex for wave-2 SSWs. One relevant question is what processes lead to the differences observed in the evolution of the vortex breakdown where the linear PV flux remains important for wave-1 events, while the nonlinear PV flux amplifies for wave-2 events. Some previous studies suggested that the pre-SSW evolution of the polar vortex is distinct between split and displacement events (Charlton and Polvani, 2007;Bancalá et al., 2012), and 470 this preconditioning could trigger the nonlinear resonance of planetary waves in the lower stratosphere, leading to the split of the polar vortex (Albers and Birner, 2014;Boljka and Birner, 2020). Here we examine the anomalies of PV and meridional wind after removing the daily climatology and zonal mean (P * * a and v * a , respectively, see Eq. (11)) for the wave-1 and wave-2 events to understand their distinct evolutions after day -10. Figure 9a-b show the spatial pattern of P * * a (shading) and v * a (green contour) at day -15 before wave-1 and wave-2 events, respectively. Both P * * a and v * a present wave-1 patterns, but the positive 475 and negative anomalies are located in different regions. The whole pattern of P * * a in wave-2 SSWs is around 60 • further east compared to wave-1 SSWs. The negative v * a is mainly located over eastern North America and the northern North Atlantic, which is important for the negative PV flux in the same region ( Figure 7a) for wave-1 SSWs, while for wave-2 SSWs the negative v * a covers all of North America. These differences in the location of P * * a (shading) and v * a (green contour) between the two types of SSW events are amplified from day -10 (Figure 9c-d). The magnitudes of P * * a and v * a in the first period from 480 day -10 to day 0 are larger than in the second period from day -25 to day -15. The negative P * * a is located more over North America, while the positive P * * a is located over the North Atlantic and Europe for wave-1 SSWs (Figure 9c). For wave-2 SSWs, the pattern of P * * a is the opposite (Figure 9d). The positive v * a extends to the full North Pacific (Figure 9c) and both P * * a and v * a maintain wave-1 structure for wave-1 SSWs. In contrast, v * a (Figures 9d) develops a wave-2 structure from day -10 onwards for wave-2 events. The weak positive P * * a over Asia (Figure 9d) further develops from day -5, resulting in a wave-2 structure 485 over mid-latitude at day 0 (not shown). The main features of nonlinear PV flux in Figure 7 can be roughly inferred by P * * a and v * a in Figure 9. We note that the composite of the nonlinear PV flux term is not equal to the direct product of the composites of P * * a and v * a , and thus some features in Figure 7 cannot be directly inferred from Figure 9. We also find that the nonlinear PV flux and the P * * a and v * a form a positive feedback from around day -10 to the onset of the wave-2 events. As the amplitude of P * * a and v * a becomes larger, the nonlinear PV flux is also amplified, particularly in the region of the negative nonlinear PV flux 490 over western North America and the North Atlantic as shown in Figure 7h. The strong negative nonlinear PV flux contributes to more negative net zonal-mean PV flux values, suggesting a zonal-mean EP flux convergence. This EP flux convergence thus further decelerates the polar vortex. According to the non-acceleration theorem (Charney and Drazin, 1961), the deceleration of the polar vortex is accompanied by stronger wave activity, which is represented by the increasing amplitude of P * * a and v * a . We also note that the spatial pattern of P * * a in Figure 9f finally leads to the split of the polar vortex in wave-2 events, with the 495 positive values corresponding to one of the daughter vortices located around 60 • W.

Conclusions
In this paper we employ a mode decomposition analysis to investigate the preconditioning of sudden stratospheric warming events. We study the (linear and nonlinear) terms in the potential vorticity equation by means of a budget analysis in order to identify the components in the first PC time series A 1 that allow us to distinguish the behavior of the polar vortex during 500 SSW events from normal winter days. Moreover, we identify characteristics of SSWs that help to identify the type of event (wave-1 vs. wave-2) during the dynamical development of SSWs. The mode decomposition analysis allows us to obtain a mode equation budget that describes the temporal evolution of the stratospheric dynamical processes that lead to the breakdown of the polar vortex. A better understanding of the vortex weakening process may help to improve the predictability of SSW events.
The rate of change of the first PC time series dA1 dt represents the evolution of the strength of the polar vortex, and we find 505 a significant increase in dA1 dt at around 25 days before the onset of SSWs. This change in dA1 dt marks the start of the vortex weakening process indicating an acceleration of the polar vortex breakdown, and is different from the evolution of dA1 dt during normal winter days. The lead time of 25 days that we identified in our analysis is far beyond the current predictability limit of SSW events (Domeisen et al., 2020b). We note that not only the composite of SSWs, but most of the individual SSW events show the increase of dA1 dt at around 20-25 days before the onset. The increase of dA1 dt is mainly due to the increase of the linear 510 PV advection term, which preconditions the weakening of the polar vortex. While recent work suggests that split SSW events are less predictable than displacement events (Taguchi, 2018;Domeisen et al., 2020b), the preconditioning by the increase in the linear PV advection and PV flux is important for both wave-1 and wave-2 SSW events at around 20 days before the onset of the event, implying a similar intrinsic predictable time scale with 20 days lead for both types of events. From around 10 days before the events, the nonlinear PV advection term increases rapidly for wave-2 SSWs, but it remains small for wave-515 1 SSWs. As the nonlinear PV advection term increases, the linear PV advection drops dramatically prior to wave-2 SSWs.
The distinct behavior of the linear and nonlinear advection terms in this 10-day period suggests that the type of SSW event could be inferred at around 10 days to one week prior to the events. Note that the type of event is determined by the larger wavenumber component of eddy heat flux in the period of -2 to 0 days before the events and hence the 10-day lead times need to be interpreted with caution. The above differences are also present in the displacement and split SSW events, but the 520 differences are somewhat smaller than those between wave-1 and wave-2 events, as not all split events are induced by wave-2 planetary waves.
Even though the contributions from linear and nonlinear PV advection terms are different in the two types of SSW events, their overall effects on dA1 dt within 10 days before the events are the same, causing dA1 dt to increase abruptly. The breakdown of the polar vortex can be divided into two periods based on the different behavior of the linear and nonlinear terms. During the 525 first period, i.e., from 25 days to two weeks before the onset of SSWs, the linear term weakens the polar vortex for both types of SSWs, and during the second period, i.e., from around 10 days before the onset date until the onset, the vortex evolution for both types of SSWs starts to diverge and the distinct vortex breakdown structures gradually develop. These two different time periods before the onset of the events are consistent with previous studies, especially for wave-2 SSWs (Labitzke, 1981;Bancalá et al., 2012;Albers and Birner, 2014), which suggested that an amplification of the wave-1 component allows the 530 wave-2 wave flux to grow and propagate more effectively into the already weakening polar vortex region.
In both the ERA-Interim reanalysis and the simplified Isca model experiments, the increase of dA1 dt is more abrupt for wave-2 SSW events than for wave-1 events and thus results in a larger dA1 dt in the 10-day period preceding the onset of the events. The abrupt changes in dA1 dt are mainly due to the exponential increase in the nonlinear PV advection term. By contrast, the linear PV advection for wave-1 SSWs increases more slowly but consistently. The rapid growth of the nonlinear process for 535 wave-2 SSWs could be related to a positive feedback between the nonlinear PV flux and the anomalies of PV and meridional wind when we tried to interpret the underlying dynamics of the increase in the nonlinear PV advection obtained from the mode equation budget. The linear and nonlinear advection terms are closely linked to the PV flux divergence, while the zonal-mean PV flux can be directly related to the zonal mean momentum budget (McIntyre and Palmer, 1983;Tung, 1986;Plumb, 2010).
The zonal-mean poleward PV flux can be further decomposed into the linear and nonlinear components, whose role in the 540 weakening of the polar vortex is similar to the effect of the PV advection terms on the increase of dA1 dt . The wave-2 spatial pattern of the linear PV flux helps to precondition the stratospheric basic state and decelerate the polar vortex in the first period of the SSW development for both types of SSWs. When the vortex weakening process evolves to the second period, the evolution of the PV flux for the two types of SSW events bifurcates as the linear and nonlinear PV flux exclusively amplify in the wave-1 and wave-2 events, respectively. This bifurcation could be due to the specific evolution of the stratospheric states in 545 the two types of events (Charlton and Polvani, 2007;Albers and Birner, 2014), which can be seen from the horizontal patterns of the PV and meridional wind zonal anomalies. Our results suggest that the high wavenumber pattern that emerged in the second period for wave-2 SSWs is closely connected to the wave-2 wave flux and could be essential to the split of the vortex.
As suggested by Aikawa et al. (2019), mode decomposition analysis allows us to investigate the contribution of each EOF mode to the breakdown of the polar vortex, and the way the associated spatial patterns play a role in the temporal evolution of 550 the first PC time series. We found that the interactions involving the low modes are the dominant contributors to the weakening of the polar vortex, especially for the increase of the linear advection term in the first period (i.e., day -25 to day -15). Further investigation of the contribution from each EOF mode to the linear and nonlinear advection terms in the mode equation budget suggests that the increase of the linear advection term in the first period is largely influenced by the second and third EOF modes, which both show a wave-1 structure. The first EOF mode only plays an important role at around one week before the of the nonlinear advection terms, the interactions amongst the first five EOF modes are important when approaching the onset of SSWs.
Even though the increase in dA1 dt and the contribution from the linear term start at around 25-20 days before the onset of SSWs in ERA-Interim, one needs to be cautious about the interpretation. The signals shown in the composite of SSWs do not 560 necessarily indicate that these signals can be used to predict each individual SSW event with lead times of 20-25 days. While most SSW events do show a consistent positive dA1 dt starting 20 days before the SSW onset, around 30% of all SSWs in ERA-Interim do not show this clear increase in dA1 dt around lead times of 20-25 days. On the other hand, the linear signals found here with lead times of 3-4 weeks may not be exclusive to SSWs. There are cases where large positive values of the first PC time series do not correspond to an SSW event, but instead to a strong deceleration event. However, dA1 dt and the contribution from 565 its linear term for SSW events are overall stronger and more persistent than that for the strong deceleration events (not shown).
What we found here suggests that the intrinsic predictability of SSWs may be longer than the current two-week practical predictability. However, more work is still needed to investigate whether the practical predictability of SSWs can actually be extended and if yes, how.
Since the signals shown here indicate that most SSWs may be predictable on subseasonal time scales, it is important to 570 understand which processes lead to the variability of the first EOF pattern and help to improve subseasonal forecast skill.
A recent study by Albers and Newman (2021) identified two modes that relate to the linear and nonlinear processes for strong downward propagating stratospheric anomalies, with one mode representing purely stratospheric processes and the other mode representing stratosphere-troposphere coupling. However, they point out that it is not clear which processes are more important for subseasonal predictability. Even though we have not connected the specific physical processes to the linear 575 and nonlinear signals that are important for the two types of SSWs in this study, the spatial patterns of the PV flux could potentially provide some hints for future investigation. For example, the wave-2 spatial patterns of the linear PV flux are relatively stationary for wave-1 SSW events. This is due to the fact that the orientation of the negative and positive anomalies does not change significantly from day -19 onwards and the orientation of the wave-1 structure also remains stationary in both PV and meridional wind anomalies. These persistent spatial patterns and the linear behaviour in the early stage of development 580 of SSWs may be related to weather phenomena, such as blocking, teleconnections, and low frequency modes in the troposphere.
For example, Smith and Kushner (2012) and Cohen and Jones (2011) suggested that displacement events are preceded by sea level pressure anomalies associated with the Siberian high which is consistent with the increase of linear vertical wave flux before the events.
In conclusion, our study finds signals that are representative of SSW events as early as 25 days preceding the events. This 585 lead time is significantly longer than the current predictability limit of SSWs. We furthermore find that mode decomposition analysis can help infer wave-1 and wave-2 events at least one week ahead of the event, which is longer than the lead times identified in previous studies (Karpechko, 2018;Taguchi, 2018;Domeisen et al., 2020b). The time scale of emergence of the distinct evolution between linear and nonlinear terms provides insights into the different dynamical processes responsible for noticeable increase in dA1 dt in the simplified GCM (Isca model) experiment, which directly indicates the weakening of the polar vortex, shows up only around 10 days before the onset of SSWs (i.e., at shorter lead times than for reanalysis), suggests that the observed atmosphere tends to be more predictable than the model, which agrees with theory Scaife and Smith, 2018). Applying the mode decomposition analysis to more complex forecasting models, i.e., S2S reforecast models (Vitart et al., 2017), to examine the predictability of SSWs will provide further insights into the dynamics of the polar vortex 595 weakening and might potentially allow for the prediction of these events beyond the current lead times.

600
In this appendix, we show the derivation procedure for obtaining the mode decomposition equation budget (Eq. (4)) and the cumulative explained variance and the power spectrum of the first 1000 EOF modes. The spatial patterns associated to the projections (U U U 1 ,U U U 2 , . . . ,U U U d ) of the wind vector daily anomalies V V V a onto the PC time series are computed by projecting V V V a onto the PC time series ({A 1 , A 2 , . . . , A d }). Note that U U U 1 ,U U U 2 , . . . ,U U U d are not necessarily orthogonal. The anomaly terms of PV and the wind vector fields can be written as Taking the inner product between Eq. (A2) and a given EOF mode E k , we obtain and δ kn = 0 for n = k), with C k being the eigenvalue of mode k, the mode equation budget is computed as which is the expression of Eq. (4) in Section 2.2.
Here we show the justification for the choice of truncation at the 1000th EOF mode and the choice of low modes. The first 615 1000 modes together explain ≈100% of the variance of the PV anomalies of all winter days in ERA-Interim. The powers of the first 25 EOF modes are concentrated in the period longer than one week. Therefore, we defined mode 1-25 as low modes. Appendix: B. Relation between the rate of change of A 1 and the PV flux In this appendix, we show the derivation for obtaining approximations of the linear and nonlinear terms (Eq. (9)) in the A 1 tendency equation using the PV flux form. Taking the inner product between Eq. (8) and E 1 , and neglecting the variation of 1 ρ θ in the inner products, we obtain the following approximation of the rate of change of A 1 : where C 1 is the eigenvalue associated with the first EOF mode of PV. Comparing the linear and nonlinear terms in Eq. (B1) with those in Eq. (A3), we can see that

625
As we mentioned in Section 2.3, given the wavenumber-0 structure of E 1 , we further approximate E 1 to be only a function of latitude. Thus, taking the inner product with E 1 can be approximated as taking a latitude-weighted integral of the meridional gradient of PV flux as demonstrated below: where a is the radius of the Earth, and φ is the latitude with φ 1 = 30 • N and φ 2 = 90 • N. Using Eq. (B3), we can obtain the 630 approximated linear and nonlinear terms from Eq. (9).
Appendix: C. EOF modes and PC time series of PV in the simplified Isca model In this appendix, we show the spatial patterns of the first 10 EOF modes of PV anomalies at 850 K and the associated first PC time series in the Isca model as described in Section 2.1. The first 10 modes together explain ≈82% of the variance of the PV anomalies of all winter days in Isca model.

635
Appendix: D. Mode decomposition budget for displacement and split SSWs in ERA-interim In this appendix, we show the A 1 mode decomposition budget for the composites of displacement and split SSW events in ERA-interim, which can be compared with Figure 3.   is the same as the legend in (a) and (b), respectively.