Influence of Northern Hemispheric Winter Warming on the Pacific Storm Track

Effect of global warming on the sub-seasonal variability of the Northern Hemispheric winter (NDJFM) Pacific storm-track (PST) activity has been investigated. Previous studies showed that the winter-averaged PST has shifted northward and intensified, which was explained in terms of energy exchange with the mean field. Effect of global warming exhibits spatio-temporal heterogeneity with predominance over the Arctic region and in the winter season. Therefore, seasonal averaging may hide important features on sub-seasonal scales. In this study, distinct sub-seasonal response in storm track activities to winter Northern Hemispheric warming is analyzed applying cyclostationary empirical orthogonal function analysis to ERA5 data. The key findings are as follows. Change in the PST is not uniform throughout the winter; the PST shifts northward in early winter (NDJ) and intensifies in late winter (FM). In early winter, the combined effect of weakened baroclinic process to the south of the climatological PST and weakened barotropic damping to the north is responsible for the northward shift. In late winter, both processes contribute to the amplification of the PST. Further, change in baroclinic energy conversion is quantitatively dominated by eddy heat flux, whereas axial tilting of eddies is primarily responsible for change in barotropic energy conversion. A close relationship between anomalous eddy heat flux and anomalous boundary heating, which is largely determined by surface turbulent heat flux, is also demonstrated.


Introduction
Weather changes we experience everyday are intimately linked with the passage of synoptic systems. Storm track is defined as the region where the activity of a synoptic system is particularly intense. Considering their extensive influences on human activities in association with extreme weather (Pfahl and Wernli 2012;Kunkel et al. 2012;Ma and Chang 2017), hydrological cycle (Fernández et al 2003;Sodemann and Stohl 2012), and large-scale circulation (Stahl et al. 2006;Lau 1978), many studies have focused on the variability of storm tracks.
Climatologically, two major storm track regions can be found in the Northern Hemisphere (NH). One is the Pacific Storm Track (PST) and the other is the Atlantic Storm Track (AST). The two storm track regions show different temporal variations with distinct three-dimensional spatial structures in the winter season. The intensity of the AST shows a single peak in January tied with strong baroclinicity. On the other hand, the PST shows relatively strong storm activities in late autumn and early spring compared to January, which is called the mid-winter suppression (Nakamura 1992). Deng and Mak (2005) showed that the sub-seasonal variation of the PST is due to a stronger effect of barotropic damping than baroclinic development in mid-winter compared to late fall/early spring. Lee et al. (2011) explained the midwinter suppression of the PST in terms of energetics. The southward shift of the Pacific jet stream during mid-winter sets apart the region of maximum meridional temperature Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s0038 2-020-05544 -4) contains supplementary material, which is available to authorized users.
* Kwang-Yul Kim kwang56@snu.ac.kr Hyung-Ju Park spalfh@snu.ac.kr 1 gradient and that of upward and poleward heat flux, which causes a weak baroclinic energy conversion (BCEC). The spatial and temporal characteristics of maximum BCEC also differ between the two storm track regions during the boreal winter. The maximum BCEC region of the PST is located in the upper troposphere (300 hPa) in late autumn/early spring, while that of the AST is located in the lower troposphere (900 hPa) in mid-winter. Storm tracks are affected by natural variability on various scales and latitudes. Madden Julian Oscillation (MJO;Madden and Julian 1971, 1972, which is associated with tropical convections in the Indian and the western Pacific, accompanies a shift of the PST on sub-seasonal time scales. When the center of MJO convection moves eastward across the eastern Indian Ocean and the western-central Pacific, the response of the PST is characterized by amplitude-varying dipole pattern of synoptic eddy kinetic energy propagating northeastward (Matthews and Kiladis 1999;Deng and Jian 2011;Lee and Lim 2012;Takahashi and Shirooka 2014;Wang et al. 2018). Furthermore, the MJO-related PST change in the northern North Pacific is affected by the tropical stratosphere, exhibiting a larger amplitude during the easterly phase of quasi-biennial oscillation (Wang et al. 2018).
Other tropical heating modifying the PST on interannual time scales is El Niño-Southern oscillation, which results in a significant equatorward and downstream shift during El Niño years and the opposite in La Niña years (Trenberth and Hurrell 1994;Straus and Shukla 1997;Zhang and Held 1999;Chang et al. 2002;Eichler and Higgins 2006). In midlatitudes, the PST is influenced by the East Asia winter monsoon (EAWM; Nakamura et al. 2002;Harnik and Chang 2004;Lee et al. 2010;Song et al. 2016). When the EAWM is strong, the subtropical Pacific jet becomes narrower and weak Pacific storm activities are observed (Harnik and Chang 2004). This strong jet-weak PST relationship is due to the modulation of baroclinicity by EAWM-induced lowerlevel cooling (Lee et al. 2010). The PST is also affected by atmospheric variability in polar regions. When the Arctic oscillation is in a strong positive phase, the AST becomes stronger and shifts northward, while the PST extends westward (Nie et al. 2008). When the North Atlantic oscillation is in a positive phase, Atlantic storm activities become stronger (Hurrell et al. 2003;Rivière and Orlanski 2007;Bader et al. 2011).
With increasing interests in climate change, many studies dealt with climatic impact on storm tracks. In a zonal mean sense, NH storm activities appear to expand northward and upward (Yin 2005;Wu et al. 2010;Chang et al. 2012;Woollings and Blackburn 2012). However, it seems that the PST and the AST respond to climate change in different ways. The PST shifts northward in response to climate change in most studies (Bengtsson et al. 2006;Wang et al. 2006;O'Gorman 2010;Bender et al. 2012;Chang et al. 2012;Wang et al. 2017). For the AST, however, some studies show weakening Zappa et al. 2013;Colle et al. 2015;Wang et al. 2017), while others show strengthening over the southern flank of the AST Harvey et al. 2015). It appears that the impact of climate change on the PST and the AST exhibits distinct spatio-temporal characteristics and mechanisms on seasonal to inter-annual scales.
The effect of warming is not uniform in terms of temporal and spatial distributions. Temperature rise is more remarkable in polar regions than in equatorial regions, and in winter rather than in summer (Serreze et al. 2009;Screen and Simmonds 2010;Kim and Kim 2018;Kim et al. 2019). Given this heterogeneity, it is expected that the storm tracks are affected non-uniformly by warming both in space and time.  reported enhanced PST in early spring and weakened AST in middle winter during the recent decades. However, detailed mechanism of how they are affected by warming on sub-seasonal scales has yet to be answered. In this study, we focus only on the PST, which earlier studies are generally in agreement to have shifted northward and intensified in response to warming.
Two methods, energetics and moist static energy (MSE) budget, have been used traditionally to analyze storm activities. Inspection of energetics allows detailed understanding of interactions between mean flow and eddies. Since it is based on the dry dynamics, however, contribution from moist processes cannot be accounted for. Nevertheless, many previous studies have shown that major change in storm characteristics coincides with change in BCEC. The MSE budget analysis allows an inspection of how heating, including radiation and turbulent heat flux at the atmospheric boundary, is balanced by horizontal transport of MSE by mean circulation and storms. Therefore, the two methods will allow complementary interpretations.
The goals of this study are as follows. We will show that change in the Northern Hemispheric PST due to warming is not uniform during the winter season but shows a noticeable inter-seasonal variation, which is characterized by a northward shift during the early winter and a strengthening during the late winter. Distinct physical processes during the two stages of winter will be investigated based on the eddy-mean flow interaction and energy budget. Change in the AST also exhibits seasonal difference with intensification in midwinter and weakening in other months. The mechanism driving the change in the AST on subseasonal scales, however, may be different from that of the PST. Therefore, we will focus only on the PST in this study. This paper is organized as follows. Data employed in the analysis is addressed in Sect. 2. A detailed description of the methodology is given in Sect. 3. Section 4 includes the results of analysis followed by summary and discussion in Sect. 5.

Data
The latest reanalysis data from the ECMWF known as ERA5 (Hersbach et al. 2020) is used for this study. Several improvements can be found in ERA5 compared to its previous generation, ERA-Interim, as detailed in Hersbach and Dee (2016). The ERA5 reanalysis data is produced at a 1-hourly time step with the original horizontal resolution of approximately 30 km and 137 vertical pressure levels. In this study, daily averaged ERA5 data is used trimmed at a horizontal resolution of 1.5° and interpolated at 37 pressure levels during the 39 years of winter (NDJFM) from 1979/1980 to 2017/2018. The resolution and period of the employed data are sufficient to achieve the goals in this study. Variables used include geopotential, air temperature, horizontal wind, vertical (pressure) velocity, and specific humidity at 37 pressure levels. Single level variables include sensible heat flux, latent heat flux from the surface and net longwave radiative flux, net shortwave radiative flux from the surface and top of the atmosphere, which are not assimilated variables but model products. Skin temperatures (SKT) and surface (2 m) air temperatures (SAT) are also used.

Energetics
Variation of storm activities is analyzed based on energetics. Storm associated variables are usually defined as deviations from mean fields or derived by using band-pass filtering. In the present study, deviations from monthly averaged fields are used to define anomalies. This definition is employed in order not to introduce temporal gap between the mean and the eddy fields. Equation governing local energetics is derived based on the quasi-geostrophic approximation (detailed derivation from the primitive equation can be found in the "Appendix"). The resulting equations are fundamentally similar to those in Cai and Mak (1990) or Deng and Mak (2005) except that diabatic forcing term is included. Eddy kinetic energy and potential energy are written respectively as where the stability parameter S is defined by Upper-case symbols denote averages for each month of daily variables and lower-case symbols represent departures except for the independent variables (coordinates). Other symbols follow the convention. Equations for eddy kinetic energy (EKE) and eddy available potential energy (EAPE) can be expressed as where Here, is 2-dimensional mean wind ( = (U, V, 0) ), is 3-dimensional anomalous wind ( = (u, v, ) ), and J is diabatic heating rate. We calculates the diabatic heating rate from the thermal dynamic energy Eq. (23e) as a departure from the adiabatic balance. The ⋅ term in (3) is barotropic energy conversion (BTEC), which represents energy conversion from mean kinetic energy (MKE) to EKE. The baroclinic energy conversion (BCEC) term ⋅ in (4) can be further decomposed into h • h and 3 • 3 , where the subscript h denotes the horizontal components and 3 the vertical component. The former represents energy conversion from mean available potential energy (MAPE) to EAPE, and the latter from EKE to EAPE, which will be referred to as BCEC1 and -BCEC2, respectively. The last term in (4) represents EAPE generation via diabatic heating ( − H ), which is referred to as diabatic energy conversion (DEC) in this study. The advection and divergence of flux terms in (3) (2) and (4) only redistribute eddy energy. We will investigate eddy activities using both EKE and EAPE. After taking time integration, (3) and (4) can be rewritten by Therefore, changes in EKE and EAPE are determined by energy conversion and redistribution during a time interval Δt . Then equation for total eddy energy (TEE) can be obtained by adding (6) and (7).
where, ΔE e = ΔK e + ΔP e . The BTEC, BCEC1, and DEC determine the generation of total eddy energy. Daily variables are used to compute each term. Then, all the terms are monthly averaged for further analysis.

Moist static energy budget
Moist static energy (MSE) budget equation is adopted to investigate the effect of global warming on eddy activities. The time integrated vertical mean MSE budget equation can be written as The MSE, m, is the sum of potential energy, internal energy, and condensational energy: (Φ + c p T + L v q) . SW and LW represent shortwave and longwave radiation trapped in the atmosphere. SH and LH denote sensible and latent heat flux from the surface. The angle brackets denote massweighted vertical integration. We decompose a variable into stationary component (overbar) and transient eddy component (prime) to write the MSE advection term as where the first three terms on the right-hand side represent eddy-induced advection of moist static energy, and the last three terms advection by stationary flow. As in the discussion of energetics, transient components are defined as deviations from monthly averages.

CSEOF analysis
The cyclostationary emperical orthogonal function (CSEOF) technique is used for analysis (Kim et al. 1996;Kim and (6) North 1997). This method has the advantage of isolating the warming signal and analyzing its sub-seasonal variability.
In CSEOF analysis, space-time target variable, T(r, t) , is decomposed as where B n (r, t) are cyclostationary loading vectors (CSLV) and T n (t) are corresponding principle component time series (PC). Unlike EOF loading vectors, CSLVs are time dependent and periodic with period d , which is called the nested period. That is, CSLVs are orthogonal to each other in space-time, and PC are uncorrelated with each other in time. One-to-one correspondences between target variable and other (called predictor) variables are established via regression analysis in CSEOF space (Kim et al. 2015). First, CSEOF analysis on another variable, P(r, t) , yields where C n (r, t) and P n (t) are the CSLV and PC of P(r, t) . Second, multivariate regression on the target PC, T n (t) , is conducted with the predictor PC, P m (t): where (n) m are regression coefficients, (n) (t) is regression error time series, and M (= 40 in this study) is the number of predictor PC used for regression. Then, regressed CSLV can be obtained as As a result of regression analysis, predictor variable can be written as Thus, target and predictor variables can be written together as which ensures that B n (r, t) and C (reg) n (r, t) evolve in a physically consistent manner. This procedure can be repeated for as many predictors as desired. In this study, SAT is used as the target variable and the nested period is

Results and discussion
4.1 Impact of warming on the winter averaged storm activities Figure 1 shows the NH winter warming signal extracted as the second CSEOF mode of SAT. The positive trend in the PC (Fig. 1c) and the positive temperature anomaly over most of the domain (Fig. 1a) indicate that regional warming signal is captured reasonably. The spatial distribution of anomalous warming and weak local cooling is consistent with that in previous studies. The conspicuous warming over the Barents-Kara (BK) Seas and relatively weak cooling over the northern Eurasia are known as the warm-Arctic cold-Eurasia (WACE) pattern (Petoukhov and Semenov 2010;Cohen et al. 2014;Kug et al. 2015;Yeo et al. 2016). The degree of warming varies from one month to another as shown in Fig. 1b. The largest warming of more than 1.5 °C is seen over the Arctic region in January and February while warming in November is relatively weak (~ 1 °C). Warming over the polar region decreases equator-to-pole temperature difference and tends to decrease MAPE. The degree of warming, therefore, is related to the magnitude of MAPE reduction. Although only a small fraction of MAPE is converted into eddy energy, it is the primary energy source for storms. Therefore, it would be of interest to examine whether the reduction of MAPE due to Arctic warming is actually coherent with changes in mid-latitude storm activities. The winter (NDJFM) averaged climatological PST and its change due to warming are analyzed to examine if the results obtained from CSEOF analysis are consistent with previous studies and understand physical processes taking place on sub-seasonal scales. To represent change in storm activities over the entire troposphere, mass-weighted vertical average from 1000 to 200 hPa is taken for EKE, EAPE, and TEE, which will be called VEKE, VEAPE, and VTEE, respectively. Figure 2 shows the climatological PST (shading) and its change due to warming (contours) with respect to (a) VEKE (b) VEAPE, and (c) VTEE. Changes due to warming are obtained from loading vectors regressed onto the SAT warming mode as explained in the method section. The location of the climatological maximum VEKE (45° N, 140° W; Fig. 2a) is further downstream of the maximum VEAPE (40° N, 170° E; Fig. 2b). The impact of warming on the PST can be characterized by changes in latitude and magnitude of maximum storm intensities. In terms of latitude, both the VEKE and VEAPE are shifted northward with positive (negative) anomalies to the north (south) of the climatological maximum. In terms of magnitude, however, EKE is strengthened since the positive anomalies, the location of which overlaps the climatological core of the PST, are more dominant than the negative anomalies (Fig. 2a). On the other hand, EAPE is weakened in low latitudes and is slightly intensified in mid-latitudes. The positive anomalies of total eddy energy ( Fig. 2c; VTEE) are mostly due to the change in VEKE, whereas the negative anomalies are due mainly to VEAPE. The winter-averaged vertical structure of the PST change is shown in Fig. 3 (color contours; primed variables). The zonal average is taken over 150° E-120° W, where change in eddy activities is clearly seen. The climatological PST distribution is shown as gray contours. The vertical structure of EKE change (color contours in the upper panels) due to warming shows that the shift is in the vertical as well as latitudinal directions, since the height of maximum positive (negative) anomaly is higher (lower) than the elevation of the climatological maximum. The poleward and upward shift of the PST is consistent with previous studies (Yin 2005;Wu et al. 2010;Chang et al. 2012;Woollings and Blackburn 2012). Change in EKE corresponds reasonably with the sum of BCEC2 ′ and BTEC ′ (Fig. 3c), except north of 60° N and in the middle troposphere near 40° N. The positive BTEC ′ in Fig. 3b is responsible for the upper-tropospheric EKE generation whereas the negative BCEC2 ′ in Fig. 3a is responsible for the lower-tropospheric EKE suppression. Considering that barotropic damping is more prominent than generation when seasonal average is taken (negative BTEC; see, for example, Chang et al. 2002 andLee et al. 2011), the positive BTEC ′ in Fig. 3b indicates a weakening of barotropic damping. Wang et al. (2017) also estimated BTEC and shows that its change favors a southward shift of the PST, which is opposite to the observed PST change. They, for that reason, concluded that the interdecadal change of the PST is not likely the result of BTEC change. This difference from our results seems to have stemmed from the difference in the definition of eddies as 2-8 day band-pass filtered anomalies. Despite that, several similar features are seen including significant positive BTEC anomalies at 30° N and 50° N and negative anomaly north of 60° N (see Fig. 6c in Wang et al. 2017).
Change in EAPE (color contours in the lower panel) is found in relatively low troposphere with a poleward shift. The positive anomaly to the north (Fig. 3f) is due to the enhanced energy supply from MAPE (positive BCEC1 ′ in Fig. 3d) and weakened energy conversion into EKE (positive -BCEC2 ′ in Fig. 3e). On the other hand, the negative anomaly to the south is due to the weakened energy supply from MAPE (negative BCEC1 ′ in Fig. 3d), which is partly cancelled by the weakened energy conversion into EKE (Fig. 3e). This result is consistent with Wang et al. (2017).
Advection, geopotential flux, and diabatic heating in Eqs.
(3) and (4) are also evaluated in Fig. S1. In the climatological average, the magnitudes of the energy redistribution terms are comparable to BCEC (Chang et al. 2002). The positive EKE anomaly north of 60° N can be understood in terms of the convergence of geopotential flux term (Fig. S1b), which is positive (negative) north (south) of 60° N. This means that storms tend to move northward more actively in a warming climate (Fig. S1b). The negative Fig. 2 The winter (NDJFM) averaged climatological storm intensity (shading; m 2 s -2 ) and anomalous storm intensity (contour) for the warming mode as represented by a VEKE, b VEAPE, and c VTEE. Contour interval for a and c is 2 m 2 s -2 from ± 2 m 2 s -2 with blue contours for negative values and red contours for positive values. Contour interval for b is 0.5 m 2 s -2 from ± 0.5 m 2 s -2 anomalous energy conversion in the middle troposphere near 40° N (Fig. 3c) is only partially offset by the positive EKE advection and geopotential flux. The DEC shows relatively large positive values near the surface at high latitudes (Fig.  S1d), which contributes to an increase in EAPE. However, contribution of the advection terms is relatively small (Fig.  S1a, and c).

Seasonal variation of the impact of warming on storm activities
As NH warming shows strong subseasonal variability (Fig. 1b), regional warming pattern near the Pacific jet entrance region also exhibits substantial seasonal variation (Fig. 4). It can be seen that cooling extends from Eurasia to the east of Japan from November to January (Fig. 4a, c), but warming is located from February to March (Fig. 4b,  c). This early winter cooling/late winter warming is seen throughout the troposphere, and is accompanied by strong seasonal variation in turbulent heat flux. Considering the strong association of jet and storms, temperature variation in the western North Pacific is expected to result in significant seasonal variation in the impacts of warming on the PST. Therefore, seasonal change in the PST within the NH winter season is investigated. Figure 5 shows the zonal mean (150° E-120° W) VEKE and VEAPE for the warming mode averaged during NDJ (blue), FM (red), and NDJFM (black). Patterns for individual months are also presented (color dashed lines). The vertical dashed lines denote the latitudes of climatological maximum VEKE and VEAPE. Consistent with earlier results, winteraveraged VEKE (black line in Fig. 5a) exhibits both poleward shift and intensification of storm activities. However, a comparison between the two stages shows that the effect of warming is not uniform through the winter. VEKE shifts poleward in the early stage (blue solid line), whereas it is intensified in the late stage (red solid line); a similar trend is also seen in monthly patterns. This seasonal difference between the two Fig. 3 The upper panels represent zonally averaged pattern (150°E-120°W) of a BCEC2 ′ , b BTEC′, and c summation of both associated with the warming mode (shading; 10 -5 W m -3 ). The color contours (blue < 0 < red, at the interval of 1 m 2 s -2 ) represent EKE for the warming mode in comparison with the climatological EKE (grey contours at the interval of 20 m 2 s -2 ). The shading in the lower panels are d BCEC1 ′ , e −BCEC2 � , and f summation of both over the same domain as the upper panels. Color contours in the lower panels represent EAPE (at the interval of 0.5 m 2 s -2 ) for the warming mode and the climatological EAPE (grey contours at the interval of 10 m 2 s -2 ). Note that the energy conversion terms in Eqs. (6) and (7) are scaled by the mean density (detailed scaling factor can be found in the "Appendix") periods is statistically significant as can also be seen in the linear trend of VTEE (see Fig. S2).  also showed that strorm track intensity in the Pacific has enhanced in early spring (particularly in February and March) during the recent decades. Detailed analysis on sub-seasonal scales, however, was not performed in earlier studies. Similar changes are found in VEAPE (Fig. 5b) except that the magnitude of the negative anomalies is larger than the positive anomalies during the early winter, implying an overall reduction in VEAPE. It is due to a significant reduction of EAPE during November in low-and high-latitude regions.
The vertical structure of EKE during the two stages are analyzed in the framework of energetics (Fig. 6). Shading in the upper (lower) panel shows change in BCEC2 (BTEC) due to warming averaged during the early stage (left column) and the late stage (right column) of winter. Zonal mean is taken over the longitude band 150°E-120°W. In the early winter, negative BCEC2 ′ (Fig. 6a) depresses EKE particularly in the middle troposphere near 40° N. BTEC ′ increases high-latitude EKE but reduces low-latitude EKE (Fig. 6c). In the late winter, intensification of EKE is consistent with the positive BTEC ′ (Fig. 6d). Therefore, BTEC ′ plays a more crucial role in changing EKE, while BCEC2 ′ plays a secondary but non-negligible role.
The vertical structure of EAPE is analyzed in Fig. 7 in a similar manner. As mentioned before in Fig. 5, the negative anomaly of EAPE in the early winter is more pronounced than the positive anomaly. Weakening of energy conversion from MAPE to EAPE is a main cause for the reduced eddy potential energy in low latitudes (Fig. 7a). Weakened energy conversion into EKE (positive -BCEC2 ′ in Fig. 7c) partly offsets the negative BCEC1 ′ in lower latitudes, but produces positive EAPE anomaly in higher latitude. In the late winter, EAPE is enhanced particularly over the latitude of the climatological maximum PST (~ 40° N). Both the energy conversion terms contribute to an enhancement of eddy potential energy with a stronger inflow of energy from MAPE and a weaker outflow to EKE. In the early stage of winter, when both EAPE and EKE are shifted to the north, mass averaged total EAPE is weakened whereas total EKE is strengthened (see Figs. 6a, c and 7a, c). Therefore, it can be said that eddy energy resides more in the form of kinetic energy and less in the form of potential energy as warming proceeds. In the late winter, on the other hand, both EKE and EAPE are enhanced (see Figs. 6b, d and 7b, d).
Total eddy energy production is primarily determined by BCEC1 and BTEC. In the early winter, negative BCEC1 ′ is largely responsible for the negative TEE′ (blue dashed contours in Fig. 8a, c) located near the southern flank of the climatological PST (Fig. 8a). On the other hand, positive TEE′ (red solid contours) located to the north is mainly due to BTEC ′ (Fig. 8c). In the late winter, both BCEC1 ′ and BTEC ′ contribute to an enhancement of TEE′. As can be seen, changes in eddy kinetic and potential energy on subseasonal scales are reasonably explained in terms of changes in energy conversion terms. Now, the question is what leads to change in the two energy conversion terms. To investigate this, component analysis of the two conversion terms is performed.

Component analysis of energy conversion terms
Component analysis of BCEC1 ′ is conducted in Fig. 9. BCEC1 is determined by multiplying eddy heat flux and temperature gradient ( h • h ). Both components are important since enhanced temperature gradient implies a more efficient energy conversion for given heat flux; heat flux Zonally averaged (150°-120° W) pattern (shading; 10 -5 W m -3 ) of BCEC2 ′ (upper panels) and BTEC ′ (lower panels) during NDJ (left column) and FM (right column) together with EKE (blue < 0 < red contours at the interval of 1 m 2 s -2 ) for the warming mode alone can not induce any energy conversion if there is no temperature gradient. Further, they are also dependent on each other; heat flux tends to relieve temperature gradient. The BCEC1 change ( h • h ) ′ due to warming can be further decomposed into The overbar and prime in (18a, 18b) denote climatological average and change due to warming, respectively. We first calculate BCEC1 from the raw data and then calculate its change due to warming as in (18a). Then the result is h due to warming are used together with their climatological averages for the estimation of BCEC1 ′ leaving out the nonlinear term (last term in 18b). The nonlinear term turns out to play only a minor role and the linear components explain majority of BCEC1 ′ . Besides, the terms associated with the change in meridional temperature gradient 2 • � 2 and poleward heat flux � 2 • 2 in (18b) turn out to be more important than their zonal counterparts ( 1 • � 1 and � 1 • 1 ). All the components needed for calculating BCEC1 ′ are shown in Fig. S3 and S4. Figure 9 shows the components of BCEC1 ′ associated with the change in meridional temperature gradient ( 2 • � 2 ; a and d) and heat flux ( � 2 • 2 ; b and e), together with the pattern of BCEC1 for the warming mode ( h • h � ; c and f). The anomalous temperature pattern Fig. 7 The zonally averaged (150°-120° W) pattern (shading; 10 -5 W m -3 ) of BCEC1 ′ (upper panels) and −BCEC2 � (lower panels) during NDJ (left column) and FM (right column) together with EAPE (blue < 0 < red contours at the interval of 0.5 m 2 s -2 ) for the warming mode is presented as contours in Fig. 9a, d. All these terms are calculated at 400-hPa level, where maximum BCEC1 ′ is observed. The horizontal distribution shows that energy conversion from MAPE to EAPE has generally weakened during NDJ (Fig. 9c) but has shifted to the north during FM (Fig. 9f), consistent with their vertical structures (Fig. 7a,  b). The component analysis reveals that the term associated with the changes in heat flux (Fig. 9b, e) is much larger and qualitatively more similar to BCEC1 ′ than that associated with changes in temperature gradient (Fig. 9a, d) in both the stages. Similar results are obtained in other levels in the troposphere as can be inferred from their vertical structures (Fig. 7a, b). Since the mean meridional temperature gradient and the stability parameter, ( ∕ p) −1 , are always negative, positive (negative) � 2 • 2 indicates poleward (equatorward) eddy heat flux. For example, negative � 2 • 2 centered over the central Pacific (170° E, 37° N) in NDJ (Fig. 9b) is accompanied by equatorward eddy heat flux. Divergence of heat flux is expected along ~ 40° N in FM. A similar change in eddy heat flux can also be found in other levels of the troposphere except near the surface.
Cooling in the early winter and warming in the late winter over the western North Pacific produce anomalous temperature gradient of opposite signs (Fig. 9a, d). In the early winter, temperature gradient south of 40° N becomes stronger, while it becomes weaker to the north, resulting in change in background baroclinicity. In the late winter, sign of anomalous temperature gradient is reversed. As can be seen in Fig. 8 The zonally averaged (150°-120° W) pattern (shading; 10 -5 W m -3 ) of BCEC1 ′ (upper panels) and BTEC (lower panels) during NDJ (left column) and FM (right column) together with TEE′ (blue < 0 < red contours at the interval of 1 m 2 s -2 ) for the warming mode Fig. 9c, ′ 2 is not significantly correlated with BCEC1 ′ in the early winter, while ′ 2 match reasonably with BCEC1 ′ in the late winter (Fig. 9f). This indicates that change in heat flux may not necessarily be associated with the change in baroclinicity. We will address this issue in more detail in the next section.
Change in BTEC is analyzed in a similar manner to BCEC1 ′ . The anomalous BTEC change, ( • ) ′ , due to warming can be further decomposed into Similar to BCEC1 ′ , we first calculate (19a) then compare with (19b) leaving out the nonlinear term ( � • � ). The x component of BTEC ′ is associated with stretching deformation ( D 1 ) of the mean field and the anisotropy of eddies ( E 1 ), and the y component shear deformation ( D 2 ) of the mean field and the axial tilting of eddies ( E 2 ) (see Mak and Cai 1989;Black and Dole 2000 for details). Climatologically, BTEC over the western to central Pacific is dominated by the term associated with shearing deformation (Lee et al. 2010).
In the warming mode, BTEC ′ ( ( • ) � ; Fig. 10c, f) is strongly correlated with the change in E 2 under shearing deformation ( E ′ 2 D 2 ; Fig. 10a, d) with the pattern correlation of 0.61 (0.68) in early (late) winter. The sum of the other linear components ( significant correlation with ( • ) � . The negative E ′ 2 D 2 in the early winter (Fig. 10a) is due to positive E ′ 2 over the region of negative shear deformation (northern flank of the jetstream; see Fig. 11c). In the late winter (see Fig. 10d), on the other hand, positive E ′ 2 D 2 is seen over the region of negative shear (northern flank of the jet) and negative E ′ 2 D 2 over the region of positive shear (southern flank of the jet), indicating that E ′ 2 is negative (Fig. 11f). While linear analysis highlights the axial tilting of eddy as a single most important parameter in the western North Pacific, BTEC change is not fully explained in the eastern part of the North Pacific (Fig. 10c, f) without the nonlinear term. In the context of linear analysis employed in the present study, it is difficult to account for the nonlinear effect in an accurate manner.
Change in E 2 anomaly arises due to change in eddy shape as well as eddy strength (i.e. EKE ). Thus, we define normalized E 2 by EKE as in (20).
The Ẽ 2 vector provides information about the eddy shape independently of the magnitude of EKE. Anomalous E 2 vector, then, can be rewritten as Change in E ′ 2 is decomposed into the terms associated with change in eddy shape ( Ẽ � 2 × − EKE ), change in EKE ( Ẽ 2 × EKE � ), and nonlinear term ( Ẽ � 2 × EKE � ). Analysis (20) Fig. 9 The 400-hPa BCEC1 ′ component (10 -4 W kg -1 ) due to (a and d) anomalous meridional temperature gradient ( h • � h ) and (b and e) anomalous meridional heat flux ( � h • h ). c, f Represent the regressed pattern of BCEC1 ′ for the warming mode. The upper panels represent the NDJ averages and the lower panel the FM averages. The color contours in a and d represent the regressed temperature pattern at the same vertical level (blue < 0 < red contours at the interval of 0.25 K). Details of calculation can be found in the main text (see Eq. 18) shows that change in EKE (Fig. 11b, e) and the nonlinear term (not shown) play a minor role and change in eddy shape (Fig. 11a, d) is a dominant factor. The BTEC ′ in Fig. 10c, f is largely owing to change in the eddy shape associated with the axial tilting of eddies ( Ẽ ′ 2 ).

Moist static energy budget
BCEC is an important component of the seasonal change in the PST under regional warming. The primary factor determining change in BCEC is shown to be eddy heat flux. In terms of time-averaged atmospheric energy budget, Fig. 10 The 350-hPa BTEC ′ component (10 -4 W kg -1 ) due to a and d axial tilting of eddy shape ( E ′ 2 D 2 ) and b and e summation of other . c, f Represent the regressed pattern of BTEC ′ for the warming mode. The upper panels represent the NDJ averages and the lower panels the FM averages. Details of calculation can be found in the main text (see Eqs. 19a,19b)

Fig. 11
The components of 350-hPa E ′ 2 (m 2 s -2 ) associated with change in eddy shape ( Ẽ � 2 × EKE ; a and d and EKE ( Ẽ 2 × EKE � ; b and e). c, f Represent the regressed pattern of E ′ 2 for the warming mode. The upper panels represent the NDJ averages and the lower panels the FM averages convergence and divergence of heat flux produced by largescale circulation and storms partly compensate heating and cooling at the boundary. Understanding how radiation or heat flux changed at the boundary due to warming and how such a change is related to heat flux will provide useful framework for diagnosing changes in the PST activities.
The MSE budget equation is adopted for this purpose, which considers not only internal heat energy and potential energy but also conversion of latent heat energy.
In the climatological MSE budget, Δ⟨m⟩ in (9) should be close to zero so that MSE advection and boundary heating terms are nearly in balance. In other words, cold advection is expected in the region of heating, and vice versa. Figure 12 shows the winter averaged climatological boundary heating (LW + SW + SH + LH) and vertically integrated advection of MSE. Much of the atmospheric heating (Fig. 12a) is confined to the adjacency of the Kuroshio Current (KC), with relatively weak cooling over the eastern Pacific and the polar region. The compensating intense negative MSE advection (Fig. 12b) along the KC and the positive MSE advection over the eastern Pacific and polar region illustrate how energy balance is achieved in the time-averaged atmosphere. The decomposition of the MSE advection into the (c) eddy and (d) stationary components shows that they are of opposite signs and tend to cancel each other. To the south of the KC axis, the negative MSE advection in Fig. 12b is caused by the transient component (Fig. 12c), partially cancelled by the stationary component (Fig. 12d). The situation is reversed to the north of the KC axis; the negative MSE advection by the stationary component is cancelled by the transient component of opposite sign. The vertically integrated MSE advection can also be interpreted as the convergence of vertically integrated horizontal MSE flux with the aid of mass conservation equation. Assuming that vertical velocity is zero at the top and bottom of the atmosphere, we have where the angle brackets denote vertical averaging. The zonally elongated convergence zone of eddy-induced MSE flux (Fig. 12c) to the north of 40° N and divergence zone to the south indicate a maximum poleward eddy heat flux near 40° N. A similar interpretation with opposite signs can be applied to the stationary component (Fig. 12d). Northward flux near 20° N, southward flux near 35° N, and northward flux near 70° N indicate the three-cell structures of largescale circulation during the boreal winter. The latitudinal structures of the MSE flux for the transient and stationary components are consistent with those in Michaud and Derome (1991).
Changes in the MSE, boundary heating, and MSE advection (MSE flux convergence) due to warming are shown in Fig. 13. Contours in (a and e) and (b and f) denote vertically integrated MSE and meridional MSE flux by the transient eddy component, respectively. Changes in the MSE advection by the stationary and all the components are also presented in Fig. S5 and S6. During the early stage of winter (left column), anomalous heating is elongated from southeast of Japan to the eastern North Pacific with a weak cooling over the northern and southern North Pacific (Fig. 13a). In the late winter, situation is reversed with anomalous heating over the Sea of Okhotsk and the southern North Pacific and cooling between them (Fig. 13e). As shown in Fig. S5 (early winter) and Fig. S6 (late winter), MSE advection roughly compensates the anomalous heating at the surface, and the residual determines the MSE change (contours in a and e). Although the patterns of MSE advection ( Fig. S5b and S6b) are somewhat noisy, they are roughly opposite in sign to the boundary heating particularly over the Sea of Okhotsk, the southeastern Japan, and the southern North Pacific in both stages of winter. Decomposition of the MSE advection reveals that the transient component over the western North Pacific (Fig. 13b, f) is positively correlated with the boundary heating, a significant fraction of which is offset by the stationary component ( Fig. S5d and S6d). MSE advection by the stationary component produces cold and dry air to the southeast of Japan and warm and moist air over the northern by c p (blue < 0 < red contours at the interval of 1 K m s -1 ). c and g Skin temperature minus 2 m air temperature (K) and d and h) saturation specific humidity at surface minus 1000-hPa specific humidity (g kg -1 ) for the warming mode. The left (right) column represents the NDJ (FM) Pacific, respectively. This change is partially compensated by advection by the transient component and boundary heating. As a result, negative eddy MSE flux is seen between Japan and the northern North Pacific (Fig. 13b). It is noted that flux of internal energy ( c p T ) dominates explaining 84% of the negative eddy MSE flux, so that a strong resemblance is seen between eddy MSE flux (contours in Fig. 13b) and heat flux (Fig. 9b). Boundary heating and anomalous heat flux by transient eddies, which largely explains the weakened BCEC1 ′ in the early winter (Fig. 9b), compensates the MSE advection by stationary component, although the causal relationship among them cannot be ascertained. A similar reasoning applies to the late winter with enhanced poleward heat flux and increased BCEC1 ′ over the western North Pacific. Since boundary forcing differs between the two periods, eddy heat flux is also expected to vary between the two periods. Turbulent heat flux (SH and LH) is found to be more important than radiative flux in the change of boundary heating in both stages of winter ( Fig. S7 and S8). Heat flux from the surface can be estimated from temperature and specific humidity differences between the surface and the air above by using the bulk formula. Figure 13 shows SKT minus SAT (c and g) and saturation specific humidity at surface obtained from the SKT based on the Clausius-Clapeyron equation minus 1000-hPa specific humidity (d and h) averaged during the early and late stages of winter. Although both SKT and SAT increase in a similar manner during winter, their difference shows distinct spatial patterns. As expected, the temperature difference (Fig. 13c, g) matches well with the sensible heat flux ( Fig. S7a and S8a). The humidity difference (Fig. 13d, h) also matches with the latent heat flux ( Fig. S7b and S8b). In both stages of winter, change in mass exchange and, to a lesser extent, heat exchange constitutes an important source of boundary heating, which is closely related to sub-seasonal variability of eddy-induced heat flux.

Summary and concluding remarks
The effect of warming on the sub-seasonal variation of the PST during the NH winter season is investigated. CSEOF analysis is used to capture the sub-seasonal variability associated with regional warming. Results in the present study are consistent with earlier studies in that the location of strong eddy activities and energy conversion is shifted northward when viewed in a winter (NDJFM) average sense. On a seasonal scale, however, the PST is shifted northward in the early winter (NDJ), while the PST is enhanced in the late winter (FM). In the early winter, energy conversion from MAPE to EAPE is weakened near the southern flank of the climatological PST. On the other hand, barotropic damping is weakened near the northern flank, resulting in a general northward shift in the eddy activities in the Pacific. In the late winter, enhanced baroclinic energy conversion as well as weakened barotropic damping contribute to the amplification of the PST. Changes in the eddy energy redistribution terms, advection and convergence of geopotential flux, are also evaluated along with the diabatic contribution. Convergence of geopotential flux contributes most among them particularly in the upper troposphere. Their magnitudes, however, are small, emphasizing the importance of baroclinic and barotropic energy conversions.
Component analysis is adopted to identify important parameters causing changes in energy conversions under regional warming. Change in energy conversion from MAPE to EAPE ( BCEC1 ′ ) is largely driven by transient eddy anomalous meridional heat flux due to warming. The anomalous eddy heat flux does not seem to be directly related to changing baroclinicity. From the MSE budget perspective, MSE advection by transient eddies is positively correlated with boundary heating but is negatively correlated with MSE advection by the stationary component in the western Pacific, where change in BCEC1 ′ is most significant. Latent heat flux constitutes the most important source for the change in boundary heating followed by sensible heat flux. Contribution of radiative flux is relatively small.
Unlike many other studies, which focused on the baroclinic process, our study, based on the quantitative analysis, also highlights the role of barotropic process in explaining the change in storm activities under warming. Exchange of kinetic energy between eddies and the jetstream (BTEC) has changed due to the axial tilting of eddies rather than the anisotropy of eddy shape or the deformation of the jetstream. While the change in axial tilting is the single most important factor for BTEC change, nonlinear interaction between anomalies of the mean field and those of eddies is not negligible particularly in the eastern part of the North Pacific. Further, it is not clear how the axial tilting of eddies changes under warming; this mechanism is yet to be answered in future studies. Although our analysis concentrates on the Pacific domain, change in storm track activities over the Atlantic may also be a topic of interest, considering the distinct characterstics of the AST.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.
Equations governing the motion are approximated as where, u, v, , , and represent zonal velocity, meridional velocity, vertical (pressure) velocity, geopotential and potential temperature, respectively, and J is diabatic heating (W kg -1 ). Let us decompose physical variables into two parts where the variables in upper case represent the long-term means and the primed variables denote anomalies from the respective means. The long-term average fields are assumed to satisfy the geostrophic and hydrostatic relationship. That is, The hydrostatic equation can be rewritten as Then, departures from the geostrophic balance satisfy We will drop the primes in the equations below, noting that the variables in lower case represent anomalies from the long-term means (upper case variables). Multiplying the zonal component of the equation of motion by u , we have which can be rewritten as In a similar manner, the meridional component of the equation motion is multiplied by v to yield Adding the two equations, we have t 1 2 where K = u 2 + v 2 ∕2 is the eddy kinetic energy, = (U, V, 0) , and = (u, v, ) . The eddy kinetic energy can be rewritten as where the subscript h denote the horizontal components. Since the geostrophic mean field is non-divergent, the second term on the right hand vanishes, so that we can write the eddy kinetic energy as Eddy available potential energy is defined as where the stability parameter S is defined by By multiplying the thermodynamic equation, we obtain By multiplying the equation above by −R( dΘ dp ) −1

, we have
The eddy kinetic energy equation and the eddy available potential energy equation can be rewritten in a succinct manner as (34) (36) P = 1 2S p By using the hydrostatic equation, we can rewrite the equations above as (52)