A Long-Term Passive Microwave Snowoff Record for the Alaska Region 1988–2016

Snowoff (SO) date—defined as the last day of observed seasonal snow cover—is an important governor of ecologic and hydrologic processes across Alaska and Arctic-Boreal landscapes; however, our understanding and capacity for the monitoring of spatial and temporal variability in the SO date is still lacking. In this study, we present a 6.25 km spatially gridded passive microwave (PMW) SO data record, complimenting current Alaskan SO records from Moderate Resolution Imaging Spectrometer (MODIS) and Landsat, but extending the SO record an additional 13 years. The PMW SO record was validated against in situ snow depth observations and showed favorable accuracy (0.66–0.92 mean correlations; 2–10 day mean absolute errors) for the major climate regions of Alaska. The PMW SO results were also within 10 days of finer spatial scale SO observational records, including Interactive Multisensor Snow and Ice Mapping System (IMS), MODIS, and Landsat, for a majority (75%) of Alaska. However, the PMW record showed a general SO delay at higher elevations and across the Alaska North Slope, and earlier SO in the Alaska interior and southwest regions relative to the other SO records. Overall, we assign an uncertainty +/−11 days to the PMW SO. The PMW SO record benefits from the near-daily temporal fidelity of underlying brightness temperature (Tb) observations and reveals a mean regional trend in earlier SO timing (−0.39 days yr−1), while significant (p < 0.1) SO trend areas encompassed 11% of the Alaska domain and ranged from −0.11 days yr−1 to −1.31 days yr−1 over the 29-year satellite record. The observed SO dates also showed anomalous early SO dates during markedly warm years. Our results clarify the pattern and rate of SO changes across Alaska, which are interactive with global warming and contributing to widespread permafrost degradation, changes in regional hydrology, ecosystems, and associated services. Our results also provide a robust means for SO monitoring from satellite PMW observations with similar precision as more traditional and finer scale observations.


Introduction
Snow is a defining regulator of the global energy budget and ecosystem function, governing both ecologic and hydrologic responses to climate variability. Warming across the high latitudes in the Northern Hemisphere has contributed to reduced snow cover [1], whereby the reduced surface albedo from less snow cover [2,3] enhances the solar energy loading that promotes additional thawing that reinforces further climate warming [1,[4][5][6]. In Alaska, the snow-free season has lengthened, prompting a cascade of environmental impacts [3], including earlier onset of snowmelt and seasonal flooding in spring; earlier and longer potential growing seasons [7,8]; changes in

Passive Microwave Record 1988-2016
We detected SO using the SSM/I(S) 19 and 37 GHz ascending orbit (~17:00 local solar time) 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 (NSIDC) [42]. This data record is a temporally extensive consortium of images from multiple sensors and platforms. The spatial resolution of the 19 and 37 GHz Tb channels used in this study are improved from the original 25-km resolution Tb format through spatially enhanced interpolation of overlapping Tb antenna patterns and regridding at 6.25 km (19 GHz) and 3.125 km (37 GHz) resolution in a northern polar EASE-grid format using the scatterometer image reconstruction (SIR) technique [43].
The SSM/I(S) Tb retrievals from different DMSP satellite platforms are highly correlated and only small improvements are made using correction techniques [44], allowing us to use overlapping observations from the DMSP F-series satellites during the SSM/I and SSMIS era to gap-fill missing coincident grid cells using a nearest-neighbor substitution. Any remaining missing days were gapfilled using a temporal linear interpolation [25]. The gap-filled daily Tb retrievals were clipped and reprojected to the Alaska Albers Conical Equal Area Projection using the North American Datum 1983 to match the same spatial extent and projection as the MODIS derived snow metrics for Alaska [13].

Deriving SO
We use the 19V-37V Tb differencing as the primary approach to detect SO [45], defined here as the last day of snow cover in spring. The resulting Tb difference is strongly positive and relatively stable prior to the seasonal onset of snowmelt indicated from in situ weather station observations

Passive Microwave Record 1988-2016
We detected SO using the SSM/I(S) 19 and 37 GHz ascending orbit (~17:00 local solar time) 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 (NSIDC) [42]. This data record is a temporally extensive consortium of images from multiple sensors and platforms. The spatial resolution of the 19 and 37 GHz Tb channels used in this study are improved from the original 25-km resolution Tb format through spatially enhanced interpolation of overlapping Tb antenna patterns and regridding at 6.25 km (19 GHz) and 3.125 km (37 GHz) resolution in a northern polar EASE-grid format using the scatterometer image reconstruction (SIR) technique [43].
The SSM/I(S) Tb retrievals from different DMSP satellite platforms are highly correlated and only small improvements are made using correction techniques [44], allowing us to use overlapping observations from the DMSP F-series satellites during the SSM/I and SSMIS era to gap-fill missing coincident grid cells using a nearest-neighbor substitution. Any remaining missing days were gap-filled using a temporal linear interpolation [25]. The gap-filled daily Tb retrievals were clipped and reprojected to the Alaska Albers Conical Equal Area Projection using the North American Datum 1983 to match the same spatial extent and projection as the MODIS derived snow metrics for Alaska [13].

Deriving SO
We use the 19V-37V Tb differencing as the primary approach to detect SO [45], defined here as the last day of snow cover in spring. The resulting Tb difference is strongly positive and relatively stable prior to the seasonal onset of snowmelt indicated from in situ weather station observations ( Figure 2). However, once the snow begins to melt the Tb difference between the two channels precipitously decreases, mirroring the snow ablation rate before reaching a seasonal minimum and becoming more temporally stable following snowoff and moving into the summer. The Tb difference reaches a seasonal minimum just after the snow has finished ablating, which we use to represent the SO condition. We observe a larger positive Tb difference between 19V and 37V GHz channels under colder winter temperatures because the emissivity is greatest at 19V relative to 37V in dry snow conditions. The Tb frequency difference is minimized as snow depth decreases and the liquid water content of the snowpack increases. Shortly after snowpack depletion the land surface is relatively wet and the microwave emissivity is generally larger at 37V than 19V GHz [46,47]; hence we define SO as the date when the 19V and 37V difference is at its lowest. colder winter temperatures because the emissivity is greatest at 19V relative to 37V in dry snow conditions. The Tb frequency difference is minimized as snow depth decreases and the liquid water content of the snowpack increases. Shortly after snowpack depletion the land surface is relatively wet and the microwave emissivity is generally larger at 37V than 19V GHz [46,47]; hence we define SO as the date when the 19V and 37V difference is at its lowest.
We use a dynamic Tb differencing approach by first calculating the average annual Tb difference value (sum) for June, which is used as a threshold to identify possible SO dates. We then classify the SO date where the daily Tb difference is less than sum and has the lowest value of all other potential selections ( Figure 2). We found that the Tb differencing approach was unable to provide reliable SO detection in coastal grid cells, which was attributed to greater water body contamination of the Tb retrievals over adjacent land areas [48]. The affected grid cells represent 2-3% of the study domain per year. For these instances we used an alternative gradient ratio polarization (GRP) algorithm [26,49] to detect SO. For the GRP approach, SO was classified for each grid cell as the day with lowest 4-day running GRP mean (rm) over each seasonal cycle (January-July).

Gridded SO Datasets
In addition to validating the PMW SO dates against SNOTEL measurements, we also compare the PMW SO with other satellite-derived Alaskan SO dates. The comparisons are quantified as We use a dynamic Tb differencing approach by first calculating the average annual Tb difference value (sum) for June, which is used as a threshold to identify possible SO dates. We then classify the SO date where the daily Tb difference is less than sum and has the lowest value of all other potential selections ( Figure 2). We found that the Tb differencing approach was unable to provide reliable SO detection in coastal grid cells, which was attributed to greater water body contamination of the Tb retrievals over adjacent land areas [48]. The affected grid cells represent 2-3% of the study domain per year. For these instances we used an alternative gradient ratio polarization (GRP) algorithm [26,49] to detect SO. For the GRP approach, SO was classified for each grid cell as the day with lowest 4-day running GRP mean (rm) over each seasonal cycle (January-July).

Gridded SO Datasets
In addition to validating the PMW SO dates against SNOTEL measurements, we also compare the PMW SO with other satellite-derived Alaskan SO dates. The comparisons are quantified as residuals Remote Sens. 2020, 12, 153 6 of 23 (days) showing the location and magnitude of PMW SO occurring earlier or later relative to the other SO records.

Alaskan Snow Metrics from MODIS
The Alaskan Snow Metrics cover Alaska, northwestern Canada, and the Russian Far East. The data originate from the MODIS Terra Snow Cover Daily L3 Global 500 m Grid data (MOD10A1) [50] from 2001-2016. The data are provided annually as 12 band GeoTIFFs-each band representing a different snow cover metric. The MODIS MOD10A1 product is a binary snow cover extent dataset derived from the NDSI. Data processing involves a 3-step filtering process to fill the data gaps, including a spatial and temporal filter, as well as a snow cycle filter [13]. In this study, we used the last day of continuous snow cover (LCLD) to represent SO and resampled each GeoTIFF image using the median LCLD value to match the same 6.25 km resolution and dimensions as the PMW-derived SO record.

Alaskan Snow Persistence
Snow persistence is defined as the transition from snow to snow-free surface conditions. The Alaskan annual snow persistence maps begin in 2001 and are roughly constrained to the Alaskan state boundaries [31]. This dataset is derived from Landsat (TM, ETM+, OLI) observations and provides a relatively high-resolution (30 m) estimation of annual SO across Alaska. After shadows and clouds are filtered from the 16-day Landsat optical reflectance data, daily snow cover extent is calculated using the NDSI. The remaining unmasked NDSI time series is then used in a regression tree model to predict SO in days of year (DOY) for each 30m grid cell [31]. The Landsat 30 m Alaska SO record used in this study extended from 2001-2017 and was resampled using the median value and returned the same 6.25 km grid projection as the PMW SO record. We also assigned a minimum +/−16 day precision threshold to the Landsat SO record based on the Landsat 16-day potential sampling frequency, which was used to assess the relative performance of the PMW SO record developed from this study.

IMS
The Interactive Multisensor Snow and Ice Mapping System (IMS) daily snow cover extent record covers different years and spatial grid formats, including 1 km (2014-present), 4 km (2004-present), and 24 km (1997-present) resolutions across the Northern Hemisphere. The IMS snow data are derived using expert interpretation of geostationary visible satellite imagery, polar orbiting multispectral satellite sensors, PMW sensors, and ground observations [51,52]. The IMS daily snow cover extent is one of the highest performing snow cover products available, likely attributed to its use of a trained analyst in combination with automated classifications and inclusion of ground observations, which decrease the influence of clouds [53]. In this study, we used the 4 km IMS snow cover extent product and calculated SO as the first day of 14 consecutive days of snow free conditions for each grid cell [30]. We assigned a +/−14 day precision threshold to the IMS derived SO record based on the 14-day moving window classification, which was used as a basis for evaluating our PMW SO record. The IMS data were also resampled to the same PMW 6.25 km resolution format using the median value. A summary of each satellite SO dataset used in this study is given in Table 1.

Ancillary Datasets
In situ snow depth observations were used to validate the PMW SO algorithm and compare uncertainty between each SO dataset. Like the IMS SO, we defined in situ SO as the first day of 14 consecutive days where snow depth was zero. We used snow depth observations from 31 SNOTEL stations across Alaska and constrained our validation and comparison period to 2004-2016, the common period that all gridded SO datasets share.
We examined the influence of three land surface properties potentially influencing the PMW SO accuracy, including surface fractional water inundation (FW), fractional forest cover (FF), and the elevation gradient within a grid cell [54,55]. We used the dynamic 5-km FW record derived from the Advanced Microwave Scanning Radiometer sensors (AMSR-E and AMSR2), calculated as 10-day FW outputs [56]; the FW record used here is derived from higher frequency (89 GHz) Tb retrievals than the PMW SO record, but has similar spatial resolution and microwave sensitivity to surface water. We calculated the annual mean FW cover during the summer months from 2002 to 2015 for each grid cell. The PMW performance was then evaluated along a regional gradient in surface water conditions based on the FW summer mean and coefficient of variation (CV) (Figure 3). Grid cells with significant open water cover (FW > 5 %) were masked from the SO validation to reduce the effect of surface water contamination of the Tb retrievals [57]. We also excluded grid cells where the FW exceeded a 70% CV threshold, indicating areas with large summer inundation variability. The total area excluded by the FW screening represented approximately 14% of the study domain ( Figure 3c). We also masked out perennial snow/ice areas identified in the Landsat snow persistence images, which represented an additional of the Alaska domain. The MODIS MOD44B V005 Vegetation Continuous Fields (VCF) product was used to represent the percent tree cover (FF) within each 6.25 km SO grid cell, as defined from the MODIS 250 m native VCF pixel resolution. The VCF has an estimated error of +/− 11.5% and is acceptable for our applications [58]. We also used the GTOPO30 digital elevation model (DEM) with 1 km native resolution to calculate the topographic complexity (TC) within each 6.25 km SO grid cell; where TC The MODIS MOD44B V005 Vegetation Continuous Fields (VCF) product was used to represent the percent tree cover (FF) within each 6.25 km SO grid cell, as defined from the MODIS 250 m native VCF pixel resolution. The VCF has an estimated error of +/−11.5% and is acceptable for our applications [58]. We also used the GTOPO30 digital elevation model (DEM) with 1 km native resolution to calculate the topographic complexity (TC) within each 6.25 km SO grid cell; where TC was calculated as the spatial standard deviation of the elevation within a cell.

SNOTEL Validation
The PMW SO results were validated against SO estimates derived from in situ daily weather station observations from Alaska NRCS SNOTEL sites (Table A1). The mean absolute error (MAE) was derived between each SNOTEL site based SO estimate and the satellite SO retrieval from the overlying PMW grid cell for each site location. Similar MAE metrics were also calculated for the MODIS, Landsat, and IMS SO datasets at each SNOTEL location. The MAE results were then used to evaluate PMW SO accuracy and performance relative to the other SO geospatial products. We also aggregated the SNOTEL annual SO results within each climate region. The SO correlation, MAE and root mean square error (RMSE) differences were then calculated for each climate region between the aggregated SNOTEL observations (independent variable) and the corresponding PMW results (dependent). The same statistics were calculated for the MODIS, Landsat, and IMS-based SO results, and used to evaluate PMW performance relative to the other datasets.

Residuals and Uncertainty
Pixel-wise residuals were calculated to quantify the sign and magnitude of differences between the PMW SO algorithm results (dependent variable) and the other gridded SO datasets (independent) from MODIS, Landsat, and IMS. The residuals were calculated over a 2004-2016 common validation period spanning the Alaska domain. The PMW SO performance was assessed through pixel-wise correlations and MAE differences relative to the other established SO geospatial records. We first calculated the per pixel correlations and MAE differences between the PMW SO results and corresponding annual SO records from MODIS, Landsat, and IMS. The mean correlations and MAE differences are then aggregated by FW, FF, and TC categories to clarify the impact of these features on PMW SO performance.

SO Trend Analysis
We performed temporal linear trend analysis for each grid cell to estimate the sign and strength of PMW SO changes over the satellite record (1988-2016). A Theil-Sen (TS) approach was used to compute the trend slope as the median value between each paired point [59]. The TS is a robust trend estimation for non-parametric time-series that do not meet the assumptions for data normality and homoscedasticity, while also being resilient to outlier effects. Significant trends were identified using the Mann-Kendall Test (MKT) and were considered significant where p < 0.1 [60].

SNOTEL Site Validation
The MAE differences between the geospatial SO records and corresponding SO estimates defined from the in situ SNOTEL site measurements are presented in Figure 4, and stratified according to the corresponding FW, FF, and TC levels represented for each site location. Overall, the PMW SO results show a similar level of accuracy as the other benchmark SO datasets relative to the Alaska SNOTEL site observations. The mean MAE for the PMW SO results was 11.4 days across the 31 SNOTEL stations, which was intermediate between IMS and Landsat-derived SO (respective mean MAE of 9.2 Remote Sens. 2020, 12, 153 9 of 23 and 10.8 days), and the MODIS SO record (mean MAE of 12.8 days). The PMW SO results were also within the estimated 14-16 day precision level of the SO benchmark observations, indicating favorable accuracy and performance. The PMW SO MAE values showed the strongest correspondence with the Landsat results (R 2 = 0.74), followed by MODIS (R 2 = 0.72) and then IMS (R 2 = 0.57). The favorable SO performance occurred despite the much coarser spatial resolution of the PMW retrievals relative to the SO benchmark records, which was attributed to strong PMW sensitivity to snow cover conditions and the near-daily temporal fidelity of the microwave retrievals. The generally consistent SO performance also indicates a common level of uncertainty among the SO records despite the different underlying sensors and methods used for the retrievals. Thus, the performance may reflect similar spatial disparity between the landscape level SO records and local SNOTEL site measurements.  There were no obvious impacts from the selected landscape heterogeneity factors on the SO MAE performance over the sparse SNOTEL locations (Figure 4). However, we compared the surface characteristics (FW, FF, TC) at each SNOTEL site location across the different SO datasets where the MAE was above or below a 14-day precision threshold. The results indicate that the PMW SO performance was optimal (MAE < 14 days) for locations where the mean FW level was low (0.94 %), and degraded (>14 days) where the FW level was larger (4.3 %). These results are generally consistent with the expected strong PMW sensitivity to surface water [56]. However, the MODIS SO record also showed a marked decrease in MAE performance from low (0.78 %) to high (3.47%) FW levels, while the Landsat and IMS records showed only marginal SO sensitivity to FW (Table A2). The PMW SO There were no obvious impacts from the selected landscape heterogeneity factors on the SO MAE performance over the sparse SNOTEL locations (Figure 4). However, we compared the surface characteristics (FW, FF, TC) at each SNOTEL site location across the different SO datasets where the MAE was above or below a 14-day precision threshold. The results indicate that the PMW SO performance was optimal (MAE < 14 days) for locations where the mean FW level was low (0.94 %), and degraded (>14 days) where the FW level was larger (4.3 %). These results are generally consistent with the expected strong PMW sensitivity to surface water [56]. However, the MODIS SO record also showed a marked decrease in MAE performance from low (0.78 %) to high (3.47%) FW levels, while the Landsat and IMS records showed only marginal SO sensitivity to FW (Table A2). The PMW SO performance was also optimal for SNOTEL locations where terrain complexity (TC) was relatively low (mean TC = 108m), and degraded at higher TC levels (112 m). The IMS and MODIS SO records also showed similar decreases in MAE performance in more complex terrain and areas of moderate FF.
The distribution of SNOTEL stations in Alaska is concentrated in the Cook Inlet and Central Interior. Four SNOTEL stations showed MAE values exceeding 21 days across each SO dataset, namely Paragon Creek, Telaquana Lake, Upper Tsaina River, and Gobblers Knob ( Figure 5). The higher levels of MAE uncertainty is consistent with greater landscape heterogeneity at these locations. Paragon Creek is located at the southern edge of the Seward Peninsula within a landscape of thermokarst lakes and a mean FW of 2%. Telaquana Lake is located at the southern end of the Alaska Range and has an FW value of 13%. The Upper Tsaina River SNOTEL station is located in the Chugach Mountains and has the largest TC value (308 m) of any other SNOTEL validation site. Lastly, Gobblers Knob is located in the south facing foothill slopes of the Brooks Range where there is likely considerable wind redistribution of snow.  The PMW SO results were compared with the other geospatial SO benchmarks and the spatially aggregated SNOTEL SO observations within three major Alaska climate zones, including North Slope, Western, and Interior regions; the resulting annual distributions are plotted for each subregion over the 2004 to 2016 study period ( Figure 6). Overall, the PMW SO results showed similar spatial distributions and interannual variability in relation to the other geospatial SO records and available SNOTEL observations. However, the PMW and IMS results were generally more similar with each other and with the SNOTEL observations, while the Landsat and MODIS SO records showed a relatively early SO bias. The SNOTEL SO observations indicated early SO conditions in The PMW SO results were compared with the other geospatial SO benchmarks and the spatially aggregated SNOTEL SO observations within three major Alaska climate zones, including North Slope, Western, and Interior regions; the resulting annual distributions are plotted for each sub-region over the 2004 to 2016 study period ( Figure 6). Overall, the PMW SO results showed similar spatial distributions and interannual variability in relation to the other geospatial SO records and available SNOTEL observations. However, the PMW and IMS results were generally more similar with each other and with the SNOTEL observations, while the Landsat and MODIS SO records showed a relatively early SO bias. The SNOTEL SO observations indicated early SO conditions in 2005 over the Alaska interior region consistent with an anomolous warm spring [61]; however, the geospatial SO records indicated more normal conditions in this region relative to the other years of record.
Remote Sens. 2020, 11, x FOR PEER REVIEW 12 of 24 geospatial SO records indicated more normal conditions in this region relative to the other years of record. Temporal correlations between the geospatial SO datasets and SNOTEL SO observations were largely significant (p < 0.1) and similar within each climate region ( Table 2). The strongest and most significant (p < 0.05) correlations occurred in the North Slope region, despite the shorter SNOTEL record, and led by the PMW and Landsat SO results (R = 0.92 and 0.95, respectively). In the interior region, the PMW SO correlation (0.72) was similar to Landsat (0.71) and intermediate between the MODIS (0.76) and IMS (0.74) results. In the Western region, the PMW results showed the lowest correlation (0.66), while the IMS record showed the best correspondence (0.72). The MAE (RMSE) uncertainty between the geospatial SO data and SNOTEL observations ranged from 2-14 (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14) days over the Alaska climate regions and were within the 14-16 day precision levels estimated for the IMS and Landsat SO records. In the interior region, the PMW and IMS SO records showed similar MAE (RMSE) differences of 10 (11) days that were larger than the MODIS and Landsat uncertainty of 5-6 (7-8) days. In contrast, the North Slope region showed lower PMW and IMS SO uncertainty of 2 and 4 (3 and 5) days compared with larger (14 and 11 day) uncertainty from the MODIS and Landsat records in this region. Temporal correlations between the geospatial SO datasets and SNOTEL SO observations were largely significant (p < 0.1) and similar within each climate region ( Table 2). The strongest and most significant (p < 0.05) correlations occurred in the North Slope region, despite the shorter SNOTEL record, and led by the PMW and Landsat SO results (R = 0.92 and 0.95, respectively). In the interior region, the PMW SO correlation (0.72) was similar to Landsat (0.71) and intermediate between the MODIS (0.76) and IMS (0.74) results. In the Western region, the PMW results showed the lowest correlation (0.66), while the IMS record showed the best correspondence (0.72). The MAE (RMSE) uncertainty between the geospatial SO data and SNOTEL observations ranged from 2-14 (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14) days over the Alaska climate regions and were within the 14-16 day precision levels estimated for the IMS and Landsat SO records. In the interior region, the PMW and IMS SO records showed similar MAE (RMSE) differences of 10 (11) days that were larger than the MODIS and Landsat uncertainty of 5-6 (7-8) days. In contrast, the North Slope region showed lower PMW and IMS SO uncertainty of 2 and 4 (3 and 5) days compared with larger (14 and 11 day) uncertainty from the MODIS and Landsat records in this region.

Residuals Distributed across Alaska
The pattern of mean residual differences (2004-2016) of the pixel-wise regression relationsips between PMW SO (dependent variable) and the other (independent) SO geospatial data over Alaska are presented in Figure 7 and summarized in Table 3. The PMW results show a similar pattern of differences with the other geospatial SO records. PMW SO timing generally occurs later than the other records over colder Alaskan tundra areas of the North Slope and northern Seward Penninsula, and in higher elevation areas of the Brooks, Wrangell, Chugach, and Alaska ranges. The PMW residuals show contrasting earlier SO conditions relative to the other benchmark records over the Alaska Interior, southern, and southwest regions. Notably, the PMW SO differences are more modest relative to the IMS and are generally higher for both Landsat and MODIS in both positive and negative directions. Overall, the PMW SO results show relatively small residual differences that are within ±10 days of the other benchmark records for a majority of the Alaska domain, including IMS (82% of study area), Landsat (70%), and MODIS (73%).

Residuals Distributed across Alaska
The pattern of mean residual differences (2004-2016) of the pixel-wise regression relationsips between PMW SO (dependent variable) and the other (independent) SO geospatial data over Alaska are presented in Figure 7 and summarized in Table 3. The PMW results show a similar pattern of differences with the other geospatial SO records. PMW SO timing generally occurs later than the other records over colder Alaskan tundra areas of the North Slope and northern Seward Penninsula, and in higher elevation areas of the Brooks, Wrangell, Chugach, and Alaska ranges. The PMW residuals show contrasting earlier SO conditions relative to the other benchmark records over the Alaska Interior, southern, and southwest regions. Notably, the PMW SO differences are more modest relative to the IMS and are generally higher for both Landsat and MODIS in both positive and negative directions. Overall, the PMW SO results show relatively small residual differences that are within ±10 days of the other benchmark records for a majority of the Alaska domain, including IMS (82% of study area), Landsat (70%), and MODIS (73%).

Assigning Uncertainty
We calculated the pixel-wise correlations and MAE differences between the PMW SO results and each geospatial SO benchmark dataset separately to evaluate the relative performance of the PMW record over the entire Alaska domain. The mean of the subsequent correlations and MAE results were then aggregated by the different landscape factors (FW, FF, TC) to evaluate the PMW SO performance and potential degradation because of greater landscape complexity relative to the finer spatial resolution SO benchmark datasets (Figure 8). The analysis here is similar to the SNOTEL validation ( Figure 4), but encompasses the entire Alaska domain rather than just the sparse SNOTEL site locations. As anticipated, we found a general decrease in PMW SO performance (MAE and correlation) with increasing FW cover. Approximately 90% of the Alaska domain showed FW levels below 6%, which contributed to favorable PMW SO average MAE (12 days) and correlation (0.53) values relative to the benchmark SO records (Figure 8a).

PMW SO Trends and Anomalies
We performed an analysis to determine the potential trends in annual SO timing over Alaska as defined from the long-term (1988-2016) PMW satellite record. The MKT results show both earlier and later SO trends across the Alaska domain (Figure 9). Approximately 51% and 42% of the domain showed trends toward, respective earlier and later SO timing, with 18% of these area trends significant at the 90% level. The later SO trends predominantly occur in the Russia Far East, North Slope, and Interior regions, while the earlier SO trends mainly occur near coastal margins. The relatively low level of significant trend areas is attributed to the large characteristic climate variability in Alaska relative to the more subtle environmental trends represented within the 29-year satellite record. For areas that are significant, the mean earlier trend was −0.39 days yr −1 and ranged between −1.31 and −0.11 days yr −1 . Likewise, the mean later trend was 0.55 days yr −1 and ranged between 0.14 Surprisingly, the PMW SO results showed better MAE and correlation performance with increasing forest cover (FF). This interaction is partially attributed to generally closer agreement among the SO datasets in forested areas. Furthermore, 90% of the Alaska domain had FF coverage of less than 40% and PMW SO average MAE differences of 11 days (r = 0.58) (Figure 8b). The PMW SO results showed degraded performance relative to the other SO benchmark datasets in areas with greater terrain complexity (TC). However, the mean MAE difference was 11 days (r = 0.55) within 90% of the Alaska domain, which coincided with TC levels of less than 209 m. Notably, the performance was degraded where TC was between 0-20 m, which was attributed to greater PMW uncertainty contributed from water contamination of the Tb and SO retrievals in low elevation coastal areas (Figure 8c). Overall, the average MAE uncertainty of the PMW SO record was 11 days for 90% of the Alaska domain despite the diversity of landscape conditions represented from the FW, FF, and TC factors. The PMW performance over Alaska relative to the other geospatial SO datasets were similar to the accuracy level defined from the sparse SNOTEL validation sites (Section 3.1). The PMW SO performance was also within the estimated 14-16 day uncertainty threshold of the IMS and Landsat SO benchmark records.

PMW SO Trends and Anomalies
We performed an analysis to determine the potential trends in annual SO timing over Alaska as defined from the long-term (1988-2016) PMW satellite record. The MKT results show both earlier and later SO trends across the Alaska domain ( Figure 9). Approximately 51% and 42% of the domain showed trends toward, respective earlier and later SO timing, with 18% of these area trends significant at the 90% level. The later SO trends predominantly occur in the Russia Far East, North Slope, and Interior regions, while the earlier SO trends mainly occur near coastal margins. The relatively low level of significant trend areas is attributed to the large characteristic climate variability in Alaska relative to the more subtle environmental trends represented within the 29-year satellite record. For areas that are significant, the mean earlier trend was −0.39 days yr −1 and ranged between −1.31 and −0.11 days yr −1 . Likewise, the mean later trend was 0.55 days yr −1 and ranged between 0.14 and 0.55 days yr −1 .
Remote Sens. 2020, 11, x FOR PEER REVIEW 16 of 24 Water years 2015 and 2016 were two of the warmest years of the 29-year satellite record and are analogous to the projected future increase in mean temperatures across Alaska [62]. We evaluated the annual SO response during these anomalous warm years as a proxy for understanding future SO behavior in a warmer climate. The pixel-wise SO anomalies were derived for water years 2015 and 2016 to quantify their annual departures from average conditions defined from the 29-year PMW record ( Figure 10). The respective average surface air temperatures from January to April were 3 °C and 6 °C above-normal over the study domain for 2015 and 2016 [62].
In 2015, 47 percent (1.6 × 10 6 km 2 ) of the domain showed anomalous early SO conditions and had an average SO day of 133, 9 days earlier than normal conditions. The anomalous early SO timing predominantly occurred over the Alaskan Arctic, but was offset by delayed SO timing over portions of the Canadian Arctic, Russian Far East, and central Alaska; where the combined area of SO delay Water years 2015 and 2016 were two of the warmest years of the 29-year satellite record and are analogous to the projected future increase in mean temperatures across Alaska [62]. We evaluated the annual SO response during these anomalous warm years as a proxy for understanding future SO behavior in a warmer climate. The pixel-wise SO anomalies were derived for water years 2015 and 2016 to quantify their annual departures from average conditions defined from the 29-year PMW record ( Figure 10). The respective average surface air temperatures from January to April were 3 • C and 6 • C above-normal over the study domain for 2015 and 2016 [62].

Discussion
The satellite PMW SO retrieval method and results described in this study performed well in relation to independent SO estimates derived from in situ SNOTEL station measurements across Alaska (r = 0.66−.92; MAE = 2-10 days). The PMW SO performance at these sites was also within ±10 days of other established geospatial SO records from the IMS, Landsat, and MODIS. These results are consistent with a prior study reporting similar PMW-derived SO accuracy over the European borealarctic [30]. Residual differences between the PMW SO results and more established geospatial SO records were also within ±10 days across a majority (≥70%) of the entire Alaska domain, and below the 14-16 day estimated precision of the Landsat and MODIS derived SO products. Regional differences between the PMW SO retrievals and other SO benchmarks were attributed to the different sensitivity of optical and microwave frequencies to snow and other landscape features, as well as the different classification algorithms, spatial footprints and temporal fidelity of the observations used for the SO retrievals.
We examined the pattern of PMW SO relative errors in relation to selected landscape features expected to contribute to PMW uncertainty, including topographic complexity (TC), fractional water (FW), and forest cover (FF) within each 6.25 km resolution PMW grid cell. Our results indicate a general decline in PMW SO performance with increasing TC, defined as the spatial standard deviation of the elevation within a grid cell. These results are consistent with the coarser footprint PMW Tb retrievals being more constrained in their ability to characterize sub-grid scale SO heterogeneity in complex landscapes relative to the finer scale optical remote sensing observations. Nevertheless, the Tb product used for this study [42,65] provides a 2-4-fold spatial resolution enhancement over other available standard PMW products for the same Tb frequencies [26,57,66], while preserving the benefit of near daily Tb temporal fidelity.
We found a smaller than expected FW influence on PMW SO uncertainty (e.g., Figure 8), despite surface water having a strong influence on microwave emissions [56]. We attributed the lower-thanexpected FW sensitivity to the screening of higher FW grid cells prior to the SO algorithm implementation, and the predominantly low FW cover of the remaining 6.25 km resolution Tb grid cells used for the SO classification. However, our results indicated greater FW cover (1% increase) during anomalous warm years relative to cool years over the domain, suggesting that the enhanced water may contribute to greater PMW SO uncertainty in warmer years. High springtime FW in Alaska may be attributed to flooding of low-lying areas enhanced by snow runoff generated during a short snowmelt period when the soil is still frozen [67] and enhanced melting of permafrost [68]. Variability in annual FW is largely a response of spatial heterogeneity in climate, surface topography and permafrost conditions, but is largely governed by seasonal temperatures in the Arctic [69,70]. In 2015, 47 percent (1.6 × 10 6 km 2 ) of the domain showed anomalous early SO conditions and had an average SO day of 133, 9 days earlier than normal conditions. The anomalous early SO timing predominantly occurred over the Alaskan Arctic, but was offset by delayed SO timing over portions of the Canadian Arctic, Russian Far East, and central Alaska; where the combined area of SO delay represented 16% (5.6 × 10 5 km 2 ) of the domain, despite generally warmer conditions.
While 2016 was also an anomalous warm year, the regional pattern of SO anomalies was quite different from the prior water year. In 2016, 65 percent (1.2 × 10 6 km 2 ) of the domain showed anomalous early SO timing, which unlike 2015, was largely distributed across central Alaska and the Russian Far East; early SO anomaly areas in 2016 had an average SO day of 129, 13 days earlier than that overall average. Unlike the prior year, only 1% of the domain showed delayed SO timing. However, the pattern of SO delay in 2016 was largely distributed over the Alaskan Arctic and central Yukon, Canada. Collectively, these results allude to a general seasonal advance in regional SO timing under projected warming over the next few decades, but with large spatial and interannual variability driven by landscape heterogeneity in precipitation and the surface energy budget. These results also indicate the important role of both temperature and precipitation, interacting with landscape features and regional ocean-atmosphere circulation patterns, in affecting SO timing [63,64].

Discussion
The satellite PMW SO retrieval method and results described in this study performed well in relation to independent SO estimates derived from in situ SNOTEL station measurements across Alaska (r = 0.66−.92; MAE = 2-10 days). The PMW SO performance at these sites was also within ±10 days of other established geospatial SO records from the IMS, Landsat, and MODIS. These results are consistent with a prior study reporting similar PMW-derived SO accuracy over the European boreal-arctic [30]. Residual differences between the PMW SO results and more established geospatial SO records were also within ±10 days across a majority (≥70%) of the entire Alaska domain, and below the 14-16 day estimated precision of the Landsat and MODIS derived SO products. Regional differences between the PMW SO retrievals and other SO benchmarks were attributed to the different sensitivity of optical and microwave frequencies to snow and other landscape features, as well as the different classification algorithms, spatial footprints and temporal fidelity of the observations used for the SO retrievals.
We examined the pattern of PMW SO relative errors in relation to selected landscape features expected to contribute to PMW uncertainty, including topographic complexity (TC), fractional water (FW), and forest cover (FF) within each 6.25 km resolution PMW grid cell. Our results indicate a general decline in PMW SO performance with increasing TC, defined as the spatial standard deviation of the elevation within a grid cell. These results are consistent with the coarser footprint PMW Tb retrievals being more constrained in their ability to characterize sub-grid scale SO heterogeneity in complex landscapes relative to the finer scale optical remote sensing observations. Nevertheless, the Tb product used for this study [42,65] provides a 2-4-fold spatial resolution enhancement over other available standard PMW products for the same Tb frequencies [26,57,66], while preserving the benefit of near daily Tb temporal fidelity.
We found a smaller than expected FW influence on PMW SO uncertainty (e.g., Figure 8), despite surface water having a strong influence on microwave emissions [56]. We attributed the lower-than-expected FW sensitivity to the screening of higher FW grid cells prior to the SO algorithm implementation, and the predominantly low FW cover of the remaining 6.25 km resolution Tb grid cells used for the SO classification. However, our results indicated greater FW cover (1% increase) during anomalous warm years relative to cool years over the domain, suggesting that the enhanced water may contribute to greater PMW SO uncertainty in warmer years. High springtime FW in Alaska may be attributed to flooding of low-lying areas enhanced by snow runoff generated during a short snowmelt period when the soil is still frozen [67] and enhanced melting of permafrost [68]. Variability in annual FW is largely a response of spatial heterogeneity in climate, surface topography and permafrost conditions, but is largely governed by seasonal temperatures in the Arctic [69,70].
Surprisingly, the PMW SO results were relatively robust to increasing FF cover and even indicated improved SO consistency at higher FF levels relative to the other benchmark observations. These results contrast with prior studies indicating elevated Tb levels over Alaska interior forest regions, relative to lower Tb values expected over snow covered areas [71]. Given that the PMW SO algorithm is governed by the 19 and 37 GHz Tb difference, the expected elevation in Tb values and reduced Tb frequency gradient would be expected to reduce the effective GRP signal-to-noise ratio over dense forests [46]. Forests are also expected to attenuate microwave emissions from surface snow cover [72], further contributing to PMW SO uncertainty. The PMW record shows relatively early SO bias over much of the Alaska interior, which may reflect uncertainty contributed from the extensive forest cover in the region. The boreal forest structure tends to be closed [73], which may contribute to a stronger signal from the underlying snowpack. A stronger-than-expected microwave snow signal, coupled with significant uncertainties from optical remote sensing of snow in forested regions [21] may explain the apparent robust PMW SO response to FF.
Other factors may contribute to PMW SO uncertainty beyond the limited landscape attributes examined in this study. For example, a PMW tendency toward earlier SO timing in southwest Alaska relative to the other SO records may reflect generally shallow seasonal snow cover and tundra shrub conditions in this region, where a shallow snowpack (i.e., <5 cm depth or <10mm SWE) is largely transparent to PMW frequencies below 25 GHz [23]; while the PMW frequencies used in this study are less sensitive to shallow snow cover, the corresponding snow signal is much stronger for optical sensing in the absence of potential degradation from overlying vegetation, low solar illumination or clouds. The PMW results indicated a general 5-10 day SO delay relative to the other geospatial SO records in colder climate areas, including the Arctic North Slope, Brooks, and Alaska Ranges. This difference is within the estimated uncertainty of the SO records from Landsat and MODIS, which may be further elevated by optical remote sensing retrieval errors under unfavorable conditions including persistent cloud cover. The MODIS LCLD has been noted to precede snow depth observations [13] and is in agreement with our results showing generally early MODIS SO bias relative to the SNOTEL SO observations in the Alaska North Slope and Western regions. Both MODIS and Landsat may show similar SO uncertainty in the region because of extensive data loss from persistent cloud cover formed under orographic conditions in mountain landscapes [74], as well as shadows created from topographic complexity [75]. The PMW regional trend analysis indicated generally earlier SO conditions over approximately 51% of the domain, but with only 9% of the negative trend areas significant (p < 0.1). These results are in general agreement with several studies showing a decline in snow cover extent and duration at higher latitudes [8, 76,77], but also more diverse and of lower magnitude than expected given the Arctic amplification and increasing temperatures [3,78]. Our PMW SO trends may be tempered by the purported global warming hiatus which extended from 1998-2013 [79,80]. The SO trends may also be influenced by the large characteristic regional climate variability superimposed on a smaller longer-term environmental trend, and represented by a limited 29-year satellite record. However, the observed annual SO anomalies generally coincide with documented warm and cool annual temperature anomalies in the region, as temperatures have been identified as a major driver of SO timing at lower elevations [76]. As temperatures across the Arctic are projected to increase into the future [3], greater variability in snow cover conditions, including SO date, are expected.

Summary and Conclusions
The SO date is an important governor of ecologic and hydrologic processes across Alaska and Arctic-Boreal landscapes, yet a better understanding of regional patterns and trends in SO timing has been limited by a regionally sparse in situ observation network and uncertainties in geospatial SO records derived from available satellite observations and models. Uncertainties in the geospatial SO records stem from multiple factors, including different retrieval methods with variable sensitivity, spatial resolution, timing and duration of observations. In this study, we presented a satellite PMW SO retrieval method that exploits the differential response between daily 19V and 37V GHz Tb channels with daily repeat and 6.25 km resolution gridding from 1988-2016 over Alaska. The resulting PMW SO data record compliments other satellite-based SO records derived from optical remote sensing, including MODIS [13] and Landsat [31], by exploiting the unique microwave sensitivity to snow cover while extending the SO record over a longer (29-year) period. The PMW SO retrievals benefit from daily satellite Tb observations at higher latitudes (>50 • N), but with relatively coarse spatial gridding compared with optical remote sensing-based products with finer (30 m to 500 m) spatial resolution but much coarser temporal sampling and compositing of the underlying observations because of orbital geometry, illumination, and cloud-atmosphere constraints.
The PMW SO record showed favorable accuracy in relation to independent in situ SO observations from Alaska SNOTEL weather stations (r = 0.66-0.92; MAE = 2-10 days). The PMW accuracy at these sites was also consistent with the performance of other geospatial SO benchmarks obtained from the IMS, MODIS, and Landsat. Analysis of regression residuals between the PMW SO retrievals and other geospatial SO benchmarks indicated relative PMW SO consistency within ±10 days of the other SO records over an average 75% of the study domain, which was also within the estimated precision of the SO benchmark observations. The PMW retrievals tended to show delayed SO timing relative to the other geospatial SO records in colder tundra and mountain areas, including the North Slope, Brooks, and Alaska Ranges; however, the PMW SO timing was comparatively early over the Alaska interior and southwest. Regional differences in SO timing were attributed to the underlying PMW and optical sensor retrievals with different sampling footprints, varying sensitivity to snow properties and potential signal degradation from cloud-atmosphere effects, surface water, vegetation cover, and terrain complexity. Despite these differences, the PMW SO performance over the entire domain was within ±11 days of the other geospatial records, which was consistent with the in situ SNOTEL validation assessment and within the estimated precision of the SO benchmark observations.
A regional trend analysis of the 29-year PMW SO record indicated generally earlier SO timing across the domain, with a stronger tendency toward early SO conditions in anomalous warm years. Hence, as the high latitudes continue to warm, uncertainty remains in our understanding of how boreal-Arctic ecosystems and hydrological function will respond to earlier SO conditions. Our results also show snow cover persistence and SO delays in many areas, even during anomalous warm years, indicating potential resilience and refugia from climate change in the region. The global coverage and