Widespread shifts in land cover and land management (LCLM) are being incentivized as tools to mitigate climate change, creating an urgent need for prognostic assessments of how LCLM impacts surface energy balance and temperature. Historically, observational studies have tended to focus on how LCLM impacts surface temperature (
Land cover and land management (LCLM) exert a powerful influence on surface and air temperature, at local to planetary scales (Diffenbaugh 2009;Bright et al. 2017;Williams et al. 2021;De Hertog et al. 2022). The influence of canopy structure and function on albedo, latent heat flux (or evapotranspiration), and sensible heat exchange can cause annually-averaged changes in surface temperature of 1-2 K or more between co-located ecosystems that experience the same macro-climate but that differ in LCLM (Juang et al. 2007;Luyssaert et al. 2014;Hemes et al. 2018;Zhang et al. 2020;Duman et al. 2021). At larger scales, plant-driven variability in latent heat, sensible heat, and momentum fluxes can cause changes in cloudiness, precipitation, and circulation patterns that feed back onto surface temperature regimes through non-local effects and "teleconnections" (Swann et al. 2012;Duveiller et al. 2018;Lague et al. 2019;Winckler et al. 2019;Cerasoli et al. 2021;De Hertog et al. 2022). Now, as scientists and societal leaders grapple with the urgent challenge of stopping climate change, the need to understand the "biophysical impacts" of LCLM is especially great (Jackson et al. 2008). "Nature-based climate solutions", which rely on LCLM shifts to enhance ecosystem CO2 sequestration and storage, are rapidly gaining support in public and private spheres (Nolan et al. 2021;Novick et al. 2022a,b;Seddon 2022). However, the holistic climate impact of these strategies remains poorly constrained, with biophysical impacts on temperature representing a major source of uncertainty (Hemes et al. 2021;Williams et al. 2021;Novick et al. 2022b). In some cases, LCLM shifts can directly reduce local temperature (Hemes et al. 2018;Zhang et al. 2020), potentially functioning as a tool for local climate adaptation in addition to global climate mitigation (Seneviratne et al. 2018).
However, in other situations, a LCLM change designed to enhance carbon uptake can increase temperature (Lee et al. 2011;Lombardozzi et al. 2018) -an adverse consequence at local scales and one that could overwhelm the net cooling from increased carbon uptake at planetary scales (Williams et al. 2021). For these reasons, considering biophysical impacts on energy budgets is critical when designing policies that incentivize widespread changes in LCLM (Duveiller et al. 2020;Anderegg 2021, Novick et al. 2022b).
These biophysical impacts originate from direct, local effects on surface temperature (Tsurf), which is arguably the most relevant driver of plant biological function (Luyssaert et al. 2014) and it is a first order driver of downstream effects on air temperature. All else being equal, an increase in albedo will tend to decrease Tsurf (Lee et al. 2011;Bright et al. 2017). Of course, in the natural world, all else is rarely equal. Ecosystems with more leaf area tend to have lower albedo, but they also tend to evaporate more water, which in isolation promotes local surface cooling (Bonan 2008). Likewise, taller ecosystems tend to be more aerodynamically rough, which increases the efficiency of heat transfer from the surface to the air, thereby reducing Tsurf (Burakowski et al. 2018). The interplay of these three mechanisms has been well-studied, at least for some LCLM shifts. We know, for example, that reforestation tends to promote local surface cooling in the tropics and much of the temperate zone (Li et al. 2015;Locatelli et al. 2015;Bright et al. 2017;Zhang et al. 2020) but can lead to local warming in boreal climates where snow-masking albedo effects dominate (Lee et al. 2011;Li et al. 2015).
Likewise, irrigation is widely recognized to have a surface cooling effect driven by increased evapotranspiration (Kueppers et al. 2007;Mueller et al. 2016;Li et al. 2020, but see Gormley-Gallagher et al. 2022).
However, with notable exceptions (Baldocchi & Ma 2013, Lobell et al. 2008;Teuling et al. 2010;Thiery et al. 2017;Seneviratne et al. 2018), most data-driven explorations have focused on longterm Tsurf impacts at the annual timescale. Systematic evaluations of how biophysical impacts differ across seasons, or at night versus during the day, are scarcer. We also lack systematic frameworks for understanding if a long-term biophysical cooling (or warming) effect persists during periods of extreme heat. These are important unknowns, because the consequences of a temperature shift depend very much on when it occurs. For example, a local cooling effect may be most important on hot summer days, when elevated temperatures are especially damaging for ecological and human systems.
Moreover, the drivers of longer-term shifts on Tsurf (e.g., albedo, roughness, and evapotranspiration) can be traced to features of canopy structure that tend to be relatively constant within a season. Thus, LCLM impacts on Tsurf at daily to weekly timescales could reflect an entirely different set of mechanisms, including stomatal responses to rising atmospheric aridity (Novick et al. 2016;Grossiord et al. 2020), or changes in incident radiation that frequently co-occur with rising Ta (Fu et al. 2022).
While Tsurf is a first-order driver of surface energy budgets and a key regulator of plant function, the near-surface air temperature (Ta,local) is an equally relevant target for climate mitigation and adaptation (Winckler et al. 2019;Novick & Katul 2020). This is the temperature experienced by humans. Moreover, since air moves, wind-driven advection is one mechanism by which land cover impacts can be propagated across the landscape (Li & Wang 2019;Barnes et al. 2023, Cohn et al. 2019). Overall, land cover effects on Ta,local have historically been challenging to evaluate. Part of problem relates to data scarcity, as no technique yet exists to directly measure near-surface Ta,local from space, and ground-based measurements are most often constrained to data from meteorological stations, which tend to be situated in grass fields. Land cover impacts on Ta,local are also conceptually more complex as they integrate a wider range of boundary-layer processes including shifts in boundary layer height (Helbig et al. 2021) and entrainment of air from aloft (Fisch et al. 2004).
Finally, whether the target is Tsurf or Ta,local, most investigations of land cover impacts on temperature parse the results on the basis of plant functional type (PFT), for example by comparing evergreen forests to deciduous forests to grasslands to croplands (Lee et al. 2011;Bright et al. 2017;Duveiller et al. 2020;Zhang et al. 2020). While the PFT-based classification is convenient, within a given PFT, impacts on Tsurf or Ta,local can vary substantially (Alkama & Cescatti 2016;Duveiller et al. 2020;Zhang et al. 2020), likely reflecting heterogeneity in climate conditions or canopy structure and function. Approaches that link changes in Tsurf and Ta.local to mappable climate and canopy characteristics would enhance our predictive understanding LCLM impacts on local temperature, within and across PFTs.
The objective of this study is to address these knowledge gaps and limitations, by exploring the coupled dynamics of Tsurf and near-surface Ta across a range of timescales, and across broad climatic and land use gradients. Our primary source of data is long-term AmeriFlux tower records, which provide information necessary quantify land-cover driven changes in Tsurf and Ta,local and attribute them to underlying mechanisms. We will meet the objective by pursuing three research questions: ; and [Q3] To what extent does land cover impact Ta,local and its coupling with Tsurf over long-and short time-scales? Our approach is designed to provide actionable guidance for understanding biophysical impacts of LCLM shifts. Towards that end, observed temperature shifts will be evaluated for how they relate to relevant climate drivers and ecosystem characteristics that can be easily mapped or estimated from remote sensing. These include albedo, evapotranspiration, and canopy height. We will also approach much of the analysis by evaluating the extent to which land cover impacts on near-surface Ta can be predicted by land cover impacts on Tsurf, since an abundance of Tsurf data is available from remote sensing (Farella et al. 2022).
Long-term eddy covariance flux tower records were acquired from the AmeriFlux network (Novick et al. 2018) for 130 sites in the United States that have adopted a CC-By-4.0 data sharing license. Metadata on elevation, International Geosphere-Programme (IGBP) plant functional type classification, and Köeppen climate zone were also acquired from AmeriFlux. The analysis is constrained to these 130 AmeriFlux sites for a few reasons. First, the majority of AmeriFlux sites have adopted an open CC-By-4.0 data sharing license. Next, our approach requires that tower data be corrected by a gridded temperature product that is relatively insensitive to fine-scale impacts of land cover (see Section 2.2.3). NOAA's DAYMET product meets this requirement but is only available for North America. Finally, our approach requires an estimate of canopy height, which frequently must be determined by searching previously published work from individual sites -a labor-intensive step that limits the number of sites that can be included in the analysis. Nonetheless, these 130 sites span a broad range of biomes, including 25 grasslands, 29 croplands or cropland mosaics, 22 deciduous forests, 35 evergreen or mixed forests (which are grouped together since there were only 3 mixed forest sites), 20 shrublands or savannahs (which are grouped together), and 8 wetlands. While "temperate" (n=55) or "warm continental" (n=24) sites are well represented in the dataset, it also includes 19 sites classified as "Dry semi-arid or arid" and 18 sites classified as "boreal," "sub-arctic,"
or "polar" (see Supplementary Information, or SI, Section S1 and Table S1 for details). However, tropical sites are not represented in this dataset.
Flux tower measurements used in this study include: latent heat flux (LE, W/m 2 ), which is an expression of evapotranspiration in units of energy, sensible heat flux (H, W/m 2 ), various proxies for temperature (K) described in more detail in the next section, incident and outgoing short-wave radiation (RSW,in and RSW,out, respectively), incident and outgoing long-wave radiation (RLW.in and RLW,out, respectively) and horizontal wind speed (WS, m/s). Analysis of LE and H was limited to midday periods (between 10:00 -14:00 hours local time) when turbulence is well developed. Thus, traditional approaches for filtering flux tower data for periods of insufficient turbulence (e.g., using friction velocity, Pastorello et al. 2020) were not applied, since they primarily screen night-time data.
Threshold filters were applied to exclude LE or H exceeding +/-1000 W/m 2 . Information on canopy height was obtained from AmeriFlux or extracted from the published literature. In a few cases, it was estimated based on the IGBP classification and site photographs. Finally, some analyses rely on an estimate of the aerodynamic conductance for heat transfer (hereafter Ga), which was estimated by inverting measured sensible heat fluxes as described in Novick & Katul (2020).
2.2 Temperature metrics. We consider multiple temperature metrics that together provide a holistic perspective of LCLM impacts on surface and air temperature (See Summary in Table 1).
The Tsurf,, which is akin to the "skin temperature," represents the projected temperature of all elements of the ecosystem that are unobscured when viewed from aloft. In dense forests, it mostly represents the temperature of sunlit canopy leaves. In sparser ecosystems, it includes both vegetation and background soil. Tsurf is sensitive to variation in radiation and features of canopy structure and function that impact albedo, evapotranspiration, and Ga. While it can be estimated from spaceborne retrievals of thermal infrared imagery (Farella et al. 2022), here it was inferred radiometrically from measurements of outgoing long-wave radiation using the Stefan-Boltzmann law: RLW,out = 𝜎𝜀 𝑠 𝑇 𝑠𝑢𝑟𝑓
, where 𝜎 is the Stefan-Boltzmann constant (5.67 x 10 -8 W/m 2 /K 4 ). The emissivity (𝜀 𝑠 ) was prescribed from an empirical relationship with albedo (𝛼) (𝜀 𝑠 = 0.99 -0.16 𝛼, Juang et al. 2007), where 𝛼 is determined from the measured 𝑅 𝑆𝑊,𝑜𝑢𝑡 /𝑅 𝑆𝑊,𝑖𝑛 at midday.
A consistent definition of near-surface air temperature that can be applied across ecosystems of varying heights has been elusive (Bonan 2008;Winckler et al. 2019, Körner & Hiltbrunner 2018). A good starting point is the air temperature measured ~2m above the ground (the so-called "screen height"). Most long-running meteorological stations measure and report air temperature in this way. Since these stations are usually situated in short grass fields over which the roughness sublayer is very thin, weather station air temperature measurements are made in the surface layer where temperature varies predictably with elevation according to similarity theory (Campbell & Norman 2000).
Outside of meteorological station networks, one of the best sources of data on Ta,local comes from flux tower networks, which provide infrastructure to enable Ta measurements above the canopies of even very tall forests. However, forested flux towers often reside entirely in the roughness sublayer, where similarity theory does not apply, and canopy structure exerts a strong influence on near-surface temperature profiles (Novick & Katul 2020). Consequently, comparing Ta measured on a forest flux tower to Ta measured on a grassland or cropland flux tower may primarily reflect the tower measurement heights, as opposed to biophysical impacts. To minimize this bias, we adopt the approach described by Novick & Katul (2020) to extrapolate air temperature to a reference elevation in the surface layer that scales with canopy height, and where canopy effects should be minimal. The approach carries some uncertainties, including in the estimation of the GA (see discussion in SI Section S2), which were reduced by filtering Ta,extrap whenever it exceeded logical thresholds (e.g. |Ta,extrap|>60 °C) or when the difference between Ta,extrap and the tower-measured air temperature (Ta,tower) exceeded 15 °C. Importantly, all analyses were also repeated using the Ta,tower as the proxy for Ta,local (as opposed to Ta,extrap, see results in the SI). Overall, conclusions do not fundamentally change whether Ta,extrap or Ta,tower is used to represent Ta,local, though LCLM effects on Ta,local are less obvious when using Ta,tower.
To understand how land cover affects local temperature at scales over which macro-climate can vary, it is necessary to have a reference against which Tsurf and Ta,local can be corrected and thus become intercomparable. Here, the "reference" is information on the macro-scale maximum and minimum near-surface Ta (hereafter Ta,macro) provided by the 1 km gridded Daymet product (Thornton et al. 2022). Daymet interpolates and extrapolates weather station data using statistical techniques that minimize the influence of heterogeneous land cover (see Barnes et al. 2023 for more detail). Estimates of daily minimum and maximum Daymet Ta were extracted for all pixels containing the flux towers used in this study, for the entire period of each flux tower record.
The locally measured maximum and minimum Tsurf and Ta,local were corrected by the Daymet Ta,macro products according to:
∆𝑇 surf,min = 𝑇 surf,min -𝑇 a,macro,min Eq. 3 ∆𝑇 a,local,min = 𝑇 a,local,min -𝑇 a,macro,min Eq. 4
These metrics function as direct indicators of the LCLM impacts on local surface and air temperature.
When evaluating the results, the difference in the ∆𝑇surf from one site to the next is more informative than the absolute value of these metrics. For example, prior work suggests (and this study will confirm) that ∆𝑇 surf,max is almost always positive (Mildrexler et al. 2011;Guo et al. 2022), meaning that during mid-day periods, the surface is usually warmer than the air. The question then becomes how much warmer is ∆𝑇 surf,max in one site compared to another.
Exploring how the temperature metrics are coordinated and coupled: Addressing Q1 requires an analysis of LCLM impacts on maximum and minimum Tsurf across seasons and over the diurnal cycle. For this, the ∆𝑇 surf,max and ∆𝑇 surf,min were determined separately for the winter (December, January, and February) and summer (June, July, August). They were then compared to each other for how they relate to ecosystem characteristics including daytime mean LE (measured at midday between 10:00 -14:00), canopy height, and aerodynamic conductance.
Next, to understand how LCLM affects the within-season dynamics of maximum Tsurf (e.g., Q2), the sensitivity of summer Tsurf,max to Ta,macro,max (hereafter dTsurf,max/dTa,macro,max) was calculated in each site as the slope of the relationship between these two variables, using all available hourly data collected from midday periods. The sensitivity of LE, H, and incident solar radiation (Rin, sum of short and longwave radiation) to maximum macro-scale air temperature (e.g., the dLE/dTa,macro,max, dH/dTa,macro,max, and dRin/dTa,macro,max) were also calculated for mid-day summer periods to explore how well they predict dTsurf,max/dTa,macro,max across sites.
The remaining analyses focus on the extent to which surface and air temperature dynamics are coupled over a range of timescales (e.g., Q3). First, the min and maximum ∆𝑇 a,local were compared to min and max ∆𝑇 surf , to understand how well seasonally-averaged LCLM impacts on Ta,local can be predicted given measured impacts on Tsurf. Next, the sensitivity (or slope) of Ta,local,max to Tsurf,max (e.g., dTa,local,max/dTsurf,max) was also calculated for each site at midday. The dTa,local/dTsurf is a particularly useful metric for examining the extent to which local surface and air temperature are coordinated within a site, and how cross-site variability in land cover affects this coupling. Following from the analysis presented in SI Section S2, beginning with the surface energy budget equation, and focusing on periods during which leaf area, height and albedo are relatively stationary, the maximum dTa,local/dTsurf can be expressed as:
where 𝜌 is air density, cp is the specific heat capacity of dry air, 𝜎 is the Stefan-Boltzmann constant, and all other terms have been previously defined. This expression suggests that the dTa,local/dTsurf should be positively related to the sensitivity of LE to surface temperature (Term 1), negatively related to the sensitivity of Rin to surface temperature (Term 2), and influenced by the reference surface temperature (here surrogated to mean summer Tsurf, Term 3). All three terms are mediated by the Ga.
The prediction that the dTa,local/dTsurf should be positively related to dLE/dTsurf is somewhat counterintuitive, given the fact that an increase in LE tends to cause a decrease in local temperature.
As explained in more detail SI Section 3, if an increase in LE occurs despite an increase in Tsurf (e.g., both dLE and dTsurf are positive), then in the absence of changes in Rin, a decrease in sensible heat flux is required to sustain energy balance. In this scenario where Tsurf is increasing, a decrease in H can only be accomplished by an even larger increase in Ta,local.
(e.g., Q1): Land cover impacts on Tsurf,max are profound. As expected, the ∆Tsurf,max is almost always positive (averaging 6.2 K in summer and 5 K in winter). In the context of the goals of this study, the more important result is that the range in ∆Tsurf,max varies by more than an order of magnitude across sites (from <0 to more than 20K), and moreover, the range approaches or exceeds 10 K in most PFTs (Figure 1a,b). In general, the summer ∆Tsurf,max is highest in grasslands and savannahs (Figure 1a, where it frequently approaches or exceeds 15-20 K), and lowest in croplands, deciduous forests, and wetlands. However, the standard deviation of summer ∆𝑇 surf,max ranges from 1.7 -6.1 K within each PFT, which is comparable to the standard deviation across all sites (≈5 K). The variability across sites remains high when the data are parsed by Köeppen climate zone (Figure 1g) or even restricted to a single PFT within a single climate zone (Figure 1h). During the dormant season, the magnitude and variability of ∆Tsurf,max are somewhat reduced (compare Figure 1a to 1b), The winter and summer ∆𝑇 surf,max are reasonably well correlated (R 2 = 0.70), but the slope defining the linear relationship between these two variables is 1.27, driven by the fact that ∆𝑇 surf,max tends to be higher in summer versus winter in many grasslands and savannahs (Figure 1c).
Land cover impacts on minimum temperatures are smaller in magnitude and less variable within and across PFTs (Figure 1d,e), and the winter and summer ∆𝑇 𝑠𝑢𝑟𝑓,𝑚𝑖𝑛 are strongly correlated across sites (R 2 = 0.78, slope = 0.94, Figure 1f). Finally, the cross-site relationship between ∆𝑇 𝑠𝑢𝑟𝑓,𝑚𝑖𝑛 and ∆𝑇 𝑠𝑢𝑟𝑓,𝑚𝑎𝑥 was very weak in both summer and winter (R 2 < 0.1) and the scatterplots for these variables again reveal the tendency for ∆𝑇 𝑠𝑢𝑟𝑓,𝑚𝑎𝑥 to greatly exceed ∆𝑇 𝑠𝑢𝑟𝑓,𝑚𝑖𝑛 (See SI Figure S1). Next, the ∆𝑇 surf,max and ∆𝑇 surf,min were explored for how they vary in relation to LE, canopy height, and albedo. During the summer, ∆𝑇 surf,max decreases as mean mid-day LE increases (Figure 2a), which is expected from first principles related to evaporative cooling. After controlling for LE, sites with higher albedo tend to have higher ∆𝑇 surf,max , particularly when mean summer daytime LE is <250 W/m 2 (Figure 2a). All else being equal, brighter sites should be relatively cool, and thus this result is somewhat surprising. It can be explained by observing that, over this same range of LE, sites with canopy height > 2 m tend to be associated with surfaces that are much cooler than the surfaces of shorter ecosystems (Figure 2b). Canopy height is directly proportional to aerodynamic conductance (R 2 = 0.57, and see SI Figure S2a), enhancing sensible heat flux and surface cooling in taller sites.
Canopy height is also related to albedo (Figure S2b), with taller sites having lower albedo. Nonetheless, the GA impacts appear to dominate; during mid-day summer periods, after controlling for evaporative effects, increases in Ga linked to increased canopy height overwhelm warming effects of lower albedo in taller (and mostly forested) sites. During winter, the relationship between ∆𝑇 surf,max and mean daytime LE is not as clear (Figure 2c,d). However, wintertime ∆𝑇 surf,max still tends to be higher in brighter ecosystems (Figure 2c), but lower in taller ecosystems (Figure 2d), and canopy height is still positively related to winter Ga (R 2 = 0.6, SI Figure S2c). Thus, the same conclusion holds -increases in Ga in taller and more aerodynamically rough stands overwhelm albedo-driven increases in ∆𝑇 surf,max , even in winter.
In contrast, minimum temperatures are not strongly influenced by cross-site differences in LE, canopy height, or albedo (Figure 2e-f). There is, however, a tendency for taller forested ecosystems to have relatively high minimum surface temperatures in both summer and winter, especially when daytime LE is relatively high (Figure 2f,h),
emerging from a non-parametric bootstrap (n=100 iterations).
The dynamics of Tsurf as summer Ta,macro rises (e.g., Q2): In summer, the change in Tsurf as the macro-scale air temperature rises (e.g., dTsurf,max/dTa,macro,max ) is highly variable across sites (Figure 3). In a plurality of sites (~52%), the Tsurf,max rises noticeably more slowly than Ta,macro (dTsurf,max/dTa,macro,max <0.9 K/K). In 27% of sites, dTsurf,max/dTa,macro,max is closer to one (e.g., 0.9 <dTsurf,max/dTa,macro,max <1 K/K). In the remaining sites (~21%), the surface temperature rises faster than the macro-scale air temperature (dTsurf,max/dTa,macro,max > 1.1). The dTsurf,max/dTa,macro,max is usually less than 1 in croplands, deciduous forests, and wetlands (Figure 3d). In grasslands, evergreen forests, and savannahs/shrublands, the dTsurf,max/dTa,macro,max is both higher and more variable across sites.
However, because dTsurf/dTa,macro tends to be <1 more often than not, in most sites, the ∆𝑇 surf,max is reduced as Ta,macro,max shifts from 298 K to 306 K (Figure 3e). The dTsurf,max/dTa,macro,max is only weakly correlated to the seasonally-averaged impacts of land cover on Tsurf across sites (∆𝑇 𝑠𝑢𝑟𝑓,𝑚𝑎𝑥 ), with R 2 = 0.06 and p = 0.01 (Figure 4a). A stronger, negative relationship is observed between the dTsurf,max/dTa,macro,max and the sensitivity of LE to rising air temperature (dLE/dTa,macro,max), with R 2 = 0.18 and p<0.0001 (Figure 4b). In other words, the dTsurf,max/dTa,macro,max is negatively associated with a higher rate of change in LE as macro-scale air temperature rises. When controlling for LE effects, dTsurf,max/dTa,macro,max is larger at sites where Rin increases more rapidly as the macro-climate warms (Figure 4b). Thus, dTsurf,max/dTa,macro,max tends to be higher where the increase in incident radiation as a function of Ta,macro (dRin/dTa,macro,max) is particularly large (compare lines in panel b).
The coupling between surface and air temperature (Q3): Land cover impacts on Ta,local are smaller than the impacts on Tsurf, but still significant. During daytime, the ∆𝑇 a,local,max averages near zero across sites in both summer and winter (Figure 5a,b). However, the variability, which is a key indicator of the influence of land cover effects, is quite large, with ∆𝑇 a,local,max as low as -5 K and more than +5 K in some sites. In contrast, at night, the ∆𝑇 a,local,min tends to be positive in most sites (Figure 5c,d), indicating the local temperature at these sites is elevated above the macro-scale temperature more often than not. This observation holds true when using tower-measured air temperature instead of extrapolated air temperature as the proxy for Ta,local, although the towermeasured temperature is generally lower and less variable than the extrapolated Ta (see SI Figure S3).
Across seasons, the land cover impacts on air temperature are correlated with the land cover impacts on surface temperature (Figure 6), but the magnitude of the impacts on Ta,local tends to be noticeably smaller when compared to LCLM impacts on Tsurf, at least during the daytime (Figure 6a,b). Finally, during summer, the ∆𝑇 a,local,max tended to decrease as both canopy height and mean summer LE increased, which is consistent with the way these variables impact Tsurf, though the correlations were not particularly strong (R 2 <0.2, relationships not shown). Winter ∆𝑇 𝑎,𝑙𝑜𝑐𝑎𝑙,𝑚𝑖𝑛 . Similar results are observed when using tower-measured air temperature as a proxy for Ta,local (see SI, Figure S3). Next, we turn our attention to the dynamic coupling between surface and air temperature in summertime. The sensitivity of local-to macro-scale air temperature (e.g., dTa,local,max/dTa,macro,max) is reasonably well coupled to the sensitivity of surface temperature to macro-scale air temperature (e.g., dTsurf,max/dTa,macro,max, R 2 =0.68, Figure 7). In the majority of croplands and many deciduous forests, the dTa,local,max/dTa,macro,max is greater than dTsurf,max/dTa,macro,max (e.g., points falling above the dotted 1:1 line in Figure 7) which means the local air temperature warms more quickly than the surface as the macro-climate warms. Nonetheless, in the vast majority of sites (including most croplands), the 20 dTa,local/dTa,macro is still less than one, meaning the local air temperature to warms less quickly than the macro-scale air temperature. Another perspective on the within-season coupling of surface and air temperature is gained by evaluating the sensitivity of Ta,local,max to Tsurf,max (e.g., dTa,local,max/dTsurf,max). Across all sites, it is less than 1.0, though it tends to be relatively high in many forests, croplands, and wetlands (Figure 8a). This is a somewhat counterintuitive result, since those sites also tend to be associated with the most pronounced surface and air cooling when averaged over the course of the entire summer (e.g., see
Figure 1). The framework presented in the SI Section S3, and discussed in Section 2.3, is useful for understanding this result. The framework predicts, and the results confirm, that that dTa,local,max/dTsurf,max should be positively related to the sensitivity of LE to increasing surface temperature (e.g., dLE/dTsurf,max, Figure 8c). Moreover, the dLE/dTsurf,max tends to be greater in forests, croplands, and wetlands (Figure 8b) -the same ecosystems with relatively high dTa,local,max/dTsurf,max.
The framework also predicts that incident radiation and the reference Tsurf may be important drivers of dTa,local,max/dTsurf,max. However, in general, these terms were small and uncorrelated with dTa,local,max/dTsurf,maxf (see SI Figure 5b,c), suggesting that evaporation dynamics are the primary driver of the coupling between local surface and air temperature as the macro-climate warms.
Ultimately, however, the influence of land cover on
Ta,local,max (e.g. ∆Ta,local,max) reflects not only the influence of dLE/dTsurf,max, but also the sensitivity of Tsurf,max to rising macro-scale air temperature (e.g. ∆Ta,local,max ∝ 𝑑𝐿𝐸 𝑑𝑇 surf,max • 𝑑𝑇 surf,max 𝑑𝑇 a,macro
). Since the latter term
) is usually much less than 1 (Figure 8a), the ∆Ta,local,max is usually more negative when the macro-scale air temperature is elevated.
The ∆Ta,local,max is enhanced when it is especially hot in only 20% of sites (e.g., sites above the 1:1 line in Figure 8d).
This study was organized around three research questions concerning how land cover impacts on surface and air temperature evolve over sub-annual and sub-daily timescales, and during periods of excessive heat. We first asked: "To what extent, and by what mechanisms, are land cover impacts on Tsurf coordinated across seasons and the diurnal cycle?" We find that land cover impacts on maximum Tsurf are very large in both summer and winter (Figure 1a,b), but they are not well-coupled to impacts on minimum temperature, which tend to be smaller and less variable (Figure 1d, e). The mechanistic underpinnings of this result are not well explained on the basis of plant functional type or climate classification, though a substantial amount of cross-site variability can be traced to variability in midday latent heat flux and canopy height (Figure 2).
We next asked: "To what extent does land cover mediate the within-season dynamics of Tsurf as macro-climate warms?" focusing specifically on mid-day periods. As macroclimate warms, the rate at which Tsurf increased varied substantially across sites (Figure 4a) in ways that were not well explained by seasonally-averaged impacts on Tsurf,max or by plant functional type (Figure 3b). However, in sites where latent heat flux increased as macro-climate warmed, the surface tended to warm more slowly (Figure 4b), such that local surface cooling benefits tend to be enhanced during hot summer days in more mesic sites.
Third, we asked: "To what extent does land cover impact Ta,local and its coupling with Tsurf over long-and short time-scales?" At seasonal timescales, we find that land cover impacts on Ta,local are smaller than, but well coupled with, the impacts on Tsurf (Figure 6). However, within a season as macroscale air temperature rises, the dynamics of Ta,local are not always coupled with the dynamics of Tsurf (Figure 8). Consequently, ecosystems that tend to promote relatively cool Ta,local most of the time may not necessarily sustain this cooling benefit during periods of extreme heat. In the rest of this discussion, these results are explored in more detail, concluding with an assessment of their practical implications.
Tsurf are profound. Land cover impacts on midday Tsurf are stark in both summer and winter, with the summer ∆𝑇 surf,max varying from as low as -2 K to more than 20 K across sites. Relatively few observational studies have explored specifically how land cover affects the maximum summer Tsurf, despite the fact that midday summer periods are arguably the time when the capacity of land cover to mediate rising temperatures is most important. The limited set of studies that do explore biophysical impacts on summertime Tsurf,max tend to confirm the possibility of land cover impacts that range from 5-10 K (Zhang et al. 2020;Guo et al. 2022) or more (Mildrexler et al. 2011).
Beyond representing a key first-order driver of downstream biophysical impacts on air temperature, Tsurf is also a particularly important driver of plant function (Luyssaert et al. 2014) that is well surrogated with the temperature of sunlit canopy leaves in closed canopies. Thus, while it is common to relate cross-site variability in processes like carbon uptake and water use to variability in air temperature (Körner & Hiltbrunner 2018), our results suggest this should be done cautiously, and motivate the use of surface temperature data instead. Moreover, averaging biophysical impacts on Tsurf over the diurnal cycle will tend to obscure the profound impacts of land cover on surface temperatures at midday, when plants are most active and most sensitive to limitations from drought and heat stress.
types" and consider non-radiative processes: In general, the surfaces of deciduous forests tend to be ~2 K cooler than the surfaces of croplands, wetlands, and evergreen forests, which in turn are 5-7 K cooler than the surfaces of grasslands, savannahs, and shrublands (Figure 1a). That forests are relatively cool is consistent with prior work demonstrating that the combination of increased evaporation and sensible heat fluxes tend to suppress forest Tsurf (Bright et al. 2017, Burakowski et al. 2018, Jackson et al. 2008, Juang et al. 2007, Zhang et al. 2020). Croplands tend to have cooler surfaces at midday than grasslands, particularly in summer, which is noteworthy since the biophysical impacts of reforestation and deforestation are often evaluated by comparing forests to just one of these two ecosystem types (Breil et al. 2020, Devaraju et al. 2018, Trail et al. 2013, Zhang et al. 2020). Our results
suggest that cropland to forest and grassland to forest conversions should be evaluated separately.
However, while some general patterns emerge across biome types, land cover impacts on land surface temperature are nonetheless quite large within PFTs (Figure 1a), within Köeppen climate zones (Figure 1g), and even with within specific PFTs in specific climate zones (Figure 1h). Thus, we also explored the extent to which ∆𝑇 𝑠𝑢𝑟𝑓,𝑚𝑎𝑥 can be explained by ecosystem characteristics (including albedo, mid-day LE and canopy height) that are first order drivers of ecosystem energy balance, and which moreover can be mapped across the landscape (Schaaf et al. 2002, Fisher et al. 2020, Simard et al. 2011). In summer, mean daytime LE was the strongest driver of cross-site variability in ∆𝑇 𝑠𝑢𝑟𝑓,𝑚𝑎𝑥 , and sites with high daytime LE had relatively cool surfaces at midday (Figure 2). Albedo effects were not immediately apparent, as they were obscured by increases in aerodynamic conductance in taller (but darker) canopies (Figure 2). Taken together, these results suggest that land cover transitions that increase LE or increase canopy height are more likely to result in local surface cooling.
Our study is limited to ecosystems in the United States, which are diverse with respect to PFT and climate classification (see section 2.1 and the SI), but are overrepresented with respect to temperate sites, and do not include tropical sites or ecosystems managed with agricultural or forestry practices more common in other countries. Thus, future work could focus more strongly on the relative importance of albedo-versus non-radiative forcings on ∆𝑇 𝑠𝑢𝑟𝑓,𝑚𝑎𝑥 across an expanded range of PFTs and climate zones. A systematic effort to collect and redistribute canopy height information from the global set of flux tower sites would enable a more globally representative analysis.
Finally, the relative importance of albedo-versus non-radiative forcings (e.g., LE and H) can shift at coarser scales. Energy consumed by LE and H can be redistributed and re-emitted in the boundary layer (Zhao & Jackson 2014, Novick et al. 2022b), such that the global effects of LCLM change on albedo may be more important than their effects on LE and H (e.g. Bala et al. 2007;Devaraju et al. 2015;Williams et al. 2021). However, it is also important to keep in mind that increased LE can lead to increased cloud formation which increases planetary albedo (Cerasoli et al. 2021) and that changes in H and LE can affect the rate at which radiation is emitted back to space (Zeng et al. 2017).
On hot summer days, the dynamics of evapotranspiration control the dynamics of land cover impacts on Tsurf: In summer, as the macro-scale air temperature rises, the rate at which
Tsurf changes (e.g., dTsurf,max/dTa,macro,max) is very diverse across sites and within PFTs (Figure 3).
Moreover, dTsurf,max/dTa,macro,max is only weakly correlated with the mean growing season ∆𝑇 surf,max (Figure 4a), which implies that the mechanisms that determine the seasonally averaged impacts of land cover on ∆𝑇 surf,max do not necessarily determine Tsurf dynamics within a season. Here, the strongest driver of dTsurf,max/dTa,macro,max across all sites was the concurrent change in LE as the macroclimate warmed (e.g., dLE/dTa,macro,max). When it is high, dTsurf,max/dTa,macro,max is relatively small (Figure 4b). A smaller but noticeable influence of variation in incident solar radiation was also detected, with relatively higher dTsurf,max/dTa,macro,max in sites where dRin/dTa,macro was relatively high.
Ultimately, because dTsurf/dTa,macro tends to be negative more often than not, in most sites, the ∆𝑇 𝑠𝑢𝑟𝑓,𝑚𝑎𝑥 actually declines as Ta,macro shifts from 298 to 306 K. This means that the potential for vegetation to confer a local surface cooling benefit is usually enhanced when it is especially hot, and especially in sites where water supply is sufficiently abundant to permit LE to increase as macro-scale air temperature rises.
Ta,local are smaller than, and coupled with, land cover impacts on Tsurf. During daytime, the magnitude of ∆𝑇 a.local,max is usually much lower than the ∆𝑇 surf,max (Figure 6 a,b), and the land cover impacts on the daytime Ta,local were almost never greater than the land cover impacts on Tsurf (e.g., only one site fell above the 1:1 line in Figure 6a or b). Thus, most of the time, a change in land use or land cover that has a cooling effect on the surface is unlikely to produce a warming effect in the air. Moreover, especially in forests, wetlands and croplands (for which ∆𝑇 a.local,max is closer to the ∆𝑇 surf,max ), the readily measurable ∆𝑇 surf,max could function as a very conservative upper bound for LCLM impacts on Ta,local. Overall, these results are broadly consistent with prior observational work showing that deforestation impacts surface temperature more than air temperature (Alkama & Cescatti 2016; Winckler et al. 2019) and that the slope of the relationship between gridded air and surface temperature products is usually less than one (Mildrexler et al. 2011). 4.5 Within the summer season, land cover impacts on air temperature are not always coupled with impacts on surface temperature. The change in near-surface air temperature as surface temperature rises during summer (dTa,local,max/dTsurf,max) is usually less than one (figure 8a).
However, it tends to be relatively high in croplands, forests and wetlands, which are the ecosystems associated with the most surface-and near-surface cooling. The mathematical framework in SI Section S3 provides a foundation for understanding this result. In particular, it reveals that in sites where changes in LE are coordinated with changes in Tsurf (e.g., dLE /dTsurf,max), as they are in many croplands, forests, and wetlands (Figure 8), the air temperature must rise more rapidly than the surface temperature to achieve a decrease in H that is necessary to close the energy budget. Overall, these results imply that some sites with relatively cool Ta,local during normal summer conditions (e.g.,
Ta,macro = 298 K) may have relatively warm Ta,local when it is especially hot (e.g., sites above the 1:1 line in Figure 7e), and vice versa.
and models. The local effects of LCLM on both surface and air temperature have important ramifications for plant biophysical function and for the potential for LCLM shifts to confer an adaptive air cooling benefit that can extend beyond ecosystem boundaries. However, a complete understanding of biophysical impacts of LCLM shifts must also consider non-local changes in circulation patterns, humidity, and cloudiness that can feed back onto both Tsurf and Ta,local (Swann et al. 2012;Pongratz et al. 2021). Understanding these non-local impacts has largely been the purview of modeling studies (Bonan 2008), which tend to report substantial non-local feedbacks (Lague et al. 2019;De Hertog et al. 2023) that can sometimes counteract local impacts (Hirsch et al. 2017;Williams et al. 2021).
However, results are sometimes contradictory from one model to the next (De Hertog et al. 2023), and models can also struggle to reproduce the local impacts with a precision necessary to be benchmarked against observations (Burakowski et al. 2018). The results presented here cannot directly bridge this gap between local and non-local effects, though they may provide new perspectives against which modeled LCLM change can be evaluated. For example, future work could explore the extent to which models capture the correlation between wintertime and summertime impacts on max and min Tsurf (e.g., Figure 1), or within-season coupling of Ta,local and Tsurf as macro-scale air temperature rises (Figure 7,8). Model-data comparisons like these may be facilitated by the vertical Ta,local extrapolation approach used here, which is thoroughly presented along with the necessary code in Novick & Katul (2020). Since the extrapolation can be applied throughout the surface layer, it can allow temperature impacts to be evaluated across models and observations at a wide range of heights. Additional observations of air temperature profiles made throughout the boundary layer would also benefit future attempts to bridge the models with observations (Helbig et al. 2021).
We conclude by offering the following practical recommendations for evaluating the biophysical impacts of managed land cover transitions on both surface and air temperature:
[1] Biophysical impacts on surface temperature and air temperature are large and should be a key consideration in policy design. The profound impacts of LCLM on surface temperature (e.g., 10 K or more) will have important ramifications for ecophysiological function that ultimately determines carbon uptake. The smaller but still sizeable impacts on air temperature are notable for their potential to be advected across the landscape and impact nearby ecosystems and communities. Finally, the large, within-site differences in surface and air temperature reported here should demotivate the use of air temperature as a factor explaining spatial variability in many ecosystem processes (e.g., photosynthesis), for which surface temperature is likely the more relevant driver.
[2] Albedo changes alone are not sufficient to understand the local impacts of managed land use and land cover transitions. Indeed, across the sites investigated here, the warming effect of reduced albedo was almost completely obscured by non-radiative processes linked to latent and sensible heat flux.
While this result may be surprising given the historic emphasis on ecosystem albedo management in the literature, it is consistent with other works reporting that non-radiative forcings can easily overwhelm albedo impacts at local scales, at least in temperate sites (Bright et al. 2017;Burakowski et al. 2018;Zhang et al. 2020).
[3] Plant functional type is not a sufficient basis on which to evaluate the biophysical impacts of land cover exchange. While our results generally confirm prior results that forests, croplands and wetlands have relatively cool surfaces, we also report extremely high variability within PFTs, and even within PFTs parsed to specific climate zones. Thus, heuristics like "forests cool the surface" are not universally true, and rich opportunities exist to expand on the results presented here that seek to link cross-site biophysical impacts to mappable ecosystem characteristics like average evaporation and canopy height.
[4] Land cover changes that increase evapotranspiration and/or that substantially increase canopy height are likely to confer a surface and air cooling benefit. This explains the tendency for forests, croplands, and wetlands to have relatively cool surfaces and to underlie relatively cool atmospheres, but also explains some of the cross-site variability in land cover impacts on surface and air temperature within these functional groups.
[5] Especially when considering the adaptation potential of LCLM shifts, it is important to focus on time periods when the negative consequences of rising temperature are most severe. LCLM impacts on both Tsurf and Ta,local were far more pronounced during the day than during the night, and were more apparent in summer than winter. Thus, studies that average these impacts over the course of an entire year, or even an entire day, will underestimate the potential for land cover to mediate micro-climate when the consequences for plant and human well-being are most stark (e.g., summer days). Finally, we note that during periods of extreme heat, the dynamics of surface and near-surface air temperature become decoupled. As a consequence, the biophysical impacts of LCLM observed during "normal" environmental conditions may not always predict impacts during more extreme conditions. Some ecosystems that tend to cool the air during normal conditions may have a diminished capacity to do so when it is very hot.