**1. Introduction**

The Northern Hemisphere stratospheric polar vortex (SPV) is characterized by a significant annual evolution or annual cycle [1]. It spins up from autumn and persists through winter into spring, at which point it breaks up and the typical wintertime subpolar westerlies are replaced by easterly summer circulation until the subsequent autumn [2,3]. The transition from a winter to a summer circulation regime is known as the stratospheric final warming (SFW) [4], which is influenced by both solar radiation and dynamics [5].

There are various indices commonly used to describe the strength of the stratospheric polar vortex. These indices include polar cap (60–90◦ N) height [6], polar cap temperature, potential vorticity on isentropic surfaces [7–9], polar mass [10], zonal-mean zonal wind at subpolar latitudes [11–13] and Northern Annular Mode (NAM) [14–17]. Baldwin and Thompson [18] have made a comprehensive comparison among those indices and found that these indices generally correspond well with each other and are capable of representing the variability in the SPV. The anomalous changes in these indicators always exhibit remarkably larger amplitudes from late fall to early spring each year and are dominated by intra-seasonal variabilities [14,15,19,20].

**Citation:** Yu, Y.; Sun, J.; Secor, M.; Cai, M.; Luo, X. A Parametric Model of Elliptic Orbits for Annual Evolutions of Northern Hemisphere Stratospheric Polar Vortex and Their Interannual Variability. *Atmosphere* **2023**, *14*, 870. https://doi.org/

Academic Editor: Yoshihiro Tomikawa

10.3390/atmos14050870

Received: 12 April 2023 Revised: 11 May 2023 Accepted: 12 May 2023 Published: 14 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The dynamic influence of the stratospheric polar vortex during the winter months is known to play a part in shaping the winter circulation patterns of the troposphere such as Arctic Oscillation and blockings [17,20–26]. As a result, the anomalous signals in the stratospheric polar vortex can exert important downward impacts on weather and climate in the extratropical troposphere. For instance, during the 1–2 months after a weaker stratospheric polar vortex event such as sudden stratospheric warming (SSW) and negative stratospheric northern annular mode event, below-normal temperatures tend to occur more frequently over the Northern Hemispheric continents [24,27–38]. In addition, the polar temperature related to the SPV in the late winter months is closely related to the formation of polar stratospheric clouds and thus affects the Arctic ozone depletion [39–42]. The interannual variations in seasonal timing and the features of stratospheric final warming have been found to affect the spring circulation anomalies, which in turn modifies the summer monsoon [4,43], and affect precipitation anomalies over South Asia in May and the temperature anomalies over Central Asia in March [44]. Later, stratospheric final warming and stronger SPV in springtime have also been found to lead to negative ice thickness anomalies in the East Siberian Sea [3]. The strength of SPV can be the sub-seasonal to seasonal forewarning of anomalous atmospheric river frequency [45]. Therefore, identifying and understanding the year-to-year varying annual cycle of SPV has potential implications for winter seasonal forecasts, as the December–February mean behavior may miss important sub-seasonal events.

Hardiman et al. [46] formulate a simple sine wave fit to the observed vacillations in the daily zonal, mean zonal wind anomalies at 10 hPa, which can be used to explain much of the sub-seasonal and interannual variability in the monthly mean vortex strength. Additionally, consistent with wave-mean flow interaction theory, their amplitude correlates positively with the magnitude of winter mean planetary wave driving. They then tested the predictive skill of this simple sine wave model throughout the winter against the persistence of the vortex strength and against a state-of-the-art seasonal forecasting system, and reported improved prediction skills of the vortex strength one month ahead.

Not only the zonal wind surrounding the SPV but also the thermal condition in the stratospheric polar region are important indicators of the SPV intensity changes. The changes in momentum and thermal fields associated with SPV variations are intimately coupled as expected from the thermal wind relationship. According to Zuev and Savelieva [47], the temperature distribution inside the SPV at a specific range of wind speed was approximately in the same temperature range, but the temperature clearly shows year-to-year differences. This study formulates a parametric ellipse orbit model for better capturing the coupled annual evolutions of both the zonal momentum of the polar jet and the total air mass in the polar cap of the stratosphere and their interannual variability.

#### **2. Data and Methods**

#### *2.1. Data*

Daily mean fields are derived from the hourly ERA5 reanalysis dataset [48], which is on 1.5◦ × 1.5◦ grids and covers 42 years (July of the previous year to June of the current year) from 1979/80 to 2020/21. We use variables including the three-dimensional temperature, geopotential height and wind fields at 37 pressure levels from 1000 hPa to 1 hPa, air pressure, 2 m air temperature and wind fields at the surface. A 7-day running mean is applied to minimize synoptic-scale perturbations.

#### *2.2. Methods*

In this study, we use the stratospheric polar mass (*M*) above the isentropic surface of 400 K inside the Arctic polar circle and mass-weighted zonally integrated zonal momentum at 60◦ N (*MU*) to represent the polar vortex and polar jet changes. The isentropic level 400 K is about 100 hPa in the subpolar latitudes. This level is chosen in our study because the stratospheric circulation signals above 400 K can actively interact with the upward propagating tropospheric waves [49–56], which in turn exert significant downward impacts on the tropospheric circulation and the weather in winter [19,57,58]. The latitude of 60◦ N is chosen to better capture the edge of the stratospheric polar vortex and the maximum of the westerly surrounding the vortex in the winter seasons. These two variables are addable and closely linked to the mass budget and angular momentum budget [10,26,59], and thus can be quantitatively attributed to the adiabatic processes mainly driven by wave dynamics and diabatic processes.

Following Pauluis et al. [60,61] and Yu et al. [10,62], we first interpolated the daily potential temperature and wind fields onto 200 equally spaced sigma (*σ*) levels from 1 to 0. The air mass between two adjacent sigma surfaces per unit area is m*<sup>σ</sup>* = Δ*σ*/g Ps, where g is the gravitational constant, Δ*σ* = 1/200 and Ps is the surface pressure. According to Yu et al. [10], the *M* can be measured by the integrated mass above 400 K isentropic level north of 60◦ N (kg), which is derived as

$$M(t) = \int\_{\pi/3}^{\pi/2} \int\_0^{2\pi} \int\_0^1 \mathbf{m}\_{\sigma}(\lambda, \ \phi, \sigma, t) \cdot H(\theta(\lambda, \ \phi, \sigma, t), 400 \text{ K}) \mathrm{d}\sigma \mathrm{Rcos}\phi \mathrm{d}\lambda \mathrm{Rd}\phi \tag{1}$$

where *λ* is the longitude and *φ* is the latitude, *R* is the Earth radius; *H* is the Heaviside function, namely *H*(*x*, *x*0) = 1 if *x* ≥ *x*0, and otherwise *H*(*x*,*x*0) = 0. In (1) and (2) below, *x = θ(λ,ϕ,σ,t)* and *x*<sup>0</sup> = 400 K.

Similarly, the intensity of the polar jet can be measured by the zonally integrated zonal momentum in the stratospheric layer above 400 K across 60◦ N (*MU*, unit: kg m s<sup>−</sup>1), which can be obtained according to

$$\text{MLI}(t) = \int\_{0}^{2\pi} \int\_{0}^{1} \mathbf{m}\_{\sigma}(\lambda, \boldsymbol{\phi} = 60 \,\text{N}, \sigma, t) \cdot \mathbf{u}(\lambda, \boldsymbol{\phi} = 60 \,\text{N}, \sigma, t) \cdot H(\boldsymbol{\theta}(\lambda, \boldsymbol{\phi} = 60 \,\text{N}, \sigma, t), 400 \,\text{K}) \mathrm{d}\sigma \text{Rcos}\boldsymbol{\phi} d\lambda \tag{2}$$

Smaller values of *M* indicate a stronger stratospheric polar vortex, while larger values of *MU* indicate a stronger polar jet and vice versa.

### **3. Results**

#### *3.1. Statistics of the SPV Indices*

From the annual evolutions (Figure 1), we see that the stratospheric mass in the polar region (*M*) and the relative zonal momentum around the polar jet (*MU*) vary generally out of phase. Namely, the *MU* increases while the *M* decreases before February and then the *MU* decreases while the *M* increases afterward. This illustrates that the Northern Hemisphere SPV and the polar jet surrounding it persist in the winter months and peak around mid-January. Radiative heating/cooling processes provide overall control of the annual cycle of *M* thermodynamically and adiabatic eddy-driven poleward mass transport processes affect the amplitude and phase of the annual cycle dynamically. The annual cycle of *MU*, besides through the thermal wind balance which is a fast process (less than a day), is controlled by planetary-scale waves that drive the poleward mass transport in the first place. Specifically, a stronger (weaker) poleward transport tends to coincide with the weakening (strengthening) of *MU*, which is about a few days prior to the increase in *M* [10]. In other words, at seasonal scales, the annual cycle of *MU* is driven by the radiative processes via thermal wind balance. However, at sub-seasonal time scales, particularly in the winter seasons, temporal evolutions of *MU* are driven by planetary-scale waves, which in turn drive the temporal evolution of *M* through eddy-driven poleward mass transport. The standard deviations of *MU* and *M* become remarkably larger in the months from December to March compared to the other months throughout the year, consistent with previous studies [14,15,19,20]. Note that the *MU* shows larger interannual variations around its climatology than the *M*. In particular, the maximum standard deviation of *MU* during the winter months is around 60% of its climatological mean, while for the *M*, its maximum standard deviation is only about 10% of its climatological mean. The seasonal changes of the coupled *M* and *MU* also exhibit remarkable interannual variation, as can be seen clearly from Figure 2. The timing of the strongest SPV, which can be represented by

the minimum *M* and maximum *MU*, ranges from December to April, during the period from 1980 to 2020.

**Figure 1.** The 42-year climatological annual evolutions (red curves) and their interannual variability (blue shadings) in *M,* (**a**) units: 1016 kg and *MU* (**b**) units: 10<sup>12</sup> kg s<sup>−</sup>1).

**Figure 2.** Annual evolutions (the ordinate) of (**a**) *M* (units: 1016 kg) and (**b**) *MU* (units: 1012 kg s−1) in the years from 1980 to 2020 (abscissa).

A closer look at the temporal evolutions of the *M* and *MU* in a few selected years (Figure 3) reveals that *MU* leads *M* by a few days, even though they are generally negatively correlated. We then investigate the lead–lag correlations of *M* and *MU* in each year and the lead days of the *MU* relative to the *M* when their maximum negative correlations are reached (Figure 4). It is seen that the maximum negative correlations occur mostly when the *MU* leads the *M*. The lead time is mainly in the range of 1–10 days with a climatological mean value of about a week. As discussed above, it is the planetary-scale wave activity that drives the temporal evolutions of *MU*, which in turn leads to changes in the *M* field with a lead-time of a few days [10]. Therefore, temporal evolutions of *M* at sub-seasonal time scales are strongly modulated by eddy-driven poleward energy transport.

In summary, the observed annual evolutions of *M* and *MU*, respectively, representing the changes in the thermal and dynamical conditions of the SPV throughout the year, have three main features: (i) the *M* and *MU* are dominated by the annual cycle with the largest interannual variability in winter; (ii) the *MU* has a larger interannual variability than the *M*; and (iii) the *MU* varies out of phase with the *M* but with *MU* leads to *M* by 1–10 days.

### *3.2. The Parametric Elliptic Orbit Model*

In this section, we will construct a parametric model of elliptic orbits to capture the coupled changes of the *M* and *MU*. To begin with, because the annual cycle of the *M* or the *MU* is close to a sine or cosine function, we can first formulate two variables, *X*, and *Y*, as

$$\begin{cases} X = X\_0 + a\cos\Theta\cos\alpha + b\sin\Theta\sin\alpha\\ Y = Y\_0 + a\cos\Theta\sin\alpha - b\sin\Theta\cos\alpha \end{cases} \tag{3}$$

where *X*<sup>0</sup> and *Y*<sup>0</sup> are the coordinates of the elliptic central point; *a* and *b* are the amplitudes; Θ is the phase angle in the range from 0 to 2*π*; and *a* corresponds to the tilting angle of

the elliptic obit with the *x*-axis ranging from −*π* to *π*. The other form of (3) is the familiar elliptic equation, namely

2

*<sup>b</sup>*<sup>2</sup> <sup>=</sup> 1 (4)

(*X* − *X*0) 2 *<sup>a</sup>*<sup>2</sup> <sup>+</sup> (*<sup>Y</sup>* <sup>−</sup> *<sup>Y</sup>*0)

**Figure 3.** Time series of *M* (red curves, units: 1016 kg) and *MU* (blue curves, units: 10<sup>12</sup> kg s−1) in four example years: (**a**) 1980, (**b**) 1990, (**c**) 2000, (**d**) 2010. The solid brown curves and light blue curves are 7-day running means of *M* (ordinate on the right-hand-side) and *MU* (ordinate on the left-hand-side), respectively; and dashed curves are original time series of *M* and *MU*, respectively. The solid magenta and blue curves in each panel (which are, respectively, identical in the four panels) are the 42-year climatology for *M* and *MU*, respectively.

**Figure 4.** Yearly time series of the lead days of *MU* with respect to *M* when their negative correlations reach maximum values in years from 1980 to 2021, derived from observational fields.

Displayed in Figure 5 is an idealized elliptic orbit with *α* = −0.25*π*. It is seen that when the Θ increases from 0 to *π*, *X* increases to a maximum while *Y* decreases to a minimum. Then, from Θ = *π* to 2*π*, *X* decreases to minimum while *Y* increases to maximum. If we consider *MU* as the X and *M* as the *Y*, such changes that accompany Θ = 0 to 2*π* are quite consistent with their annual cycle from the previous year, from summer to winter and then to the concurrent summer.

**Figure 5.** (**a**) Orbits of two variables, *X* and *Y*, satisfy the idealized parametric ellipse orbit model, with the tilting angle of the ellipse with the horizontal axis (*α*) equaling −0.25 *π*. (**b**) Lead days of *X* with respect to *Y* when their negative correlations reach maximum values as a function of the ratio of *b* to a and the tilting angle *α*.

According to the parametric elliptic orbit model, the maximum negative correlations between *X* and *Y* can occur at different lead-lag times, which highly depend on the ratio b/a and the angle *<sup>α</sup>*. For the *<sup>α</sup>* in the range of 0~−*<sup>π</sup>* <sup>2</sup> , a larger *b*/*a* would yield a longer lead time for the *X* relative to *Y* to achieve their maximum negative correlation. Contrarily, for the *<sup>α</sup>* in the range of <sup>−</sup>*<sup>π</sup>* <sup>2</sup> ∼ −*π*, a smaller *b*/*a* would yield a longer lead time for the *X* relative to *Y*.

We apply a MATLAB function, "fit\_ellipse" to a pair of daily time series of *MU* and *M* from 1 July to 30 June of the next year to determine the parameters *X*0, *Y*0, *a*, *b* and *α* with <sup>Θ</sup> <sup>=</sup> <sup>2</sup>*π<sup>t</sup>* <sup>365</sup> , where *<sup>t</sup>* starts on 1 July as *<sup>t</sup>* = 1 and ends on June 30 as *<sup>t</sup>* = 365 (29 February in a leap year is excluded in our analysis). The information from "fit\_ellipse" can be found at https://www.mathworks.com/matlabcentral/fileexchange/3215-fit\_ellipse (accessed on 15 January 2023). To better capture the observed lead time information of *MU* with respect to *M*, we introduce an auxiliary equation when applying the "fit\_ellipse" function such that the constructed time series of *MU* and *M* from the fitted elliptic orbit meet the same lead time of the observed *MU* with respect to the observed *M*. As illustrated in Figure 5, when the value of *α* is fixed, the lead-lag days of maximum correlation and the ratio b/a have a one-to-one correspondence relation. The climatological mean annual cycles of *MU* and *M* and their elliptic orbit are displayed in Figure 6. The close resemblance between the fitted orbit elliptic (dashed magenta) and the scattering plot of the observed *MU* and *M* (colored dots) indicates that an elliptic orbit is capable of capturing the annual evolutions of *MU* and *M*. In particular, the fitted *MU* also varies out of phase with the fitted *M*.

**Figure 6.** Scatter plot of the climatological mean *M* versus *MU* (colored dots) in the period 1980–2021 and its fitted ellipse orbit (dashed black curve).

#### *3.3. Year-to-Year Variations of the Elliptic Orbits for the Annual Evolutions of M and MU*

In this section, we will examine the effectiveness of the parametric elliptic orbit model in capturing the inter-annual variations in the annual evolutions of the *MU* and *M*. The yearly time series of the parameters derived from the ellipse orbit fitting function are shown in Figure 7. It is seen that each of these yearly time series exhibits a strong yearto-year variability. In particular, they have pronounced differences to their counterparts derived from the climatological orbit (dashed magenta), indicating a strong year-to-year variation in the elliptic orbits. The generally larger values of *X*<sup>0</sup> and smaller values of *Y*<sup>0</sup> than their counterparts for the climatological elliptic orbit are observed, and the general out-of-phase relationship between *MU* and *M* is also present in the year-to-year variation in their annual means. The generally larger values of a and b than their climatological counterparts indicate that the amplitudes of the seasonal cycles of individual years are stronger than those of the climatological annual cycle. The generally larger values of b/a and α than their climatological counterparts indicate that the lead days of *MU* with respect to *M* in individual years are longer than for their climatological counterpart.

**Figure 7.** Yearly time series of the parameters of the fitted elliptic orbit model for *M* and *MU*. (**a**) The parameter *X*0, (**b**) the parameter *Y*0, (**c**) the parameter *a*, (**d**) the parameter *b*, (**e**) the parameter *α*, and (**f**) the ratio of *b* to *a*. The horizontal dashed red lines are the values of the parameters for the fitted orbit of the climatological mean annual cycles of *MU* and *M*.

It is seen from Figure 8 that the fitted elliptic orbits of individual years exhibit similar patterns to the yearly scatter plots of *MU* versus *M* in terms of both timing and intensity. In particular, the fitted orbits capture faithfully the interannual variations in extreme values of both *MU* and *M* in the winter seasons. For example, the fitted orbits capture the weaker SPV in the winters of 1982/83, 2001/02, 2002/03, 2003/04, 2006/07 and 2013/2014, which correspond to the panels in Figure 8 (or the abscissa's tick marks in Figure 9), labeled as 1982, 2001, 2003, 2006 and 2013, respectively. The year-to-year annual evolutions of *MU* and *M* obtained from the fitted elliptic orbits are shown in Figure 9. In comparison with Figure 2, it can be seen that the fitted elliptic orbits can capture both the annual evolutions and their interannual variability in *MU* and *M*, as well as the generally negative correlation between *MU* and *M* at both intraseasonal and interannual scales.

**Figure 8.** Scatter plot (color dots) of the *M* and *MU* in all years from 1980 to 2021 and the fitted elliptic orbit (black orbit).

Next, we quantitatively compare the extremes of *MU* and *M* in the winter seasons and their annual-cycle amplitudes obtained from the fitted orbits with the original time series of *MU* and *M*. Figure 10a,b shows that the yearly extreme values of both *MU* and *M* in winter seasons obtained from the fitted orbits have a high positive correlation (exceeding 0.94) with their counterparts obtained from the 31-day running mean time series of observed *MU* and *M*. The slightly larger values of the maximum *MU* and minimum *M* indicate that the fitted orbits tend to overestimate the strength of SPV in the winter seasons. The year-to-year variations in the amplitude of the annual cycle of *MU* and *M* (defined as the differences between the maximum and minimum values in the annual evolutions of *MU* and *M*) are displayed in Figure 10c,d, where the orbit-fitted amplitudes are compared with the observed ones. The figures show that the fitted orbits can almost perfectly capture the observed year-to-year variations in the annual cycle amplitudes for both *MU* and *M*, as is evident from the positive correlation exceeding 0.9. Next, we examine if the parametric elliptic orbit model can faithfully reflect the coupling between *MU* and *M* at the interannual time scale, as it does at the annual time scales (e.g., Figure 9 versus Figure 2). Figure 10e shows the scatter plot of the yearly time series of maximum *MU* and minimum *M*, which jointly measure the year-to-year variations in the strength of the SPV in the winter seasons. It is seen that the correlation between the yearly time series of *MU* and *M* derived from the fitted orbits (red circles) is identical to that derived from the observations (blue circles), equaling −0.78. Therefore, the observed negative correlation between the interannual variability in *MU* and that in *M* is well captured by the fitted orbits. Figure 10f shows the scatter plot of the yearly time series of the amplitude of the annual evolutions of *MU* and *M*. The observations show that the annual evolution amplitude of *MU* of individual years is positively and strongly correlated with that of *M* (about 0.69) and so are the fitted orbits (0.72) and the observation (0.69). It can be summarized from Figures 8–10 that both the annual evolutions and their yearly variations in *MU* and *M* are well coupled with significant correlations (about or larger than 0.7), which can be captured with almost identical correlations by the orbit fitting method, with five parameters varying yearly, as illustrated in Figure 7.

**Figure 10.** Scatter plot of the orbit-fitted extremes (red circles), and the observed extremes of maximum value of *MU* and minimum value of *M* in each winter season. (**a**) The maximum of *MU* from orbit fitting (ordinate) versus the observation (abscissa), (**b**) the same as (a) but for the minimum of *M*, (**c**) the same as (a) but for the yearly difference between maximum and minimum values of *MU*, (**d**) the same as (a) but for the difference between maximum and minimum values of *M*, (**e**) the maximum of *MU* (abscissa) versus the minimum of *M* (ordinate) for observation (blue circles) and fitted orbits (red circles), and (**f**) the same as (**e**) but for the yearly differences between maximum and minimum values of *MU* and *M*.

#### **4. Conclusions**

The year-to-year varying annual evolutions of the stratospheric polar vortex (SPV) have an important downward impact on the weather and climate from winter to summer and thus potential implications for seasonal forecasts. This study considers the necessity of investigating both the thermal and dynamical conditions of SPV by jointly considering the daily time series of the zonally integrated mass-weighted zonal momentum at 60◦ N (*MU*) and the total air mass above the isentropic surface of 400 K (*M*) over the latitude band of 60–90◦ N. The annual evolutions of *MU* are characterized by a generally out-of-phase relationship with the annual evolutions of *M*, and the changes in *MU* tend to lead the changes in *M* by 1–10 days. Moreover, their annual evolutions have significant interannual variations, which also exhibit a strong negative correlation between *MU* and *M*.

By constructing a parametric elliptic orbit model to fit the daily time series of *MU* and *M* in each year, we can closely reproduce the year-to-year variations in the annual evolutions of the observed *MU* and *M* jointly from the year-to-year variations in their elliptic orbits, including their amplitude, extremes in winter seasons, as well the strong negative correlation between *MU* and *M* and their amplitudes in the annual cycle at the interannual time scales. The findings of this study may pave a new way for short-time climate forecasts of the annual evolutions of SPV, including its temporal evolutions over individual cold seasons (including the spring and fall seasons). In particular, one would be able to predict timings of the minimum *MU* and maximum *M* in a given year by predicting the five parameters (*X*0, *Y*0, *a*, *b*, and *α*) of the corresponding yearly elliptic orbit. Because the timings of minimum *MU* and maximum *M* correspond closely to the timings of SSW events, one could in turn predict the timings of the high probability occurrence of SSW events in winters. Such a circulation condition during this timing would provide a favorable background for the break-up of the stratospheric polar vortex, which may yield a higher probability of the occurrence of SSW, which will be a topic of our future studies.

We here wish to add that this study is focused on the seasonal cycle of the stratospheric polar vortex in the Northern Hemisphere. In the Southern Hemisphere, the annual cycle of the stratospheric polar vortex also plays an important role in extreme surface weather conditions and natural hazards, as it may raise the risk of increased rainfall in the latitudinal band of 35–50◦ S [63–65]. According to previous studies, the Southern Hemisphere stratospheric polar vortex tends to be steadier and varies at a longer period [9,66] and break-up events occur less frequently [65]. It is expected that the annual cycle of the stratospheric polar vortex as well as extreme events such as SSW in the Southern Hemisphere could be also well captured by our parametric elliptic orbit model, which will be verified in our future work.

**Author Contributions:** Conceptualization, M.C. and Y.Y.; methodology, Y.Y.; validation, J.S., M.S. and Y.Y.; formal analysis, J.S. and M.S.; investigation, M.C.; data curation, J.S.; writing—original draft preparation, Y.Y.; writing—review and editing, M.C., J.S., M.S. and X.L.; visualization, J.S.; supervision, M.C.; project administration, M.C.; funding acquisition, M.C. and Y.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was jointly funded in part by the National Science Foundation of China, grant number 42075052, and the Natural Science Foundation of Jiangsu Province, grant number BK20211288, Climate Program Office of National Oceanic and Atmospheric Administration (NA20OAR4310380), and the National Science Foundation (AGS-2032542).

**Data Availability Statement:** The ERA5 datasets used in this study are freely available on the official website. (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressurelevels?tab=overview; https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-singlelevels?tab=overview; accessed on 1 November 2022).

**Acknowledgments:** We would like to thank the three anonymous reviewers for their constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.
