Spatiotemporal Evolution of Heat Wave Severity and Coverage Across the United States

Heat waves have pronounced impacts on human health, ecosystems, and society. Heat waves have become more frequent and intense globally and are likely to intensify further in a warming climate. Across the United States there is a warming trend in average surface temperatures, but concordant increase in heat wave severity appears absent. Limitations in heat waves studies may be responsible for limited detection of a heat wave warming signal. We track daily spatiotemporal evolution of heat waves using geometric concepts and clustering algorithms to investigate how heat manifests on the land surface. We develop a spatial metric combining heat wave frequency, magnitude, duration, and areal extent. We find mixed trends in some individual heat wave characteristics across the United States during 1981–2018. However, exploration of the spatiotemporal evolution of combined heat wave characteristics shows considerable increases during this period and indicates a substantial increase in heat wave hazard across the United States.


Introduction
Extreme heat is detrimental to human health as it exacerbates preexisting medical conditions and causes mortality rates to increase (Ahmadalipour et al., 2019;Ahmadalipour & Moradkhani, 2018;Curriero et al., 2002;Ebi, 2008;Gasparrini et al., 2015). Heat waves also have numerous adverse ecosystem and societal impacts through reducing agriculture and vegetation productivity (Ciais et al., 2005;Mishra & Cherkauer, 2010), degrading air quality (Shaposhnikov et al., 2014;Vautard et al., 2005), inducing infrastructure failures (Garcia-Herrera et al., 2010;McEvoy et al., 2012), and increasing energy use (Miller et al., 2008). Heat waves have also been linked to enlargement of the size of wildfires causing severe vegetation and economic impacts (Zamuda et al., 2013). Heat waves are reportedly occurring more frequently across the globe, and under a warming climate they are projected to increase in frequency, magnitude, and duration (Barriopedro et al., 2011;Coumou & Rahmstorf, 2012;Diffenbaugh & Ashfaq, 2010;IPCC, 2014;Seneviratne et al., 2014). A recent study found that almost half of the world's population would likely be exposed to lethal heat wave events annually by 2100 (Mora et al., 2017). However, at the subcontinental scale of the United States there has been limited observational evidence of increased heat wave severity (Mazdiyasni & AghaKouchak, 2015;Schoof et al., 2017).
Positive warming trends in average surface temperature have been observed over the majority of the United States in recent decades, but evidence of increases in heat wave magnitude, frequency, and areal extent are much less pronounced (Grose et al., 2017;Peterson et al., 2013;Smith et al., 2013). Previous studies have been limited in their exploration of the spatiotemporal dynamics of heat waves as they have aggregated events both seasonally and regionally and, in general, employed heat wave indices that take into account only one single characteristic, that is, frequency, magnitude, or duration (Perkins, 2015). These limitations have led to a wide range of conclusions regarding heat wave trends across the United States and hinder our ability to understand the dynamic mechanisms responsible for heat waves (Chen & Li, 2017;Grotjahn et al., 2016). Further, while detection and attribution studies have found a slight anthropogenic influence on the severity of heat waves, the human influence on recent U.S. heat waves has not been sufficient to exceed modes of natural variability (Dole et al., 2014;Mazdiyasni et al., 2019;Min et al., 2013;Peterson et al., 2013).
To improve our understanding of heat wave response to a warming climate we use a novel combination of simple geometric concepts, a clustering algorithm, and indicator kriging to track the spatial evolution of individual heat wave events daily. This methodology allows us to develop a new spatial severity metric for heat waves that combines frequency, magnitude, duration, and areal extent to give a more complete characterization of the heat wave hazard.

Data
Given that surface station air temperature and humidity series often suffer from inhomogeneities combined with our focus on spatial analysis we utilize the gridded parameter-elevation regression on independent slopes model (PRISM) daily temperature dataset for 1981-2018 (4-km spatial resolution). PRISM combines station observations with a digital elevation model to construct daily gridded data for the contiguous United States (Daly et al., 2008;Thornton et al., 2014). PRISM has been highly quality controlled and is widely used in climatological studies. PRISM has also been validated against four reanalysis products with widespread agreement across several heat wave metrics (Schoof et al., 2017).
Heat waves can be defined as a sequence of days with temperature above a predefined high percentile threshold, typically between the 90th and 99th percentiles of daily maximum or mean temperature distribution (Perkins, 2015). Here we use the 95th percentile of local (grid cell) daily maximum temperature of the 1981-2018 climatology for each summer month (May-September) to define events. A variable monthly threshold is used to account for the seasonal cycle of heat. However, our methodology could be applied to any other high threshold or variable of interest. For example, consideration of the modifying influence of humidity on temperature is left for future work. The use of locally derived and monthly variable thresholds is necessitated by the clear evidence of acclimatization to heat extremes and local climatological variability that together negate the use of absolute and temporally/geographically static thresholds (Hondula et al., 2015). Both a minimum consecutive duration of 3 days above the 95th percentile threshold, and a three-day temporal independence to define separation between events, is set given epidemiological evidence of required multiday heat persistence and little to no heat wave effect after four or more days post event (Curriero et al., 2002;Perkins, 2015).

Methods
A regional analysis is undertaken where climate regions of the Continental United States are defined using National Oceanic and Atmospheric Administration's (NOAA) climatically consistent regions (Karl & Koss, 1984). Heat waves are tracked separately within each region, and results are presented regionally and aggregated across the entire Continental United States. The geometric concepts of convex hull, area, clusters, and solid area are used to track and quantify the evolution of heat waves at the surface (Keellings et al., 2018;Mcgarigal, 2015;Zick & Matyas, 2016) (Table 1). This approach allows for tracking of the development of the heatwave over the landscape, rather than solely looking at the size of area affected. First, heat waves are identified in each grid cell, and thus, a binary image of heat waves is created for each summertime day. A hierarchical clustering algorithm then identifies clusters of heat and delimits individual heat wave events based on Euclidean distance between clusters (Murtagh & Legendre, 2014). Heat waves can be considered to be generated at the synoptic scale (1,000 km) (Orlanski, 1975;Sheridan & Kalkstein, 2004). However, argument could be made that the generating scale of heat waves may vary regionally. Therefore, we regionally quantify the spatial length scale of heat waves, defined as the distance at which grid cells identified as in heat wave are no longer substantially correlated within that region. This is accomplished by employing indicator semivariograms similar to Touma et al. (2018) in their exploration of length scales of extreme precipitation events. Indicator semivariograms of grid cell pairings in which at least one cell is in heat wave are first created (i.e., pairs both in heat wave have value of 0 and dichotomous pairs have value 0.5), and then the mean of these values at lagged distances is taken. An exponential semivariogram model is then fitted, and the model's range (distance at which the model asymptotes and spatial autocorrelation ceases) then represents the heat wave spatial length scale. This is repeated on every identified heat wave day, and the median value is used as the regional heat wave spatial length scale. Clusters of grid cells in heat wave are identified as belonging to the same event if the Euclidean distance between them is less than or equal to the regionally defined heat wave spatial length scale. A convex hull is then drawn around all heat wave clusters belonging to the same event within each region on each day in order to delimit the boundary of the heat wave and to allow for it to be tracked from 1 day to the next based on spatial overlap. Each convex hull is given a unique identifier and is tracked daily from genesis through lysis based on spatial overlap with previously delineated convex hulls up to 3 days prior. Subsequently, the solid area of heat wave within the convex hull is calculated. The use of this spatial delineation and tracking methodology is essential as multiple smaller events, separated by the heat wave spatial length scale within a region, may subsequently merge into a single larger event or vice versa.
By delineating heat waves in this manner, with essentially no preconceived definition of severity and duration, we are able to determine event frequency (in each summer and/or regionally) and event duration. Event magnitude or intensity is also of great interest. To capture heat wave event magnitude and to inform the binary threshold exceedance, the average threshold exceedance among grid cells experiencing heat wave conditions on each day in each event is recorded (Schoof et al., 2017). As the heat wave hazard may be considered multidimensional including frequency, magnitude, duration, and areal extent we propose the Heat Severity and Coverage Index (HSCI) that, when used in conjunction with the methodology outlined here, gives a single value indicative of all combined dimensions of the heat wave hazard. The HSCI is conceptually similar to the Drought Severity and Coverage Index employed by the U.S. Drought Monitor in that it provides a convenient way to combine severity (e.g., frequency, magnitude, and duration) with areal extent and to aggregate regionally (Svoboda et al., 2002). The HSCI is defined as where m i is the average threshold exceedance (magnitude in degrees Celsius above threshold) and a i is the areal proportion of the region in heat wave, both of which are calculated for each day of the heat wave n.  Thus, the HSCI captures magnitude, areal extent, and event duration as the daily product of magnitude and areal extent is summed through the lifespan of the heat wave. The accumulated HSCI (AHSCI) can be calculated seasonally by summing all individual heat wave events from that season, thus capturing seasonal frequency, and can also be aggregated regionally. A conceptual example of the heat wave delineation methodology is illustrated in Figure 1 for an event in the south climate region during July/August 2012.

Results
We first explored the relationship between heat wave magnitude and area as well as duration and area. As heat wave events are tracked daily, we compared the maximum magnitude and duration of each event with its maximum areal extent. Figure 2 shows the significant positive correlation between these variables and indicates that events with a larger areal extent are likely to have a higher magnitude and last for a longer duration. Seasonal total heat wave area and seasonal heat wave magnitude were then examined regionally. Figures S1 through S4 show the respective seasonal total heat wave area, seasonal maximum heat wave magnitude, seasonal heat wave frequency, and seasonal mean heat wave duration series for each of nine climate regions, see inset map in Figure 3. Few of the nine regions show significant upward trends in any of these individual heat wave characteristics, but all regions exhibit pronounced interannual variability.
Next we calculated the AHSCI for each season and in each region. Figure 3 shows series of the HSCI aggregated seasonally and to each climate region. As each series exhibits substantial interannual variability, they have been smoothed using an exponentially weighted moving average in order to better visualize trends (Holt, 2004). The unsmoothed series for each region are shown in Figure S5 with significant linear trends indicated. The majority of regional AHSCI series show upward trends toward increased values of AHSCI with the exception of the Southeast. Most regions exhibit higher AHSCI values during the most recent two decades with relatively large increases in the South, Southwest, West, Central, and Northeast.
The observed AHCSI for the continental United States (Figure 4) shows pronounced interannual variability, and it captures the extreme heat wave years of 1988, 2003, and 2012 which are known to be highly impactful historical events (Lott & Ross, 2006;Lyon & Dole, 1995;NOAA, 2015). These events are estimated to have resulted in 30.2-57.7 billion USD damages with 454 direct and 5,000-10,000 excess deaths (1988), 6.2-7.5 billion USD damages with 35 direct and unknown excess deaths (2003), and 31.5-38.2 billion USD damages with 123 direct and unknown excess deaths (2012) (NOAA, 2015). Through combining all dimensions of the heat wave hazard including frequency, magnitude, duration, and areal extent the AHSCI captures these historic extreme events. Examining heat waves in this manner also indicates that extreme heat events have become more likely through the record as evidenced by the statistically significant (at the 0.05 significance level) upward trend in AHSCI (9.5 decade −1 ) for the continental United States.
However, this conclusion cannot be reached by examining heat wave hazards individually (see Figures S1 through S4), as while some previous studies have found regional upward trends in single heat wave characteristics such as frequency, areal extent, or magnitude, there has been only suggestive evidence of a slight increase in heat wave intensity (Habeeb et al., 2015;Mazdiyasni & AghaKouchak, 2015;Peterson et al., 2013;Russo et al., 2014;Schoof et al., 2017;Smith et al., 2013). This study is the first to examine heat waves holistically and to find a significant upward trend using a heat wave index that combines frequency, magnitude, duration, and areal extent.

Conclusions
This new methodology focuses on exploration of heat wave events spatially and allows the spatiotemporal evolution of individual events to be tracked. Previous studies have used a multitude of heat wave indices to measure events, but these indices often only capture individual characteristics such as magnitude or frequency, or confound frequency, magnitude, and duration by simply performing counts above a threshold, or by aggregating to the monthly or annual scales (Perkins, 2015). There has also been a lack of consideration of areal extent in heat wave metrics and when it is considered it is generally aggregated at the event or region scale. Schoetter et al. (2015) examined cumulative severity by employing an index of mean intensity, duration, and extent, across western Europe, but this was again aggregated to the event scale rather than tracked daily. The lack of consideration and understanding of spatiotemporal patterns of heat waves and their drivers has hampered planning and mitigation efforts for future public health impacts and response during heat wave events (Hattis et al., 2012) as fine scale spatial variation is important for determining health outcomes associated with heat waves (Kloog et al., 2015;Lee et al., 2016). These studies highlight the importance of studying the spatial characteristics of heat waves; to that end, the methodology outlined in this paper illustrates an approach to examine the spatial evolution of heat waves and shows statistical changes in holistic heat wave characteristics beyond those achieved with previously used methods. In addition, the spatial methodology allows for future exploration of additional spatial characteristics such as fragmentation and connectivity of heat within each event as well as the possibility to elucidate relationships between heat wave characteristics and geophysical variables.
Our findings raise the question of why the integration of individual heat wave characteristics, through the methods outlined here, reveals significant upward trends in heat waves while trends in those individual characteristics are limited? A possible explanation for this is changes in the persistence of heat within events but not necessarily accompanying changes in outright magnitude or duration. Our methodology tracks events daily, and the HSCI is computed as a sum across the entire event meaning that it captures not just the event's statistics such as mean or maximum magnitude but the accumulation of heat magnitude throughout the event.
This study is not without limitations in its methodology and results. As outlined in previous sections there are multiple possible definitions of heat waves and, while here we employ a commonly used approach with applications across societal and ecosystem impacts, we do not consider humidity which is known to decrease human comfort and intensify heat wave health impacts. A further limitation of the study is the approach of tracking heat waves separately within each climate region, that is, if a heat wave spans multiple regions it is tracked as a separate event in each region and constrained by the region's boundary. We employ this regional analysis in order to explore regional variations in heat wave characteristics, in recognition that heat waves may behave differently across climate regions and as a necessity of reducing computational expense in the application of the clustering algorithm across approximately 900,000 grid cells.
of natural modes of variability and human-induced climate change (DeGaetano & Allen, 2002;Gaffen & Ross, 1998;Meehl & Tebaldi, 2004). Furthermore, synoptic systems are associated with heat wave occurrence and of the synoptic weather types moist tropical and dry tropical are the most frequently correlated with heat waves (Sheridan & Kalkstein, 2004). It is interesting to note that across the United States the frequency of such weather types has increased in several regions (Senkbeil et al., 2017;Vanos et al., 2015).
General circulation models have projected more frequent, more intense, and longer duration extreme temperature events in North America across all Intergovernmental Panel on Climate Change scenarios throughout the rest of this century (Kunkel et al., 2013;Meehl & Tebaldi, 2004;Sillmann et al., 2013). With this in mind, greater clarity on the evolution of heat waves is essential as, despite the pervasive use of cooling technology in American homes, extreme heat continues to be linked to more annual fatalities than any other form of extreme weather in the United States and such impacts are likely to intensify in the future (Davis & Gertler, 2015;NWS, 2016). In the context of these other studies, our findings indicate a substantial increase in the heat wave burden across the United States during the last four decades that may be a presage of continued increases in heat waves to come.

Author Contributions
DK designed the study and performed the analysis. DK and HM wrote the manuscript. Both authors interpreted the results and improved the manuscript.