Introduction

Despite the profound impacts of midlatitude weather extremes, the mechanisms and underlying causes for these events in a changing climate remain an active topic of debate. Regional extremes of surface temperature or air pollution are frequently linked to large-scale atmospheric dynamics that feature the meandering of midlatitude westerly jets (e.g., Buehler et al. 2011;Pfahl and Wernli 2012;Screen and Simmonds 2014;Martineau et al. 2017;Sun et al. 2019). Under climate warming, it was hypothesized that Arctic amplification or the enhanced surface warming in the Arctic compared to lower latitudes would weaken the midlatitude jets, causing more frequent jet meanders and associated extremes (e.g., Francis andVavrus 2012, 2015). However, not only the observed trends in jet waviness are obfuscated by the internal variability of the climate system (e.g., Blackport and Screen 2020;Nie et al. 2023) but also the weakened midlatitude jets are offset by the zonal wind acceleration due to tropical upper-tropospheric warming (e.g., Shaw et al. 2016;Chen et al. 2020;Nie et al. 2022). The net effect yields a small decrease in jet waviness toward the end of the twenty-first century, albeit with large uncertainties (Cattiaux et al. 2016; Barnes and Polvani 2015;Peings et al. 2017Peings et al. , 2018;;Woollings et al. 2018;Nie et al. 2023). Recently, Chen et al. (2022) developed a tracer-based wave activity for quantifying jet meandering, and Nie et al. (2023) showcased that waviness subject to changes in time mean zonal winds can be constrained by a potential vorticity (PV)-like passive tracer in an advection-diffusion model. In this paper, we will show the eddy mixing time scale as a key measure of the Northern Hemisphere winter circulation waviness.

Midlatitude weather systems are typically understood in terms of generation, propagation, and dissipation of baroclinic eddies. Baroclinic instability of the atmosphere converts the available potential energy (APE) of the atmosphere to eddy kinetic energy (EKE) of synoptic storms, which, in turn, propagate eastward in the form of Rossby wave packets due to their faster group speed than individual phase speeds (e.g., Simmons and Hoskins 1978;Orlanski and Katzfey 1991;Chang and Yu 1999;Chang 1999;Chang et al. 2002;Wirth et al. 2018). In the conservative limit, the local rate of change in wave activity is equal to the convergence of wave activity flux, and in the steady state, the wave activity flux diverges away from the region of wave generation and converges into the region of wave dissipation (e.g., Edmon et al. 1980;Plumb 1985Plumb , 1986)). While EKE and wave activity are intricately related, the spatial distributions of the two quantities, however, are largely distinct. Compared with storm tracks or the regions of large EKE (e.g., Chang et al. 2002), the climatological pattern in local wave activity (LWA) is characterized by local maxima to the north of North Pacific and Atlantic storm tracks or near the exits of storm tracks (e.g., G. Chen et al. 2015;Huang and Nakamura 2017;Martineau et al. 2017). Indeed, the LWA pattern is similar to that of blocking frequency [see the review by Woollings et al. (2018)] or the frequency of high-latitude Rossby wave breaking (Strong and Magnusdottir 2008;Michel and Rivière 2011;Barnes and Hartmann 2012;Martineau et al. 2017). This spatial separation between EKE and LWA may be attributed to different stages of the baroclinic eddy life cycle (e.g., Simmons and Hoskins 1978), in which EKE is closely related to the baroclinic growth of eddies, but blocking or Rossby wave breaking occurs during the barotropic decay of eddies. This paper is focused on the "waviness" of the westerly flow that is more relevant to blocking or wave breaking.

Numerous diagnostics have been developed to quantify the waviness (e.g., meridional displacement or sinuosity) of circumglobal contours in 500-hPa geopotential height (Z 500 ) (e.g., Francis and Vavrus 2012;Hassanzadeh et al. 2014;Francis and Vavrus 2015;Cattiaux et al. 2016). From the perspective of stagnant weather regimes, the primary circulation regimes of the Western Hemisphere and the European sector are characterized by blocking over the Alaska and Greenland regions, respectively (e.g., Cheng and Wallace 1993;Bao and Wallace 2015;Lee et al. 2019), and the cold spells over East Asia are closely related to blocks over the Ural Mountains (e.g., Luo et al. 2018;Wang et al. 2021). This emphasis on persistent circulation waviness is consistent with the detection of blocking, which identifies persistent anomalies in Z 500 (e.g., Dole 1986;Nakamura et al. 1997) or the latitudinal reversal of Z 500 or PV (e.g., Tibaldi and Molteni 1990;Pelly and Hoskins 2003). Importantly, Huang and Nakamura (2016) developed a theoretical framework for the budget of LWA from Kelvin's circulation theorem, which provides dynamical insights on the variability of midlatitude storm tracks (Huang and Nakamura 2017) and the onset of blocking events (Nakamura and Huang 2018;Paradise et al. 2019). To better compare LWA with other diagnostics of Z 500 , G. Chen et al. (2015) proposed the use of the LWA diagnostic for the waviness of daily Z 500 . As Z 500 is related to the lower-tropospheric temperature through the hypsometric equation, this can be thought of as the wave activity of lower-level temperature. Recently, Geen et al. (2023) examined an idealized atmospheric model in response to polar amplification and concluded that the changes in the geometric structure of Z 500 are consistent among different waviness metrics.

The behavior of jet meandering can also be understood through the lens of a PV-like passive tracer (see discussions about a passive tracer versus PV in section 6) in an advection-diffusion model (Chen et al. 2022;Nie et al. 2023). Eddy diffusivity, often involving Rossby wave breaking, is a fundamental property of the atmosphere that redistributes dynamical tracers such as PV or potential temperature (e.g., Green 1970;Plumb 1979;Stone and Yao 1987;Held 1999). Nakamura et al. (1997) showed that an idealized model of "contour advection" for a PV-like tracer can disentangle eddy-mean flow interactions during the blocking evolution. Chen et al. (2020) found that modifying the zonal mean zonal wind for the advection of vorticity and temperature in an atmospheric dynamic core can isolate eddy feedback. Woollings et al. (2023) proposed that large-scale polar flow can be conceptualized as a mix of geostrophic turbulence and Rossby wave propagation. Compared to the LWA budget that emphasizes the jet's capacity for blocking formation (Nakamura and Huang 2018), this diffusive viewpoint characterizes a strong zonal jet as a "mixing barrier" (e.g., Haynes and Shuckburgh 2000;Allen and Nakamura 2001;Chen and Plumb 2014;Abalos et al. 2016) or a barrier to blocking. This is also evident for the leading mode of jet variability, featuring weaker eddy mixing at positive zonal wind anomalies and vice versa (Nie et al. 2014;Burrows et al. 2017).

This paper highlights the time scale of jet meandering based on the theory of eddy diffusivity scaling. This is illustrated for a simple case below and more realistic Rossby wave packets in section 3. For a zonally symmetric background flow u, a linear theory of eddy mixing can be formulated by considering a single wave solution of the form exp[ik(x 2 ct)], where x is the zonal position, t is the time, k is the zonal wavenumber, and c is the eddy phase speed. If the intrinsic frequency of Rossby waves is denoted as v i , which describes the eddy frequency in the frame of reference moving with the mean flow u, then eddy diffusivity can be written as [i.e., rewrite Eq. (1) of Chen et al. (2022) in the formfoot_0 of D 5 tEKE and see the derivation in their appendix]

Here, EKE 5 yfoot_1 , t is the decorrelation time scale of eddies, and t e is referred to as the effective time scale of eddies resulting from nonzero intrinsic frequency v i , which can be thought of as eddy persistence in the reference frame that moves at the mean flow speed.

The effective eddy time scale in Eq. ( 1) has two important limits. (i) In the limit where v 2 i 5 0 or the eddy phase speed is equal to the mean zonal wind speed (c 5 u), the effective time scale is reduced to t e ' t. This corresponds to the regime of isotropic turbulence, in which eddy diffusivity scales as eddy characteristic speed (EKE 1/2 ) times mixing length (tEKE 1/2 ) (e.g., Held 1999). (ii) In the limit where v 2 i 5 (c 2 u) 2 k 2 . . t 22 or the mean zonal speed is much greater than eddy phase speed, the effective time scale is t e ' v 22 i t 21 , , t. This indicates that a strong zonal jet suppresses the effective eddy time scale compared with isotropic diffusivity. This may be regarded as a wave regime rather than a turbulence regime because a strong jet acts as a waveguide for zonal Rossby wave propagation. Chen et al. (2022) and Nie et al. (2023) found that this linear theory helps to understand the waviness of the PV-like tracer in an advection-diffusion model in response to changes in time mean zonal winds. This paper will expand this single-wavenumber theory to realistic wave packets, with a focus on the eddy mixing time scale.

The paper is organized as follows. In section 2, we introduce the reanalysis and model data and outline the diagnostics of waviness for Z 500 and a PV-like tracer. Section 3 describes the effective time scale for circulation waviness from the eddy mixing theory, which is derived in the appendix. In section 4, we use the eddy mixing theory to explain the spatial sensitivity of effective time scale to time mean zonal winds. This theory is further exemplified in section 5 for the interannual variability of LWA in the reanalysis and the climate response of LWA in climate models. Section 6 offers brief discussions on the causality for LWA in kinematic and dynamical models. Section 7 concludes the paper.

Data and methods

a. Reanalysis data and model simulations

We use 6-hourly data from the fifth major global reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ERA5) (Hersbach et al. 2020). The ERA5 data are analyzed at the resolution of 18 longitude 3 18 latitude for the winters [December-February (DJF)] in 1980-2020. To be consistent with climate model outputs, the 6-hourly geopotential height is averaged to daily data before computing wave activity. The meridional wind at 1200 UTC is analyzed for the lagged correlation analysis of eddies.

We also employ initial-condition large ensembles from four coupled atmosphere-ocean models: CESM1 (40 members) (Kay et al. 2015), CanESM2 (50 members) (Kirchmeier-Young et al. 2017), GFDL CM3 (20 members) (Sun et al. 2018), and CESM2 (50 members) (Rodgers et al. 2021). The ensemble members of each climate model are driven by the same climate forcings but with only small differences in initial conditions, and thus, the differences between ensemble members are attributed to the internal variability of the climate system. To evaluate the projected future changes at the end of the twentyfirst century, we calculate the difference between the historical period and the future period (2080-99).

b. Diagnostics of circulation waviness

The diagnostics of jet waviness in the literature have been designed to characterize various aspects of the westerly flow (e.g., Francis and Vavrus 2012;Hassanzadeh et al. 2014;Francis and Vavrus 2015;Cattiaux et al. 2016;G. Chen et al. 2015;Huang and Nakamura 2016). As PV is conserved following the Lagrangian motion of an air parcel, we focus on the dynamics-based metrics related to the geometric characteristics of PV.

1) LWA

The LWA diagnostic was introduced by Huang and Nakamura (2016) to quantify the wave activity of a PV snapshot through Kelvin's circulation theorem. To better compare with other metrics of circulation waviness in the literature (e.g., Francis and Vavrus 2012;Hassanzadeh et al. 2014;Cattiaux et al. 2016), G. Chen et al. (2015) noted that Z 500 is approximately linearly related to PV at 500 hPa and demonstrated that LWA could quantify the waviness of daily Z 500 . In comparison with other methods that are designed for the waviness over preferred regions, this method produces an instantaneous two-dimensional (longitude 3 latitude) map of circulation waviness. Specifically,

Here, l is the longitude, f is the latitude, and a is Earth's radius. The term f e (Z 500 ) is the equivalent latitude of the geopotential height contour with the value Z 500 . The term ẑ(l, f) 5 z(l, f) 2 Z 500 (f e ) represents the eddy perturbation to the basic state of geopotential height Z 500 (f e ).

2) TRACER WAVINESS Similar to PV, the spatial distributions of chemical tracers in the atmosphere and oceans are largely regulated by transport and mixing in the flow (e.g., Green 1970;Plumb 1979;Stone and Yao 1987;Held 1999;R. Chen et al. 2015). Since a passive tracer advected by observed winds exhibits a similar geometric structure to PV in observations, another measure of the PV geometry is through the Lagrangian conservation of a PV-like passive tracer that is advected by observed winds (e.g., Haynes and Shuckburgh 2000;Allen and Nakamura 2001;Abalos et al. 2016). As the tracer geometry is only a function of prescribed advecting velocities after some initial adjustment in the tracer field, the tracer model enables us to attribute the changes in tracer geometry to the changes in advecting winds. Using an advection-diffusion model of passive tracers driven by 250-hPa reanalysis winds, Chen et al. (2022) and Nie et al. (2023) show that the waviness of a passive tracer, after mapping from tracer to PV (i.e., coordinate transformation), can describe the spatial patterns of circulation waviness and their changes under climate change.

In contrast to the meridional PV gradient that is sustained by differential heating between the equator and poles, the meridional gradient of a passive tracer decays in time. The tracer wave activity itself is largely affected by the vanishing tracer gradient, and thus, it is not a good metric for circulation waviness. Rather than mapping from tracer to PV [i.e., a transformation that replaces the tracer gradient with the PV gradient as in Eq. ( 9) of Chen et al. ( 2022)], we define tracer waviness as the tracer wave activity divided by the tracer gradient (Nie et al. 2023).

q cos fdf :

(3)

Here, Q tracer denotes the PV-like tracer from an advectiondiffusion model [see Chen et al. (2022) for details of the tracer model and references therein], f e (Q tracer ) is the equivalent latitude for the tracer contour with the value Q tracer , and q(l, f) 5 q tracer (l, f) 2 Q tracer (f e ) describes the deviation of the tracer field from the basic state Q tracer (f e ) due to horizontal advection.

For small-amplitude waves, q(l, f) ' 2 fQ tracer /f e , where f is the displacement of the tracer contour away from latitude f 5 f e ; substituting this into Eq. ( 3) yields W tracer ' f2 . This ascertains that W tracer describes the geometric property of a tracer field. Importantly, the tracer model does not rely on the assumptions of the linear theory in Eq. ( 1). Generally, LWA Z 500 and W tracer represent different but interconnected aspects of the flow, and thus, it is plausible that a scaling relationship exists for their spatial pattern and variability. To the extent that they are affected by the same synoptic weather systems, the fractional changes of LWA at 500 hPa may be related to the tracer waviness as

where overbars describe the time means in the historical period and d denotes the deviations from time means. The deviations can be either interannual variability or long-term trends.

Diffusivity scaling and effective time scale for circulation waviness

While the tracer model can simulate the spatial distributions of waviness and their response to climate change, it does not offer a simple explanation for the eddy mixing that controls tracer geometry. We next consider a multiwavenumber theory of eddy diffusivity e.g., (e.g., Davis 1991; R. Chen et al. 2015), which is more appropriate for Rossby wave packets than the single-wavenumber theory in Eq. ( 1). For simplicity, we consider the time mean flow u(x, y) and eddy meridional wind y (x, y, t), where overbars are the time means and primes denote the anomalies. The time-mean flow is assumed to vary slowly on the spatial scale larger than the eddy length scale. While this assumption is not exactly satisfied, the scale separation between the mean flow and eddies is typically used in a linear theory. We also compare this linear theory with the tracer waviness in the advection-diffusion model.

We perform a Galilean transformation from longitudinal position x to a new x L coordinate that moves at the zonal wind speed u (see the schematic in Fig. 1a). This is denoted as x L (x, t) 5 x 2 ut, where x is the longitude, t is the time, and x L is the longitudinal position in the new frame of reference. Physically, the moving coordinate will produce a Doppler shift in eddy frequency, and the effects on the spectral power of eddies are discussed in the appendix. Then, the eddy diffusivity for a passive tracer can be derived as [see Eq. (A4) in the appendix] (e.g., Davis 1991; R. Chen et al. 2015)

where the eddy autocorrelation in the moving coordinate x L is evaluated at lag time Dt as

Note for the time lag Dt 5 0, x 5 x L , and EKE 5 y 2 (x) 5 y 2 (x L ). Furthermore, we will approximate Eq. ( 6) in the form of a damped oscillatorfoot_4 , with frequency v i and decay time scale t:

Substituting Eq. ( 7) into Eq. ( 5) yields

This recovers the single-wavenumber diffusivity in Eq. ( 1), but v i and t are estimated from a wave packet using the eddy autocorrelation in the moving x L coordinate. We note that the intrinsic frequency of a Rossby wave is negative, and the damped oscillator can only estimate the magnitude of v i . Similar to Eq. ( 4), we can estimate the fractional changes of LWA at 500 hPa in relation to eddy mixing from Eq. ( 5).

500 dt e t e 1 dEKE EKE 5 LWA Z 500 dt e du du t e 1 dEKE EKE :

Here, d denotes either interannual variability or long-term trends. We have made an approximation dt e 5 (dt e /du)du, where the sensitivity parameter dt e /du denotes the dependence of effective eddy time scale t e on the local zonal wind speed u.

More specifically, the diffusivity scaling for waviness is calculated as follows (see Fig. 1):

1) We first compute the one-point correlation map between the two time series y (x, t) and y (x 1 Dx, t 1 Dt) for each grid point at 250 hPa over the historical period (i.e., no coordinate transformation), denoted as r(Dx, Dt), where Dx is the longitudinal shift and Dt is the time lag. This yields a two-dimensional (2D) correlation matrix of eddy meridional wind as a function of Dx and Dt. In practice, we have set an upper limit for time lag as Dt max 5 min(10 days, {pa cos f/[u(l, f)]} 3 (1 day/86 400 s)), and the second limit in the brackets sets the maximum longitudinal shift to half a latitudinal circle for the time mean zonal speed u. As the eddy correlation map is calculated from daily data at 1200 UTC, the temporal resolution of the correlation map is 1 day. To facilitate the coordinate transformation in step 2, the temporal resolution is increased from 1 to 0.1 days by linear interpolation. 2) We next conduct a coordinate transformation from x to

x L . For a given Dt at each grid point, we linearly interpolate eddy correlation from the 2D map as a function of Dx and Dt to the straight line (uDt, Dt) (see the black line in Fig. 1b). This yields r Dt L 5 r(uDt, Dt), which represents the eddy autocorrelation at lag time Dt in the coordinate that moves at the zonal speed u [see Eq. ( 6)]. Physically speaking, this may be understood as the composite of decaying eddies with the phase speed of u, which may be thought of as the critical layer of Rossby waves in the linear theory (i.e., c 5 u).

3) Rather than integrating Eq. ( 5) directly, we find the best fit of r Dt L to a damped oscillator with frequency v i and decay time scale t [see Eq. ( 7)], denoted as rDt L . This is achieved by minimizing (r Dt L 2 rDt L ) 2 using MATLAB's optimization toolbox (see the black and red lines in Fig. 1c). This produces the climatological values of v i , t, and effective time scale t e . 4) To estimate the interannual anomalies or trends in t e in Eq. ( 9), we assume that the 2D correlation map r(Dx, Dt) in step 1 is fixed and that t e is only affected by the coordinate transformation in step 2 due to the anomalies in u.

That is, dt e 5 (dt e /du)du, with a sensitivity parameter dt e /du. The parameter dt e /du 5 dt e /du may be interpreted as differences in eddy time scale dt e between the composites of decaying eddies with the phase speeds u and u 1 du. This sensitivity parameter is estimated from the correlation map r(Dx, Dt), and the analysis procedure is described in detail in section 4 with respect to Fig. 4.

It is noteworthy that t e bears a simple physical meaning of eddy persistence in the reference frame that moves at the time mean zonal speed. For example, if the autocorrelation function in Eq. ( 5) decays with the e-folding time scale t as r Dt L 5 exp(2t 21 |Dt|), it follows that t e 5 t. Thus, a larger eddy time scale t e implies more persistent meridional meandering. For a typical Rossby wave packet, the mean wind speed is faster than the eddy phase speed (see the schematic in Fig. 1b). This suggests that for a Rossby wave packet, the slowdown of the time mean zonal speed, all else being equal, would enhance the eddy persistence in the meridional wind.

The spatial pattern of waviness time scale and its sensitivity to zonal wind

In this section, we use the eddy mixing theory [Eq. ( 5)] to understand the spatial distribution of circulation waviness and its sensitivity to zonal wind changes. Figure 2 displays the onepoint correlation maps of eddy meridional wind at the North Atlantic and Pacific jet cores, in comparison with the highlatitude regions over Greenland and Alaska where atmospheric blocking occurs frequently. The 2D eddy correlation at the jet cores exhibits the well-known eastward propagation of a Rossby wave packet, whereas the eddy correlation map around the high-latitude blocking is quasi-stationary, with a slow westward retrogradation. As the dispersion of a Rossby wave packet is influenced by mean flow, the seasonal mean zonal wind is marked in Fig. 2 by black lines. The interpolation of 2D eddy correlation along a black line gives the autocorrelation function in the frame of reference with respect to the seasonal mean zonal wind, which may be regarded as the composite of eddies with phase speeds equal to the time mean wind speed.

We further quantify the different characteristics of eddy autocorrelations (Figs. 2b,d,f,h) by the parameters of a damped oscillator [Eq. ( 7)]. Notably, the eddy decay time scale is greater at the two high-latitude grid points (t 5 2.5-3 days) than at the jet cores (t ; 0.9 days), but the intrinsic frequency is much greater at the jets (v i 5 4-6 day 21 ) than at the two polar points (v i ; 0.25 day 21 ). This produces v 2 i t 2 ; 13 2 29 at the jet cores, whereas v 2 i t 2 5 0:4 2 0:6 at two polar locations. The product v i t is also proportional to the ratio between the decay time scale of the wave envelope and the wave period. As such, the wave packets that exhibit multiple peaks (Figs. 2a-d) give v i t . 1, and the one-point correlation maps dominated by an exponential decay (Figs. 2e-h) yield v i t , 1. Using t e 5 t/(1 1 v 2 i t 2 ), the effective time scale is slightly smaller than the decay time scale in the high latitudes (t e ; 1.9 days) but becomes close to zero at the jet cores (t e 5 0.03-0.06 days). In other words, eddies in the high latitudes are more persistent than those at the jet cores. And the downstream development of anomalies at the jet cores does not imply more persistence relative to the mean jet. The decrease in t e /t by a factor of 10 from the polar regions to the jet cores is largely due to the increased intrinsic frequency in the jet regions.

Figure 3 compares the spatial patterns of EKE and LWA. Note that EKE is calculated as y 2 without any temporal filtering here, although the temporal filtering to synoptic time scales yields a similar pattern. Over the North Pacific and North Atlantic Oceans, high EKE values indicate intense storm activity developing downstream of the North Pacific and North Atlantic jet streams, respectively (Fig. 3a), consistent with different measures of storm tracks (e.g., Chang et al. 2002). By contrast, large values of LWA are observed near Alaska or Greenland to the north of the North Pacific or North Atlantic jet streams (Fig. 3b; also see G. Chen et al. 2015). These regions with large LWA magnitude roughly coincide with the primary centers of action for atmospheric blocking (e.g., Woollings et al. 2018).

Given that fast meridional wind speed or enhanced EKE is expected to produce a sizeable meridional displacement of the westerly flow, it is imperative to explain why large LWA tends to occur to the north of jet streams with muted EKE rather than at the jet latitudes with active eddies. This can be reconciled by examining the spatial distribution of eddy mixing, which approximately reproduces the regions of large LWA over Alaska and Greenland, as well as to some extent over western Europe (Fig. 3d). As eddy diffusivity scales as the product between EKE and effective eddy time scale [Eq. ( 5)], differences between EKE and LWA can be attributed to the spatial variation in effective eddy time scale. The examples in Fig. 2 provide evidence that the spatial pattern of effective time scale can be thought of as a transition from a wave propagation regime along the jet, which features high frequency and t e /t ' v 22 i t 22 , , 1, to a turbulence regime in high latitudes that is characterized by low frequency and t e /t ' 1. The regions with high values in frequency overlap the regions of strong zonal jets (Fig. 3c), which is expected from a Doppler shift of the dispersion relationship, |v i | 5 (u 2 c)k. As a result, while EKE is large at the jet latitude, larger values in eddy diffusivity or waviness are found in the polar regions, associated with the large effective time scale that is a combination of 6)] and the best fit (red) to a damped oscillator [see Eq. ( 7)]. Estimates of intrinsic frequency v i , decay time scale t, and effective eddy time scale t e for the damped oscillator are listed at the top of each subpanel. The four selected points are (a),(b) the Atlantic jet at 388N, 748W, (c),(d) Pacific jet at 338N, 1438E, (e),(f) Greenland at 748N, 778W, and (g),(h) Alaska at 588N, 1488E. These locations are marked by * in Fig. 3. large decay time scale t and small frequency v i . Thus, we conclude that the increased jet meandering over Alaska or Greenland is related to enhanced eddy persistence rather than eddy intensity.

More importantly, the eddy mixing theory can quantify the sensitivity of eddy time scale to changes in time mean zonal wind speed. Using Greenland as an example, the autocorrelation function along the direction of increasing zonal wind speed exhibits a robust decrease in eddy persistence (from black to green lines in Figs. 4a,b). This is depicted schematically in Fig. 1b that the direction of the faster mean zonal wind will be further away from the direction of eddy phase speed. This is consistent with the notion that blocking or wave breaking is less likely to occur in a stronger background zonal flow. This sensitivity agrees with the linear theory in Eq. ( 1), suggesting that an increase in u results in a decrease in effective time scale. Moreover, we have quantified the sensitivity parameter as dt e /du 5 dt e /du, where dt e is the response of eddy time scale to a change of zonal wind by du for the 2D eddy correlation map at each location, as exemplified in Figs. 4a and 4b. Specifically, we use a range of du from 1 to 10 m s 21 with an increment of 1 m s 21 , and dt e /du is obtained from the least squares fit over the range of dt e in response to zonal wind changes. This represents the average sensitivity for the range of increases in zonal wind from 1 to 10 m s 21 . Figure 4c shows that dt e /du is generally negative, indicating that the composites of eddies with faster zonal phase speeds tend to be less persistent in time. Interestingly, the sensitivity dt e /du exhibits a high spatial resemblance to t e , indicating the regions with large eddy time scales are more sensitive to changes in the time mean zonal wind. b) LWA Z500 (shading; 310 8 m 2 ) and Z500 (black contour; CI: 100 m). (c) Intrinsic frequency v i (shading; day 21 ) and zonal wind (contour; CI: 10 m s 21 ). (d) Effective eddy time scale t e 5 t/(1 1 v 2 i t 2 ) (contour; CI: 0.5 days) and eddy diffusivity D L (shading; 310 6 m 2 s 21 ). The damped oscillator parameters (i.e., v i , t, and t e ) are estimated from local autocorrelation functions as in Fig. 2.

The spatial sensitivity of LWA to zonal wind perturbations is demonstrated in Fig. 5 using a tracer model that is driven by reanalysis winds. Details of the tracer model are described in Chen et al. (2022), and the sensitivity is also shown in Fig. S5 in Nie et al. (2023). The idealized zonal wind perturbations have a Gaussian shape with a maximum speed of 5 m s 21 . For the zonal wind perturbations that decelerate the jet streams over the North Pacific or North Atlantic, the response of LWA is nearly zero. In contrast, a large increase in LWA is seen for the slowdown of the westerly wind at the poleward flank of the jets. These results agree with the spatial sensitivity of dt e /du in Fig. 4d. We have also predicted the LWA response to the zonal wind perturbations directly from Eq. ( 9). The similarity between the LWA responses from the tracer model (Figs. 5e-h) and eddy mixing theory (Figs. 5i-l) corroborate the linear theory, albeit the linear assumptions made in deriving eddy diffusivity.

While the signals downstream of the changes in zonal wind are consistent between the eddy mixing theory and tracer model (Figs. 4h,i), the upstream signals differ, likely due to the assumptions in the linear theory. Moreover, we have focused on the symmetric response to positive and negative zonal wind anomalies. The asymmetric response to zonal wind deceleration may be relevant to the abrupt onset of blocking, although this is out of the scope of the linear theory in this study.

Explaining the interannual variability and climate response of waviness

Given the spatial sensitivity of dt e /du, we can estimate the changes in LWA to different zonal wind anomalies using Eq. ( 9). Here, we will briefly present (i) the interannual variability of waviness associated with the negative phase of the Arctic Oscillation (AO) in the reanalysis examined in Chen et al. ( 2022) and (ii) the response of waviness to climate change in large-ensemble climate simulations analyzed in Nie et al. (2023). The comparison between the eddy mixing theory and the tracer model helps assess the accuracy of the eddy mixing theory.

a. Interannual variability of waviness associated with the AO

Figure 6 shows the composite spatial distributions of DJF anomalies of eddy activity in the negative phase of the AO. Consistent with the equatorward shift of the North Atlantic jet, the Atlantic storm tracks are displaced southward (Fig. 6a). In contrast, LWA displays a dipolar pattern of opposite sign to zonal wind changes over the North Atlantic (Fig. 6b), consistent with the slower and blocked zonal flow observed over the Atlantic Ocean in the negative phase of the AO (Woollings et al. 2008;Cattiaux et al. 2016). A similar anticorrelation between zonal wind and LWA is found over East Siberia/Alaska. Chen et al. (2022) found that the anomalous LWA in the negative phase of AO can be roughly simulated in a tracer model with prescribed changes in advecting winds, which is displayed in Fig. 6c. Here, we further show that the anomalous LWA pattern over the North Atlantic can be largely explained by changes in eddy mixing and to a lesser extent over the North Pacific (Fig. 6d). Moreover, Eq. ( 9) indicates that the changes in waviness can be separated into changes due to effective time scale and changes due to EKE. The changes in effective time scale are, in turn, anticorrelated with the changes in zonal wind, due to the negative sign of dt e /du. Similar to the tracer model simulations in Chen et al. (2022), the diffusivity-scaling-based estimate suggests that effective time scale due to zonal wind changes dominate the effects on LWA over EKE (Figs. 6e,f).

The domain-averaged time series of LWA are evaluated over Greenland or Alaska regions. The tracer model can simulate the interannual variability in LWA, with a higher correlation over Greenland (r 5 0.77) than Alaska (r 5 0.56) (Figs. 7a,b,e,f). Note that the correlation coefficients with LWA based on Z 500 here are slightly lower than those based on the 250-hPa PV in Fig. 8 of Chen et al. (2022), as the tracer driven by 250-hPa winds is expected to be more correlated with the 250-hPa PV rather Z 500 . These temporal variations in LWA are, to a large extent, captured by the predictions from the eddy mixing theory. The correlation coefficients indicate that about 60% of the LWA variance over Greenland and 30% of the variance over Alaska are explained by the respective variance in effective time scale due to zonal wind changes (Figs. 7c,d,g,h).

b. Climate response in waviness in climate simulations

We next switch to the climate response in waviness in largeensemble climate simulations. Figure 8 displays the future changes in LWA between the historical period and the future period (2080-99), using initial-condition large ensembles from four climate models (i.e., CESM1, CanESM2, GFDL CM3, and CESM2). As noted by Nie et al. (2023), the four large ensembles disagree on the future changes in waviness over the North Pacific sector at the end of the twenty-first century, with a decrease in GFDL CM3 and CESM2 and an increase in CanESM2, but the models predict a consistent weakening in LWA over the North Atlantic sector. Using the climate response in zonal wind, the different signs and magnitudes of future changes over the Pacific and Atlantic sectors are largely reproduced by both the advection-diffusion model and eddy mixing theory. The pattern correlations between the LWA results are r 5 0.43, 0.7, 0.54, and 0.62 for tracer calculations and r 5 0.52, 0.74, 0.54, and 0.66 for diffusivity scaling.

Given that changes in EKE are not incorporated in either tracer calculation or diffusivity scaling, this indicates anomalous zonal wind advection alone can capture the large-scale feature of the varied model response in LWA to climate change. As discussed in Nie et al. (2023), these changes may be attributed to a consistent poleward shift in the North Pacific jet (with the exception of CanESM2) and a robust extension of the North Atlantic jet into Europe. From Fig. 3 in Harvey et al. (2020), those regions correspond to consistent increases in EKE, which indicates that the changes in EKE may counteract the changes in effective time scale and result in a smaller net change in LWA.

Finally, the estimates from either the tracer model or eddy mixing theory tend to be more zonally oriented while LWA in large ensembles is more regional. This may be attributed to a dynamical inconsistency in eddy-mean flow interactions and jet meandering (Haynes et al. 2007;Valva and Nakamura 2021), as neither the passive tracer nor diffusivity scaling has included the feedback on transient eddies (e.g., EKE).

Causality for circulation waviness in kinematic and dynamically consistent flow

In this section, we reiterate the causality discussed in the paper based on kinematic models, in which passive tracers are advected by prescribed velocities, and after a spin-up period, the geometry of tracer field is only a function of prescribed advecting velocities (Chen et al. 2022;Nie et al. 2023). This is motivated by previous studies in which the tracer model is advected by observed winds and initialized with a zonally symmetric passive tracer increasing with latitude like PV, and after being integrated for about a month, the geometric structure of the tracer captures many fine-scale features of PV that are independent of initial tracer profiles (e.g., Haynes and Shuckburgh 2000;Allen and Nakamura 2001;Abalos et al. 2016). In this kinematic framework, the waviness of passive tracers or eddy mixing results solely from the property of the advecting flow (e.g., EKE, the intrinsic frequency of Rossby waves or their decay rate). We further argue that the behavior of PV contours, in the absence of diabatic heating or friction, can be largely understood through this PV-like tracer advected by the observed winds. Specifically, tracer and PV contours corresponding to the same equivalent latitude display similar geometric structures, which may be quantified by a transformation from tracer to PV. As the advecting velocities of the tracer can be directly modified in the advection-diffusion model, we can evaluate the effect of the time mean flow on the geometry of a PV-like tracer or the causality between the two in a modeling framework.

In a dynamically consistent flow, PV is related to the advecting flow through the invertibility principle (e.g., Hoskins et al. 1985), and thus, PV actively modulates the advecting flow. Large PV gradients are typically associated with small PV excursions; small PV gradients allow the PV contours to meander strongly, yet the kinematic model does not directly consider this dynamical feedback of PV on eddy mixing (e.g., Haynes et al. 2007). While the kinematic perspective can describe the causal relationship between the advecting flow and the waviness of passive tracers, additional feedback should be considered for the two-way interactions between PV and the advecting flow. One possible approach is to extend the tracer model to a global circulation model in which the horizontal winds are nudged to reanalysis winds. In this type of nudging experiment, the simulated horizontal winds are not necessarily consistent with the PV invertibility, but they offer dynamical insights into the two-way coupling between the mean flow and planetary waves (e.g., Hitchcock and Simpson 2014;Ding et al. 2023). The analyses in section 5 may be further examined by exploring different nudging strategies for the impact of changes in the time mean wind on circulation and extreme events in a more realistic modeling framework.

Conclusions

Using an advection-diffusion model for PV-like passive tracers, previous studies have found that variations in circulation waviness can be largely constrained by anomalous seasonal mean zonal winds (Chen et al. 2022;Nie et al. 2023). This eddy diffusive perspective characterizes a strong zonal jet as a "mixing barrier" (e.g., Haynes and Shuckburgh 2000;Allen and Nakamura 2001;Chen and Plumb 2014;Abalos et al. 2016) or as a barrier to blocking. In this paper, we emphasize the eddy mixing time scale for passive tracers as a key measure of jet waviness rather than EKE. Based on the onepoint correlation map of anomalous meridional wind, we estimate the intrinsic frequency and decorrelation time scale of eddies for the Rossby wave packet, with the former characterizing the frequency of Rossby wave propagation and the latter characterizing the decay of wave amplitude. The combination of the two is quantified by a damped oscillator that describes the time scale of jet waviness [Eq. ( 8)].

The spatial pattern of the waviness time scale can largely capture the spatial distribution of atmospheric waviness as depicted by LWA (Fig. 3). In the high latitudes where zonal winds are weak, the long eddy mixing time scales coincide with the large values in LWA. In the subtropical regions characterized by strong westerlies, the intrinsic frequency is large, resulting in suppressed eddy mixing time scales and LWA despite the large values in EKE. Moreover, the eddy time scale can roughly capture the interannual variability of LWA in observations and the response of LWA to climate change in climate model simulations. As the intrinsic frequency is largely modulated by the time mean zonal winds, changes in intrinsic frequency provide a physical mechanism for how the zonal winds impact regional circulation waviness in a changing climate.

It is noteworthy that the time scale of circulation waviness (Fig. 3d) is shorter than the typical time scale of blocking events, but LWA encapsulates the Rossby wave events of all time scales, and this is consistent with the time scale of temperature fluctuations. Using an equator-to-pole temperature gradient dT/dy 5 6 3 10 26 K m 21 and the daily fluctuation of temperature T 5 3 K, this gives the magnitude of meridional displacement in temperature contours as h 5 T /(dT/dy) 5 5 3 10 5 m. If the magnitude of wind is estimated as y 5 5 m s 21 , the time scale of waviness in temperature is t 5 h/y ' 1 day, consistent with the estimate in Geen et al. (2023) (their Fig. 3). Future research can extend the analysis from the average of all wave events to extreme events by focusing on the Rossby wave packets during blocking events and their persistence.

computing [ ỹ (k, v) ỹ * (k, v)] v5uk requires an integral of the spectral power ỹ (k, v) ỹ * (k, v) along constant phase speed c 5 u. Alternatively, we employ the eddy covariance form in Eq. (A4), which computes eddy autocovariance in the moving coordinate at speed u.