The 2021 Western North American Heatwave and Its Subseasonal Predictions

Record‐breaking above‐normal temperatures were observed across western North America in June–July 2021. In this study, our ability to predict this heatwave 2–3 weeks in advance is assessed based on 10 Subseasonal‐to‐Seasonal prediction models. It is found that the above normal temperature in Western Canada during June 28–July 4 was predicted by most models as early as June 10. However, for the forecasts initialized earlier than June 17, not a single ensemble member of all the models is able to capture the magnitude of the observed temperature anomaly. We identify two important processes: an upper tropospheric wave train associated with the boreal summer intraseasonal oscillation in Southeast Asia and an anomalous North Pacific atmospheric river leading to high moisture conditions. Most models are able to predict the wave train across the North Pacific. A realistic representation of moisture transport and its pattern appears crucial for the extended‐range forecast of this heatwave.


Data and Methodology
One set of the data used in this study is the ERA5 reanalysis (Hersbach et al., 2020). Variables analyzed include 2-m air temperature (T2m), 200-hPa and 500-hPa geopotential height (Z200 and Z500), 850-hPa zonal wind (U850), total column water vapor (TCWV), vertically integrated vapor transport (IVT), and mean large-scale and convective precipitation rates. The daily averaged satellite-observed outgoing longwave radiation (OLR) data from the National Oceanic and Atmospheric Administration (NOAA) polar-orbiting series of satellites (Liebmann & Smith, 1996) are also analyzed. Weekly anomalies are calculated as departures from the climatology of the same week over the 30-year period from 1991 to 2020.
For the S2S forecast, we use the real-time and reforecast data from the archive database of the S2S prediction project (Vitart & Robertson, 2018;Vitart et al., 2017). Data from 10 operational centers are analyzed. Table S1 in Supporting Information S1 lists the models and their description along with the related forecast and reforecast information. Variables analyzed are T2m, Z500, Z200, OLR and U850. Total column water (TCW) is provided in several models, which is also analyzed. TCW includes cloud water and cloud ice in some models, and thus may be slightly different from TCWV. For the heatwave, we focus on the forecast target week of June 28-4 July 2021. We verify forecasts initialized on three dates: 24th, seventeenth and tenth of June, with lead time of 5-11, 12-18 and 19-25 days. The model climatology is calculated using the available reforecast data initialized each year on dates that are close to those of the real-time forecast for each model, and thus varies between models and is usually different from the one calculated from the ERA5 reanalysis. These differences in climatology are not expected to affect our results.
To represent the boreal summer intraseasonal oscillation (BSISO) in the East Asian-western North Pacific (EAWNP) region, the EAWNP ISO index defined in Lin (2014) is used, which is calculated as principal components (PCs) of the two leading empirical orthogonal function (EOF) modes of the EOF analysis on the combined anomaly fields of 90°-150°E zonally averaged U850 and OLR from 10°S to 40°N. This EAWNP ISO index has been applied in previous studies to analyze the impact of atmospheric subseasonal variability on persistent heavy rainfall in South China, tropical cyclone genesis over the western North Pacific, and Mei-Yu Onset (e.g., Gao et al., 2016;You et al., 2018;Yao et al., 2019). Here, we look at the association between tropical convection and atmospheric circulation preceding the June 2021 heatwave, and compare with the composite response to the EAWNP ISO pattern. An agreement would indicate that the tropical convection anomaly in this particular event was associated with the EAWNP ISO and its related teleconnection was not likely caused by chance. Different BSISO indices have been designed in previous studies to represent the northward BSISO propagation (e.g., Kikuchi et al., 2012;Lee et al., 2013). The results presented here are not expected to be sensitive to the BSISO index definition.
The simple general circulation model (SGCM) as described in Hall (2000) is applied to study the atmospheric response to a specified tropical heating pattern. The SGCM is a global spectral model with a resolution of T31 and 10 vertical levels, based on the primitive equation model of Hoskins and Simmons (1975). The model solves equations of vorticity, divergence, surface pressure and temperature with a scale-selective dissipation and a level-dependent linear damping on temperature and momentum. Previous studies show that this SGCM simulates well atmospheric response to specified tropical thermal forcing (e.g., Hall, 2000;Hall & Derome, 2000;Lin & Brunet, 2018;Lin et al., 2010). In this study, we use the approach described in Hall and Derome (2000) to calculate the linear response to an imposed heating perturbation. The basic state used is the model climate of a perpetual early summer (May-June-July) SGCM integration of 3,600 days under a time-independent forcing that is calculated from the reanalysis daily data.

S2S Prediction of the Heatwave
Shown in Figures 1a and 1b are the observed T2m anomalies averaged during the weeks of June 21-27 and June 28-4 July 2021, based on the ERA-5 reanalysis. The maximum strength of the heatwave occurred in the week of June 28-July 4. Above normal T2m is seen over a large area of western North America, with a maximum weekly mean anomaly over 10°C in British Columbia and Washington State.
To see the evolution of atmospheric circulation associated with the heatwave, Figures 1c and 1d show the weekly averaged 500-hPa geopotential height (Z500) anomalies for June 21-27, and June 28-July 4, respectively. During  Figure 1c), a wave train is seen in the North Pacific originating from the tropical western Pacific, with a positive geopotential height anomaly center located over the northeast North Pacific near the west coast of Canada. This positive height anomaly center moves eastward onto the North America continent, expands and intensifies in the following week, providing a condition favorable for the high temperature ( Figure 1d).
We focus on the S2S model forecast for the week of the maximum heatwave, June 28-July 4. The ensemble mean forecasts of 2-m temperature anomaly for this target week by the 10 S2S models with the initial condition dates of 24th, seventeenth, and tenth of June 2021 are calculated, which correspond to lead time of 5-11, 12-18 and 19-25 days, respectively. Shown in Figures 2a-2c are the ensemble mean T2m anomalies predicted by the ECMWF model. The forecasts of T2m by the other nine models are presented in Figures S1-S3 in Supporting Information S1. It is seen that from the initial condition of June 24th almost all the models are able to predict the T2m anomaly distribution and amplitude of the observed heatwave in western North America. Thus, the heatwave can be accurately predicted with a lead time of 5-11 days. For the initial condition dates of June 17th and tenth, most models predict warmer than normal T2m in the heatwave area, but the amplitude is much weaker than the observation. Shown in Figures 2g-2i is the predicted T2m anomaly averaged over the western North American region (green rectangle in Figure 1b) for the target week of June 28 to July 4 by individual ensemble members as well as the ensemble mean. For the forecast initialized on June 24 (Figure 2g), all ensemble members of all the models predict above normal temperature, and the ensemble mean T2m anomaly of most models has a magnitude close to the observed value (the blue line). For the forecasts initialized on June 17 and 10 (Figures 2h Figure 1. (a) and (b) Observed weekly mean T2m anomaly for June 21-27 and June 28-4 July 2021, respectively, based on ERA5 reanalysis. (c) and (d) Observed weekly mean Z500 anomaly for June 21-27 and June 28-4 July 2021, respectively, based on ERA5 reanalysis. The anomaly is calculated as departure from the climatology of the same week over 30 years of 1991-2020. The green rectangle in (b) is the area used to calculate area mean T2m anomaly (Tav) for the heatwave. and 2i), most models are still able to produce statistically significant ensemble mean above normal T2m forecast. However, not a single ensemble member of all these models predict a T2m anomaly as strong as the observations. This indicates that the current S2S models have enormous difficulties to predict the strength of such heatwaves more than two weeks in advance.
The ensemble mean forecast of Z500 anomaly for the week of June 28-July 4 starting from the 24th, seventeenth and tenth of June 2021 are shown in Figures 2d-2f for the ECMWF model and Figures S4-S6 in Supporting Information S1 for the other models. It can be seen that that all the S2S models predict well the distribution and strength of the observed Z500 anomaly (Figure 2d and Figure S4 in Supporting Information S1). Specifically, the positive Z500 anomaly over the continental western North America is predicted by these models. For the forecast initialized on June 17, however, the predicted positive Z500 anomaly center is located over the eastern North Pacific Ocean, which is common to all the models (Figure 2e and Figure S5 in Supporting Information S1). This westward shift of the positive Z500 anomaly compared to the observation leads to discrepancies in the heatwave forecast. To further demonstrate the impact of the positive Z500 anomaly over western Canada, we compare the ensemble members of the June 17 forecast that produce area-averaged T2m anomaly (Tav) greater than 4°C for the target week with the other ensemble members. Indeed, those ensemble members producing a Tav greater than 4°C also produce positive Z500 anomalies over the continental western North America, whereas the other ensemble members generate positive Z500 anomalies over the eastern North Pacific Ocean (Figures S7a and S7b in Supporting Information S1).

Mechanisms Contributing to the Heatwave
This heat episode was extremely rare. It was found in a recent attribution study that this event was virtually impossible without human-caused climate change (World Weather Attribution, 2021). In addition to climate change, several processes were found to contribute to this heatwave, including strong rainfall in China in late June which disturbed the jet stream (Sullivan & Malik, 2021), persistent atmospheric blocking and heat dome over western North America (e.g., World Weather Attribution, 2021; Sullivan & Malik, 2021). However, how the jet stream was disturbed and what led to the formation of the blocking and heat dome are unclear. Here we propose two mechanisms for this heatwave event, i.e., teleconnection related to the East Asia -western Pacific summer monsoon variability and the effect of moisture transport. It is important for S2S models to capture these processes in order to make reliable heatwave predictions.

Teleconnection
As seen in Figure 1c, in the week preceding the strong heatwave, there was a Rossby wave train across the North Pacific originating from the subtropical western Pacific. An element of this wave train, the positive geopotential height anomaly near the west coast of North America later developed and led to a major blocking system over western North America (Figure 1d). How this wave train was initiated is thus of great interest.
Tropical diabatic heating anomalies, such as those associated with El Niño-Southern Oscillation (ENSO) and the Madden-Julian Oscillation (MJO; e.g., Madden & Julian, 1994;Zhang, 2005), have been observed to influence middle and high latitude weather and climate through teleconnections (e.g., Horel & Wallace, 1981;Lin et al., 2010). To see if a tropical forcing was involved in the June 2021 heatwave event, Figure 3a presents OLR anomalies in the tropical Indian Ocean and western North Pacific region during the week of June 21-27, 2021. Positive OLR anomalies are observed in the EAWNP region near 90°-150°E along about 10°N, which correspond to suppressed convection and reduced diabatic heating.
The above observed OLR anomaly is likely associated with BSISO, which has a strong component of northward propagation (e.g., Yasunari, 1979;Lau & Chan, 1986), in contrast to the MJO that mainly propagates eastward along the equator. Here we use the EAWNP ISO index as defined in Lin (2014) to represent the BSISO in the EAWNP region. Eight phases of the EAWNP ISO correspond to different latitudinal locations of the ISO convection ( Figure S8a in Supporting Information S1). During the second half of June 2021, the EAWNP ISO was strong (many days with amplitude greater than 2) and in Phase 7 ( Figure S8b in Supporting Information S1). The OLR anomaly during June 21-27, 2021 (Figure 3a) does look similar to that of EAWNP ISO Phase 7. Historically, during and following Phase 7 of the EAWNP ISO, western North America experiences above normal air temperature (Figures 3b and 3c), and in the upper troposphere a wave train across the North Pacific can be observed (Figures 3d and 3e) with positive geopotential height anomalies near the west coast of North America. Such features are similar to what was seen for the present heatwave (Figure 1c).
The tropical teleconnection in the S2S models is assessed by comparing Z500 anomalies of June 28-July 4 between ensemble members that successfully predicted EAWNP ISO Phase 7 during June 21-27 with those that predicted non-phase 7. As shown in Figure S9 in Supporting Information S1 for the ECMWF model initialized on the 17 June 2021, the positive Z500 anomaly in the northeast Pacific for Phase 7 members is much stronger than that for the non-phase 7 members, indicating that the reduced convection activity in the EAWNP region during that period was contributing to the amplification of the Z500 anomaly near the west coast of North America.
To further confirm that the North Pacific wave train is induced by the tropical diabatic heating anomaly of EAWNP ISO Phase 7, a numerical experiment is conducted using the SGCM (e.g., Hall & Derome, 2000;Lin et al., 2010) to find the linear response to an imposed thermal anomaly. With a 3-dimensional basic state close to the May-June-July climatology, a negative diabatic heating centered at 120°E, 10°N ( Figure S10 in Supporting Information S1) which represents the effect of EAWNP ISO Phase 7 is added to the temperature equation. The cooling perturbation has a vertical profile with maximum near 350 hPa with a vertically averaged heating rate at the center of −2.5 K/day, representing reduced deep convection, and has an elliptical form in the horizontal with a semi-major axe of 28° of longitude and a semi-minor axe of 11° of latitude. The 200-hPa geopotential height anomaly response at day 10 and day 15 indeed shows a Rossby wave train across the North Pacific, with a positive anomaly center in the eastern North Pacific near the Canadian coast (Figures 3f and 3g).

Moisture Contribution
The tropical heating anomaly of EAWNP ISO alone is not able to cause the strong heatwave, as the tropical heating induced positive tropospheric geopotential height anomalies are mainly located over the ocean (Figures 3f  and 3g). There must have been other mechanisms which made the positive geopotential height anomaly center move over the continent and intensify as observed in Figure 1d.
Here we look at the contribution of atmospheric moisture and analyze the forecast TCW output available in several S2S models. For the forecast initialized on June 24, there is higher ensemble mean forecast TCW anomaly during June 28-July 4 over western North America than the forecast initialized on June 17 ( Figure S11 in Supporting Information S1), indicating that a more accurate forecast of the heatwave is associated with a higher content of moisture. The shorter lead-time (initialized on June 24) forecast of TCW for the week of June 28-July 4 is much more successful in terms of reproducing the observed TCWC anomaly than the longer lead-time (initialized on June 17) forecast. Then, we compare the ensemble members of the June 17 forecast that produce Tav >4°C during the week of June 28-July 4 with ensemble members predicting Tav <4°C. As is evident from Figure  S12 in Supporting Information S1, the TCWV is also better represented in the warmer forecasts than colder forecasts, and ensemble members predicting higher T2m have higher atmospheric total moisture in western North America. This highlights the importance of moisture for this heatwave and its prediction. During the week of heatwave, in spite of high moisture content there was small amount of precipitation in the heatwave region (see Figure S13b in Supporting Information S1) and thus likely few clouds to reflect the solar radiation. The abundant moisture in the air could act as a greenhouse gas, absorbing longwave radiation, which further warmed the air that had already warmed up as a result of descending motion in the high pressure dominated region. It is likely that the moisture effect also helped the positive geopotential height anomaly center of the wave train move to the continent and intensify.
One may argue that the increased atmospheric moisture in the warmest ensemble members may be a consequence of the members being particularly warm, instead of a cause for the warm temperature. However, most models have the same soil moisture initial condition for all the members, excluding the possibility that different evaporation might occur for these ensemble members.
One question that remains is where the moisture in western Canada during that period came from. An analysis of atmospheric water vapor transport shows that it was caused by a strong atmospheric river (AR; e.g., Ralph et al., 2018) around June 25 (Figure S13a in Supporting Information S1), which ran over the Coast Mountains and penetrated into the continent. This AR produced precipitation only on the windward side of the mountains Most S2S models were able to predict the AR across the North Pacific several days in advance, as demonstrated in Figure 4 for the forecast TCW anomaly during 21-27 June 2021 initialized on June 17. However, the models failed to capture the observed high value of water vapor penetrated into the continent. This is likely an important reason why the models were unable to predict the amplitude of this extreme heatwave more than two weeks in advance.

Summary and Conclusions
Accurate S2S predictions of extreme weather events are of great societal and economic values. In this study, we evaluate the ability of the state-of-the-art operational S2S systems in predicting the record-breaking heatwave that occurred in June-July of 2021 over western North America. The result shows that most of the S2S models were able to predict the above-normal temperature 2-3 weeks in advance, but the magnitude of the heatwave was severely underestimated.
Two dynamical processes are found important for this heatwave. First, during late June 2021, weak monsoon precipitation in Southeast Asia associated with the summer intraseasonal oscillation induced a Rossby wave train across the North Pacific, leading to a positive geopotential height anomaly center near the west coast of North America. This positive geopotential height anomaly intensified and moved eastward onto the North American continent, providing a condition favorable for the heatwave. Most S2S models were able to predict this wave train, but failed to capture the eastward movement of the positive geopotential height center onto the North American continent. Second, a significant atmospheric river was formed during the period preceding the heatwave that transported substantial moisture into the western North American region. The S2S models demonstrated little success in predicting the penetration of moisture into the North American continent, which is likely the main reason for the failure of the forecast for the heatwave strength. To improve the prediction skill for extreme heatwaves, it is crucial to better represent the moisture transport in the S2S models.

Data Availability Statement
This work is based on S2S data. S2S is a joint initiative of the World Weather Research Programme (WWRP) and the World Climate Research Programme (WCRP). The S2S database is hosted at ECMWF as an extension of the TIGGE database, at https://apps.ecmwf.int/datasets/data/s2s.