A physical analysis of summertime North American heatwaves

This study examines the dominant heatwave variability over North America (NA), extracted from an empirical orthogonal function (EOF) analysis of summertime monthly warm extreme index anomalies over 1959–2021. The principal mode features a dipole structure with a large area of anomaly over northwestern NA and an anomaly of opposite sign over the southern U.S. The corresponding principal component is associated with a large-scale atmospheric wave train extending from the North Pacific to North America (NP-NA) and a northeastward injection of moisture from the subtropical western Pacific towards western NA, which are key factors in supporting the NA heatwave variability. The NP-NA wave train can be systematically reinforced and supported by synoptic-scale eddies, and may also be forced by an anomalous convection over the tropical-subtropical western Pacific. Surface radiation heating directly contributes to surface temperature anomalies and is dominated by anomalous downwelling shortwave and longwave radiations. In association with a positive phase of the heatwave variability, the NP-NA wave train brings an anticyclonic anomaly over northern NA, leading to anomalous descent, reduced total cloud cover and below-normal precipitation and surface relative humidity over northern NA. Over northwestern NA, the anomalous subsidence causes air to warm through compression. Reduced cloud cover results in increased downward shortwave radiation that is a key contributor to surface radiation heating. In addition, increase in vertically integrated water vapour through the moisture injection from the North Pacific collocates with tropospheric warming. The atmosphere has more water vapor holding capability and acts as a greenhouse gas to absorb longwave radiation, leading to increased downward longwave radiation that is the second major contributor to surface radiation heating. Processes with circulation and surface radiation anomalies of opposite signs will likewise lead to the negative heatwave variability.


Introduction
Climatic extremes have significant impacts on natural systems and human life and have been extensively studied, including characterizing internal/natural and projected changes in climate extremes and exploring their physical mechanisms (e.g., Meehl and Tebaldi 2004;van Oldenborgh et al. 2015;Sillmann et al. 2017;IPCC 2013IPCC , 2021. Changes in atmospheric circulation are responsible for changes in temperature, moisture, and precipitation, which significantly affect regional climate and are associated with extreme weather conditions leading to severe economic and social damages. Therefore, extreme temperature changes such as severe summer heatwaves and cold winters have been attributed to changes in regional and large-scale atmospheric circulations (e.g., Coumou et al. 2012;Wu et al. 2012;Horton et al. 2015;Wolter et al. 2015;Grotjahn et al. 2016;Harnik et al. 2016;Müller et al. 2020;Suarez-Gutierrez et al. 2020), especially atmospheric teleconnections (e.g., Teng et al. 2013;Johnson et al 2018;Yu et al. 2018Yu et al. , 2019Yu et al. , 2020, and regional-and global-scale thermodynamic processes (e.g., Campbell and Vonder Haar 1997;Durre and Wallace 2001;Loikith and Broccoli 2012;Kamae et al. 2014;Pithan and Mauritsen 2014;Holmes et al. 2015;Schneider et al. 2015;Yu and Zhang 2015;Coumou et al. 2018;Gross et al. 2020). Extreme temperature anomalies have also been related to 1 3 high frequency climate activities, such as cyclone/anticyclone tracks and blocking activity (e.g., Favre and Gershunov 2006;Parker et al. 2014;Chang et al. 2016;Sillmann et al. 2017;Kirshbaum et al. 2017;Schaller et al. 2018;Yu et al. 2022).
Focusing on recent extreme heatwaves, record-breaking above-normal temperatures were observed in western North America (NA) from late June to early July of 2021. Lin et al. (2022) and Mo et al. (2022) explored the physical mechanism and subseasonal prediction of this extreme heatwave event. A northward propagation of the Southeast Asian summer monsoon (SEASM) in late June of 2021 was found to cause precipitation and convection anomalies in the subtropical western Pacific. The precipitation anomaly acted as an anomalous diabatic heating to trigger an extratropical Rossby wave train propagating from the western Pacific, through the North Pacific towards North America, leading to synoptic conditions favorable for the extreme heatwave in western NA. In particular, the anomalous wave train brought a strong high-pressure known as heat dome (NOAA 2022) over western NA, under which the warm air was trapped and the associated subsidence caused further warming by compression. Meanwhile, there was a northeastward advance of a strong atmospheric river (AR, Ralph et al. 2004Ralph et al. , 2018 across the North Pacific towards western NA. Strong ARs can pump moisture out of tropical regions and incorporate mid-latitude moisture sources along their paths (Dettinger et al. 2015;Eiras-Barca et al. 2017;Mo et al. 2022). The AR injected moisture in western NA could act as a greenhouse gas to trap solar radiation that further warmed the lower atmosphere. Through the mechanisms of the atmospheric teleconnection and AR contribution, an extreme heatwave occurred across western NA in June-July, 2021 Mo et al. 2022).
Does this particular extreme heatwave in western NA indicate an intrinsic pattern of NA heatwave variability? Are the atmospheric teleconnection and ARs from the North Pacific to North America key factors supporting historical extreme heatwaves in NA? Furthermore, what are the physical processes responsible for the large-scale atmospheric circulation and, more directly, surface temperature anomalies associated with historical heatwaves in NA? These issues are of fundamental importance to North American climate variability and prediction. In this study, we address these questions from a climatological perspective, by analyzing the variability of monthly extreme surface temperatures and the associated atmospheric and oceanic fields over the past 60 years.
In the following section, we outline the data and methods employed in this study. We then examine the leading summertime heatwave variability over NA in Sect. 3. Section 4 presents physical diagnostics of the leading heatwave variability, including the heatwave associated anomalies of large-scale circulation and moisture; we also analyze the maintenance of the circulation and surface energy budget anomalies to help understand the generation and maintenance of NA heatwaves. Results are summarized and further discussed in Sect. 5.

Data and methodology
(a) Reanalysis data, observations, and climate indices This study is mainly based on June-September (JJAS) atmospheric variables extracted from the fifth generation of atmospheric reanalysis (ERA5, Hersbach et al. 2020) of the European Centre for Medium-Range Weather Forecasts (ECMWF) over the period from 1959 to 2021. We use hourly surface air temperature from ERA5 to create a warm extreme index (TX90) with the percentage of time when daily maximum temperature is above its 90th percentile. TX90 is derived using 1981-2010 as the base period and applying a 5-day running window, following Donat et al. (2013). This extreme index is then analyzed on a monthly basis (e.g., Sillmann et al. 2013;Yu et al. 2020). We use 6-hourly wind data to calculate the monthly synoptic eddy forcing defined in Eq. (2) below. In addition, the monthly ERA5 reanalysis fields we employed include mean sea level pressure (MSLP), heat fluxes at the surface, air temperature, geopotential and wind velocities in the troposphere, surface relative humidity, precipitation, total column vertically-integrated water vapour (TCWV), and total cloud cover (TCC). These variables are analyzed on 2.5° × 2.5° grids. Note that TCWV is the same as integrated water vapour (IWV) used in AR analysis, or precipitable water (PW) used in weather forecasting. It is closely correlated to the vertically-integrated horizontal water vapor fluxes, and therefore is often used as a proxy to determine the position and width of ARs (Ralph et al. 2004).
We also use observed monthly TX90 index extracted from the updated Hadley Centre extreme (HadEX2) dataset (Donat et al. 2013) from 1959 to 2010. The observed extreme index is available over land and on 2.5° × 3.75° (latitude-longitude) grids. In addition, we employ monthly sea surface temperature (SST) over 1959-2021, obtained from the Extended Reconstructed Sea Surface Temperature version 5 (ERSSTv5) dataset (Huang et al. 2017). SST data are interpolated from 2.0° × 2.0° grids to 2.5° × 2.5° grids using a bilinear interpolation. We also use monthly outgoing longwave radiation (OLR) from NCAR archives (https:// www. ncei. noaa. gov/ produ cts/ clima te-data-recor ds/ outgo ing-longw ave-radia tion-month ly), with gaps filled with a temporal and spatial interpolation (Liebmann and Smith 1996), on 2.5° × 2.5° grids over 1979-2021. Atmospheric pattern indices of the Pacific-North American (PNA, Wallace and Gutzler 1981), East Pacific-North Pacific (EP-NP, Bell and Janowiak 1995), West Pacific (WP, Wallace and Gutzler 1981), North Atlantic Oscillation (NAO, e.g., Hurrell et al. 2003), and Arctic Oscillation (AO, Thompson and Wallace 1998) from 1959 to 2021 are obtained from the Climate Prediction Center (CPC, http:// www. cpc. ncep. noaa. gov/ data/ indic es). Except AO, these patterns are identified from a rotated empirical orthogonal function (REOF) analysis of monthly mean normalized 500-hPa height anomalies over the northern extratropics (poleward of 20°N), based on the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis (Kistler et al. 2001). The AO is identified from an empirical orthogonal function (EOF) analysis of monthly mean 1000-hPa height anomalies over the northern extratropics. The Asian-Bering-North American (ABNA, Yu et al. 2016) pattern index is also calculated based on the ERA5 monthly data. This index is constructed by an additive combination of three regionally averaged geopotential anomalies over the action centers of ABNA (North Asia, Bering Sea and Strait, and North America), using the normalized 500-hPa geopotential field after linearly removing the PNA pattern contribution (Yu et al. 2016(Yu et al. , 2018. In addition, the monthly oceanic Niño3.4 index from CPC is used to characterize the tropical El Niño-Southern Oscillation (ENSO) variability.

(b) Data processing and diagnostic methods
We use JJAS monthly data to examine the variability and formation of NA heatwaves. Monthly anomalies are calculated relative to the 30-year respective monthly mean climatology over 1981-2010 after removing the secular linear trend over the whole period considered. We perform an EOF analysis to identify the principal mode of area-weighted monthly TX90 anomalies over NA. The corresponding leading principal component (PC1) is then employed as an index of the heatwave variability. Relationships between the PC1 index and variables of interest are quantified through correlation and regression analyses. The statistical significance of a correlation is determined by a Student's t test, where the effective degree of freedom is estimated by considering the autocorrelation of the time series analyzed (Bretherton et al. 1999).
To explore the maintenance and driving mechanisms of the circulation anomalies associated with the heatwave variability, we examine the PC1 associated dynamical anomalies of wave activity flux and transient eddy forcing in the upper troposphere. The horizontal wave activity flux ( WAF ) vector is calculated following Takaya and Nakamura (2001) as follows: where g = (u g , v g ) and are horizontal geostrophic winds and stream function, respectively, which are calculated from the geopotential field. The double overlines and prime denote the 30-year JJAS climatological mean over 1981-2010 and its anomaly, respectively. The subscripts represent partial derivatives. WAF tends to be parallel to the local group velocity of stationary Rossby wave, which illustrates sources and sinks of wave activity and reveals the atmospheric wave dispersion (Takaya and Nakamura 2001).
In the upper troposphere, synoptic-scale eddy forcing is crucial in reinforcing the anomalous circulation (e.g., Lau 1988;Kug and Jin 2009). The synoptic-scale vorticity forcing (F v ), representing the geopotential tendency, can be expressed as, where f is the Coriolis parameter, V(u,v) the wind vector including zonal and meridional components, and the relative vorticity, following standard conventions. The asterisk denotes a 2-8 day bandpass filtered synoptic perturbation and the overbar indicates a monthly average.
We also examine the PC1 associated thermodynamic anomalies of temperature advection and surface energy budget to aid in understanding the formation mechanism of surface temperature anomalies. The anomalous horizontal temperature advection (F adv ) in the lower troposphere, dominated by wind variations, can be expressed as, where V' is the anomalous horizontal wind velocity, and T denotes the 30-year JJAS climatological mean temperature over 1981-2010.
The surface energy budget over land can be written as (e.g., Yu and Boer 2002;Zhang et al. 2011), where C s is the surface layer heat capacity, T s the surface skin temperature, R d s ( R u s ) the surface downward (upward) shortwave radiation, R d l ( R u l ) the surface downward (upward) longwave radiation, LH (SH) the latent (sensible) heat flux at the surface, F g the ground energy flux, and F m the surface heat flux used for melt. R = R d s + R u s + R d l + R u l denotes the surface net radiation, and B = LH + SH is the surface turbulent heat flux. In addition, the surface radiation heating is defined as the sum of surface net shortwave radiation and downward longwave radiation, which indicates how cloud feedback and water vapour contribute to surface temperature anomalies (Zhang et al. 2011). All fluxes are defined to be positive downward to warm the surface. For monthly and JJAS averages, the left-hand side term of Eq. (4) is much smaller than heat flux terms and is not considered (e.g., Yu and Zhang 2015). In addition, F g and F m are smaller than other heat flux terms and are also not discussed.

The leading heatwave variability over North America
The principal mode (EOF1) of JJAS monthly TX90 anomalies over 1959-2021 from ERA5 accounts for 25.6% of the total variance over NA and is well separated from subsequent EOFs according to the criterion of North et al. (1982). EOF1 is dominated by a dipole pattern of the TX90 anomaly featuring a large area of anomalies in northwestern NA, with the center of action located over British Columbia, Alberta and Washington State, accompanied by slightly stronger opposite-sign anomalies in the southern U.S., with the action center over Arkansas, Louisiana and Mississippi (Fig. 1,. Similar results can be obtained using the observational data. EOF1 from the observed monthly TX90 anomalies over 1959-2010 accounts for 20.4% of the total variance ( Fig. 1, top-right) and is also separated from subsequent EOFs. EOF1 patterns from the two datasets bear close spatial resemblance, with differences mainly seen in the magnitude of action centers. The corresponding principal components (PC1s) are then used to characterize the variability of the leading warm extreme pattern. Variations of PC1s from the observation and ERA5 are also similar in the overlapping period ( Fig. 1, bottom), especially after 1979 probably due to inclusion of satellite data post 1979 in the reanalysis. The two PC1s are highly correlated, with a correlation coefficient of 0.92 (0.95) for the monthly series over 1959-2010 (1979-2010). This implies that the two PC1 series have about 85% (90%) variance in common over the 52 (32) summers. Similarities of both spatial and temporal  1959-2010 (blue) patterns in the observation and reanalysis suggest that ERA5 can well represent the dominant TX90 variability over NA. The TX90 anomaly over the northwestern part of NA in EOF1 (Fig. 1, top panels) resembles the record-breaking heatwave distribution over western NA in June-July 2021 ( Fig. 1 in Lin et al. 2022). The PC1 index from ERA5 shows pronounced anomalies of 1.64 in June and 2.83 in July 2021, the latter of which is the maximum over 1959-2021 (Fig. 1,  bottom). By contrast, the predominantly positive TX90 anomalies in negative PC1 extreme months, such as July 1980 and June-August 2010, occur in the southern U.S. (not shown). In addition, the amplitude of PC1 increases in recent summers. The standard deviations of the PC1 index are 0.76 over 1959-2000 and 1.14 over 2001-2021, reflecting a stronger fluctuation of the dipole pattern in recent summers. For example, high PC1 values in the summers of 2007 and 2009 also correspond to remarkable heatwaves recorded over northwestern NA (e.g., Goodrich et al. 2011). Hence, the PC1 from ERA5 that includes the up-to-date summer series over 1959-2021 is used as an index of the heatwave variability over NA in the following calculation.

Physical diagnostics of the heatwave variability
To explore the generation and maintenance of the NA heatwave variability, we examine the PC1 associated anomalies of large-scale atmospheric circulation, surface energy balance, and relevant dynamical and thermodynamic processes.
(a) Large-scale circulation and moisture anomalies We first calculate correlations between the PC1 index and several indices of prominent modes of climate variability in the Northern Hemisphere, which are associated with monthly/seasonal mean surface temperature anomalies over NA as demonstrated in previous studies (e.g., Higgins et al. 2002;Yu et al. 2019). PC1 correlates with the NAO and AO indices, with correlation coefficients of 0.34 and 0.24 over the 63 summers, respectively, statistically significant at the 5% level (Table 1). This indicates that the PC1 associated atmospheric circulation anomalies bear some resemblance to the NAO and AO patterns. In association with a positive phase of the PC1 index, the 500-hPa geopotential (Φ 500 ) and MSLP anomalies reveal below-normal geopotentials and pressure across the high latitude North Atlantic and above-normal anomalies over the central North Atlantic (Fig. 2), which are signatures of the NAO pattern and the AO pattern over the North Atlantic sector (e.g., Hurrell et al. 2003;Thompson and Wallace 1998). PC1 is also marginally correlated with the EP-NP index at the 5% level (Table 1), which is likely due to the significant anticyclonic circulation anomaly over western NA (Fig. 2) that is one of three action centers of the EP-NP pattern (e.g., Bell and Janowiak 1995). By contrast, the correlations between the PC1 and PNA, ABNA, WP and Niño3.4 indices are low ( Table 1).
The PC1 associated Φ 500 anomalies in the northern midlatitudes are dominated by a North Pacific to North American (NP-NA) wave train, featured by a sequence of high and low height centers arcing from the subtropical western Pacific, across the North Pacific and north-central NA, and then to the southeastern U.S. (Fig. 2, top). In particular, in association with a positive phase of the PC1 index, a pronounced anticyclone anomaly can be seen in northern and central NA along the wave train path. The center of action of this anticyclone anomaly is located over British Columbia, Alberta and Washington, collocating with the corresponding center of the EOF1 pattern of TX90 anomalies over NA (Fig. 1, top-left). In addition, the NP-NA wave train is clearly evident in both Φ 500 and MSLP anomalies associated with the PC1 index, with differences mainly in magnitude of the anomalies, indicating an equivalent barotropic structure of the large-scale circulation anomalies in the troposphere. Figure 3 displays the PC1 associated precipitation and total column vertically-integrated water vapour anomalies. In association with a positive phase of the PC1 index, significant precipitation increases appear over the subtropical western Pacific, the Sea of Okhotsk, southeastern Alaska and the northeastern Pacific, the southern U.S., and southern Greenland. By contrast, significant precipitation decreases are mainly seen over South Asia, the tropical western Pacific, the Yellow Sea, the central North Pacific, central NA, and parts of the tropical and central North Atlantic (Fig. 3, top). The reduced precipitation and convection in South Asiawestern Pacific indicates weakened Southeast Asian summer monsoon activity. This represents an anomalously negative diabatic heating in the tropical region. The TCWV anomaly is broadly similar to the precipitation anomaly over the Table 1 Correlation coefficients between the PC1 index and atmospheric and oceanic indices, based on JJAS monthly series over 1959-2021 The effective sample size is found to be at least 165 for the 252-month time series in these calculations. The 5% significance level of a correlation is r > 0. 15  northern Pacific (Fig. 3, bottom). The area-weighted pattern correlation between the two anomalous fields is 0.61 over the northern (0-90 o N) Pacific. In addition, it is interesting to see an anomalous TCWV path across the mid-latitude North Pacific, indicating a northeastward injection of moisture from the subtropical western Pacific towards western NA. By contrast, the TCWV anomalies are mostly positive over western NA and negative over eastern NA, which differs from the precipitation anomalies over NA. In particular, many areas of increased (reduced) TCWV collocate with reduced (increased) precipitation over the NA continent (Fig. 3). This suggests that the anomalous precipitation is not only related to the water vapor anomaly, but also to other mechanisms such as those responsible for inducing adiabatic vertical motions and moisture convergence (e.g., Zhang et al. 2011;Mo et al. 2021). Previous observations and climate model results also suggested that deficits in soil moisture would lead to more frequent and severe hot temperatures in summer, and vice versa. The above-normal (below-normal) precipitation often increases (decreases) the soil moisture, which in turn favors less (more) heatwaves (e.g., Dai et al. 1998;Alexander 2010;Wu et al. 2012). In addition, both precipitation and TCWV anomalies are weak upstream over most of the Eurasian continent. The strongest water vapour transport (AR-like feature) often penetrates through a cyclone-anticyclone couplet. The Rossby wave train consists of alternative cyclones and anticyclones, which forms a favorable weather pattern for the northeastward moisture transport from the subtropical western Pacific towards the west coast of NA (e.g., Mo et al. 2022). Hence, the NP-NA wave train pattern and the northeastward transport of moisture are correlated. To quantitatively analyze the linkage between the TCWV anomaly pattern and the NP-NA teleconnection, we have projected the Φ 500 anomalies onto the NP-NA pattern (Fig. 2, top) over (20-80° N, 120° E-90° W) to create a NP-NA pattern index, and projected the TCWV anomalies onto the moisture pattern (Fig. 3, bottom) over the same domain as a TCWV pattern index. The two indices obtained are not sensitive to a slight change of the projected domain. Figure 4 shows the normalized JJAS monthly time series of the projected Φ 500 and TCWV pattern indices, superimposed on the PC1 index, over 1959-2021. The correlation between PC1 and the NP-NA (TCWV) pattern index is 0.71 (0.64), implying Black dots indicate anomalies that are significantly correlated with the PC1 index at the 5% level the good relations between the NA heatwave variability and the two patterns. The correlation between the NP-NA and TCWV indices is 0.75, significant at the 1% level. It also indicates that the two series have over 55% variance in common. In addition, most of the extreme values in the NP-NA index are well corresponding to the extremes in the TCWV pattern index (Fig. 4), indicating that the significant correlation between them is not due to a few extreme events.
The PC1 associated atmospheric circulation and water vapor anomalies across the North Pacific towards NA, as well as the precipitation anomalies over NA, bear resemblance to those occurred in the extreme heatwave event in June-July 2021 (see Figs. 1 and 4 in Lin et al. 2022, andFigs. 1, 2 in Mo et al. 2022). This indicates that the NP-NA atmospheric teleconnection and moisture anomalies can be seen in both the 2021 heatwave event and variations of NA heatwaves from a climatological point of view, and suggests that they are key factors in supporting the NA heatwave. Nevertheless, it is worth noting that the geopotential anomaly associated with the PC1 index covers a wider area of north-central NA compared to the 2021 heatwave event. As will be discussed below in Sect. 4.c, a combination of above-normal (below-normal) temperature and below-normal (above-normal) precipitation tends to be a consequence of a positive (negative) geopotential anomaly over NA.
(b) Maintenance of the circulation anomalies Similar geopotential anomalies associated with the PC1 index are apparent at 200-hPa and 500-hPa (c.f. contours in the top panel of Fig. 5 with Fig. 2), with differences mainly in magnitude of the anomalies. Diagnosis of the wave activity flux (WAF, Eq. (1)) for stationary Rossby waves (vectors in the top panel of Fig. 5) reveals the propagation of the NP-NA wave train. The WAF analysis also suggests that the main source of wave activity associated with the NA heatwave variability may be traced back to the western North Pacific.
The centers of action of the NP-NA wave train are located near the mid-latitude jet stream and storm tracks, suggesting that the generation of this wave train may be related to local dynamical processes such as transient eddy variations and jet instabilities. Interaction of synoptic eddies and low-frequency airflow is critical for the formation of atmospheric circulation variability in mid-latitudes (e.g., Trenberth and Hurrell 1994;Branstator 1995;Kug and Jin 2009). In particular, synoptic vorticity forcing can reinforce anomalous circulation in the upper troposphere (e.g., Lau 1988; Kug and Jin 2009;Choi et al. 2010;Yu and Lin 2013). Figure 5 (bottom) displays the geopotential and synoptic eddy vorticity forcing (F v , defined in Eq. (2)) anomalies at 200-hPa regressed upon the PC1 index. In association with the NA heatwave variability, dominant anticyclonic circulation forcings can be seen over the northwestern Pacific, NA, and the mid-latitude North Atlantic, which coincide well with the corresponding positive Φ 200 anomalies. Meanwhile, dominant cyclonic forcings are apparent over northeastern Asia, the northeastern Pacific, and the high-latitude North Atlantic that collocate with the corresponding negative Φ 200 anomalies. The collocation of the anomalous anticyclonic (cyclonic) vorticity forcing with the positive (negative) geopotential anomaly indicates that synoptic eddies are systematically reinforcing and supporting the anomalous wave train extending from the North Pacific towards NA.
The tropical-subtropical heating anomaly may act as a forcing source for atmospheric circulation anomalies. Previous studies indicated that the atmospheric circulation is especially sensitive to perturbations in the tropical-subtropical western Pacific (e.g., Simmons et al. 1983;Palmer and Mansfield 1984;Trenberth et al. 1998;Barsugli et al. 2002;Lin 2009). Significant precipitation anomalies associated with the heatwave variability appear in the tropical-subtropical western Pacific (Fig. 3, top), which is also evident in the PC1 associated OLR anomalies over 1979-2021 (not shown). This indicates that an anomalous convection over the tropical-subtropical western Pacific may act as a forcing on the NP-NA wave train. Influences of the precipitationinduced cooling anomaly in the Philippine Sea and heating anomaly in the subtropical SEASM region on the NP-NA teleconnection have been confirmed in numerical experiments by Lin et al. (2022) and Qian et al. (2022). On the other hand, the PC1 associated SST anomalies are weak in the tropics (Fig. 6, top), consistent with the low correlation between the PC1 and Niño3.4 indices noted above. SST is dominated by anomalies in the mid-latitude North Pacific and the tropical-subtropical North Atlantic. Strong low-level wind anomalies associated the PC1 index are also seen in the North Pacific (Fig. 6, bottom and Fig. 2, bottom) and are well collocated with the corresponding action centers of SST anomalies. This implies that the anomalous atmospheric circulation largely drives the SST in the mid-latitude North Pacific, in agreement with previous studies (e.g., Cayan 1992;Kushnir et al. 2002). Figure 6 (bottom) further reveals that the thermal advection (F adv , Eq. (3)) anomalies at 850-hPa in association with the PC1 index are weak over land, especially over northwestern NA. The weak thermal advection anomaly in the lower troposphere would be expected in summer. It also suggests that surface temperature anomalies over NA in summer are likely determined by local thermodynamic processes, through surface heat flux anomalies.
(c) Surface heat flux anomalies. Figure 7 displays the anomalies of surface energy budget (Eq. (4)) over NA associated with the PC1 index. Table 2 lists the spatial standard deviations (< X + > 1/2 ) of the anomalous TX90 (EOF1) and anomalies of various terms in the energy budget calculated over NA, as well as the spatial covariances (< TX90 + , X + >) and correlations (r(TX90 + , X + )) between the anomalies of TX90 and surface energy balance components. Here, we define X = < X > + X + with In association with the NA heatwave variability, the anomalous radiative heat flux is dominated by surface shortwave downward ( R d s ) and longwave upward ( R u l ) radiation anomalies (Fig. 7), which is apparent in the comparison of both the spatial standard deviations of various radiation components over NA and the spatial covariances between the anomalies of TX90 and surface radiation terms (Table 2). R d s ( R u l ) anomalies are highly and positively (negatively) correlated with TX90 anomalies, with a pattern correlation of 0.86 (-0.95) over NA. This indicates that the shortwave downward and longwave upward radiation anomalies are structurally similar to the TX90 anomaly. However, the anomalous R d s tends to support and maintain the TX90 anomaly, while the anomalous R u l damps the TX90 anomaly. Notably, R u l only represents the radiative response to anomalous surface temperatures following the Stefan-Boltzmann law (e.g., Peixoto and Oort 1992). Surface longwave downward ( R d l ) anomalies are weaker than R d s and R u l anomalies, but still contribute to supporting the TX90 anomaly (Fig. 7). The magnitude of the spatial covariance between the anomalies of TX90 and R d l (5.18%Wm −2 ) is about half of those between the  anomalies of TX90 and R d s (11.79%Wm −2 ) and between TX90 and R u l (−10.27%Wm −2 ) ( Table 2). In addition, surface shortwave upward ( R u s ) anomalies are much weaker compared to other radiation anomalies (Fig. 7), and are only weakly and negatively correlated with the TX90 anomaly ( Table 2).
The anomalous surface turbulent heat flux (Fig. 7) is dominated by the sensible heat flux (SH) anomaly, which shows a higher spatial standard deviation and covariance with the TX90 anomaly than the latent heat flux (LH). SH anomalies tend to damp the TX90 anomaly over NA (Fig. 7 and Table 2). The magnitude of the spatial covariance between the anomalies of TX90 and SH (-5.72%Wm −2 ) is about half of the strongest covariance between TX90 and R d s (Table 2). In addition, the anomalies of ground heat flux F g and F m are small and negligible (not shown). Overall, the surface net radiation anomaly (R, Fig. 7) supports the TX90 anomaly (Table 2) and is primarily compensated by the surface turbulent heat flux anomaly (B, Fig. 7).

(d) Contributors of the surface radiation heating
Surface radiation heating directly contributes to surface temperature changes. In association with the PC1 index, Q s is dominated by downwelling shortwave and longwave radiation ( R d s and R d l ) anomalies, with little contributions from the anomalous upwelling shortwave radiation R u s . Over NA, large Q s anomalies collocate well with large TX90 anomalies (cf. Fig. 7, bottom-left with Fig. 1 top-left). The spatial covariance between the anomalies of TX90 and Q s is 15.94%Wm −2 and the pattern correlation between them is high of 0.96 (Table 2), indicating that the TX90 anomaly over NA is largely determined by surface radiation heating.
The atmospheric circulation anomaly would modify the hydrological cycle, including the radiative effects of water vapor, clouds, and changes in surface properties such as soil moisture (e.g., Zhang et al. 2011). In association with a positive phase of the PC1 index, the NP-NA wave train brings positive (negative) Φ 500 anomalies over the northern (southern) part of NA, with the corresponding action centers over northwestern NA and the southeastern U.S. The circulation anomalies result in anomalous descent (positive 500-hPa vertical velocity ω 500 anomalies), reduced total cloud cover and below-normal precipitation and surface relative humidity (RH 1000 , relative humidity at 1000-hPa) over northern NA, and opposite anomalies over southern NA (Fig. 8). The pattern correlation between the anomalies of Φ 500 and TCC (precipitation) is -0.76 (-0.65) over NA (Table 3). Over northwestern NA, the anomalous subsidence causes air to warm through compression. In addition, decreases (increases) in cloud cover directly lead to increases (decreases) of downward shortwave radiation at the surface (Figs. 7 and 8). The pattern correlation between the anomalies of TCC and R d s over NA is strong of − 0.97 (Table 3). R d s anomaly caused by TCC changes is a key contributor to surface radiation heating.
On the other hand, increases (decreases) of total column vertically-integrated water vapour over NA are generally consistent with vertically-averaged tropospheric warming (cooling), indicated by the anomalies of the mid-tropospheric Φ 500 and thickness (Ф 250 -Ф 1000 ) (Fig. 8). The pattern correlation between the anomalies of TCWV and Φ 500 is 0.62 over NA (Table 3). Warming (cooling) regions of the troposphere have more (less) water vapor holding capability, which acts as a greenhouse gas to absorb (release) longwave radiation. This is a process known as water vapor feedback (Ingram 2010;Held and Shell 2012), which leads to downward longwave radiation anomalies at the surface. The pattern correlation between the anomalies of TCWV and R d l is 0.89 (Table 3). R d l anomaly caused by TCWV changes is the second major contributor to surface radiation heating.

Summary and discussion
This study analyzes the North American summertime heatwave variability from a climatological perspective, using the ERA5 reanalysis over 1959-2021 and observational data over 1959-2010. To aid in unraveling the generation and maintenance of the heatwave variability, we examine the heatwave associated anomalies of large-scale atmospheric circulation and surface energy balance, and the relevant dynamical and thermodynamic processes.
The principal mode of warm extreme anomalies, extracted from an EOF analysis of JJAS monthly TX90 anomalies over NA, is employed to characterize the dominant heatwave variability. EOF1 features a dipole pattern with a large area of anomaly over northwestern NA and an anomaly of opposite sign over the southern U.S. EOF1 over the northwestern part of NA resembles the record-breaking heatwave over western NA in June-July 2021. The heatwave variability, represented by the corresponding PC1 index, is associated with a large-scale atmospheric wave train extending from the North Pacific to North America, together with a northeastward injection of moisture from the subtropical western Pacific towards western NA. The NP-NA atmospheric teleconnection and moisture anomalies are key factors in supporting the NA heatwave variability.
The generation of the NP-NA wave train is related to internal dynamical processes such as transient eddy variations and jet instabilities. In particular, synoptic-scale eddies can systematically reinforce and support the wave train. The anomalous convection over the tropical-subtropical western Pacific, which is related to a northward propagation of the Fig. 8 Anomalies of total cloud cover TCC (interval 1.0%), precipitation (interval 0.2 mmday −1 ), relative humidity RH 1000 (interval 1.0%), vertical velocity ω 500 (interval 5.0 × 10 -5 hPas −1 ), total column vertically-integrated water vapour TCWV (interval 0.2 Kgm −2 ), geo-potential Φ 500 (interval 50 m 2 s −2 ), and thickness Φ 250 -Φ 1000 (interval 100 m 2 s −2 ) over North America regressed on the PC1 index. Black dots indicate anomalies that are significantly correlated with the PC1 index at the 5% level Southeast Asian summer monsoon system, may also act as a forcing to trigger the NP-NA teleconnection. Surface heat fluxes associated with the heatwave variability are dominated by downward shortwave and upward longwave radiation anomalies. Surface radiation heating directly contributes to surface temperature anomalies and is dominated by anomalous downward shortwave and longwave radiations. In association with a positive phase of the PC1 index, the NP-NA wave train brings positive (negative) geopotential anomalies over the northern (southern) part of NA, leading to anomalous descent, reduced total cloud cover and below-normal precipitation over northern NA, and opposite anomalies over southern NA. Over northwestern NA, the anomalous subsidence causes air to warm through compression. Reduced cloud cover results in increased downward shortwave radiation. In addition, the increase in vertically-integrated water vapour over the total column is consistent with vertically averaged tropospheric warming over northwestern NA, where the troposphere has more water vapor holding capability and acts as a greenhouse gas to absorb longwave radiation, leading to increased downward longwave radiation.
The dominant mode of TX90 anomalies over NA, as well as its associated atmospheric wave train and moisture anomalies from the North Pacific to North America, are also clearly evident in the NCEP/NCAR reanalysis over the same period (not shown). This increases confidence in the results reported here, since ERA5 and NCEP/NCAR reanalyses are generated from two different climate centers. Wu et al. (2012) identified the interdecadal and interannual modes of the NA heatwave frequency variability. The interannual mode shares some features with the dipole pattern identified here, especially for the anomalies over northwestern NA and the southern U.S. In addition, Teng et al. (2013) found that the U.S. heatwave anomaly is associated with a zonal wavenumber-5 pattern. The anomalous atmospheric planetary wave resembles the NP-NA teleconnection, while its action center over NA is located in the central U.S. The anomalous circulation center coincides with the dominant U.S. heatwave anomaly, which is also consistent with the dynamical structure over northwestern NA reported here. As a matter of fact, the U.S. heatwave anomaly and its associated circulation pattern identified by Teng et al. (2013) bear close resemblance to the second mode (EOF2) of the TX90 anomalies over NA and the corresponding circulation anomalies, respectively (not shown). The EOF2 pattern accounts for 10.5% of the total variance over NA for the period from 1959 to 2021. However, it is worth noting that there is no universal definition of a heatwave (e.g., Peterson et al. 2008;Fischer and Schär 2010;Wu et al. 2012). Different definitions highlight various aspects of heatwaves, such as the heatwave amplitude, frequency and duration. Hence, inconsistency remains in results of the NA heatwave variability, which may stem from various definitions of heatwaves, as well as different analysis methods and datasets used.
We have also analyzed the relationship between the PC1 index and the Southeast Asian Summer Monsoon index. The SEASM index used is defined as an area-averaged JJAS dynamical normalized seasonality at 850 hPa over the Southeast Asia (2.5-20° N, 70-110° E) (Li and Zeng 2002), downloaded from Prof. Jianping Li's homepage (http:// lijia nping. cn/ dct/ page/ 65544). The correlation between the two indices is positive but low (0.11) over 1959-2020. The positive correlation may be attributed to some similarities of the precipitation anomalies over the tropical western Pacific in association with the PC1 and SEASM indices (cf. Fig. 3 here with Fig. 2 of Li and Zeng 2002). Nevertheless, the PC1 variability cannot be fully captured by the SEASM index employed. Notably, inconsistency exists between the monsoon circulation and convection indices, especially for the SEASM (e.g., Wang and Fan 1999). Additionally, the NP-NA teleconnection can be influenced by cooling anomaly in the Philippine Sea and heating anomaly in the subtropical SEASM region Qian et al. 2022), which indicates differences of the relevant regional convection anomalies in the tropical-subtropical western Pacific. Hence, it remains to be explored whether there exists a robust relationship between the NA heatwave variability and East Asian summer monsoon. In addition, this study only explores the internal climate variability of NA heatwaves. It is noteworthy that external forcings, especially the anthropogenic climate change, also play important roles in shaping climate extremes including heatwaves (e.g., IPCC 2021). The relative role of external forcing and internal variability on the NA heatwave and its associated physical processes merit further investigation. 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/.