Albedo-Induced Global Warming Impact at Multiple Temporal Scales within an Upper Midwest USA Watershed

: Land surface albedo is a signiﬁcant regulator of climate. Changes in land use worldwide have greatly reshaped landscapes in the recent decades. Deforestation, agricultural development, and urban expansion alter land surface albedo, each with unique inﬂuences on shortwave radiative forcing and global warming impact (GWI). Here, we characterize the changes in landscape albedo-induced GWI (GWI ∆ α ) at multiple temporal scales, with a special focus on the seasonal and monthly GWI ∆ α over a 19-year period for different land cover types in ﬁve ecoregions within a watershed in the upper Midwest USA. The results show that land cover changes from the original forest exhibited a net cooling effect, with contributions of annual GWI ∆ α varying by cover type and ecoregion. Seasonal and monthly variations of the GWI ∆ α showed unique trends over the 19-year period and contributed differently to the total GWI ∆ α . Cropland contributed most to cooling the local climate, with seasonal and monthly offsets of 18% and 83%, respectively, of the annual greenhouse gas emissions of maize ﬁelds in the same area. Urban areas exhibited both cooling and warming effects. Cropland and urban areas showed signiﬁcantly different seasonal GWI ∆ α at some ecoregions. The landscape composition of the ﬁve ecoregions could cause different net landscape GWI ∆ α .


Introduction
Surface albedo-the amount of solar radiation reflected by a surface relative to the total incident solar radiation-is a fundamental component of the Earth's surface energy balance [1,2].Unlike greenhouse gases (GHGs), which regulate climate by the interception of longwave radiation that affects the Earth's radiation balance, the warming or cooling effects of surface albedo are directly due to instantaneous changes in the amounts of shortwave radiation reflected to outer space.Changes in land use worldwide have greatly reshaped landscapes in recent decades at increasing rates since the Industrial Revolution [3][4][5][6].Earth's surface albedo has changed accordingly, resulting in alterations of the Earth's radiation balance that are partially responsible for the changing climate.
Globally, deforestation, agricultural development-including forest and grassland conversion-and urban expansion are major sources for albedo change [6,7], which, in turn, can directly affect the Earth's radiation balance.Imbalances due to albedo changes are Land 2022, 11, 283 2 of 19 described by the albedo-induced radiative forcing (RF ∆α , W m −2 )-changes in the fraction of solar radiation reflected back to the atmosphere from the Earth's surface [8].For example, according to the Intergovernmental Panel on Climate Change (IPCC; [9]), the RF of wellmixed GHGs has a warming effect equivalent to ~+2.83 W m −2 , while the RF ∆α due to land use and land cover change (LULCC) has a cooling effect equivalent to ~−0.15 W m −2 .In other words, albedo changes have offset ~5% of the energy imbalance caused by well-mixed GHGs, with the offsets varying substantially by region.However, this offset is a global average in reference to LULCC, with dominant changes from forest to non-forest since 1750.Therefore, the local contributions of RF ∆α due to LULCC are unknown and might play an important role in the overall global average of climate regulation effects.
The present scientific understanding of the forcing effects of albedo due to LULCC is ranked as medium-low relative to the rich scientific evidence of the forcing effects of GHGs [1].To cross-examine effects with the climate impact of other GHGs (i.e., biogeochemical GWI), albedo-induced warming or cooling can be converted into equivalents of carbon-dioxide (CO 2eq ) and/or carbon (C eq ) atmospheric radiative forcing via the concept of albedo-induced global warming potential (GWP ∆α , kg CO 2eq m −2 yr −1 ; [9,10]) metric-hereinafter referred to as global warming impact (GWI ∆α ), to be in line with our previous study [11].For example, Houspanossian et al. [12] found that conversion from forests to croplands, from forests to pastures, and from pastures to croplands in dry subtropical forests of South America offset 12-27 Mg C eq ha −1 during a 12-year period, or from 15% to 55% of the total C emissions due to deforestation.In Europe, Carrer et al. [13] reported that inclusion of cover crops in annual cropping systems could have cooling effects equivalent to a mitigation of −0.03 Mg C eq ha −1 yr −1 , while Lugato et al. [14] showed that such mitigation potential due to the inclusion of cover crops could be substantially enhanced by growing high-albedo chlorophyll-deficient cover crops.In southwest Michigan, USA, Chen et al. [15] estimated that land conversion from forest to maize (Zea mays L.) can provide a cooling equivalent to a mitigation of −0.043 Mg C eq ha −1 yr −1 due to a 0.051 (i.e., 5.1%) increase in albedo.At watershed scale, Sciusco et al. [11] demonstrated that altered landscapes could produce cooling effects relative to the intact, native, late successional forests typical of pre-European settlement and contribute a range of −0.1 to −0.4 Mg C eq ha −1 yr −1 , which is the same order of magnitude of biogeochemical GWI emissions due to many crop management components [16,17].
Despite the potential importance of albedo modification strategies for regional to global climate mitigation by IPCC [18] and numerous discussions [19,20] on intentionally increasing surface albedo to cool the Earth, little effort has been made to understand changes in landscape RF ∆α or GWI ∆α in the context of landscape composition at broader temporal scales and for multiple anthropogenic LULCC [11,21,22].A critical unknown is how different cover types contribute to total landscape GWI ∆α at different times of the year.Over the long term (i.e., years to decades), little is known about whether the intra-annual variations of landscape GWI ∆α are significant.
Here, we build on Sciusco et al. [11] to estimate the contributions of the GWI ∆α to the landscape warming or cooling effects at seasonal and monthly timescales over 19 years for multiple ecoregion subtypes.We quantify different cover types in different ecoregions of the upper Midwest USA watershed by estimating (1) the monthly and seasonal contributions to the total landscape cooling or warming; (2) the variations of GWI ∆α contributions by cover type, ecoregion, and year; and (3) the magnitude of cooling or warming effects due to land cover change relative to mature forest cover.

Study Area and Landscape Composition
The Kalamazoo River Watershed (5621 km 2 ; Figure 1) is located in southwest Michigan, USA, and includes portions of 10 counties: Allegan, Barry, Calhoun, Eaton, Hillsdale, Jackson, Kalamazoo, Kent, Ottawa, and Van Buren.The mean annual temperature (1981-2010) is 9.9 • C, and the average annual precipitation is 900 mm evenly distributed throughout the year [23].The dominant cover type prior to European settlement in the early 1800 s was eastern broadleaf deciduous forest [24], with scattered patches of tallgrass prairie, oak savanna, lakes, and wetlands [25].The dominant land cover includes cultivated crops, successional forest stands, pasture-hay grasslands, and two urban areas (i.e., Kalamazoo and Battle Creek).Medium to coarse texture soils and mesic climate allow the continuous recharge of groundwater [26].

Study Area and Landscape Composition
The Kalamazoo River Watershed (5621 km 2 ; Figure 1) is located in southwest Michigan, USA, and includes portions of 10 counties: Allegan, Barry, Calhoun, Eaton, Hillsdale, Jackson, Kalamazoo, Kent, Ottawa, and Van Buren.The mean annual temperature (1981-2010) is 9.9 °C, and the average annual precipitation is 900 mm evenly distributed throughout the year [23].The dominant cover type prior to European settlement in the early 1800 s was eastern broadleaf deciduous forest [24], with scattered patches of tallgrass prairie, oak savanna, lakes, and wetlands [25].The dominant land cover includes cultivated crops, successional forest stands, pasture-hay grasslands, and two urban areas (i.e., Kalamazoo and Battle Creek).Medium to coarse texture soils and mesic climate allow the continuous recharge of groundwater [26].Within the watershed, there are five United States Environmental Protection Agency (US EPA) Level IV ecoregions (Figure 1).Level IV ecoregions have the finest resolution and exhibit unique physiographic, geologic, pedologic, botanic, hydrologic, and climatic characteristics [27].The five ecoregions studied here are Battle Creek Outwash Plain (56b), Michigan Lake Plain (56d), Lake Michigan Moraines (56f), Lansing Loamy Plain (56g), and Interlobate Dead Ice Moraines (56h).For further details, see the US EPA [28].
In this study, we used the National Land Cover Database (NLCD; [29,30]), which provides nine land cover types (barren, cropland, forest, grassland, pasture, shrubland, urban, water, and wetland) at 30 × 30 m spatial resolution, with an overall accuracy ranging between 80% and 90%, in the central and western U.S. [29].Land cover classifications from NLCD are available for the years 2001, 2004, 2006, 2008, 2011, 2013, and 2016.Because land cover maps from NLCD are not annually provided like albedo data (see next section), we assumed land cover is similar to the prior year for years where NLCD data are not available.For example, the land cover map for 2002 was assumed to Within the watershed, there are five United States Environmental Protection Agency (US EPA) Level IV ecoregions (Figure 1).Level IV ecoregions have the finest resolution and exhibit unique physiographic, geologic, pedologic, botanic, hydrologic, and climatic characteristics [27].The five ecoregions studied here are Battle Creek Outwash Plain (56b), Michigan Lake Plain (56d), Lake Michigan Moraines (56f), Lansing Loamy Plain (56g), and Interlobate Dead Ice Moraines (56h).For further details, see the US EPA [28].
In this study, we used the National Land Cover Database (NLCD; [29,30]), which provides nine land cover types (barren, cropland, forest, grassland, pasture, shrubland, urban, water, and wetland) at 30 × 30 m spatial resolution, with an overall accuracy ranging between 80% and 90%, in the central and western U.S. [29].Land cover classifications from NLCD are available for the years 2001, 2004, 2006, 2008, 2011, 2013, and 2016.Because land cover maps from NLCD are not annually provided like albedo data (see next section), we assumed land cover is similar to the prior year for years where NLCD data are not available.For example, the land cover map for 2002 was assumed to be the same as that of 2001, and for 2017-2019, we assumed the land cover had no significant changes.
Growing season and monthly albedos (α gs and α mo , respectively) at 10:30 a.m.local time were derived by stacking (i.e., median image composite) the daily images into growing seasons or months by year.Specifically, α gs accounted for 19 composites (2001-2019), while α mo accounted for 11 composites (January-December, over the 19-year period, less March, which we removed because few images were available and likely due to high cloud cover).We applied the same methodology as Jeong et al. [33] and Sciusco et al. [11] to identify the growing season (roughly from March to November) for each year.Briefly, we used the enhanced vegetation index (EVI; see Appendix A) to identify the growing season by detecting the EVI inflection points (i.e., the dates) when maximum and minimum change rate in greenness occurred over the entire watershed (which was roughly between March and November for the 19-year period).
We employed the Google Earth Engine platform [34] to analyze and process all datasets (see Appendix A).We then performed a zonal statistical analysis within ArcMap (v.10.6) to calculate the proportion of NLCD cover types within each MODIS pixel, and to extract α gs and α mo values by pixel before statistical analysis in RStudio v.1.2.5033 [35].

Albedo-Induced Global Warming Impact (GWI ∆α )
We employed the linear downscaling approach of Chen et al. [36] to estimate surface albedo of cover type i (α si ) at each MODIS pixel.For a MODIS pixel, αsi is considered as the sum of surface albedo of cover type i (α si ) within each MODIS pixel, as follows [36]: where αsi is the estimated surface albedo of cover type i for a time-period (t) (i.e., GS: t = 2001-2019 and monthly: t = January-December, less March), k i is the proportion (0-1) of cover type i in each MODIS pixel, and ε t represents model residuals.
The calculation of landscape albedo-induced radiative forcing (RF ∆α ) and global warming impact (GWI ∆α ) is based on the change in surface albedo due to land cover conversion.This is normally considered as the surface albedo difference (∆α s ) between a cover type i and the native vegetation cover type (i.e., the reference) [11].Here, forest is the dominant land cover type prior to European settlement and serves as our reference [24].Thus, the surface albedo difference of cover type i (∆α si ) for growing season and monthly periods and for each ecoregion at 10:30 a.m.local time is calculated as: where αsi and αsf are the estimated surface albedos of cover type i and of the reference forest f for growing season and monthly periods and for each ecoregion.We calculated ∆α si only when the proportion of a cover type i was ≥80% of the 500 × 500 m MODIS pixel.Conversions to barren, grassland, and shrubland were excluded, as they have a small percentage of the total study area, as were current water and wetland covers.Consequently, only cropland, forest, pasture, and urban covers were considered in this study.
We then used ∆α si to calculate the landscape RF ∆α as follows [13,37,38]: where RF ∆α (W m −2 ) is the landscape albedo-induced radiative forcing at the top-ofatmosphere at 10:30 a.m.local time, and SW in , K T , and ∆α si are the incident shortwave radiation at the surface, the clearness index (for more details, see Equations (A1)-(A5) in Appendix A), and the surface albedo difference between a cover type i and the reference forest, respectively, for the growing season and monthly periods and for each ecoregion (Equation ( 2)).The incident shortwave radiation at the surface (SW in ) and the clearness index (K T ) were derived from the solar and meteorological dataset NASA POWER [39] at daily time-step for multiple locations (i.e., five Level IV ecoregions) within the Kalamazoo River Watershed.We then averaged SW in and K T values to match the 19 growing seasons and the 11 months.Positive or negative values of RF ∆α indicate warming and cooling effects, respectively.Lastly, we calculated the landscape GWI ∆α as follows [13,37,38]: where  [13,41,42] and TH is fixed at 100 years (i.e., the number of time steps the GWI is then divided by) [43,44].With Equation (4), we calculate the equivalent RF ∆α that a unit area of A would have at the global scale.Positive or negative values of the GWI ∆α indicate effects equivalent to CO 2 emission or mitigation, respectively.Here, we report the results of landscape seasonal and monthly GWI ∆α (GWI ∆αgs and GWI ∆αmo , respectively) expressed with units of Mg C eq ha −1 gs −1 and Mg C eq ha −1 mo −1 , respectively, which refer to 10:30 a.m.local time.

Contributions of Land Cover Change to GWI ∆α
We performed a nested analysis of variance (ANOVA) with repeated measurements (i.e., growing season and monthly periods) to quantify the contribution to landscape GWI ∆α by cover type and ecoregion for the growing season and monthly periods.Specifically, we looked at contributions among and within the five ecoregions, with two linear models: GW I ∆α (t) = cover type(t) where ecoregion refers to the five Level IV ecoregions (i.e., 56b, 56d, 56f, 56g, and 56h), and cover type refers to the three cover types (i.e., cropland, pasture, and urban) used to determine the albedo difference from the reference forest (Equation ( 2)) for growing season and monthly periods at each ecoregion.For further details about the ANOVA analysis, see Appendix A.
The average GWI ∆α showed an overall cooling effect (Figure 3a,b) for most cover types, with the exception of urban, which showed neutral effects.Overall, cropland cover type had seasonal and monthly average cooling effects equivalent to −0.35 ± 0.05 Mg C eq ha −1 gs −1 and −0.68 ± 0.61 Mg C eq ha −1 mo −1 , respectively (Figure 3a,b).The highest seasonal and monthly cooling effects reached −0.47 Mg C eq ha −1 gs −1 (in 2015 for Ecoregions 56b, 56d, 56f, and 56h; Tables S3.1-3.3, and 3.5 and −2.15 Mg C eq ha −1 mo −1 (in February for Ecoregion 56b; Table S4.1), for the two periods, respectively.These cooling effects represented ~26% and ~68% more than the seasonal and monthly annual averages, respectively.On the other hand, urban cover was the only cover type showing neutral effects (i.e., seasonal and monthly average effects equivalent to −0.001 ± 0.034 Mg C eq ha −1 gs −1 and −0.158 ± 0.256 Mg C eq ha −1 mo −1 , respectively; Figure 3a,b).S4.1), for the two periods, respectively.These cooling effects represented ~26% and ~68% more than the seasonal and monthly annual averages, respectively.On the other hand, urban cover was the only cover type showing neutral effects (i.e., seasonal and monthly average effects equivalent to −0.001 ± 0.034 Mg Ceq ha −1 gs −1 and −0.158 ± 0.256 Mg Ceq ha −1 mo −1 , respectively; Figure 3a,b).S4.1), for the two periods, respectively.These cooling effects represented ~26% and ~68% more than the seasonal and monthly annual averages, respectively.On the other hand, urban cover was the only cover type showing neutral effects (i.e., seasonal and monthly average effects equivalent to −0.001 ± 0.034 Mg Ceq ha −1 gs −1 and −0.158 ± 0.256 Mg Ceq ha −1 mo −1 , respectively; Figure 3a,b).All three cover types showed similar trends in GWI ∆αgs 2001-2019 across the five ecoregions (Tables S3.1-3.5 and Figure 4a) with similar deviations from the GWI ∆αgs mean.
The growing season cooling effects (GWI ∆αgs ) ranged between −0.39 and −0.47 Mg C eq ha −1 gs −1 for cropland, between −0.25 and −0.31 Mg C eq ha −1 gs −1 for pasture, and between −0.01 and −0.11Mg C eq ha −1 gs −1 for urban.Urban was the only cover type that also showed warming effects, with GWI ∆αgs raging between 0.01 and 0.05 Mg C eq ha −1 gs −1 .The inter-monthly variation of GWI ∆αmo for cropland, pasture, and urban cover showed similar trends with higher cooling effects in January, February, and December (Tables S4.1-4.5 and Figure 4b).Among the three cover types, the highest cooling occurred in cropland (−0.18 to −2.15 Mg C eq ha −1 mo −1 ), followed by pasture (−0.13 to −1.64 Mg C eq ha −1 mo −1 ) and urban (−0.02 to −0.470 Mg C eq ha −1 mo −1 ) cover.For the three cover types, GWI ∆αmo was relatively constant from April to November.However, urban had slightly bell-shaped trends with small warming effects in June and July (0.01 to 0.1 Mg C eq ha −1 mo −1 ; Tables S4.1-4.5 and Figure 4b), and cropland had an inverted bell-shaped trend with slight rises in June and October.Lastly, pasture (in Ecoregions 56b, 56f, and 56g) had relatively constant trends.The variation among ecoregions in GWI ∆αgs was significant (p < 0.001) by ecoregion, cover type, and their interactions (ANOVA model in Equation ( 5); Table S5), while the variation in GWI ∆αmo was significant (p < 0.001) only by cover type.Neither years nor months were significant for GWI ∆α variations.Most of the variation in GWI ∆αgs was explained by cover type (η 2 = ~95%), followed by the interaction between ecoregions and cover type (η 2 = ~52%) and ecoregions (η 2 = 32%).In comparison, the variation in GWI ∆αmo was almost equally explained by ecoregion, cover type, and their interactions, although only cover type was significant (η 2 = ~24% at p < 0.001).The variation in both GWI ∆αgs and GWI ∆αmo within ecoregions (Equation ( 6); Table S5), however, was significant (p < 0.001) by cover type, which explained more of the variation in GWI ∆αgs (η 2 = 99%) than in GWI ∆αmo (η 2 = 65%).S3.1-3.5 and Figure 4a) with similar deviations from the GWIΔαgs mean.The growing season cooling effects (GWIΔαgs) ranged between −0.39 and −0.47 Mg Ceq ha −1 gs −1 for cropland, between −0.25 and −0.31 Mg Ceq ha −1 gs −1 for pasture, and between −0.01 and −0.11Mg Ceq ha −1 gs −1 for urban.Urban was the only cover type that also showed warming effects, with GWIΔαgs raging between 0.01 and 0.05 Mg Ceq ha −1 gs −1 .The intermonthly variation of GWIΔαmo for cropland, pasture, and urban cover showed similar A post-hoc Tukey test analysis (Figure 5a,b) showed that, within each ecoregion, the least square means (LSMs) of GWI ∆αgs had low variability and were significantly different among the cover types (Figure 5a).The LSMs for GWI ∆αmo were more variable, and many cover types had statistically similar means (Figure 5b).Among ecoregions, the LSMs of cropland GWI ∆αgs at Ecoregion 56g were significantly different from those at Ecoregions 56b and 56d, while the LSMs at urban GWI ∆αgs at Ecoregion 56g were significantly different from those at Ecoregion 56d.However, no significant differences in their LSMs for GWI ∆αmo were observed (Figure 5b).
A post-hoc Tukey test analysis (Figure 5a,b) showed that, within each ecoregion, the least square means (LSMs) of GWIΔαgs had low variability and were significantly different among the cover types (Figure 5a).The LSMs for GWIΔαmo were more variable, and many cover types had statistically similar means (Figure 5b).Among ecoregions, the LSMs of cropland GWIΔαgs at Ecoregion 56g were significantly different from those at Ecoregions 56b and 56d, while the LSMs at urban GWIΔαgs at Ecoregion 56g were significantly different from those at Ecoregion 56d.However, no significant differences in their LSMs for GWIΔαmo were observed (Figure 5b).The overall GWIΔα contribution from different seasons and months varied by cover type, and it was exclusively higher during the non-growing season (NGS) than during the growing season (GS) months for all ecoregions (Table S(6) and Figure 6), with the NGS months being characterized by only cooling effects (Table S(7)).As a general trend, the highest contributions were in February and the lowest in October.During the NGS, urban (at all ecoregions) contributed the most to the total cooling effect (between 18% and 31%), followed by pasture (at Ecoregions 56b, 56f, and 56g; contribution between 14% and 25%), Lowercase letters indicate differences among cover types within the five ecoregions, while uppercase letters indicate differences of same cover type among the five ecoregions.The among ecoregions analysis only considered cropland and urban covers (i.e., the cover types that were in every ecoregion).
The overall GWI ∆α contribution from different seasons and months varied by cover type, and it was exclusively higher during the non-growing season (NGS) than during the growing season (GS) months for all ecoregions (Table S6 and Figure 6), with the NGS months being characterized by only cooling effects (Table S7).As a general trend, the highest contributions were in February and the lowest in October.During the NGS, urban (at all ecoregions) contributed the most to the total cooling effect (between 18% and 31%), followed by pasture (at Ecoregions 56b, 56f, and 56g; contribution between 14% and 25%), and cropland (at all ecoregions; contribution between 14% and 24%).It is worth noting that during GS months, no cover type had a contribution >8% (i.e., 1/12 of the annual total).Nevertheless, climate regulations of urban (for all ecoregions) were close to the overall mean value of 8% in April, while cropland (at all ecoregions) and pasture (at Ecoregions 56b, 56f, and 56g) were close in April and May.
and cropland (at all ecoregions; contribution between 14% and 24%).It is worth noting that during GS months, no cover type had a contribution >8% (i.e., 1/12 of the annual total).Nevertheless, climate regulations of urban (for all ecoregions) were close to the overall mean value of 8% in April, while cropland (at all ecoregions) and pasture (at Ecoregions 56b, 56f, and 56g) were close in April and May.

Discussion
Our results showed that the albedo-induced global warming impact (GWIΔα) accounted for a significant portion of the climate cooling effects (i.e., Ceq mitigation) due to land cover changes and landscape composition.Individual contributions varied by cover type, ecoregion, and season/month, with cropland showing the highest cooling effects, followed by pasture.Urban showed both cooling and warming effects, the latter of which occurred only during growing season months.Overall, the cooling effects of the monthly GWIΔα were higher than the seasonal ones, most likely due to the substantial influence of snow cover on land surface albedo (i.e., high presence of snow during the winter months; [45]) combined with the effect of management practices on different vegetation surfaces.For example, snowfalls on harvested crop fields create a highly reflective layer, whereas on forested fields, snowfalls tend to be masked by the tree structures and reflect less.For

Discussion
Our results showed that the albedo-induced global warming impact (GWI ∆α ) accounted for a significant portion of the climate cooling effects (i.e., C eq mitigation) due to land cover changes and landscape composition.Individual contributions varied by cover type, ecoregion, and season/month, with cropland showing the highest cooling effects, followed by pasture.Urban showed both cooling and warming effects, the latter of which occurred only during growing season months.Overall, the cooling effects of the monthly GWI ∆α were higher than the seasonal ones, most likely due to the substantial influence of snow cover on land surface albedo (i.e., high presence of snow during the winter months; [45]) combined with the effect of management practices on different vegetation surfaces.For example, snowfalls on harvested crop fields create a highly reflective layer, whereas on forested fields, snowfalls tend to be masked by the tree structures and reflect less.For the same reason, seasonal analysis showed that the cooling contributions during the non-growing season months were higher than during the growing season months.Overall, our results seem to be promising in the context of climate regulation of albedo changes due to land cover changes and landscape composition.Nevertheless, several assumptions and limitations to our study could benefit GWI computations elsewhere.

Cooling Effects
During the 19-year study period, the highest cropland cooling effect of the seasonal and monthly GWI ∆α was equivalent to −0.47 Mg C eq ha −1 gs −1 in 2015 for Ecoregions 56b, 56d, 56f, and 56h, and −2.15 Mg C eq ha −1 mo −1 in February for Ecoregion 56b, respectively.In comparison, Abraha et al. [46], accounting for GHGs using whole-system lifecycle analysis, found emissions of 2.6 Mg C eq ha −1 yr −1 over 8 years in Conservation Reserve Program (CRP) grasslands converted to maize.Thus, the seasonal and monthly maximum albedoinduced cooling effects from cropland represent (i.e., offset) 18% and ~83%, respectively, of the annual C eq over the 8 years from CRP grasslands converted to maize fields.The total highest albedo cooling effect due to both the seasonal and monthly GWI ∆α would completely offset the annual emissions due to the grassland converted to maize over the 8 years.Moreover, the abovementioned cropland cooling effects are more than enough to offset the annual net biogeochemical GWI-i.e., warming effects due to the net contributions from CO 2 , methane (CH 4 ), and nitrous-oxide (N 2 O) at 0.31 Mg C eq ha −1 yr −1 -produced by annual crop systems (i.e., maize-soybean-wheat rotation) under conventional tillage management of the same area [47].

Variable Effects
Urban areas appeared to have either cooling or warming effects depending on the temporal scale examined.Unlike cropland, which had cooling effects during summer months mostly due to changes in vegetation over the growing season [48][49][50][51][52], urban cover had warming effects over the growing season months and cooling effects during the rest of the year.At Ecoregion 56g, the highest seasonal and monthly warming effects, estimated at 0.05 Mg C eq ha −1 gs −1 in 2002 and 2012, and ~0.1 Mg C eq ha −1 mo −1 in June and July, are equivalent to ~29% and 45%, respectively, of the annual net ecosystem production (NEP) of deciduous forest stands in northern Michigan [53].On the other hand, the maximum monthly cooling effect of urban landscapes at Ecoregion 56h was equivalent to −0.7 Mg C eq ha −1 mo −1 (in February).In comparison, Xu et al. [38] estimated a C eq offset between −2.2 Mg C eq ha −1 and −4.4 Mg C eq ha −1 induced by an increase in pavement albedo of 0.01 (i.e., 1%) across two major US cities over a 50-year period.Such results are important, considering that the watershed includes two major urban centers, Kalamazoo and Battle Creek, with a total population of >500,000 people in 2010 [36].More so, as previous studies predicted [54], Michigan urban areas is anticipated to increase >50% by 2030.

Intra-and Inter-Annual Variability of Albedo and GWI ∆α
Despite our expectation that the intra-and inter-annual variability of surface albedo would vary due to seasonality and climatic conditions [11], we did not find significant differences in inter-annual variation of ∆α nor GWI ∆α .For example, each cover type showed unique inter-annual trends that appeared to be similar during the seasonal and monthly periods.However, we found that during the growing season, the GWI ∆α of cropland and urban cover types was not the same in every ecoregion, which emphasizes that cover types may have different contributions depending on the location.In other words, changes in landscape composition in the five ecoregions could cause different net landscape GWI ∆α .For example, contrasting landscape compositions among the five ecoregions led to different cumulative cooling effects during the growing season over 19 years.This value varied from −6.89 Mg C eq ha −1 gs −1 to −11.36 Mg C eq ha −1 gs −1 at Ecoregions 56d and 56b, respectively.At monthly scales, landscape composition produced cumulative cooling effects between −9.14 Mg C eq ha −1 mo −1 and −14.89Mg C eq ha −1 mo −1 at Ecoregions 56h and 56b, respectively.
Our results also suggest that a total forest loss of ~680 ha due to the conversion to cropland and urban (i.e., the main cover type classes within the watershed) during 2001−2016 led to seasonal and monthly cooling effects at the watershed scale that, on average, were equivalent to ~−179 Mg C eq gs −1 and ~−374 Mg C eq mo −1 , respectively.These amounts equal approximately 90 ha and 190 ha, respectively, of mature forest net carbon sequestration in the same region, assuming an average value of 2 Mg C ha −1 yr −1 .

Seasonal Percent Contributions to the Total Cooling and Warming
In line with other studies [55], the largest contribution to the overall total seasonal GWI ∆α came from the non-growing season months, during which all the cover types exhibited cooling effects that varied in magnitude depending on the ecoregion.Once again, urban was the only cover type that contributed to warming effects in the growing season, generally following a decreasing trend going from June to September.Such results echo the need reported in previous work to advance research on the importance of surface albedo modification within urban components (e.g., pavements, roofs, walls) as climate regulation strategy to resolve the urban energy budget and energy demand [38,56].However, it should be noted that urban areas are composed of infrastructure (e.g., roofs, walls, pavements, etc.) and other cover types (e.g., trees, grasses, bare soil, water bodies, etc.) with varying albedo contributions, which were not considered in this study (see next section for more details).Nevertheless, the seasonal analysis clearly confirmed that the contributions to the total landscape cooling or warming effect varied by ecoregion.We found that albedo climate benefits either contributed net cooling/warming or a net neutral effect of the at landscape scale depending on the ecoregion.

Assumptions and Limitations of the Study
Several assumptions and limitations in our study could benefit computations of albedoinduced global warming impact (GWI ∆α ) elsewhere.The first assumption is related to the choice of the time horizon (TH, see Equation ( 4)) fixed at a 100-year period.The choice of either short or long time horizons can either over-or de-accentuate GWI ∆α values [57].Specifically, by keeping TH fixed at 100 years, we assume that the land cover composition of the study area will remain the same for the next 100 years, although it is likely that the land cover over the next 100 years will be very different.However, by setting TH = 100, we aligned our study with the Kyoto Protocol [11,44], and hence with the IPCC protocols.
There is also uncertainty associated with the datasets employed in this study.For example, the MODIS MCD43 albedo product has a pixel with a nominal spatial resolution of 500 × 500 m, which has been shown not to properly match the effective spatial resolution (usually much higher than the nominal one [58,59]).However, previous attempts by researchers to analyze the effective spatial resolution of the MODIS albedo product [60,61] were limited to a single homogeneous area.For areas characterized by substantial land surface heterogeneity, similar to the one presented in this work, the effective spatial representativeness of the pixel is hard to determine [58].Regarding land cover classification assumptions, we acknowledge that NLCD data are not annually available like albedo data.To overcome this limitation, when NLCD data were not available for a particular year, we assumed land cover was similar to the previous year available.Although the inter-annual LULCC for those years may not be captured, in this way, we were able to carry on a longer timeseries analysis of albedo which otherwise would have been limited to only 7 years (i.e., 2001, 2004, 2006, 2008, 2011, 2013, and 2016).A second assumption regards our analysis of urban areas.The NLCD provides four sub-classes of impervious surfaces (i.e., <20%, between 20% and 49%, between 50% and 79%, and between 80% and 100% of the total cover) composing the developed (i.e., urban) class.However, in this study, we aggregated the four sub-classes into one single "urban" land cover class.We are aware that by doing so, we may have obtained less than precise values of albedo change in urban areas.Other studies [62] have demonstrated that urban heat island intensity (UHII) is highly sensitive to the spatial context of urban areas (i.e., urban, rural, and their combination).Despite this, the focus of our study was not to investigate the importance of urban albedo modifications in the context of UHII effects.We acknowledge that a clear definition of urban and rural contexts would benefit future investigations of albedo changes in urban areas.
We investigated the contribution of landscape composition due to land transformation on GWI ∆α in the context of climate change mitigations by considering the forest cover type as the reference for the entire study area [24].However, without other references, our calculations cannot estimate the forests' contribution to GWI ∆α .This hinders a comprehensive synthesis of the total climate emission or mitigation of the watershed, considering that the low albedo of forests contributes to climate warming [63].Moreover, regarding indirect biophysical effects, forests' role in mitigating climate change is multifaceted.Recent studies have demonstrated how re-/afforestation strategies can increase low-level cloud cover formation, which, depending on the forest type, results in cooling effects [64].Analyzing land transformation with reference to forests is only one method.For example, other studies and policymakers have focused on land transformation in the context of bioenergy conversions [55] and land management practices to compare landscape dynamics to agriculture [65][66][67].
We also acknowledge that the instantaneous albedo values at 10:30 a.m.local time are likely different from the daily averages.In this study, we utilized albedo estimated by the MODIS BRDF function, which is the composite of a 16-day period [68], with albedo values from a single snapshot at 10:30 a.m.local time (e.g., MODIS Terra morning overpassing time).However, there is increasing evidence showing diurnal variations of albedo [68,69] under different sky conditions.Nevertheless, the MODIS MCD43 albedo product is soundly validated [70][71][72], as well as widely accepted for retrieving albedo from other remote sensing products [32,[73][74][75][76][77] and presenting overall good accuracy compared to in situ daily averages [70].
Lastly, we considered the growing season and monthly albedo at 10:30 a.m.local time by computing the median composite (see Appendix A.1), which prevented us from accounting for the effects of land surface characteristics (i.e., vegetation properties, such as leaf area index, and landscape heterogeneity) on the spatiotemporal variation of albedo among and within patches of the same type.This represents a limitation, as during the growing season, vegetation cover and canopy structure and albedo are negatively correlated due to the varying capacity of the canopy to absorb incoming solar radiation [78].Our cover type categories did not reflect these differences, so future efforts will be needed to quantify such differences, including the use of other remote sensing metrics and instantaneous measurements [79].

Conclusions
Albedo-induced global warming impact (GWI ∆α ) accounted for significant climate benefits (i.e., C eq mitigation) due to land cover changes and structural variations across the landscape.Looking at individual cover types, the climate benefits were higher in cropland, with seasonal and monthly offsets of 18% and 83%, respectively, of the annual greenhouse gas emissions of maize fields in the same area.The second-highest benefits were found in pasture lands.However, urban showed near-neutral albedo climate benefits.Notably, the overall change in landscape composition within the five ecoregions caused different net landscape C eq mitigations.Seasonal climate benefits ranged from −6.89 Mg C eq ha −1 gs −1 at Ecoregion 56d to −11.36 Mg C eq ha −1 gs −1 at Ecoregion 56b, and monthly climate benefits ranged from −9.14 Mg C eq ha −1 mo −1 at Ecoregion 56h to −14.89 Mg C eq ha −1 mo −1 at Ecoregion 56b.Changes in albedo due to land cover changes and landscape composition are of fundamental importance in the context of landscape climate regulation.We must couple these estimates with biogeochemical (i.e., GHGs) climate benefits to increase our understanding of the magnitude that various human-induced mechanisms contribute to climate change.Table S5.Nested analysis of variance (ANOVA with repeated measurements), among and within ecoregions, based on the linear models (see Equations ( 5) and ( 6) in the main text).Dependent variables: GWIΔαgs (Mg Ceq ha -1 gs -1 ) and GWIΔαmo (Mg Ceq ha -1 gs -1 ).

Figure 1 .
Figure 1.The study area of the five United States Environmental Protection Agency Level IV ecoregions and the nine National Land Cover Database cover type classes within the Kalamazoo River Watershed in 2001.(Map projection: WGS-84 UTM Zone 16 N).

Figure 1 .
Figure 1.The study area of the five United States Environmental Protection Agency Level IV ecoregions and the nine National Land Cover Database cover type classes within the Kalamazoo River Watershed in 2001.(Map projection: WGS-84 UTM Zone 16 N).

Figure 2 .
Figure 2. Albedo difference (∆α) between a cover type i and the forest for (a) growing season (∆α gs ) and (b) monthly (∆α mo ) periods, within the five Level IV ecoregions and the entire Kalamazoo River Watershed 2001-2019.Annotations indicate mean (±standard deviation) of ∆α gs and ∆α mo .
Land 2022, 11, x FOR PEER REVIEW 8 of 20 All three cover types showed similar trends in GWIΔαgs 2001-2019 across the five ecoregions (Tables

Figure 4 .
Figure 4. Albedo-induced global warming impact (GWIΔα) for a given cover type within the five Level IV ecoregions.Panels (a,b) represent the GWIΔα for the growing season (GWIΔαgs; 2001-2019) and monthly (GWIΔαmo; January-December, less March) periods, respectively.Positive and negative values of GWIΔα indicate warming and cooling effects, respectively, equivalent to Carbon (Ceq) emission and mitigation, respectively.

Figure 4 .
Figure 4. Albedo-induced global warming impact (GWI ∆α ) for a given cover type within the five Level IV ecoregions.Panels (a,b) represent the GWI ∆α for the growing season (GWI ∆αgs ; 2001-2019) and monthly (GWI ∆αmo ; January-December, less March) periods, respectively.Positive and negative values of GWI ∆α indicate warming and cooling effects, respectively, equivalent to Carbon (C eq ) emission and mitigation, respectively.

Figure 5 .
Figure 5.The least square means (LSMs) multi-comparison analysis of albedo-induced global warming impact (GWIΔα) for a given cover type across the five Level IV ecoregions for the (a) growing season (GWIΔαgs; 2001-2019) and (b) monthly (GWIΔαmo; January-December) periods.Whiskers represent the lower and upper limits of the 95% family-wise confidence level of the LSMs.Bars sharing the same letters are not significantly different according to the post-hoc Tukey test analysis.Lowercase letters indicate differences among cover types within the five ecoregions, while uppercase letters indicate differences of same cover type among the five ecoregions.The among ecoregions analysis only considered cropland and urban covers (i.e., the cover types that were in every ecoregion).

Figure 5 .
Figure 5.The least square means (LSMs) multi-comparison analysis of albedo-induced global warming impact (GWI ∆α ) for a given cover type across the five Level IV ecoregions for the (a) growing season (GWI ∆αgs ; 2001-2019) and (b) monthly (GWI ∆αmo ; January-December) periods.Whiskers represent the lower and upper limits of the 95% family-wise confidence level of the LSMs.Bars sharing the same letters are not significantly different according to the post-hoc Tukey test analysis.Lowercase letters indicate differences among cover types within the five ecoregions, while uppercase letters indicate differences of same cover type among the five ecoregions.The among ecoregions analysis only considered cropland and urban covers (i.e., the cover types that were in every ecoregion).

Figure 6 .
Figure 6.Percent contribution of albedo-induced global warming impact (GWIΔα) to cooling (values in bold) or warming effects by season and month periods for major cover types in the five Level IV ecoregions.The horizontal dashed line represents the average contribution ~8% (i.e., 1/12 of the annual total), and the solid vertical lines separate the non-growing season (NGS) from the growing season (GS) months.Values of March were missing and gap-filled as the mean GWIΔα of February and April.

Figure 6 .
Figure 6.Percent contribution of albedo-induced global warming impact (GWI ∆α ) to cooling (values in bold) or warming effects by season and month periods for major cover types in the five Level IV ecoregions.The horizontal dashed line represents the average contribution ~8% (i.e., 1/12 of the annual total), and the solid vertical lines separate the non-growing season (NGS) from the growing season (GS) months.Values of March were missing and gap-filled as the mean GWI ∆α of February and April.

Figure S1 .
Figure S1.Albedo-induced radiative forcing (RFΔα) for the five Level IV ecoregions.Panels a,b represent the RFΔα for the growing season (RFΔαgs; 2001-2019) and monthly (RFΔαmo; January-December, less March) periods, respectively.Positive and negative values of RFΔαgs and RFΔαmo indicate warming and cooling effects, respectively, during 19 growing seasons and 11 months over the 19-year period, respectively.

Table 1 .
Land cover composition in ha (%) of the Kalamazoo River Watershed, and net gain (+) and loss (−) for each cover type during the period 2001-2016.

Table 2 .
Pivot table showing the land cover conversion (ha) of each cover type during the period 2001-2016 across the Kalamazoo River Watershed.Bold values indicate the main land cover conversions, while "*" indicates no land cover conversion.

Table S2 .
Mean (± standard deviation) for α and α 2001-2019 for the four cover types.Average values (X ) across the four cover types (X cover type) and for the 19-years growing season (X gs) and the 11-month (X mo) periods are also shown.

Table S4 . 1 .
Monthly albedo-induced global warming impact (GWIΔαmo, Mg Ceq ha -1 mo -1 ) January-December for Ecoregion 56b.Cumulative GWIΔαmo value for the entire ecoregion over the 11-month period is underlined.Positive and negative values of GWIΔαmo indicate monthly warming and cooling effects, respectively, equivalent to Carbon (Ceq) emission and mitigation, respectively, during the 19-year period.

Table S4 . 2 .
Monthly albedo-induced global warming impact (GWIΔαmo, Mg Ceq ha -1 mo -1 ) January-December for Ecoregion 56d.Cumulative GWIΔαmo value for the entire ecoregion over the 11-month period is underlined.Positive and negative values of GWIΔαmo indicate monthly warming and cooling effects, respectively, equivalent to Carbon (Ceq) emission and mitigation, respectively, during the 19-year period.

Table S4 . 3 .
Monthly albedo-induced global warming impact (GWIΔαmo, Mg Ceq ha -1 mo -1 ) January-December for Ecoregion 56f.Cumulative GWIΔαmo value for the entire ecoregion over the 11-month period is underlined.Positive and negative values of GWIΔαmo indicate monthly warming and cooling effects, respectively, equivalent to Carbon (Ceq) emission and mitigation, respectively, during the 19-year period.

Table S4 . 4 .
Monthly albedo-induced global warming impact (GWIΔαgs, Mg Ceq ha -1 mo -1 ) January-December for Ecoregion 56g.Cumulative GWIΔαmo value for the entire ecoregion over the 11-month period is underlined.Positive and negative values of GWIΔαmo indicate monthly warming and cooling effects, respectively, equivalent to Carbon (Ceq) emission and mitigation, respectively, during the 19-year period.

Table S4 . 5 .
Monthly albedo-induced global warming impact (GWIΔαmo, Mg Ceq ha -1 mo -1 ) January-December for Ecoregion 56h.Cumulative GWIΔαmo value for the entire ecoregion over the 11-month period is underlined.Positive and negative values of GWIΔαmo indicate monthly warming and cooling effects, respectively, equivalent to Carbon (Ceq) emission and mitigation, respectively, during the 19-year period.
2: generalized eta squares indicating the variance in the dependent variable (i.e., GWIΔα) accounted for by the independent variables (i.e., ecoregion, cover type, and their interactions)." a ": η 2 values obtained from the Mauchly's test of sphericity by applying the Greenhouse-Geisser correction. η

Table S6 .
Overall percent contribution and percent contribution higher than the average (>8%; i.e., 1/12 th of the total, i.e., equal contribution across the 12-month period) of albedo-induced global warming impact (GWIΔα) by season and month in the five Level IV ecoregions.Cooling effects showing the contributions higher than the average are underlined.

Table S7 .
Percent contribution higher than the average (>8%; i.e., 1/12 th of the total, i.e., equal contribution across the 12-month period) of albedo-induced global warming impact (GWIΔα) by season and month for cropland, pasture, and urban in the five Level IV ecoregions.

Table S8 . 1 .
Growing season albedo-induced radiative forcing (RFΔαgs, W m -2 ) 2001-2019 for Ecoregion 56b.Cumulative RFΔαgs value for the entire ecoregion over the 19-year period is underlined.Positive and negative values of RFΔαgs indicate growing season warming and cooling effects, respectively.

Table S8 . 2 .
Growing season albedo-induced radiative forcing (RFΔαgs, W m -2 ) 2001-2019 for Ecoregion 56d.Cumulative RFΔαgs value for the entire ecoregion over the 19-year period is underlined.Positive and negative values of RFΔαgs indicate growing season warming and cooling effects, respectively.

Table S8 . 3 .
Growing season albedo-induced radiative forcing (RFΔαgs, W m -2 ) 2001-2019 for Ecoregion 56f.Cumulative RFΔαgs value for the entire ecoregion over the 19-year period is underlined.Positive and negative values of RFΔαgs indicate growing season warming and cooling effects, respectively.

Table S9 . 1 .
Monthly albedo-induced radiative forcing (RFΔαmo, W m -2 ) January-December for Ecoregion 56b.Cumulative RFΔαmo value for the entire ecoregion over the 11-month period is underlined.Positive and negative values of RFΔαmo indicate monthly warming and cooling effects, respectively.

Table S9 . 2 .
Monthly albedo-induced radiative forcing (RFΔαmo, W m -2 ) January-December for Ecoregion 56d.Cumulative RFΔαmo value for the entire ecoregion over the 11-month period is underlined.Positive and negative values of RFΔαmo indicate monthly warming and cooling effects, respectively.

Table S9 . 3 .
Monthly albedo-induced radiative forcing (RFΔαmo, W m -2 ) January-December for Ecoregion 56f.Cumulative RFΔαmo value for the entire ecoregion over the 11-month period is underlined.Positive and negative values of RFΔαmo indicate monthly warming and cooling effects, respectively.

Table S9 . 5 .
Monthly albedo-induced radiative forcing (RFΔαmo, W m -2 ) January-December for Ecoregion 56h.Cumulative RFΔαmo value for the entire ecoregion over the 11-month period is underlined.Positive and negative values of RFΔαmo indicate monthly warming and cooling effects, respectively.