Forest Fire Effects on Landscape Snow Albedo Recovery and Decay

: Surface snow albedo (SSA) darkens immediately following a forest ﬁre, while landscape snow albedo (LSA) brightens as more of the snow-covered surface becomes visible under the charred canopy. The duration and variability of the post-ﬁre snow albedo recovery process remain unknown beyond a few years following the ﬁre. We evaluated the temporal variability of post-ﬁre snow albedo recovery relative to burn severity across a chronosequence of eight burned forests burned from 2000 to 2019, using pre- and post-ﬁre daily, seasonal, and annual landscape snow albedo data derived from the Moderate Resolution Imaging Spectroradiometer (MOD10A1). Post-ﬁre annual LSA increased by 21% the ﬁrst year following the ﬁre and increased continually by 33% on average across all eight forest ﬁres and burn severity classiﬁcations over the period of record (18 years following a ﬁre). Post-ﬁre LSA measurements increased by 63% and 53% in high and moderate burn severity areas over ten years following ﬁre. While minimum and maximum snow albedo values increased relative to annual post-ﬁre LSA recovery, daily snow albedo decay following fresh snowfall accelerated following forest ﬁre during the snowmelt period. Snow albedo recovery over 10 years following ﬁre did not resemble the antecedent pre-ﬁre unburned forest but more resembled open meadows. The degradation of forest canopy structure is the key driver underlying the paradox of the post-ﬁre snow albedo change (SSA vs. LSA). severity, demonstrating the dependence of LSA on forest density. Post-ﬁre LSA persistently increased by 1.3 per year following a ﬁre, while spring LSA measurements increased by 1.0 per year following a ﬁre event, suggesting a potential inﬂuence of post-ﬁre radiative forcing on surface snow albedo during the springtime snowmelt phase.


Introduction
Seasonal snowpack is a critical water resource and contributes over half the total runoff and as much as 70% of streamflow in mountainous areas of the western US [1]. Snowpack acts as a seasonal high mountain water reservoir, recharging aquifers and nourishing the rivers, streams, and tributaries into the summer and fall following snowmelt [1][2][3], providing important moisture for soils and vegetation. A warming climate has reduced snowpack volume [4][5][6], advanced the timing of spring snowmelt [5,7,8], limited downstream water availability [1,7,9], and ultimately reduced summer soil moisture [10]. Moreover, due to climate warming, forest fires are increasing in intensity, extent, frequency, and duration across the western US, particularly during early snowmelt years [11][12][13].
In the western US, most forest fires occur in the seasonal snow zone [14]. Forest fires leave behind charred forests, which stand blackened for decades shedding black carbon (BC) and burned woody debris (BWD) into the snowpack below for many winters following fire. BC and BWD concentrate on the snowpack surface during the melt, visibly darkening the snowpack and lowering the surface snow albedo (SSA) [14,15]. Forest fires drive increases in net snowpack shortwave radiation as a result of the reduction in postfire SSA in combination with the more open post-fire forest canopy, which allows more incoming solar radiation incidents on the snowpack surface [15,16]. Darkening of SSA due Temporal variability of post-fire LSA depends on the intra-annual seasonality of LSA, the delayed mortality of burned forests, and forest regeneration over the years following the fire. Disturbance history, drought stress, pathogens such as bark beetle, competition, season, and soil type all interact and influence post-fire tree mortality [40][41][42], making delayed tree mortality particularly challenging to predict [43]. Patterns of fireinduced tree mortality and post-fire vegetation re-development affect the chance of forest regeneration and growth [44]. Recovery of post-fire landscape snow albedo to pre-fire conditions depends primarily on the complex interactions driving forest regeneration to pre-fire antecedent conditions [45], and the long-term effects of fire on radiative forcing and climate continue to be heavily studied. In many high-severity fires, forested areas cannot recover to their pre-fire conditions, remaining a heterogeneous mix of shrubland [46]. This alters the fire risk of the areas in that high severity burns are more likely in the less fire tolerant conditions. This pattern of high severity fires leading to future high severity re-burns has been seen in multiple areas across the Western US [47]. Fire-induced tree mortality and delayed tree mortality are complex issues caused by direct and indirect factors that require further investigation [48].
It is essential to characterize forest fire effects on LSA recovery and decay to understand the variability and magnitude of long-term implications of forest fires on climate, hydrology, and water resource availability. No previous study has evaluated the temporal variability of forest fire effects on landscape snow albedo relative to year-since-fire or burn severity. Remote sensing imagery provides the only spatially distributed long-term data Spatial variability of post-fire LSA immediately following fire and post-fire LSA recovery over time (Figure 1), depends on antecedent (pre-fire) forest composition and structure, as well as burn severity and related canopy mortality over time directly, due to the proportion of visible snow-covered surface, and indirectly, through canopy shadowing and darkening snow surfaces [39]. As forest fires increase in extent, severity, duration, and frequency, understanding these forest fire effects on snow albedo will be even more crucial to understanding its impact on overall climate cooling.
Temporal variability of post-fire LSA depends on the intra-annual seasonality of LSA, the delayed mortality of burned forests, and forest regeneration over the years following the fire. Disturbance history, drought stress, pathogens such as bark beetle, competition, season, and soil type all interact and influence post-fire tree mortality [40][41][42], making delayed tree mortality particularly challenging to predict [43]. Patterns of fire-induced tree mortality and post-fire vegetation re-development affect the chance of forest regeneration and growth [44]. Recovery of post-fire landscape snow albedo to pre-fire conditions depends primarily on the complex interactions driving forest regeneration to pre-fire antecedent conditions [45], and the long-term effects of fire on radiative forcing and climate continue to be heavily studied. In many high-severity fires, forested areas cannot recover to their pre-fire conditions, remaining a heterogeneous mix of shrubland [46]. This alters the fire risk of the areas in that high severity burns are more likely in the less fire tolerant conditions. This pattern of high severity fires leading to future high severity re-burns has been seen in multiple areas across the Western US [47]. Fire-induced tree mortality and delayed tree mortality are complex issues caused by direct and indirect factors that require further investigation [48].
It is essential to characterize forest fire effects on LSA recovery and decay to understand the variability and magnitude of long-term implications of forest fires on climate, hydrology, and water resource availability. No previous study has evaluated the temporal variability Remote Sens. 2022, 14, 4079 4 of 17 of forest fire effects on landscape snow albedo relative to year-since-fire or burn severity. Remote sensing imagery provides the only spatially distributed long-term data available to investigate the spatial and temporal variability of forest fire effects on LSA. Using daily global snow albedo data from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument from 2000 to 2018, we evaluated the temporal variability of preand post-fire landscape snow albedo recovery and decay relative to burn severity, tree mortality, and years-since-fire. Quantifying the temporal variability of forest fire effects on snow albedo will assist in understanding the hydrologic and climatic implications of this feedback loop and how the future water vulnerability and availability will be impacted by increasing forest fires and decreasing snow-water storage under future warming climate conditions [49][50][51].

Materials and Methods
Our approach combined remote sensing, downscaling, and statistical analysis to evaluate pre-vs. post-fire variability in LSA across a chronosequence of eight mixed severity burned forests in western Wyoming. We evaluated MODIS-derived winter and spring LSA relative to burn severity for years following the fire. To isolate forest fire effects on landscape snow albedo and omit interannual variability of snow accumulation and melt patterns, we normalized snow albedo data from within the burn perimeter with data from the 5-km buffered area outside the burn perimeter. Buffer areas outside the burn perimeter accounted for interannual variability and adjusted for normal year-to-year fluctuations in snow accumulation and melt patterns.

Study Area
We examined post-fire snow albedo change in a chronosequence of eight burned forests in the seasonal snow zone of western Wyoming, which burned from 2000 to 2018 with mixed burn severity ( Figure 2). The study area extends around Jackson, Wyoming, and across the Triple Divide, the headwaters of three major river systems, the Colorado, Columbia, and Missouri Rivers. The elevation of Jackson Hole Airport is 1966 m. The climate in Jackson, Wyoming, is cold and temperate, with an average annual temperature of 3.6 • C and annual precipitation of 851 mm. The cold season lasts approximately December through March each year, where December is the coldest with an average temperature of −8.8 • C and experiences the most precipitation with an average of 95 mm (climate-data.org, accessed on 10 November 2020).
This area is vital because it supplies much of the western US with snow-water resource storage. While as little as 37% of the entire precipitation in the western US falls as snow [1], across this mountainous region, approximately 60-80% of precipitation falls as snow [52]. The headwaters of these river systems begin high in the Rocky Mountains, where snow accounts for a large portion of the overall water accumulation that flows downstream into the Pacific Ocean, Gulf of California, and the Gulf of Mexico. The area comprises pine-dominated forests, predominately lodgepole pine (Pinus contorta) and whitebark pine (Pinus albicaulis), with a history of regular and frequent forest fires across the season snow zone.
The chronosequence of burned forests of progressing ages provides a study area with varying temporal features of post-fire impacts on LSA while reducing the spatial variability between sites due to a relatively similar forest structure and composition and topographic complexity across the study area. The eight fires evaluated in this study occurred between 2000 and 2018 in and near the Hoback River Basin. They contained continuous and mixed patches of unburned, low, moderate, and high burn severity areas (Table 1).  The chronosequence of burned forests of progressing ages provides a study area with varying temporal features of post-fire impacts on LSA while reducing the spatial variability between sites due to a relatively similar forest structure and composition and topographic complexity across the study area. The eight fires evaluated in this study occurred between 2000 and 2018 in and near the Hoback River Basin. They contained continuous and mixed patches of unburned, low, moderate, and high burn severity areas (Table 1).

Remote Sensing Data
Daily LSA was derived from MODIS MOD10A1 dataset at 500-m spatial resolution; the MODIS/Terra Snow Cover Daily L3 Global 500 m SIN Grid, Version 6 [53]. Winter (January and February) and spring (March and April) snow albedo data were obtained from the two spatial tiles that included the chronosequence of eight burned forests from 2000 to 2019. We restricted data acquisition to these dates to evaluate the winter period of snow accumulation and the spring period of spring snow ablation. The seasonal snow albedo is typically brighter during winter and darker during spring due to grain growth and LAP concentration on the surface during snowmelt, ultimately driving seasonal variability in the radiative heat forcing on snow [54,55]. LSA values in the MOD10A1 dataset were calculated as the surface reflectance within snow-covered areas integrated from narrowband albedo to broadband albedo after atmospheric and forest reflectance properties [56]. Unlike SSA, measured or modeled at the snow surface, LSA is measured from a birds-eye-view on a satellite-based instrument that integrates variability across the surface into a broadscale landscape snow albedo value.
Yearly forest density was derived from the percent tree cover parameter in the MODIS MOD44B dataset at a 250-m spatial resolution; the MODIS/Terra Vegetation Continuous Fields (VCF) Yearly L3 Global 250 m SIN Grid, Version 6 [57]. Forest fire perimeters and gridded burn severity data at a 30-m spatial resolution were collected from the MTBS project [43]. Within the burn perimeters, burn severity was classified into high burn severity, moderate burn severity, and low/unburned forests from the differenced relativized differenced normalized burn ratio (RdNBR) in MTBS [43]. To distinguish open meadow from unburned forest within the burn perimeters and 5-km buffer areas, the National Landcover Database (NLCD) landcover classifications were used to differentiate unburned forest and open areas at a 30-m spatial resolution. The NLCD serves as the Landsat-based, 30-m resolution land cover database for the United States and provides spatial reference and data for land surface characteristics [58]. Excluding bodies of water, perennial ice, and snow cover, each landcover classification was categorized as either forest or open area. All areas classified as deciduous, mixed, and coniferous forests were categorized as forests. All areas classified as developed, barren, shrubland, herbaceous, planted/cultivated, and wetland were categorized as open areas.

MODIS Downscale
Broad-scale remote sensing data such as from MODIS, including snow albedo at a 500 m and forest density at a 250 m spatial resolution, was too coarse relative to the actual spatial variability to accurately evaluate the influence of burn severity on snow albedo following a forest fire. In order to reduce the large uncertainties persistent with MODIS snow albedo data primarily due to the coarse spatial resolution and mixed pixels, we downscaled the broad-scale snow albedo and forest density data using the 30 m burn severity data to~150 m spatial resolution. We used a resampling method as the base procedure in performing an effective downscale of the MODIS data [59,60]. The MOD10A1 and MOD44B datasets were resampled via bilinear interpolation from the Landsat-derived burn severity classifications using the values of the four nearest cell centers to determine the value of the resampled MOD10A1 and MOD44B datasets. We performed this procedure on every daily snow albedo and annual forest density measurement within each forest fire area and burn severity classification. Particularly in areas of complex landcover patchiness, the discrepancies in the spatial scale of the processes influencing snow albedo to the scale of the measurements introduce uncertainty even after being transformed via downscaling methods [61].

SNOTEL Data and Days since Snowfall
Daily snow depth (cm) data were retrieved from snow telemetry (SNOTEL) stations using the Natural Resources Conservation Services (NRCS) Data Retrieval Tool from 1 January 2001-2 April 2021. The SNOTEL network is an automated system of snowpack and related micro-meteorological sensors operated by the NRCS of the United States Department of Agriculture. SNOTEL sites throughout the Western US and Alaska typically collect continuous snow water equivalent (SWE), total precipitation, snow depth, and air temperature data [52]. The meteorological data was first retrieved by performing a query using Google Earth Engine. The SNOTEL stations contained in or near the Hoback River Basin were used in this study, including Blind Bull Summit (353) Pass (822), and Triple Peak (831). For the entire period of record from the MODIS data, days since snowfall was derived by calculating the number of days since a daily snow depth accumulation greater than 5 cm.

LSA Recovery and Decay Statistical Analysis
Average winter, spring, and annual snow albedo values and standard deviation were calculated for each burn severity classification for each fire. In total, we evaluated a total of 24,605 snow albedo values spread across eight fires and four burn severity classifications; 13,674 pre-fire values and 10,931 post-fire values, and an additional total of 11,744 daily snow albedo values spread across the unburned forest and open meadow in the 5-km buffered areas; 6502 pre-fire values and 5242 post-fire values. Average winter, spring, and annual snow albedo metrics within each fire were normalized by the average winter, spring, and annual snow albedo values from within the 5-km buffered area. The buffer forest area snow albedo metrics were used to normalize the high severity, moderate severity, and unburned forest areas in the burned forest, while the buffer open area was used to normalize the open area within the burn perimeter. Average annual forest density was calculated for each burn severity classification for each fire for the entire period of record from 2000-2019. Post-fire percent differences in snow albedo and forest density were calculated for all burn severity classifications and years since the fire. Normalized pre-vs. post-fire LSA comparisons were compared for winter, spring, and annual snow albedo values for each burned forest in the chronosequence using Welch's two-sample t-test.
Multivariate linear regression models were used to characterize annual normalized LSA (nLSA) in the burn perimeter as a function of forest density and years-since-fire, relative to burn severity classification. Spline regression models were developed using the splines package in R to investigate key stages and potential thresholds in the post-fire snow albedo recovery. Each spline model function was partitioned with a 75-25 split, such that 75% of the data was used to train the models, and 25% of the data was used to test each model. Manual knots were identified to minimize the residual sum of squares on validation data and used in the spline regression model to represent the major stages in post-fire LSA recovery. To characterize post-fire snow albedo decay functions for snow accumulation and snow ablation periods, winter and spring daily snow albedo values were plotted over average days since snowfall for the study region. Exponential decay functions were used to model the post-fire snow albedo decay relative to the key post-fire snow albedo recovery stages. Exponential decay models were computed using a generalized linear model with a Gaussian log-link. ArcGIS Pro 2.7.0 and R Studio software version 1.2.5 [62,63] were used for all data analysis. Statistical relationships were tested with a significance level of 0.05.

Annual Landscape Snow Albedo Recovery following Forest Fire
Annual post-fire LSA increased the first winter following forest fire by 21% on average (range of 15-28%, mean increase relative to pre-fire LSA = 6.6, sd = 12.8, p-value < 0.0001). Across all eight fires, LSA brightened by 33% on average (range of 8-39% increase) over the 18 years following the fire and the period of record in this study ( Table 2, mean increase relative to pre-fire LSA = 10.1, sd = 11.9, p-value < 0.0001). Annual LSA recovery occurred over the 10 years following the fire, while annual LSA continued to brighten across all eight burned forests (Figure 3). After a decade following the fire, no significant change in annual post-fire LSA recovery was observed.
Post-fire LSA during winter was brighter by 19% overall than post-fire LSA during spring (Figure 4; winter mean = 54.4, sd = 9.0; spring mean = 44.0, sd = 7.5). The difference between winter and spring post-fire LSA increased with burn severity, demonstrating the dependence of LSA on forest density. Post-fire winter LSA persistently increased by 1.3 per year following a fire, while spring LSA measurements increased by 1.0 per year following a fire event, suggesting a potential influence of post-fire radiative forcing on surface snow albedo during the springtime snowmelt phase. Table 2. Pre-fire and post-fire LSA statistics for all eight separate fires, including Pre-fire count (n) of MODIS cells used in the analysis, pre-fire mean LSA, pre-fire standard deviation (sd) of LSA, post-fire count (n) of MODIS cells used in the analysis, post-fire mean LSA, post-fire standard deviation (sd) of LSA, the mean difference, and post-fire total percent increase for each fire. The sample size was determined by the number of cloudless (50% or less) MODIS scenes available since 2000 relative to the year the fire occurred. The significance reported with the mean difference was calculated using Welch's two-sample t-test, *** = p-value < 0.001.  Table 2. Pre-fire and post-fire LSA statistics for all eight separate fires, including Pre-fire count (n) of MODIS cells used in the analysis, pre-fire mean LSA, pre-fire standard deviation (sd) of LSA, post-fire count (n) of MODIS cells used in the analysis, post-fire mean LSA, post-fire standard deviation (sd) of LSA, the mean difference, and post-fire total percent increase for each fire. The sample size was determined by the number of cloudless (50% or less) MODIS scenes available since 2000 relative to the year the fire occurred. The significance reported with the mean difference was calculated using Welch's two-sample t-test, *** = p-value < 0.001. Post-fire LSA during winter was brighter by 19% overall than post-fire LSA during spring (Figure 4; winter mean = 54.4, sd = 9.0; spring mean = 44.0, sd = 7.5). The difference between winter and spring post-fire LSA increased with burn severity, demonstrating the dependence of LSA on forest density. Post-fire winter LSA persistently increased by 1.3 per year following a fire, while spring LSA measurements increased by 1.0 per year following a fire event, suggesting a potential influence of post-fire radiative forcing on surface snow albedo during the springtime snowmelt phase.   Areas with high and moderate burn severity experienced the greatest change in prefire vs. post-fire normalized LSA, undergoing a 44% and 36% increase, respectively, over the period of record (Table 3; p-value < 0.001). Within the burn perimeter, unburned forest and open meadow areas increased by 32% and 22%, respectively, over the period of record (Table 3; p-value < 0.001). Annual LSA values were not significantly different between the forested burn severity classes within the burned perimeter during the pre-fire or post-fire periods. Post-fire comparisons within the burn perimeter demonstrate that all burn severity classes brighten after the fire, and the forest classes, including "unburned", "high severity", and "moderate severity" classes, were similar. In contrast, the open meadow within the burn perimeter was consistently brighter following the fire. Overall, brightening occurred over a decade following the fire until LSA in all forest classes within the burn perimeter resembled the LSA in the open meadow ( Figure S1). Areas with high and moderate burn severity experienced the greatest change in prefire vs. post-fire normalized LSA, undergoing a 44% and 36% increase, respectively, over the period of record (Table 3; p-value < 0.001). Within the burn perimeter, unburned forest and open meadow areas increased by 32% and 22%, respectively, over the period of record (Table 3; p-value < 0.001). Annual LSA values were not significantly different between the forested burn severity classes within the burned perimeter during the pre-fire or post-fire periods. Post-fire comparisons within the burn perimeter demonstrate that all burn severity classes brighten after the fire, and the forest classes, including "unburned", "high severity", and "moderate severity" classes, were similar. In contrast, the open meadow within the burn perimeter was consistently brighter following the fire. Overall, brightening occurred over a decade following the fire until LSA in all forest classes within the burn perimeter resembled the LSA in the open meadow ( Figure S1). Table 3. Post-fire comparisons of percent LSA increase relative to burn severity classification, and pre-vs. post-fire comparisons of percent LSA increase relative to burn severity. Results from the daily normalized LSA Tukey Honestly Significant Difference test based on the piecewise comparisons of burn severity classification in post-fire comparisons and pre-vs. post-fire comparisons. *** = p-value < 0.001.

Post-Fire Annual Landscape Snow Albedo Recovery
Post-fire landscape snow albedo increased persistently for more than ten years following a fire as a function of forest density, years-since-fire, and burn severity (p < 0.0001, R 2 = 0.68; nLSA = 12.2 + 1 years-since-fire − 0.4 × forest density + 12.2 × high + 0.2 × moderate + 0.7 × unburned − 10.1 × open). Spline regression models demonstrated the rate of post-fire LSA recovery occurred in four main stages, including (1) an immediate post-fire increase in snow albedo (0-1 years), (2) an early post-fire period of slow snow albedo recovery (1-6 years), (3) an intermediate post-fire period of faster snow albedo recovery (6-10 years), and (4) a "final" post-fire period of snow albedo stabilization after 10 years following the fire. The models associated with these main recovery periods take the form of: fire LSA recovery occurred in four main stages, including (1) an immediate post-fire increase in snow albedo (0-1 years), (2) an early post-fire period of slow snow albedo recovery (1-6 years), (3) an intermediate post-fire period of faster snow albedo recovery (6-10 years), and (4) a "final" post-fire period of snow albedo stabilization after 10 years following the fire. The models associated with these main recovery periods take the form of: where: x = years since fire, a & b = model coefficients.

Post-Fire Burn Severity and Annual Landscape Snow Albedo Recovery
Landscape snow albedo across the burned forest area showed similar long-term albedo values and overall recovery patterns regardless of the burn severity classification ( Figure 5). Specifically, LSA in high and medium severity burned forests showed no significant difference from one another for up to 15 years following the fire. Within the fire perimeters, LSA in low severity burned forests was lower than LSA in high severity burned forests only in the first 6 years following fire and showed no significant difference from moderate severity burned regions over the entire post-fire period. Open meadow areas within the fire perimeters were consistently brighter than forested areas in low, medium, and high severity burned areas over the entire post-fire period. Over the decade observed following the fire, LSA in high and moderate severity burned forests converged on LSA in the open meadow ( Figure 5). High and moderate severity burned forests were much darker immediately following forest fire than open meadows outside the fire perimeter. In contrast, over the post-fire snow albedo recovery period 10 years following the fire, the LSA in the high and moderate severity burned forests transitioned to near zero difference with open meadows outside the burn perimeter ( Figure S1).

Post-Fire Daily Landscape Snow Albedo Decay
Post-fire daily landscape snow albedo decay (over days since snowfall) occurred more rapidly than in the pre-fire period during spring ablation (Figure 6, pre-fire vs. post-fire during accumulation, n = 12,537, Wald Z = 49.82, p-value ≤ 0.0001; pre-fire vs. post-fire during ablation, n = 12,066, Wald Z = 24.92, p-value ≤ 0.0001). Post-fire daily LSA decay demonstrated steeper exponential decay coefficients than pre-fire daily LSA decay during the snowmelt period (pre-fire accumulation (LSA = 0.397e −0.00705 days ), post-fire accumulation (LSA = 0.535e −0.0067 days ), pre-fire ablation (LSA = 0.283e −0.00408 days ), post-fire ablation (LSA = 0.392e −0.00957 days ). Post-fire daily LSA decay models were different during the initial and final snow albedo decay recovery periods but not during the intermediate recovery period. Therefore, they were aggregated into one primary decay model for the pre-fire and post-fire periods. Pre-and post-fire LSA maximum and minimum values, derived from LSA decay models, demonstrate how daily LSA decay relates to decadal scale post-fire LSA recovery. Post-fire LSA maximum values were brighter following fire across LSA decay during accumulation. However, post-fire LSA maximum values converge during ablation on the lower values of pre-fire LSA after weeks of post-fire LSA decay. Lower overall maximum and minimum daily LSA values relative to snow albedo decay were observed in the initial stage of post-fire LSA recovery than in the final stage ( Figures 5 and S1), reflective of the broader post-fire landscape snow albedo recovery over time.

Discussion
Across the chronosequence of eight burned forests in western Wyoming, LSA immediately brightened the first year following the fire and continued to increase for at least ten years following the fire. In the burn perimeters, a mosaic of post-fire forest structure was associated with heterogeneous antecedent conditions and post-fire burn severity. Within the burn perimeters, the forested areas, including high and moderate burn sever- Figure 6. Pre-and post-fire landscape snow albedo decay during accumulation and ablation derived from MODIS data over days since snowfall > 5 cm derived from the average snowfall of all SNOTEL sites within the study domain. Box plots show the distribution of snow albedo pixels within burned forests for each 4-day period following a fresh snowfall.

Discussion
Across the chronosequence of eight burned forests in western Wyoming, LSA immediately brightened the first year following the fire and continued to increase for at least ten years following the fire. In the burn perimeters, a mosaic of post-fire forest structure was associated with heterogeneous antecedent conditions and post-fire burn severity. Within the burn perimeters, the forested areas, including high and moderate burn severities and unburned forests, experienced a more significant change in post-fire LSA than open areas. Over a decade post-fire, snow albedo changes do not demonstrate significant recovery to the pre-fire unburned forest but more resemble an open meadow.

Post-Fire Landscape Snow Albedo Recovery
Like surface snow albedo, landscape snow albedo recovery occurs over a decade following a forest fire, with the degradation of forest canopy structure as the key driver underlying the paradox of the post-fire snow albedo change. The immediate removal in canopy density caused an immediate increase in exposure of the snow surfaces and higher post-fire LSA values. Over the years following a forest fire, continuing forest mortality and shedding of the forest canopy reduced forest density [64], leading to the recovery in post-fire LSA over a decade following fire. The post-fire LSA recovery stabilized after ten years-since-fire as LSA flattened out, suggesting a decade-long post-fire LSA recovery period ( Figure 4). Across all eight fires, forest density was at its lowest approximately fifteen years following a fire before showing indications of potential forest recovery.
Four main stages of post-fire LSA recovery occurred across the chronosequence of burned forests, including (1) an immediate post-fire increase in snow albedo (0-1 years), (2) an early post-fire period of slow snow albedo recovery (1-6 years), (3) an intermediate post-fire period of faster snow albedo recovery (6-10 years), and (4) a "final" post-fire period of stabilization after 10 years following fire ( Figure 4b). Spline regression model results provide an understanding of the stages and duration of post-fire snow albedo recovery which should inform hydro-climatological modeling.
Over the decade of observed post-fire recovery, the difference in mean LSA between burned forests and open meadows shrinks in magnitude over years-since-fire indicating that the mean LSA of burned forests becomes more similar to the mean LSA of unburned open meadows as the post-fire forest recovers ( Figure 5). Conversely, the difference in mean LSA between burned and unburned forests grow in magnitude over the years since the fire, indicating that the mean LSA of burned forests becomes less similar to the mean LSA of unburned forests over the years since the fire ( Figure S1). Together, these results show that the LSA of snowpack in burned forests does not recover back to the LSA of unburned forests after 10 years following fire and instead recovers to the LSA of open meadows 10 years post-fire.

Uncertainties in Post-Fire Landscape Snow Albedo Recovery
We quantified the spatio-temporal relationship between post-fire LSA recovery and forest density. However, many complex measures influence recovery after ten years following a fire. In many areas throughout the seasonal snow zone, post-fire climate conditions are likely to become increasingly unfavorable to tree regeneration, even if seed sources are nearby [65,66]. High severity burned areas are less likely to see complete forest regeneration, with probable conversion to non-forest landscape [65,66]. Studies have suggested the increase in LSA following a forest fire may be significant enough to neutralize the initial carbon release caused by the fire and thus may not necessarily accelerate climate warming to a significant degree [31,34]. There are negative hydrologic implications of a forested area being converted to a non-forested area. Post-fire landscape heterogeneity continues to increase with delayed tree mortality and regeneration rates, and it is unclear when or if LSA in a burned forest returns to the levels of pre-fire conditions. While signs of regeneration were visible in the MOD44B data approximately fifteen years following a fire, there may have been a delayed response of decreasing post-fire LSA back to pre-fire conditions due to MOD44B measuring on an annual temporal scale. Vegetation regeneration and regrowth were most likely prevalent and visible in the late spring-and summertime. Past studies with MODIS products suggest that the difference between pre-fire and post-fire spring-and summertime surface albedo increases with time during the first five to eight years following a fire [31,45]. In snow-dominated, mountainous areas, late-spring snow cover may influence the appearance of vegetation recovery.
To remove the interannual variability of snow conditions from the post-fire snow albedo recovery patterns, we normalized LSA within the burned areas by the LSA within the 5 km buffer area surrounding the burned areas. Open areas in the buffer zones were brighter even before the fire occurred than those in the burned perimeter, possibly due to developed areas or high shrublands in the buffer zones that are different. For this reason, and also due to the mixed pixels inherent in coarse resolution MODIS data, we observed negative nLSA values in the burned forest prior to forest fire occurrence (Figure 4).
There was an observed period approximately from seven-to ten years following the fire where the LSA in the high burn, moderate burn, and unburned forest areas became brighter than the LSA in the buffer open area. This was unexpected, as the buffer open area was predicted to represent the brightest areas of measured LSA. However, when the buffer areas were constructed from the NLCD dataset, multiple landcover identifications classified as an open area might have made it significantly darker than an open meadow in the wilderness. For example, developed areas, which are typically impervious and dark surfaces, and shrublands, which were described as areas of two classes, one of which being shrubs no higher than three meters tall. Both landcover types were classified as buffer open areas in the study, while they have a probability of darkening the LSA more than an open meadow. For this reason, it was possible for the LSA in the burned perimeter to appear brighter in the measured data than that of the buffer open area.

Post-Fire Landscape Snow Albedo Decay
Post-fire snow albedo was higher overall than pre-fire snow albedo following fresh snowfall during both the accumulation and ablation periods ( Figure 6). During accumulation, pre-and post-fire snow albedo decay rates were similar, but during ablation, post-fire snow albedo decay was much faster than pre-fire snow albedo decay. In early winter, as snow accumulates, solar incidence angle is low, and temperature gradients within the snowpack are typically low, so snow albedo decay is slower. While in early spring, as the snow melts, the solar incidence angle is higher, and temperature gradients in the snowpack are greater. In contrast, the post-fire snowpack surface is darkened by black carbon and burned woody debris, lowering snow albedo overall and accelerating snow albedo decay rates. The spatio-temporal variability in post-fire LSA decay observed in this study using coarse resolution remote sensing data aligns with previous research which used in situ SSA measurements [15]. These post-fire LSA decay functions may be useful in parameterizing the spatio-temporal variability of snow albedo in hydrologic modeling.

Uncertainties in Post-Fire Landscape Snow Albedo Decay
The spatial heterogeneity caused by forest fires introduces new complexities influencing snow accumulation and ablation in burned forests. Throughout this study, we found the variability of daily post-fire LSA measurements increased compared to the daily pre-fire LSA. The spatial scale of MOD10A1 made it difficult to differentiate the burn severity classifications in this matter. Therefore, we could not statistically differentiate between the post-fire LSA of high burn, moderate burn, and unburned forest areas. A contributing factor in this may have been the incorporation of low-burned areas into our unburned classification, although low-burned areas likely produce little forest damage, mostly affecting the understory. However, MOD10A1 captured the increased effect forest fires have on springtime snowmelt vs. winter accumulation.
At the large spatial resolution presented by the MOD10A1 pixels, key snow-forest processes are occurring at the sub-pixel resolution that influenced the results of LSA that we most likely were not able to account for. This was the main purpose of downscaling the MODIS resolution. Pre-fire forest conditions are generally more homogenous, in that forested areas are dominated by vegetation and open areas are more dominated by shrubland or grassland. However, the landscape mosaic becomes much patchier and heterogenous following a forest fire, increasing the likelihood that landcover will be misclassified at the sub-pixel level. For instance, MODIS pixels assigned to a high burn severity designation most likely also contain patches of low or moderate burn severity. Therefore, the snow albedo calculated from coarse spatial resolution MODIS pixels over high burn severity area is underrepresented and lower than it typically would be when determined by the finer spatial resolution of an instrument like Landsat-8 [67].

Conclusions
Post-fire landscape snow albedo increases up to 10 years following a fire over four main stages of recovery due to decreasing forest density and increasing exposure of the snowpack. After the initial stage of immediate brightening, post-fire LSA continues to increase. Still, the brightening is greatest between five and ten years following a fire event, possibly due to the transition between initial and delayed tree mortality. After ten years post-fire, post-fire LSA recovery stabilized where the burned areas were more like an open area than the unburned forest before the fire occurrence.
As climate warming continues, seasonal snowpacks will continue to decline, resulting in earlier snowmelt and lower annual landscape albedo. Moreover, due to warming, forest fires will continue to increase in intensity, size, frequency, and duration, impacting the environment and human life [68,69]. Particularly in the western US, the increase in forest fires and rising temperatures will continue to deplete seasonal snowpack, making water availability an increasing issue. Post-fire snow albedo recovery is critical to improving models estimating snow albedo-climate feedbacks, snow-water storage, water resource availability, and vulnerability.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/rs14164079/s1, Figure S1: Tukey plot shows significance in pairwise differences in mean landscape snow albedo data between burned forest (high, moderate, low burn severity classifications) relative to open meadow in the burned regions and forest and open meadow in the buffer outside the burned regions.