Snow Phenology and Hydrologic Timing in the Yukon River Basin, AK, USA

: The Yukon River basin encompasses over 832,000 km 2 of boreal Arctic Alaska and north-west Canada, providing a major transportation corridor and multiple natural resources to regional communities. The river seasonal hydrology is deﬁned by a long winter frozen season and a snowmelt-driven spring ﬂood pulse. Capabilities for accurate monitoring and forecasting of the annual spring freshet and river ice breakup (RIB) in the Yukon and other northern rivers is limited, but critical for understanding hydrologic processes related to snow, and for assessing ﬂood-related risks to regional communities. We developed a regional snow phenology record using satellite passive microwave remote sensing to elucidate interactions between the timing of upland snowmelt and the downstream spring ﬂood pulse and RIB in the Yukon. The seasonal snow metrics included annual Main Melt Onset Date (MMOD), Snowoff (SO) and Snowmelt Duration (SMD) derived from multifrequency (18.7 and 36.5 GHz) daily brightness temperatures and a physically-based Gradient Ratio Polarization (GRP) retrieval algorithm. The resulting snow phenology record extends over a 29-year period (1988–2016) with 6.25 km grid resolution. The MMOD retrievals showed good agreement with similar snow metrics derived from in situ weather station measurements of snowpack water equivalence (r = 0.48, bias = − 3.63 days) and surface air temperatures (r = 0.69, bias = 1 day). The MMOD and SO impact on the spring freshet was investigated by comparing areal quantiles of the remotely sensed snow metrics with measured streamﬂow quantiles over selected sub-basins. The SO 50% quantile showed the strongest (p < 0.1) correspondence with the measured spring ﬂood pulse at Stevens Village (r = 0.71) and Pilot (r = 0.63) river gaging stations, representing two major Yukon sub-basins. MMOD quantiles indicating 20% and 50% of a catchment under active snowmelt corresponded favorably with downstream RIB (r = 0.61) from 19 river observation stations spanning a range of Yukon sub-basins; these results also revealed a 14–27 day lag between MMOD and subsequent RIB. Together, the satellite based MMOD and SO metrics show potential value for regional monitoring and forecasting of the spring ﬂood pulse and RIB timing in the Yukon and other boreal Arctic basins.


Introduction
Annual hydrologic variability in snowmelt dominated basins is reflected in snow cover dynamics [1]. Several studies in the western United States found warming temperatures promoted an earlier snowmelt onset and subsequent runoff [2][3][4]. The earlier onset of snowmelt can slow the rate of melt due to less early season incoming solar radiation [5] while proportionally reducing annual streamflow [1]. However, seasonal trends toward earlier snowmelt and longer melt durations and the associated impacts on snowmelt driven flooding and river ice break up (RIB) patterns remain relatively uncertain [6,7].

Study Domain
The Yukon River traverses east to west and, along with its tributaries, constitutes one of North America's largest river basins. This region experiences six to nine months of snow cover annually, and spring snowmelt runoff is the main hydrologic contribution to the system [10]. The YRB has a mean annual discharge of 6400 m 3 s −1 , a drainage area exceeding 853,300 km 2 [8,29] and covers 10 degrees of latitude from 59 • N to 69 • N, extending into the Canadian Yukon and British Columbia to the east, and reaching the Alaska Bering Sea coast to the west. The diverse topography, with a median elevation of 617 m and extending from sea level to the highest elevations of the Brooks (2735 m) and Alaska (6190 m) Ranges, encompasses a diversity of northern boreal, arctic, alpine and maritime biomes. Evergreen needleleaf forests are the dominant vegetation cover (54%) followed by broadleaf deciduous forests (9%) covering the valley bottoms and into the mid-elevations. The Yukon Delta and higher elevations have tall and low shrubs (9%) mixed with some dry and wet herbaceous (9%) tundra as the dominant plant community. Permafrost is present to a large extent in the YRB, and comprises several types including sporadic (14%), discontinuous (46%) and continuous (16%) and moderately thick to thin permafrost (24%) [32].
We constrained our analysis of the YRB to three major catchments defined by the location of reliable long-term gauging stations located on main stem of the Yukon River. Catchment distributions within the YRB are shown in Figure 1, including Eagle (287,800 km 2 ), Stevens Village (500,968 km 2 ) and Pilot (817,961 km 2 ), in order of catchment size.

Passive Microwave Satellite Record 1988-2016
The detection of MMOD and SO used the 19 and 37 GHz afternoon Tb retrievals at horizontal (H) and vertical (V) polarizations from the MEaSUREs Calibrated Enhanced-Resolution Passive Microwave Daily EASE-Grid 2.0 Brightness Temperature ESDR, available at the National Snow and Ice Data Center [33].This multidecadal data record represents Tb retrieval records calibrated across multiple sensors and platforms for different frequencies and polarizations from the NOAA DMSP Special Sensor Microwave/Imager

Passive Microwave Satellite Record 1988-2016
The detection of MMOD and SO used the 19 and 37 GHz afternoon Tb retrievals at horizontal (H) and vertical (V) polarizations from the MEaSUREs Calibrated Enhanced-Resolution Passive Microwave Daily EASE-Grid 2.0 Brightness Temperature ESDR, available at the National Snow and Ice Data Center [33].This multidecadal data record represents Tb retrieval records calibrated across multiple sensors and platforms for different frequencies and polarizations from the NOAA DMSP Special Sensor Microwave/Imager (SSM/I) and the Special Sensor Microwave Imager/Sounder (SSMIS) [34]. Each platform has several sensors. When using SSM/I, we used DMSP-F8 (1987)(1988)(1989)(1990)(1991), DMSP-F10 (1992DMSP-F10 ( -1997, and DMSP-F13 (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005). For missing temporal observations, we gap-filled using overlapping sensors. The sampling resolution of the combined 19 and 37 GHz Tb retrievals are~25 km or coarser; however, the MEaSUREs products used for this study were processed using the scatterometer image reconstruction approach to obtain an enhanced spatial grid resolution of 6.25 km (19 GHz) and 3.125 km (37 GHz) from the overlapping Tb antenna patterns [34,35]. To establish a continuous record, 3.125 km 37 GHz retrievals were resampled to 6.25 km using a nearest neighbor interpolation, and missing grid cells were gap-filled using a temporal linear interpolation of adjacent Tb retrievals [22].

Other Ancillary Datasets
We used composited weekly FW time series from a regional AMSR (Advanced Microwave Scanning Radiometer) FW record [36] to derive a mean summer (JJA) fractional water (FW) map for the AMSR period of record (2002 to 2016). The resulting static FW map was used to define where water bodies persist across the landscape and areas where the PMW land retrievals may be influenced by surface water contamination. We used glacier outlines acquired from the Global Land Ice Measurements from Space (GLIMS) program to identify and mask out SO dates and SMD grids from statistical analysis for regions in Alaska that have snow and ice year-round. We also used the Fractional Forest (FF) cover (%) from the MODIS (Moderate Resolution Imaging Spectroradiometer) MOD44B land cover product [37] to assess the influence of vegetation cover on the relative accuracy of the snowmelt metrics. The 250 m resolution MOD44B data were reprojected to match the 6.25 km PMW grid using nearest-neighbor resampling.
Monthly gridded surface air temperatures were obtained for the YRB from the Weather Research and Forecasting (WRF) Reanalysis, downscaled to 20 km using ERA-Interim historical reanalysis data  acquired from the Scenarios Network for Alaska and Arctic Planning (SNAP) [38]. We used nearest-neighbor interpolation to resample monthly average temperatures to the 6.25 km PMW grid. Resampled annual monthly temperature grids were then used to produce correlations with the annual snowmelt metrics. A complete list of the primary datasets used in this study is found in Table 1.

Deriving MMOD
We used a Tb spectral gradient ratio polarization (GRP) [24,39] and a Tb differencing approach [40] to detect MMOD (henceforth we denote the PMW derived snow metrics as MMOD PMW , SO PMW , and SMD PMW ). The GRP algorithm detects the MMOD PMW within a given grid cell and annual (water year) time series when the GRP running mean (rm) is less than the difference between the water year's average winter (January-March) GRP (win) and an input melt parameter threshold, param ( Figure 2). Here, param is a dynamic variable that ranges from 0.2 to 0.6 and determines how low the GRP must be relative to win for the MMOD PMW to be detected. The duration of the rm and param are spatially dependent on an ancillary static surface FW map. In general, we found the GRP snow signal to be degraded for grid cells with higher FW cover. However, the GRP derived MMOD PMW performance was improved by assigning a higher param value for cells with higher water coverage below a 39% FW threshold, while the param value was set at 0.2 for higher FW levels above this threshold. For certain years and locations, the GRP algorithm was unable to detect the MMOD PMW and for these grid cells we used the Tb difference between the 19V and 37V channels for the MMOD PMW classification [22,40]. between m and TbDIFF exceeds THold for four or more days, this indicates the MMOD (Figure 2b). More information on this algorithm is outlined in [22]. ; four-day GRP running mean (rm), (orange); mean GRP value from 1 January to 1 March (win) (black); threshold value derived from subtracting the winter mean GRP (win) and FW melt threshold parameter (param) (horizontal blue line). (b) MMOD detection algorithm outlined in [22]. (c) Daily average air temperature and snow depth.
SOPMW was derived using the 19V-37V Tb differencing approach, which exploits varying sensitivity of the different Tb frequencies to surface scattering, wetness and dielectric properties between snow-covered and fully ablated landscapes [25]. The Tb difference is positive and relatively stable prior to the seasonal onset of snowmelt, but as the snow begins to ablate the Tb difference precipitously decreases, mirroring the snow ablation rate before reaching a seasonal minimum. The Tb difference reaches a seasonal minimum just after the snowpack has largely ablated, which we use to represent the SOPMW condition. SMDPMW is then derived by taking the difference between the SOPMW and MMODPMW, resulting in the number of days required for the snowpack to melt out for a given grid cell. ; four-day GRP running mean (rm), (orange); mean GRP value from 1 January to 1 March (win) (black); threshold value derived from subtracting the winter mean GRP (win) and FW melt threshold parameter (param) (horizontal blue line). (b) MMOD detection algorithm outlined in [22]. (c) Daily average air temperature and snow depth.

MMOD Validation and Evaluation
The Tb difference algorithm requires the derivation of three variables, including the daily difference between 19V-37V (Tb DIFF ), previous three-day average of Tb DIFF (m), and the product of m multiplied by an empirical constant, 0.35 (TH old ). When the difference between m and Tb DIFF exceeds TH old for four or more days, this indicates the MMOD (Figure 2b). More information on this algorithm is outlined in [22]. SO PMW was derived using the 19V-37V Tb differencing approach, which exploits varying sensitivity of the different Tb frequencies to surface scattering, wetness and dielectric properties between snow-covered and fully ablated landscapes [25]. The Tb difference is positive and relatively stable prior to the seasonal onset of snowmelt, but as the snow begins to ablate the Tb difference precipitously decreases, mirroring the snow ablation rate before reaching a seasonal minimum. The Tb difference reaches a seasonal minimum just after the snowpack has largely ablated, which we use to represent the SO PMW condition. SMD PMW is then derived by taking the difference between the SO PMW and MMOD PMW , resulting in the number of days required for the snowpack to melt out for a given grid cell.

MMOD Validation and Evaluation
The MMOD PMW metric was validated using daily Snow Water Equivalent (SWE) and air temperature measurements at 20 Alaska SNOw TELemetry (SNOTEL) network stations (Appendix A). At each SNOTEL station, we calculated the site-level MMOD using both SWE (MMOD SWE ) and temperature (MMOD T ) measurements and compared these measurements to the collocated satellite based MMOD PMW retrieval. Each local station MMOD was derived using an eight-day forward moving window temporal mean of the daily SWE measurements to extract the date where SWE was at its peak, with the assumption that depletion after this date corresponds with a melting snowpack ( Figure 2). Each alternative temperature derived MMOD T was calculated using an eight-day rolling mean of daily air temperature measurements at each station, where the MMOD was defined as the first date where the mean air temperature exceeded 0 • C. Both the SWE and air temperature derived MMOD definitions infer the shift in the seasonal energy budget that initializes snowmelt [41]. The relative accuracy of the MMOD PMW retrievals was assessed against the SNOTEL site-based MMOD observations using bias and correlation as measures of performance. The MMOD PMW record was also compared against a spatially continuous annual record of primary spring thaw timing derived from a daily Freeze/Thaw (FT) classification of SSM/I 37V GHz Tb retrievals spanning the same domain and multiyear record as the current study but derived at a coarser 25 km spatial resolution [23].

Snow and Streamflow Indices in the Yukon River Basin
The differences in timing between upland snowmelt and basin streamflow were examined across a selection of regional catchments within the larger YRB. Daily streamflow data were obtained from the Global Runoff Database Centre (GRDC) for three major river gaging stations located along the Yukon River main stem at Eagle, AK (64.79, −141.20), Stevens Village, AK (65.88, −149.72) and Pilot Station, AK (61.93, −162.88); the three stations are located at the outlets of the major YRB catchments, representing respective drainage areas of 287,800 km 2 , 500,969 km 2 , and 817,962 km 2 . Except for the years 1997-2002 at Pilot, all stations had complete streamflow data records spanning the study period .
Streamflow timing indices were represented as Q20, Q50, and Q80 terms calculated as the respective days of the year (DOY) when 20, 50 and 80% of the total annual (WY) flow passed the station gage location [2,3]. We also extracted the DOY representing the peak annual discharge at each station ( Figure 3). Snowmelt indices were also represented from the PMW record as quantiles but calculated as the DOY when 20, 50 and 80% of the cumulative catchment area [%] above each streamflow station was either under active snowmelt or snowoff conditions. We represented the cumulative area of snowmelt contributing to the measured discharge for each basin as the difference between the cumulative catchment area under respective melt and snowoff conditions (denoted as contributing in Figure 4). Discharge and snowmelt indices were extracted for each year and station location. Least-squares linear regression was used to quantify the relationship between the streamflow (dependent) and snowmelt (independent) indices to determine how much streamflow variability could be explained by the level of snowmelt activity. The resulting regressions were also used to identify any apparent lead time between the PMW derived catchment snow metrics and the subsequent downstream flood pulse indicated from the catchment discharge measurements. Remote Sens. 2021, 13, x 8 of 25

Snow and River Ice Breakup in the Yukon River Basin
We next examined the interaction between the timing of upland snowmelt and seasonal river ice breakup within the YRB. RIB observations were acquired for 19 station locations ( Figure 1) along the Yukon River main stem from the National Weather Service's Historical River Observations Database (www.weather.gov/aprfc/rivobs, Accessed: November 1, 2020). The number of annual RIB observations ranged from 18 to 28 years with a mean record length of 24 years from 1988 to 2016. All river ice observation locations used in this study are listed in Appendix B.
We again used linear regression analysis to quantify the relationship between annual RIB date and the snowmelt indices within each catchment. However, unlike the discharge analysis involving three major sub-basins, the river ice analysis encompassed a larger number of catchments associated with the more extensive RIB station network. Similar snow metric quantiles were calculated for the upstream catchments associated with each RIB station location. For snowmelt metrics that were statistically significant (p < 0.05), we then calculated the average annual difference (i.e. temporal lag) between the mean timing (DOY) of a given catchment snowmelt metric and the associated RIB date.

MMOD Classification Accuracy
The MMOD PMW results were generally consistent with the primary spring thaw onset indicated from the coarser (25 km resolution) FT-ESDR product [23]. In both datasets, MMOD/Spring Onset typically occurs later in the Alaska North Slope and at higher elevations ( Figure 5). The mean MMOD PMW in the YRB from 1988-2016 was DOY 113 (±11 days; temporal SD), while the FT-ESDR Spring Onset mean was DOY 108 (± 7 days) for the same period. Thus, MMOD PMW typically occurs about five days later than the FT-ESDR spring onset and has higher spatial (SD) heterogeneity attributed to the spatially enhanced Tb record used for the MMOD PMW retrievals. The MMOD PMW pattern also appears to better preserve the influence of the regional topography and land cover on spring melt timing, as described in the following sections.

Snow and River Ice Breakup in the Yukon River Basin
We next examined the interaction between the timing of upland snowmelt and seasonal river ice breakup within the YRB. RIB observations were acquired for 19 station locations ( Figure 1) along the Yukon River main stem from the National Weather Service's Historical River Observations Database (www.weather.gov/aprfc/rivobs, Accessed: November 1, 2020). The number of annual RIB observations ranged from 18 to 28 years with a mean record length of 24 years from 1988 to 2016. All river ice observation locations used in this study are listed in Appendix B.
We again used linear regression analysis to quantify the relationship between annual RIB date and the snowmelt indices within each catchment. However, unlike the discharge analysis involving three major sub-basins, the river ice analysis encompassed a larger number of catchments associated with the more extensive RIB station network. Similar snow metric quantiles were calculated for the upstream catchments associated with each RIB station location. For snowmelt metrics that were statistically significant (p < 0.05), we then calculated the average annual difference (i.e. temporal lag) between the mean timing (DOY) of a given catchment snowmelt metric and the associated RIB date.

MMOD Classification Accuracy
The MMODPMW results were generally consistent with the primary spring thaw onset indicated from the coarser (25 km resolution) FT-ESDR product [23]. In both datasets, MMOD/Spring Onset typically occurs later in the Alaska North Slope and at higher elevations ( Figure 5). The mean MMODPMW in the YRB from 1988-2016 was DOY 113 ( ± 11 days; temporal SD), while the FT-ESDR Spring Onset mean was DOY 108 ( ± 7 days) for the same period. Thus, MMODPMW typically occurs about five days later than the FT-ESDR spring onset and has higher spatial (SD) heterogeneity attributed to the spatially enhanced Tb record used for the MMODPMW retrievals. The MMODPMW pattern also appears to better preserve the influence of the regional topography and land cover on spring melt timing, as described in the following sections.  The mean correlation between the satellite MMOD PMW retrievals and in situ MMOD estimates derived from SNOTEL station SWE measurements (MMOD SWE ) was 0.49 with a bias of −3.63 days, indicating an earlier MMOD PMW . The mean correlation between MMOD PMW and the in situ air temperature measurement-based MMOD estimates (MMOD T ) was stronger (0.69) and showed a smaller (~1 day) bias. Aggregation of the correlation results by snow cover classification [39] showed generally stronger MMOD PMW mean correspondence in the colder Tundra/Taiga snow regime with the respective SNOTEL SWE and air temperature-based MMOD observations (0.62 and 0.81). The MMOD mean correlations were generally lower for Alpine/Prairie (0.31 and 0.61) and Maritime (0.35 and 0.41) snow regimes.
For each of the three snow classification zones, we investigated the influence of different landscape factors on the relative bias between the MMOD PMW retrievals and the SNOTEL MMOD SWE,T observations, including FW cover; fractional forest (FF) cover; terrain aspect and elevation, and topographic roughness index (TRI). In the Tundra/Taiga region, the MMOD PMW and MMOD SWE bias was moderately correlated with FF (r = −0.61, p < 0.1) and moderately correlated with FW (r = −0.86, p < 0.1) within a surrounding grid cell. In contrast, the MMOD PMW and MMOD T bias corresponded more strongly with terrain aspect northness (r = −0.87, p < 0.1). In the Alpine/Prairie region, MMOD PMW and MMOD SWE biases were strongly influenced by elevation (r = 0.92, p < 0.1) and TRI (r = 0.94, p < 0.1), whereas the MMOD PMW and MMOD T biases were only influenced by FW (r = 0.99, p < 0.1). In the Maritime region, MMOD PMW and MMOD SWE biases were influenced by FW (r = −0.74, p < 0.1), whereas the MMOD PMW and MMOD T biases were influenced by FW (r = −0.77, p < 0.1) and aspect northness (r = −0.68, p < 0.1). These results indicate that enhanced surface moisture influences the agreement between the MMOD PMW and SNOTEL observations, as the 37 GHz Tb observations are strongly sensitive to surface water within the satellite footprint [23]. The influence of FF cover particularly in the boreal interior regions was demonstrated, as microwave emissions from surface snow cover can be adversely affected by the overlying forest cover at both 19 and 37 GHz frequencies [42].

Snow Metric Distribution in the YRB
The average snow metric spatial distribution (1988-2016) indicated a topographic influence with generally earlier (later) MMOD PMW and SO PMW at lower (higher) elevations. The PMW record also showed significant interannual variability in the spring snow metrics, indicated by the extensive early spring onset during the exceptionally warm year of 2016 relative to a more climatological normal year in 2001. The SMD PMW distribution showed relatively less spatial and annual variability but a longer duration of spring snowmelt in the YRB upper headwaters and lower delta regions ( Figure 6). The spatial distribution of correlations between the snow metrics and annual May temperatures are shown in Figure 7. Here, we quantified the relationship between annual snow metrics and temperature using monthly aggregated air temperatures from the downscaled (20 km resolution) WRF Reanalysis (Table 2). Overall, all spring snow metrics showed generally significant but variable correlations with May temperatures in both sign On average, MMOD PMW in the YRB ranged from DOY 82 at lower elevations to DOY 153 at higher elevations, with a regional mean of DOY 113 ± 12. Alaska experienced record-setting warmth during the 2015/16 snow season (October-April), with statewide temperatures 4 • C above the mean [43]. During this year, the average MMOD PMW in the YRB was DOY 101 ± 16, which was 12 days earlier than the long-term mean (1988-2016) and 22 days earlier than in 2001.
SO PMW in the YRB on average ranged from DOY 116 to 162 between lower and higher elevations, with a regional mean value of DOY 138 ± 9. Like MMOD PMW , SO PMW showed a regional mean of DOY 127 ± 13 during the exceptionally warm year in 2016, which was approximately 11 days earlier than normal relative to the long-term mean (1988-2016), and 20 days earlier than in 2001.
The average spring SMD PMW in the YRB ranged spatially from 16 to 42 days, with a regional average of 25 ± 4 days over the long-term record. SMD PMW appeared to have a less distinct spatial distribution relative to MMOD PMW and SO PMW , but with longer durations in the YRB upper headwater and delta regions. A longer SMD PMW occurred during the unprecedented warm year in 2016 (26 ± 8 days) relative to 2001 (24 ± 9 days). The longer duration in 2016 is in line with other studies that identified a lengthening snowmelt season under warming conditions [44,45].
The spatial distribution of correlations between the snow metrics and annual May temperatures are shown in Figure 7. Here, we quantified the relationship between annual snow metrics and temperature using monthly aggregated air temperatures from the downscaled (20 km resolution) WRF Reanalysis (Table 2). Overall, all spring snow metrics showed generally significant but variable correlations with May temperatures in both sign and magnitude. MMOD PMW and temperature regressions conducted within 100 m elevational bins showed lower elevations having the greatest percentage of significant (p < 0.1) grid cells, ranging from 51-60% over the 0-700 m elevational range, with moderate correlations (−0.41 to −0.44). The strongest correlations occurred at higher elevations (1200-1600 m) and ranged from −0.45 to −0.47, although the area of significant grid cells was lower (30%). Overall, MMOD PMW had a significant relationship with temperature for 42% of the YRB, respectively.

Interaction between Snow Metrics and Discharge in the YRB
We iterated through 36 regressions for Eagle, Stevens Village and Pilot sub-basins. Regressions included quantiles (Q20, Q50, Q80) for each snow metric (MMOD, SO, and SMD) defined as the percent area of a catchment regressed against streamflow quantiles (Q20, Q50, Q80, and Peak), defined as the percent of cumulative annual flow at each gaging station. At the Eagle station, statistically significant (p < 0.1) regressions were observed only for the timing of MMODPMW quantiles (Q20MMOD, Q50MMOD, and Q80MMOD) and streamflow Q20Flow. Regressions were also significant for the same SO quantiles (Q20SO, Q50SO, and Q80SO) and streamflow Q20Flow (Appendix C). For MMOD, Q20MMOD had the highest correlation with the observed Q20Flow (r = 0.5) followed by Q50Flow (0.48) and Q80Flow (0.47). SOPMW regressions were higher relative to MMODPMW, with Q20SO having a moderately strong correlation of 0.61 with Q20Flow, suggesting a stronger relationship between SOPMW The SO PMW correlations with May temperatures showed significant relationships over 60 to 85 percent of the area at elevations below 1000 m in the YRB with correlations ranging from −0.40 to −0.57. The number of grid cells with significant SO PMW and spring temperature relationships were lower at higher elevations (1100-3000 m), with correlations ranging from −0.30 to −0.52. The SO PMW relationship with temperature was also predominantly negative over 87% of the YRB, indicating generally earlier (delayed) snowpack depletion in warmer (cooler) years. The SMD PMW relationship with May temperatures was relatively spatially complex due to the variable influence of temperature on MMOD PMW and SO PMW . SMD PMW was significantly correlated with temperature over 75% of the YRB, but with both positive and negative relationships for 22% and 53% of the domain, respectively. The positive SMD PMW temperature response was predominantly located at higher elevations (>800 m) characterized by deeper snowpack conditions [46].

Interaction between Snow Metrics and Discharge in the YRB
We iterated through 36 regressions for Eagle, Stevens Village and Pilot sub-basins. Regressions included quantiles (Q20, Q50, Q80) for each snow metric (MMOD, SO, and SMD) defined as the percent area of a catchment regressed against streamflow quantiles (Q20, Q50, Q80, and Peak), defined as the percent of cumulative annual flow at each gaging station. At the Eagle station, statistically significant (p < 0.1) regressions were observed only for the timing of MMOD PMW quantiles (Q20 MMOD , Q50 MMOD , and Q80 MMOD ) and streamflow Q20 Flow . Regressions were also significant for the same SO quantiles (Q20 SO , Q50 SO , and Q80 SO ) and streamflow Q20 Flow (Appendix C). For MMOD, Q20 MMOD had the highest correlation with the observed Q20 Flow (r = 0.5) followed by Q50 Flow (0.48) and Q80 Flow (0.47). SO PMW regressions were higher relative to MMOD PMW , with Q20 SO having a moderately strong correlation of 0.61 with Q20 Flow , suggesting a stronger relationship between SO PMW and streamflow compared to MMOD PMW at Eagle (Table 2).
Moving downriver from Eagle to Stevens Village, we found overall relationships strengthened between the MMOD PMW and SO PMW quantiles, and streamflow. At Stevens Village, the regressions showed significant relationships between the snow metric quantiles and Q20 Flow in addition to peak flow. The timing of Q80 MMOD correlated well with the timing of the hydrologic peak (r = 0.71). Moreover, the timing of Q50 SO and Q80 SO were also strongly correlated with the timing of the hydrologic peak, with respective correlations of 0.71 and 0.66. We also found a strong relationship between Q20 SO and the Q20 Flow (r = 0.68) at Stevens Village.
At the Pilot gaging station, Q20 MMOD no longer showed a significant relationship with streamflow. However, Q50 MMOD and Q80 MMOD were still moderately correlated with the timing of the streamflow peak, with respective correlations of 0.42 and 0.58. The SO PMW quantiles and the streamflow quantiles remained stronger than for MMOD PMW , with Q20 SO and Q50 SO having correlations of 0.52 and 0.63 with the timing of the peak flow. Overall, Q20 SO had the highest correlation with Q20 Flow at the Eagle station. However, for Stevens Village and Pilot stations, Q50 SO had the strongest relationship with the timing of the streamflow peak. For the 1988 to 2016 record, Q50 SO occurred, on average, 16 ± 16 days earlier than the streamflow peak at Eagle; the lag time was less at Stevens Village (9 ± 10 days) and longer at the Pilot (17 ± 9 days) station ( Figure 8).
quantiles and the streamflow quantiles remained stronger than for MMODPMW, with Q20SO and Q50SO having correlations of 0.52 and 0.63 with the timing of the peak flow. Overall, Q20SO had the highest correlation with Q20Flow at the Eagle station. However, for Stevens Village and Pilot stations, Q50SO had the strongest relationship with the timing of the streamflow peak. For the 1988 to 2016 record, Q50SO occurred, on average, 16 ± 16 days earlier than the streamflow peak at Eagle; the lag time was less at Stevens Village (9 ± 10 days) and longer at the Pilot (17 ± 9 days) station ( Figure 8). Overall, there is an initial influence of snow processes on discharge at the higher reaches, represented by Eagle. Snowmelt had a stronger influence on discharge in the middle YRB reaches at Stevens Village. Moving to Pilot, larger MMOD quantiles (Q50, Q80) influenced the peak flow but to a lesser degree relative to Q20SO and Q50SO.

Interaction between Snow Metrics and River Ice Breakup in the YRB
We examined the interaction between the satellite PMW-derived annual snow metrics and observed RIB dates using snow metric pixel values at 19 RIB measurement locations along the Koyukuk and Yukon rivers. Significant relationships (p < 0.1) between MMODPMW and RIB were identified at 11 of the 19 locations (Figure 9). The strongest relationships were found at Beaver (r = 0.75), Nulato (r = 0.62) and Allakaket (r = 0.6) (Appendix D). The remaining eight RIB locations had modest correlations with MMODPMW averaging 0.42. For these 11 RIB locations, the MMODPMW occurred 26 days earlier than the RIB. SMDPMW was statistically significant at only four RIB locations, including Russian Mountain (r = 0.52), Circle (r= 0.42), Holy Cross (r = 0.4) and Kaltag (r = 0.36) ( Table 3). Overall, there is an initial influence of snow processes on discharge at the higher reaches, represented by Eagle. Snowmelt had a stronger influence on discharge in the middle YRB reaches at Stevens Village. Moving to Pilot, larger MMOD quantiles (Q50, Q80) influenced the peak flow but to a lesser degree relative to Q20 SO and Q50 SO .

Interaction between Snow Metrics and River Ice Breakup in the YRB
We examined the interaction between the satellite PMW-derived annual snow metrics and observed RIB dates using snow metric pixel values at 19 RIB measurement locations along the Koyukuk and Yukon rivers. Significant relationships (p < 0.1) between MMOD PMW and RIB were identified at 11 of the 19 locations (Figure 9). The strongest relationships were found at Beaver (r = 0.75), Nulato (r = 0.62) and Allakaket (r = 0.6) (Appendix D). The remaining eight RIB locations had modest correlations with MMOD PMW averaging 0.42. For these 11 RIB locations, the MMOD PMW occurred 26 days earlier than the RIB. SMD PMW was statistically significant at only four RIB locations, including Russian Mountain (r = 0.52), Circle (r = 0.42), Holy Cross (r = 0.4) and Kaltag (r = 0.36) ( Table 3).  Q20SO was the strongest correlated snow metric to RIB and significant at all RIB locations, with a mean correlation of 0.77. There were also strong interactions with the other SO quantiles, but correlations diminished as quantiles increased. However, we again found the annual mean difference between Q20SO and RIB to occur within a day of each other, while the other SO quantiles occurred, on average, after RIB. Hence, the potential for the satellite-derived SO metric in forecasting RIB is limited.  The SO PMW snow metric was significantly correlated with RIB timing at all 19 observation locations, with an average correlation of 0.58. Nine of the 19 locations were exceptionally significant (p < 0.001), had a mean correlation of 0.71 and ranged between 0.63 (Circle and Eagle) and 0.81 (Beaver). On average for the 19 observation locations, SO PMW occurred less than a day before RIB, although for some sites and years SO PMW preceded RIB by up to nine days (Holy Cross) or followed RIB by up to 11 days (Eagle).
We next examined relations between the timing of the satellite PMW-derived snow metric quantiles and observed RIB on the premise that river ice conditions are responsive to the snowmelt runoff pulse contributed from the surrounding drainage basin. The PMW snow metric quantiles generally showed a stronger relationship with RIB than the streamflow quantiles. RIB timing was significantly correlated with Q20 MMOD (p < 0.01) and Q50 MMOD (p < 0.05) in the surrounding catchments at all 19 RIB observation locations. The correlation between Q20 MMOD and RIB averaged 0.61 and ranged from 0.48 (Galena) to 0.77 (Bettles) across the 19 RIB observation sites. At the 19 RIB locations, Q20 MMOD occurred an average of 27 days before the RIB. At the Bettles station, Q20 MMOD occurred only eight days before RIB on average, while Q20 MMOD preceded RIB by an average of 35 days at the Pilot station. The correlation between Q50 MMOD and RIB was slightly lower than for Q20 MMOD but still significant, ranging from 0.37 (Eagle) to 0.79 (Allakaket). As expected, the temporal window between Q50 MMOD and RIB dates was narrower than for Q20 MMOD , with Q50 MMOD preceding RIB by an average of only 16 days throughout the YRB (Table 4). Q20 SO was the strongest correlated snow metric to RIB and significant at all RIB locations, with a mean correlation of 0.77. There were also strong interactions with the other SO quantiles, but correlations diminished as quantiles increased. However, we again found the annual mean difference between Q20 SO and RIB to occur within a day of each other, while the other SO quantiles occurred, on average, after RIB. Hence, the potential for the satellite-derived SO metric in forecasting RIB is limited.
In addition to quantiles, we also derived a Contributing Peak metric defined as the date at which the most grid cells in each catchment had begun melting, as indicated from the MMOD PMW retrievals, but where the snowpack had not yet fully depleted (identified by SO PMW ). We found the Contributing Peak to be a strong indicator of RIB, showing significant (p < 0.05) relationships at all RIB locations, and with a mean correlation of 0.67. The Contributing Peak occurred after the RIB date at only five locations (Allakaket, Bettles, Dawson, Eagle, and Hughes). For the remaining 14 locations, the Contributing Peak occurred an average of seven days before RIB.

MMOD Algorithm Performance
The MMOD PMW retrieval method performed favorably in relation to independent SWE and temperature measurement-based MMOD estimates from in situ SNOTEL sites distributed across Alaska. There was an overall stronger mean relationship between MMOD PMW and MMOD T (r = 0.69, bias = −1 day) than with MMOD SWE (r = 0.49, bias = −3.63 days). The stronger temperature relationship was attributed to generally greater spatial and temporal heterogeneity in the in situ SWE measurements relative to air temperature, which can propagate to greater uncertainty in identifying the MMOD signal. The satellite retrievals may also falsely classify MMOD PMW in response to early season temporary melt events in lieu of the actual seasonal melt signal [22,47]. Our results identified the weakest correspondence with air temperature in the maritime regions and a strong correlation between the temperature bias and FW (r = −0.77, p < 0.1). This FW interaction in the maritime region suggests that enhanced regional moisture, such as rain-on-snow, can adversely affect the GRP algorithm performance [24]. Future work should address a more comprehensive approach to distinguish between rain-on-snow and early season melt events from the MMOD [24,46].
While the MMOD PMW retrievals showed generally less correspondence with MMOD SWE than MMOD T , we found SWE bias to be modestly correlated to FF in the Tundra/Taiga (r = −0.61, p < 0.1) and FW in both the Tundra/Taiga (r = −0.86, p < 0.1) and Alpine/Prairie (r = −0.74, p < 0.1). These results suggest MMOD PMW to be most strongly sensitive to the change in surface wetness conditions during seasonal transitions. Further, vegetation biomass contribution to the landscape freeze-thaw (FT) signal can inflate MMOD PMW errors where the timing of the vegetation canopy seasonal FT transition differs from the surrounding snow cover [42]. In more densely vegetated areas, we would expect to see much less direct snow signal and more vegetation FT signal [28,47]. Vegetation cover can also influence snowpack spatial heterogeneity and representativeness of in situ SNOTEL SWE measurements relative to the surrounding satellite footprint; thus, contributing to larger satellite-site differences [23].

Changes in Snowmelt Properties
We initially examined long-term trends in snowmelt properties across the YRB from 1988-2016 but found no temporally or spatially consistent patterns despite an overall declining trend in snow cover extent across the Northern Hemisphere [48,49]. Yet, our results indicate that both MMOD PMW and SO PMW occur earlier during anomalously warm years and, conversely, occur later during cooler years. This identified interaction between temperature and snow metrics is very important, as Alaska has experienced several annual temperature records broken over the last few years [43]. With anticipated warming, future lines of research should focus on how these snow metrics can inform the intensity of oncoming wildfire seasons, wildlife movements and ecosystem productivity.

Snow and Hydrologic Interactions
We found significant relationships (p < 0.01) between the satellite-derived snowmelt quantiles and the spring flood pulse indicated from in situ streamflow measurements within the major YRB catchments. Yet, the timing of Q50 SO and Q80 SO generally showed the strongest correspondence with the streamflow quantiles (Q20, Q50) but occurred after these events, degrading their forecast potential. The day of peak flow is used here to represent the spring flood pulse in the YRB and was strongly correlated to Q50 SO at Stevens Village (r = 0.71) and Pilot (r = 0.63). The peak flow consistently occurred after the Q50 SO , on average 16, 9, and 17 days earlier than the peak flow for Eagle, Stevens Village, and Pilot, respectively. We would anticipate a smaller lag time at Eagle; however, the associated stream gage station has relatively high interannual variability and showed a weak correlation to Q20 SO (r = 0.21). Regardless, other studies have placed greater importance on MMOD in understanding the timing of the Spring flood pulse [8,50], although our results indicate that SO quantiles are stronger predictors of peak flow in the YRB. Additional research can benefit by exploring a lower quantile to identify further relationships between snow metrics and their forecasting potential.
A novel component of this study was the derivation of the Contributing Peak metric, derived using MMOD and SO. This metric effectively defines the area that is actively melting, or not melting, across a basin. Of the 19 RIB locations, strong relationships between RIB and Contributing Peak were identified at 11 stations, with a mean correlation of 0.67. At the 11 locations, the Contributing Peak occurred seven days before the RIB on average. Hence, the strong correlation and occurrence before RIB, indicate the Contributing Peak as having forecasting potential. There were five stations (Allakaket, Bettles, Dawson, Eagle, and Hughes) where RIB occurred before the day of Contributing Peak. These stations all exist at the higher elevations, where we also noted RIB occurs earlier relative to the lower elevations. This likely indicates other driver, other than snowmelt, governing RIB.
Our ground observations identified an earlier RIB date in the upper reaches of the Yukon, relative to the lower reaches. These results are in line with other large rivers including the Mackenzie and Yenisey. However, for large high-latitude rivers flowing north, latitude plays a prominent role in an earlier RIB date in southern reaches relative to northern. As the YRB trends from east to west, the influence of latitude is minimized in favor of river channel characteristics such as slope, curvature and other channel-scale factors [51]. The amount of SWE and rate of warming/melting in the catchment, freeze/thaw status and moisture holding capacity of underlying soils also influences the magnitude and rate of melt pulse to rivers affecting the RIB [51,52]. Our satellite-derived snow metrics did not fully capture these additional factors, which may degrade the observed correspondence with discharge and RIB.
The satellite-derived snow metrics developed here have several intrinsic limitations. The moderate spatial resolution (6.25 km) captures the overall latitudinal and altitudinal patterns across the domain. However, a significant amount of landscape heterogeneity exists within each pixel, which influences algorithm performance. Currently, the snow metrics capture the timing of snow processes rather than amount of melt. Hence, the snow metrics lack a representation of SWE and snowmelt rates controlling runoff magnitude, routing and soil storage influencing discharge. The strength of the snow metrics resides in the use of daily Tb observations and their strong sensitivity to the rapid changes in surface wetness from snowmelt.

Conclusions
The seasonal and spatial timing in snowmelt properties are important controls on ecosystem and socioeconomic processes across Arctic boreal landscapes, including the YRB. However, capabilities for regional monitoring of snow phenology and its influence on the spring flood pulse and RIB are constrained by the sparse in situ station observation network and limited satellite resources. In this study, we presented a new satellite PMW MMOD retrieval method and dataset that exploits the differential response between 19V and 37V GHz Tb channels with daily repeat and 6.25 km resolution gridding from 1988 to 2016 over the YRB. The developed algorithm is physically-based, drawing from the established GRP method and incorporating a dynamic detection threshold adjusted according to surface FW cover. The MMOD PMW results compared favorably with the regional pattern and annual variability in the primary spring thaw signal from the FT-ESDR occurring about five days later. The MMOD PMW retrievals also showed favorable agreement with independent in situ SWE (r = 0.49, bias = −3.6 days) and air temperature (r = 0.69, bias = 1 day) measurement based MMOD estimates from the regional weather station network. The regional pattern in MMOD PMW performance was influenced by one or more landscape factors, including the prevailing climate and snow type, terrain complexity, forest cover (FF) and the fractional open water body cover (FW) within the satellite footprint.
The resulting satellite snow metrics developed in this study showed that MMOD extends from early March into early June across the YRB. Melt begins at lower elevations and reaches of the Yukon before progressing inland and into higher elevations of the upper YRB reaches. SO timing followed a similar pattern to MMOD, but with a delayed response, where SO generally occurred in March in the YRB coastal areas and lower elevations, and extended as late as July in the higher elevations and upper reaches of the basin. The SMD was generally longer where MMOD and SO occurred earlier. Conversely, SMD was the shortest (<20 days) at higher elevations and in more interior regions. The timing of all snow metrics was sensitive to relatively warm and cool years. For example, both MMOD and SO occurred approximately 12 days and 11 days earlier during the anomalous warm year in 2016, respectively.
Positive and significant regressions were identified between the spring flood pulse and both MMOD and SO. The timing of Q50 SO was found to be the best indicator of peak discharge at the Stevens Village and Pilot stations. Q50 SO also occurred before peak discharge at both locations, making it one of the better predictors of the spring flood pulse. Overall, regressions were stronger between the snow metrics and RIB, with Q20 SO having the best correlation at all stations. However, for all stations, on average, Q20 SO occurred within one day of the RIB, limiting its potential utility as a forecast indicator. The Q20 MMOD and Q50 MMOD metrics indicated stronger potential utility as early predictors of seasonal ice breakup in the YRB, showing strong correlations with RIB and preceding seasonal ice breakup by approximately 35 days and 16 days, respectively. Yet, the Contributing Peak possessed the strongest correlations and occurred seven days on average before RIB at 11 of the 19 locations.
Anomalously warm years like 2016, when the cold season temperatures were + 4 • C above normal, coincided with MMOD and SO timing that was more than 10 days earlier than normal. This interaction indicates that the snow metrics accurately reflect seasonal surface air temperature conditions represented from regional reanalysis data, which is congruent with other studies [51,53]. Hence, projected regional warming trends of approximately 1-3 • C in spring air temperatures by midcentury [9,54] are expected to promote generally earlier onset of MMOD and SO across the YRB, along with generally earlier ice breakup and earlier onset of the spring flood pulse. Such changes in snow and hydrologic processes will affect fluvial dynamics and fisheries, as well as terrestrial processes like spring carbon uptake and wildfires, making further research at the intersection of snow and ecologic processes of high value.
Author Contributions: C.G.P., P.B.K. and J.S.K. designed the study; C.G.P., P.B.K. and J.D., analyzed the data; all authors contributed to the writing and editing. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
We thank three anonymous reviewers for their thoughtful and constructive comments.

Conflicts of Interest:
The authors declare no conflict of interest.   Quantile derived from streamflow Q20 SO SO Q20