## Abstract

The increasing frequency of very high summertime temperatures has motivated growing interest in the processes determining the probability distribution of surface temperature over land. Here, we show that on monthly time scales, temperature anomalies can be modeled as linear responses to fluctuations in shortwave radiation and precipitation. Our model contains only three adjustable parameters, and, surprisingly, these can be taken as constant across the globe, notwithstanding large spatial variability in topography, vegetation, and hydrological processes. Using observations of shortwave radiation and precipitation from 2000 to 2017, the model accurately reproduces the observed pattern of temperature variance throughout the Northern Hemisphere midlatitudes. In addition, the variance in latent heat flux estimated by the model agrees well with the few long-term records that are available in the central United States. As an application of the model, we investigate the changes in the variance of monthly averaged surface temperature that might be expected due to anthropogenic climate change. We find that a climatic warming of 4°C causes a 10% increase in temperature variance in parts of North America.

## 1. Introduction

As the global mean temperature rises, the impacts of extreme summertime temperatures are expected to grow more severe (Kirtman et al. 2013). An increase in seasonally averaged temperatures has already led to more record-breaking high temperatures around the world (Rhines and Huybers 2013; Diffenbaugh et al. 2017), but whether anthropogenic climate change will drive changes in the probability distribution of summertime temperatures is an unsolved question. Some authors argue that only a shift in the underlying temperature probability distribution could generate such extreme events as the record-breaking 2010 European heat wave (Schär et al. 2004; Barriopedro et al. 2011). Others argue that, while unprecedented in the past 600 years, this extreme seasonal anomaly can be attributed to gradual warming without invoking a changed probability distribution (Tingley and Huybers 2013). Climate models project increased summertime temperature variance in simulations where greenhouse gas emissions are not curtailed, but underlying biases in model representations of the historical record along with highly parameterized relationships between land surface fluxes, soil moisture, and 2-m air temperature make the accuracy of such predictions difficult to evaluate (Seneviratne et al. 2010; Mueller and Seneviratne 2014). Given the complexity of global climate models and their large biases in contemporary summertime temperature variance (Vargas Zeppetello et al. 2020), developing alternative modeling approaches is a key strategy for understanding and improving the projections of future temperature variance in the context of anthropogenic climate change.

In this paper, we develop a diagnostic model for summertime temperature variance written in terms of monthly anomalies in shortwave radiation and precipitation. Since observations of shortwave radiation, precipitation, and temperature are available from 2000 to 2017, we can use the model to estimate temperature variance during this period and compare our estimate to the observed values. This diagnostic model is unique in that it depends only on well-observed quantities; other novel modeling approaches (e.g., Schwingshackl et al. 2018; Vargas Zeppetello et al. 2020) use correlation analysis to relate temperature variability to various components of the surface energy budget available only within a climate modeling context. In contrast, our diagnostic model can be used to estimate summertime temperature variance in any dataset that includes information on monthly radiation and precipitation variability.

First, we develop the diagnostic model from considerations of the land surface energy and water budgets. Next, we use the model to estimate the variance in monthly average temperature during summer and compare it to the observed temperature variance over North America from 2000 to 2017. We also show the diagnostic model’s estimate of latent heat flux variance and compare it to that from the longest observed record of latent heat flux in the United States. Finally, we use the model to explore how changes in various climatological mean variables expected due to increasing anthropogenic greenhouse gas emissions may affect summertime temperature variance. A summary and conclusions follow.

## 2. The diagnostic model

On monthly time scales, we assume that the land surface energy and water budgets are in equilibrium. These two equilibrium budget equations are

All terms in Eq. (1) are given in W m^{−2}, while all terms in Eq. (2) are given in kg H_{2}O m^{−2} s^{−1}. The term $F$ is the net shortwave radiation incident at the land surface, while *F*_{LW} is the net surface longwave radiation flux. The terms *LE* and *H* are the turbulent fluxes of latent and sensible heat respectively, while *G* is the flux of energy downward into the soil column. Also, $P$ is precipitation, *R* and *I* are the surface runoff and infiltration moisture fluxes respectively, *E* is the net evapotranspiration, and *L* is the latent enthalpy of vaporization. Primed quantities denote departures from respective monthly means for June, July, and August over the 2000–17 period. For example, the net shortwave radiation anomaly in the *i*th month of the *j*th year can be written as

where $Fi\xaf$ denotes the climatological monthly mean for the *i*th month over the 2000–17 period, $Fi\xaf=(1/18)\u2211jFi,j$. Using this definition of a monthly anomaly, Fig. 1a shows the summertime 2-m air temperature variance $\sigma 2\u2061(T)=(1/N)\u2211i,j\u2061(Ti,j\u2032)2$ between 2000 and 2017 computed from the interpolated weather station data presented in Willmott and Matsuura (2001). Since our record comprises 18 sets of three summer months, *N* = 18 × 3 = 54. In the remainder of this paper, we will drop the subscripts on the anomaly terms for clarity.

### a. Linear response surface fluxes

The sum of monthly net longwave, sensible heat, and ground heat flux anomalies is assumed to be linearly proportional to temperature perturbations, thus

Here, *ν* (in W m^{−2} K^{−1}) is the reciprocal of the response of 2-m air temperature *T*′ to a radiative forcing $F\u2032$ in the absence of evapotranspiration anomalies [see Eq. (1)]. While it is common practice to write anomalies in upward surface longwave radiation, sensible heat flux, and ground heat flux as linear functions of near-surface temperature, recent studies (e.g., Schwingshackl et al. 2018; Vargas Zeppetello et al. 2019b) have demonstrated that downward longwave radiation anomalies can also be written as linear functions of surface temperature anomalies.

The use of 2-m air temperature, rather than land surface temperature, introduces some ambiguity into the value of *ν*, as sensible heat flux has been shown to impact land surface and 2-m air temperature differently on submonthly time scales (Gallego-Elvira et al. 2016; Panwar et al. 2019). However, the extremely high correlation between monthly anomalies in land surface and 2-m air temperature found in reanalysis and global climate models suggests treating the land–atmosphere interface as a phantom layer that is as equally well characterized by the surface temperature as it is by the 2-m air temperature (Vargas Zeppetello et al. 2020). Because reliable observations of 2-m air temperature span the globe, we will accept the ambiguity associated with the *ν* parameter in order to derive a closed model for temperature variance that can be compared to observations.

Similar to the fluxes in Eq. (4), the sum of runoff and infiltration anomalies is assumed to be linearly proportional to soil moisture perturbations, thus:

In Eq. (5), “fractional surface saturation” *m* is a unitless number between zero and one that designates the fraction of available pore space in the evapotranspiration-accessible portion of the soil column that is occupied by liquid water. The relationship between fractional surface saturation and a particular mass flux of liquid water depends on soil properties that are not known at continental scales. Despite this knowledge barrier, hydrological models of varying complexity feature high correlations between runoff, infiltration, and incident precipitation on seasonal time scales (Haddeland et al. 2011; Oki and Kim 2016). To ensure proper scaling between these three mass fluxes and the fractional surface saturation, at each grid cell we set the “surface moisture capacity” *μ* (in kg m^{−2} s^{−1}) to be

where $\sigma \u2061(P)$ is the local summertime standard deviation in monthly averaged precipitation and *α* is a unitless parameter that controls the mass of liquid water required to effectively change the soil’s fractional saturation *m*. High *α* could correspond, for example, to a shallow soil column where infiltration syphons off a large portion of incident precipitation without changing the soil’s effective saturation. In addition to the depth of the evapotranspiration-accessible portion of the soil column, *α* is dependent on the soil type, porosity, and topography. Since observations of $\sigma \u2061(P)$ exist (see section 2c), *α* is the real tunable parameter embedded in Eq. (5), but we will retain *μ* in the rest of our equations for simplicity.

In wet regions, where mean precipitation and $\sigma \u2061(P)$ are both generally high, Eqs. (5) and (6) imply that anomalies in runoff and infiltration are generally large for a given $P\u2032$. Higher runoff and infiltration fluxes in wet regions have been inferred from satellite observations of soil moisture (McColl et al. 2019) and in soil chamber experiments dating back to the earliest examinations of flow through porous media (Richards 1928). In dry regions where $\sigma \u2061(P)$ is generally small, Eqs. (5) and (6) imply that soil moisture fluctuations are not associated with large changes in the lateral or vertical liquid water fluxes. This is consistent with modeling and catchment analyses that show low runoff fractions in dry soils (Haddeland et al. 2011).

### b. Evapotranspiration

We write total evapotranspiration as

In Eq. (7), *ρ*_{a} (kg air m^{−3}) is the density of air, *r*_{s} (s m^{−1}) is the “bulk surface resistance” parameter, and *V* (kg H_{2}O kg air^{−1}) is the atmospheric water vapor demand written as the difference *q*_{s}(*T*) − *q*, where *q*_{s} is the saturation specific humidity at the 2-m air temperature *T*, and *q* is the boundary layer specific humidity. The first-order terms in a Taylor expansion of Eq. (7) are

where barred terms indicate summertime mean values. In Eq. (8), *γ* [kg H_{2}O (kg air)^{−1} K^{−1}] is a linearized value of the derivative *dq*_{s}/*dT* given by the Clausius–Clapeyron equation evaluated at the summertime mean 2-m air temperature $T\xaf$. We have assumed that anomalies in specific humidity are small compared to those in *q*_{s} on monthly time scales—as is true in limited observations available for comparison (van Heerwaarden et al. 2010; Vargas Zeppetello et al. 2019a).

### c. Diagnostic equations

where we have substituted

as the climatological mean potential evapotranspiration, or the mean evapotranspiration $E\xaf$ expected for $m\xaf=1$, or saturated soils. Note that *β* increases exponentially with $T\xaf$ due to the dependence of saturation specific humidity on the Clausius–Clapeyron relationship. Combining Eq. (9) with Eqs. (1) and (8), we obtain

where *ζ* = (1 + *μ*/*β*)^{−1} ⊂ [0, 1] is a dryness index and Γ^{−1} is the “moist surface climate sensitivity”:

Maps of Γ and *ζ* are shown in Fig. 2. The inverse of the first term in Eq. (12), *v*^{−1} (K W^{−1} m^{2}), can be considered the “dry surface climate sensitivity” [see Eq. (4)]. The second term in Eq. (12) is associated with the climatological mean latent heat flux [see Eq. (8)] that renders drier soils (smaller $m\xaf$ and larger *ζ*) more sensitive to forcing than wet soils.

We have now obtained an equation for *T*′ in terms of monthly anomalies $F\u2032$, $P\u2032$, and climatological values $m\xaf$ and $V\xaf$ that manifest themselves in Eq. (11) through the Γ and *ζ* terms. Observations of climatological fractional surface soil moisture $m\xaf$ and atmospheric water vapor demand $V\xaf$ are not available at the continent scale. For climatological soil moisture $m\xaf$, we normalize the summertime precipitation climatology $P\u2061(x)\xaf$ at every grid cell by the maximum monthly averaged precipitation over North America max $\u2061(P\xaf)$; hence $m\u2061(x)\xaf=P\u2061(x)\xaf/max\u2061(P\xaf)$.^{1} Finally, for each point in space, we calculate climatological atmospheric water vapor demand from the observed summertime mean surface temperature, and the summertime mean specific humidity from ERA5 reanalysis’ lowest atmospheric layer (Copernicus Climate Change Service 2017).

Using Eq. (11), we arrive at our equation for temperature variance written as the sum of three terms:

We will refer to $\sigma 2\u2061(F)$, $\sigma 2\u2061(LP)$, and $F\u2032LP\u2032\xaf$ as the “forcing components,” each in W^{2} m^{−4}. Figures 1b–d shows the observed forcing components for the 2000–17 period. Note that everywhere $F\u2032LP\u2032\xaf<0$ due to the negative correlation between shortwave radiation and precipitation. Observations of summertime radiation are from the CERES satellite (CERES Science Team 2000), while observations of precipitation are from interpolated weather station data presented in Willmott and Matsuura (2001). Equation (13) can be used to estimate temperature variance in any dataset that includes surface shortwave radiation and precipitation; the appendix contains maps of the summertime temperature variance (Fig. A1) and forcing components (Figs. A2–A4) in the Northern Hemisphere midlatitudes from the observations, ERA5, and an ensemble of global climate model simulations that we will analyze to complement our analyses presented in section 3.

We can use Eqs. (8)–(11) to write an equation for the variance in latent heat flux. Written in terms of the forcing components:

Diagnostic Eqs. (13) and (14) constitute our model and contain three tunable parameters *ν* [Eq. (4)], *μ* [Eq. (5)], and *r*_{s} [Eq. (7)]. Observations of these quantities do not exist; we assume that the three parameters are constant across the model domain. By tuning the model parameters to obtain good overall agreement with the observed temperature variance over the North American domain, we obtain *ν* = 14 W m^{−2} K^{−1}, $\mu =5\sigma \u2061(P)kgm\u22122s\u22121$, and *r*_{s} = 75 m^{−1} s. The values of Γ and *ζ* depend on these parameters and the summertime mean values of $m\xaf$ and $V\xaf$.

## 3. Results

### a. Temperature variance

We begin with an evaluation of Eq. (13): can this diagnostic model accurately estimate the summertime temperature variance in the midlatitudes? The assumption of constant parameter values ignores regional variations in roughness, vegetation phenology, and other surface properties. Nonetheless, in Fig. 3 we show that over the North American domain the temperature variance estimated by the model using forcing components from observations (shown in Figs. 1b–d) agrees well with the observed pattern of temperature variance. Figure 3a shows a map of the ratio of the temperature variance estimated by the model using the observed forcings to the observed temperature variance *σ*^{2}(*T*)_{Model}/*σ*^{2}(*T*)_{Obs.}, while Fig. 3b shows histograms of this variance ratio for various sets of grid cells described below. The structure of Eq. (13) allows us to decompose the temperature variance into three terms associated with each forcing component. Maps of the variance computed according to Eq. (13) and each term associated with this decomposition are shown in Fig. 4. One notable aspect of this decomposition is that the component of temperature variance associated with precipitation forcing is similar to that found in Koster et al. (2015) using different methods.

The model agrees well with the observed pattern of temperature variance over North America (the training dataset); the solid black histogram is centered on unity, indicating good overall agreement, and the diagnostic model’s estimated temperature variance is off by more than a factor of 2 in fewer than 15% of grid cells. While the three spatially invariant parameters listed above were tuned to achieve this agreement, extending the model throughout the Northern Hemisphere midlatitudes (between 30° and 75°N) also shows a good fit to observations (the dashed black histogram in Fig. 3b); 78% of grid cells are within a factor of 2 of the observed temperature variance and the distribution is centered on unity with no obvious bias toward over or underestimating temperature variance. However, there are several specific regions where the model’s basic assumptions are violated, resulting in spurious estimates.

The model tends to overestimate temperature variance in high-elevation regions. This is likely due to our simple formulation of runoff flux, which does not account for topography and thus results in an overestimation of temperature variance associated with precipitation. Removing high-elevation regions from the set of points analyzed in Fig. 3b, we obtain the blue histogram that has notably fewer points with overestimates of temperature variance. Problems associated with our runoff parameterization are also likely responsible for the overestimate of temperature variance around the Gulf of Mexico. Monthly averaging in the Gulf Coast means that short time scale variability associated with tropical cyclones is not explicitly accounted for and creates a high bias in soil moisture variability that is reflected in the overestimate of temperature variance (see Fig. 3a). With a more sophisticated runoff parameterization, the large precipitation anomalies (primarily associated with tropical cyclones) would have less influence on soil moisture due to high runoff rates driven by these systems.

Underestimates of temperature variance are also present in our evaluation exercise, particularly in coastal regions such as the west coast of the United States and far northern Canada. In these regions, remote forcing due to thermal advection (and possibly associated with summer sea ice variability at high latitudes) plays a significant role in the temperature variance (Holmes et al. 2016). The red histogram in Fig. 3b shows the ratio of variances in grid cells where the observed temperature variance exceeds 2°C^{2}. Since many of these grid cells are in the high latitudes (see Fig. 1a), where summertime sea ice variability and thermal advection likely contribute to high temperature variance, we attribute the systemic underestimate of temperature variance in this subset of grid boxes to the absence of thermal advection in our model formulation. If we only consider high temperature variance regions, south of 50°N, however, our model is a fair fit to observations: 70% of points lie within 50% of the observed temperature variance. This suggests that even in high temperature variance regions in the continental United States, our model captures the essential physics governing the distribution of land surface temperatures. While some studies (e.g., Schneider et al. 2015) have noted the important role of thermal advection in regulating midlattitude temperature variability in wintertime, our work aligns with the results of Holmes et al. (2016), who argue that in most continental regions in summertime, thermal advection is not a major contributor to local summertime temperature variance. Another region where the model underestimates temperature variance is the Ohio Valley (south of the Great Lakes). Mueller et al. (2017) quantified the influence of agriculture on temperature extremes in this region; the neglect of vegetation changes during summertime may contribute to the model’s underestimate of temperature variance.

The histograms detailed above testify to the skill of our model [Eq. (13)] to simulate the observed temperature variance. Analogous maps and histograms that show the diagnostic model’s ability to reproduce the temperature variance in ERA5 and an ensemble of global climate models are shown in the appendix (Figs. A5 and A6).

The temperature variance and forcings in the three datasets shown in Figs. A1–A4 (observations, ERA5, and the global climate model ensemble) differ substantially from one another, but the same values of three tunable parameters (bulk surface resistance *r*_{s}, dry surface climate sensitivity *ν*, and surface moisture capacity *μ*) are used in each realization of the model and lead to accurate estimates of each product’s temperature variance. This suggests that the model’s skill does not stem from overtuning; Eq. (13) captures the basic physics underlying the relationships between temperature variance and the surface budgets of energy and moisture. Also, with the exception of the mean surface moisture $m\xaf$, variations in soil hydrological parameters and exchange coefficients are unimportant in capturing the spatial pattern and amplitude of monthly temperature variance. This result is reminiscent of a recent study by Byrne and O’Gorman (2018) that accurately captured observed trends in land surface relative humidity using only a single land surface dryness parameter.

### b. Latent heat flux variance

We have shown that our model reproduces the observed temperature variance given the observed forcing. How does the pattern of temperature variance relate to the variance in latent heat flux over the same period? Figure 5 shows the latent heat flux variance estimated by Eq. (14). Similar to the temperature variance, the latent heat flux variance can be decomposed into three forcing terms. While each forcing component makes a comparable contribution to the temperature variance, precipitation forcing is overwhelmingly responsible for the variance in latent heat flux, leading to a maximum in latent heat flux variance in the central United States. Maps of each component are shown in Fig. 6.

Direct global observations of latent heat flux variance do not exist. Long records of latent heat flux covering the entire 2000–17 period are only available from three energy balance Bowen ratio stations at three Atmospheric Radiation Measurement sites within the Southern Great Plains (SGP) field campaign (Xie and Gaustad 1993). By averaging the output from the three sites and removing monthly means in a manner analogous to Eq. (3), we estimate the latent heat flux variance in the SGP domain (shown by the red circle in Fig. 5) as 283 W^{2} m^{−4}. This compares favorably to the latent heat flux variance calculated with Eq. (14) using the observed forcing (see inset in Fig. 5). For comparison, the latent heat flux variance over the same period and region from ERA5, a new reanalysis product, is shown in the inset to Fig. 5; it is 35% less than the variance observed in the SGP record. That our model reproduces the observed latent heat flux variance at the SGP is independent evidence that the parameter values (chosen to match our model to the observed temperature variance) are physically realistic.

### c. Temperature variance in a changing climate

The plasticity of the diagnostic model allows us to explore the sensitivity of temperature variance to changes in selected state variables. We specify changes in state variables that are expected with a 4°C global average warming such as a decrease in land relative humidity of 4% (Byrne and O’Gorman 2016) and a decrease in soil moisture of 10% (Berg et al. 2017). We emphasize that our results should not be interpreted as projections because we use spatially uniform mean-state changes in this analysis. Rather, our results serve to highlight the importance of different aspects of land–atmosphere interaction on local temperature probability distributions.

Increasing the climatological mean surface temperature by 4°C while holding relative humidity constant increases climatological mean atmospheric water vapor demand $V\xaf$ by way of increasing the surface *q*_{s} (and its derivative *γ* = *dq*_{s}/*dT*) relative to *q*. As a result, the temperature variance increases everywhere (Fig. 7a) due to an increase in *β*, which increases the moist surface climate sensitivity Γ^{−1} [see Eqs. (10)–(12)]. Including a 4% reduction in relative humidity amplifies the increase in temperature variance associated with a 4°C warming (cf. Figs. 7a,b) by further increasing the atmospheric water vapor demand $V\xaf$. Importantly, these changes are largest in wet regions in the eastern United States, where evaporative cooling is more sensitive to background atmospheric demand $V\xaf$ than to the availability of soil moisture $m\xaf$ (Seneviratne et al. 2010; Vargas Zeppetello et al. 2019a). In these regions, the combination of a 4°C increase in temperature and a 4% decrease in relative humidity causes temperature variance to increase by up to 25% (Fig. 7b).

Climate models simulate a decrease in climatological surface soil moisture $m\xaf$ associated with increased summertime evapotranspiration in greenhouse warming scenarios (Berg et al. 2017). Figure 7c shows the change in temperature variance associated with a 10% reduction in soil moisture (a typical value for the end of this century in the RCP8.5 emission scenarios). A reduction in mean soil moisture increases temperature variance by reducing the latent heat flux anomalies associated with a nominal temperature perturbation, thus increasing the moist surface climate sensitivity Γ^{−1} [see Eq. (12)]. The increase in temperature variance due to a climatological mean soil drying of 10% is less than that due to a 4°C mean warming if the projected decrease in relative humidity is included (cf. Figs. 7b,c). The regions with the largest increases in temperature variance are coincident with maximum values of $m\xaf$, where the 10% reduction in soil moisture constitutes the maximum absolute decrease in soil moisture.

Colman (2015) showed that anthropogenic climate forcing is projected to cause a reduction in low clouds throughout the midlatitudes. Under the assumption that a decrease in average cloudiness will be accompanied by a decrease in the shortwave radiation variance at Earth’s surface, Fig. 7d shows the impact of a 10% decrease in $\sigma 2\u2061(F)$ at each grid cell in the domain. This results in a *σ*^{2}(*T*) reduction everywhere, but the signal is particularly strong in the high latitudes, where radiative forcing is largest, and in the dry American West, where radiative forcing is the largest contributor to temperature variance due to the lack of available soil moisture for evapotranspiration (see Fig. 4).

Our results suggest that a mean warming combined with an atmospheric drying could dramatically amplify temperature variance—particularly in wet regions where soil moisture is plentiful and evapotranspiration responds primarily to evaporative demand (Fig. 7b). Swann et al. (2016) demonstrated that, on a global scale, CO_{2} climate forcing in an ensemble of models may not translate into meaningful changes in mean evapotranspiration because plants compensate for the enhanced atmospheric demand by closing their stomata in response to higher levels of atmospheric CO_{2}. In terms of our diagnostic model, this corresponds to an increase in bulk surface resistance *r*_{s}. Figure 7e shows the impact on temperature variance of a 10% reduction in surface resistance everywhere in the domain (*r*_{s} = 82.5 vs *r*_{s} = 75). As expected, the impact of this change is limited to the regions where latent heat flux variance is largest (see Fig. 5). Interestingly, the sign of the change is negative: increasing surface resistance leads to *decreasing* summertime temperature variance. This change is somewhat surprising because increasing surface resistance causes greater mean temperatures. However, increasing surface resistance also reduces evapotranspiration sensitivity to any forcing at all time scales. Thus, increasing surface resistance damps temperature variance in the central United States, where we have shown that latent heat flux variations contribute most to overall temperature variance. This effect is not trivial; Fig. 7f shows that combining the impacts of reduced soil moisture and increased surface resistance results in a muted change over the central United States due to the tug of war between the two large but opposing impacts on latent heat flux in this region.

## 4. Conclusions

In this paper, we use the land surface budgets of energy and liquid water to derive a simple diagnostic model of the variance in monthly averaged temperature and latent heat flux during summertime. The variances are linear functions of local shortwave radiation variance, precipitation variance, and the covariance between monthly anomalies in shortwave radiation and precipitation. Our approach is uniquely independent of any particular modeling platform; we can use observations, reanalysis, and global climate model output to estimate a particular dataset’s temperature variance. Using the observed forcing values, our model reproduces the summertime temperature variance in observations, ERA5, and a global climate model ensemble for the 2000–17 period with considerable skill. We also provide an estimate of summertime latent heat flux variance over North America. Unfortunately, an observed record of latent heat flux that is long enough to estimate the monthly variance is only available at one location: the ARM Southern Great Plains site. Our model’s estimate for the latent heat flux variance in this region agrees favorably with the observed record and is closer to the observed Southern Great Plains site value than a state-of-the-art reanalysis product (ERA5). One drawback of the model is that it does not include advective terms or feedbacks between the land surface and atmospheric forcings. Future work in land–atmosphere interaction will focus on short time scale feedbacks that could complicate the relationship between atmospheric forcings and surface fluxes that were assumed independent in our model.

We used the model to explore the impacts of climate change on surface temperature variance. We find that, in isolation, 4°C warming causes an increase in temperature variance over all of North America. Importantly, the increase in temperature variance is extremely sensitive to changes in relative humidity of the near-surface air; including a 4% decrease in relative humidity (projected to occur with a 4°C increase in temperatures; see Byrne and O’Gorman 2016) increases the temperature variance by up to 25% in some climatologically wet regions of the United States compared to the 2000–18 value. When changes in both temperature and relative humidity are accounted for, temperature variance increases by more than 20% in some parts of the eastern United States. We also find that the increase in temperature variance associated with the projected decreases in soil moisture could be compensated for by increased stomatal resistance in response to higher atmospheric CO_{2}, but this estimate is poorly constrained due to the complexities associated with quantifying the stomatal response to carbon emissions.

In simulations of the twentieth century, climate models display large biases in summertime temperature variance on monthly time scales (Vargas Zeppetello et al. 2020). These biases raise serious questions concerning the reliability of the model-projected changes in temperature variance driven by anthropogenic climate change. The present study presents a tool that can be used to identify the sources of bias in simulated temperature variance by quantitatively assessing the relative contributions of biases in the forcings (shortwave and precipitation variance) and biases in the mean state (temperature, humidity, and soil moisture) to biases in the simulated temperature variance. In this way, the model can be used to identify pathways for improving the models, which presumably will lead to more reliable projections of changes in temperature variance driven by climate change. The biases in climate models have led investigators evaluating the impact of climate change on food production and human heat stress to make the conservative assumption that the variability of summertime temperature will remain constant as the climate warms (e.g., Tigchelaar et al. 2018). Our results suggests that this assumption ignores potentially dramatic increases in temperature variance associated with climate change and provides a framework for quantifying how large those changes are likely to be so that they can be integrated more easily with impact assessments.

## Acknowledgments

This paper contains data from the Copernicus Climate Change Service, the CERES satellite, the NOAA ESRL, and the CMIP5 data archive. We are grateful to all who made these data available to the public. We thank three anonymous reviewers, Steve Warren, and the EcoClimate group at the University of Washington for helpful feedback on the manuscript. LRVZ was supported by a NSF GRFP Fellowship and DSB was supported by a grant from the Tamaki Foundation.

### APPENDIX

#### Comparison between Observations, Reanalysis, and Global Climate Models

We present global maps of midlatitude Northern Hemisphere temperature variance (Fig. A1) and forcings (Figs. A2–A4 ) for observations, ERA5, and a global climate model ensemble (GCME) taken from the CMIP5 data archive. ERA5 runs from 1979 to 2019; for consistency with the available observations we used the 2000–17 period to calculate all quantities relevant for the diagnostic model. The GCME is composed of scenarios where the atmospheric CO_{2} concentration rises by 1% each year from a preindustrial control baseline in six models listed in Table A1; we use the first 18 years of these runs. For the GCME temperature variance, we calculated an average of each forcing from each of the six models considered, then input them into the diagnostic model to obtain the estimate of temperature variance. We then compare this estimate to the average of the temperature variances simulated by the six models.

The temperature variance and radiative forcing in ERA5 are nearly identical to the observations, while the GCME has high biases in both quantities. We have identified biases in incident shortwave radiation variability as an important contributor to climate model temperature variance biases in other work (Vargas Zeppetello et al. 2020). Despite the biases in shortwave radiation variance, the precipitation forcing and covariance terms from the GCME compare favorably to those from both ERA5 and the observations.

In Figs. A5 and A6, we show the skill of the diagnostic model as the ratio of temperature variance estimated from the diagnostic model driven by the relevant forcings from observations, reanalysis, and the GCME to the temperature variance found in those three datasets. We use the same values for the parameters *ν*, *μ*, and *r*_{s} in each realization and obtain good agreement between the model estimate and the temperature variance in all three products.

## REFERENCES

*Science*

*Geophys. Res. Lett.*

*J. Climate*

*Proc. Natl. Acad. Sci. USA*

*J. Geophys. Res. Atmos.*

*Proc. Natl. Acad. Sci. USA*

*Geophys. Res. Lett.*

*J. Hydrometeor.*

*J. Climate*

*Climate Change 2013: The Physical Science Basis*, T. F. Stocker et al., Eds., Cambridge University Press, 953–1028.

*Front. Earth Sci.*

*J. Hydrometeor.*

*Geophys. Res. Lett.*

*J. Climate*

*Terrestrial Water Cycle and Climate Change: Natural and Human-Induced Impacts*, Q. Tang and T. Oki, Eds., John Wiley & Sons, 1–16, https://doi.org/10.1002/9781118971772.ch1.

*Geophys. Res. Lett.*

*Proc. Natl. Acad. Sci. USA*

*J. Agric. Res.*

*Nature*

*J. Climate*

*Geophys. Res. Lett.*

*Earth-Sci. Rev.*

_{2}reduce estimates of climate impacts on drought severity

*Proc. Natl. Acad. Sci. USA*

*Proc. Natl. Acad. Sci. USA*

*Nature*

*J. Hydrometeor.*

*J. Climate*

*Geophys. Res. Lett.*

*J. Climate*

## Footnotes

^{a}

Current affiliation: Department of Earth and Space Science, University of Washington, Seattle, Washington.

^{1}

We note that mean soil moisture does indeed scale linearly with mean precipitation in both ERA5 and the global climate model ensemble we analyzed in this work. The scaling factor we chose was $max\u2061(P\xaf)$, but the results were insensitive to changes in this value.