Landscape pattern is critical for the moderation of thermal extremes

Abstract

Temperature is highly variable across space and time at multiple scales, shapes landscape pattern, and dictates ecological processes. While our knowledge of ecological phenomena is vast relative to many landscape metrics, thermal patterns which shape landscape mosaics are largely unknown. To address this disconnect, we investigated the thermal landscape by measuring black bulb temperature (T_bb) at intervals as small as 15 min across 3 yr in a mixed‐grass shrub vegetation community. We found that the thermal landscape was highly heterogeneous displaying a prevalence for thermal extremes (i.e., T_bb > 50°C) and that T_bb was driven by the synergism of environmental, terrain, and vegetation factors. Specifically, variation of T_bb on the landscape was best predicted by the inclusion of ambient temperature (T_air), solar radiation (S_rad), low woody cover, and tall woody cover as variables. Moreover, models of single vegetation parameters (i.e., bare ground, low woody, or tall woody cover) each had greater relative importance than those containing a single terrain variable (i.e., slope or aspect) based on AIC, providing evidence that vegetation is a key driver of T_bb on the landscape. Within the thermally heterogeneous landscape, tall woody cover moderated T_bb by 10°C more than bare ground, herbaceous, or low woody cover during peak diurnal heating (14:00), and was the only cover type that remained <50°C on average. Given that tall woody cover comprises only about 7% of the landscape in our study, these findings have direct conservation implications for species inhabiting shrub communities, specifically that the distribution of tall woody cover is a spatially limited but key predictor of potential thermal refugia on the landscape. Our findings also demonstrate that local interactions between vegetation and temperature can create thermal patterns that shape dynamic landscape mosaics across space and time. Furthermore, we show that structural heterogeneity can maximize thermal complexity across landscapes which can provide greater potential thermal options for organisms. However, our modeled climate projections suggest that far greater thermal extremes will be possible across increasingly larger swaths of the landscape in the future, making assessments and quantifications of thermal landscapes increasingly critical.

Introduction

Abiotic and biotic components of landscape structure (e.g., vegetation and terrain features) shape landscape pattern and, in turn, drive ecological processes (Turner et al. ). For example, the juxtaposition of thermal patterns across space and time influences ecosystem function by dictating decomposition, nutrient cycling, evaporation, and heat flux (Rosenberg et al. , Hobbie ). Fine‐scale thermal environments also regulate the life stages, daily activity, and thermoregulation of organisms (Sartorius et al. , Angilletta , Sears et al. ). Consequently, both regional distributions of plants and animals, as well as fine‐scale constraints on the movement and resource use of individual organisms, are simultaneously dictated by thermal patterns (Peterman and Semlitsch , Lawson et al. ). Despite influencing a variety of ecological phenomena, as well as being subject to scale dependency and high variability across space and time, temperature remains an understudied component of landscape pattern in ecology (Saunders et al. , Faye et al. ).

The distributions of possible temperatures across landscapes are considerably broader than the thermal tolerances of most organisms (Gilchrist ). As a result, the spatial variation of microclimate (i.e., thermal heterogeneity) dictates how and when organisms utilize landscape patches based on their thermal preferences and limits (Campbell and Norman , Angilletta ). Microclimate is formed by the collective effects of physical factors such as ambient temperature, solar radiation, and wind that determine the thermal conditions near ground level (Porter and Gates , Rosenberg et al. ). Importantly, the survival of some individuals is often contingent upon the immediate accessibility of thermally moderated microclimates (i.e., refuges) during extreme heat events (Suggitt et al. ). Further, species persistence requires regional landscape features that shelter populations from climatic extremes (i.e., refugia; Ashcroft ). Identifying both local thermal refuges and broadscale thermal refugia requires an understanding of the capacity of landscapes to buffer against thermal extremes, and accordingly has become a topic of significant importance recently in ecological studies (Keppel and Wardell‐Johnson , Moritz and Agudo , Hovick et al. ). For many ecosystems, basic questions still remain regarding how much of the landscape may serve as thermal refugia (Keppel et al. ), as well as thresholds at which previously thermally suitable patches become unsuitable at different timescales (i.e., hourly, daily; Tracy and Christian ).

Ecosystem function and biological diversity are widely accepted as being affected by heterogeneity (Christensen , Wiens , Fuhlendorf et al. ). Patchiness of vegetation structure remains the most commonly investigated facet of heterogeneity; however, other fundamental aspects such as microclimate have received less attention (Limb et al. ). In addition to terrain features (i.e., slope and aspect) which influence thermal mosaics (Bennie et al. ), the spatiotemporal variation of microclimate across landscapes is strongly influenced by the interaction between temperature and vegetation structure and composition (Saunders et al. , Jenerette et al. , Schut et al. ). Consequently, the resulting configuration of thermal environments dictates the options available to organisms, especially during extreme temperatures (Guthery et al. , Carroll et al. , b). For example, discrete woody plants can substantially modulate fine‐scale thermal environments through blocking solar radiation (Vetaas ) and provide critical thermal refuge for taxa, ranging from ectothermic reptiles (Attum and Eason , Attum et al. ) to endothermic birds (Guthery , McKechnie et al. , Carroll et al. , b) and mammals (Chappell and Bartholomew , van Beest et al. ). However, the common practice of examining thermal landscapes at ≥1‐km scales can overlook fine‐scale variation in thermal and ecological patterns resulting from site‐specific factors (Suggitt et al. , Logan et al. , Potter et al. ). Overcoming this discrepancy is especially important given that determining the scales at which organisms make selection choices is a key question in landscape ecology (Wiens , Graf et al. , Jackson and Fahrig ).

Understanding the thermal heterogeneity present across landscapes is critical now and will become more critical in the future due to climate change. Increases in both the frequency and intensity of extreme temperature events due to climate change (IPCC ) are predicted to alter the pattern juxtaposition of thermal environments across landscapes (Opdam and Wascher ). Already, temperature increases and thermal extremes have been implicated in the reduced fitness and survival of a wide range of species (Sinervo et al. , du Plessis et al. , Moses et al. ), and these negative impacts are predicted to worsen (Huey et al. , Sinervo et al. , du Plessis et al. , Moses et al. , Cunningham et al. ). For example, broadscale increases in thermal extremes have already been shown to translate to fine scales, causing cascading effects within lizard species whereby constraints on activity (e.g., growth, reproduction, feeding) have led to actual species extinctions (Sinervo et al. ). Nonetheless, information on how broadscale shifts in regional climates manifest at local scales is scarce (Parmesan , Varner and Dearing ), but will be integral for assisting ecologists in disentangling the complex ways that these changes will occur (Helmuth et al. , Angilletta ). Moreover, an understanding of the location and timing of climate change impacts on ecological communities remains largely unclear (Helmuth et al. ), especially regarding the role of landscapes as moderators of thermal extremes (Sears et al. , Potter et al. ).

The purpose of this study was to investigate the relationship between local variation in thermal environments and landscape structure (abiotic and biotic). In addition, we aimed to assess the thermal landscape mosaic of a mixed‐grass shrub ecosystem in order to determine how much of the landscape acts as thermally moderated refugia during summer heat and how the spatial extent of potential refugia constricts with increasing temperatures over short timescales (i.e., within‐day variation). To accomplish this, we assessed site‐specific factors influencing spatial and temporal variability in landscape configuration (i.e., temperature and vegetation structure) under a range of thermal conditions, and upscaled those measurements to elucidate landscape patterns (Liang and Schwartz , Guillevic et al. ). Finally, we modeled how thermal patterns may be altered under future climate change projections to provide a linkage between complex thermal landscapes and anticipated changing environmental conditions.

Methods

Study area

Our study site is located on the 7956‐ha Packsaddle Wildlife Management Area (WMA) managed by the Oklahoma Department of Wildlife Conservation in western Oklahoma. Precipitation ranges from 241.81 to 746.6 mm, with an average of 554.4 mm/yr from 1994 to 2013 for this location (Arnett Oklahoma Mesonet Site; Oklahoma Mesonet , b). Prone to high‐heat events, summer temperatures in the region can reach an average of 25 d/yr >37.8°C (Arndt ).

The vegetation in this landscape consists of shrubs and mixed‐grasses and is dominated by sand shinnery oak (Quercus havardii). In addition to sand shinnery oak, other shrubs include sand sagebrush (Artemisia filifolia), sand plum (Prunus angustifolia), and aromatic sumac (Rhus aromatica; DeMaso et al. , Vermeire and Wester ). Herbaceous plants include little bluestem (Schizachyrium scoparium), big bluestem (Andropogon gerardii), sideoats grama (Bouteloua curtipendula), blue grama (Bouteloua gracilis), western ragweed (Ambrosia psilostachya), Texas croton (Croton texensis), and prairie sunflower (Helianthus petiolaris; DeMaso et al. , Peterson and Boyd ). Tall woody cover makes up only about 7% of the landscape and primarily consists of sand shinnery oak mottes (≥2 m in height), but also hackberry (Celtis occidentalis), soap berry (Sapindus drummondii), black locust (Robinia pseudoacacia), and cottonwood (Populus deltoides). The terrain on the study area is generally flat to rolling.

Data collection

Data were collected during the summers of 2012–2014. To quantify patterns in thermal environments (i.e., thermal heterogeneity) and to assess the distribution of temperatures under thermal extremes, we measured black bulb temperature (T_bb). By combining ambient temperature, solar radiation, and wind effects into a single numeric index, T_bb provides a proxy for environmental conditions and can be used to evaluate variation in microclimates across the landscape (Porter and Gates , Campbell and Norman ). It is closer to the thermal conditions that an organism experiences than ambient temperature (Helmuth et al. ). We used steel spheres painted flat black (hereafter, black bulbs) and placed them at ground level (101.6 mm diameter; 20 gauge thickness) to measure T_bb (Guthery et al. , Allred et al. , Hovick et al. ). Each black bulb was fitted with a temperature probe suspended in the center and connected to a HOBO U12 data logger (Onset Corporation, Bourn, Massachusetts, USA) which recorded T_bb.

At three on‐site meteorological stations situated 2 m above ground level, average hourly ambient temperature (T_air) and solar radiation (S_air) were measured to compare broad environmental conditions with site‐specific T_bb. Meteorological stations were situated along an adjacent east–west orientation <7 km apart and at similar slope (1.54–3.26°), elevation (658–691 m), and vegetation type (low herbaceous cover). T_bb was paired with T_air and S_rad measured from the nearest meteorological station for analysis. T_bb was sampled under a similar range of T_air at bare ground (19.48–40.50°C), herbaceous (16.70–41.88°C), low woody (19.66–41.88°C), and tall woody (19.91–40.84°C) cover types (Table ).

Results of multiple regression models of black bulb temperature (T_bb) modeled as a function of ambient temperature (T_air), solar radiation (S_rad), slope, eastness aspect (East), and northness aspect (North), as well as bare ground (Bare), low woody (Low), and tall woody (Tall) cover, and interactions

Candidate model	K	AIC	ΔAIC
Tair + Srad + Low + Tall	5	28,007.03	0.0
Tair + Srad + North + Bare + Low + Tall	7	28,010.49	3.46
Tair + Srad + Slope +East+ Low + Tall	7	28,011.06	4.03
Tair + Srad + Slope + Low + Tall	6	28,013.30	6.27
Tair + Srad + Slope + North + Slope × North + Bare + Low + Tall	9	28,015.76	8.73
Tair + Srad + Slope + East + Slope × East + Low + Tall	8	28,016.32	9.29
Tair + Srad + Slope + East + Slope × East + Bare + Low + Tall	9	28,017.57	10.54
Tair + Srad + Tair × Srad + East + Bare + Low + Tall	8	28,022.86	15.83
Tair + Srad + Tair × Srad + Low + Tall	6	28,023.88	16.85
Tair + Srad + Tair × Srad + Bare + Low + Tall	7	28,025.17	18.14
Tair + Srad + Srad × Low + Srad × Tall	5	28,025.56	18.53
Tair + Srad + Tair × Srad + North + Bare + Low + Tall	8	28,027.34	20.31
Tair + Srad + Tair × Srad + Slope + Bare + Low + Tall	8	28,031.45	24.42
Tair + Srad + Tair × Srad + Slope + North + Slope × North + Bare + Low + Tall	10	28,032.61	25.58

Notes

Models include T_bb measured during peak diurnal heating periods (11:00–15:00 h) at the Packsaddle WMA, Oklahoma, USA, 2012–2014 (n = 4632). Models were ranked by AIC value.

Because biophysical and ecological detail is reduced when studies rely solely on broadscale climate parameters to generalize fine‐scale thermal environments (Sears et al. , Duarte et al. ), we quantified site‐specific T_bb within diurnal periods. To assess spatial variation in thermal patterns across the landscape, we used a stratified random sampling approach to distribute sampling locations proportional to available vegetation cover classes and across the distribution of slope and aspect features within the study area as mapped using ArcGIS 10.3 (Environmental Systems Research Institute [ESRI], Redlands, California, USA). Vegetation cover was assessed using an IKONOS multispectral image with 2 m spatial resolution collected during summer 2013. Specifically, a supervised classification using 125 points with known vegetation cover was used to create and train a thematic layer of cover classes on the study area. The resulting map was further refined through groundtruthing at 215 points prior to and following map creation. Cover classes were categorized as bare ground, herbaceous, low woody, and tall woody (≥2 m). To evenly distribute points across the available slope and aspect features on the study area, we calculated spatial topographic variables using ArcGIS 10.3 (ESRI). Specifically, we used a 1/3‐arc second Nation Elevation Dataset (10 m resolution; https://gdg.sc.egov.usda.gov/) Digital Elevation Map (DEM) to calculate aspect (0–360°) and slope (0–90°). The range of terrain aspect at sampling points (0–337°) was similar to that available on the landscape (0–358°). Terrain slope at sampling points ranged from 0° to 23.30°, and similarly, 99.8% of this flat‐to‐rolling landscape had a slope of less than 24° (87% and 62% were less than 7° and 1.5°, respectively). Therefore, sampling sites were adequately distributed across the landscape to capture variability in both vegetation and topographic features.

To account for structural heterogeneity and to assess fine‐scale thermal variation, we quantified thermal environments using sampling arrays consisting of three black bulbs connected to a single data logger. Arrays were arranged so that one black bulb was positioned at each random point and also at 2 and 4 m from each random point in a randomly determined cardinal direction to assess fine‐scale variation in T_bb. Thermal sampling was conducted in June, July, and August to optimize our chances of capturing a full range of summer conditions, given that temperatures can be highly variable during summer, especially in the Southern Great Plains of North America. To avoid systematically biasing our thermal sampling, we sampled vegetation types proportional to their availability continuously as the study progressed so as not to over‐ or underrepresent certain vegetation types during hotter or cooler periods (Table ). Because our main objective of the study was to rigorously examine differences in the thermal landscape under a host of temperatures, we sampled from the same pool of stratified random points resulting in 318 sampling events over 3 yr. Of those, we sampled 107 in 2012, 127 in 2013, and 84 in 2014 (each sampling location was a cluster of three black bulb thermometers). We measured T_bb at 15‐min intervals for 24 h, and T_bb measurements were averaged across each hourly period (Logan et al. ) to allow for comparisons with hourly averages of T_air recorded at local meteorological stations. This sampling effort provided 91,584 sampling occasions and 22,896 possible hourly averages. Due to rodent damage (i.e., cords being chewed), 21,988 hourly averages were used in our analysis for 24‐h time periods. Rodent damage was not restricted to any one portion of the study area or vegetation type, and although we took great precautions to avoid it (i.e., cables were buried and hidden under the ground surface), its periodic occurrence was mostly unavoidable given the rigorous sampling approach that we employed (i.e., sampling at hundreds of locations across the landscape). Hourly averages for 24‐h periods (n = 21,988), hours experiencing T_air ≥ 25°C (n = 17,491), and diurnal periods of high heat (11:00–15:00 h; n = 4632) were included in the analyses conducted in this study.

To investigate biotic factors influencing site‐ and scale‐specific thermal heterogeneity and buffering, we examined vegetation structure at each sampling location. Cover classes at each sampling point were classified as bare ground (i.e., no vegetation), herbaceous, low woody, and tall woody cover (<2 m). In eight compass directions (cardinal and subcardinal), overhead vegetation obstruction was also measured by aiming a digital carpenter's level affixed to a 2‐m pole at the top of the nearest vegetation and recording the angle (Kopp et al. ).

Analyses

We used generalized linear mixed‐effects models using an AIC model selection framework to assess the influence of environmental variables (T_air, S_rad, and their interaction) during the hottest portions of the day (11:00–15:00 h), terrain features (i.e., slope, aspect, and their interaction), and vegetation cover (i.e., bare ground, herbaceous, low woody, and tall woody cover and their interaction with S_rad) as fixed effects on T_bb as the dependent variable. Sampling location was included as a random effect. Aspect data were circular and thus were arcsine‐transformed prior to analysis resulting in values represented as eastness and northness (Roberts , Nadyeina et al. ). We also assessed the relative importance of each terrain and vegetation variable singularly in models that each contained T_air, and S_rad given that these variables are known to explain a large portion of the variation in T_bb (Hovick et al. ). Models were ranked by AIC value and those within ΔAIC = 2 were considered to have similar explanatory power (Burnham and Anderson ). We estimated model fit by calculating the variance explained by fixed effects (i.e., marginal R²) and the variance explained by fixed and random effects (Nakagawa and Schielzeth 2013: equations 29 and 30) using the MuMIn package in R (Barton , R Development Core Team ). Possible collinearity among variables was assessed with a Pearson correlation test prior to analysis, and correlation was ≤0.40 for all variables (range: −0.008 to 0.40) except for low woody and herbaceous cover (0.78). Consequently, herbaceous cover was not included as a variable in candidate models. Additionally, because T_bb is assumed to be linear function of T_air (Gunderson and Leal , Logan et al. ), we used simple linear models of T_bb as a function of T_air to descriptively depict T_bb occurring among each cover type during comparatively warm periods, when T_air ≥ 25°C, which was the main focus of this study (Hovick et al. ).

To determine whether rates of T_bb modulation varied among cover classes, we calculated the difference between T_bb measurements and the simultaneously recorded T_air value (T_bb − T_air) for each cover class. This standardized difference depicts to what degree cover‐specific T_bb is modulated for a given T_air experienced broadly on the study area. We also investigated the role of vegetation in modulating T_bb using one‐way analyses of variance (ANOVAs; Zar ) with T_bb measurements as a dependent variable among cover class (bare ground, herbaceous, low woody, and tall woody) as an independent variable. In a separate ANOVA, we tested for differences in angle of obstruction (dependent variable) among cover classes as an independent variable. Multiple comparisons were made using a Tukey multiple comparison test (Zar ) following significant ANOVA. Differences were deemed significant at the P < 0.05 level.

Mean (±SE) hourly T_bb recorded during 24‐h daily periods (00:00–23:00 h) within each cover class were calculated to examine temporal variability in T_bb among cover types and to assess thermal heterogeneity throughout diurnal periods. To assess spatial variation in T_bb among random point, 2‐m, 4‐m, and landscape sites, we calculated variances within each site category and compared T_bb measurements using Welch's t tests for each site combination. Along with space and time, we also recorded T_bb across a wide range of T_air conditions and, thus, were able to examine the distribution of T_bb in discrete categories as T_air increased. To assess shifts in T_bb concomitant to changes in T_air, we calculated proportions of T_bb experienced on the landscape across the ranges of ambient temperature (T_air) measurements sampled.

Finally, we predicted what T_bb will occur across the landscape given future climate projections. To account for possible nonuniform rates of increase in T_bb among cover classes, we used the resulting linear equations and slope coefficients from the simple linear models to model projections of T_bb. Specifically, we calculated the T_air that was experienced at 14:00 h (i.e., the hottest period of the day) using the mean T_bb recorded and the specific linear equation for each cover type. We then added the projected ambient temperature increase to the resulting T_air value and solved for the new projected T_bb under end‐of‐century low (2.7°C increase in T_air)‐ and high‐emission (4.6°C increase in T_air) climate change scenarios (www.climatewizard.org; Girvetz et al. ). For our projections of future T_bb, we averaged available climate scenarios for the study area as predicted by both climate change scenarios. Discrete values of current and future T_bb modeled for low‐ and high‐emission climate scenarios were assigned to pixel values in ArcGIS 10.3 for each cover type to display potential changes in the thermal landscape. T_bb was mapped based on natural breaking points in observed T_bb values (<42°C, ≥42°C to <50°C, and ≥50°C) and because these intervals also hold biological relevance given that temperatures >50°C can denature proteins and prohibit biotic processes needed for life (Calder and King , Larcher ).

Results

We found that this shrub‐dominated landscape was composed of a mosaic of thermally heterogeneous T_bb that ranged from 27° to 79°C (Fig. ). When T_air exceeded 39°C, we observed that T_bb was capable of exceeding 70°C, with a large portion of the landscape reaching >50°C (Fig. ). Of the suite of multiple regression candidate models, the model that included T_air, S_rad, low woody cover, and tall woody cover received the greatest support based on AIC ranking (Table ). The effect of each parameter included in the top model was significant (P < 0.05). All other candidate models had ΔAIC > 2 indicating a lack of relative support (Burnham and Anderson ). The global and null models had ΔAIC values of 57.85 and 371.31, respectively. Additionally, low woody and tall woody cover were the only vegetation or terrain variables that were included each of the top models (Table ). The marginal and conditional coefficients of determination (R²) were 0.41 and 0.67, for the top model, respectively. Moreover, marginal and conditional R² values were 0.42 and 0.68 for the global model, respectively, and 0.30 and 0.65 for the null model, respectively. Among the models that included T_air, S_rad, and their interaction along with one single terrain (slope or aspect) or vegetation variable (cover class), the model containing low woody cover received greater support than models for any other variable, followed by tall woody cover (Table ). Single‐variable models containing any one of the vegetation parameters (bare ground, low woody, or tall woody), as well as the null model, each received greater model support than those containing any one terrain variable (slope, sin aspect, or cosine aspect; Table ). However, single‐variable models had very low relative support among the full suite of candidate models.

Distribution of black bulb temperatures (Tbb) observed at ambient temperature Tair ≥ 25°C (i.e., warmer focal periods) at the Packsaddle WMA, Oklahoma, USA, 2012–2014 (n = 17,491).

Models including environmental drivers of T_bb (i.e., T_air, S_rad, and their interaction) and single terrain and vegetation variables to assess their relative importance as predictors of T_bb during peak diurnal heating periods (11:00–15:00 h) at the Packsaddle WMA, Oklahoma, USA, 2012–2014 (n = 4632)

Candidate model	K	AIC	ΔAIC
Tair + Srad + Low Woody	4	28,131.46	0.0
Tair + Srad + Tall Woody	4	28,136.33	4.9
Tair + Srad + Bare	4	28,227.70	96.2
Tair + Srad	3	28,229.87	98.4
Tair + Srad + Eastness	4	28,230.54	99.1
Tair + Srad + Northness	4	28,230.62	99.2
Tair + Srad + Slope	4	28,235.73	104.3

Note

Models were ranked by AIC value.

Although we observed that many factors (i.e., abiotic and biotic) in combination drive T_bb, specific patch characteristic (e.g., cover class) was associated with substantial differences in observed T_bb across the landscape. Linear models demonstrate that T_bb among tall woody cover remained cooler than in other cover types (Fig. ). Additionally, T_bb among cover classes were significantly different (F_3,21864 = 57.9, P < 0.0001). Tall woody cover comprises a small portion of the landscape (~7%) yet provided mean T_bb that was more than 4°C cooler than low woody, bare ground, and herbaceous cover for all measurements (Table ) and more than 10°C cooler during peak diurnal heating (14:00 h; Fig. ). Interestingly, mean T_bb for bare ground was similar to herbaceous cover and low woody (P > 0.05), but the T_bb for herbaceous cover differed from low woody cover (P < 0.05; Table ). For example, we observed maximum T_bb > 70°C in low woody and bare ground cover types and T_bb approaching nearly 80°C in herbaceous cover types, but in woody cover types, T_bb peaked at only 63°C (Table ). Angle of obstruction was 5, 3, and 0.3 times greater than at bare ground, herbaceous, or low woody cover types, respectively (F_3,953 = 391, P < 0.0001). Thus, increasing angle of obstruction generally corresponded to the coolest conditions which were associated with tall woody cover.

Relationship of black bulb temperature (Tbb) as a function of ambient temperature (Tair) ≥25°C (i.e., warmer focal periods) measured at bare ground, herbaceous, low woody, and tall woody cover types at the Packsaddle WMA, Oklahoma, USA, 2012–2014 (n = 17,491).

Range of ambient temperature (T_air) and black bulb temperature (T_bb) sampled from 0:00 to 23:00 h in bare ground, herbaceous, low woody, and tall woody cover types at the Packsaddle WMA, Oklahoma, USA, 2012–2014 (n = 21,988)

Vegetation type	Tair range (°C)	Tbb range (°C)	Tbb, mean (±SE)
Bare	17.21–40.04	17.10–72.55	34.83 (±0.30)a,b
Herbaceous	16.70–41.60	17.32–79.56	35.17 (±0.12)a
Low Woody	16.70–41.60	14.80–75.81	34.15 (±0.14)b
Tall Woody	16.70–40.84	15.44–63.36	30.34 (±0.26)c

⁴Different superscript letters denote significant differences (P > 0.05; Tukey's multiple comparisons).

Daily variation (i.e., 24‐h period) in mean black bulb temperature (Tbb; ±SE) distribution experienced at bare ground, herbaceous, low woody, and tall woody cover types at the Packsaddle WMA, Oklahoma, USA, 2012–2014 (n = 21,988).

Temporally, T_bb did not appear to differ between cover classes from 20:00 to 06:00 h, but exhibited considerable disparity throughout the rest of the day (07:00–19:00 h; Fig. ). We observed substantial heterogeneity during midday and afternoon periods (12:00–14:00 h; Fig. ), and tall woody cover exhibited the most moderated thermal environments across the full range of temporal variability. Specifically, diurnal distributions of T_bb demonstrated that during midday, tall woody cover provided microsites that were on average 10.22°C, 10.74°C, and 12.7°C less than low woody, bare, and herbaceous cover types, respectively (Fig. ). Spatially, we found that variances were similar among all spatial scales that we assessed (range: 137.8–140.8), especially at the finest point scale (137.8) and broadest landscape scale (139.5). Additionally, we found that local‐scale T_bb did not differ between point and 4‐m sites (t = −0.147, df = 11,612.3, P = 0.88), point and 2‐m sites (t = −0.165, df = 11,638.7, P = 0.87), and 2‐m and 4‐m sites (t = −0.312, df = 11,630.4, P = 0.76).

Standardized differences (T_bb − T_air) in the magnitude of thermal environments at tall woody cover types were 0.81°C compared to 5.71°C at low woody, 6.77°C at bare, and 6.95°C at herbaceous cover types (F_3,17433 = 149.4, P < 0.0001; Fig. ). Therefore, tall woody cover moderated T_bb relative to T_air by 7.05, 8.35, and 8.59 times more than low woody, bare, and herbaceous cover types, respectively (Fig. ). The greater thermal buffering capacity (i.e., reduced T_bb) observed at tall woody cover types corresponded to greater angle of obstruction with tall woody cover than other vegetation cover types.

Standardized differences between black bulb temperature (Tbb) and ambient temperature (Tair; i.e., Tbb − Tair; ±SE) depicting the magnitude of relative thermal buffering at bare ground, herbaceous, low woody, and tall woody cover types (n = 17,491) at Tair ≥ 25°C (i.e., warmer focal periods) at the Packsaddle WMA, Oklahoma, USA, 2012–2014. Different letters denote statistical significance (P < 0.05).

The variability of T_bb was evident across a wide range of T_air conditions ranging from 20° to >40°C, and comparatively moderated thermal environments became increasingly scarce as T_air increased (Table ). For example, 86.8% and 98.4% of T_bb measurements exceeded 40°C at T_air ranging from 35° to <40°C and ≥40°C, respectively (Table ). Furthermore, the prevalence of thermal extremes was demonstrated by our finding that 49% and 73.3% of T_bb measurements exceeded 50°C at T_air ranging from 35° to <40°C and ≥40°C, respectively (Table ).

Proportion of black bulb temperature (T_bb) experienced on the landscape across the ranges of ambient temperature (T_air) measurements (°C) sampled during 2012–2014 at the Packsaddle WMA, Oklahoma, USA, 2012–2014 (n = 21,988)

T bb	T air
20° to <25°	25° to <30°	30° to <35°	35° to <40°	≥40°
<20°	14.00	0.69	0.00	0.00	0.00
20° to <25°	77.60	33.02	0.48	0.00	0.00
25° to <30°	7.61	46.80	11.28	0.17	0.00
30° to <35°	0.67	9.11	21.05	2.60	0.00
35° to <40°	0.150	5.87	16.09	10.61	1.60
40° to <45°	0.060	3.02	17.07	16.71	8.82
45° to <50°	0.000	1.13	16.97	20.91	16.31
≥50°	0.000	0.31	17.06	49.00	73.26

Under observed conditions, we found that during peak heating (14:00 h), mean T_bb averaged within bare ground, herbaceous, and low woody cover (93% of the total landscape) exceeded 50°C while T_bb in tall woody cover (7% of the total landscape) was 41.5°C, substantially buffering against conditions prevalent on the surrounding landscape (Fig. ). When modeling future conditions under the low‐emission scenario, tall woody cover shifted to 45.8°C, while the remainder of the landscape remained >50°C, reaching up to 57.5°C for low woody, 58.5°C for bare ground, and 60.8°C for herbaceous cover types. Under the modeled high‐emission scenario, the entire landscape exceeded 50°C, with T_bb for tall woody, low woody, bare ground, and herbaceous cover types being 55.0°C, 61.2°C, 62.9°C, and 65.5°C, respectively (Fig. ).

Thermal characterization of current (A), and future thermal conditions under low (B)‐ and high (C)‐emission scenarios by 2080 at 14:00 h at the Packsaddle WMA, Oklahoma, USA, 2012–2014. Mean Tbb < 42°C (blue), ≥42° to <50°C (yellow), and ≥50°C (red) shown.

Discussion

We found that T_bb across the landscape was driven by many environmental, terrain, and vegetation factors in concert across a mixed‐grass shrub ecosystem. Accordingly, we observed a high degree of thermal heterogeneity, which resulted in a diversity of potential thermal choices for organisms. Our finding that vegetation, specifically low and tall woody cover, best predicted T_bb provided evidence that vegetation is a key driver of T_bb on the landscape that we studied. The thermal landscape was also prone to extreme T_bb during bouts of high ambient temperature and solar radiation, creating a scarcity of thermally moderated microsites. Compared to other cover classes, tall woody cover substantially buffered extreme temperatures further supporting the hypothesis that landscapes can serve as moderators of thermal extremes (Suggitt et al. , Hovick et al. ). These findings also substantiate the premise that the interaction between vegetation and temperature can be informative for making assessments of fine‐scale thermal environments (Saunders et al. ), especially in shrub landscapes. Broadly, our results also demonstrate that thermal landscapes should be viewed as scale‐dependent mosaics that exhibit variability across space and time resulting from the synergistic relationship of both abiotic and biotic factors.

Interestingly, our observation that the thermal landscape exhibited the potential for experiencing T_bb exceeding 70°C was similar to observations from studies in other regions, including the tallgrass prairie in Oklahoma where T_bb exceeded 70°C and Mescalero sand dunes in New Mexico where T_bb exceeded 65°C (Hovick et al. , Sears et al. ; Fig. ). Moreover, during extreme heat (i.e., T_air ≥ 40°C), 98.4% and 73.3% of T_bb measurements exceeded 40°C and 50°C, respectively. Given that temperatures >50°C can result in protein denaturation and the inhibition of biotic processes (Calder and King , Larcher ), thermally buffered sites serve a critical landscape function, especially for species that require regular use of thermal refuges (Lagarde et al. , Shi et al. ). In the face of thermal extremes during peak diurnal heating (14:00 h), tall woody cover moderated T_bb by more than 10°C compared to other vegetation classes (e.g., bare ground, herbaceous, and low woody), demonstrating that a comparatively small segment (7%) of the landscape substantially buffered microclimate. In addition, tall woody cover buffered T_bb seven times more than bare ground, herbaceous, and low woody cover at T_air > 25°C, demonstrating the degree to which vegetation components can dictate thermal environments. These findings have major ecological implications given that in order for organisms to obtain suitable microclimates, such microclimates must be available and accessible within the range of conditions occurring on the landscape (Faye et al. ). If thermally moderated sites are juxtaposed so that they are distant from one another, the potential for ecological traps could be high for species seeking thermal refuge during the heat of the day or during extreme heat events. For example, due to their limited mobility, Moorish tortoises (Testudo graeca) in Morocco require thermal refuges within 0.5 km of their diurnal position during hot summer conditions to avoid heat stress (Moulherat et al. ). Characterizations of thermal landscapes such as those in this study combined with behavioral and physiological data could serve to answer questions about how environmental conditions may induce ecological traps for species now and in the future.

Although tall woody cover that had dense canopies (i.e., high angle of obstructions) was a significant moderator of site‐specific microclimate, it only composed a small proportion of the total landscape (7%). Hence, any practice that reduces or prevents tall woody cover in regions where it is comparatively scarce would likely act to increasingly homogenize the thermal landscape toward more extreme conditions during high heat. However, landscape components that moderate thermal extremes may not always align with the other habitat needs of specific species. This disconnect can occur locally due to predation risks or foraging trade‐offs that could potentially override thermal considerations (Amo et al. , Levy et al. , Cunningham et al. ) or broadly given that structural features on the landscape that moderate extremes may not coincide with the habitat suitability thresholds for sustaining populations. Despite providing thermal buffering and promoting thermal heterogeneity, increases in woody plant abundance pose a major threat to many grassland bird species by altering the structure and composition of critical habitat space (Coppedge et al. , Grant et al. ). As a result, woody plant encroachment in the Great Plains of the United States has led to substantial declines in species diversity and composition for grassland and shrubland birds (Engle et al. ), potentially creating a management conundrum that presents conservationist with a difficult challenge. Therefore, although characterizing thermal landscapes is critical for understanding the environments that organisms inhabit now and in the future, it will also be necessary that efforts to maximize the availability of thermally buffered sites consider the physiological tolerances of imperiled species and incur a minimal amount of degradation to overall habitat quality.

Along with exhibiting spatial variability, as indicated by the similarity of variances among each scale measured (i.e., point, 2 m, 4 m, and landscape), our results also agreed with previous studies that found that temperature is also subject to high temporal variability (Kustas et al. , Raney et al. ). We observed that thermal conditions were spatially homogenous from 20:00 to 06:00 h indicating that organisms would face similar thermal conditions regardless of their position or patch selection on the landscape. However, thermal heterogeneity increasingly diverged among vegetation cover types throughout the day, peaking at 14:00 h providing further evidence that thermal landscapes, and the occurrence of thermal extremes, are temporally dynamic. Our findings of thermal variability at different times of the day were not unexpected; however, the differences (up to >10°C) that we observed between tall woody cover and all other cover types were stark. During this period, organisms likely face a habitat selection conundrum because although the variability among potential thermal choices increases, available thermal environments become increasingly extreme (Carroll et al. 2015b). Therefore, these results also show the significance of acknowledging temporal variation in thermal patterns in assessments of ecological communities, especially given that short bouts (i.e., hours) of extreme heat may be more impactful or detrimental to species than increases in seasonal or yearly averages (Gilbert et al. , Faye et al. ).

The influence of climate change on biota remains uncertain yet holds a high potential for shifting the distributions of species and resulting in other “ecological surprises” (Williams and Jackson ). Therefore, it is critical that we better understand how climate change will alter potential risks to thermal exposure within ecological communities (Parmesan ). Our models of T_bb indicate that microclimates will continue to vary as ambient temperature (T_air) increases due to climate change, but will all be subject to far greater extremes. Under low‐emission climate scenarios, tall woody cover remained the only cover type with T_bb < 50°C (45.8°C), while under high‐emission scenarios, the entire landscape experienced T_bb > 50°C representing highly extreme conditions from a biological standpoint (Calder and King , Larcher ). Even though organisms may actually experience more moderate thermal conditions than those indicated by black bulb temperatures, such scenarios would also likely reduce the capacity of the landscape to provide refuge or refugia for species that are sensitive to thermal extremes. For example, in the absence of tall woody structure, northern bobwhites (Colinus virginianus) inhabiting mixed‐grass shrublands on the subtropical and semiarid periphery of their North American distribution would likely not persist during hot summer conditions due to excessive heat and solar radiation exposure (Guthery ). However, if increases in environmental heat loads reach a point where they overwhelm the ability of microsites to moderate T_bb within thermally tolerable limits for organisms, then the presence of previously viable thermal refuge may be rendered ineffective.

A key goal of landscape ecology is to investigate the structural features and patch configurations that ultimately create mosaic patterns (Wiens and Milne ). This assessment of current and future thermal environments provides a step toward understanding the role of abiotic factors as drivers of ecological complexity and increases our understanding of thermal landscapes now and in the future. Our findings reinforce that thermal landscapes should be viewed as dynamic systems that are subject to diverse spatiotemporal variation (Saunders et al. ) and in mixed‐grass shrub systems can be largely a by‐product of the interaction between temperature and vegetation structure. Specifically, tall woody cover was the primary component for the maximization of thermal buffering across the landscape by providing shade which can be at premium on the landscape during large parts of the day. Thus, in the absence of tall woody cover, thermal patterns would likely shift toward increasingly greater and more spatially homogenized extremes. In mixed‐grass shrub systems, vegetation can provide thermal heterogeneity and thermal buffering (Attum and Eason , Carroll et al. , b, this study), and therefore, conservation practices should focus on preserving both vertical and horizontal vegetation structural heterogeneity to create complex thermal landscapes. Doing so will maximize the capacity of landscapes to support a diversity of species and their wide‐ranging thermal tolerances. However, our models of the changes in microclimatic conditions associated with future climate change suggest that landscapes will be subject to far greater thermal extremes across increasingly larger swaths of the landscape, subsequently increasing the need to understand the drivers of thermal patterns multiple scales.

Acknowledgments

Funding was provided by the Pittman‐Robertson Federal Aid to Wildlife Restoration Act under project W‐161‐R (F11AF00069) of the Oklahoma Department of Wildlife Conservation and Oklahoma State University, administered through the Oklahoma Cooperative Fish and Wildlife Research Unit (Oklahoma Department of Wildlife Conservation, Oklahoma State University, U.S. Geological Survey, U.S. Fish and Wildlife Service, and the Wildlife Management Institute cooperating). We received further support from the Oklahoma Agricultural Experiment Station. Portions of this manuscript result from work supported by the National Science Foundation under Grant No. IIA‐1301789 and Grant No. 1413900. We thank Oklahoma Department of Wildlife Conservation personnel Scott Parry, Jeremiah Zurenda, and Todd Atha for their efforts in managing the Packsaddle WMA. We also thank many hardworking research technicians for their assistance in data collection. The capture and handling procedures used in this study were approved by The Institutional Animal Care and Use Committee at Oklahoma State University (IACUC; Protocol No.: # AG 11‐22).