Estimating Meltwater Drainage Onset Timing and Duration of Landfast Ice in the Canadian Arctic Archipelago Using AMSR-E Passive Microwave Data

: Meltwater drainage onset (DO) timing and drainage duration (DD) related to snowmelt-water redistribution are both important for understanding not only the Arctic energy and heat budgets but also the salt / heat balance of the mixed layer in the ocean and sea-ice ecosystem. We present DO and DD as determined from the time series of Advanced Microwave Scanning Radiometer-Earth observing system (AMSR-E) melt pond fraction (MPF) estimates in an area with Canadian landfast ice. To address the lack of evaluation on a day-by-day basis for the AMSR-E MPF estimate, we ﬁrst compared AMSR-E MPF with the daily Medium Resolution Imaging Spectrometer (MERIS) MPF. The AMSR-E MPF estimate correlates signiﬁcantly with the MERIS MPF ( r = 0.73–0.83). The estimate has a product quality similar to the MERIS MPF only when the albedo is around 0.5–0.7 and a positive bias of up to 10% in areas with an albedo of 0.7–0.9, including melting snow. The DO / DD estimates are determined by using a polynomial regression curve ﬁtted on the time series of the AMSR-E MPF. The DOs / DDs from time series of the AMSR-E and MERIS MPFs are compared, revealing consistency in both DD and DO. The DO timing from 2006 to 2011 is correlated with melt onset timing. To the best of our knowledge, our study provides the ﬁrst large-scale information on both DO timing and DD.


Introduction
The surface melt of Arctic sea ice begins in early summer, with liquid water pools (called melt ponds) forming on the ice.These melt ponds reduce the regional albedo, enhance heat uptake on the sea ice [1], and have an impact on simulations of an accelerating under-ice melt [2,3].Schröder et al. [4] found that an accurate accounting of melt ponds during the melt season leads to better predictions of sea ice extent in the fall.Thus, a better understanding of the temporal and spatial variability in melt pond fraction (MPF) facilitates an investigation of the function of melt ponds in the Arctic and in the earth's climate system [5][6][7].
A critical component of the time series of the MPF is the existence of pathways for meltwater drainage to sea level through the sea ice.The MPF is largely controlled by a hydraulic balance of meltwater inflows and outflows [8,9] (Figure 1).The first part of the melt season is characterized by a dominant lateral meltwater flow, widespread pond coverage, and a rapid increase in MPF [8].Later, in the second part of the melt season, meltwater drainage occurs mainly by the lateral flow of meltwater across the ice surface, toward macroscopic flaws such as brine drainage channels, ice floe edges, and cracks.Then the snow cover melts completely, concurrent with a reduction in MPF and enhancement of albedo [8,9].In the third part of the melt season, sea ice has high permeability and open macroscopic flaws, which accelerate ice floe breakup and disintegration [8].The permeability of multiyear ice (MYI) open macroscopic flaws, which accelerate ice floe breakup and disintegration [8].The permeability of multiyear ice (MYI) is lower than that of first-year ice (FYI), and the low permeability of MYI helps to retain meltwater on the sea ice's surface [10,11].In the last part of the melt season, the surface of the melt pond freezes in cases where sea ice remains throughout the melt season [8].Meltwater drainage (particularly on MYI) does not always occur or slowly occurs, meaning that meltwater remains on the surface throughout the melt season, that there may be no sharp boundary for onset and offset, and that the meltwater surface eventually freezes in the fall freeze-up period (Figure 1b).The first and second parts of the melt season, characterized as having a very low albedo and a near peak of annual insolation, have a particularly large impact on the surface energy balance [9].Then, sea ice melting is promoted by a positive permeability/heat-flux feedback, resulting in abundant meltwater drainage into highly permeable ice, particularly FYI [8,9].Meltwater of low salinity and low nutrient concentrations is also retained under ice or in bottom depressions, forming under-ice freshwater lenses (i.e., under-ice melt ponds).The redistribution of meltwater affects the salt and heat balance at the layer where fresh meltwater and saline water mix [12], and the supply of atmospheric pollutants into the ocean and the sea-ice ecosystem [13].Understanding the drainage onset (DO) timing and drainage duration (DD) related to meltwater redistribution in the early melt season will help to model changes in physical and biological coupling within the Arctic sea ice ecosystem.
An understanding of DO and DD is often achieved through both airborne and in situ measurements [8][9][10][11]14], although this is difficult at larger scales for temporal and spatial limitations.Satellite passive microwave and optical remote sensing is useful for monitoring the entire polar regions daily [15,16], although optical remote sensing is limited by cloud cover when observing sea ice conditions.Although the timing of melt onset (MO) has been derived using passive microwave remote sensing [17][18][19][20], little is known about large-scale DO timing and DD in the melt season.
This study aims to derive this DO and DD using the time series of MPF retrieved data from the Advanced Microwave Scanning Radiometer-Earth observing system (AMSR-E MPF), an estimate proposed by Tanaka et al. [21] (see Supplementary Text S1).The AMSR-E MPF estimate is valid in the Beaufort Sea with pack ice and Canadian Arctic Archipelago (CAA) areas with landfast ice.To determine the DO and DD, only the CAA region with landfast ice was chosen for this study since the same sea ice needs to be observed during the melt season.To avoid the contamination of open water in the AMSR-E MPF estimate, we used a totally ice-covered area (10/10 ice concentration), without open water contamination in the AMSR-E MPF estimate until the breakup of landfast ice using ice charts, available from 2006 to present (Section 2.1).To address the lack of evaluation on a day-by-day basis for the AMSR-E MPF estimate (see Supplementary Text S2), we first compare the AMSR-E MPF to the daily MPF product from the Medium Resolution Imaging Spectrometer (MERIS MPF) and then discuss the biases related to melting snow in the early melt season.Results from the six-year period The first and second parts of the melt season, characterized as having a very low albedo and a near peak of annual insolation, have a particularly large impact on the surface energy balance [9].Then, sea ice melting is promoted by a positive permeability/heat-flux feedback, resulting in abundant meltwater drainage into highly permeable ice, particularly FYI [8,9].Meltwater of low salinity and low nutrient concentrations is also retained under ice or in bottom depressions, forming under-ice freshwater lenses (i.e., under-ice melt ponds).The redistribution of meltwater affects the salt and heat balance at the layer where fresh meltwater and saline water mix [12], and the supply of atmospheric pollutants into the ocean and the sea-ice ecosystem [13].Understanding the drainage onset (DO) timing and drainage duration (DD) related to meltwater redistribution in the early melt season will help to model changes in physical and biological coupling within the Arctic sea ice ecosystem.
An understanding of DO and DD is often achieved through both airborne and in situ measurements [8][9][10][11]14], although this is difficult at larger scales for temporal and spatial limitations.Satellite passive microwave and optical remote sensing is useful for monitoring the entire polar regions daily [15,16], although optical remote sensing is limited by cloud cover when observing sea ice conditions.Although the timing of melt onset (MO) has been derived using passive microwave remote sensing [17][18][19][20], little is known about large-scale DO timing and DD in the melt season.
This study aims to derive this DO and DD using the time series of MPF retrieved data from the Advanced Microwave Scanning Radiometer-Earth observing system (AMSR-E MPF), an estimate proposed by Tanaka et al. [21] (see Supplementary Text S1).The AMSR-E MPF estimate is valid in the Beaufort Sea with pack ice and Canadian Arctic Archipelago (CAA) areas with landfast ice.To determine the DO and DD, only the CAA region with landfast ice was chosen for this study since the same sea ice needs to be observed during the melt season.To avoid the contamination of open water in the AMSR-E MPF estimate, we used a totally ice-covered area (10/10 ice concentration), without open water contamination in the AMSR-E MPF estimate until the breakup of landfast ice using ice charts, available from 2006 to present (Section 2.1).To address the lack of evaluation on a day-by-day basis for the AMSR-E MPF estimate (see Supplementary Text S2), we first compare the AMSR-E MPF to the daily MPF product from the Medium Resolution Imaging Spectrometer (MERIS MPF) and then discuss the biases related to melting snow in the early melt season.Results from the six-year period of 2006-2011 are used for both evaluation and examination.Data from the AMSR2 sensor, which superseded the AMSR-E sensor, on the Global Change Observation Mission 1st-Water (GCOM-W1) satellite will extend the continuous record of the MPF, DO, and DD from 2012 to present, when the MPF retrieved from AMSR2 data is evaluated in a future paper.

Canadian Ice Service Regional Ice Charts
Figure 2 is a map of our study area, including the CAA and its subregions.This also shows the extent of AMSR-E data and the products of MPF and spectral sea ice albedo from MERIS data, as discussed in the next subsection.To obtain total ice concentration data in the CAA, we used the Canada Ice Service (CIS) weekly ice charts in Sea Ice Geo Referenced Information and Data (SIGRID)-3 format [22], including total sea ice concentration from 2006 to 2011 [23].The ice charts were created based on the Manual of Standard Procedures for Observing and Reporting Ice Conditions (MANICE) [24] and were prepared by following internationally recommended terminology and symbology established by the World Meteorological Organization (WMO).The ice charts from 2006 to the present were created manually by the CIS based on various satellite data, aerial reconnaissance, ship reports, and operational model results and then made available through the National Snow and Ice Data Center (NSIDC).These data are provided in six different regions, but this study used only charts from the Western and Eastern Arctic regions, as shown in Figure 2. SIGRID-3 data are in ESRI shapefile format, following the JCOMM convention [22].Supplementary Text S3 summarizes the main sources of error in the ice charts described by Tivy et al. [25].
superseded the AMSR-E sensor, on the Global Change Observation Mission 1st-Water (GCOM-W1) satellite will extend the continuous record of the MPF, DO, and DD from 2012 to present, when the MPF retrieved from AMSR2 data is evaluated in a future paper.

Canadian Ice Service Regional Ice Charts
Figure 2 is a map of our study area, including the CAA and its subregions.This also shows the extent of AMSR-E data and the products of MPF and spectral sea ice albedo from MERIS data, as discussed in the next subsection.To obtain total ice concentration data in the CAA, we used the Canada Ice Service (CIS) weekly ice charts in Sea Ice Geo Referenced Information and Data (SIGRID)-3 format [22], including total sea ice concentration from 2006 to 2011 [23].The ice charts were created based on the Manual of Standard Procedures for Observing and Reporting Ice Conditions (MANICE) [24] and were prepared by following internationally recommended terminology and symbology established by the World Meteorological Organization (WMO).The ice charts from 2006 to the present were created manually by the CIS based on various satellite data, aerial reconnaissance, ship reports, and operational model results and then made available through the National Snow and Ice Data Center (NSIDC).These data are provided in six different regions, but this study used only charts from the Western and Eastern Arctic regions, as shown in Figure 2. SIGRID-3 data are in ESRI shapefile format, following the JCOMM convention [22].Supplementary Text S3 summarizes the main sources of error in the ice charts described by Tivy et al. [25].
The ice charts used by Tanaka et al. [21] were paper maps, using areas with "9+/10" and "10/10".Each area is defined as having a sea ice concentration exceeding 95%, which was manually extracted.The ice charts in this study provide three sea ice concentration ranges of "9-10/10", "9+/10", and "10/10" as one code.To avoid the contamination of open water in the AMSR-E MPF estimate, we used only the totally ice-covered areas with landfast ice (10/10 sea ice concentration).The ice charts were resampled into 25 km × 25 km NSIDC stereographic projection grids (NSIDC grids).The ice charts used by Tanaka et al. [21] were paper maps, using areas with "9+/10" and "10/10".Each area is defined as having a sea ice concentration exceeding 95%, which was manually extracted.The ice charts in this study provide three sea ice concentration ranges of "9-10/10", "9+/10", and "10/10" as one code.To avoid the contamination of open water in the AMSR-E MPF estimate, we used only the totally ice-covered areas with landfast ice (10/10 sea ice concentration).The ice charts were resampled into 25 km × 25 km NSIDC stereographic projection grids (NSIDC grids).

Satellite Data Product
We used the daily brightness temperature (T B ) in the AMSR-E/Aqua Daily L3 25 km polar stereographic product provided by the NSIDC for the AMSR-E MPF estimate.The AMSR-E operated from 2002 to 2011 and was a passive microwave radiometer on NASA's Aqua satellite, equipped with 12 channels at 6.9, 10.7, 18.7, 23.8, 36.5, and 89.0 GHz, including both vertical (V) and horizontal (H) polarizations.It provided data in both ascending and descending (morning and evening) orbits.The spatial resolutions of individual frequencies ranged from 4 km × 6 km at 89.0 GHz to 43 km × 75 km at 6.9 GHz.In the product, the resampling interval is 25 km × 25 km at all frequencies.We used the daily average, T B , and only the H polarization channel at 6.9 GHz and the V polarization channel at 89 GHz.The input data for the latest version 3 T B product are the AMSR-E L2A Global Swath Spatially Resampled Brightness Temperatures product, including land contamination.The contamination impact on T B (particularly at 6.9 GHz) near shore results from the fact that T B is higher over land than sea.Tanaka et al. [21] found that an increase in the T B at 6.9 GHz means a decrease in the MPF.Therefore, land contamination may lead to a lower MPF near the shore.We used the DMSP SSM/I-SSMIS Daily Polar Gridded Brightness Temperatures Version 5 product provided by the NSIDC, to determine the freeze onset (FO) date.We used only the T B s from 2006 to 2011, with V polarization at 19 and 37 GHz.
For evaluating the AMSR-E MPF from 2006 to 2011, we used the evaluation product, which consists of daily averaged MPFs, and broadband albedo, using the 12.5 km NSIDC grid with the criterion of more than 50% valid pixels (both spatially and temporally).This period corresponds to the availability period of the ice chart data (Section 2.1).The product was made based on the swath MERIS Level 1b data with 1 km × 1 km spatial resolution over the Arctic summers from 2002 to 2011.The MPF data were obtained by the melt pond detector (MPD) algorithm [16].The 12.5 km NSIDC grid is projected onto the 25 km NSIDC grid, which averages by using four of the 12.5 km NSIDC grid cells.All comparisons were made using the 25 km NSIDC grid.Supplementary Text S4 summarizes the characteristics of the products of MPF and spectral sea ice albedo from MERIS data [26,27].

Estimating Meltwater Drainage Onset Timing and Duration
As an example, we show how to determine the DO timing and DD by using the six points plotted in Figure 2. Figure 3 shows the time series of the AMSR-E MPF and daily average 2 m air temperature in the ERA5 dataset provided by the European Centre for Medium-Range Weather Forecasts (ECMWF, https://www.ecmwf.int/).The MO date is obtained by using the dynamic threshold variability method based on 36.5 GHz vertically polarized T B data [20].The FO date is retrieved from the melt/freeze-up algorithm, using 19 and 37 GHz vertically polarized T B data [18].As previously stated, the MPF after MO rapidly increases, and later, when the MPF decreases temporally, is when an increase in ice permeability or meltwater drainage occurs (Figure 1).This feature was also found in the time series of the AMSR-E MPF (black circles in Figure 3d-f) and then used to estimate the most likely DO timing and DD. Figure 4 is a flowchart of the method for estimating DO timing and DD, relative to the following five steps: (1) The time series of the MPF is created from the AMSR-E MPF estimation method.The time series in Figure 3 includes MPF estimates outside the conditions when melt ponds are expected to occur, i.e., prior to MO: the gray circles in Figure 3 after the FYI completely melts or breaks up and there is open water (gray x's in Figure 3), and after the FO.These data points outside the possible melt season are problematic and are not to be included; we only used the realistic time series over 10/10 ice concentration and between the MO and FO dates, as the black circles in Figure 3 show.If the threshold of the MO is not determined due to the high inter-quartile range in the histogram of the potential MO [20], a data point of more than 0 • C air temperature is used for determining a realistic time series.
(2) The time series can be expressed accurately as a third or fourth order function, as shown in Figure 1.We assume a third or fourth (k = 3 or 4) order polynomial regression in one variable, (4) The maximum and minimum values are detected on the k = 3 or 4 order regression curves.In [8,9], the first maximum was taken as the DO date, and the first minimum was taken as the end-of-drainage (ED) date (i.e., the onset of melt pond evolution).Following these, we defined the first maximum as the DO and the first minimum as the ED posterior to the maximum, on the curves.If there are no maximum values, the dates are not estimated for that cell (Case 2).
(5) Finally, DD is calculated from the difference between ED and DO. is an explanatory variable.To examine the significance of b1 to bk, we determine pk-values for all coefficients of the k = 3 or 4 (p3 and p4, respectively).If p3 > 0.05 and p4 > 0.05, the grid cell is defined as an unfitted case (Case 1).
(3) The time series is fitted with k = 3 if p3 < 0.05 and p4 > 0.05, and with k = 4 if p3 > 0.05 and p4 < 0.05.If both p3 and p4 are less than 0.05, k is determined by using the degree of freedom adjusted coefficient of determination (Rk).
(4) The maximum and minimum values are detected on the k = 3 or 4 order regression curves.In [8,9], the first maximum was taken as the DO date, and the first minimum was taken as the endof-drainage (ED) date (i.e., the onset of melt pond evolution).Following these, we defined the first maximum as the DO and the first minimum as the ED posterior to the maximum, on the curves.If there are no maximum values, the dates are not estimated for that cell (Case 2).
(5) Finally, DD is calculated from the difference between ED and DO.   Figure 5 shows a histogram and the spatial distributions for each case.Case 1 is 20%-35% of all grid cells and is distributed in the West Arctic Waterway, Baffin Inlet, and East Parry Channel.This is because a passive microwave signal has difficulty distinguishing between saline seawater and freshwater in ponds in cold-water regions when the first-year ice is completely melting or breaking up (Figure 3a) and low significance of the k associated with high variability of the time series (Figure 3b).Case 2 is around 10%.Case 3 (k = 3 regression) is 20%-40%, and the realistic time series is broken by contamination of open water after the ED timing.Case 4 (k = 4 regression) is 30%-55%.The time series in Figure 3e (similar to that in Figure 3d) incorporates k = 4 regression.This case needs a higher k because the variability of the time series is higher for Figure 3e than Figure 3d.The difference in variability between Case 3 and Case 4 areas may be due to melting snow coverage or land contamination on the TB (particularly for the coarse spatial resolution at 6.9 GHz).  Figure 5 shows a histogram and the spatial distributions for each case.Case 1 is 20%-35% of all grid cells and is distributed in the West Arctic Waterway, Baffin Inlet, and East Parry Channel.This is because a passive microwave signal has difficulty distinguishing between saline seawater and freshwater in ponds in cold-water regions when the first-year ice is completely melting or breaking up (Figure 3a) and low significance of the k associated with high variability of the time series (Figure 3b).Case 2 is around 10%.Case 3 (k = 3 regression) is 20%-40%, and the realistic time series is broken by contamination of open water after the ED timing.Case 4 (k = 4 regression) is 30%-55%.The time series in Figure 3e (similar to that in Figure 3d) incorporates k = 4 regression.This case needs a higher k because the variability of the time series is higher for Figure 3e than Figure 3d.The difference in variability between Case 3 and Case 4 areas may be due to melting snow coverage or land contamination on the T B (particularly for the coarse spatial resolution at 6.9 GHz).  Figure 5 shows a histogram and the spatial distributions for each case.Case 1 is 20%-35% of all grid cells and is distributed in the West Arctic Waterway, Baffin Inlet, and East Parry Channel.This is because a passive microwave signal has difficulty distinguishing between saline seawater and freshwater in ponds in cold-water regions when the first-year ice is completely melting or breaking up (Figure 3a) and low significance of the k associated with high variability of the time series (Figure 3b).Case 2 is around 10%.Case 3 (k = 3 regression) is 20%-40%, and the realistic time series is broken by contamination of open water after the ED timing.Case 4 (k = 4 regression) is 30%-55%.The time series in Figure 3e (similar to that in Figure 3d) incorporates k = 4 regression.This case needs a higher k because the variability of the time series is higher for Figure 3e than Figure 3d.The difference in variability between Case 3 and Case 4 areas may be due to melting snow coverage or land contamination on the TB (particularly for the coarse spatial resolution at 6.9 GHz).

Results and Discussion
Figure 6a-f shows the relationship between the MERIS and AMSR-E MPFs from 2006 to 2011.We show only the relationship over 10/10 ice concentration and from the MO date to the FO date.The data of AMSR-E MPFs strongly correlated with those of MERIS MPFs with r = 0.73−0.83and a root mean square error (RMSE) of 8-9%.The AMSR-E MPF estimates in areas less than 10%-20% and ranging from 20% to 50% overestimate and underestimate the MERIS MPF estimates, respectively.The lower MPFs were biased toward the positive, and the higher MPFs were biased toward the negative, meaning that the negative and positive biases largely cancel each other out.This leads to the low bias of ±2% from 2006 to 2011.The biases strongly correlated with an average albedo with r = 0.68−0.96at the 99% confidence level (Figure 6g).

Results and Discussion
Figure 6a-f shows the relationship between the MERIS and AMSR-E MPFs from 2006 to 2011.We show only the relationship over 10/10 ice concentration and from the MO date to the FO date.The data of AMSR-E MPFs strongly correlated with those of MERIS MPFs with r = 0.73−0.83and a root mean square error (RMSE) of 8-9%.The AMSR-E MPF estimates in areas less than 10%-20% and ranging from 20% to 50% overestimate and underestimate the MERIS MPF estimates, respectively.The lower MPFs were biased toward the positive, and the higher MPFs were biased toward the negative, meaning that the negative and positive biases largely cancel each other out.This leads to the low bias of ±2% from 2006 to 2011.The biases strongly correlated with an average albedo with r = 0.68−0.96at the 99% confidence level (Figure 6g).Microwave radiometric properties are influenced by sea ice condition between the melt ponds and by melting snow outside melt ponds (e.g., prior to).We now discuss this effect on the AMSR-E MPF estimate.Different ice conditions during the melt season are often summarized as broadband albedos [28][29][30], so we determined three categories according to the broadband albedo obtained from the products of MPF and spectral sea ice albedo from MERIS data, as follows: melting snow between the ponds or melting snow outside ponds for an albedo > 0.7; bare ice including light ponds for an albedo of 0.6-0.7;and ponded ice for an albedo of 0.3-0.6.Figure 6g shows the bias at each albedo with 0.02 bin size from 2006 to 2011.The areas over 10/10 ice concentration and from the MO to FO Microwave radiometric properties are influenced by sea ice condition between the melt ponds and by melting snow outside melt ponds (e.g., prior to).We now discuss this effect on the AMSR-E MPF estimate.Different ice conditions during the melt season are often summarized as broadband albedos [28][29][30], so we determined three categories according to the broadband albedo obtained from the products of MPF and spectral sea ice albedo from MERIS data, as follows: melting snow between the ponds or melting snow outside ponds for an albedo > 0.7; bare ice including light ponds for an albedo of 0.6-0.7;and ponded ice for an albedo of 0.3-0.6.Figure 6g shows the bias at each albedo with 0.02 bin size from 2006 to 2011.The areas over 10/10 ice concentration and from the MO to FO dates are used to obtain the average bias at each albedo.The biases for an albedo between 0.4 and 0.5 were negative and close to zero around albedo of 0.5.At an albedo between 0.5 and 0.7, the biases were around zero and constant, and they shifted from negative to positive at an albedo over 0.7.The largest positive/negative biases were at an albedo of 0.7-0.9/around0.4.We found that the AMSR-E MPFs in areas with an albedo of 0.7-0.9(including melting snow) are 10% larger than the MERIS MPFs, and that bare ice shows little bias compared to the MERIS MPF.T B decreases at 6.9 GHz and increases at 89.0 GHz, with an increase in the MPF, meaning that the gradient ratio GR 06H-89V decreases with an increase in the MPF [21].When dry snow starts to melt, the T B also decreases around 6 GHz and increases around 89 GHz.Therefore, we consider that an increase of snow wetness on sea ice probably led to a positive bias (an increase of pseudo-AMSR-E MPF).After all the snow was gone, there were few biases in the area with an albedo of 0.5-0.7.This suggests that the AMSR-E MPF product provides a quality product, similar to the MERIS MPF, only when the albedo is around 0.5-0.7.The negative biases to 20% were apparent around an albedo of 0.4-0.5.The AMSR-E MPF estimate is mostly useful after snow cover has melted completely.In contrast, the reflectance of wet ice is very close to that of a melt pond [16], so this too may lead to an increase of pseudo-MERIS MPF.We can conclude that an increase in wet ice leads to negative biases.
The DOs/DDs from the time series of the AMSR-E and MERIS MPFs are compared.Here we used 1954 datapoints on a number of AMSR-E and co-located MERIS grid cells.Figure 7 shows the differences between DOs/DDs from the time series of the AMSR-E and MERIS MPFs.These differences show a normal distribution.The mode and average of the DO/DD differences were 2.0/-2.0days and 2.1/-3.4days, respectively, with a standard deviation of 6.6/6.0 days, indicating that DO and DD are consistent.dates are used to obtain the average bias at each albedo.The biases for an albedo between 0.4 and 0.5 were negative and close to zero around albedo of 0.5.At an albedo between 0.5 and 0.7, the biases were around zero and constant, and they shifted from negative to positive at an albedo over 0.7.The largest positive/negative biases were at an albedo of 0.7-0.9/around0.4.We found that the AMSR-E MPFs in areas with an albedo of 0.7-0.9(including melting snow) are 10% larger than the MERIS MPFs, and that bare ice shows little bias compared to the MERIS MPF.TB decreases at 6.9 GHz and increases at 89.0 GHz, with an increase in the MPF, meaning that the gradient ratio GR06H-89V decreases with an increase in the MPF [21].When dry snow starts to melt, the TB also decreases around 6 GHz and increases around 89 GHz [31].Therefore, we consider that an increase of snow wetness on sea ice probably led to a positive bias (an increase of pseudo-AMSR-E MPF).After all the snow was gone, there were few biases in the area with an albedo of 0.5-0.7.This suggests that the AMSR-E MPF product provides a quality product, similar to the MERIS MPF, only when the albedo is around 0.5-0.7.The negative biases to 20% were apparent around an albedo of 0.4-0.5.The AMSR-E MPF estimate is mostly useful after snow cover has melted completely.In contrast, the reflectance of wet ice is very close to that of a melt pond [16], so this too may lead to an increase of pseudo-MERIS MPF.
We can conclude that an increase in wet ice leads to negative biases.The DOs/DDs from the time series of the AMSR-E and MERIS MPFs are compared.Here we used 1954 datapoints on a number of AMSR-E and co-located MERIS grid cells.Figure 7 shows the differences between DOs/DDs from the time series of the AMSR-E and MERIS MPFs.These differences show a normal distribution.The mode and average of the DO/DD differences were 2.0/-2.0days and 2.1/-3.4days, respectively, with a standard deviation of 6.6/6.0 days, indicating that DO and DD are consistent.We show maps of the MO timing (Supplementary Figure S1), DO timing (Figure 8), and DD (Figure 9), from 2006 to 2011.Results from the CAA and seven subregions are summarized by averaging the 2006-2011 data (Supplementary Table S1).The average date of DO in the CAA was 177 DOY, with a standard deviation of 9 days.On average from 2006 to 2011, the earliest and the latest DO dates were 171 DOY in the West Arctic Waterway and 182 DOY in the QEI, respectively (Figure 8a).The standard deviations of DO timing were higher in the QEI and West Parry Channel, and lower in other subregions (Figure 8b).The determinations of DO timing and DD in the East Parry Channel were limited by the contamination of open water in the AMSR-E MPF estimate from earlier melt seasons.In each year, contrasts in the DO timing between the QEI and the other six regions in 2006, 2009, and 2010 (Figure 8c,f,g) were higher than those in 2007, 2008, and 2011 (Figure 8d,e,h).The characteristics of both high and low contrasts were found in the map of MO timing (Supplementary Figure S1).We compared the averaged MO timings to the averaged DO timings.The value of r for MO and DO indicated a strong relationship (r = 0.77, p-value < 0.001).The scatterplot for the East Parry Channel in the relationship between the averaged MO timing and the averaged DO timing was We show maps of the MO timing (Supplementary Figure S1), DO timing (Figure 8), and DD (Figure 9), from 2006 to 2011.Results from the CAA and seven subregions are summarized by averaging the 2006-2011 data (Supplementary Table S1).The average date of DO in the CAA was 177 DOY, with a standard deviation of 9 days.On average from 2006 to 2011, the earliest and the latest DO dates were 171 DOY in the West Arctic Waterway and 182 DOY in the QEI, respectively (Figure 8a).The standard deviations of DO timing were higher in the QEI and West Parry Channel, and lower in other subregions (Figure 8b).According to [8][9][10][11], FYI and MYI have higher and lower permeabilities, respectively, and active horizontal flow, particularly in level ice, during low ice permeability, leads to meltwater into macroscopic pathways, such as seal breathing holes and cracks.Therefore, the variability of DO timing and DD in each year may be associated with sea ice type and roughness.To explain the relationship between the variability and sea ice type or roughness, time series of AMSR2 MPF from 2012 to present is required for long-term records of the MPF, DO timing, and DD.

Conclusion
This study leads to following results: (1) We evaluated the AMSR-E MPF estimate [21] by using daily average TB for the AMSR-E L3

Figure 1 .
Figure 1.Conceptual sketch of the temporal evolution of melt pond fraction (MPF) (solid curve) and meltwater loss at macroscopic holes (dashed curve) in.(a) Drainage onset (DO) is determined as the first maximum.End of drainage (ED) is determined as the first minimum in (b).

Figure 1 .
Figure 1.Conceptual sketch of the temporal evolution of melt pond fraction (MPF) (solid curve) and meltwater loss at macroscopic holes (dashed curve) in.(a) Drainage onset (DO) is determined as the first maximum.End of drainage (ED) is determined as the first minimum in (b).

Figure 2 .
Figure 2. Locations of the Canadian Arctic Archipelago and its subregions.Triangle, black x, red x, red circle, red square, red hexagon, and black diamond show monitoring points in 2007, as shown in Figure 3.

Figure 2 .
Figure 2. Locations of the Canadian Arctic Archipelago and its subregions.Triangle, black x, red x, red circle, red square, red hexagon, and black diamond show monitoring points in 2007, as shown in Figure 3.
where y (MPF) is the predicted outcome of the polynomial regression.The regression coefficients are b 1 to b k for each order, where b 0 is the y-intercept, and x (day of year, DOY) is an explanatory variable.To examine the significance of b 1 to b k , we determine p k -values for all coefficients of the k = 3 or 4 (p 3 and p 4 , respectively).If p 3 > 0.05 and p 4 > 0.05, the grid cell is defined as an unfitted case (Case 1).(3)The time series is fitted with k = 3 if p 3 < 0.05 and p 4 > 0.05, and with k = 4 if p 3 > 0.05 and p 4 < 0.05.If both p 3 and p 4 are less than 0.05, k is determined by using the degree of freedom adjusted coefficient of determination (R k ).

Figure 3 .
Figure 3.Time series of the Advanced Microwave Scanning Radiometer-Earth observing system (AMSR-E) melt pond fractions (MPFs) and 2 m air temperature from the ERA5 dataset.This figure shows six sets of monitoring points for 2007 in (a) Case 1 (Baffin Inlet), (b) Case 1 (QEI), (c) Case 2 (Baffin Inlet), (d) Case 3 (West Arctic Waterway), (e) Case 4 (Baffin Inlet), and (f) Case 4 (QEI), shown in Figure 2. Black dots are datapoints for estimating drainage onset (DO) and drainage duration (DD).Gray circles and x's are before melt onset (MO) and include open water contamination under 9/10 sea ice concentration, respectively.Red curves show k-order polynomial regression.Solid gray curves show the 2 m air temperature.Vertical dashed lines show MO (back), DO (green), and DD (blue).Vertical solid black lines show freeze onset.Horizontal solid black lines show 0 o C air temperature.

Figure 3 .
Figure 3.Time series of the Advanced Microwave Scanning Radiometer-Earth observing system (AMSR-E) melt pond fractions (MPFs) and 2 m air temperature from the ERA5 dataset.This figure shows six sets of monitoring points for 2007 in (a) Case 1 (Baffin Inlet), (b) Case 1 (QEI), (c) Case 2 (Baffin Inlet), (d) Case 3 (West Arctic Waterway), (e) Case 4 (Baffin Inlet), and (f) Case 4 (QEI), shown in Figure 2. Black dots are datapoints for estimating drainage onset (DO) and drainage duration (DD).Gray circles and x's are before melt onset (MO) and include open water contamination under 9/10 sea ice concentration, respectively.Red curves show k-order polynomial regression.Solid gray curves show the 2 m air temperature.Vertical dashed lines show MO (back), DO (green), and DD (blue).Vertical solid black lines show freeze onset.Horizontal solid black lines show 0 • C air temperature.

Figure 4 .
Figure 4. Flowchart of determination of drainage onset (DO), end of drainage (ED), and drainage duration (DD) estimates.

Figure 5 .
Figure 5. (a) Histogram and (b) spatial distributions for each case in the Canadian Arctic Archipelago (CAA), from 2006 to 2011.

Figure 4 .
Figure 4. Flowchart of determination of drainage onset (DO), end of drainage (ED), and drainage duration (DD) estimates.

Figure 4 .
Figure 4. Flowchart of determination of drainage onset (DO), end of drainage (ED), and drainage duration (DD) estimates.

Figure 5 .
Figure 5. (a) Histogram and (b) spatial distributions for each case in the Canadian Arctic Archipelago (CAA), from 2006 to 2011.

Figure 5 .
Figure 5. (a) Histogram and (b) spatial distributions for each case in the Canadian Arctic Archipelago (CAA), from 2006 to 2011.

Figure 6 .
Figure 6.(a-f) Relationship between the Medium Resolution Imaging Spectrometer (MERIS) and AMSR-E MPFs.The red line: linear regression, color bar: number of grid cells (N), r: correlation coefficient, RMSE: root mean square error.(g) Average bias between AMSR-E and MERIS MPFs at each albedo with 0.02 bin size.The albedo was obtained from the products of MPF and spectral sea ice albedo from MERIS data.The error bars are standard deviation of the bias.

Figure 6 .
Figure 6.(a-f) Relationship between the Medium Resolution Imaging Spectrometer (MERIS) and AMSR-E MPFs.The red line: linear regression, color bar: number of grid cells (N), r: correlation coefficient, RMSE: root mean square error.(g) Average bias between AMSR-E and MERIS MPFs at each albedo with 0.02 bin size.The albedo was obtained from the products of MPF and spectral sea ice albedo from MERIS data.The error bars are standard deviation of the bias.

Figure 7 .
Figure 7. DO and DD differences between AMSR-E and MERIS.SD: standard deviation.
The determinations of DO timing and DD in the East Parry Channel were limited by the contamination of open water in the AMSR-E MPF estimate from earlier melt seasons.In each year, contrasts in the DO timing between the QEI and the other six regions in 2006, 2009, and 2010 (Figure 8c,f,g) were higher than those in 2007, 2008, and 2011 (Figure 8d,e,h).The characteristics of both high and low contrasts were found in the map of MO timing (Supplementary Figure S1).We compared the averaged MO timings to the averaged DO timings.The value of r for MO and DO indicated a strong relationship (r = 0.77, p-value < 0.001).The scatterplot for the East Parry Channel in the relationship between the averaged MO timing and the averaged DO timing was far from the linear regression line (not shown).This may limit determining the DO timing.The correlation increases to r = 0.91 (p-value < 0.001) without any points in the East Parry Channel, so that the average DO timing follows the average MO timing in the subregions.The determination of DDs in the East Parry Channel and Baffin Inlet was also limited by the contamination of open water.We exclude data from the East Parry Channel and Baffin Inlet.The average DD in the CAA was 13 days with a standard deviation of nine days (Figure 9a,b).The shortest and longest DDs were 11 days in the West Parry Channel and 17 days in the QEI, respectively.Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 12 far from the linear regression line (not shown).This may limit determining the DO timing.The correlation increases to r = 0.91 (p-value < 0.001) without any points in the East Parry Channel, so that the average DO timing follows the average MO timing in the subregions.The determination of DDs in the East Parry Channel and Baffin Inlet was also limited by the contamination of open water.We exclude data from the East Parry Channel and Baffin Inlet.The average DD in the CAA was 13 days with a standard deviation of nine days (Figure 9a,b).The shortest and longest DDs were 11 days in the West Parry Channel and 17 days in the QEI, respectively.

Figure 8 .
Figure 8. Spatial distributions of drainage onset (DO) in the Canadian Arctic Archipelago, from 2006 to 2011.Black points are Case 2.

Figure 8 . 12 Figure 9 .
Figure 8. Spatial distributions of drainage onset (DO) in the Canadian Arctic Archipelago, from 2006 to 2011.Black points are Case 2. Remote Sens. 11, x FOR PEER REVIEW 10 of 12