A Century of Forest Regrowth and Snow Loss Alters the Cooling Effect of Historical Land Use in the Upper Midwest

Anthropogenic land use affects climate by altering the energy balance of the Earth’s surface. In temperate regions, cooling from increased albedo is a common result of historical land-use change. However, this albedo cooling effect is dependent mainly on the exposure of snow cover following forest canopy removal and may change over time due to simultaneous changes in both land cover and snow cover. In this paper, we combine modern remote sensing data and historical records, incorporating over 100 years of realized land use and climatic change into an empirical assessment of centennial-scale surface forcings in the Upper Midwestern USA. We show that, although increases in surface albedo cooled through strong negative shortwave forcings, those forcings were reduced over time by a combination of forest regrowth and snow-cover loss. Deforestation cooled strongly (− 5.3 Wm−2) and mainly in winter, while composition shift cooled less strongly (− 3.03 Wm−2) and mainly in summer. Combined, changes in albedo due to deforestation, shifts in species composition, and the return of historical forest cover resulted in − 2.81 Wm−2 of regional radiative cooling, 55% less than full deforestation. Forcings due to changing vegetation were further reduced by 0.32 Wm−2 of warming from a shortened snow-covered season and a thinning of seasonal snowpack. Our findings suggest that accounting for long-term changes in land cover and snow cover reduces the estimated cooling impact of deforestation, with implications for long-term land-use planning.


INTRODUCTION
Anthropogenic land-cover change may exacerbate or ameliorate climate warming through the alteration of energy fluxes others 1988, 2002;Foster and others 2003;Feddema and others 2005;Foley and others 2005;Bonan 2008;McAlpine and others 2010;Mahmood and others 2014). Since 1750, land-cover change has cooled climate by an estimated -0.15 (± 0.1) Wm -2 through increases in reflected shortwave radiation, mostly due to deforestation in areas with seasonal snow cover (Myhre and others 2013;He and others 2014). This radiative cooling effect has led some to argue that reducing forest cover through avoided reforestation or shortened rotation lengths may offset climate warming from greenhouse gasses (Bala and others 2007;Bonan 2008;Anderson and others 2011;Arora and Montenegro 2011;Lutz and Howarth 2015;Schwaab and others 2015;Lutz and others 2016). However, it is not clear that albedo will remain a dominant biophysical force as land use and forest succession continue to develop against the backdrop of a warming world.
The radiative impact of land-cover change depends on both reflected shortwave radiation, which is controlled by albedo, and outgoing longwave radiation, which is controlled mostly by surface temperature. However, both albedo and surface temperature also affect energy balance through local, non-radiative effects. Increases in albedo cool both within the atmospheric boundary layer (locally) and at the top-of-atmosphere (globally) by increasing reflected shortwave radiation, while increases in surface temperature warm within the atmospheric boundary layer by altering turbulent fluxes but cool at the top-of-atmosphere by increasing outgoing longwave radiation and enhancing planetary heat loss (Jin and Dickinson 2010; Peng and others 2014; Ghimire and others 2014). In temperate zones, deforestation-the replacement of forested cover by agricultural or urban land-typically increases albedo in winter, where newly exposed snow increases shortwave reflectance by up to about 80% (Bonan and others 1992;Chapin and others 2005;Bala and others 2007;Davin and Noblet-Doucoudre 2010). In contrast, temperate deforestation typically increases surface temperature in summer, when decreases in leaf area and roughness reduce evaporative cooling and inhibit heat transfer to the atmosphere (Wickham and others 2013;Zhang and others 2014;Rigden and Li 2017). Locally, warming from changes in turbulent fluxes offsets cooling from albedo and dominates local energy balance (Bright and others 2017;Burakowski and others 2018). At the top-of-atmosphere, however, increases in surface temperature increase outgoing longwave radiation from the surface, augmenting radiative cooling from albedo. (Oleson and others 2004;Bonan 2008). Thus, shifts in albedo or surface temperature from vegetation change both affect energy balance, but differ in the seasonality and radiative properties of their effects.
The consequences of past vegetation and climatic change are still evolving today. Much of the temperate northern hemisphere, including Russia, China, Northern Europe, and the USA, was subject to near-universal deforestation followed by either continued agricultural land use or secondary reforestation by deciduous trees (Pfaff 2000;Rudel and others 2005;Lambin and Meyfroidt 2011;Meyfroidt and Lambin 2011;PengSong and others 2018). Simultaneously, climate warming over the past 100 years decreased the length of the snowcover season by 2 weeks and its extent by 5% (Brown 2002. Retrieved from National Snow and Ice Data Center: http://nsidc.org/data/G02131; Sultaire and others 2016). These long-term changes in snow and land cover have taken decades to develop, but have important impacts on climate. As snow cover declines over time, for example, winter shortwave forcing decreases, increasing the relative importance of warmer summer temperatures (Oleson and others 2004;Bonan 2008; Swann and others 2010; for example, Houspanossian and others 2013; others 2009, 2016). Similarly, where initially deforested areas regrow, shifts in species composition may preserve winter shortwave cooling while reducing summer warming (Klingaman and others 2008;Burakowski and others 2016). Together, changes in snow and land cover over decades produce novel patterns of climate forcing that are poorly captured through short-term observations alone (Naudts and others 2016). Long-term data are therefore needed to understand the short-and long-term consequences of forest management decisions for a changing and uncertain future (Betts 2007;Chapin and others 2008;McAlpine and others 2010;Loveland and Mahmood 2014;Mahmood and others 2010).
As a region where typical land-use and climatic patterns are captured in the historical record, the Upper Midwest is a useful case study for long-term change in land cover and snow cover. Much like other parts of the temperate Northern Hemisphere, the Upper Midwestern landscape is a product of land-use legacies that stem from nineteenth-century deforestation, when low-albedo evergreen forests were replaced by higher-albedo crops and deciduous forests (Schulte and others 2007;Rhemtulla and others 2009;Nowacki and Abrams 2015). Over the following approximately 150 years, climate warming decreased snow cover, reducing the albedo of the land surface. Uniquely, these changes were captured through historical vegetation surveys and weather records in the Upper Midwest, making detailed analysis of changes in land cover and snow possible (Peterson and Vose 1997;Kassulke 2009). However, the net impact of land and snow-cover change remains uncertain, as past work has been limited relatively coarse 'potential natural vegetation' type estimates of land cover (for example , Bonan 1999;Matthews and others 2003;Betts and others 2007, but see Steyaert andKnox 2008 andothers 2016 for examples utilizing demographic data), focused mainly on shortwave radiative forcing (for example, Schwaab and others 2015;Lutz and Howarth 2015;Lutz and others 2016) or failed to account for long-term changes in snow cover.
Using newly gridded high-resolution historical survey data and a novel empirical approach, we attempt to resolve these gaps. We focus mainly on top-of-atmosphere radiative effects, which are controlled by properties directly measured by satellites and apply at a global scale. We hypothesize that land-cover change has produced an overall radiative cooling effect in Upper Midwest, but that reduced snow cover and natural forest regrowth have reversed some of the cooling of initial deforestation. To test this idea, we combine modern remote sensing data and historical records to create an empirical analysis of land-cover change that explicitly addresses long-term trends in land use and snow cover.

Region of Interest
The Upper Midwest includes land between 40°and 50°north latitude and 83°and 97°west longitude within the state borders of Michigan, Wisconsin, and Minnesota. This region covers 645,285 km 2 and encompasses 31% of the US Census Bureau designated Midwestern region (US Census Bureau 2010). The climate is continental, with cooler temperatures in the northern part of the region and warmer temperatures to the south. Precipitation is greater in the eastern part of the region than in the west. Summers are warm to hot, and winters are below freezing, with variable snow cover that generally increases with latitude. The growing season ranges from 130 to 170 days in the north to 200-220 days in the south (IIASA/FAO 2012, retrieved from http://gaez.fao.org/Main.html). Regional terrain is glaciated and flat to rolling, with less than 550 m difference between its highest and lowest points.
Modern vegetation communities in the Upper Midwest have been strongly influenced by human land use (Schulte and others 2007;Nowacki and Abrams 2008;Hanberry 2013). Historically, this region encompassed a structured 'tension zone' in which southern savannas and deciduous forests graded into northern mixed and evergreen forests as conditions became wetter and cooler (Curtis 1959). Today, the northern part of the region consists of homogenized post-logging mixed and deciduous forests, while the southern part of the region consists of post-logging deciduous forests and open agricultural land. The persistence of deciduous forests in the non-agricultural portion of the landscape was likely due to successional dynamics, continued logging, altered fire regimes, and/or post-logging ecosystem feedbacks (Frelich 2002).

Methods Summary
We evaluated the long-term impact of land-cover change on surface top-of-atmosphere (TOA) radiative forcing by reconstructing albedo and surface temperature for the pre-Euro-American settlement landscape. Radiative forcing was used to quantify changes in planetary energy balance, with negative forcings indicating a loss of energy from the climate system (that is, climate cooling) and positive forcings indicating a gain of energy in the climate system (that is, climate warming). We projected MODIS-derived, land-cover-specific seasonal cycles of albedo and surface temperatures to the pre-settlement landscape using historical vegetation cover data and empirical models. These empirical models quantified the average decadal seasonal cycle of surface temperature or albedo based on land-cover type, leaf area, snow, and geographic information. We also quantified TOA albedo and surface-temperature forcings associated with changes in snow-cover seasonality using historical and modern snow-cover datasets within the region of interest. Comparing modern and historical land-cover maps of albedo and surface temperature allowed for quantification of radiative forcing from land-and snow-cover changes over the past century. Details on the datasets and methods are described below.

Land-Cover Datasets
Modern land cover (Figure 1e) was characterized with MODerate resolution Imaging Spectrometer (MODIS) data (Friedl and others 2010) using the International Geosphere-Biosphere Program (IGBP) classifications included in the 1000 m MCD12Q1 product (Loveland and others 2000). We limited analysis to the six land-cover classes most common in the region of interest: evergreen needleleaf forest, deciduous broadleaf forest, mixed forest, vegetation mosaic, cropland ('open land'), and urban cover. The first five of these six categories were present in the historical reconstruction, with urban land cover present on the modern landscape only.
Historical vegetation (Figure 1a) was characterized using a Public Land Survey (PLS)-based gridded data product created by Paciorek and others 2016 and Goring and others 2016. The PLS was carried out prior to the first large waves of Euro-American settler colonization ($ 1830-1870) and provides a record of 'pre-settlement' species composition and stem density (Bourdo 1956;Cole and others 1998;Kassulke 2009). The 8-km statistically estimated product used here takes advantage of the observational, fine-scale nature of PLS data while minimizing the known biases associated with past survey practices (Schulte and Mladenoff 2001;Liu and others 2011;Paciorek and others 2016;Goring and others 2016).
Historical vegetation was categorized into landcover classes using a modified IGBP DISCover classification scheme (Loveland and others 2000) to match the modern MODIS land-cover classes. The IGBP scheme was based on percent canopy cover and foliage seasonality, which can be approximated by the stem density and species designations present in the statistically reconstructed PLS data. We used Goring and others' density designations for prairie (< 0.5 stems ha -1 ), savanna (0.5-47 stems ha -1 ), and forest (> 47 stems ha -1 ) while using IGBP criteria for composition dominance (60% cover). Areas with more than 47 stems ha -1 and greater than 60% deciduous stems were classified as 'deciduous broadleaf forest.' Similarly, areas with more than 47 stems ha -1 and greater than 60% evergreen stems were classified as 'evergreen needleleaf forest. Areas with more than 47 stems ha -1 , but less than 60% presence of evergreen or deciduous taxa were classified as 'mixed forest.' Areas with 0.5-47 stems ha -1 were classified as 'vegetation mosaic,' and areas with fewer than 0.5 stems ha -1 were classified as 'open land.'

MODIS Surface Datasets
Land cover was linked to land-cover-specific seasonal cycles of albedo and surface temperature using decadal averages (2002-2011) of MODIS version 5 land surface data (Table 1, Figure 1b   Decadal averages were used to smooth the impact of interannual variability on seasonal cycles. We used four MODIS-derived products: albedo, snow presence, leaf area index (LAI), and land surface temperature. All products were 8-day composites (46 'MODIS dates' per year) retrieved from the Land Process Distributed Active Archive Center (LP-DAAC) and processed using the MODIS Reproject Tool (Dwyer and Schmidt 2006). Quality control and gap filling were carried out prior to analysis for all products except snow probability (see supplemental material). First, we removed MODIS fill value pixels from the dataset. Fill value pixels are pixels for which no acceptable values could be retrieved by the MODIS satellite and algorithms, and are typically indicated by numbers outside the normal data range for the product, for example, 9999. To fill the resulting gaps, we used a modified Savitsky-Golay filter. The modified Savitsky-Golay filter first applied a linear interpolation to missing data and then smoothed the gap-filled signal using convolution.

Albedo
Albedo (Figure 1b) was the first of two surface properties of interest in this paper due to its control over shortwave forcing. Albedo for each land-cover class was quantified with the MCD43A3 Albedo product (Schaaf and others 2002). This product included both black sky albedo (direct light) and white sky albedo (diffuse light). Blue sky (actual) albedo was calculated as the weighted sum of black sky albedo and white sky albedo after accounting for atmospheric direct and diffuse light ratios. Consequently, black sky albedo (BSA) was multiplied by a standard 0.8 (that is, 80% direct light), whereas white sky albedo (WSA) was multiplied by a standard 0.2 (that is, 20% diffuse light) (Lewis and Barnsley 1994): We used blue sky albedo in our analyses.

Snow Probability
Snow probability (Figure 1c) was one of two variables used to predict land-cover-specific albedo. Snow probability influenced albedo through its high reflectivity under sparse leaf cover (LAI < 1.0). We quantified snow probability using the snow presence layer of the MCD43B2 Albedo Quality product. The snow presence layer was a binary product with a one indicating snow-covered albedo and a zero indicating snow-free albedo. The decadal average of this product for each composite date resulted in a prob-ability of snow cover for each pixel and date ranging from 0 to 1. We chose a measure of snow probability because MODIS products do not provide snow quantity (for example, depth or snow water equivalent), and because surface albedo increases only marginally once a continuous snowpack has been established.

Leaf Area Index
Leaf area index (Figure 1d), defined as leaf area per unit ground area (m 2 m -2 ), was the second of two variables used to predict land-cover-specific albedo. Surface albedo was more representative of the canopy at high LAI (LAI > 1), and more representative of the sub-canopy, snow, or soil at low LAI (LAI < 1). LAI was quantified with the MCD15A2 LAI product (Myneni and others 2002). Although this product has known negative biases for evergreen pixels in low leaf area canopies and in the presence of snow (Tian and others 2004; Serbin and others 2013), fully evergreen pixels were uncommon on both the historical (13%) and modern (< 2%) landscapes, and most evergreen pixels had an LAI greater than 1. Additionally, the effects of any biases in LAI applied equally to historical and modern reconstructions, reducing bias in calculated albedo change.

Surface Temperature
Surface temperature (Figure 1f) was the second of two surface properties of interest in this paper due to its control over longwave forcing. Surface temperature for each land-cover class was quantified using the MOD11A3 Surface-Temperature product (Wan and Dozier 1996;Wan 2008). Although cloud contamination and emissivity uncertainties in mixed cover (for example, vegetation mosaic) are known sources of error for this product, it is especially robust over the cropland and closed forests that make up the majority of the region of interest (Strahler and others 1994;Wan 2008). Average surface temperature for each composite date was calculated using a weighted average of day-and nighttime surface temperature, with weights determined by the number of daylight hours for each date.

Model Construction
We developed land-cover specific seasonal cycles of albedo and surface temperature using empirically derived locally weighted scatterplot smoothing (LOWESS) regressions for each MODIS acquisition date (6 vegetation types * 46 dates = 276 regres-sions) (Figure 1h-i, Figure S1, Table S1-S2). LOWESS regressions consisted of two input variables and a response (that is, albedo or land surface temperature). The LOWESS procedure fit a response surface based on these three-dimensional data points, utilizing a user-defined proportion of the data (span) around each point to predict a local value. In contrast to more traditional regressions, LOWESS allows the model to respond flexibly to multi-dimensional changes in the input data while still granting the user some control over model behavior (see supplemental materials). For albedo, we used MODIS LAI and snow probability as predictors, and MODIS blue sky albedo as a response (Figure 1h). For surface temperature, we used latitude and longitude as predictors, and MODIS surface temperature as response (Figure 1i). Longitude and latitude were chosen as inputs for surface-temperature models because they captured broad climatic patterns that influence surface temperature.

Modeled Projections of Albedo and Surface-Temperature
Response to Land-Cover Change LOWESS models projected albedo and surface temperature onto the historical and modern landscape ( Figure 1j). Here, we focused solely on landcover changes, holding snow cover and spatial air temperature variability constant between the historic and modern time periods. Forcings associated with snow-cover change were calculated separately and described in section 'Snow-Cover Shift.' For historical projections, we used the empirical LOWESS models with PLS historical land-and modern snow-cover to derive historical seasonal cycles of albedo and surface temperature. For modern projections, we used the empirical LOW-ESS models with modern land-and snow-cover to derive modern seasonal cycles of albedo and surface temperature. Modern projections from the LOWESS models were compared to upscaled MODIS data for model validation. Validation revealed high agreement across the region for all vegetation types, indicating the robustness of the LOWESS models in projecting seasonal cycles of albedo and surface temperature to land-cover type.

Radiative Forcing of Land-Cover Change
Albedo Forcing Change in albedo was calculated by subtracting historically reconstructed albedo values from modern values on a pixel-by-pixel basis ( Figure 1m): Radiative forcing for albedo was calculated using radiative kernels. Radiative kernels are created by repeatedly perturbing a surface property in a climate model and 'observing' the top-of-atmosphere response. They account for both the decoupling between surface and planetary albedo and albedo feedbacks. The radiative kernel approach is commonly used and preferred to physically based forcing calculations because it captures not only the first-order effect of such a perturbation, but also its short-term climatic feedbacks. We calculated monthly radiative forcing by multiplying Da values by European Centre HAMburg model 6 (ECHAM6) pre-industrial TOA radiative kernels for albedo (Block and Mauritsen 2013): where RF aTOA is the top-of-atmosphere albedo forcing for a given month m (Wm -2 ), Da m is the change in albedo (unitless; %) for month m, and R m is the radiative kernel value (Wm -2 ) for month m.
Surface-Temperature Forcing Change in land surface temperature was calculated by subtracting historically reconstructed surface-temperature values from modern values on a pixel-by-pixel basis ( Figure 1n): Radiative forcing due to changes in surface temperature was also calculated using radiative kernels. We calculated monthly radiative forcing by multiplying DT values by ECHAM6 pre-industrial TOA radiative kernels for surface temperature (Block and Mauritsen 2013): where RF TTOA is the top-of-atmosphere surface temperature forcing for a given month m (Wm -2 ), DT m is the change in land surface temperature (K) for month m, and R m is the radiative kernel value (Wm -2 K -1 ) for month m.

Snow-Cover Shift
Snow Data Snow changes and the associated forcings were estimated independently of landcover forcings by holding historical land-cover constant and applying snow changes from the presettlement to modern time period (Figure 1l). We used monthly 0.5 9 0.5-degree resolution snow water equivalent (SWE) from the terrestrial water budget product created by Willmott and Matsuura (2012) to quantify seasonal changes in snow cover since the historical time period (Figure 1g). The SWE product was derived from monthly Global Historical Climatology Network version 2 (GHCN2) temperature and precipitation data using a modified Thornthwaite water budget procedure (Willmott and others 1985). The GHCN2 record covers several thousand global stations with high station density in the Upper Midwest even in the early nineteenth century. We used SWE from 1900 to 1910 (earliest years available) to characterize historical snow cover, and from 2000 to 2010 to characterize modern snow cover. We used the SWE dataset as a proxy for snow cover due to a lack of historical spatially resolved snow-cover products for the pre-settlement period. We justify this approach by modifying the dataset to focus on snowcover seasonality with long-term decadal and large spatial scale averages. Long-term decadal averages minimized the impact of interannual variability in snow cover, corresponding with our treatment of MODIS land-cover products, while large spatial scale averages minimized the small-scale impacts of SWE that arise from forest cover, wind transport, and snow age effects on albedo and surface temperature.
Albedo Shift We used SWE seasonality to translate changes in the snow-free season into albedo seasonality ( Figure S2). For both the modern and historical landscape, SWE seasonality was established for each pixel by dividing monthly SWE values by a seasonal maximum, yielding monthly values ranging from 0 (minimum snow cover) to 1 (maximum snow cover). MODIS albedo data were shifted on a pixel-by-pixel basis for each month in proportion to the change in SWE seasonality (historical-modern) using Eq. 6: where a Shift(m) is the seasonally shifted 'historical' albedo, a Modern(m) is modern albedo, S Historical(m) and S Modern(m) are normalized historical and modern SWE for a given month, ranging from 0 to 1, (max(a Modern )min(a Modern ) is the seasonal range of modern albedo, and subscript (m) indicates a value that varies by month.
The first term in Eq. 6, a Modern(m) , defined the modern albedo to be shifted, while the second and third terms in Eq. 6 defined the snow-related albedo shift. The second term, (S Historical(m) -S Modern(m) ), quantifies the change in snow seasonality, while the third term, (max(a Modern )-min(a Modern )), is the range of albedo in the pixel to be shifted. We multiply (S Historical(m) -S Modern(m) ) by (max(a Modern )min(a Modern )) because we assume that albedo is at its maximum value at 100% of peak snow and at its minimum at 0% of peak snow.
For example, for a pixel with a Modern(m) = 0.5, S Historical(m) = 0.9 (90% of maximum snow), S Modern(m) = 0.6 (60% of maximum snow due to accelerated melting), max(a Modern ) = 0.8 and min(a Modern ) = 0.2, the albedo shift is (0.30) * (0.6) = 0.18 and the historically reconstructed value is 0.5 + 0.18 = 0.71. We apply Eq. 6 to spring months only because changes in SWE correspond most closely to changes in snow cover in that season. Radiative forcings for the effect of snow on albedo were calculated as for land cover, by multiplying the change in albedo by radiative kernels (Figure 1o).

Surface-Temperature Shift
We used the modern sensitivity of LST to SWE to quantify changes in surface temperature attributable to changes in snow cover ( Figure S3). Through our LOWESS projections of LST to the modern landscape, we observed a systematic bias in the model and observation residuals that significantly covaried with SWE. The slope of this covariation (m) was used as the sensitivity of LST to SWE and emerged through snow water equivalents' impact on sensible heat transfer and latent heat (Cohen and Rind 1991). The sensitivity of LST to SWE was used in conjunction with changes in SWE across time (modern-historical) to estimate the effect of changes in snow on LST: where LST Shift (°C) is the shifted 'historical' surface temperature, LST Modern (°C) is modern LST, DSWE (mm) is the change in SWE from the historical to modern landscape, and m is the slope of the LST residual-SWE relationship (DLST DSWE -1 ) (Figure S3). Thus, we applied the modern spatial LST sensitivity to SWE to the historical landscape to calculate LST changes associated with historical snow-cover change. Radiative forcings for the effect of snow on surface temperature were calculated as for land-cover, by multiplying the change surface temperature by radiative kernels (Figure 1o).

Land-Cover and Snow Change
Historical land cover in the Midwest was mostly forested, with deciduous, mixed, or evergreen tree coverage representing 70% of the historical landscape (1830-1870) (Figure 2A). Deciduous forests represented 34% of historical land cover, mixed forests 23%, and evergreen forests 13%. Evergreen and mixed forests dominated in Michigan and Northern Wisconsin, whereas deciduous forests dominated in central Wisconsin and southern Minnesota. A transition between forest and prairie occurred in western Minnesota and Southern Wisconsin, with evergreen and deciduous forests thinning to vegetation mosaic and open land to the south and west.
Modern land cover was less forested and more deciduous than past land cover ( Figure 2B). Deforestation, the persistent replacement of his-torical forests by open land, vegetation mosaic, or urban land cover, or replacement of mosaic with open land, occurred in 45% of the region, with open land more than doubling. Composition shifts, the replacement of historical evergreen or mixed forests by deciduous forests, and the replacement of evergreen forests by mixed forest, occurred in 18% of the region (25% of historically forested land). No detectable change in forest composition occurred in $ 15% of the region (22% of historically forested land). Historically-open land was replaced with forests or urban land cover in less than 1% of the region ( Figure 6). Combined, deforestation and composition shifts were responsible for 93% of land-cover change.
Modern snow cover (in terms of SWE) was consistently less than past snow cover ( Figure 2C, D). Modern SWE was -6.5 mm less than historic SWE on average, with the greatest decreases concentrated in the melt season (Feb-Apr). Decreases in SWE for the northernmost third of the region were often an order of magnitude greater than those in the southernmost third ( Figure 2C, D),

Model-Data Fit
Modeled reconstructions of modern albedo and surface temperature matched MODIS data reasonably well. The albedo reconstruction fell along a 1:1 line, with clusters likely a result of the forced discretization associated with vegetation-specific LOWESS regressions (annual R 2 = 0.86, Figure S4). We observed no systematic trend in the distribution of grid cells that diverged most from MODIS data, suggesting that albedo may simply be atypical for the assigned vegetation class in some locations. Models did not show overall bias in estimating albedo, with mean residuals < ± 0.6% in all seasons except winter (± 2%). The surface-temperature reconstruction matched MODIS data well, with tight clustering around a 1:1 line and mean residuals < ± 0.2°C in all seasons (annual R 2 = 0.98, Figure S4).
Deforestation (67%) produced strong negative radiative forcings in winter but much weaker radiative forcings in summer ( Figure 5A). In contrast, shifts to deciduous taxa (26%) produced its strongest radiative forcings in summer. Regional radiative forcing was driven mainly by deciduousto-crop, deciduous-to-mosaic, and mosaic-to-crop transitions ( Figure 6).

Snow-Cover Change and Radiative Forcing
Historical snow cover lasted 17 days longer than modern snow cover on average. Changes in snow seasonality of a month or more in the northern third of the region drove a -1% decrease in spring albedo (Figure 7). Snow shortwave forcing was positive and strongest in early spring ( Figure 6), with average spring forcings of +1.43 Wm -2 (95% CI -2.14, +5.01), and average annual forcings of +0.36 Wm -2 (95% CI -0.54, +1.25). Changes in snow seasonality and snow depth al- tered LST in both winter and spring, reducing LST by -0.1°C in winter and -0.09°C in spring. LST longwave forcing was small but greatest in late winter and early spring, with average winter forcings of -0.04 Wm -2 (95% CI -0.13, +0.05), average spring forcings of -0.08 Wm 2 (95% CI -0.25, +0.09), and average annual forcings of -0.04 Wm -2 (95% CI -0.16, +0.03).

Net Forcing of Land-Cover Change and Snow
Regional forcing was most strongly driven by -2.81 Wm -2 of negative forcing from land-cover change. However, these forcings were offset by +0.32 Wm -2 of net forcing from snow-cover loss (Figure 8). Region-wide, the interaction of landcover and snow-cover changes resulted in total forcings with a clear seasonality. In winter, negative shortwave forcings from deforestation and composition shifts resulted in a land-cover forcing of -3.29 Wm -2 that was augmented by weakly negative snow forcings. In spring, positive snow forcings offset almost half of decreasingly negative shortwave forcing from deforestation, resulting in a net forcing of -1.98 Wm -2 . In summer, moderate shortwave forcings from deforestation and composition shift dominated, resulting in a net forcing of -2.48 Wm -2 . Finally, in fall, negative shortwave forcings from deforestation and composition shift decreased slightly, resulting in a net forcing of -2.12 Wm -2 .
Annually, snow forcings offset 11% of landcover forcing (Figure 8). Offsetting effects reduced the net forcing of human land use from -2.81 Wm -2 considering land-cover forcings alone to -2.49 Wm -2 with the effects of snow-cover change included.

DISCUSSION
Over much of the world, including the Upper Midwest, the modern landscape is a legacy of landuse choices made decades to centuries ago. The consequences these long-term patterns have long been recognized for the carbon budget, but less so for the energy budget despite the fact that changes in surface energy balance can produce climate forcings as strong as those produced by greenhouse gasses (Claussen and others 2001;Bonan 2008;Pongratz and others 2010;Anderson and others 2011;Avila and others 2012). Land-cover change in historically deforested zones is often assumed to consist mainly of persistent, near total deforestation (for example, Betts and others 2007; Burakowski and others 2018), especially at the regional and global scales most relevant to climate predictions. This assumption is inconsistent with our empirically based assessment of centennial-scale landcover surface forcings. We showed that the climatic cooling effect of total deforestation is reduced by long-term changes in land cover (for example, secondary forest regrowth) and shrinking snow cover (Figure 8), making a case for greater nuance in current thinking about the climatic impact of deforestation in temperate latitudes.
We found that deforestation indeed made up a plurality of land-cover conversions, but that over 40% of initially deforested land has returned to a forested state (Figure 2A, B). Of these regrown forests, half experienced a shift in composition toward more deciduous species. Spatially heterogeneous forest regrowth and composition shift were detectable in our historical dataset, but have been difficult to incorporate into land surface models or coarse-scale empirical estimates (Pacala and Deutschman, 1995;Pappas and others 2015). As forests continue to regrow, composition shifts may become the dominant long-term effect of historical deforestation, increasing the need to detect and quantify these transitions (PengSong and others 2018).  Most estimates of climate forcings by deforestation are primarily based on observations made shortly after clearing, but such short-term focus leads to an overestimation of the radiative impact of deforestation. Land-cover heterogeneity resulted in a regional forcing that was less negative than that of total deforestation by 55% (applying deforestation forcings to all historically forested areas), a finding in agreement with model experiments by Burakowski and others (2016) who also find that reforestation has a warming effect. Composition shift resulted in weak top-of-atmosphere radiative cooling, mainly in summer, while deforestation resulted in strong top-of-atmosphere radiative cooling, mainly in winter ( Figure 5). Other vegetation conversions created additional seasonal and annual patterns of forcings: weak thermal and radiative warming for shifts to evergreen species and strong thermal warming for transitions to urban land cover. Within each conversion, the magnitude of forcing also varied among vegetation transitions (Figure 6). Evergreen-to-crop transitions, for example, created nearly 6 times the forcing of deciduous-to-mosaic or mosaic-to-crop transitions. In aggregate, diverse vegetation changes resulted in calculated forcings that were less than some coarse-scale estimates (for example, Oleson and others 2004 [their Figure 10]), highlighting the importance of including long-term landscape processes in regional radiative forcing estimates (Euskirchen and others 2016;Thom and others 2017).
Snow-cover change further reduced the radiative forcing of initial deforestation. The snow-cover period declined by about 2 weeks and spring snowcover extent declined by -28% across the Upper Midwestern region over the past 100 years (Groisman and others 2004), indicating that increases in snow melt exceeded observed increases in lake-effect snowfall (Burnett and others 2003). Despite widespread acknowledgment of declining snow cover however, few studies have incorporated it into radiative forcing assessments (Pitman and others 2011; but see others 2009, 2016;Li and others 2018). In areas where snow cover has declined, deforestation produced weaker cooling due to a loss of albedo-induced shortwave forcings (Figure 7). These declines will likely continue in the future. Assuming historical sensitivity of snow-cover forcings to temperature (+0.38Wm -2 forcing/1°C warming; +0.32Wm -2 calculated snow forcing/0.85°C global temperature increase since 1850 = + 0.38 Wm -2°C-1 ), continued snow melt in the region will offset more than half of the remaining cooling from land-cover change by 2100. This estimate is likely to be conservative, as it does not account for the potential acceleration of snowmelt due to snow-albedo feedbacks (Qu and Hall 2007). Thus, we show that declining snow cover has already reduced the cooling impact of land-cover change and will reduce it further in the future.
Planetary longwave cooling was negligible compared to shortwave cooling, but the effect of both snow-cover and land-cover changes is also felt through a change in surface temperature itself (Davin and Noblet-Ducoudre 2010). Both deforestation and composition shift resulted in minor (< + 0.5°C) thermal warming on an annual basis ( Figure 5). However, deforestation amplified seasonal temperature extremes by producing strong summer warming and strong winter cooling, whereas composition shift muted seasonal temperature extremes by producing summer cooling and winter warming. Locally, warmer growing season temperatures in deforested areas may exacerbate the detrimental effects of climate warming (for example, temperature-driven drought) despite radiative cooling. Composition shift, in contrast, may mitigate those effects. Thus, estimates of radiative cooling must be accompanied by an evaluation of differences in energy redistribution at a local level.
In deforested areas, warmer surface temperatures during periods of negative radiative forcing imply that these landscapes have a reduced capacity to redistribute heat. When a surface is absorbing less solar energy (negative forcing), cooler surface temperatures would be expected. We observe the opposite for deforestation, which experiences radiative cooling but thermal warming in summer. Changes in surface temperature due to non-radiative effects, typically controlled by evapotranspiration and turbulent heat transfer, have been increasingly explored in recent work (Rocha and Shaver 2011;Zhao and Jackson 2014;Luyssaert and others 2014;Li and others 2015;Bright and others 2017;Burakowski and others 2018). In the absence of observational information on historical turbulent fluxes, it is difficult for us to calculate such metrics for the deep past. However, loss of energy redistribution with deforestation would be consistent with the average energy redistribution factor values described by Bright and others (2017, [their Table S1]) for forested and open land. This reduced capacity to dissipate heat may lead to greater local temperature increases in response to future greenhouse warming in deforested areas.
Shifts to deciduous taxa, in contrast, experienced less warming in the presence of negative forcing, except in early spring and fall. Early spring warming in evergreen-to-deciduous transitions is consistent with findings by Zhao and others (2014) and may be a result of leaf phenology. In early spring and fall, deciduous forests are leafless, reducing surface roughness and the corresponding turbulent redistribution of energy. During the remainder of the year, however, composition shifts cool when radiative forcing is negative. This effect is most pronounced in midsummer, implying that areas experiencing evergreen-to-deciduous shifts have a greater ability to dissipate energy in that season, and may be more robust to higher summer temperatures under future greenhouse warming.
As Earth's climate continues to warm, snowcover declines will diminish the impact of radiative forcings, while increasing summer temperatures will enhance the importance of energy redistribution. Such shifts have implications for future landuse decisions. In the absence of forest preservation initiatives, climate warming may support the northward expansion of agriculture, increasing deforested area (King and others 2018). Because deforestation cools through snow-mediated forcings, snow-cover loss will reduce the cooling impact of deforestation. Alternatively, reductions in timber harvest due to economic forces or forest conservation initiatives may stabilize or increase forest area . If the goal of subsequent forest management is to restore historical forests, regrowth will reverse some of the past century and a half of land-cover change, moving land-cover forcing toward zero. However, several bodies of evidence (including this paper) point to a potentially irreversible shift toward deciduous secondary forest due to altered disturbance regimes (Schulte and others 2007;Nowacki and Abrams 2008;Rhemtulla and others 2009;Hanberry 2013). Because shifts to deciduous forest cool radiatively through snow-independent forcings and cool within the atmospheric boundary layer in the hottest part of the year, composition shifts have the potential to maintain their cooling role well into the future. Although avoided reforestation or short rotation lengths may cool climate in the short-term, long-term changes in land cover and climate predict a future where passive reforestation to deciduous secondary forest may maximize climatic benefits.

CONCLUSIONS
Reforestation in historically deforested areas is common in much of the temperate northern hemisphere and remains ongoing in more recently deforested areas around the world . Because today's deforestation will have effects that persist for decades to centuries, it was critical to assess the long-term consequences of deforestation against the backdrop of a changing climate. To meet this challenge, we took advantage of a unique dataset, collected both prior to and a century after forest clearing in the Upper Midwest, to assess the long-term effects of deforestation. We show that historical deforestation ultimately results a multitude of land-cover transitions, each with different effects on the energy budget. In conjunction with decreases in snow cover, this mixture of land-cover transitions reduced the cooling effect of deforestation, challenging the notion that temperate deforestation ultimately results in strong albedo-mediated cooling.
In particular, we focused on 'composition shift,' a land-cover transition in which regrown forests are more deciduous than the historical forests. Composition shift occurred in a majority of reforested land, implying that it may be a common long-term ecosystem response to the dual pressures of land use and climatic change (Schulte and others 2007;Goring and others 2016). Composition shift affected the energy budget differently than deforestation, in both the magnitude and seasonality of radiative forcing. Because radiative forcing affects surface conditions both directly and through climatic feedbacks, failing to account for an ultimate transition to 'composition shifted' forest may introduce significant error in long-term predictions (Burakowski and others 2016). Although our results pertain most specifically to areas where regrown forests are deciduous and snow cover is declining, similarly 'nontraditional' long-term land-cover shifts-post-communist agricultural abandonment in Eastern Europe, for example-exist worldwide (Kummerle and others 2009; Duan and others 2011). We show more generally, then, that long-term changes in land cover can have unexpected impacts on surface energy balance that may be overlooked if short-term observations are used to make long-term predictions.

AC KNOWLEDGEMENT S
This material is based upon work supported by the National Science Foundation PalEON MacroSystems Biology under grant DEB #1241874. This work is also supported by the National Science Foundation grant number #1556772 to the University of Notre Dame.

OPEN ACCESS
This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4 .0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.