Spatio-temporal patterns and trends in MODIS-retrieved radiative forcing by snow impurities over the Western US from 2001 to 2022

The seasonal mountain snowpack of the Western US (WUS) is a key water resource to millions of people and an important component of the regional climate system. Impurities at the snow surface can affect snowmelt timing and rate through snow radiative forcing (RF), resulting in earlier streamflow, snow disappearance, and less water availability in dry months. Predicting the locations, timing, and intensity of impurities is challenging, and little is known concerning whether snow RF has changed over recent decades. Here we analyzed the relative magnitude and spatio-temporal variability of snow RF across the WUS at three spatial scales (pixel, watershed, regional) using remotely sensed RF from spatially and temporally complete (STC) MODIS data sets (STC-MODIS Snow Covered Area and Grain Size/MODIS Dust Radiative Forcing on Snow) from 2001 to 2022. To quantify snow RF impacts, we calculated a pixel-integrated metric over each snowmelt season (1st March–30th June) in all 22 years. We tested for long-term trend significance with the Mann–Kendall test and trend magnitude with Theil–Sen’s slope. Mean snow RF was highest in the Upper Colorado region, but notable in less-studied regions, including the Great Basin and Pacific Northwest. Watersheds with high snow RF also tended to have high spatial and temporal variability in RF, and these tended to be near arid regions. Snow RF trends were largely absent; only a small percent of mountain ecoregions (0.03%–8%) had significant trends, and these were typically decreasing trends. All mountain ecoregions exhibited a net decline in snow RF. While the spatial extent of significant RF trends was minimal, we found declining trends most frequently in the Sierra Nevada, North Cascades, and Canadian Rockies, and increasing trends in the Idaho Batholith. This study establishes a two-decade chronology of snow impurities in the WUS, helping inform where and when RF impacts on snowmelt may need to be considered in hydrologic models and regional hydroclimate studies.


Introduction
Mountain snowpack of the Western US (WUS) is a key part of the regional hydroclimate system and acts as a primary water resource to millions of people and a diverse natural ecosystem.Seasonal snow cover dominates the hydrology of this region, providing drinking water, hydropower, and irrigation to an increasingly water stressed environment (Li et al 2017).Impurities at the snow surface, termed light absorbing particles (LAPs), are impactful to seasonal snowmelt in tandem with rising anthropogenic activity and climate change, as these together alter the timing and rate of snowmelt.Drivers of snowmelt must be accurately predicted and characterized to optimize management of this hydroclimate system.When LAPs are not considered, operational hydrologic models will underestimate melt rates (Skiles et al 2012), negatively affecting the reliability of operational streamflow forecasting (Bryant et al 2013).Furthermore, to understand changing drivers of snowmelt in sensitive mountain systems, it is critical to characterize LAP distributions and trends for historical reanalysis and future projections of snowpack and the climate system (Skiles et al 2018, Réveillet et al 2022).
LAPs absorb incoming solar radiation strongly in the visible range of the electromagnetic spectrum, thereby darkening the snow surface (i.e.reducing snow albedo) and absorbing more energy that warms and melts the snowpack (if isothermal).In areas with significant LAP concentrations, this process, known as snow radiative forcing (RF), is a key driver of the snow surface energy balance, snow temperature, early melt onset, increased melt rates, and decreased snow persistence in spring and summer (Painter et al 2012, Skiles et al 2018).Direct effects of LAPs to the snow come from surface darkening, while indirect effects are seen in accelerated snow grain growth and exposure of darker snow-free land, further reducing albedo through positive feedbacks (Painter et al 2012).This potential to impact snowmelt rates raises questions about the role of LAPs in regional snowmelt spatiotemporal patterns and trends, in tandem with the effects of climate change (Deems et al 2013).
Dust, black carbon (BC), microbial growth, and potentially microplastics are the primary types of LAPs responsible for driving increased melt (Skiles and Painter 2019).Dust is the most common atmospheric aerosol by mass, and can be found in snow cover downwind of arid and semi-arid environments having undergone disturbances resulting in a decline of natural biophysical soil crusts (Painter et al 2007, 2010, Mauro et al 2015).Mineral dust concentrations are typically higher than that of BC in seasonal snow (Skiles and Painter 2018), but BC is a much more effective absorber of incoming solar radiation in the visible wavelengths (0.38-0.7 nm) (Wiscombe and Warren 1980).BC is created by incomplete fuel combustion, and is associated with areas of high population or wildfires.Snow downwind of urban areas can have high BC concentrations (Ming et al 2009), proven to significantly reduce snow albedo (Hadley and Kirchstetter 2012).Both BC and dust may be present in a snowpack, and recent studies have presented analyses quantifying their relative and combined impacts on snow RF (compare Skiles andPainter 2018, Gleason et al 2022).Microbial presence on snow can also have substantial RF effects.Recent research in the Cascades demonstrates that snow algae can increase RF during peak bloom in mid-latitude snow through their snow-darkening effects (Healy and Khan 2023).Biotic presence in mountain snowpack may become a growing issue in the future, as the habitat extent of snow algae may expand with a warming climate, lengthened melt season, and increased nutrient arrival at high elevations (Benning et al 2014).Finally, microplastics (particles ⩽ 5 mm) have recently been detected in snowpack in increasing quantities, but their associated effects on RF are still being studied (Reynolds et al 2022).
The effects of snow RF on snowmelt and streamflow timing manifest in unique ways for each region and microclimate (e.g. total snowfall accumulated in winter, springtime solar radiation regime).In the Upper Colorado Basin, springtime episodes of rapid dust deposition can magnify snow RF and snowmelt rates; in extreme conditions, these may advance snowpack disappearance by over 50 d (Skiles et al 2012), alter the rate and timing of seasonal streamflow (Painter et al 2018), and reduce annual streamflow volume an estimated 5% due to increased evapotranspiration losses with a longer snow-free season (Painter et al 2010).Recent dust-on-snow in the Upper Colorado is impacting snowmelt more than recent changes in climate (Skiles et al 2012, Deems et al 2013).Dust is expected to increase with further aridification of the WUS as a result of ongoing climate change, increasing LAPs at the snow surface, and contributing to the positive feedback loop of climate change (Mahowald et al 2010, Skiles et al 2018).In High Mountain Asia, dust has been found to substantially contribute to snowmelt and rising snowlines with climate warming (Sarangi et al 2020).In European mountains, dust and BC have been shown to advance peak streamflow by 10-15 d and shorten the snow season by an average of 17 d (Réveillet et al 2022).
LAP's have been reported in WUS snowpacks with a high level of spatiotemporal variability (Doherty et al 2016, Skiles et al 2018, Cui et al 2021).Each region/subregion of the WUS is uniquely impacted by LAP's; sources (e.g.dust events, fuel emissions, wildfires, etc.) are variable, along with entrainment mechanisms, physical geographic characteristics, and meteorological patterns.Controls of LAP deposition vary in space and time, and shift with land use, soil and vegetation states, and anthropogenic climate, resulting in low predictability of LAP deposition and seasonal RF impacts (Painter et al 2007, 2010, Mahowald et al 2010, Skiles et al 2018).Much inconsistency is observed as the number and intensity of dust deposition events greatly varies regionally and with each season (Skiles et al 2015).
The current understanding of LAP concentrations and snow RF in snowpack is based on spatially limited and few snow measurements, often collected from areas of relatively high atmospheric deposition (Skiles et al 2015, Chellman et al 2017, Skiles and Painter 2018).Additional understanding of LAP impacts on snowmelt has come with numerical modeling studies (e.g.Oaida et al 2015).Recent work has utilized satellite remote sensing to assess snow RF from LAPs over a 16 year period across the Northern Hemisphere (e.g.Cui et al 2021Cui et al , 2023)), which have revealed differences at sub-continental scales.No study has yet established and compared the geographic distribution of remotely sensed snow RF during the snowmelt season over the WUS across a range of spatial scales (e.g.pixel to large basin) and subregions.While estimates of snow RF have been reported (Hadley et al 2010, Huang et al 2022, Painter et al 2012, Skiles et al 2012, Seidel et al 2016), these have generally been in confined geographic areas (plot to watershed scale) and for short records.Little is known about longer term interannual variability and possible trends in snow RF over recent decades and at different scales throughout the WUS.Cui et al (2023) show no significant trend in snow RF over middle North America (which includes WUS), but there is recent evidence of increased LAP deposition in some basins (Clow et al 2016).Whether and where trends in snow RF exist at more local or watershed scales throughout the full WUS domain is unknown.
Remotely sensed data from the Moderate Resolution Imaging Spectroradiometer (MODIS) enables detection and quantification of RF due to LAPs in snow over vast geographic areas.MODIS is ideal for snow RF characterization, with a long-standing data record (Terra launched in 1999), frequent temporal resolution (1 d over the WUS), and moderate spatial resolution (463 m at nadir).Advancements to snow remote sensing and analysis provide the opportunity to characterize spatio-temporal variations in snow RF.Algorithms under continual refinement have been developed which allow detection of LAP presence and estimates of RF impacts to snow, such as the MODIS Snow Covered Area and Grain Size (MODSCAG) and MODIS Dust Radiative Forcing on Snow (MODDRFS) algorithms employed here (Painter et al 2009(Painter et al , 2012)).MODSCAG and MODDRFS data are further refined by using a time series of the data to account for clouds and off-nadir viewing issues, producing a spatially and temporally complete (STC) daily time series (Rittger et al 2020).Cui et al (2023) also utilized MODIS to assess spatio-temporal patterns in snow RF, but at a much coarser spatial resolution (0.05 deg.) and based on a different snow albedo product.Their study showed little to no snow remaining in the WUS by April, however it is known that the critical spring snowmelt season in a typical year starts between March and May and can last through June or later (e.g.Trujillo and Molotch 2014).MODIS snow detection methods typically use a band ratio but this approach has degraded performance in the late melt season relative to spectral mixture analysis that is employed by MODSCAG (Rittger et al 2013).Hence, patterns and trends in snow RF during the spring snowmelt season have yet to be established with MODIS over this domain.
In this study, we use remotely sensed data (i.e.fractional snow covered area (fSCA) from STC-MODSCAG and snow RF from STC-MODDRFS) to characterize spatio-temporal patterns of LAP presence at the snow surface from 2001 to 2022 during the annual snowmelt season (March-June) in the WUS.At three distinct spatial scales (regional, large river basin, and MODIS pixel scales), we investigate (1) prevailing geographic patterns of mean snow RF across multiple spatial scales, (2) the interannual variability of snow RF at the watershed scale, and (3) the locations and spatial scales at which trends emerge in the WUS.We focus on characterizing snow RF, as MODIS data are unable to distinguish the mechanism of LAP deposition (i.e.wet or dry events) or LAP type.The latter capability is possible with other recently launched hyperspectral missions (e.g.EnMAP launched in 2022, Guanter et al 2015), however these new missions lack the record length to establish long-term variations and trends in drivers of snow RF.A strength of the MODIS data is that they can document how snow RF has varied within and between regions and over recent decades, which will inform hydroclimate research and allow watershed management to better understand the extent to which current snowmelt is affected by impurities.In section 2, we describe the study areas, remote sensing data, and analysis approach.Results are highlighted in section 3 and discussed and interpreted in section 4, and we make conclusions in section 5.

Study area
Our analysis spanned the WUS, with the eastern boundary in Eastern Colorado being the farthest extent for which STC-MODSCAG/MODDRFS data are currently processed over the contiguous United States.This domain comprises a wide variety of climates, including coastal rainforest, Mediterranean, midlatitude desert, semiarid steppe, alpine, and humid continental.Variable snow conditions are associated with these environments, with boreal forest, montane forest, tundra, maritime, prairie, and ephemeral snow classes being found (Sturm and Liston 2021).
US watersheds are delineated by Hydrologic Unit Code (HUC), a hierarchical classification system created by the United States Geological Survey and based on surface hydrologic features in a standard geographic framework.This system divides the contiguous US into 21 regions (2-digit or 'HUC2') and 222 subregions (4-digit, 'HUC4' , here also termed 'watersheds') (Seaber et al 1987).This study focuses on 9 regions (HUC2 level) and 98 nested subregions (HUC4 level) in the WUS (figure 1), areas where the STC-MODSCAG/MODDRFS data are available.Snow RF was assessed at three spatial scales: per region (HUC2), sub-region/watershed (HUC4), and MODIS pixel levels (463 m).HUC2s average ∼4.6 × 10 5 km 2 in size and give boundaries of major geographic areas, whereas HUC4s average ∼2.7 × 10 4 km 2 in size and delineate large river basins (Seaber et al 1987).
Running the analysis at these three scales allowed us to investigate and contextualize the spatiotemporal variations of snow RF, and at what spatial scales they emerge from broad regional scales down to local levels.In addition, pixel trends were analyzed within 12 Level III Western Mountain Ecoregions (Omernik and Griffith 2014), due to the expected prevalence of LAPs in higher elevation mountain snowpack.

MODIS remotely sensed snow cover and RF
This study used daily STC snow RF gridded at 463 m from the MODSCAG (Dozier et al 2008, Painter et al 2009, Rittger et al 2020) and Dust Radiative Forcing of Snow (MODDRFS; Painter et al 2012) products derived from surface reflectance data collected by NASA's MODIS onboard the Terra satellite.MODIS collects daily information of Earth's surface reflectance at 36 non-continuous spectral bands ranging from 0.4 to 14.0 µm.
MODSCAG is created from a physically based algorithm which uses spectral mixture analysis to simultaneously identify subpixel fractional snow-covered area and estimates of grain size (Painter et al 2009).This model has been validated over the Rocky Mountains, high plains of Colorado, Sierra Nevadas, and Himalayas using Landsat 5, 7, and 8 fSCA (Painter et al 2009, Rittger et al 2013, 2021), fine-scale snow mapping with airborne lidar (Stillinger et al 2023), and in situ albedo measurements (Painter et al 2009).Using multiple endmember spectral mixture analysis, the MODSCAG algorithm leverages the unique spectral reflectance characteristics of various land surface types.Spectral reflectance retrievals from MODIS/Terra (MOD09GA), are used as input.Endmember reflectance values are combined by a series of simultaneous linear equations, where the relative contribution of each represents its fractional area coverage (Painter et al 2009).Both the fraction of each pixel that is covered by snow, as well as the associated grain size is simultaneously estimated, producing the fSCA data used later (see below).
The MODDRFS algorithm infers per-pixel snow RF by impurities at the snow surface including dust, BC, algae, and forest litter using MODIS surface reflectance data (Terra MOD09GA) and MODSCAG fSCA estimates (Painter et al 2012).This is done from a coupled radiative transfer model that distinguishes between clean snow and impurity-laden snow by inferring spectral differences according to snow grain size matching using the optical grain radius (OGR) in the near infrared and shortwave infrared, wavelengths which are less affected by impurity absorption.MODDRFS estimates the OGR and identifies the clean snow hemispherical-directional reflectance factor (HDRF-MOD09 retrieval), which varies with grain size, through the normalized difference grain size index (NDGSI).Estimates of instantaneous RF are produced using the difference in snow albedo (α) between clean snow from the NDGSI and observed snow in pixels with >95% fSCA.This is done by modeling α of snow both in the presence and absence of LAPs, then integrating the enhanced absorption of solar radiation over the range of snow reflectance (0.35-2.5 µm) (Skiles et al 2018).Note that we do not model the daily total RF, which would require additional data and assumptions about subdaily shortwave radiation variations (e.g.Cui et al 2023).
We utilize both MODSCAG and MODDRFS outputs to create STC products including the fSCA and RF analyzed herein.Our resulting STC data are minimally influenced by cloud cover, which frequently causes errors in snow identification, as it employs a snow-cloud identification method that utilizes a time series of fractional vegetation and rock land-surface data as well as persistence filters to flag snow-cloud identification errors and improve accuracy, reducing bias by 20% over other methods (Rittger et al 2020).MODIS data are impacted by off-nadir viewing, which STC accounts for with an interpolation weighting scheme based on the ground pixel size (Dozier et al 2008).Broadband α of clean snow is derived from the STC-MODSCAG grain size and solar illumination angle of each pixel.This offers a theoretical value for anticipated α under conditions absent of snow impurities (Painter et al 2009).Snow albedo is reduced in the visible part of the spectrum, with a 4%-6% RMSE and negligible bias relative to observations (Bair et al 2019).All data were produced and made available by Snow Today at the National Snow and Ice Data Center (NSIDC: Rittger and Raleigh 2022).

Analysis
We used the STC-MODSCAG/MODDRFS data (see section 2.2) to quantify snow RF while accounting for the snow covered area affected by the impurities.More specifically, this was done on a daily basis (1 March to 30 June each year) by multiplying the fSCA of each pixel area (m 2 ) by the instantaneous RF measured (W m −2 ), resulting in a metric that measures the instantaneous power (kW) due to snow impurities: (1) where: RF pixel (i) = instantaneous snow radiative forcing (kW) for the pixel on day i RF(i) = instantaneous snow radiative forcing for a unit area (W m −2 ) from MODDRFS on day i fSCA(i) = fractional snow covered area from MODSCAG on day i A = MODIS pixel area = 463 m * 463 m = 214 400 m 2 .This metric weighted the instantaneous RF according to fSCA in order to give a proportionally higher influence to pixels that have greater fSCA because both the amount of snow cover in a pixel and the instantaneous RF contribute to the actual impact of snow impurities on areal snowmelt.For instance, a pixel with low fSCA but high snow RF might have limited impacts due to the constrained amount of horizontal snow cover.Both fSCA and RF are dynamic in space and time, and this metric allows us to integrate both variables.Additionally, we note that fSCA and RF are not completely independent of each other in a single season.For instance, high RF acts to increase snowmelt rates which can decrease fSCA.New snowstorms can increase fSCA while decreasing RF when burying snow impurities.
In calculating the RF pixel metric, we masked out some pixels where uncertainty was greater; namely, pixels at lower elevations, with low fSCA, with significant forest cover, and/or that were glaciated.More specifically, only non-glaciated pixels above 800 m elevation that have <25% forest cover and >15% fSCA were included in the analysis.Lower elevations were filtered out due to prevalence of ephemeral snow in those zones.Recent validation with 3 m airborne observations showed high accuracy above 10% fSCA (Stillinger et al 2023).Previous studies have shown increasing uncertainty in fSCA with increasing forest cover (e.g.Raleigh et al 2013) with highest retrieval uncertainty somewhere between 66% and 81% forest cover (Rittger et al 2020).Our decision to mask out forests is further justified because the snow surface receives less shortwave radiation due to canopy shading in intact forests, thereby minimizing the impact of LAPs in subcanopy areas (Maurer and Bowling 2014).We used the Hansen et al (2013) global forest cover dataset, taking the maximum fractional tree cover at each pixel across the 2000 and 2015 datasets, masking out any burned forests as well as healthy forests.We employed viewable snow fraction data, without canopy corrections, to ensure analysis on the original observations.Viewable snow fraction is the fraction of a given area of the earth's horizontal surface that is covered by snow as seen from the satellite, not accounting for snow hidden by vegetation.Additionally, we masked out glaciated pixels using a glacier mask derived from Global Land Ice Measurements from Space data, as melt-water on the surface of glaciers can induce uncertainty to remote sensing of snow, and some glaciers may have been snow covered in some months of this study but not others (Kargel et al 2005, Mauro andFugazza 2022).
Seasonal summary statistics of RF pixel (equation ( 1)) were calculated and assessed at the three spatial scales of our analysis (pixel, HUC4, and HUC2) in each of the 22 study years.Using pixels that were not masked out (see above), we calculated a single average RF power value for the ablation season, which we defined as the 122 d spanning 1st March-30th June, yielding annual 22 values from 2001 to 2022.HUC2 and HUC4 masks were used to summarize statistics on long-term mean and spatial standard deviation in RF and fSCA from the STC-MODSCAG/MODDRFS data.Pixels along watershed boundaries were included when >50% of the pixel area fell within the HUC2 or HUC4 watershed.Although fewer years of observations are available than what is considered a typical climatology (30+ years), the 22 year dataset provides a substantial collection of information from which to derive contemporary mean conditions and assess year-to-year variations over two decades.
With a baseline established for average RF conditions per HUC4 and HUC2 unit area, we examined the interannual fluctuations since 2001 by calculating standard deviation in annual RF.We also generated individual year-by-year maps of anomalous RF conditions per HUC4, with each annual RF transformed with a z-score normalization based on the 22 year mean and standard deviation of each HUC4.In this interannual analysis, we assess the regional consistency of anomalies in annual RF and identify years with high and low anomalies in snow RF.
To detect interannual trends per HUC2, HUC4, and pixel, we used a Mann-Kendall significance test with year as the independent variable and mean RF as the dependent variable.This is a robust, nonparametric test used to detect monotonic trends within a time series (Mann 1945, Kendall 1975).Thiel-Sen's slope was used to determine the magnitude of the trend on pixels only, a method which calculates a median slope between data values.This suited the needs of this study due to the relatively short record, as Sen's slope is not easily influenced by outliers (Sen 1968).

Mean snow RF spatial patterns
Mean pixel-scale RF exhibited complex spatial variability across the WUS (figure 2).Generally, areas of higher elevation exhibited elevated snow RF (r = 0.34, p < 0.001).This elevation effect may be partially explained by an increase in solar radiation, both spatially (i.e. less atmosphere at higher elevations) and temporally, as higher elevation snow cover persists longer into spring and summer (r = 0.43, p < 0.001) when potential solar radiation is maximized.Mountain ranges exhibited greater RF including the San Juan and Parks Ranges of the Colorado Rockies, Wasatch and Uinta mountains in Utah, and the central-to-southern Sierra Nevadas in California.Additionally, the mountains of south-central Idaho, western Wyoming, and peaks of the Cascades showed significant RF but were less ubiquitous.
Spatially expansive coverage of RF values were located centrally in the WUS study area, notably in the lowlands of the Great Basin (Nevada, Utah, Wyoming, and southern Idaho).The plains to the east of the Rockies in Montana, Wyoming, and Nebraska also showed broad coverage of low RF values, but these were generally lower than the values seen in the Great Basin and intermountain zones.
At the HUC4 sub-basin level (figure 3(a)), the greatest values of snow RF were found in watersheds of the Colorado Rockies in the Upper Colorado region, where multiple HUC4s had mean RF values in the range of 120-132 kW.Notable RF intensity ranging from 105 to 116 kW was also seen in the following sub-basins: San Joaquin (105 kW) of the Central Sierra Nevadas in California, Willamette (116 kW) in the Oregon Cascades, Black Rock Desert-Humboldt (105 kW) in northern Nevada, Upper Snake (105 kW) in Idaho and Wyoming, and Escalante Desert-Sevier Lake (110 kW) of Utah.A mix of average RF values ranging from 8 to 100 kW were found in watersheds across the remaining WUS.The lowest RF values were found in the eastern drainages from the Rocky Mountains, with similar values in watersheds across all latitudes (3-15 kW).Of the 98 HUC4s in the WUS, we analyzed a total of 77.The other 21 HUC4s were removed from the analysis because they were biased by spurious RF hotspots (e.g.salt flats) or did not have sufficient snow cover to be included (e.g.southwestern California, southern Arizona and New Mexico, and on the eastern edge of the study area; plotted as white in figure 3).
Within sub-basin variability in RF was high in many basins (figure 3(b)) but was highest in watersheds of the Upper Colorado and the Great Basin (145-180 kW).Watersheds of New Mexico and southern Colorado (Rio Grande region) showed substantial variability across space.Watersheds west of the Cascade crest in the Pacific Northwest (PNW) showed moderate variability (70-100 kW), as well as watersheds in northern Montana and the western Dakotas.Lower spatial variability (0-30 kW) was seen in the watersheds to the east of the Rockies.At the HUC2 level, long-term mean snow RF exhibited distinct regional differences (figure 4).The Upper Colorado region exhibited the greatest median RF (119 kW), followed by the PNW and Great Basin (both 89.7 kW), California (65.7 kW), and the Rio Grande (62.5 kW).Lower RF medians were seen in the Missouri (34.0 kW), Arkansas-White-Red (27.8 kW), and Lower Colorado (19.3 kW).The Texas-Gulf region (not shown in figure 4) had isolated and high RF values that were deemed spurious (e.g.dry lake beds misidentified as snow), and was therefore excluded from further analysis.

Interannual variability
We summarized the HUC2-scale interannual variability in average snow RF as shown in the summarized boxplot distributions in figure 4. Based on the interquartile range (IQR), the greatest interannual variability in RF was found in the Rio Grande region (25 kW), followed by the Great Basin (22 kW), then the Upper Colorado (17 kW).Moderate interannual variability was found in the PNW (14 kW).The IQR of the remaining regions displayed lower temporal variability (8-11 kW), with Arkansas-White-Red, Missouri, and Lower Colorado showing the lowest values.Generally, the spread of values was commensurate with the overall RF measurements taken from each region (greater RF values lended to greater year-to-year variability).
At the HUC4 scale, the greatest interannual variability was found in the Lower Colorado (35-45 kW), Upper-Colorado (22-27 kW), Great Basin (11-23 kW), and Eastern PNW (17-21 kW) areas (figure 5).The California region included two HUC4s that exhibited high interannual variability: the North Lahontan just  Examining annual RF at the HUC4 level (figure 6), we found that anomalous high or low snow RF conditions tended to be grouped together regionally, although much complexity existed in both space and time.Notable high RF years across the full domain included 2008 while lower RF years included 2015 and 2019.Most years indicated that over the entire WUS, localized areas of both above and below average conditions existed in different watersheds within the ablation season.One year diverged from this pattern; 2020 was the year when RF impacts were most spatially limited in the 22 year record, as 85% of analyzed HUCs displayed below average RF values during that year.

Snow RF trends
A relatively small area (0.03%-8%) of analyzed snow exhibited a significant snow RF trend at the pixel level, and for these areas the trends tended to be decreasing (figure 7).There were some coherent areas with statistically significant (p < 0.05) RF trends (both increasing and decreasing), with much complexity across the WUS study area (figure 7(a)).Over the WUS, most pixel-scale snow RF trends ranged between −0.5 and +0.4 MW yr −1 , with outlying min and max values ranging from −1.1 to +0.9 MW yr −1 .
We examined the percentage of analyzed pixels with either a significantly decreasing or significantly increasing trend per mountain ecoregion to characterize the unique pattern and distribution of RF trends.We detected trends in each ecoregion (figure 7(b)).In general, decreasing trends were more common than increasing trends, with the North Cascades showing the greatest percent decline in RF (8%), followed by the Canadian Rockies (7%) and Sierras (6.5%) (figures 8(a) and (c)).The Southern Rockies (figure 8(d)), usually exhibited declining trends where trends were significant, except for a minor cluster of increasing trends observed towards the eastern side of the San Juan mountains.The Arizona/New Mexico mountains showed very few pixels with a trend (0.03% pixels decreasing).The greatest increase in snow RF was found in the Idaho Batholith ecoregion (4.5% pixels increasing) (figure 8(b)).The Idaho Batholith region was also unique in that it had a balance in the land area with increasing and decreasing RF trends (figure 7(b)).The two most northern mountain ecoregions-the Canadian Rockies and North Cascades-showed the next two most significant increases in RF trend by area (∼2.3%).All other ecoregions rarely exhibited (<1%) increasing trends in RF.
Because ecoregions are variable in overall size and analyzable snow cover (e.g.due to differences in forest cover, elevation distributions, see section 2.3), these percentages did not represent the absolute land area with a trend.For example, when comparing the North Cascades, (8% decline) and Sierra Nevada (6.5% decline), we found that a greater number of pixels had a significantly declining trend in the Sierras, although this was the third highest percentage by mountain ecoregion.We also noted elevated decreasing trends in RF around the periphery of forests.This may be an artifact of the data analysis (e.g.some forests being included in MODIS pixels at off-nadir views, Rittger et al 2020).To qualitatively investigate the possible connections between trends in SCA and trends in RF, at both the watershed and pixel scales we visually compared trend maps of each variable.We found that spatial and temporal patterns of snow cover did not closely align with that of RF trends and patterns (not shown).
At the HUC4 scale, we found that 7 of the 77 analyzed watersheds demonstrated a statistically significant (p < 0.05) trend in RF (no figures shown).The HUC4s with a decreasing trend were: Northern-Mojave Mono Lake in southeastern California, Bear in northern Utah and southeastern Idaho, Upper Colorado-Dirty Devil of southern Utah, White-Yampa of northern Colorado, Lower Canadian, Upper Canadian, and North Canadian.The latter three were located in northern New Mexico and Texas, on the plains to the east of the Sangre de Cristo Mountains and San Juan Mountains.This cluster is within the Arkansas-White-Red HUC2.
The RF trend analysis at the HUC2 scale revealed that only one region had a significant (p < 0.05) trend in RF (no figures shown).Specifically, the Arkansas-White-Red had a slightly decreasing trend in snow RF.While this HUC2 includes a vast area of land (642 000 km 2 ), most observable snow RF measurements were contributed from confined areas at the eastern foothills and plains of the Colorado Rockies and New Mexico Sangre de Cristo Mountains east of the San Juans.

Spatiotemporal patterns by region
This study confirms the expectation that snow RF in the Upper Colorado basin is significantly higher than other regions in the WUS (figures 2-4).These patterns align with what is reported in the literature, as multiple studies occurring in the southwestern Colorado Rockies have found significant dust impacts on snow (Painter et al 2007, 2012, 2018, Skiles et al 2012, 2015) as well as BC impacts on snow in the Rocky Mountains (e.g.Gleason et al 2022).We find elevated RF in the San Juan and Parks ranges of Colorado, the first high-altitude contact area for predominantly southwesterly winds transporting dust from the southern Colorado Plateau (Skiles et al 2015).Major dust-on-snow impacts have also been reported in western areas of the Upper Colorado basin, such as in the Utah Wasatch mountains (Lang et al 2023), which are also shown in the analyzed MODIS data (figure 2).
This study also highlights notable RF in the Great Basin (tied for second highest median pixel RF by HUC2), which has seen less attention in the literature than the Upper Colorado Basin.The Great Basin is an arid environment with minimal forest cover where much of the area is covered by snow in the spring season.
Here, mountain snowpack tends to retain snow while lower elevations likely sublimate or melt faster, Cascades; greatest percent decrease in RF (8% of pixels analyzed).(b) Idaho Batholith; greatest percent increase in RF (4.5% of pixels analyzed).(c) Sierra Nevada; most coherent decline in RF over space, and the 3rd highest percentage of pixels with a declining trend (6.5%).(d) Southern Rockies; nearly unanimous decline in snow RF (4.5% of analyzed pixels exhibit a decrease, and <0.05% exhibit an increasing trend).
resulting in plentiful snow next to snow-free soil.It is possible that dust and/or BC is accumulating consistently over this large area, making the effects of RF more subtle while still having a significant influence on driving seasonal melt.
Similar to the Great Basin, the PNW had a high mean RF (particularly in the central Idaho mountains and in Washington and Oregon at high elevations along the western Cascades), but in contrast there was less analyzed area and it was more discontinuous spatially (figure 2).Part of this result is due to our masking of pixels with >25% canopy cover and elevations <800 m (section 2.3).For the Cascades and Olympic ranges in particular, this constrained the analysis to just the higher peaks below glaciated areas (figure 2).A potential contributing factor to the high RF values in the high Cascades may be that these areas can accumulate a deep and dense snowpack in winter, which may persist late into summer when insolation values are greatest.Additionally, previous studies have reported BC and mineral dust in high elevation areas in the Cascade and Olympic Ranges (Doherty et al 2014, Kaspari et al 2015), as well as across Idaho (Doherty et al 2016) where we see moderate-to-high snow RF (figure 3(a)) and the most notable increasing trends in snow RF (figure 8(b)).Further inland locations in the PNW (e.g.Wallowas and Idaho Mountains) are located in more arid regions adjacent to agricultural and range lands, and recorded high RF values in the MODIS data (figure 2).
Relatively moderate snow RF was found in the Rio Grande and California regions.Much of the Rio Grande region does not experience consistent snow cover, but the northern portion extends up into the San Juan mountains and San Luis Valley of Colorado, which highly impact the values reported here.In California, the primary snow RF values were measured in the High Sierras, with some secondary contribution from the southern Cascades.Previous studies have shown high impacts on snow albedo from BC and pumice deposition in the Eastern Sierra (Sterle et al 2013).
The remaining regions (Arkansas-White-Red, Missouri, and Lower Colorado) had the lowest relative snow RF in the WUS.The Arkansas-White-Red region has comparatively little snow cover on average; the western corner includes the eastern edge of the Colorado and New Mexico Rockies, responsible for the RF values here.This is where we find multiple HUC4 watersheds with statistically significant declining trends in RF over the study period.The Missouri region is characterized by a large area of analyzable snow, but very low values of RF.The Lower Colorado receives little persisting snow on average and the second lowest mean pixel RF.

Regional characteristics and possible drivers of snow RF
The primary scope of the current study is to establish geographic and annual patterns of snow RF.At the pixel scale, the spatial patterns in snow RF are linked to elevation (figure 2) and locations where snow persists longer into summer when solar radiation is more intense, in agreement with other studies (e.g.Gleason et al 2022, Réveillet et al 2022).However, as noted by Réveillet et al (2022), these elevation-dependent patterns are superimposed on more complex regional patterns in LAP source (e.g.dust vs. BC) and deposition dynamics, as well as solar radiation regime and snow persistence patterns due to spatiotemporal variations in snowfall accumulation.Additional work is needed for more direct attribution to explain drivers such as climate and land use, and to more directly connect snow RF patterns to LAP sources and deposition events (e.g.back trajectory analyses, Skiles et al 2015).
To interpret our results and shed light on possible drivers of snow RF at the HUC4 level, we assess a hypothesis that the observed spatial and temporal variations in RF are linked to the type of impurities present in a watershed as well as climate (snowfall and solar radiation).We first classify the HUC4 basins using kmeans clustering based on their long-term mean (figure 3(a)), within-basin spatial standard deviation (figure 3(b)), and annual standard deviation (figure 5) in snow RF.We then examine the characteristics of these five clusters, namely: (1) mean annual snowfall derived from PRISM (Daly et al 1994; assuming daily precipitation is snowfall when daily temperature <0 • C), (2) March-June solar radiation from PRISM (Norm91m product, Rupp et al 2022), and surface concentrations of (3) dust and (4) BC over January-June from the MERRA-2 aerosol reanalysis (M2T1NXAER v. 5.12.4,Buchard et al 2017, Randles et al 2017).All variables are the average for each HUC4.We then rescale all four variables with the z-score normalization to focus on relative spatial patterns of each across the WUS domain.Because the dust and BC data each spanned multiple orders of magnitude, we also apply a log base 10 transformation before applying the z-score.Spatial maps of these four variables show higher snowfall in the Rockies and Northwest, higher solar and dust toward the southwest, and higher BC along the west coast and near population centers (figure 9).
Figures 10(a)-(c) shows the clustering analysis.Snow RF is highest and most variable in both space and time in cluster 5, which includes the Upper Colorado, northern Great Basin, and southeastern PNW regions (Idaho).Cluster 5 is also characterized by relatively high snowfall (figure 10(d)) and high solar radiation (figure 10(e)), with greater dust than BC relative to the rest of the WUS (figures 10(f) and (g)).This suggests that watersheds which experience greater variability of snow RF in both space and time may be associated with dust over other LAP types.The high interannual variability in snow RF in cluster 5 (figure 10(b)) is consistent with the episodic nature of dust storms, which vary in number and severity from year to year (e.g.Skiles et al 2015), depending on antecedent winter moisture conditions, surface crusting, and vegetative effects (Brazel and Nickling 1986).Generally, greater spatial variability (within-watershed) is exhibited in watersheds located within or downwind of more arid environments (figure 10), which are known source regions for dust.
According to the MERRA-2 aerosol reanalysis, the highest surface concentrations of dust and BC do not coincide with the areas with highest snow RF (figure 10), which indicates that other factors (e.g.snowfall and snow duration, climate, and LAP transport) are important.For instance, the highest relative dust concentrations are found in the arid southwest (e.g.cluster 2; figure 10), which has experienced recent prolonged droughts (2002( -2003( , 2007( -2009( , and 2012( -2016( ) (MacDonald 2010)).However, these areas with relatively high dust concentration have lower snow RF, likely due to the low snowfall in these areas (figures 10(d) and (f)), which inherently limits the snow season length and the amount of additional radiation that may be absorbed by impurities.In contrast, the highest snow RF area (cluster 5) has moderate dust concentrations (in relative terms), which when coupled with high snowfall and high solar radiation may explain the high snow RF (figures 10(d)-(f)).Note that there may be limitations in the source data (MERRA-2 reanalysis), such as the spatial resolution.More work is needed to validate its accuracy in these regions.
Moderate snow RF is found in the other clusters, but these vary in terms of their spatial and temporal variability.For instance, cluster 4 (e.g.western PNW, southern Sierra) has lower spatial variability than cluster 3 (figure 10(a)) and lower temporal variability than cluster 5 (figure 10(b)).Cluster 4 is associated with high snowfall (figure 10(d)) and the highest BC concentration across the WUS (figure 10(g)), but the lowest dust concentrations (figure 10(f)) and lowest solar radiation (figure 10(e)) due to clouds and latitude.This suggests that areas in the WUS with lower year-to-year variability and lower spatial variability of RF may be associated with anthropogenic sources like BC.With the exception of unusual societal disruptions such as the Covid-19 pandemic and its possible impacts on snow impurities in 2020 (figure 6 and Bair et al 2021), anthropogenic activity is fairly consistent from year-to-year, so associated impurity sources might be emitted with more temporal consistency compared to dust entrainment in arid environments.Long-range, anthropogenic impurity sources are known to have a pattern of consistent, frequent transport (VanCuren and Cahill 2002).Pacific Coast adjacent watersheds including mountainous areas (Cascades and Southern Sierra Nevada) may possibly be influenced by trans-Pacific particulate movement, BC emissions from nearby large urban centers along the west coast, or dust entrainment from extensive farming activities in the Central and Willamette Valleys, positioned upwind of the mountain ranges (Hadley et al 2007, Vicars andSickman 2011).
The four factors shown in figures 9 and 10 (i.e.snowfall, solar radiation, dust, and BC) are significant predictors of mean snow RF (figure 3(a)).A multiple linear regression model developed with these four predictors explained 50% of the variance in the long-term mean snow RF across HUC4s (see supplementary document).The most impactful predictor of snow RF was total snowfall, followed by solar, then dust, and finally BC.However, a caveat is that there is moderate collinearity between these predictors.For instance, the spatial patterns of solar radiation and dust concentration have similarity (see figures 9(b) and (c)).More mechanistic studies are needed to refine understanding of the drivers in snow RF.

Limitations & outstanding research needs
The study was limited in several respects, and future research should address outstanding issues not addressed here.First, we did not analyze all landscapes, including post-wildfire areas and glaciers, focusing on mid-to-high elevation open areas over the WUS.Post-wildfire areas have been shown to have high snow RF, and the ecoregions of the WUS have experienced high rates of wildfires (e.g.Gleason et al 2019, Kampf et al 2022).While a suite of vegetation and soil are included in MODSCAG, burned vegetation is not.We did not analyze these areas so as to avoid conflating changes in snow RF with changes in land cover.However, the fSCA, snow albedo, and snow duration data have separately been analyzed for a few small fires in California which identified snow susceptibility from melt with both decreased snow albedo and canopy cover (shading) resulting in 50% less snow cover in 2022 and 50 fewer snow cover days compared to 2013, a year with similar snow accumulation (Hatchett et al 2023).In addition to forests, we also masked out glaciated areas due to challenges with discerning RF by particulates over bare ice with MODIS retrievals.However, glaciers are also impacted by impurities and absorb solar radiation through all of summer, and play an important role in cryosphere-hydrosphere interactions.Impacts from wildfire areas and glacial zones are therefore not quantified here, and both of these factors necessitate additional study.Future work should also examine finer-scale (e.g.watershed) patterns and trends in snow RF beyond the WUS, as increased RF by contaminants is a known issue with mountain snowpack globally (Yasunari et al 2010, Malmros et al 2018, Sarangi et al 2019, 2020).At sub-continental scales, Cui et al (2023) found significant decreasing trends in snow RF over a 16 year period in two of six large study areas, but more local assessments are needed to inform hydroclimate studies.
Another limitation is that our summarized data, while suitable for this study, is somewhat temporally coarse, aggregated from 1st March through 30th June to produce a single estimate of snow RF per season.We do not assess when changes occur within each melt season, which may provide insights into the nature of specific contamination events (e.g.dust storms, wet vs. dry deposition) as well as evidence of sources.Future work may involve investigating RF development in greater temporal detail within the snowmelt season, for example identifying and analyzing events relative to the background.
Future research should investigate the interplay between the timing and number of spring snowstorms, which we did not assess here.Spring snowfall impacts snow RF with more frequent storms refreshing snow albedo and lowering RF; this acts as a limiting factor of total RF that can be produced.LAP deposition events may occur during or between storms, and thus LAP particles may be buried in the snow and have a delayed impact, which amplifies as they gradually concentrate at the surface during the ablation season (Sterle et al 2013).More persistent snow (e.g.due to a higher snow accumulation in winter and persisting into the late spring-early summer) has a higher likelihood of surface accumulation of impurities, while also being present under greater insolation values thereby increasing the snow RF (e.g.Seidel et al 2016).Past research has high snow accumulation areas may have the highest sensitivity to RF impacts (Deems et al 2013), which aligns with our interpretation of our results (section 4.2).The interplay between snow RF and snowfall impacts each region uniquely, as some mountain ranges tend to have more late season storms.Interannual variations in snow accumulation are thus partially responsible for the patterns of year-to-year variability in snow RF (figure 6), but direct anthropogenic influences play a role as well.For instance, 2020 had the most widespread, below-average RF conditions across the study area, where nearly 85% of watersheds displayed below average conditions (figure 6).Other studies have found more dramatic reductions in snow RF in the Indus and Himalaya regions during this same time frame, directly linked to decreased pollution from the COVID-19 lockdowns (Bair et al 2021).
More research is also needed to understand the reasons for divergent trends in snow RF and the connectivity of specific LAP sources and the snow RF patterns, which could be achieved with techniques such as back trajectory calculations.One complexity is that our annual snow RF metric is influenced by multiple factors including the snow darkening due to LAPs and the persistence of seasonal snow cover.The latter is impacted by both winter/spring snow accumulation and is also impacted by LAPs.A decreasing trend in snow RF may signify a reduction in LAPs and/or a reduction in snow persistence/duration (e.g.due to less winter snowfall).As an example of the first effect (changes in LAPs), in California, updated emission standards and restrictions on diesel engines and biomass burning have significantly reduced atmospheric BC concentrations since the 1960s (CARB 2015, Kirchstetter et al 2017).This aligns with the tendency for declining snow RF (6.5% of the pixels) across the Sierra Nevada mountain range reported here, however other sources of BC (e.g.trans-Pacific) need to also be accounted for to better understand the declining trend in the Sierras.Past research has reported that both BC and dust affect mountain snowpack across Idaho (Doherty et al 2016) where we see the most notable increasing trends.To the second effect (changes in snow persistence/duration), large interannual variations in snow accumulation (e.g.Siler et al 2019) and snow cover make it challenging to detect trends in snow persistence (Bormann et al 2018), which can impact the snow RF trends.
Finally, we note that there is a need for longer-term monitoring of snow RF.The MODIS record does not have enough years of data available to meet the requirements of a more conventional climatology (e.g.30+ years ideal for climate purposes).Some regions show quite large interannual variability (figures 6 and 7) which can mask possible trends in a shorter record.In order to gain more robust insights into the spatio-temporal patterns examined here, we suggest future work extend the current analysis as more years of data become available.MODIS is nearing its end of operation, but the Visible Infrared Imaging Radiometer Suite onboard the Suomi NPP and NOAA-20 satellites can map similar snow properties (Rittger et al 2021) and holds the potential to continue the record going forward.

Conclusions and outlook
This study offers a first-of-its-kind documentation and characterization of the spatio-temporal patterns and trends of snow RF by LAPs across multiple spatial scales in the WUS for the snowmelt season (March-June) over the period of 2001-2022.This information enables us to better understand patterns of snow albedo, radiative balance, and melt processes, all climate-sensitive variables that are critical to monitor for water management and regional hydroclimate.Future work should pursue more direct, mechanistic attribution studies (e.g. with back trajectory analyses) to investigate LAP sources and to more fully contextualize the patterns and trends in snow RF described here.
At three distinct spatial scales, we use remote sensing data spanning 2001-2022 to establish baseline conditions, examine interannual variability, and conduct trend analyses of RF on snow by LAPs.We confirm that the Upper Colorado region has the highest consistent values of snow RF in the WUS.However, other regions (e.g.Great Basin, PNW) display relatively high snow RF and have not been the subject of as many RF studies.
The temporal patterns of snow RF are largely stable over the last two decades for the mountain snowpack of the WUS, with over 90% of the analyzed area having no significant trends in snow RF since 2001.Where there are trends, these tend to be declining RF trends and for a small amount of mountainous area (up to 8%).Areas with declining trends are found in all mountain ecoregions, such as in the North Cascades, Sierra Nevadas, and throughout the Rocky Mountains.Some small and coherent areas with increasing trends were found in central Idaho.The Arkansas-White-Red was the only HUC2 region with a significant trend in snow RF (decreasing).Seven HUC4 watersheds exhibited a significant decreasing trend in RF.Whether these decreasing trends in snow RF are due to decreases in snowfall and snow cover duration (e.g.climate change) and/or decreases in LAPs in snow is an open question that merits more detailed analysis in future work.
Increasing aridity as a result of climate change leads to increased dust emissions (Skiles et al 2018).Both rising temperatures and aridification will work in tandem to decrease snow persistence and water availability globally, although the relative contributions of each are yet to be fully understood.Having less available water to the environment in the dry months contributes to arid, dusty conditions, hence creating a positive feedback loop of early snow disappearance and less water in a warming climate.The predicted effects of climate change will continue to shift environmental conditions in the WUS in complex ways, likely altering snow accumulation and melt dynamics but also aridity, wind erosion, and consequent dust generation while regional populations and industries continue to grow.This may have policy implications, considering that population and anthropogenic activity of this region have steadily grown over the past decades.With innovative emissions and agricultural management policies, population growth and development can take place sustainably and dryland erosion can be minimized.In turn, this may decrease impacts on the scarce and overallocated water resources of the west.
Our results help highlight potential hotspots for LAP loading to mountain snowpack that may lead to drought conditions and water scarcity.Future research should be conducted in these LAP hotspots to more directly characterize the prevalent LAP types and source regions, thereby allowing the implementation of effective and directed mitigation strategies.Our results suggest that LAP abatement efforts may be an important component of climate adaptation strategies in order to reduce scenarios of drought and water scarcity.When impurities at the snow surface are not considered, operational hydrologic models will underestimate melt rates (Skiles et al 2012), adversely affecting the reliability of streamflow forecasting (Bryant et al 2013).This work informs where and when RF impacts on snowmelt may need to be considered in hydrologic models.Additionally, it provides essential context for understanding patterns and trends in snowmelt, which influence, and are influenced by climate change (Deems et al 2013).

Figure 2 .
Figure 2. Mean RF by snow impurities over the Western US, 2001-2022.Note that pixels with forest cover, glaciers, insufficient snow cover, or at low elevations were not analyzed and masked out (see section 2.3).Note also that the dataset analyzed did not include values in areas outside the US (e.g.parts of Pacific Northwest HUC2 in Canada).

Figure 3 .
Figure 3. (a) Average RF of snow per HUC4 sub-basin over Western US from 2001 to 2022.(b) Snow RF variability within-basin, as measured by the spatial standard deviation averaged over the snowmelt season.

Figure 4 .
Figure 4. Boxplot summarizing mean annual (2001-2022) snow RF averaged over each HUC2 region, organized in descending order by median.Wider boxes represent greater interannual variability while more constrained boxes represent more consistent RF conditions from year-to-year.The color-coded inset map shows the locations of the HUC2 regions.

Figure 5 .
Figure 5. Temporal standard deviation of ablation season RF values per WUS sub-basin.Note that some basins had values exceeding the highest colorbar value.

Figure 6 .
Figure 6.Annual RF anomalies from the 2001 to 2022 mean at the HUC4 scale.

Figure 7 .
Figure 7. Pixels exhibiting a significant (p < 0.05) trend in snow RF conditions, per ecoregion.(a) Spatial map of the significant trends, and (b) percentage of each mountain ecoregion with increasing or decreasing trends in snow RF.

Figure 8 .
Figure 8. Ecoregion trends in snow RF (MW yr −1 ) per pixel from 2001-2022, zoomed to show notable results.(a) NorthCascades; greatest percent decrease in RF (8% of pixels analyzed).(b) Idaho Batholith; greatest percent increase in RF (4.5% of pixels analyzed).(c) Sierra Nevada; most coherent decline in RF over space, and the 3rd highest percentage of pixels with a declining trend (6.5%).(d) Southern Rockies; nearly unanimous decline in snow RF (4.5% of analyzed pixels exhibit a decrease, and <0.05% exhibit an increasing trend).

9.
Spatial maps of long-term (a) annual snowfall, (b) solar radiation over March-June, (c) dust surface concentration over January-June, and (d) BC surface concentration over January-June.All characteristics are summarized at the HUC4 level and transformed into z-scores to emphasize differences from the westwide distribution (see main text).Snowfall and radiation are from PRISM climate data while the dust and BC concentration data are from the MERRA-2 aerosol reanalysis.

Figure 10 .
Figure 10.Comparisons of long-term mean snow RF versus (a) spatial standard deviation in RF and (b) annual standard deviation in RF across the HUC4 watersheds.Watersheds are classified into five groups based on kmeans clustering (with plus signs indicating cluster centroid), with (c) of the clusters mapped.The cluster characteristics are summarized with the boxplots in terms of z-score normalized (d) annual snowfall, (e) spring solar radiation, (f) dust surface concentration, and (g) BC surface concentration.See figure 9 for the spatial distributions of these characteristics.