We used the daily-mean sea surface height anomaly (SSHA) with a 0.25° × 0.25° resolution provided by the Archiving, Validation, and Interpretation of Satellite Oceanographic (AVISO) (Ducet et al., 2000). The sensible and latent heat fluxes are from the Objectively Analyzed Air–Sea Fluxes (OAFlux) conducted by the Woods Hole Oceanographic Institution, available in daily mean on a 1° × 1° grid (Yu and Weller, 2007). To obtain the anomalous atmospheric circulation related to the KE fluctuation, the reanalysis states, including SST, air temperature, geopotential height, snow cover, precipitation, and horizontal winds from the fifth major global ReAnalysis data from the European Centre for Medium-Range Weather Forecasts (ERA5), are analyzed in this study. The ERA5 dataset has hourly intervals, with a horizontal resolution of 31 km and 137 levels in the vertical (https://www.ecmwf.int/en/newsletter/147/news/era5-reanalysis-production). The present study focuses on the wintertime spanning from 1 December to 28 February (DJF) of the following year. Here, the first winter of 1993 denotes December 1993 and January and February 1994, and a total of 27 winters from 1993 to 2019 are analyzed in this study. To eliminate disturbances from the climate trend, all the data are linearly detrended before proceeding with further analysis. Since the KE index has non-negligible autocorrelations, a Student’s t-test with an effective degree of freedom proposed by Bretherton et al. (1999) is used to evaluate the significance of variable anomalies.

Following Qiu et al. (2014), the KE index (KEI), defined as the 1-year lowpass-filtered SSHA averaged in 31°–36°N, 140°–165°E, is employed to characterize the interannual KE dynamical state. Previous studies have extensively used this identification method to fully capture the KE variation feature (e.g., Sun et al., 2022; Usui and Ogawa, 2022).

The trough-tracking algorithm proposed by Sen et al. (2019) is applied to extract the MedT axis in a region of 23°–57°N and 0°–70°E where the MedT is mostly observed. First, we locate the minimum z₅₀₀ at a specific latitude of 39.75°N between 0°E and 70°E. Then, the sequential minimum z₅₀₀ to the south and north is detected in a zonal length of 12° with the minimum at the previous latitude as the center. If no minimum can be located or the detected point is out of the domain interested in a direction, the algorithm is determined in this direction. Finally, we connect all the points selected to obtain the MedT axis.

We conducted the empirical orthogonal function (EOF) technique in the MedT region (25°–55°N, 0°–65°E) to extract the dominant modes and principal components (PCs) of the 500-hPa geopotential height variability of the MedT. The EOF analysis has been widely used in obtaining the dominant modes of atmospheric fields of various monsoon systems, such as the East Asian winter monsoon (Chen et al., 2014; Shi, 2021), EAT (Wang et al., 2009; Leung and Zhou, 2015; Chen et al., 2018), and MedT (Sen et al., 2019).

As the upper troposphere is likely to be equivalent barotropic over the midlatitudes, the atmospheric anomalous processes are explained through the barotropic vorticity equation. The anomalous Rossby wave source (RWS) is diagnosed as follows (Sardeshmukh and Hoskins, 1988):

where V _χ is the divergent winds and ξ is the absolute vorticity. The overbars and primes denote the wind-mean states and perturbations from the mean states, respectively.

In the middle and upper troposphere, the wave activity flux (WAF) is employed to calculate the wave propagation of stationary eddies on the zonally varying basic flow (Takaya and Nakamura, 1997; Takaya and Nakamura, 2001). The phase-independent WAF, which parallels its local group velocity, is used to depict the Rossby wave activity associated with the KE fluctuation. The formula is as follows:

Here, (u,v) is the winter-mean horizontal wind over the period from 1993 to 2019, | U → | is the speed for climatological horizontal flow, ψ' is the geostrophic stream function perturbation, f is the Coriolis parameter, R is the gas constant, p is the air pressure, and σ = ( R T ¯ / C p p ) − d T ¯ / d p is based on the temperature (T ) and the specific heat at constant pressure (C _p ). The subscripts x , y , and p indicate the partial derivatives along the zonal, meridional, and vertical directions, respectively.

The maximum Eady growth rate (EGR), proposed by Hoskins and Valdes (1990), is used to estimate the atmospheric baroclinicity and determine whether the atmospheric dynamical state is conducive to the weather system origination. We calculated the EGR at the 300-hPa level using the equation of

with f being the Coriolis parameter, N being the Brunt–Väisälä frequency, U being the horizontal wind speed, and z being the vertical height.

To effectively extract the synoptic-scale transient eddies, the meridional wind variance v'v' at the 300-hPa level subject to a 2.5–6-day bandpass 31-point filter (Sun and Zhang, 1992) is applied to represent the storm track.

To evaluate the geopotential height anomaly generated by low-frequency atmospheric waves and high-frequency eddies associated with the KE variation, we applied the equation given by Lau and Nath (2014) to calculate the Rossby wave-induced and transient eddy-induced geopotential height tendency at the 500-hPa level. The equation is shown as follows:

where Z is the monthly geopotential height, f is the Coriolis parameter, g is the gravitational acceleration setting to 9.80 m/s², and V' and ζ‘ are the wind and relatively vorticity disturbances, respectively, subjected to an 8–90-day Butterworth bandpass filtering for Rossby wave and a 2.5–6-day bandpass filtering for synoptic eddies. The contributions of Rossby waves and transient eddies to the observed geopotential height tendency can be quantified by comparing their magnitudes and spatial patterns. This method has been extensively used in calculating the effects of Rossby waves and transient eddies on the geopotential height anomaly (Song et al., 2016; Sun et al., 2022).

Over midlatitudes, the feedback response of the atmosphere to the extratropical forcing in the ocean and land is difficult to evaluate, because the atmospheric response, to a large extent, is overwhelmed by the tropical remote forcing (i.e., ENSO) and atmospheric internal variability. To extract the midlatitude oceanic feedback on the atmosphere, we follow the lag regression approach of Révelard et al. (2016), which separates the slow atmospheric response to the SST from the rapid atmospheric internal variability. Hence, any atmospheric variable x(t) consists of three components, as follows:

Here, αKEI(t−2) and βe(t−1) represent the quasi-equilibrium atmospheric responses to the KE 2 months earlier and ENSO 1 month earlier, respectively, with and β as the feedback parameters, KEI as the KE index described above, and e as the Niño 3.4 index (defined as the 5-month running mean SST over the region of 5°S–5°N, 170°–120°W). n(t) represents the atmospheric white noise. It generally takes 1–2 months for the extratropical SST-forced baroclinic signal to fully transform into the large-scale equivalent barotropic circulation and reach its maximum amplitude via the transient eddy–mean flow interactions (Liu et al., 2007; Ferreira and Frankignoul, 2008; Liu et al., 2008; Révelard et al., 2016). Indeed, the KE fingerprint on the large-scale atmospheric circulation has some persistence as the atmosphere lags the ocean by 1 to 6 months, but the strongest imprint emerges when the delay is set to be 2 months (not shown). In addition, we have assumed a 1-month delay for ENSO signals reaching the Mediterranean region, considering the maximum ENSO teleconnection patterns. For the above reasons, the atmosphere response times to the KE and ENSO are regarded as 2 and 1 months, respectively, in our paper, consistent with the lag time applied in Révelard et al. (2016) and Sun et al. (2022).

To access the feedback parameter α , first, we eliminate the ENSO signal, yielding the ENSO-filtered atmospheric variables and KEI , expressed as follows:

Here, A=x(t),e(t−1)/e(t−1),e(t−1) , and B=KEI(t), e(t+1)/e(t+1), e(t+1), with p,q representing the covariance of p and q . Then, given that the KE variability has no correlation with the atmospheric intrinsic variability, we substitute the and KEI in Equation (5) with Equations (6) and (7), leading to the following:

After multiplying K E I ( t − 2 ) ˜ on two sides of Equation (8), we take the ensemble average and obtain as follows:

The feedback parameter α for the atmospheric variable shows similar spatial patterns with the composite differences between the KE stable period (SP) and unstable period (USP) ( Supplementary Figures ), validating the effectiveness of the lag regression approach that this paper applied.

Since this paper focuses on the atmospheric feedback in winter (DJF) to the KE variability 2 months earlier, Figure 1 highlights the time series of KEI from October to December (OND) and the corresponding SSHA patterns for the two opposite KE states. As depicted in Figure 1A , the KEI exhibits a feature of bimodal interannual fluctuation, with a low-frequency shift between positive and negative phases. The OND SSHA pattern, with correspondence of KEI larger than 0.5, is characterized by a prominent positive SSHA anomaly to the south of the KE axis and a weaker negative SSHA to the north, suggesting a stronger and steadier KE jet, hereafter KE SP (OND in 2002, 2004, 2010–2014, 2016, and 2018–2019) ( Figure 1B ). The KE USP is selected when the OND KEI is less than −0.5 (OND in 1993–1997 and 2006–2008), which shows a basically opposite SSHA pattern with smaller amplitude compared to that in the KE SP ( Figure 1C ). The KE variation feature during 1993–2012 above is consistent with Qiu et al. (2014) and sustains a stable dynamical state from 2012 to 2019 (Usui and Ogawa, 2022).

Figure 2 shows the spatial distribution of year-to-year winter MedT axes and corresponding z₅₀₀ along the trough lines for KE SP and USP. The MedT location is highly unstable zonally, shifting in the east–west direction substantially from year to year. When the KE is in its stable (unstable) state, the MedT axis is situated to the west (east) of its climatological location, with stronger (weaker) intensity along the trough line. On average, the longitude and z₅₀₀ differences of MedT between the KE SP and USP are approximately −12° and −15 gpm, respectively, pointing to a close connection between KE and MedT dynamical states. To further investigate this relation, Figure 3 displays the KE-induced atmospheric anomalies derived from the lag regression approach. Indeed, the KE exerts significant impacts on the atmospheric circulation in the Euro-Mediterranean region. An anomalous cyclone center of −18 gpm, slightly located to the west of the MedT climatological position, is developed in the middle troposphere over the Euro-Mediterranean region ( Figure 3B ), along with accelerated (decelerated) westerly winds to its south (north) ( Figure 3D ). The atmospheric anomaly on the 250-hPa level bears similar distributions to that on the 500-hPa level ( Figures 3A, C ), which implies that the atmospheric response to the KE fluctuation is equivalent to barotropic. Therefore, the MedT tends to be stronger and locates farther west in the KE SP, and the opposite for the USP.

By performing EOF analysis, the first three leading modes of MedT variability and their relations with KEI are given in Figure 4 . Here, the area of 25°–55°N, 0°–65°E is chosen for EOF analysis, taking into account the region where the MedT is mostly detected and the location of 500-hPa geopotential height anomaly in Figure 3C . Our results are broadly consistent with Sen et al. (2019), except that the order of the first and second modes is contrary to their results, probably due to the different periods and regions we focused on. The first eigenvector of detrended and normalized z₅₀₀ fields during 1993–2019, which contributes to 31.86% of the total variance, depicts that this mode corresponds to a pronounced anomaly extending northeastward from the eastern Mediterranean Sea to the Caucasus Mountains ( Figure 4A ). The correlation coefficient between the first PC time series (PC1) and MedT intensity index (the average z₅₀₀ in the box of 30°–50°N and 20°–45°E) proposed by Sen et al. (2019) is 0.85 (p< 0.001), revealing that the leading EOF mode can well represent the MedT intensity variation. Figure 4B shows there exists a relatively weak co-variation between the PC1 and KEI. After applying a Savitzky–Golay filter to extract the lowpass MedT signal, the relation between the PC1 and KEI is strengthened, with the correlation coefficient increasing to 0.57 (p< 0.01), demonstrating that the KE fluctuation has some potential influences on the MedT intensity. The second PC (PC2) time series is found to be highly correlated with the MedT zonal displacement index (the z₅₀₀ difference between 39°N, 10°E and 39°N, 50°E) developed by Sen et al. (2019), with a correlation coefficient of 0.86 (p< 0.001). Therefore, the second EOF mode, which accounts for 31.37% of the total variance, mainly explains the zonal displacement variability of MedT ( Figure 4C ), also having linkage with the lowpass KEI (r = 0.40, p< 0.1) ( Figure 4D ). The third EOF mode, which is responsible for 13.68% of the total variance, is hard to be identified as far as its physical implication, but it shows a similar spatial pattern with the oscillation between the eastern and western Mediterranean (Xoplaki et al., 2004) ( Figure 4E ). The correlation between the third PC (PC3) time series and KEI is weak and insignificant (r = 0.28, p > 0.1) ( Figure 4F ). In summary, there exists a significant relationship between the variability of the KE dynamical state and the intensity/zonal displacement of MedT, which raises a new question: how does the KE fluctuation affect the MedT intensity/zonal displacement?

Some studies suggested that the KE interannual variability triggers the mesoscale SST anomalies in the KE region, thus modulating the atmospheric large-scale circulation by generating the atmospheric Rossby wave activity, transient eddy–mean flow interaction, and thermal winds (Ma X. et al., 2015; O’Reilly and Czaja, 2015; Révelard et al., 2016; Lei et al., 2020; Sun et al., 2022). In our study, the Rossby wave and eddy–mean flow interaction are likely to be the possible underlying mechanisms. Therefore, the aforementioned two physical processes will be analyzed in the following subsections separately.

Figure 5 shows the estimated responses of SST and THF to the KE variability in the cold season. When the KE is stable, the SST decreases (increases) by 0.8°C (0.4°C) to the west (east) of 155°E in the KE region. The THF anomaly bears some resemblances to the SST anomaly, showing a corresponding decreased (increased) THF of −30 W m⁻² (20 W m⁻²), which is equivalent to nearly 10% of the climatological THF in wintertime. Such a large THF anomaly is supposed to exert pronounced impacts on the atmosphere.

The KE fluctuation exerts prominent influences on the surface parameters in the Euro-Mediterranean region, as depicted in Figure 10 . In the winter with a stable KE, the 2-m air temperature shows a dipole-like pattern, with colder air temperature over the Euro-Mediterranean and North Africa and warmer air temperature over northeastern Europe and western Central Asia, consistent with the stronger and western-positioned MedT. The coldest temperature anomaly of −0.8°C is situated over the Balkan Peninsula and Romania, which is upstream of the MedT axis in the KE SP ( Figure 10A ). The change of snow cover displays basically the same pattern as air temperature, with snow accumulating (diminishing) to the west (east) of the climatological MedT position ( Figure 10B ). The convective precipitation in the northeastern Mediterranean Sea and large-scale precipitation in the Anatolian Peninsula are enhanced ( Figures 10C, D ), mainly due to the intensified water vapor transport driven by the accelerated westerlies from the Mediterranean Sea and stronger upward motion resulting from the deeper and western-positioned MedT (not shown). The aforementioned climate variability is highly correlated with the intensity and zonal displacement changes of MedT, implying a strong linkage between the KE fluctuation and the Euro-Mediterranean climate.