The critical role of humidity in modeling summer electricity demand across the United States

Cooling demand is projected to increase under climate change. However, most of the existing projections are based on rising air temperatures alone, ignoring that rising temperatures are associated with increased humidity; a lethal combination that could significantly increase morbidity and mortality rates during extreme heat events. We bridge this gap by identifying the key measures of heat stress, considering both air temperature and near-surface humidity, in characterizing the climate sensitivity of electricity demand at a national scale. Here we show that in many of the high energy consuming states, such as California and Texas, projections based on air temperature alone underestimates cooling demand by as much as 10–15% under both present and future climate scenarios. Our results establish that air temperature is a necessary but not sufficient variable for adequately characterizing the climate sensitivity of cooling load, and that near-surface humidity plays an equally important role.

A ccurate predictions of demand are a key challenge in electricity adequacy planning. Climate, technological, and socioeconomic factors are commonly used in predictive models of electricity demand 1,2 to ensure reliable planning and operation in the electricity sector by adequately balancing supply and demand. However, more frequent and intense climate extremes such as sustained heat waves 3,4 cause unanticipated changes in load 5 , challenging the reliability of electricity demand predictions. This poses a significant risk to the resilient operation of power systems 6 . Specifically, unanticipated higher demand for space conditioning and refrigeration during heat waves has led to rolling outages with serious socio-economic and public health consequences [7][8][9][10][11][12] . In light of the expected increase in frequency and intensity of climate extremes 3 , climate-induced outages are an increasing risk to the resilient operation of the electric power infrastructure [13][14][15][16][17] .
Quantifying the climate-sensitive portion of the electricity demand, referred to as the climate-demand nexus, hinges on effective characterizations of the key measures of heat stress. Heat stress measures refer to a combination of temperature and humidity that capture human discomfort levels and its ability to efficiently dissipate heat to avoid life-threatening conditions 9,10,18,19 . Heat stress is generally associated with high morbidity and mortality rates during heat waves. Despite the rich literature in electricity demand prediction 2,20-23 , little prior work has focused on exploring the climate sensitivity of demand with respect to different measures of heat beyond air temperature. While climate science research establishes air temperature as an incomplete measure of the surface heat content 9,10,18,24,25 , the majority of the existing research on climate-demand nexus use air temperature -or features derived from air temperature, such as cooling and heating degree days-as key predictors 2,26-28 . We address this gap by leveraging advanced data analytics to characterize the climate sensitivity of residential electricity demand as a function of a range of heat stress measures beyond air temperature, including dew point temperature, discomfort index (DI), heat index (HI), humidex, wet bulb temperature, and wet bulb globe temperature (see Methods). These measures are specifically chosen because they have been established in the literature as most effective in capturing heat stress 29 . However, no consensus has yet been reached on which measures are most predictive of the climate-demand nexus. We isolate the climate sensitivity of residential energy demand (see Methods), which has been established as the most climate-sensitive sector 15,16,[30][31][32][33][34] . Thus, our results primarily reflect changes in demand due to climate variability. We focus on summer months-that is, May, June, July, and August-as the majority of the US states peak during summer, and climate extremes such as heat waves occur in the summer 35 .
The central thesis in this paper is that air temperature alone underestimates the climate-sensitive portion of energy demand for cooling. To test this hypothesis, we first comprehensively assess the role of air temperature in capturing the climate sensitivity of summer electricity load at a US national scale; then, we identify the measures of heat stress that are most predictive of the climate-demand nexus and finally we quantify the likely underestimation of models based on air temperature alone under current and future climate scenarios.
We use monthly aggregated, state-level electricity consumption data across the contiguous United States, extracted from the Energy Information Administration (EIA) public reports over the years of 1990-2016 as well as population data from the US Census Bureau, and harmonized climate data from NCEP North American Regional Reanalysis (NARR) 36 . The electricity consumption data is carefully adjusted to remove trends associated with demographic and technological changes to isolate the climate effects on energy demand 37 (see Methods).
To quantify the effectiveness of air temperature alone in explaining the climate sensitivity of residential sector energy load, we develop two sets of models at the state level: air-temperatureonly models, which only use surface air temperature as a predictor, and selected-features models, where data-driven variable selection is performed to identify the key measures of heat stress considered in this study. Comparing the results of the airtemperature-only and selected-features models, we determine the influence of each measure of heat stress on the overall climate sensitivity of demand. The predictive models are developed using a state-of-the-art, stochastic, non-parametric Bayesian ensembleof-trees algorithm 38 (see Methods), which has been shown to outperform other climate-demand nexus models in terms of predictive accuracy 1,14,15,31,39 .
We conclude that the air-temperature-only models underestimate residential energy consumption, specially for future climate scenarios. With that, we show that near-surface humidity has an equally important role in electricity prediction as air temperature.

Results
Climate features identification. Our results reveal substantial improvements in prediction accuracy in 35 states for the cooling demand as a function of both temperature and humidity over the air-temperature-only model. Fig. 1a illustrates the key predictors of the climate-sensitive portion of electricity demand across the 48 US states as identified by the selected-features model. A darker shade represents states with higher energy consumption and the pie charts depict the selected climate variables for each state. We find significant variability in the identified key measures of heat stress across the United States, illustrating the importance of conducting region-specific analysis in characterizing the climate-energy-demand nexus. The most energyintensive states such as California and Texas, as well as many of the summer-peaking southern states (Supplementary Fig. 1 and Supplementary Note 1), benefit from the inclusion of measures of heat stress beyond air temperature-such as dew point temperature, wet bulb temperature, wet bulb globe temperature, and HI (see Methods)-in estimating the climate sensitivity of demand. This underscores the importance of accounting for near-surface humidity combined with temperature to better model human comfort levels and response to heat, which translates to residential cooling demand.
The highlighted states in Fig. 1b represent the extent in which the selected-features models outperformed air-temperature-only models in terms of predictive accuracy. Comparing the airtemperature-only models and selected-features models reveal that 35 states benefit from including measures of heat stress beyond air temperature, based on the comparison of root mean square error (RMSE) improvement in prediction accuracy (Fig. 1b). In addition to comparing model performances for the average consumption levels (i.e., 50th quantile predictions, Fig. 1b), we also compare the 90th quantile predictions (Supplementary Figs. 2 and 3). In this way, the effects of extreme heat events can be deduced from the higher quantiles of the aggregated electricity demand data. In other words, since heat waves typically coincide with the energy peaking summer months, the results for the 90th quantile RMSE improvement ( Supplementary Fig. 3) reflect months with more intense and frequent heat events. We observe significant model improvements, in terms of out-of-sample RMSE, for a number of high-energy-consuming states of the order of 8% (7.7%) in California, 8.5% (21.1%) in Texas, 8.6% (26.1%) in Illinois, and 7.1% (9%) in North Carolina, for the 50th quantile predictions (and 90th quantiles predictions-Supplementary Fig. 3). The observed improvements are substantial. For instance, in a state such as Texas, the selected-features model's average improvement of 8.5% in a given summer month such as August 2016 is equivalent to 1,498,968 MWh of energy consumption, which could sustain the residents of Austin-TX for more than 4 months [40][41][42][43] . In the same year and month, the 8% median increase in the state of California could sustain almost 1.5 million of Californians households 44 .
Sensitivity analysis under future climate scenarios. To quantify the extent to which the air-temperature-only models underestimate the climate sensitivity of demand under different climate scenarios, we use projected climate data extracted from five widely used Global Circulation Models (GCMs).
The GCM-extracted climate data are used to train state-level air-temperature-only and selected-features models to predict for the reference period of 1991-2010, and the future period of 2031-2050 (Methods). Since the objective here is to evaluate the changes in demand pattern as a result of future climate change (i.e., to quantify the relative increase over the periods, as opposed to absolute values), reference period data from the GCM is used instead of the historical data to remove bias induced by the use of data from different sources in the sensitivity analysis. In other words, input data generated from the same source (i.e., the GCMs) over the two time periods of reference (1991-2010) and future (2031-2050) are used in the developed predictive framework to track the relative increase in the climate-sensitive portion of demand as estimated by the air-temperature-only versus selected-features models.
The 1991-2010 reference period was chosen to match the period of the observed data as closely as possible and, at the same time, the same length as the future period (i.e., 20 years). While the common practice in climate science research is making   Fig. 1 Key climate predictors of cooling demand. a States are shaded according to total energy consumption intensity (scale bar goes from 1 to 48, with darker tones representing higher intensity states in terms of energy consumption). The pie charts illustrate the state-level key measures of heat stress for predicting the climate-sensitive portion of electricity demand. The majority of states have up to two most important features. It is important to point out that climate sensitivity is not an exclusively geographical characteristic, and considering the population density weight method for the variables calculation (see Methods), it is possible to expect neighboring states to not necessarily have the same selected climate variables. Non-summer-peaking states with many selected features (e.g., Washington (WA) and Oregon (OR)) reveal relative climate insensitivity in comparison (Supplementary Note 1). Note that the top energy-consuming states of California (CA) and Texas (TX) do not include air temperature as the key predictors of the climate-sensitive portion of demand. b Highlighted states show percentage improvement in predictive accuracy when comparing the selected-features models to the air-temperatureonly models (comparison based on out-of-sample RMSE). Source data are provided as a Source Data file. projections until 2100, that timeline is not appropriate in the energy sector, given the life span of the existing energy infrastructure [45][46][47] . The future period of 2031-2050 was chosen since energy infrastructure planning has the 2050 year as a comparable reference target in most energy projection and planning reports 48,49 .
The selected-features models project a significantly higher increase in the climate-sensitive portion of demand over time compared to the air-temperature-only models, when comparing the relative increase of each model (Fig. 2). The projections based on the selected-features models for the high-energy-consuming states of Texas, California, and Florida show a relative demand increase of about 12%, 8% and 10%, respectively, over the reference period (Fig. 2).   This comparison reveals that the underestimation in energy demand attributable to the air-temperature-only models is substantially more significant for the projections of higher demand values that are typically associated with periods of intense heat stress months-where the resilient operation of the grid is most challenged and unexpected demand spikes lead to supply inadequacy and thereby increase the risks of morbidity and mortality.

Discussion
The existing projections of climate-induced demand increase based on rising air temperatures alone ignore the fact that rising temperatures are associated with increased humidity-a lethal combination that increase morbidity and mortality in the absence of adequate cooling capacity during extreme heat events. We show that air temperature is a necessary but not sufficient measure to characterize residential space cooling demand during summer times. Humidity levels are also critical in capturing what truly reflects heat sensation for the human body 9 . Ignoring the role of humidity leads to underestimating the climate sensitivity of demand, challenging the resilient operation of power systems -especially under future warming scenarios where summertime energy production will be further constrained 50 . Inadequate supply to meet rising demand will have significant socioeconomic costs due to adverse health effects-particularly among the most vulnerable population-and warrants rapid and costly investments in energy infrastructure expansion as well as adaptation measures 26 .
We propose a data-driven framework to quantify the extent of underestimation of the climate-sensitive portion of cooling load, attributable to the use of air temperature alone as a measure of heat stress. Our results reveal a significant increase in predictive accuracy (of the order of 8% for high-energy-consuming states such as California and Texas) in characterizing the climatesensitive portion of demand through a more holistic inclusion of measures of heat stress. The results based on projecting demand under future climate scenarios are consistent, in that models based on air temperature alone show a systematic underestimation of the climate-sensitive portion of demand. The underestimation is particularly pronounced at the upper tail of the demand distribution, suggesting that the projected increase would be almost twice as large when considering measures of heat stress beyond air temperature in high-energy states such as California, Florida, and Texas.
It is important to note that our results pertain to only the climate-sensitive portion of the residential cooling demand. To project changes in total residential cooling demand under future climate scenarios, additional information regarding economic and population growth as well as technology and demographic changes need to be considered. What we wish to highlight here is that, under similar cooling technology, the existing projections based on air temperature alone substantially underestimate the anticipated demand increase under future warming scenarios, which could lead to inadequate and ineffective investments in capacity expansion and demand response programs as well as adaptation measures. Furthermore, while this study is focused on the residential energy demand, the methodology can be extended to other climate-sensitive sectors 51 .

Methods
Observed electricity consumption data. Monthly electricity consumption data was extracted from the EIA 52 over the years of 1990-2016 at a state level, with a total of 68 months per state (4 months per year for a total of 17 years). The source file includes electricity consumption data in Megawatt hours (converted to kWh), separated by sector: residential, commercial, industrial, transportation and others. The electricity consumption of the residential sector was the focus of the analysis presented in this paper. The residential electricity consumption data was detrended to isolate the climate effects on energy demand from techno-demographic changes such that increases in energy were not attributable to non-climatic factors such as technology shifts and population increase 16,37 . The residential energy consumption data in each state was first normalized using state-level population data-from the US Census Bureau 41 -to obtain a per-capita electricity demand, and then  where the total years, n years , range from 1990 to 2016, m denotes the month, and y the year. Then, an adjustment factor was calculated for each year as the sum of the percapita consumption per month divided by the yearly consumption E(y). This process was repeated for all states as follows: Finally, the adjusted energy consumption was calculated by dividing the monthly consumption by the adjustment factor: The detrended monthly electricity consumption data (described above) is referred to as observed electricity consumption data. This is a widely used method for trend-adjusting monthly aggregated regional data. For a comparative assessment of different de-trending methodologies please refer to Bessec and Fouquau 53 .
Observed climate data. The observed climate data was acquired at a sub-daily time scale for the period of 1990-2016 from the NCEP NARR 36,54,55 , aggregated at a monthly level to match the temporal scale of electricity consumption data, and weighted by population so as to give a greater weight to areas with higher population density when aggregating to the state level. The input climate data include air temperature, dew point temperature, DI, HI, humidex, wet bulb temperature, and wet bulb globe temperature. All climate variables use a combination of temperature, humidity, and pressure.
Dew point refers to the temperature in which air is saturated with water vapor, calculated in this paper using the equation below: where t d is the dew point temperature and t is air temperature in Celsius, and RH stands for relative humidity in percent (same RH as in the other equations in this section) 56 . Wet bulb temperature is the lowest temperature to which air can be cooled by water evaporation at a constant pressure. In this paper, it follows the equation: where π ¼ ðp=p 0 Þ 1=λ is the Exner function, used for scaling, λ is the inverse of the Poisson constant for dry air (3.504), p pressure (mb), p 0 the reference pressure (1000 mb); C is the freezing temperature in Kelvin, and T E the equivalent temperature, which is the moist potential temperature scaled by λ. For reviewing the complete derivation of the equation, refer to the Appendix A of refs. 29,57 . DI was developed in the 1950s to calibrate air conditioners and further adapted by the Israeli Defense Force as a main measure of heat stress: where T W stands for wet bulb temperature, RH is percentage relative humidity, and T C is the temperature in Celsius. DI is unitless 29,58,59 . Wet bulb globe temperature (sWBGT) is the heat stress measure in direct sunlight: where e RH = (RH/100)e sPa ; wet bulb temperature used here is unitless and widely used. e RH represents the vapor pressure in pascals and is calculated from RH in percentage and saturated vapor pressure 29,60 .
Humidex was developed for the Meteorological Service of Canada. It is a unitless 29,61 index, aiming to explain what the temperature feels like for the human body: HI is also called apparent temperature, describing what the temperature feels like to the human body when relative humidity is combined with air temperature: HI ¼ À 42:379 þ 2:04901523T F þ 10:14333127RH À 0:22475541T F RH À 6:83783x10 À3 T 2 F À 5:481717x10 À2 RH 2 þ 1:22874x10 À3 T 2 F RH þ 8:5282x10 À4 T F RH 2 À 1:99x10 À6 T 2 F RH; where (T F ) denotes the air temperature, and HI are measured in degrees Fahrenheit 29,62 . The data are considered under the emission scenario that has an end-of-century radiative forcing equal to 8.5 W m −2 , a Representative Concentration Pathway that is characterized by high greenhouse emission levels (RCP8.5) 15,64 . This global data is approximated by a 0.5°spatial resolution (~50 km) 65 and is processed from 1950 to 2099. The same variables were calculated for the projected data, namely, dew point temperature, DI, HI, humidex, wet bulb temperature, and wet bulb globe temperature. Two periods were extracted from the projected data, namely, the reference period of 1991-2010, and the future period of 2031-2050. There were two distinct stages in model development and analysis, namely, predictive model development using the observed data and sensitivity analysis using the projected data ( Supplementary Fig. 5). In the predictive model development stage, electricity consumption data, climate data, and population data were aggregated and normalized at the state level, as explained earlier, yielding monthly observed data for 48 states. Using these data sets, the development of airtemperature-only models and selected-features models involved three steps: for the selected-features model, variable selection was performed to identify the key predictors of cooling demand (details in the next section: BART (Bayesian Additive Regression Trees) algorithm and modeling process); then, the model hyperparameters were tuned (through cross-validation), and finally a ten-fold crossvalidation was performed to assess the predictive accuracy of the final models. For the air-temperature-only models, the same process was followed with the exception of the variable selection step, as there was only one independent variable. After developing these statistical models, the results were processed with the collection of in-sample (training set, 90%) and out-of-sample (test set, 10%) R 2 and RMSE values. The comparative assessments between the models were conducted based on the out-of-sample values of errors (RMSE results shown on Supplementary Figs. 2 and 3 and R 2 results shown on Supplementary Figs. 6 and 7).
In the sensitivity analysis stage, projected climate data from the five GCMs over the reference period (1991-2010) and future period (2031-2050) were used as inputs to the previously developed predictive models. This generated four sets of data, namely, air temperature, reference: estimates from the air-temperature-only model using the reference period data; air temperature, future: estimates from the air-temperature-only model using the future period; selected features, reference: estimates from the selected-features model for the reference period; and finally, the selected features, future, estimates from the selected-features model for the future period. This allowed for calculating the delta change over time between the same models (i.e., calculating the difference between selected-features models over the periods of 2031-2050 and selected 1991-2010). Similar analysis was also conducted on the 90th quantile of data, as electricity inadequacy during extreme heat episodes is typically associated with the tail of the demand distribution.
BART algorithm and modeling process. We leveraged a non-parametric Bayesian ensemble-of-trees algorithm to characterize the climate sensitivity of electricity consumption, since the algorithm was shown to outperform other climate-demand nexus models in terms of predictive accuracy 1,14,15,31,38,39 . The independent variables in the development of the BART models were heat stress measures, and the response variable was state-level electricity consumption. Models were trained and tested using a 10-fold cross-validation technique 66 .
The non-parametric Bayesian sum-of-trees model can be represented using the equation below 38,67 : Here, Y denotes the response variable, g denotes the regression tree function, T j denotes individual tree structures, M j denotes a set of parameters (e.g., mean values) associated with the tree nodes. BART is the sum of m total trees. The decision rules associated with each tree involve binary splitting of the predictor space. So for a given T and M, g(x; T, M) represents a tree with parameters M assigning the binary splits for the terminal nodes of T. ε is the stochastic irreducible error, assumed to be normally distributed with a constant variance 38 .
The BART algorithm uses a prior over the sum-of-trees parameters to regularize the trees, that is, to weaken the effect of individual trees with the goal for better prediction when adding all decision trees. A prior is imposed on all parameters, (T j , M j ) [∀j = 1, . . . , m] and σ. This regularization is made to keep individual trees effects from being unduly influential over the entire sum of trees. Without a regularization, large tree components would overwhelm the sum of trees and would diminish the benefits of leveraging additive tree structures. The prior regularization can be mathematically represented as: where tree components T j and M j are assumed to be independent of each other and of σ.
As mentioned earlier, to assess the role of air temperature in capturing the climate sensitivity of electricity consumption compared to other measures of heat stress, two models were developed for each state: the air-temperature-only model, and the selected-features model. Air-temperature-only models were constructed with a single input variable of surface temperature and electric demand as the response variable. To develop the selected features a data-driven variable selection algorithm was harnessed. Variable selection: specifically, the variable selection algorithm involved identifying the variables most frequently used as decision/ splitting rules in the sum-of-tree model 38 . The calculation was then averaged over 100 model iterations to achieve stable estimates throughout the sum-of-trees distribution. This is a reasonably robust approach for variable selection using the Bayesian ensemble-of-trees algorithm 38,67 .
Pseudo-code. Feature selection is performed for 48 states before running BART and out-of-sample RMSE is averaged over 20 models (cross validation performed at line 3). Code was performed in parallel cores.

Data availability
All data sets used in this study are publicly available from the referenced sources.