Satellite Microwave Assessment of Northern Hemisphere Lake Ice Phenology from 2002 to 2015

A new automated method for satellite assessment of seasonal lake ice phenology at 5-km resolution was developed for all lake pixels (water coverage ≥ 90%) in the Northern Hemisphere using 36.5 GHz, H-polarized brightness 10 temperature (Tb) observations from the Advanced Microwave Scanning Radiometer (AMSR-E/2) sensors. The lake phenology metrics include seasonal timing and duration of annual ice cover. A Moving t-Test (MTT) algorithm allows for automated lake ice retrievals with daily temporal fidelity and 5-km resolution gridding. The resulting ice phenology record shows strong agreement with available ground-based observations from the Global Lake and River Ice Phenology Database (95.4% temporal agreement), and favourable correlations (R) with alternative ice phenology records from the Interactive 15 Multisensor Snow and Ice Mapping System (R = 0.84 for water clear of ice [WCI] dates; R = 0.41 for complete freeze over [CFO] dates) and Canadian Ice Service (R = 0.86 for WCI dates; R = 0.69 for CFO dates). Analysis of the resulting 12-year (2002-2015) AMSR ice record indicates increasingly shorter ice-cover duration for 43 out of 71 (60.6%) Northern Hemisphere lakes examined, with significant (p < 0.05) regional trends toward earlier ice melting for only five lakes. Higher latitude lakes reveal more widespread and larger trends toward shorter ice cover duration than lower latitude lakes, consistent 20 with enhanced polar warming. This study documents a new satellite-based approach for rapid assessment and regional monitoring of seasonal ice cover changes over large lakes, with resulting accuracy suitable for global change studies.


Introduction
Ice phenology describes the seasonal cycle of lake ice cover and encompasses freeze-up and breakup periods and ice cover duration (Duguay et al., 2015a).Freeze-up corresponds to the time period between the beginning of ice formation and the formation of a complete sheet of ice; breakup involves the time period between the onset of spring melt and the complete disappearance of ice from the lake surface (Kang et al., 2012).These ice phenology variables are key metrics sensitive to weather and climate conditions and influence lakeatmosphere interactions and hydrological and ecological processes in high-latitude and high-altitude regions (Duguay et al., 2006(Duguay et al., , 2012(Duguay et al., , 2015a;;Mishra et al., 2011).By insulating lake water from the overlying atmosphere and minimizing water and atmosphere heat and gas exchanges, lake ice has a controlling influence on water-column oxygen concentration, water temperature, and the composition and abundance of aquatic species (Livingstone, 1997;Bengtsson and Herschy, 2012;Kang et al., 2012;Wrona et al., 2016).In addition to the impacts on aquatic life, the formation and disappearance of lake ice also has a significant influence on the spread of man-made pollutants such as perfluorinated chemicals (PFCs) (Veillette et al., 2012;Wrona et al., 2016).The extent and duration of lake ice also affect human activities, including hydroelectric power generation, navigation Published by Copernicus Publications on behalf of the European Geosciences Union.and winter transportation, and production and distribution of food and water (Schröter et al., 2005;Weyhenmeyer et al., 2011).Moreover, lake ice phenology is closely coupled with atmospheric heat fluxes (Latifovic and Pouliot, 2007;Park et al., 2016) and sensitive to the alteration of weather patterns under projected global warming (Magnuson et al., 2000).
Accurate and consistent records of lake ice phenology especially in data-sparse regions, including much of the pan-Arctic and Qinghai-Tibetan Plateau regions, provide valuable information for monitoring global change impacts on high-latitude and high-altitude environments (Magnuson et al., 2000).Previous studies have documented significantly earlier ice breakup between the 1950s and 2000s for lakes in Canada (Duguay et al., 2006;Latifovic and Pouliot, 2007;Prowse et al., 2011) and decreasing ice cover duration of Eurasia lakes during the last few decades (Vuglinsky and Gronskaya, 2006;Karetnikov and Naumenko, 2008;Prowse et al., 2011).Shorter ice cover seasons may promote greater CH 4 emissions from northern lakes (Greene et al., 2014), which could reinforce further climate warming due to the role of CH 4 as a potent greenhouse gas.Despite a general tendency for later freezing and earlier breakup in the Northern Hemisphere (Magnuson et al., 2000), various tendencies including earlier ice formation and later ice breakup over specific lakes and time periods may exist.For example, observations from satellite altimetry and radiometry over 1992-2004 for Lake Baikal showed a tendency for colder winters, with earlier ice formation, later ice breakup, and ice duration increase (Kouraev et al., 2007a, b).
A historic lake ice phenology database has been assembled from long-term ground-based observations across the northern domain (Magnuson et al., 2000;Benson et al., 2012); however, the number of monitoring sites is extremely sparse, with variable observational recording periods and methods, which limits capabilities for regional assessment and monitoring of environmental changes.Data acquired from spaceborne optical-thermal infrared (TIR) and microwave sensors have also been applied for monitoring river and lake ice phenology at regional and global scales (Chu and Lindenschmidt, 2016).Optical-TIR remote sensing can provide accurate estimation of land surface temperature (LST) and classification of land cover types at relatively fine spatial resolutions (∼ 10s to 100s of meters), while LST (1 km resolution) and snow cover (500 m resolution) products derived from MODIS (Moderate Resolution Imaging Spectroradiometer) have been used to infer lake ice conditions (Nonaka et al., 2007;Hall et al., 2010;Kheyrollah Pour et al., 2012).Time series of AVHRR (Advanced Very High Resolution Radiometer) imagery have also been used to classify Canadian lake ice phenology events with relatively high accuracy (Latifovic and Pouliot, 2007).However, regional monitoring of lake ice dynamics from satellite optical-TIR sensors is constrained by signal degradation and data loss stemming from seasonal reductions in solar illumination at higher latitudes and persistent cloud cover, smoke, and other atmospheric aerosol contamination (Maslanik and Barry, 1987;Jeffries et al., 2005;Helfrich et al., 2007;Kang et al., 2012).
Satellite microwave remote sensing at lower frequencies (∼ < 89 GHz) is relatively insensitive to solar illumination and atmosphere constraints, while current microwave radiometers onboard polar orbiting satellites provide frequent (∼ daily) observations spanning northern (≥ 45 • N) land areas.The active and passive microwave retrievals are also highly sensitive to the large contrast in surface dielectric properties between open water and ice cover over large lakes.Despite successful applications using active microwave remote sensing in lake ice retrieval (Leconte and Klassen, 1991;Nolan et al., 2002;Howell et al., 2009;Geldsetzer et al., 2010), capabilities for global lake ice monitoring from satellite radar sensors have been constrained by limited global coverage and temporal frequency of observations.Alternatively, spaceborne microwave radiometers have provided brightness temperature (T b ) observations since 1978 with relatively high temporal fidelity (∼ 1-2 days) especially at higher (≥ 45 • N) latitudes.Frequent microwave radiometer data acquisition and complete time series of images are valuable to ice phenology studies and also supportive to improving numerical weather model predictions (Helfrich et al., 2007) and timely monitoring of lake ice events, including transient ice disturbances (Jeffries et al., 2005).Despite having relatively coarser spatial resolution retrievals than optical-TIR sensors, the capability for consistent daily lake ice monitoring available from passive microwave observations provides added precision for delineating lake ice phenology trends, which may be much smaller than year-to-year ice cover variability.The satellite T b retrievals are capable of detecting lake ice phenology events coinciding with large changes in surface emissivity, but the passive microwave retrievals are constrained by a generally coarser spatial resolution than radar and optical-TIR sensors.Despite these limitations, ice freeze-up and breakup events for Great Slave Lake (GSL) were monitored using a threshold-based method for SSM/I (Special Sensor Microwave Imager) observations at 85 GHz (Walker and Davey, 1993;Ménard et al., 2002).Recently, H-polarized AMSR-E (Advanced Microwave Scanning Radiometer for EOS) T b observations at 19 GHz were analyzed to determine ice phenology for GSL and Great Bear Lake (GBL), the two largest lakes in northern Canada (Kang et al., 2012(Kang et al., , 2014)).Similar T b records from SSM/I and SMMR (Scanning Multichannel Microwave Radiometer) were also used to monitor lake ice phenology for Nam Co Lake (Ke et al., 2013) and Qinghai Lake (Che et al., 2009) within the high-elevation Qinghai-Tibetan Plateau.Previous studies based on satellite passive microwave remote sensing have mainly focused on one or two lakes using empirical algorithm thresholds developed for specific study areas.There is a great need to develop a universal lake ice detection algorithm and establish a consistent ice phenology database covering major lakes at the global scale for climate impact assessments such as those published by the Intergovernmental Panel and Climate Change (Stocker et al., 2013).
In this study, we present a new automated method to derive lake ice phenology using 36.5 GHz H-polarized satellite radiometric T b measurements from AMSR-E and AMSR2 (Advanced Microwave Scanning Radiometer 2).The algorithm is used to produce daily lake ice maps with 5 km gridding.The resulting AMSR-E/2 Lake Ice Phenology (LIP) record encompasses all 5 km by 5 km lake pixels (water coverage ≥ 90 %) within the Northern Hemisphere (≥ 0 • N) and spans more than 12 years of observations including both AMSR-E (June 2002 to September 2011) and AMSR2 (June 2012 to December 2015) satellite sensor records.Here we present a detailed methods description and evaluation of the LIP product against other independent observations and alternative lake ice products.A trend analysis is also conducted to characterize recent regional LIP changes over the study period.

Study domain
This study utilizes satellite passive microwave remote sensing to detect lake ice changes for 5 km lake pixels in the Northern Hemisphere, with a particular focus on lake ice phenology in the mid-and high latitudes (≥ 30 • N).The domain (Fig. 1) includes the high northern pan-Arctic region and high-altitude Qinghai-Tibetan Plateau, which are data-sparse but strongly sensitive to the Arctic amplification effect (Serreze and Francis, 2006;Woo et al., 2007) and/or elevation-dependent warming (Wang et al., 2011;Mountain Research Initiative EDW Working Group, 2015).Both regions are also characterized by cold climate conditions with extensive winter ice cover.The resulting domain includes three sets of lakes for algorithm evaluation and lake ice phenology analysis.Among the lakes analyzed, four are represented in the Global Lake and River Ice Phenology Database (GLRIPD) (Benson and Magnuson, 2000) and were used to evaluate the LIP estimates on a per pixel basis against available ground-based observations; the four GLRIPD lakes evaluated include Lake Superior in the USA and Canada and Lake Oulujarvi, Lake Haukivesi, and Lake Paijanne in Finland.In addition, 12 North American lakes (GBL, GSL, Smallwood Lake, Nettiling Lake, Dubawnt Lake, Amadjuak Lake, Wollaston Lake, Baker Lake, Kasba Lake, Lesser Slave Lake, and Peter Pond Lake in Canada and Red Lake in USA) that experience annual breakup and freeze-up events were selected for lake-wide intercomparisons between the LIP metrics derived from this study and alternative lake ice products from the Interactive Multisensor Snow and Ice Mapping System (IMS) (Helfrich et al., 2007; http://www.natice.noaa.gov/ims/)and the Canadian Ice Service (CIS) (Howell et al., 2009).Finally, regional LIP trends were assessed over the 12-year (2002-2015) satellite record for 71 Northern Hemisphere lakes identified in the Global Lakes and Wetlands Database (GLWD) (Lehner and Döll, 2004).The lakes selected represent approximately 23 % (297 044 km 2 ) of the estimated total surface area of large lakes (area ≥ 50 km 2 ) within the domain (Lehner and Döll, 2004).The 71 lakes were selected on the basis of having (a) at least one 5 km lake pixel with 100 % water coverage and located outside of a 5 km land buffer zone and (b) all pixels representing the lake having at least 20 days with full ice coverage and 20 days with open water.Criterion (a) was used to reduce potential contamination from adjacent land areas since the native resolutions of 36.5 GHz observations are approximately 14 km × 8 km for AMSR-E and 12 km × 7 km for AMSR2, respectively (Imaoka et al., 2010).Criterion (b) emphasizes lakes having extended ice and open water seasons rather than those with a temporary ice cover or short open water phase.The 20-day minimum duration was set according to the predefined subsample sizes of our algorithm (Sect.2.3).

Datasets used for algorithm development
The lake ice detection algorithm developed in this study relies primarily on 36.5 GHz H-polarized T b retrievals from AMSR-E and AMSR2.The AMSR-E sensor was operational on the NASA Aqua satellite from June 2002 to October 2011 and provided twice-daily measurements of global microwave emissions over land with descending/ascending orbital equatorial crossings at 01:30/13:30 local time and vertically (V) and horizontally (H) polarized T b retrievals at six frequencies (6.9,10.7,18.7,23.8,36.5,89.0 GHz).For this study, we used the AMSR-E ascending 36.5 GHz orbital swath data at the native footprint resolution of approximately 14 km × 8 km (Kawanishi et al., 2003).After the cessation of AMSR-E operations on 4 October 2011, its successor AMSR2 was launched on 18 May 2012 onboard the sunsynchronous JAXA GCOM-W1 satellite.AMSR2 is similar to AMSR-E in sensor configuration, including frequencies, incidence angle, and orbital equatorial crossing time.Major AMSR2 advancements over AMSR-E include an additional frequency at 7.3 GHz designed for mitigating radio frequency interference (RFI) and a larger (2.0 m diameter) main reflector for enhanced spatial resolution.The AMSR2 L1R (version 1.2) resampled ascending swath 36.5 GHz T b retrievals at approximately 12 km × 7 km resolution were used for this study.The uncalibrated AMSR2 T b retrievals were estimated to be positively biased against AMSR-E by ∼ 1.3 K (Du et al., 2014).However, the sensor inconsistency is expected to have minimal impacts on our algorithm, which relies on T b time series change signal detection rather than T b absolute accuracy. .Three sets of lakes selected for the LIP analysis in the Northern Hemisphere mid-and high latitudes (≥ 30 • N).The purple star symbols denote lakes used for evaluating LIP retrievals on a per pixel basis against GLRIPD ground-based observations, including Lake Superior in the USA and Canada and Lake Oulujarvi, Lake Haukivesi, and Lake Paijanne in Finland.The red star symbols denote 12 lakes (Great Bear Lake, Great Slave Lake, Smallwood Lake, Nettiling Lake, Dubawnt Lake, Amadjuak Lake, Wollaston Lake, Baker Lake, Kasba Lake, Lesser Slave Lake, and Peter Pond Lake in Canada and Red Lake in USA) used for the lake-wide comparisons between the LIP results from this study and other regional lake ice phenology records (IMS, CIS).The 71 lakes selected for assessing LIP trends over the 12-year satellite record are in bright blue.
The Level 1 GLWD (Lehner and Döll, 2004) comprises the 3067 largest lakes (area ≥ 50 km 2 ) and 654 largest reservoirs (storage capacity ≥ 0.5 km 3 ) worldwide and was used for identifying Northern Hemisphere water bodies and the 71 large lakes used for the LIP assessment (Fig. 1).We also used the MODIS 250 m land-water mask (MOD44W) for calculating the proportional water coverage of 5 km resolution pixels within lake areas identified by the GLWD (Carroll et al., 2009).

Datasets used for algorithm evaluation
Four lake ice phenology databases were used to evaluate the LIP retrievals, including (a) the GLRIPD (Benson and Magnuson, 2000), (b) the National Oceanic and Atmospheric Administration (NOAA) IMS 4 km daily snow and ice product (Helfrich et al., 2007; http://www.natice.noaa.gov/ims/),(c) the CIS lake-wide ice product (Howell et al., 2009), and (d) MODIS quick-look images for GBL downloaded from the Geographic Information Network of Alaska (http: //www.gina.alaska.edu).
The GLRIPD contains descriptive ice cover data for 865 lakes and rivers in the Northern Hemisphere (Benson and Magnuson, 2000).The GLRIPD includes groundbased (lakeshore) observations that were used for evaluating the corresponding LIP results for the targeted lakes.The GLRIPD records the first date when the water body was observed to be completely ice covered and the date when the last ice breakup was observed before the summer open water phase for each year of record (Benson and Magnuson, 2000).For evaluating LIP results representing 5 km lake dominant pixels (water coverage ≥ 90 %), the lake was assumed completely covered with ice for the period between the first date with complete ice cover and last ice breakup date as recorded in the GLRIPD, while lakes were classified as open water condition for other dates within each annual cycle.Only four lakes were selected for the LIP comparisons due to a predom-The Cryosphere, 11, 47-63, 2017 inance of ice observations from smaller lakes in the GLRIPD database.The temporal coverage of GLRIPD observations for the four lakes that overlaps with the AMSR-E/2 record extends from 2002 to 2007 for Lake Superior and 2003 to 2007 for Lake Oulujarvi, Lake Haukivesi, and Lake Paijanne.
The NOAA IMS daily snow and ice product provides snow and ice cover extent information derived from ground observations and an extensive variety of satellite observations, including AVHRR, GOES (Geostationary Operational Environmental Satellite), SSM/I, and AMSU (Advanced Microwave Sounding Unit) (Helfrich et al., 2007).The CIS lake ice product estimates lake ice cover fraction in tenths (0: open water; 10: complete ice cover) for nearly 140 lakes across Canada and the northern USA from visual interpretation of 1.1 km resolution NOAA AVHRR and 100 m resolution RADARSAT ScanSAR imagery (Howell et al., 2009), and MODIS (250-500 m) and Visible/Infrared Imaging Radiometer Suite (VIIRS) I-band (375 m) observations.The CIS product provides a single lake-wide value per lake on a weekly basis.Both the 4 km IMS grid products (year 2004-2015) and CIS data (year 2002-2015) were used for lakewide comparisons against the resulting LIP retrievals.
MODIS quick-look images in true-color composites (bands 1, 4, 3 in RGB) were selected for qualitative visual comparisons with the LIP results.The MODIS images were acquired over the breakup season (2012-2013) with clearsky conditions on 22, 27 June and 8 July 2013 and extensive cloud cover on 5 July 2013.The quick-look products were provided at 250 m resolution in Albers equal-area conic projection.
In addition, ERA-Interim (Dee et al., 2011) quarter-degree reanalysis surface air temperature (SAT) data was analyzed for evaluating LIP trends over the 71 Northern Hemisphere lakes selected (Fig. 1).ERA-Interim is a global atmospheric model data reanalysis produced by the European Centre for Medium-Range Weather Forecasts, and the data assimilation system used to produce ERA-Interim is based on a 2006 release of the IFS (Cy31r2) including a four-dimensional variational analysis (4D-Var) with a 12 h analysis window (Dee et al., 2011).Daily average SAT over the spring (MAM) and fall (SON) seasons of years 2002 to 2015 was extracted for the quarter-degree grids encompassing the lake centers.

Data processing
To derive the LIP estimates, AMSR-E/2 36.5 GHz orbital swath T b data were spatially resampled to a 5 km resolution polar EASE-Grid (version 2) format using an inverse distance squared weighting method (Ashcroft and Wentz, 2000;Brodzik et al., 2012Brodzik et al., , 2014)).It is worth noting that the T b spatial gridding is posted at 5 km resolution while the original 36.5 GHz AMSR-E/2 observations have coarser native sensor footprints (∼ 12 km for AMSR-E and 9 km for AMSR2).The finer grid spacing is intended to facilitate product comparisons and analyses with other alternative lake products de-rived at similar resolutions, including the NOAA IMS 4 km daily snow and ice product and a land surface fractional open water cover dataset derived from AMSR-E/2 at 5 km resolution (Du et al., 2016).Before carrying out the T b gridding process, an additional altitude correction was made to the AMSR-E data by considering the actual surface of the Earth instead of that of an ideal Earth ellipsoid.The same altitude correction was used for the AMSR2 L1R data (Maeda et al., 2016).According to Maeda et al. (2016), an altitude of 3000 m leads to about 4 km displacement of AMSR2 T b geolocation.Thus the altitude correction is a necessary prerequisite to ensure reliable analysis of AMSR-E/2 lake ice phenology retrievals at higher elevations, including the Qinghai-Tibetan Plateau.
The finer-resolution MOD44W static open water maps were aggregated to the same 5 km resolution polar EASEgrid 2.0 projection format as LIP and used with the GLWD to identify dominant lake pixels (water coverage ≥ 90 %) where the AMSR-E/2 lake ice detection was made.The 250 m resolution MODIS quick-look images were re-projected to the EASE-grid 2.0 projection for visual comparisons with the LIP results.

Algorithm theoretical basis
Accurate modeling of satellite observed microwave emissions from lakes is complex and requires good understanding of microwave scattering and emitting mechanisms from atmosphere and lake elements.Microwave emissions from a non-scattering atmosphere are governed by both air temperature and atmosphere optical thickness, which is approximately the sum of the optical thickness of oxygen, cloud liquid water, and atmospheric water vapor (Wang and Tedesco, 2007;Du et al., 2015).Microwave emissions from a lake with an upper layer that may consist of water, ice, and snow are determined by a number of factors; these factors include lake surface roughness, water dielectric properties mainly affected by water salinity and temperature, ice thickness and dielectric properties, and snow cover dielectric properties mainly controlled by snow density and wetness, snow particle size, and stratification of snow and ice layers (Du et al., 2010;Lemmetyinen et al., 2010Lemmetyinen et al., , 2011)).Despite the complexity of the lake emission problem, sharp changes in satellite microwave T b observations at multiple frequencies are evident during the transitions between lake freeze-up and breakup periods.For example, previous studies showed low T b measurements (< 150 K) from H-polarized 37 GHz SMMR data over low-emissivity open water regions of the Great Lakes and Gulf of Mexico, contrasting with much higher T b values (> 215 K) over western Lake Superior under frozen conditions due to the high emissivity of lake ice (Ferraro et al., 1986).Similarly, the H-polarized emissivity at 35 GHz and 50 • incidence angle is approximately 0.356 for a calm and unfrozen lake at 0-8 • C and is well below the emissivity (> 0.610) of different types of snow and ice (Mätzler, www.the-cryosphere.net/11/47/2017/The Cryosphere, 11, 47-63, 2017Cryosphere, 11, 47-63, 1994)).These studies suggest a very large Ka-band T b difference (> 60 K) between a lake at 0 • with no ice and 100 % ice coverage.The timing of ice formation and disappearance can therefore be determined by the large characteristic T b changes indicated from satellite passive microwave observations (Walker and Davey, 1993;Che et al., 2009;Kang et al., 2012).

Algorithm development
For identifying freeze-up and breakup events, a moving t test method (MTT) was introduced to detect abrupt temporal changes in the H-polarized 36.5 GHz T b observations from AMSR-E and AMSR2.Selection of 36.5 GHz T b observations from other available AMSR-E/2 frequencies represents a compromise between finer spatial resolutions gained from higher frequencies and less sensitivity to potential atmosphere contamination available from lower T b frequency observations.Moreover, H-polarization T b retrievals were used instead of V-polarization data due to their reported higher sensitivity to lake freeze-up/breakup signals (Kang et al., 2012).The detailed lake ice detection method used in this study is described below.

Step 1: detection of abrupt changing point
The MTT method was initially developed for detecting abrupt climate changes by examining whether the difference between the mean values of two subsamples is statistically significant (Jiang and You, 1996;Xiao and Li, 2007).As detailed in the literature (Xiao and Li, 2007), for a time series with n elements, a t test is made at each point x k for evaluating the difference of the two subsets x k1 and x k2 (n 1 ≤ k ≤ n − n 2 ; n 1 , n 2 are the subsample sizes) before and after x k .The t statistic is defined as where n 1 +n 2 −2 , x k2 and x k1 are the mean values, and s 2 k1 and s 2 k2 are the variances for the two subsets, respectively.Given a significance level α, x k is determined as an abrupt changing point if |t| ≥ t α .In this study, we define α as 0.005 and temporal subsample sizes as n 1 = n 2 = 20 days.The 20-day requirement is set for excluding potentially dynamic T b changes caused by short-term weather events such as storms.

Step 2: determining reference T b values for lake ice conditions
For a group of detected changing points sequenced from p to q, the mean T b values x p1 and x q2 as defined in Eq. ( 1) are representative of the satellite observations over the stable stages before and after the changing period, respectively.
For a lake experiencing a complete annual freeze-up/breakup cycle, at least two groups of seasonal changing points can be defined.Besides the sharp T b increases as lake water freezes, the melting of dry snow overlying lake ice to wet snow can induce further increases in the observed T b since microwave emissions from wet snow are close to that of a blackbody (Ulaby et al., 1986).Therefore, assuming x p1 is always smaller than x q2 , we define the lowest x p1 of all changing groups as the reference T b for lake water and x q2 from the same group as the reference T b for lake ice.Lake ice conditions for a given date i can thus be determined as where T b i is the T b for date i, T b threshold = (x p1 +x q2 )/2.0, and x p1 − x q2 is required to be larger than 30 K since liquidice phase changes of lake water can lead to large T b changes exceeding 60 K as introduced in Sect.2.2.

Step 3: deriving lake ice status
Based on Eq. ( 2), lake ice status is first derived for each point i in the T b time series where n 1 ≤ i ≤ n − n 2 using a temporally smoothed T b i defined as the mean T b within the range The use of a smoothed T b i minimizes the impact of high temporal frequency events in the time series while emphasizing lower frequency lake icecovered and ice-free signals.Thus for point j , whose temporally adjacent points have different lake ice status, the refined lake ice detection is carried out using Eq.(2) for each observed T b value within the range [j − n 1 /2, j + n 2 /2].For running the algorithm, missing daily T b retrievals were obtained through temporal linear interpolation of adjacent successful T b retrievals acquired from the same ascending orbits (Kim et al., 2012).However, only the lake ice detection results corresponding to the actual satellite observations were output for further analysis.The above lake ice detection process was carried out for each T b time series from AMSR-E and AMSR2 separately because of the 7-month gap (4 October 2011-18 May 2012) in the observation records between the two sensors.The above MTT algorithm was applied to all 5 km pixels with dominant (≥ 90 %) open water coverage within the Northern Hemisphere domain on a daily basis to generate the AMSR-E/2 LIP dataset describing lake ice conditions.The dominant (≥ 90 %) open water coverage criterion is set to include lake pixels while reducing potential contamination from adjacent land areas.

Evaluation of lake ice phenology retrievals
The resulting LIP retrievals were evaluated against other available lake ice databases, including the GLRIPD, IMS, and CIS products.The remotely sensed lake ice phenology variable definitions from this study are summarized in Ta- ble 1 relative to the other lake ice observational data records used for LIP validation on a per pixel basis and for entire lakes (Kang et al., 2012;Duguay et al., 2015a).
The LIP-derived annual CFO (complete freeze over) and WCI (water clear of ice) dates for the 12 selected lakes were also compared with alternative IMS and CIS ice products.Different from the dominant open water coverage (≥ 90 %) requirement set for generating the LIP database, only lake pixels with complete (100 %) open water coverage and outside of the 5 km land buffer zone were considered in the IMS and LIP comparisons; this same criterion was set for the lake-wide comparisons to minimize potential contamination from adjacent land areas since the native 36.5 GHz AMSR-E/2 footprint ranges from approximately 9 to 12 km.Considering possible retrieval uncertainties, the CFO-WCI dates derived from the LIP and IMS datasets were slightly adjusted from the definitions in Table 1 and were determined as the dates when most (99.5 % for this study) lake pixels were identified as ice/water.The CIS CFO dates were determined when the reported lake ice fraction was 9, followed by changes from 9 to 10, and WCI dates were derived when the reported lake ice fraction was 1 followed by changes from 1 to 0. The derived CIS CFO-WCI dates are comparable with corresponding LIP results, excluding waters adjacent to land such as part of the eastern arm of the GSL where a high concentration of islands exists and ice formation/melting timing was found to be different from other GSL areas (Howell et al., 2009;Kang et al., 2012).

Analysis of lake ice phenology changes
Based on the LIP database covering the AMSR-E (June 2002-October 2011) and AMSR2 (June 2012-December 2015) observation periods, we selected 71 lakes in the Northern Hemisphere from 250 of the world's largest lakes (including both natural and artificial lakes), as described in the GLWD, to analyze potential lake ice phenology trends, including WCI date, CFO date, and annual ICDe.In order to assess the pattern of recent Northern Hemisphere lake ice phenology changes, a temporal trend analysis was performed on the 12-year LIP record for each of the 71 lakes.The assumption of independent observations was first determined using a correlogram (Noguchi et al., 2011).For ice phenology time series without significant autocorrelation detected, the magnitude and significance of temporal trends were tested using the non-parametric Mann-Kendall and Sen's methods (Sen, 1968;Duguay et al., 2006).Alternatively, for a time series with persistent serial correlation, additional prewhitening approaches were applied (Zhang et al., 2000).For evaluating LIP-derived lake phenology, a similar temporal trend analysis was also carried out on the ERA-Interim daily average SAT over the lakes for the spring and fall seasons from 2002 to 2015.

LIP comparisons with GLRIPD lake observations
The lake ice status derived from the LIP and GLRIPD records are plotted for the selected large lake validation sites (Fig. 2), including Lake Superior (a), Lake Oulujarvi (b), Lake Haukivesi (c), and Lake Paijanne (d), along with the daily ascending AMSR-E/2 T b retrievals.The LIP results show generally strong agreement with the GLRIPD site observations of lake ice conditions for the four lakes examined, with overall retrieval accuracy of 95.4 %.The lake ice/water retrieval error at the beginning of the record for Lake Haukivesi (Fig. 2c) may be caused by partial melting www.the-cryosphere.net/11/47/2017/The Cryosphere, 11, 47-63, 2017 Figure 2. Comparison of lake ice status for Lake Superior, USA and Canada (a), Lake Oulujarvi, Finland (b), Lake Haukivesi, Finland (c), and Lake Paijanne, Finland (d), derived from the Global Lake and River Ice Phenology Database (GLRIPD) (blue dots) and AMSR Lake Ice Phenology retrievals (LIP) (red dots).The AMSR-E/2 36.5 GHz H-polarized daily T b retrievals used in the LIP algorithm are also plotted for reference (black line).The blue/red dots represent GLRIPD/LIP-derived ice conditions as indicated by their y axis positions for the dates described by their x axis coordinates. of lake ice in January and February 2003 that resulted in low AMSR-E/2 T b observations.While the AMSR-E/2 T b observations show dynamic daily fluctuations due to changing water and atmosphere properties (Sect.2.2), lake freeze-up and breakup events constitute the dominant factors affecting seasonal T b changes.The effects of higher temporal frequency T b variations are minimized in the LIP algorithm by the predefined 20-day subsample sizes (Sect.2.3), which represent a compromise between the algorithm's capability in capturing shorter-term lake ice formation or melt events, and potentially degraded lake ice/water seasonal retrieval accuracy.

LIP comparisons against MODIS imagery, IMS and CIS products
An example visual comparison between the LIP results and MODIS quick-look imagery (Fig. 3) shows the 2013 spring ice breakup process over the GBL.In this example, both datasets show a general onset of lake ice breakup on 22 June (Fig. 3a), similar spatial ice distribution patterns on 27 June (Fig. 3b), and ice-free conditions on 8 July (Fig. 3d).Despite extensive cloud presence in the MODIS image for 5 July, both MODIS and LIP show remaining ice cover on the western edge of GBL (Fig. 3c).A few remaining pixels along the GBL coast line were identified as ice covered in the LIP results for 8 July (Fig. 3d); the apparent LIP retrieval error is attributed to land contamination, while the affected pixels are within the 5 km land buffer zone (Sect.2.4) and excluded in the final LIP product.The LIP products were compared against similar lake ice phenology metrics from the IMS and CIS datasets for the 12 North American study lakes.For GBL and GSL, the LIP products agreed well with the CIS records; temporal corre-lations (R value) of 0.90 and 0.85 were observed for the WCI dates for GBL (Fig. 4a) and GSL (Fig. 4b), respectively, while correlations of 0.90 and 0.72 were determined for GBL (Fig. 4c) and GSL (Fig. 4d) CFO dates.The LIP results were also strongly correlated with the IMS record on derived WCI dates for both lakes (R = 0.94 for GBL and R = 0.91 for GSL) (Fig. 4a, b).However, lower correspondence was found between LIP and IMS CFO dates, with respective correlations of 0.54 and 0.63 for GBL and GSL (Fig. 4c, d).These results indicate that the LIP-derived lake ice phenology variables are generally consistent with the IMS and CIS records for the two lakes examined, with generally higher (lower) correspondence for WCI (CFO) dates.For GBL, the LIP estimated CFO dates occur similar to the CIS records (0-day difference) and later than the IMS records by about 3 days; the LIP WCI dates occur earlier than the CIS and IMS records by about 2 and 3 days, respectively.For GSL, the LIP record also shows earlier (later) CFO than the CIS (IMS) records by about 2 (3) days and earlier WCI dates than both CIS and IMS by about 6 and 5 days, respectively.The intercomparisons between CIS and IMS show average 5-and 0-day differences in respective CFO and WCI dates for their overlapping period from 2004 to 2015 for GBL; corresponding differences for GSL are 5 and 1 day, respectively.The differences between the LIP and CIS/IMS metrics are of similar magnitude as the differences between the CIS and IMS metrics.
The LIP comparison results for all 12 study lakes are summarized in Table 2. Similar to the comparisons for GBL and GSL, the LIP results are strongly correlated with both CIS and IMS records for WCI dates, with respective average temporal correlations (R value) of 0.86 and 0.84.For CFO www.the-cryosphere.net/11/47/2017/The Cryosphere, 11, 47-63, 2017  dates, the average correlation between LIP and CIS results are also strong (R = 0.69) while only moderate correlation (R = 0.41) was found between the LIP and IMS results.The LIP estimated CFO dates tend to occur earlier than the CIS record by about 6 days and later than the IMS record by about 1 day.The LIP record also shows earlier WCI dates than the CIS and IMS records by about 7 and 2 days, respectively.

Analysis of LIP lake ice phenology changes
The magnitude and direction of LIP trends were calculated for the 71 Northern Hemisphere study lakes for the 12-year AMSR-E/2 record.Among all 71 lakes, 43 (60.6 %) show declining trends in ICDe, indicating an increasingly shorter ice cover season, while the other lakes show either increasing or minimal change in annual ice cover (Fig. 5).However, no observed ice trends are statistically significant (p ≥ 0.05).The lack of significant trends is attributed to large yearly variability (±13.3 day) in average ICDe and a relatively short (12-year) LIP observation record.The changing trends also demonstrate a latitudinal pattern, as 81.0% of the lakes (17 out of 21) at higher latitudes (> 60 • N) show declining ICDe trends while only 45.0 % (9 out of 20) of lower-latitude (< 50 • N) lakes show a similar trend.The observed changes in ICDe are the net result of changes in fall CFO and spring WCI dates.A tendency toward earlier WCI dates was found for 40 lakes, including 5 lakes (Lake Vygozero, Lake Barun-Torey, Lake Segozero, Novosibirsk Reservoir in Russia and Lake Teshekpuk in USA) with significant LIP trends.However, no lakes showed significant trends toward later spring breakup.Similar to the ICDe analysis, most high-latitude lakes (81.0 % of lakes above 60 • N) show earlier spring thaw trends, while only 45.0 % of lower-latitude (< 50 • N) lakes show similar trends.A tendency toward delayed CFO was found for 35 of the 71 lakes (49.3 %) examined, but no trends are statistically significant (p ≥ 0.05).Lake Bosten in China was the only lake with a significant trend toward earlier freeze-up.There was no clear relationship between changes in lake CFO dates with latitude.
Similar analysis of quarter-degree ERA-Interim SAT over the study lakes indicates a much stronger warming trend in spring (0.073 • C year −1 ) than fall (0.023 • C year −1 ).Moreover, similar to the latitudinal pattern shown in the LIP-based analysis, the SAT increase in the spring is positively correlated with latitude (R = 0.33; p = 0.005) indicating greater warming during the study period at higher latitudes, while no SAT correlation with latitude is found for the fall (R ∼ 0.0).

Discussion
We developed a new satellite approach for regionally consistent classification of ice phenology for large lakes in the Northern Hemisphere from the AMSR-E/2 sensor record.We used similar 36.5 GHz H-polarization daily brightness temperature retrievals from AMSR-E and AMSR2 sensor records with 5 km posted spatial resolution.The resulting LIP record documents the timing and duration of seasonal lake phenology events of 5 km lake pixels (water coverage ≥ 90 %) over the 12-year AMSR-E/2 record.The LIP results showed strong agreement with GLRIPD site observations from four lakes, with agreement ranging from 92.4 to 98.7 %.Differences between the LIP and GLRIPD results can be attributed to several factors.First, each database has a different definition of lake ice conditions; lake ice coverage determined by satellite microwave sensors is dependent on ice thickness, which may vary from the ice detection approach used by on-site observers or observed from optical sensors.According to the literature (Hall et al., 1981), lake ice thickness and T b are linearly related for multiple frequencies (from 5 to 37 GHz).The reported maximum microwave penetration depths of fresh lake ice at 37 GHz range from 0.70 to 1.4 m, depending on ice temperatures (Chang et al., 1997;Surdyk, 2002;Kang et al., 2014).This implies that the formation of thin ice, resulting in relatively small T b increases, may not be detectable using the defined T b thresholds in the LIP algorithm.Differences between the LIP and GLRIPD results may also reflect spatial inconsistencies in lake observation area between the ground-based lakeshore observations and the coarser satellite footprint.Thus, the lake area observed on-site may not completely overlap with the AMSR-E/2 lake pixel used for the LIP classification.The GLRIPD also does not provide explicit descriptions of lake ice status for the period between the first date when the water body was completely ice covered and the date when the last ice breakup occurred; thus, short-term events such as temporary ice melting or formation may not be recorded in the GLRIPD.For example, though identified as ice covered in the GLRIPD, Lake Oulujarvi was more likely to have thawed on 9 January 2007 since the low T b (178.7 K) observation is more characteristic of open water emissions (Fig. 2b).The satellite T b observations at 36.5 GHz are also affected by other factors than surface freeze-up/breakup transitions, including changes in atmosphere water vapor and cloud liquid water (e.g., Sect.2.2).For example, the LIP detected iceon conditions for Lake Haukivesi Finland in mid-summer (30 July 2004) (Fig. 2c), which is likely incorrect and may be due to increased atmosphere water vapor concentrations under warm summer conditions, resulting in a large T b increase similar to a seasonal freeze-up event.
In the lake-based comparisons for the 12 lakes examined, including GSL and GBL, the LIP results show strong correspondence with the CIS product for both CFO and WCI dates and similar high correlations with the IMS results for WCI dates; however, the LIP WCI (CFO) dates differ by approximately 7 (6) days from the CIS and 2 (1) days from the IMS.These differences are attributed to the different sensor spatial/temporal resolutions and retrieval methods associated with the different products.As described in Sect.2.1.3,the CIS product is derived for individual lakes from visual interpretation of imagery from optical and SAR sensors and has a ±1-week accuracy due to the weekly product reporting.Both CIS and IMS products rely partially on observations from optical sensors such as AVHRR and their accuracy is influenced by adverse weather conditions, including the presence of cloud cover.IMS-derived lake ice products have been widely used in monitoring global climate change (Duguay et al., 2013(Duguay et al., , 2014(Duguay et al., , 2015b)); however, the IMS detected freeze onset was found to be too early for some lakes in northern Québec, presumably due to misclassification by inclusion of coarse-resolution satellite passive microwave observations during periods of prolonged cloud cover (Brown and Duguay, 2012;Brown et al., 2014); this may cause the low correlations between LIP and IMS CFO dates, as well as a delayed LIP CFO bias relative to IMS.
In addition, the relatively coarse spatial resolution of AMSR-E/2 observations limits capabilities for resolving lake ice conditions of finer-scale water bodies.Spaceborne and airborne optical-TIR and radar sensors are capable of improved delineation of smaller lakes and rivers (Chu et al., 2016), but at the expense of degraded temporal fidelity for regional and global applications.
As a proxy indicator of climate variability and change (Duguay et al., 2006), lake ice phenology variables and their changing trends are important for monitoring and understanding climate change and its feedbacks.As described above, only 5 of the 71 lakes examined showed statistically significant trends towards earlier WCI dates while no lake showed a significant later CFO trend.Earlier ice breakup events are signs of warmer spring conditions, which promote melting and breakup of lake ice and a lower surface albedo that absorbs more incoming solar radiation and further intensifies the rate of ice melt (Mishra et al., 2011).These results are also consistent with previous studies over Canada that found a general trend toward earlier springs and WCI dates particularly over western Canada but little change in isotherm and CFO dates in fall (Duguay et al., 2006).Our results also indicate that lakes at higher latitudes are more likely to experience trends toward earlier spring ice breakup and shorter ICDe, which is consistent with enhanced warming trends at higher latitudes (Solomon et al., 2007;Deutsch et al., 2008).The above ice phenology trends coincide with regional SAT trends from ERA-Interim that show an average spring warming rate that is more than triple that of fall, as well as stronger warming trends for higher latitudes.
Though the LIP lake ice phenology trends are generally consistent with regional climate warming (Magnuson et al., 2000;Solomon et al., 2007), further analysis based on a longer period of record is needed for distinguishing longterm climate trends from large interannual variability and periodic climate cycles, including the North Atlantic Oscillation (NAO), El-Niño Southern Oscillation (ENSO), and Pacific Decadal Oscillation (PDO) (Mishra et al., 2011).

Conclusions
Lake ice phenology is strongly influenced by variations in air temperature, while consistent long-term records of lake ice changes provide a sensitive climate change indicator (Magnuson et al., 2000;Weyhenmeyer et al., 2011).Continuous and accurate monitoring of lake ice dynamics is greatly needed for studies of global change and for monitoring lake ice impacts on ecosystems and infrastructure, especially for high-latitude and high-altitude regions.In this study, we developed a new automated algorithm for consistent daily retrieval of lake ice conditions over the Northern Hemisphere using similar 36.5 GHz H-polarized T b observations from AMSR-E and AMSR2 sensor records.The resulting 5 km resolution lake ice phenology record allows for daily monitoring of lake ice conditions without being significantly degraded by variations in solar illumination or cloud and atmosphere contamination effects.In particular, the LIP record distinguished 71 large lakes that satisfied 20-day minimum ICDe and open water season algorithm thresholds; these lakes represent approximately 23 % of the total surface area of large lakes (area ≥ 50 km 2 ) within the Northern Hemisphere domain.Smaller water bodies were excluded from the lake-wide analysis if the lakes had no pixels with complete (100 %) open water coverage outside of a 5 km land buffer zone.The relatively coarse spatial resolution of AMSR-E/2 observations limits capabilities for resolving finer-scale water bodies, while the conservative lake selection criterion minimizes potential land contamination effects.The LIPderived lake ice conditions were found to be largely consistent with GLRIPD ground-based observations, with an average agreement of 95.4 % for the four lakes examined.The LIP record also showed favorable correspondence with other lake CFO and WCI assessments defined from the CIS and IMS products for 12 large study lakes.The LIP, CIS, and IMS differences were attributed to the different data sources and methods used to construct the different products, including differences in spatial and temporal resolutions of observations, and distinct nature of optical and microwave remote sensing.Though the design of the LIP algorithm, including the MTT method, helps to identify lake breakup/freeze-up events, while minimizing other T b disturbances from shortterm weather events, atmosphere effects can still lead to retrieval errors, especially from persistent high atmosphere water vapor concentrations over high-latitude lakes in the summer.Based on the LIP record from 2002 to 2015, significant earlier melting of lake ice cover was detected for 5 of the 71 lakes examined in the Northern Hemisphere, while lakes at higher latitudes show a more evident warming trend toward earlier ice breakup and shorter ICDe than those at lower latitudes.As the operations from AMSR2 and similar sensors continue, the MTT algorithm will allow for automated retrieval and consistent monitoring of ice conditions for large Northern Hemisphere lakes into the future.

Data availability
The AMSR-E/2 derived Northern Hemisphere LIP record described in this study is publicly available through the following link http://files.ntsg.umt.edu/data/AMSRE2_LAKE_ICE_PHEN (Du et al., 2017).

Figure 1
Figure1.Three sets of lakes selected for the LIP analysis in the Northern Hemisphere mid-and high latitudes (≥ 30 • N).The purple star symbols denote lakes used for evaluating LIP retrievals on a per pixel basis against GLRIPD ground-based observations, including Lake Superior in the USA and Canada and Lake Oulujarvi, Lake Haukivesi, and Lake Paijanne in Finland.The red star symbols denote 12 lakes (Great Bear Lake, Great Slave Lake, Smallwood Lake, Nettiling Lake, Dubawnt Lake, Amadjuak Lake, Wollaston Lake, Baker Lake, Kasba Lake, Lesser Slave Lake, and Peter Pond Lake in Canada and Red Lake in USA) used for the lake-wide comparisons between the LIP results from this study and other regional lake ice phenology records (IMS, CIS).The 71 lakes selected for assessing LIP trends over the 12-year satellite record are in bright blue.

Figure 3 .
Figure 3. Comparisons between MODIS quick-look images (left column) and AMSR-E/2 LIP results (right column) for Great Bear Lake (GBL) on 22 June (a), 27 June (b), 5 July (c), and 8 July 2013 (d).The images are in the EASE-GRID version 2 polar projection format, consistent with the underlying AMSR-E/2 gridded T b dataset used for the LIP classification.

Figure 4 .
Figure 4. Comparisons of water clear of ice (WCI) dates for (a) Great Bear Lake and (b) Great Slave Lake and complete freeze over (CFO) dates for (c) Great Bear Lake and (d) Great Slave Lake, derived from the AMSR-E/2 Lake Ice Phenology (LIP) dataset developed in this study.The LIP results are compared against similar metrics derived from the NOAA/IMS (IMS) and Canadian Ice Service (CIS).Missing LIP data from 2011 to 2012 denote the period between the end of AMSR-E operations and the start of the AMSR2 record.

Figure 5 .
Figure 5. Changing trends of (a) ice cover duration (ICDe), (b) water clear of ice (WCI) dates, and (c) complete freeze over (CFO) dates of 71 lakes for the period 2002-2015.Lake changing trends are shown by bar symbols whose heights are proportional to the trend magnitudes; the significant trend lakes are marked by yellow stars, while purple triangles denote lakes where no trend was detected (rate of change is ∼ 0.0 day year −1 ).

Table 1 .
Definitions of remotely sensed lake ice phenology variables from this study in relation to the other lake ice observational datasets used for the LIP validation assessment on a per pixel basis and for entire lakes.

Table 2 .
Summary of the comparison results for water clear of ice (WCI) and complete freeze over (CFO) dates derived from AMSR Lake Ice Phenology (LIP), Canadian Ice Service (CIS) datasets for the period 2002-2015, and the NOAA/IMS (IMS) dataset for the period 2004-2015.R LIP,CIS/IMS denotes the correlation coefficient between the LIP and CIS/IMS datasets; D LIP,CIS/IMS is the average difference (unit: day) in WCI or CFO dates calculated by LIP minus CIS/IMS.LIP,CIS R LIP,IMS D LIP,CIS D LIP,IMS R LIP,CIS R LIP,IMS D LIP,CIS D LIP,IMS