Rapid attribution analysis of the extraordinary heatwave on the Pacific Coast of the US and Canada June 2021

. Towards the end of June 2021, temperature records were broken by several degrees Celsius in several cities in the Pacific northwest areas of the U.S. and Canada, leading to spikes in sudden deaths, and sharp increases in emergency calls and hospital visits for heat-related illnesses. Here we present a multi-model, multi-method attribution analysis to investigate to what extent human-induced climate change has influenced the probability and intensity of extreme heatwaves in this region. Based on observations, modelling and a classical statistical approach, the occurrence of a heatwave defined as the maximum 5 daily temperatures (TXx) observed in the area 45 ◦ N–52 ◦ N, 119 ◦ W–123 ◦ W, was found to be virtually impossible without


Introduction
During the last days of June 2021, Pacific northwest areas of the U.S. and Canada experienced temperatures never previously observed, with temperature records broken in multiple cities by several degrees Celsius. Temperatures far above 40 • C (104 20 • F) occurred on Sunday 27 to Tuesday 29 June (Figs 1a,b for Monday) in the Pacific northwest areas of the U.S. and western Provinces of Canada, with the maximum warmth moving from the western to the eastern part of the domain from Monday to Tuesday. The anomalies relative to the daily maximum temperature climatology for the time of year reached 16 • C to 20 • C (Figs 1c,d). It is noteworthy that these record temperatures occurred one whole month before the climatologically warmest part of the year (end of July, early August), making them particularly exceptional. Even compared to the average annual maximum 25 temperatures in other years independent of the month considered, the recent event exceeds those temperatures by about 5 • C ( Figure 2). Records were shattered in a very large area, including the village of Lytton where a new all-time Canadian temperature record of 49.6 • C was set on June 29, and where wildfires spread on the following day.
Given that the observed temperatures were so far outside historical experiences and occurred in a region with only about 50% household air conditioning penetration, large impacts on health were expected. In British Columbia, Canada, there were 30 an estimated 619 heat-related excess deaths putting the event among the deadliest weather related events in Canada 1 . There were an estimated 196 extra deaths in Washington and 100 in Oregon. 2 3 Below we investigate the role of human-induced climate change in contributing to the likelihood and intensity of this extreme heatwave, following an established approach to multi-model multi-method extreme event attribution (Philip et al., 2020;. We focus the analysis on the maximum temperatures in the region where most people were affected 35 by the heat (45 • N-52 • N, 119 • W-123 • W) including the cities of Seattle, Portland, and Vancouver. While the extreme 1 https://www2.gov.bc.ca/assets/gov/birth-adoption-death-marriage-and-divorce/deaths/coroners-service/death-review-panel/extreme_heat_death_ review_panel_report.pdf 2 https://www.opb.org/article/2021/08/06/oregon-june-heat-wave-deaths-names-revealed-medical-examiner/ 3 https://doh.wa.gov/emergencies/be-prepared-be-safe/severe-weather-and-natural-disasters/hot-weather-safety/heat-wave-2021 heat was an important driver of the observed impacts, it is important to note that these impacts strongly depend on exposure, vulnerability, and other climatological variables beyond temperature. In addition to the attribution of the extreme temperatures, we qualitatively assess whether meteorological drivers and antecedent conditions played an important role in the observed extreme temperatures in Section 6. The black box indicates the study region. Source: GHCN-D

Event definition
For the attribution analysis, a definition of the event is required. As there are several ways to define an event, this step requires some expert judgement. Within this study, we will analyse the daily maximum temperatures, which characterised the event and dominated headlines in the large number of media reports describing the heatwave and its associated impacts. We therefore define the event based on the annual maximum of daily maximum temperature, TXx. Here, we first average over the region 45 and then take the annual maximum. Other options for variables that could be selected for analysis include e.g. 3-day averaged temperatures, as there is some evidence that longer time scales better describe the health impacts (e.g., D'Ippoliti et al., 2010) or high minimum temperatures which also have strong impacts on human health. However, TXx is a standard heat impact index and thus the results can easily be compared to other studies. We intentionally focus on one event definition to keep this analysis succinct and its results easy to communicate, choosing TXx, which not only characterises the extreme character of the 50 event but is also readily available in climate models allowing us to use a large range of different models. Recognising that the WMO standard definition of a heatwave is a multi-day measure associated with persistent heat, we additionally considered the annual maxima of 5-day and 10-day averages of daily maximum temperatures, TX5x and TX10x, for observations only. Trends in these time series, as seen by the intensity change, turn out to be very similar to the results for TXx presented in Sect. 3.
For the spatial scale of the event we consider the area 45 • N-52 • N, 119 • W-123 • W. This covers the more populated region 55 around Portland, Seattle and Vancouver that were impacted heavily by the heat (with a total population of over 9.4 million in their combined metropolitan areas), but excludes the rainforest to the west and arid areas to the east. Note that this spatial event definition is based on the expected and reported human impacts rather than on the meteorological extremity. Besides the analysis for this region, we also analysed long and homogeneous observational records of three stations in Portland, Seattle and Vancouver.  Figure 3 shows the observed trends in TXx in the GHCN-D dataset over 1900-2019. The stations were selected on the basis of a long time series of at least 50 years of data, and were required to be at least 2 • apart. The trend is defined as the regression on the global mean temperature, so the numbers represent how much slower or faster the temperature has changed compared to the global mean temperature. Individual stations with different trends than nearby stations usually have inhomogeneities in 65 the observational method or local environment. There are large positive trends in heatwaves in Europe. These are not yet fully understood or adequately represented in climate models . The negative trends in eastern North America and parts of California are well-understood to be the result of land use changes, irrigation and changes in agricultural practice (Donat et al., 2016(Donat et al., , 2017Thiery et al., 2017;Cowan et al., 2020 western-north-american-extreme-heat-virtually-impossible-without-human-caused-climate-change/ used ERA5 extended by the ECMWF operational analysis and the ECMWF forecast. The differences in the heatwave amplitude between the rapid analysis and the updated analysis presented here are minor and do not affect the rounded estimate of the return period used in this study. Therefore, the analysis results do not require updating.

Previous trends in heatwaves
Temperature observations were used to assess probability ratios and return periods associated with the event for three major 80 cities in the study area: Portland, Seattle, and Vancouver. Observing sites were chosen based on (i) the availability of long homogenised historical records, (ii) their ability to represent the severity of the event by avoiding sites within the proximity of large water bodies, and (iii) their representativeness for populous areas of each city to better illuminate impact on inhabitants.
For Portland, the Portland International Airport National Weather Service station was used, which has continuous observations over 1938-2021. The airport is located close to the city centre, adjacent to the Columbia River. The river's influence is 85 thought to be small. For Seattle, Seattle-Tacoma International Airport was chosen, which has almost continuous observations between 1948 and 2021, making it one of the stations with the longest records in the Seattle area. This location is further inland and thus avoids the influence of Lake Washington that affects downtown Seattle. Two long records exist adjacent to downtown Vancouver, but they are both very exposed to the Georgia Strait that influenced observations due to local onshore flow during the peak of the event. Thus the time series from a site further inland at New Westminster was selected, which starts in 1875 but 90 contains data gaps in 1882-1893, 1928, and 1980-1993. The data for Portland International Airport and Seattle-Tacoma International Airport were gathered from the Global Historical Climatology Network Daily (GHCND; Menne et al., 2012) while data for New Westminster were gathered from the Adjusted Homogenized Canadian Climate Dataset (AHCCD) for daily temperature (Vincent et al., 2020). This station's record is a composite of data from three locations in two nearby cities as location changes took place in 1966and 1980. From 1874 to 1966, the station operated at an elevation of 118 m above sea level near the centre of New Westminster. In 1966, the station was moved about 2 km east and to an elevation of 18 m. The portion of the homogenised record from 1980 onward is from Pitt Meadows, BC, located about 14 km east of the previous location at an elevation of 5 m. Using a composite station is non-ideal given the potential influence of local micro-climatic effects and particularly the increasing distance from the Strait of Georgia which exhibits a cooling effect to sites in its vicinity. Use of data from this composite site may increase the uncertainty 100 of our analysis, but given the magnitude of the signal and the consistency of results among the datasets presented here (and analysis of other temperature datasets from BC, not shown here), we do not expect this to substantially affect our results. The AHCCD dataset is updated annually and ends in 2020. Data for 2021 were appended from unhomogenized recent records from Environment and Climate Change Canada. Overlapping data for 2020 were compared between the two sources and found to be identical except for several duplicate/missing observations. Such duplicate or missing data would not cause error in the present 105 analysis because the records are complete for June 2021.
As a measure of anthropogenic climate change we use the global mean surface temperature (GMST), where GMST is taken from the National Aeronautics and Space Administration (NASA) Goddard Institute for Space Science (GISS) surface temperature analysis (Hansen et al., 2010;Lenssen et al., 2019, GISTEMP,). We apply a 4-yr running mean low-pass filter to suppress the influence of ENSO and winter variability at high northern latitudes as these are unforced variations.

Model and experiment descriptions
In this study a variety of climate model simulations are analysed and will be described in the following paragraphs.
Model simulations from the 6th Coupled Model Intercomparison Project (CMIP6; Eyring et al., 2016)  In addition an ensemble of extended historical simulations from the IPSL-CM6A-LR model is used (see Boucher et al., 2020, for a description of the model), which follows the CMIP6 protocol (Eyring et al., 2016). It is composed of 32 members, and 120 the simulations cover the historical period (1850-2014). It has been extended until 2029 using all forcings from the SSP2-4.5 scenario, except for the ozone concentration which has been kept constant at its 2014 climatology, as it was not available at the time the extensions were generated. This ensemble is used to explore the influence of internal variability.
Furthermore we use the GFDL-CM2.5/FLOR model (Vecchi et al., 2014), which is a fully coupled climate model developed at the Geophysical Fluid Dynamics Laboratory (GFDL). While the ocean and ice components have a horizontal resolution of 125 only 1 • , the resolution of atmosphere and land components is about 50 km and might therefore provide a better simulation than coarser model simulations of certain extreme weather events (Baldwin et al., 2019). The data used in this study cover the period from 1860 to 2100, and combine historical and RCP4.5 experiments driven by transient radiative forcings from CMIP5 (Taylor et al., 2011).
We also examine five ensemble members of the AMIP experiment (1871-2019) from the GFDL-AM2.5C360 (Yang et al., simulations resulting in a sample size of 3823 years.

Statistical methods
A more detailed description of the statistical methods is given in (Philip et al., 2020;. Here we give a description of the most important aspects.

145
As discussed in Section 1.1, we analyse the annual maximum of daily maximum temperatures (TXx) averaged over 45 • N- Initially, we analyse reanalysis data and station data from sites with long data records. Next, we analyse To statistically model the selected event, we use a GEV distribution (see e.g., Coles, 2001) that shifts with GMST, i.e., the location parameter has a term proportional to GMST and the scale and shape parameters are assumed constant. Uncertainties corresponding to the statistical-model uncertainty, are obtained using a non-parametric bootstrap procedure. With this GEV 160 distribution, first the PR and intensity change are calculated from observations, as well as the return period in the current climate. Next, the return period is used as a threshold to specify the event magnitude for the models. For this return period, the PRs and intensity changes between 2021 and the counterfactual climate are calculated from different models. This is, however, only done for models that pass our validation tests on the seasonal cycle, the spatial pattern of the climatology, and the scale and shape parameters of the GEV distribution, see Section 4. Finally, both observational and model results are synthesised into 165 a consistent attribution statement, see Section 5. For models (except IPSL-CM6A-LR and CAM5.1), we additionally analyse the PR between the current climate and a future climate at +2 • C above the 1850-1900 reference, which is equivalent to +0.8 • C above the current climate of 2021. For this analysis of future change we use model data up to about 2050 or when the model GMST reaches +0.8 • C compared to now.
The CMIP6 data are analysed using the same statistical models as the main method. However, for practical reasons, the   is 5.7 • C above the previous record of 34.0 • C. This extremely large increase leads to difficulties in the statistical analysis described below. There are two possible sources of this extreme jump in peak temperatures. The first is that an event with very low probability occurred -the statistical equivalent of "really bad luck". An event with very low probability could also have occurred in the pre-industrial climate but its amplitude would have been increased by climate change in the current climate which already includes about 1.2 • C of global warming. The second option is that strong nonlinear interactions and feedbacks 180 took place in this event, amplifying the intensity of the extreme, This would mean that climate change is acting to exacerbate extreme heatwaves more rapidly than implied by a GEV with a location parameter which scales linearly with GMST. In this case, the event would not belong to the "same population" as the other ones and we would not expect the method applied here to be successful. This second possibility requires further investigation. While we keep this possibility in mind, we assume the first option within this study. In Figure 5a we show the seasonal cycle of the daily maximum temperature averaged over the index region and in Figure 5b we show the spatial pattern of TXx averaged over many years at each grid point individually. These two metrics are used in the model validation procedure, see Section 4. Figure 6 shows the analysis of the gridded (ERA5) data using the ordinary extreme value analysis applied for attribution 190 studies within WWA. This approach excludes data of the extreme event of interest from the statistical fit to avoid selection bias, especially when the study area has been defined based on the extent of the extreme event. Figure 6a shows the observed TXx as a function of smoothed GMST and shows that the value observed in 2021 is far outside the range of any values observed to date. The distribution of TXx including data up to 2020 is described very well by a GEV distribution that has linearly warmed at a rate about twice as fast as the GMST, see Figure 6b. The warming rate is consistent with expectations, 195 as summer temperatures over continents increase faster than the global mean. The GEV fit has a negative shape parameter ξ, An alternative to the standard approach for which no information of the event under study is used (to avoid a selection bias),

200
is to use the information that it actually happened, yet without including the value observed in 2021 in the fit. Specifically, we again assume that the data up to 2020 can be described by a GEV distribution with constant scale and shape parameters, but we reject all GEV models in which the upper bound is below the value observed in 2021. In other words, we enforce that the fit parameters are within a subset of parameters that are compatible with the 2021 event. The result is shown in Figure 7.
While the distribution now includes the 2021 event, the fit to the data up to and including the year 2020 is noticeably worse than when not taking 2021 into account. The return period for the 2021 event under these assumptions still has a lower bound of 10,000 years in the current climate. The fit differs from the previous one mainly in the shape parameter, which is now much less negative (about −0.2 instead of −0.4). This shifts the upper bound to higher values. The fit also gives a somewhat higher trend parameter.
The third possibility is to fit the GEV distribution over all available data, including observations from 2021. This approach 210 implicitly assumes that the 2021 event is drawn from the same population. We typically do not make this assumption in cases where we intentionally select a specific region in order to maximise the extremity (i.e. return period) of the event to avoid a selection bias. Note however that we may have overestimated the return period by excluding the event, a question left for future investigations. In intentionally choosing the region with the largest (rarest) return period for the 2021 event, the extreme value is drawn from a different distribution and cannot be included. However, this is only partly the case here. We did choose the 215 general region because the temperatures were exceptional there, however, we based the choice of subregion for the analysis on population density and type of terrain -parameters that are independent of the heatwave. The benefit of this approach is that it uses all available information. With this approach we still assume this was an event happening by chance, that is, the behaviour is in line with that of a chaotic deterministic system and by chance we observe a low-probability event in this short time series. Given the extremity of the event and the relatively low number of data, a robust GEV fit is hard to obtain, and the appropri-220 ateness of the method is difficult to assess. However, the application of this classical method in this case is interesting provided we keep in mind the assumptions we make. While we acknowledge that none of the three possibilities to fit a GEV distribution is fully satisfying, we decided on using this third approach to estimate the return period leading to an estimate of 1,000 years (95% CI >100 yr). Follow-up research will be necessary to investigate the potential reasons for this exceptional event and the consequences for assumptions for these fits (see also the discussion in Sections 6 and 8). Also, further research is needed into

Analysis of temperature in Portland, Seattle and Vancouver
To represent Portland we chose the International Airport station, which is located on the northern edge of the city and has been collecting data since April 1938; the data are in the GHCN-D v2 database. We fit a GEV distribution to this data, including 2021 ( Figure 9, lower panels). It gives a return period of 700 yr for the 2021 event with a lower bound of 70 yr. For the PR we can only give a lower bound (6), since the best estimate is infinite. This corresponds to an increase in TXx of 3.4 • C with a large uncertainty of 0.3 to 5.3 • C. The large uncertainties are due to the somewhat shorter time series and large variability at this station.

240
In Seattle, the only station with a sufficiently long time series that includes 2021 is Seattle-Tacoma International Airport. It is located ∼15 km south of the city but has similar terrain, without the proximity to water of the city itself. The previous record was 39.4 • C in 2009, and in 2021 it reached 42.2 • C. This is still a large increase of 2.8 • C over the previous record. The event was also not quite as improbable, with a return period of 300 yr (lower bound 40 yr) in the current climate ( Figure 10). The PR is again infinite with a lower bound of 7, and the increase in temperature from a late nineteenth century climate is 3.8 • C (0.7 245 • C to 5.7 • C).

Model evaluation and analysis
In this section we show the results of the model validation. The validation criteria assess the similarity between the modelled and observed seasonal cycle, the spatial pattern of the climatology, and the scale and shape parameters of the GEV distribution.
The assessment results in a label "good", "reasonable" or "bad", according to the criteria defined in Ciavarella et al. (2021). In this study, we use models that are labelled "good" or "reasonable". However, if five or more models are classified as "good" 255 within a particular model set such as the CMIP6 models, then we do not include all of the "reasonable" models but only those that pass the specific test on fit parameters as "good". Table 1 shows the model validation results. The full table including the models that did not pass the validation tests is given in Table 3. In total 21 models and a combined 224 ensemble members passed the validation tests.
Next, we show probability ratios and change in intensity ∆T for models that pass the validation tests and we also include 260 the threshold values for a 1-in-1000 year event (Table 2). Results are given both for changes in the current climate (1.2 • C) (c) Figure 11. as Figure 9 but for the station data at New Westminster. Data source: AHCCD.
compared to the past (pre-industrial conditions) and, when available, for a climate at +2 • C of global warming above preindustrial climate compared with the current climate. The results are visualised in Section 5.

Hazard synthesis
In Sections 3 and 4 we calculated the probability ratio as well as the change in magnitude of the event in the observations and the  of an unweighted average of observations and the model average. This will partly correct for a common model bias, if the observations are reliable, and is indicated by the white box accompanying the magenta bar in the synthesis figures. Figure 12 shows the synthesis results for the current vs. past climate; the results for the future vs. current climate are presented in Figure 13. Where the results for the probability ratio do not give a finite number, we replace them by 10000 to 280 allow all models to be included in the synthesis analysis. This means that the reported synthesised probability ratio gives a more conservative, lower value. For the intensity change we report the weighted synthesis value. For the probability ratio we can only give a lower estimate of the range. Generally, we do not see any consistent departures in the model results that can be traced back to experiment differences, except that model ensembles which consist of many members tend to have smaller uncertainties as expected.

285
Results for current vs. past climate, i.e., for 1. which is the scenario closest to current emission levels, be reached as early as the 2040s (Lee et al., 2021).

The broader context of the heatwave
In the previous section we summarised and synthesised trends in TXx that were detected in observations and attributed to climate change using model data. In this section we provide some context to the analysed heatwave event by evaluating 295 the assumption that this heatwave occurred in this location by chance and by discussing factors that possibly influenced the extremity of the event, being the specific meteorological conditions and dynamics, preceding dryness -which can amplify temperature during heatwaves and reduce evaporation, and the ENSO and PDO modes of of natural variability that are relevant for this region.

Probability of a chance event 300
In Section 3 we offered two explanations for having such a record breaking heatwave as that studied here; it could either occur by chance (a low-probability event) or, for instance, nonlinear effects that have not been observed at this location before could have made such an extreme heatwave possible. Here we provide some context to the first option by providing an estimate of how many heatwaves that are characterised by a return period of 1000 years at their given location we can expect on the entire globe. In this study our analysis focuses on the area 45 • N-52 • N, 119 • W-123 • W which was strongly affected by the heatwave.

305
However, the entire area affected by heat is larger -about 1500km x 1500km which is about 1.5% of the land area of the world.
Assuming that the event occurred just by chance also over the entire area affected by the heatwave and that the return period is similar to the 1000 years obtained for the study area, we can roughly estimate the global return period, i.e. the worldwide Heatwaves today compared to the pre-industrial climate two decades. Further research on this and other exceptional heatwaves is needed to determine whether this estimate is indeed realistic, i.e., whether or not we should reject the assumption that this heatwave occurred by chance at this location.

Meteorological analysis and dynamics
The evolution of this event can be explained by a confluence of meso-and synoptic-scale dynamical features, potentially including antecedent low soil moisture conditions and anomalously high column specific humidity that are a hallmarks of 320 extreme heat in western North America (Stewart et al., 2017;Bumbaco et al., 2013). At the synoptic scale, an omega-block developed over the study area beginning at roughly 00UTC on June 25th centred at ∼125 • W, 52 • N, which then slowly progressed eastward over subsequent days. This ridge featured a maximal 500 hpa geopotential height of ∼5980 m, which is unprecedented for this area of western North America for the period from 1948 through to June 2021 at least (Figure 14).
Despite being a record, this extreme high pressure system -a feature sometimes called a "Heat dome" -is not that 325 anomalous given the long-term trend in 500 hPa driven by thermal expansion (Christidis and Stott, 2015). Also, comparing The circulation pattern itself also appears typical for hot summertime temperatures: using analogues of 500 hPa and a pattern correlation metric to compare fields, we find that about 1% of June and July circulation patterns, defined as the 500 330 hPa geopotential height pattern within [160 • W-110 • W; 35 • N-65 • N] in previous years have an anomaly correlation larger than 0.8 with the 28 June pattern. This degree of correlation is typical among days with this type of blocking pattern during the months of June and July. Roughly one third of June and July geopotential height fields have 1% or fewer analogues with an anomaly correlation larger than 0.8. We also find that this fraction does not change when restricting the analogues search within 3 distinct time periods between 1948 and 2020. We conclude that the 28 June circulation is probably not exceptional, 335 while temperatures associated with it were.
At the mesoscale, high solar irradiance during the longest days of the year and strong subsidence increased near-surface air temperatures during the event. As is typical for summer heatwaves in the region (Brewer et al., 2012(Brewer et al., , 2013, a mesoscale thermal trough developed over western Oregon by 00UTC on the 28th June. This feature migrated northward reaching the northern tip of Washington State by 00UTC on the 29th. Further offshore, a small cut-off low travelled southwest to northeast 340 around the synoptic-scale trough that made up the west arm of the omega block. The pressure gradients associated with the thermal trough and the cut-off low promoted moderate E-SE flow in the northern and eastern sectors of the feature and S-SW flow to the south. Near-surface winds with easterly components crossed the Cascade Range of Washington and Oregon and the southern Coast Mountains of British Columbia where they were lighter but sufficient to displace cooler marine air. The difference in elevation on the west and east sides of the mountain ranges contributed to more adiabatic heating than cooling, 345 which helped drive the warmest temperatures observed in the event along the foot of the west slope of these mountains, near or slightly above sea level. These dynamics are illustrated in Figure 15. By 12UTC on the 29th of June 2021, the southwestern portion of the study area was under the influence of southerly to southwesterly near-surface flows that advected marine air and forced marked cooling. Unfortunately, winds associated with this transition intensified a wildfire that quickly consumed the town of Lytton, BC where Canada's nationwide all time high temperature was set just a day before.

350
There is no scientific consensus whether blocking events are made more severe or persistent because of Arctic amplification or other mechanisms (i.e. Tang et al., 2014;Barnes and Screen, 2015;Vavrus, 2018). We contend that Arctic sea ice was unlikely to have played a large role in this event largely due to the timing. In early summer, Arctic sea ice remains extensive but continues to melt, thus keeping near surface temperatures near 0 • C. This causes summer trends in near-surface temperatures over the Arctic ocean to be smaller than for the midlatitudes. During the months prior to the event, the sea ice extent was below 355 the 1981-2010 mean, but was similar to values observed from 2011 to 2020 (Fetterer et al., 2017 (updated daily). Instead, Arctic Amplification in summer is characterised by strong warming over high-latitude land areas (as can clearly be seen in Figure   16) and this warming signal reaches into the upper-troposphere. This enhanced warming is likely related to strong downward trends in early summer snow cover. There is evidence, from observations (Coumou et al., 2015;Chang et al., 2016), climate models (Harvey et al., 2020;Lehmann et al., 2014) and paleo-proxies (Routson et al., 2019), that this enhanced warming over 360 high latitudes leads to a weakening of the jet and storm tracks in summer. This weakening could favour more persistent weather conditions (Pfleiderer et al., 2019;Kornhuber and Tamarin-Brodsky, 2021). Regional-scale interactions between loss of snow cover and low soil moisture associated with earlier snowmelt and rapid springtime soil moisture drying, may have had an enhanced warming impact into early summer in the Arctic. At mid-atmospheric levels there is some amplification remaining due to the winter season ( Figure 16), but at the jet level (∼250 hPa) the usual increase of the thermal gradient due to tropical 365 upper tropospheric warming is advected North by the Hadley circulation (Haarsma et al., 2013). The final effect on the jet stream is therefore a competition between factors enhancing and decreasing the temperature gradient.

Drought
An additional feature of the event is the very dry antecedent conditions that may have contributed to the observed extreme temperatures through reduced latent cooling from inhibited evapotranspiration. Low soil moisture conditions can lead to a strong 370 amplification of temperature during heatwaves, including nonlinear effects (Seneviratne et al., 2010;Mueller and Seneviratne, 2012;Hauser et al., 2016;Wehrli et al., 2019). In addition, low spring snow level conditions can also further amplify this feed- back (Hall et al., 2008). In this section we briefly explore whether precipitation anomalies and evapotranspiration measures could have played a role in the extreme heat via feedbacks related to soil moisture conditions.  The available moisture is also influenced by evapotranspiration, which depends strongly on temperature, radiation and available atmospheric moisture. Evaporation in the study area was below normal in the ERA5 reanalysis from March and became more negatively anomalous until May (not shown). During the event in late June, there was progressive soil desiccation, creating ideal conditions for strong negative evaporation anomalies. On the other hand, surface net radiation was high, especially before peak temperatures were reached (June 26-28). As a consequence, evaporation deficits during the heat build-up were 385 rather slight compared to most days in May and the first half of June. Together with the extreme near-surface temperatures that suppressed surface heating due to an already hot surface, this resulted in only moderately positive surface sensible heat flux anomalies.
Satellite-based measurements of surface soil moisture based on microwave remote sensing from the European Space Agency (ESA) Climate Change Initiative (CCI) provided by the Copernicus service suggest that surface soil moisture was below 390 normal in the region since the beginning of April and that the anomalous conditions persisted until June (https://dataviewer. geo.tuwien.ac.at/?state=88bf0c), in agreement with the decreased precipitation and somewhat below normal evapotranspiration in the ERA5 reanalysis.

Influence of modes of natural variability
The El Niño Southern Oscillation is the dominant source of interannual variability in the region through the Pacific North Altogether, external modes of variability appear to have played little to no role in the formation of the event.
7 Vulnerability and exposure The Pacific Northwest region is not accustomed to very hot temperatures such as those experienced during the June 2021 405 heatwave. Heatwaves are one of the deadliest natural hazards, resulting in high excess mortality through direct impacts of heat (e.g. heat stroke) and by exacerbating pre-existing medical conditions linked to respiratory and cardiovascular issues (Haines et al., 2006;Ebi et al., 2021). In addition to at least 815 excess heat-related deaths, there was a significant increase in emergency department visits. 4 5 On June 28, 2021 alone, there were 1,038 heat-related emergency department visits in the U.S. The June 2021 heatwave also affected critical infrastructure such as roads and rail and caused power outages, agricultural impacts, and forced many businesses and schools to close. 7 8 Rapid snowmelt in BC caused water levels to rise, leading to 415 evacuation orders north of Vancouver. 9 Furthermore, in some places, wildfires, the risk of which has increased due to climate change in this region (Kirchmeier-Young et al., 2019), started and quickly spread, requiring entire towns to evacuate. 10 The co-occurence of such events may result in compound risks, for example when households are advised to shut windows to keep outdoor wildfire smoke from getting inside while simultaneously being threatened by high indoor temperatures when lacking air conditioning.

420
Timely warnings were issued throughout the region by the US National Weather Service, Environment and Climate Change Canada, and local governments. British Columbia has a "Municipal Heat Response Planning" summary review that gathers information on heat response plans throughout the province, including responses such as increasing access to cooling facilities and distribution of drinking water. In long-term strategies, changes to the built environment are emphasised (Lubik et al., 2017).
Not all municipalities throughout the Pacific Northwest and BC have formalised heat response plans, and others have limited 425 planning, thought to be due to low heat risk perceptions throughout the area, as well as a lack of local data for risk assessments (Lubik et al., 2017).
The extremely high temperatures that occurred in this heat episode meant that everyone was vulnerable to its effects if exposed for a long enough period of time. This record-breaking extreme event has been analysed under the assumption that it was simply a low probability random event. A rudimentary calculation looking into the probability of a random event of similar extent and severity to occur anywhere over the earth's land area, gave an estimated chance of order 1-in-15 years, which at first impression makes a random event 455 seem plausible, but this should be more thoroughly investigated.
The alternative is that nonlinear interactions and feedbacks occurred, which amplified the intensity of this extreme, placing it in a different population of heatwave events with different (and possibly unknown) statistics. We briefly considered dynamical and hydrological (drought) mechanisms and modes of natural variability that could have had an amplifying role. The conditions in the preceding months were dry but not extremely anomalous. The circulation itself was highly anomalous but not exceptional 460 enough to explain the record breaking heat alone, however, local topography and preceding dryness may have amplified the associated temperatures. Also, it cannot be excluded that dynamical mechanisms (Arctic Amplification) at work influenced the persistence of blocking conditions.
Further research is planned to investigate whether these or other feedbacks were operating in this exceptional event, and whether those feedbacks are related to human-induced climate change and if they increase the frequency beyond that expected 465 for random events of such extreme temperatures. Also, further research is needed to overcome the known limitations of standard GEV analysis on annual maxima with short records and very extreme values.
Whether or not local or dynamical feedbacks are responsible for amplifying the extreme temperatures in this particular event, this study shows that the human-induced warming that has occurred since pre-industrial conditions does make extreme events like this possible in the current climate and study region, and many times more likely than in the pre-industrial era.

470
Adaptation measures therefore need to be much more ambitious and take account of the rising risk of heatwaves around the world. Although this extreme heat event is still rare in today's climate, the analysis above shows that the frequency will increase with further warming. Deaths from extreme heat can be dramatically reduced with adequate preparedness action -a number of adaptation and risk management priorities are becoming clear: it is crucial that local governments and their emergency management partners establish heat action plans to ensure well coordinated response actions during an extreme heat event -475 tailored to high-risk groups (Ebi, 2019). Heatwave early warning systems also need to be improved, which includes tailoring messages to inform and motivate vulnerable groups, as well as providing tiered warnings that take into account vulnerable groups may have lower thresholds for risk (Hess and Ebi, 2016). In other words, it is important to start to warn the most vulnerable as temperatures start to rise even though the general population is not yet acutely at risk. In cases where heat action plans and heat early warning systems are already robust, it is important that they are reviewed and updated to capture the 480 implications of rising risks -every five years or more often (Hess and Ebi, 2016). Further, heatwave early warning systems should undergo stress tests to evaluate their robustness to temperature extremes beyond recent experience and to identify modifications to ensure continued effectiveness in a changing climate (Ebi et al., 2018).