Remote Sensing of Environmental Changes in Cold Regions: Methods, Achievements and Challenges

: Cold regions, including high-latitude and high-altitude landscapes, are experiencing profound environmental changes driven by global warming. With the advance of earth observation technology, remote sensing has become increasingly important for detecting, monitoring, and understanding environmental changes over vast and remote regions. This paper provides an overview of recent achievements, challenges, and opportunities for land remote sensing of cold regions by (a) summarizing the physical principles and methods in remote sensing of selected key variables related to ice, snow, permafrost, water bodies, and vegetation; (b) highlighting recent environmental nonstationarity occurring in the Arctic, Tibetan Plateau, and Antarctica as detected from satellite observations; (c) discussing the limits of available remote sensing data and approaches for regional monitoring; and (d) exploring new opportunities from next-generation satellite missions and emerging methods for accurate, timely, and multi-scale mapping of cold regions.


Introduction
Cold regions, including high-latitude and high-altitude landscapes, are experiencing climate warming with amplification at roughly twice the global rate for the Arctic region (>60 • N) [1,2] and Tibetan Plateau (TP) [3,4].Cold regions are typically characterized by the presence of permafrost, extensive snow and ice cover, and rich reserves of stored soil organic carbon in the northern regions and TP [5,6].Ecosystems within these regions are highly vulnerable to changes resulting from the rapid destabilization and melting of ice above the 0 • C isotherm [7,8], lengthening the annual non-frozen season [9,10], and thawing of carbon-rich permafrost soils [11].
Recent region-wide warming trends have altered vegetation, and interactions and feedbacks between the water, energy, and carbon cycles [12][13][14]; these changes have also resulted in a myriad of impacts to landscape function and ecosystem services [15].Warmer temperatures have reduced the duration of seasonal snow and ice cover over land, ocean, and inland water bodies [16][17][18][19].Longer snow and ice-free seasons have led to lower surface albedos and greater net energy loading, reinforcing regional warming trends [20][21][22].The alteration of seasonal snow and ice cover has also altered wildlife habitats and human mobility, including degrading the stability of snow and ice cover for winter travel [23].
Permafrost soils are estimated to store up to 1,600 billion tonnes of soil carbon, representing roughly twice the amount of carbon stored in the atmosphere [11].Warming soils have promoted permafrost degradation and active layer deepening, enhancing the mobilization and potential transfer of soil carbon to the atmosphere [24,25].Ground surface deformation from degrading permafrost has also increased the risk of damage to human infrastructure, including roads, pipelines, and buildings [26].Changes in permafrost properties greatly impact the surface water budget because the soil ice layers form a relatively impermeable barrier to soil drainage [27].Surface subsidence into the water table driven by the thawing of ice rich soils has increased surface water inundation and lake expansion in continuous permafrost areas (where more than 80% of the ground is underlain by permafrost).In contrast, extensive draining of lakes and wetlands has been observed in more degraded permafrost areas [28][29][30].
Satellite observations have indicated vegetation greening over northern latitude tundra, attributed to enhanced vegetation growth from a longer frost-free season, contrasting with vegetation browning in boreal forest and some tundra regions, that may result from greater drought stress due to warmer temperatures and a longer frost-free season [31].Boreal forests have also been affected by an increase in the frequency and severity of wildfires exacerbated by warmer and drier conditions, and insect-related disturbances [32,33].The net effect of these changes is a complex snow/ice, vegetation, soil, and wetland mosaic where the terrestrial water, carbon, and energy cycles are strongly coupled and interactive with the climate.
Remote sensing provides an unprecedented approach for characterizing the timing, magnitude, and patterns of environmental changes.This is especially advantageous for geographically remote cold regions, where site observations are often spatially sparse and temporally limited.The multi-scale nature of remote sensing also provides insight into the often emergent spatial and temporal patterns and properties of ecosystems that may not be fully identified nor understood when approached from the perspective of a local region.Earth parameter data records derived from optical-infrared (optical-IR) and microwave satellite observations spanning multiple decades are particularly valuable in distinguishing large characteristic natural climate variability from more subtle environmental trends in cold regions [14,19,34,35] and for the detection of local to regional disturbances [36,37].Finer spatial-resolution optical-IR and active microwave sensors are essential for distinguishing the complex spatial heterogeneity in permafrost landscapes [38][39][40] and for near-real time applications such as monitoring river ice jams [41] and glacial lake outburst flooding [42].
This paper provides an overview of recent progress and prospects in remote sensing of cold regions by first reviewing general principles and methods in measuring a selection of key environmental variables, including glacier ice; snow; surface water bodies; permafrost and surface deformation; vegetation and terrestrial carbon process.We then summarize recent environmental changes documented by the remote sensing data record.Finally, we explore opportunities for leveraging available and future satellite missions, and integrating emerging remote sensing and big data techniques to establish a next-generation monitoring system for cold regions.

Principles and Methods
Satellite optical-IR and microwave sensors (Supplementary Table S1) have provided complementary observations of cold regions since the 1970s.In general, optical-IR sensors are well suited for mapping environmental variables over heterogeneous landscapes due to their relatively high-resolution (sub-meter to 1 km) imaging capability, though the signal-to-noise ratio of the observations may be degraded by cloud-atmosphere aerosol contamination, low solar elevation, and long periods of seasonal darkness at higher latitudes.Microwave remote sensing is less affected by atmospheric conditions and provides earth observations day-or-night under nearly all-weather conditions [43].The microwave penetration ability is generally superior to optical-IR wavelengths and depends on sensor frequency and landscape conditions.Lower-frequency (e.g.1-2GHz or L-band, 0.3-1GHz or P-band) observations provide better measurements of forest biomass and enhanced soil sensitivity under low to moderate vegetation and snow cover, while higher frequency (e.g., 12-18GHz or Ku-, 26.5-40 GHz or Ka-band) signals are more suitable to detect sparsely vegetated soil, snow properties (e.g., snow depth, surface roughness, stratification, and microstructures) [44,45] and ephemeral surface freeze-thaw (FT) conditions [38,[46][47][48].Among microwave sensors, satellite synthetic aperture radar (SAR) measures backscatter signals at relatively high spatial resolution (1-100s m), though the utility of operational SARs has been constrained by limited data access, incomplete global coverage and low temporal sampling.Alternatively, satellite microwave scatterometers and radiometers provide global coverage and frequent sampling (i.e., every 1-3 days) valuable for monitoring environmental dynamics over large regions, but at relatively coarse spatial resolution (~5-36 km).
A variety of sensor configurations and remote sensing techniques have been applied for monitoring cold regions, based on radiative transfer theory and the unique spectral signatures of various target variables.These approaches are summarized in the following subsections for selected variables where remote sensing has been used to document significant environmental changes attributed to global warming.

Glacier Mass and Movement
Glaciers are slow moving masses of ice formed over time by the accumulation and compaction of snow, holding 75% of Earth's freshwater [49] and 10% of the global land area, including most of Greenland and Antarctica [50].Glacier mass balance is highly sensitive to climate change and controls a glacier's long-term behavior and evolution.A glacier flows under its own weight due to the pull of gravity, and thus transports ice mass to lower altitudes.Remote sensing is the only practical approach for inferring glacier mass and movement over large regions.
Glacier mass is commonly estimated through independent or combined gravimetry and altimetry measurements [51].Satellite and aircraft-based gravimetry measurements have been widely used in glacier mass change assessments [52].Glacier mass can also be indirectly estimated through glacier area and thickness measurements.Glacier area change events, such as ice calving, can be precisely detected using satellite images acquired over different times [53].Glacier thickness change caused by ice melting or accumulation can be measured via geodetic approaches, including point measurements from altimetry and digital elevation models (DEMs) derived from photogrammetry or interferometric SAR (InSAR) techniques (Supplementary Table S2).Satellite laser altimeter measurements can achieve decimeter to centimeter accuracy levels; for example, the Geoscience Laser Altimeter System (GLAS) onboard the Ice, Cloud, and land Elevation Satellite (ICESat) provided decimeter-accuracy elevation data with a 70-m ground footprint over the global ice sheets [54].The ICESat-2 satellite was launched in late 2018 and has a significantly improved laser system providing observations with enhanced spatial resolution, temporal sampling, and measurement accuracy [55].Alternatively, satellite photogrammetry using optical-IR (e.g., ASTER, IKONOS) and SAR/InSAR (e.g., ERS-1/2 and ENVISAT) measurements have been successful in providing glacier raster DEMs [56,57].In particular, the SAR Interferometer Radar Altimeter (SIRAL) onboard Cryosat-2 is capable of measuring changes in the thickness of both sea ice and land ice under three different measurement modes (low resolution, SAR and InSAR) [58].
Glacier movement can be detected from repeat-pass satellite images using feature tracking and Differential InSAR (DInSAR) techniques.The feature tracking method identifies and matches ice surface features from satellite images and calculates the moving distance of the features over different acquisition times.The methods have been applied to both optical-IR and SAR image series, including Landsat Operational Land Imager (OLI) [59], Moderate Resolution Imaging Spectroradiometer (MODIS) [60], and ERS-1/2 SAR [61] sensors.DInSAR uses repeat-pass SAR imagery to calculate glacier motion velocity after removing the topography signals from the sensor interferograms using an external DEM or a combination of interferograms [62,63].The accuracy of the feature-tracking method can be within the sub-pixel level, while that of DInSAR is up to half of the radar wavelength.

Lake Ice Cover
Ice cover plays an important role in lake-atmosphere interactions at high latitudes.The presence (or absence) and extent/concentration of ice cover on large lakes has a significant impact on regional weather and climate (e.g., lake-effect snowfall, thermal moderation effect) [64][65][66][67][68]. Ice cover (extent) and ice thickness have recently been identified as Essential Climate Variables by the Global Climate Observing System (GCOS) of the World Meteorological Organization [69].Both ice extent, from which ice phenology (i.e., ice dates during freeze-up and break-up, and ice cover duration) can be determined, and ice thickness are sensitive indicators of climate change [70,71].Not identified by GCOS is the bedfast ice regime of shallow Arctic/sub-Arctic lakes (less than about 3-m).Such lakes are widespread across permafrost regions of Alaska, Northern Canada, and Siberia.Determining if and when entire lakes or lake sections become bedfast (i.e., ice cover is thick enough to reach lake bottoms) in winter has been shown to also be relevant for climate monitoring [72,73].Winters with a larger (smaller) fraction of bedfast ice are generally indicative of colder (warmer) air temperature and/or lower (higher) on-ice snow depth conditions which can lead to thicker (thinner) ice.Considering the sparse distribution of weather stations in northern high latitudes, whose temperature measurements are not representative for large areas, satellite remote sensing provides an alternative to measure regional ice cover extent (phenology), ice thickness, and bedfast ice as summarized below.
Ice cover extent and phenology-Satellite remote sensing has assumed a greater role in lake ice observations in recent years due to the dramatic reduction in ground-based observational recordings and the availability of increasingly longer satellite time series, particularly from the 2000s onward [74].Ice cover extent products are either generated manually, largely from visual interpretation of multi-source/frequency satellite imagery such as the National Oceanic and Atmospheric Administration (NOAA) National Ice Center Interactive Multisensor Snow and Ice Mapping System (IMS).The IMS products are produced manually through assimilation of various sources of data, including polar-orbit and geostationary satellite imagery and in situ data.In some cases, automated algorithms are applied to these data to facilitate analysis.The IMS products are available at various resolutions (1 km, 4 km, and 24 km) [75].Ice phenology dates (freeze-up/ice-on and break-up/ice-off dates) and ice cover duration can also be derived from the IMS products [74,76].MODIS 500-m snow (MOD10A1/MYD10A1) and 250-m surface reflectance (MOD09GQ) products have been used in a few recent studies, alone or in combination with each other and 1-km MODIS (MOD11A1/MYD11A1) Land Surface Temperature (LST), to derive ice dates (start and end of break-up and freeze-up dates) and their associated trends (2001-2017 or shorter) [77][78][79][80][81][82].Approaches that use top-of-atmosphere or surface reflectance (e.g., MODIS near-infrared and red bands) are based on threshold values where ice is determined to be present/absent above/below a certain value.High solar zenith angle, which is important during freeze-up for high-latitude lakes, and cloud cover are two factors that affect the quality of lake ice products.Hence, research has also focused on developing ice phenology retrieval algorithms from passive and active microwave observations.At passive microwave frequencies (e.g., 18-37 GHz) used for satellite remote sensing of lake ice cover, nadir emissivity from open water is low (ε = 0.443-0.504at 24 GHz) compared to that of ice (ε = 0.858-0.908at 24 GHz) [83].This makes the determination of the timing of ice formation and decay on lakes feasible from brightness temperature (Tb) measurements.The emissivity of ice, and therefore Tb, further increases during ice formation, as the influence of radiometrically cold water under the ice cover decreases with ice thickening [84].Kang et al. (2012) found AMSR-E (Advanced Microwave Scanning Radiometer for Earth Observing System) 18.7 GHz (H-pol) Tb data (interpolated onto a 10-km grid) to be the most suitable for estimating ice dates (freeze-onset, ice-on dates, melt-onset, and ice-off dates), as well as the duration of ice cover and ice-free seasons using a thresholding approach.Du et al. [17] further demonstrated that ice dates could be derived from AMSR-E and AMSR2 (Advanced Microwave Scanning Radiometer 2) 36.5 GHz (H-pol) Tb (re-gridded at 5-km) data using a moving t-test algorithm.Derived ice dates compared favorably with those obtained from ground-based observations and other satellite products such as IMS.
Threshold-based and semi-automated (region-based segmentation followed by manual labelling of ice/open water) approaches have also been developed to generate ice cover extent and phenology products using SAR data.Wang et al. [85] evaluated the semi-automated segmentation algorithm "glocal" Iterative Region Growing with Semantics (IRGS) [86] for lake ice classification using dual polarized (HH and HV) RADARSAT-2 imagery acquired over Lake Erie.Their analysis showed that the algorithm could provide reliable discrimination between ice and open water with high overall classification accuracy (90.4%) when compared to Great Lakes image analysis charts from the Canadian Ice Service.Murfitt et al. [87] developed a threshold-based approach for estimating ice phenology events for mid-latitude lakes in Central Ontario by tracking the temporal evolution in backscatter from HH-polarization RADARSAT-2 imagery (2008-2017).The authors reported mean absolute errors of 2.5-10 days for freeze events and 1.5-7.1 days for water clear-of-ice when compared to MODIS imagery.The method was also successful in detecting multiple freeze (high backscatter) and melt (low backscatter) events throughout the ice season.By combining acquisitions from ENVISAT Advanced SAR (ASAR) wide swath and RADARSAT-2 ScanSAR data, Surdu et al. [88] showed the advantage of more frequent sampling (i.e., every 2-5 days over the 2005-2011 study period), but also the need for sensor incidence angle correction for more precise ice phenology detection from backscatter thresholds over Alaskan North Slope lakes.
Ice thickness-Field measurements of ice thickness are spatially and temporally sparse in cold regions.Recent investigations have developed approaches to estimate ice thickness from passive microwave, active microwave (altimetry and SAR) and thermal remote sensing data.Kang et al. [84] showed that the temporal evolution of Tb measurements from AMSR-E at 10.7 GHz and 18.7 GHz frequency (V polarization) during the ice growth season on Great Bear Lake (GBL) and Great Slave Lake (GSL), Canada, is strongly related to ice thickness.The authors proposed simple linear regression equations to estimate ice thickness for the lakes using 18.7 GHz V-pol data (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009), while the estimated ice thicknesses compared favorably with in situ measurements (Mean Bias Error, MBE, 6 cm; Root Mean Square Error, RMSE, 19 cm).Beckers et al. [89] explored waveforms from CryoSat-2 Ku-band radar altimetry to estimate ice thickness also on GBL and GSL.Their study obtained ice thickness estimates with RMSE of 33 cm when compared to in-situ measurements obtained at GSL. Murfitt et al. [90] evaluated RADARSAT-2 data for estimating lake ice thickness in Central Ontario, Canada.They reported RMSE values of 11.7 cm and attributed the uncertainty to unexplored questions about scattering mechanisms and the interaction of the radar signal with lake ice having complex structure within the ice layer and at the ice-water interface.In addition to the radar-based investigations, lake surface (ice/snow) temperature observations from MODIS have also been evaluated for estimating lake ice thickness [91].Using heat balance terms derived from the Canadian Lake Ice Model [92], the authors retrieved ice thicknesses up to 1.7 m from MODIS with RMSE of 17 cm and MBE of 7 cm when compared to field measurements acquired on GSL and Baker Lake, Canada.Work on the estimation of ice thickness from satellite remote sensing is still in its infancy.Biases in retrievals are relatively large from some spaceborne instruments, while the remote sensing time series are generally too short to analyze climate trends in satellite-derived ice thickness.
Bedfast ice-Radar remote sensing allows for distinguishing lakes with bedfast (grounded) ice due to the difference in backscatter intensities between floating ice (generally higher backscatter return) and bedfast ice (lower return) [93].Recent analyses of polarimetric SAR (X-, C-, and L-band) satellite and ground-based scatterometer (Ku-and X-band) measurements, supported by radiative transfer modeling experiments, have revealed that the high backscatter of floating ice on shallow Arctic lakes is from the ice-water interface (due to appreciable surface roughness or preferentially oriented ice facets), dominated by single-bounce scattering [94,95].Areas of bedfast and floating lake ice are monitored/mapped from SAR using image thresholding [73] or region-based segmentation approaches [96].Analyses of C-band SAR time series (ERS-1/2, RADARSAT-1/2, ENVISAT ASAR, and Sentinel-1) have been used to document trends and variability in bedfast ice across Alaska, over the last 20-25 years [73,88].Antonova et al. [97] have also shown the potential of a unique time series of three-year repeat-pass TerraSAR-X imagery with higher temporal (11 days) and spatial (10 m) resolutions than available in past studies for monitoring both bedfast ice and lake ice phenology in the Lena River Delta, Siberia.The authors also analyzed an 11-day sequential interferometric coherence time series from TerraSAR-X as a supplementary approach for bedfast ice monitoring.Coherence time series have been found to detect most areas of bedfast ice as well as spring snow/ice melt onset.

Remote Sensing of Snow
Snow and glaciers provide one-sixth of the world's population with fresh drinking water, and seasonal snow is the main fresh-water source at mid-latitudes [98].Snow is also a crucial factor controlling the seasonal radiation balance of the land surface, and a sensitive indicator of global climate change.Snow measurement is essential to snowmelt driven runoff predictions, water resources management, flood control, and climate change studies [99].Key snow properties derived from remote sensing include snow cover area or extent, structure (e.g., depth, density), and water equivalent.

Snow Cover Area
Snow cover area has been estimated using satellite optical sensors such as Landsat TM, Aqua/Terra MODIS, and NOAA AVHRR (advanced very-high-resolution radiometer).Snow cover can be identified under clear-sky conditions using the Normalized Difference Snow Index (NDSI), which exploits the contrasting reflectance of snow in the visible and short-wave infrared bands [100,101].Utilizing ancillary spectral indices such as Normalized Difference Vegetation Index (NDVI) and Normalized Difference Forest Snow Index (NDFSI) helps incorporate vegetation information in snow detection and improves performance in mapping forest snow cover [102,103].Radiation transfer models can also be used for improved snow mapping over forest areas [104].Recent snow mapping efforts have focused on generating long-term snow cover products using observations from multiple satellite sensors [105], high-spatial resolution snow mapping using Sentinel-2 and Landsat optical imagery [106,107] and machine-learning techniques [108,109].In addition, satellite sub-pixel snow cover is valuable for improved estimation of snowmelt runoff and understanding energy exchanges between the land surface and atmosphere [110].Two main approaches to derive sub-pixel snow fraction include empirical linear regression of NDSI [111] and spectral mixture analysis [112][113][114][115][116].
Cloud contamination can significantly limit the signal quality of snow property detections made by satellite optical-IR remote sensing.Daily composites of half-hourly to hourly observations from geostationary satellites enable optimized snow detection while suppressing cloud contamination effects.Automated snow mapping has been achieved for a variety of geostationary satellites including GOES (Geostationary Operational Environmental Satellite) over North America [117]; Meteosat Second Generation (MSG) satellites over Europe [118,119]; and Multifunctional Transport Satellites (MTSAT)-2, Himawari-8, and Feng Yun (FY)-2 over Asia [120,121].A linear interpolation method [117] has been used to derive fractional snow cover from FY-2 satellite observations [122].Cloud-free snow cover products have also been derived using combined observations from polar-orbiting and geostationary satellites [121], additional observations from satellite microwave sensors [123][124][125][126]; and through processing of satellite data using advanced spatiotemporal filters and interpolation techniques [127,128].Optical-IR sensors have limitations in distinguishing between dry and wet snow, which can be effectively addressed by incorporating SAR observations.

Snow Water Equivalent
Snow water equivalent (SWE) describes the amount of water contained in the snowpack when completely melted.To estimate SWE from satellite microwave observations, the scattering and emission contributions from intervening atmosphere, snow surface, snowpack, and underlying soil and snow-soil interactions need to be distinguished and accounted for, which can be done by exploiting the frequency-dependent sensitivity of microwaves to land surface components [129,130] (Figure 1).Snow water equivalent (SWE) describes the amount of water contained in the snowpack when completely melted.To estimate SWE from satellite microwave observations, the scattering and emission contributions from intervening atmosphere, snow surface, snowpack, and underlying soil and snow-soil interactions need to be distinguished and accounted for, which can be done by exploiting the frequency-dependent sensitivity of microwaves to land surface components [129,130] (Figure 1).Satellite passive microwave SWE estimation relies on multi-frequency observations from SMMR (Scanning Multichannel Microwave Radiometer), SSM/I (Special Sensor Microwave/Imager), SSMIS (Special Sensor Microwave Imager Sounder), AMSR-E/2, and FY-3 MWRI (Microwave Radiation Imager) [125] sensors.The associated SWE algorithms include (a) static [130] and dynamic semiempirical algorithms [131,132]; (b) iterative algorithms [133]; (c) physically based statistical algorithms [134]; (d) probabilistic approaches [135]; (e) machine learning methods [136,137]; and (f) data assimilation methods [138,139].
The satellite SAR SWE retrieval algorithms can be grouped into two categories: (a) physical inversion algorithms and (b) interferometry methods.By utilizing the frequency-dependent sensitivity to snow and underlying soil properties, combined multi-frequency SAR observations (e.g., X-and Ku-band) are capable of SWE retrievals as demonstrated by model simulations and field experiments [140][141][142][143].The uncertainties related to snow density, ice microstructure, snow layer stratification, vegetation, and terrain effects are the main issues affecting the performance of both passive and active microwave snow retrieval algorithms [144][145][146].Alternatively, SAR interferometry techniques show promise in overcoming many of the above difficulties by utilizing interferometry phase difference information for SWE estimation [147,148].

Landscape Freeze/Thaw States
The landscape FT status is closely linked to ecosystem carbon, water and energy exchanges, snow melt dynamics, and permafrost extent and stability [21,34,[149][150].Global FT observational data records spanning the modern satellite era have been used to document environmental trends from global warming, including earlier and longer non-frozen seasons as a driver of northern vegetation greening, increased trends in damaging frost events in early spring [9], degrading permafrost, and an earlier spring flood pulse across the pan-Arctic [151].By availing of the high sensitivity of microwaves to significant dielectric shifts between solid and liquid water phase transitions, the FT signals can be captured by classifying Tb or radar backscatter coefficient (σ 0 ) changes relative to frozen or non-frozen reference conditions [35].
The satellite SAR SWE retrieval algorithms can be grouped into two categories: (a) physical inversion algorithms and (b) interferometry methods.By utilizing the frequency-dependent sensitivity to snow and underlying soil properties, combined multi-frequency SAR observations (e.g., Xand Ku-band) are capable of SWE retrievals as demonstrated by model simulations and field experiments [140][141][142][143].The uncertainties related to snow density, ice microstructure, snow layer stratification, vegetation, and terrain effects are the main issues affecting the performance of both passive and active microwave snow retrieval algorithms [144][145][146].Alternatively, SAR interferometry techniques show promise in overcoming many of the above difficulties by utilizing interferometry phase difference information for SWE estimation [147,148].

Landscape Freeze/Thaw States
The landscape FT status is closely linked to ecosystem carbon, water and energy exchanges, snow melt dynamics, and permafrost extent and stability [21,34,149,150].Global FT observational data records spanning the modern satellite era have been used to document environmental trends from global warming, including earlier and longer non-frozen seasons as a driver of northern vegetation greening, increased trends in damaging frost events in early spring [9], degrading permafrost, and an earlier spring flood pulse across the pan-Arctic [151].By availing of the high sensitivity of microwaves to significant dielectric shifts between solid and liquid water phase transitions, the FT signals can be captured by classifying Tb or radar backscatter coefficient (σ 0 ) changes relative to frozen or non-frozen reference conditions [35].
FT algorithms applied to both active and passive microwave observations include threshold-based methods [35,38,152], change-detection approaches [48,153], and multi-channel data fusion or machinelearning methods [154][155][156][157][158]. The threshold-based methods determine FT conditions by comparing the satellite observations with reference Tb or σ 0 values representative of seasonal frozen and non-frozen conditions.Such approaches are robust and relatively easy to implement for operational FT retrievals.More sophisticated multi-channel combinations, decision tree, and probabilistic model methods can effectively distinguish FT conditions from precipitation events in sparsely vegetated and drier climate zones.As a recent development, data fusion, and machine learning methods show promise in providing potentially enhanced FT retrievals by exploiting massive archives of satellite observations and ancillary data [157][158][159].
Multiple global FT data records have been developed using observations from satellite microwave radiometers and scatterometers, including SMMR and SSM/I[S] [160], AMSR-E/2 [161], Aquarius [162], SMOS (Soil Moisture and Ocean Salinity) [46], SMAP [163,164], and ASCAT [165].Long-term (>39-year) global daily FT data records have been developed using similar overlapping 37 GHz Tb retrievals from SMMR and SSM/I[S] sensors with moderate (~25km) spatial resolution [160].Finer (~12km) resolution FT data have also been developed using calibrated 36.5 GHz orbital swath Tb records from the AMSR-E/2 sensors [161,166].The SMAP mission provides an operational FT data record derived from L-band (1.4 GHz) Tb retrievals with global coverage, 1-3-day temporal fidelity, and 9-km and 36-km resolution gridding [163, 164,167].Example maps of the observed non-frozen days in 2017 are shown (Figure 2) from three different FT data products and operational satellite sensors, including SMAP [167], SSM/I [160], and AMSR2 [161].The SSM/I and AMSR2 FT records are derived from vertically polarized Tb retrievals and a modified single channel algorithm.The SMAP FT record is derived using a dual algorithm approach including a normalized polarization ratio (NPR) of vertically and horizontally polarized Tb differences from NPR reference states, and a single channel algorithm where conditions are unfavorable for the NPR.All of these records show similar FT regional patterns, including generally fewer frost days at lower latitudes and altitudes, and in coastal areas relative to higher latitude, alpine, and inland areas.However, the SMAP and AMSR2 products show generally enhanced spatial delineation of FT patterns due to the respective finer 9-km and 6-km resolution gridding of these products, relative to the 25-km resolution SSM/I global grid product.The SMAP L-band products also have greater soil FT sensitivity than the K-band retrievals from AMSR2 or SSM/I [152].
In addition, the SMAP radar produced operational FT retrievals over northern (≥45 • N) land areas with ~3-km resolution and 1-3-day fidelity until the radar transmitter failed in July, 2015; however, these data are overlapping with SMAP radiometer based FT records which share the same satellite antenna receiver [152].FT retrievals have also been acquired from L-band Global Navigation Satellite System (GNSS) signals captured from the SMAP radar receiver, which has provided kilometer scale observations of FT seasonal transitions.
FT regional patterns, including generally fewer frost days at lower latitudes and altitudes, and in coastal areas relative to higher latitude, alpine, and inland areas.However, the SMAP and AMSR2 products show generally enhanced spatial delineation of FT patterns due to the respective finer 9-km and 6-km resolution gridding of these products, relative to the 25-km resolution SSM/I global grid product.The SMAP L-band products also have greater soil FT sensitivity than the K-band retrievals from AMSR2 or SSM/I [152].

Surface Deformation
In cold regions, climate change could significantly affect surface morphology (e.g., permafrost melting due to global warming), which may pose a threat to 70% of current infrastructure in Arctic permafrost regions by 2050 (Figure 3) [26].Remote sensing technologies could offer a useful platform to better understand these geomorphic changes, especially those that guarantee high-resolution topography analysis [168].
Remote Sens. 2019, 10, x FOR PEER REVIEW 9 of 36 In addition, the SMAP radar produced operational FT retrievals over northern (≥45˚N) land areas with ~3-km resolution and 1-3-day fidelity until the radar transmitter failed in July, 2015; however, these data are overlapping with SMAP radiometer based FT records which share the same satellite antenna receiver [152].FT retrievals have also been acquired from L-band Global Navigation Satellite System (GNSS) signals captured from the SMAP radar receiver, which has provided kilometer scale observations of FT seasonal transitions.

Surface Deformation
In cold regions, climate change could significantly affect surface morphology (e.g., permafrost melting due to global warming), which may pose a threat to 70% of current infrastructure in Arctic permafrost regions by 2050 (Figure 3) [26].Remote sensing technologies could offer a useful platform to better understand these geomorphic changes, especially those that guarantee high-resolution topography analysis [168].Liu et al. [169] used spaceborne InSAR data to map surface subsidence trends at a thermokarst landform located near Deadhorse on the North Slope of Alaska.The motivation of this study was the fact that the intrinsic dynamic thermokarst process of surface subsidence remains a challenge to quantify and is seldom examined using remote sensing methods.Subsequent InSAR analysis using Phased Array type L-band SAR (PALSAR) images revealed localized thermokarst subsidence of a few cm yr −1 between 2006 and 2010.Luo et al. [170] used terrestrial laser scanning for quantifying surface deformation pertaining to underlying hydrological-thermal processes affecting active layer FT conditions in the Tibetan Autonomous Region (China).The Terrestrial Laser Scanner and six Trimble 5700 GNSS systems were deployed in the region between May 2014 and October 2015, where the site was monitored four times.The results indicated that as air temperature and precipitation increase with climate warming the active layer will become more unstable, exacerbating slope instability as phase changes (thawing and freezing) occur.Jorgenson and Grosse [171] summarized recent developments (2010-2015) in remote sensing applications to detect and monitor landscape changes in permafrost regions, analyzing surface temperatures, snow cover, topography, surface water, vegetation cover, and disturbances from fire and human activities.According to this review, repeated light detection and ranging (LIDAR), InSAR, and airborne geophysics will be key tools for monitoring permafrost changes (topographic and subsurface) in Arctic and boreal regions.Arenson et al. [172] stressed the fact that in situ monitoring of periglacial dynamics is essential for the study of periglacial morphology and the design of mitigation and adaptation measures for infrastructure in permafrost zones.The application of structure-from-motion photogrammetric techniques is relatively low-cost and easy to use (e.g., see [173] for details), and provide capabilities for multi-temporal surveys of surface deformation.Meng et al. [174] used X-band SAR Interferometry for the detection of surface deformation in the Sichuan-Tibet Grid Connection Project Area (China).In this area, landslides, and debris flows triggered by climate change are becoming a major threat.Surface deformation time series observations were obtained through sequential TerraSAR X-band images.The analysis suggested that the deformation rate tends to increase dramatically during the late spring and late autumn, but with reduced deformation during colder winter months.Stettner et al. [175] discussed the capability of high spatiotemporal resolution X-band microwave satellite data (obtained from the TerraSAR-X satellite) to quantify cliff-top erosion of an ice-rich permafrost riverbank in the central Lena Delta.The results indicated continuous erosion from June to September in 2014 and 2015 with no significant seasonality across the thawing season.The authors identified X-band backscatter time series as a useful tool complementing optical remote sensing and in situ monitoring for rapid analysis of tundra permafrost erosion along riverbanks and coastal areas.Chen et al. [176] used Persistent Scatterer Interferometry (PSI) to map and quantify permafrost thaw subsidence in the Qinghai-Tibet Plateau.According to the authors, the PSI approach is less affected by temporal or geometric decorrelation, while their results indicated that permafrost areas near gullies are more vulnerable to gradual thawing and degradation.

Remote Sensing of Water Bodies
Surface water (SW) over inland areas is a key component of the water and energy cycles, covering about three percent (4.46 million km 2 ) of Earth's land [177].The dynamics of SW in cold land regions are closely linked to terrestrial water storage changes [178], flood and drought events [179], seasonal thawing and the spring flood pulse [180], microtopography, underlying geology and permafrost conditions [181].SW changes are also occurring in Arctic-boreal wetlands, lakes, rivers, and streams as permafrost degrades with regional climate warming [29,180,182]; surface subsidence during the initial stages of permafrost degradation leads to increased inundation, while later stages of permafrost thaw lead to surface drying and reduced wetland extent as drainage pathways increase [27].The emerging glacier and thermokarst lakes formed as ice melts have a strong climate feedback and may increase regional hazards from outburst flooding [183,184].
Clear and calm water appears dark and is readily distinguished from surrounding land features using optical-IR, microwave radiometer, and mono-static radar sensors.Optical-IR satellite sensors such as PlanetScope multispectral cameras, MODIS, AVHRR, and Landsat provide potential daily to 16-day repeat global observations of SW cover at moderate to high-resolution (3-1000 m), while screening and temporal compositing of the data to reduce the influence of cloud and atmosphere contamination may degrade temporal fidelity [177,181,185].
Potential drawbacks from optical-IR sensors can be partially overcome by active [186] and passive [178,[187][188][189] microwave remote sensing.Daily satellite passive microwave observations were used to monitor spatial variability and multi-year trends in surface inundation in permafrost affected regions [29,188].Lower frequency (e.g., L-band) microwave retrievals have shown greater sensitivity and detection of surface water even under low to moderate vegetation cover [190].The synergistic use of optical-IR and microwave observations, and ancillary hydrologic information has shown promise in producing optimum SW mapping results in terms of accuracy, and spatial and temporal resolution [190][191][192][193][194][195].

Remote Sensing of Terrestrial Ecosystems
Vegetation growth in cold regions is limited by multiple environmental constraints including low temperatures, frozen sub-surface soils in permafrost affected terrain, low light levels in winter and shoulder seasons, water stress in summer, and nutrient limitations due to the very slow release of plant available nitrogen and phosphorus from seasonally frozen or inundated soils [196,197].Although ground-based measurements are needed to inform local-scale investigations, remote sensing is especially important for the cold regions because it provides spatially and temporally resolved data over broader scales, allowing for synoptic investigations of vegetation ranging from individual plots to regional and global extents.

Vegetation Mapping
Land cover type is a general term encompassing a range of important information about ecosystems, including biotic and abiotic properties related to vegetation, energy balance, and carbon exchanges.Satellite remote sensing has been used throughout the cold regions to provide various maps of land cover that group vegetation according to geobotanical themes, including physiology.The spatial patterns in the vegetation maps can provide useful insight into how microclimates, soil type, and hydrology, disturbance and plant succession contribute to variations in plant community characteristics at landscape to regional levels.These maps are also used to parameterize ecosystem process models by means of parameter look-up tables aggregated according to generalized plant functional types [198,199].
Classification algorithms trained on satellite optical-IR spectral information have been used to map land cover and vegetation type over large areas.For example, the United States Geological Survey (USGS) provides 30-m land cover data for the state of Alaska for years 2001 and 2011, using the C5 decision-tree classifier applied to Landsat TM and ETM+ (Enhanced Thematic Mapper) imagery (http://www.mrlc.gov).For the Anderson Level II classification (19 classes, defining different types of forest, shrub, herbaceous, and wetland), the overall accuracy ranged from 59% to 76%, depending on the definition of agreement with reference data [200].The Earth Observation for Sustainable Development of Forests (EOSD) project used an unsupervised k-means clustering approach with Landsat ETM+ data to map land cover across Canada at 25-m resolution for the year 2000, using a detailed 23-class system [201].In another study, phenological data derived from Landsat 8 NDVI timeseries, along with other inputs, were used to classify land cover across ice-free portions of Greenland at 30-m resolution [202].Mapped classes included fen, dry heath and grasslands, wet heath, and copse and tall shrubs, with an overall classification accuracy of 89%.Other studies have used higher resolution imagery to map vegetation communities in heterogenous tundra landscapes (e.g., using IKONOS imagery [203]).For sites across the North Slope of Alaska, tundra vegetation communities were shown to be separable based on visible wavelengths collected through field spectroscopy, and vegetation type was mapped at 2-m resolution with ~70% accuracy using linear discriminant analysis applied to WorldView-2 data [40].
In addition to optical-IR data, airborne and satellite SAR data can be effective for separating land cover classes due to its sensitivity to vegetation structure and water content [204][205][206].A supervised classification approach was applied to airborne, multifrequency polarimetric SAR imagery to map functional vegetation types at 30-m resolution across a boreal forest area [204].Classes such as jack pine, black spruce, and trembling aspen were mapped with high accuracy (> 90%).Radar data can also be used to distinguish different types of wetland vegetation, as higher frequency microwave data can detect flooded short vegetation (e.g., fens and bogs), while lower frequency data are better at detecting flooded tall vegetation (e.g., forests and swamps) [206].The fusion of microwave and optical data through traditional supervised classification (e.g., maximum-likelihood classifier) [205] and machine-learning approaches [206,207] enables improved land cover classification over the high-latitudes.Deep-learning approaches also show promise in utilizing semantic information of satellite imagery for classifying complex ecosystems [208].

Vegetation Growth and Photosynthetic Carbon Assimilation
A number of remote sensing vegetation indices such as NDVI, leaf area index (LAI), and enhanced vegetation index (EVI) have been developed to characterize vegetation properties related to photosynthesis on a per-pixel basis [209].In Arctic-boreal regions, multispectral satellite data have been used to quantify variables related to vegetation growth [210], biomass [211,212], and carbon fluxes [213].Satellite microwave systems can also provide information about vegetation growth and carbon assimilation.For example, vegetation optical depth (VOD) derived from satellite microwave Tb observations is sensitive to vegetation water content and provides information on both canopy biomass (photosynthetic and non-photosynthetic) phenology and drought stress [214][215][216][217][218]. VOD has also been used to monitor vegetation growth and recovery after fires in boreal forests [219].
Much of the remote sensing-based work for plant growth and carbon cycling has focused on modeling gross primary productivity (GPP), which represents carbon biomass created by plants through the process of photosynthesis over a given length of time and space (e.g., m 2 day −1 ).GPP also represents the amount of atmospheric CO 2 sequestered by plants within biomass.GPP models are often based on a light-use efficiency (LUE) framework, which estimates GPP as a function of APAR, the fraction of absorbed photosynthetically active radiation (PAR), and a photosynthetic efficiency parameter which describes the rate at which absorbed radiation is used for carbon fixation [220].To model APAR, inputs of LAI and fPAR (the fraction of canopy absorbed PAR) are needed, which can be modeled using spectral reflectance data combined with radiative transfer algorithms [221].Alternatively, NDVI can be used as a proxy for fPAR [222].The efficiency parameter is more difficult to model and can vary based on vegetation type, water limitations, and light conditions [220].However, GPP has also been shown to be directly related to the EVI, which tends to be less saturated than NDVI over higher canopy densities and less sensitive to soil background noise [220].Relationships between EVI and GPP were shown to vary among sites, with correlations generally stronger for deciduous sites than for evergreen sites [220].Another important and newer proxy for GPP is solar-induced fluorescence (SIF) [223].SIF quantifies the amount of light reemitted by chlorophyll molecules as a byproduct of photosynthesis, and has been shown to be directly proportional to GPP [224].Satellite-based SIF (e.g., from the Global Ozone Monitoring Experiment 2 [225] and Orbiting Carbon Observatory-2 [226]), allows for global monitoring of terrestrial GPP and the carbon cycle.However, the SIF signal has relatively small magnitude and generally requires a large sensor footprint and coarse temporal compositing of the data to obtain an adequate signal-to-noise ratio relative to NDVI and EVI observations [227].
Aboveground biomass (AGB) is an important component of global carbon accounting.Remote sensing techniques are used to estimate AGB by capturing both vertical structure and spatial variability of canopies.For example, AGB has been estimated from LIDAR-derived canopy height and vertical structure metrics [228,229].Biomass was mapped for the circumboreal zone using multi-step modeling that combined field-based biomass data, airborne and satellite LIDAR data [230,231], and similarly at the regional scale for boreal forest in Québec [232].Other studies demonstrated the potential of low-frequency (e.g., L-and P-band) polarimetric SAR/InSAR data for estimating AGB in boreal forests and over complex terrain [233][234][235][236][237].

Changes and Trends
Climate warming in cold regions has been altering the phenology of snow, lake, and river ice, and vegetation, causing accelerated ice melting in polar regions, and leading to complex wetting and drying trends, and greening and browning patterns in northern high latitudes and TP.Multi-decade remote sensing observations are essential in documenting the environmental changes and revealing the long-term trends  The changes in the annual non-frozen season length reflect the overall warming trend in the climate system.The timing and duration of the non-frozen season is an important factor affecting water, carbon, and energy budgets in cold land areas.Recent trends toward earlier and longer nonfrozen seasons coincide with global warming and have been shown to be a major driver of northern vegetation greening, active layer deepening and permafrost degradation, enhanced evapotranspiration, earlier snowmelt onset, and associated changes in terrestrial water and energy budgets [21,34,149].For example, satellite observations indicate that the snow end date in spring advanced by 5.11 days from 2001 to 2014 in the high northern latitudes (52-75°N) [238], along with shorter lake ice cover duration at higher latitudes from 2002 to 2015 [17].The melting of permafrost ice provides the water supply to thermokarst lakes [239] and leads to surface water expansion detected within continuous and discontinuous permafrost zones [29].
Satellite remote sensing observations have been used to detect changes in warming Arcticboreal ecosystems.Increasing temperatures and atmospheric CO2 concentrations can impact the production, dynamics, and composition of vegetation, as well as soil moisture and other soil properties [240].One of the most pronounced ongoing changes is tundra shrub expansion [203].Changes to vegetation are likely to occur around vegetation ecotones [203,240].Much effort has been given to identifying regions of vegetation production increase (greening) or decrease (browning).A 10-year time series of monthly average AVHRR NDVI indicated an increase in photosynthetic activity in the northern high latitudes from 1981 to 1991 driven by earlier growing season onset [241].A subsequent study using a longer AVHRR time series indicated an increase in tundra vegetation The changes in the annual non-frozen season length reflect the overall warming trend in the climate system.The timing and duration of the non-frozen season is an important factor affecting water, carbon, and energy budgets in cold land areas.Recent trends toward earlier and longer non-frozen seasons coincide with global warming and have been shown to be a major driver of northern vegetation greening, active layer deepening and permafrost degradation, enhanced evapotranspiration, earlier snowmelt onset, and associated changes in terrestrial water and energy budgets [21,34,149].For example, satellite observations indicate that the snow end date in spring advanced by 5.11 days from 2001 to 2014 in the high northern latitudes (52-75 • N) [238], along with shorter lake ice cover duration at higher latitudes from 2002 to 2015 [17].The melting of permafrost ice provides the water supply to thermokarst lakes [239] and leads to surface water expansion detected within continuous and discontinuous permafrost zones [29].
Satellite remote sensing observations have been used to detect changes in warming Arctic-boreal ecosystems.Increasing temperatures and atmospheric CO 2 concentrations can impact the production, dynamics, and composition of vegetation, as well as soil moisture and other soil properties [240].One of the most pronounced ongoing changes is tundra shrub expansion [203].Changes to vegetation are likely to occur around vegetation ecotones [203,240].Much effort has been given to identifying regions of vegetation production increase (greening) or decrease (browning).A 10-year time series of monthly average AVHRR NDVI indicated an increase in photosynthetic activity in the northern high latitudes from 1981 to 1991 driven by earlier growing season onset [241].A subsequent study using a longer AVHRR time series indicated an increase in tundra vegetation growth in boreal North America, largely driven by longer growing seasons, but decreasing growth in interior forests [242].Recent work used time series of Landsat vegetation indices in Arctic sites, with results indicating more areas with increasing growth than decreasing growth [243].Analysis of a 28-year Landsat NDVI record indicated that greening and browning of Canadian boreal forests were largely driven by disturbances from wild fire and insect damages; and, in forests not affected by disturbance, climate changes were associated with both areas of greening and browning [240].Regarding phenological changes, a 30-year Landsat record was used to detect an earlier/heterogeneous leaf emergence trend in temperate/boreal deciduous forests [197], while a 33-year AVHRR NDVI time series revealed regional trends toward earlier growing season onset, later growing season end, and longer growing season duration over the high northern latitudes [31].Based on recent satellite optical and microwave observations for years 2003 through 2017 over the high northern latitudes, the mean summer (JJA) MODIS NDVI record (MYD13A1.006;[244]) showed similar spatial patterns with the AMSR-E/2 VOD record [216], despite the different spatial resolutions and underlying physics of the observations (Figure 5a,b).Major greening and biomass growth trends are found in northern taiga and tundra regions, where both NDVI and VOD show significant increases (p-value < 0.05) (Figure 5c,d).However, declining NDVI and VOD trends are also widespread, indicating decreasing productivity.The declining trend areas largely occur in boreal forest but are generally less significant than positive trend areas.

Antarctic and Greenland Ice
The Antarctic and Greenland ice sheets are the largest ice bodies on Earth and are significantly affected by changing air temperatures and solar radiation.If completely melted, the polar ice sheets

Antarctic and Greenland Ice
The Antarctic and Greenland ice sheets are the largest ice bodies on Earth and are significantly affected by changing air temperatures and solar radiation.If completely melted, the polar ice sheets would raise global sea level by 70 m [245].Recent studies show that the magnitude of recent melting of the Greenland ice sheet is exceptional over at least the last 350 years [246].Greenland's ice is melting so fast that it could become a major contributor to sea-level rise within two decades.Greenland ice loss mainly occurred in the southeast and northwest margins of the ice sheet in the 2000-2010 period; while the largest sustained acceleration (~10 years) in ice loss was detected in southwest Greenland from GRACE (Gravity Recovery and Climate Experiment) observations [247].The overall transformation of ice into liquid water appears to be accelerating and Greenland loses an average of 270 billion tons of ice each year [248].
The recent loss of continental ice includes both the northern and southern hemispheres.An estimate of the mass balance of the entire Antarctic ice sheet over a 25-year record (1992 to 2017) shows that the Antarctic Peninsula, the smallest ice sheet in Antarctica, has lost an average of 20 Gigatonnes (Gt) of ice per year.The loss rate increased during the study period especially after the year 2000.The West Antarctic Ice Sheet lost 53 ± 29 Gt yr −1 from 1992-1997, and the loss rate accelerated to 159 ± 26 Gt yr −1 from 2012-2017.The East Antarctic Ice Sheet is relatively stable, with small gains over the study period [249].The changes of polar mass balance are associated with snow and ice surface darkening [250], warmer atmosphere and ice surface conditions [251][252][253], and increased surface melt duration and extent [253][254][255] as observed by satellite optical and microwave sensors.
The velocity of ice flow in the Antarctic has been closely monitored using optical and radar remote sensing due to its importance in determining ice discharge and sea level rise [60,256].The velocity map derived from the MODIS-based Mosaic of Antarctica data showed the general flow patterns of glaciers and ice sheets moving from interior Antarctica toward the ocean for the periods from 2003-2004 and 2008-2009 (Figure 6).Continuous monitoring of the widespread ice flow over the entire continent is highly needed for improving our understanding of ice sheet dynamics and evolution in a warming climate [60,256].
Remote Sens. 2019, 10, x FOR PEER REVIEW 15 of 36 year 2000.The West Antarctic Ice Sheet lost 53±29 Gt yr −1 from 1992-1997, and the loss rate accelerated to 159±26 Gt yr −1 from 2012-2017.The East Antarctic Ice Sheet is relatively stable, with small gains over the study period [249].The changes of polar mass balance are associated with snow and ice surface darkening [250], warmer atmosphere and ice surface conditions [251][252][253], and increased surface melt duration and extent [253][254][255] as observed by satellite optical and microwave sensors.
The velocity of ice flow in the Antarctic has been closely monitored using optical and radar remote sensing due to its importance in determining ice discharge and sea level rise [60,256].The velocity map derived from the MODIS-based Mosaic of Antarctica data showed the general flow patterns of glaciers and ice sheets moving from interior Antarctica toward the ocean for the periods from 2003-2004 and 2008-2009 (Figure 6).Continuous monitoring of the widespread ice flow over the entire continent is highly needed for improving our understanding of ice sheet dynamics and evolution in a warming climate [60,256].

Tibetan Plateau
The Tibetan Plateau (TP), the most extensive highland in the world, has an area of approximately 2.5 × 106 km 2 and an average elevation of over 4000 m.The TP also has the Earth's largest storage of ice outside of the north and south polar regions.The climate of the TP is changing rapidly with temperatures warming at a rate of around 0.36 °C/decade [3], which is twice the mean global trend

Tibetan Plateau
The Tibetan Plateau (TP), the most extensive highland in the world, has an area of approximately 2.5 × 106 km 2 and an average elevation of over 4000 m.The TP also has the Earth's largest storage of ice outside of the north and south polar regions.The climate of the TP is changing rapidly with temperatures warming at a rate of around 0.36 • C/decade [3], which is twice the mean global trend [4].Consequently, the impacts of climate change on the TP environment are most pronounced, leading to earlier onset of seasonal thawing, accelerating glacier melting, permafrost degradation, and complex changes in snow, lakes, and vegetation.
Glacier melting has been observed over the TP and larger High Mountain Asia region.Kaab et al.Global and localized satellite snow products have been used for studying environmental changes in the TP, including snow impacts on the regional water cycle, ecosystems, and atmospheric circulation [261].
Based on the MODIS snow cover product (MOD10A2) and observations from 37 meteorological stations, a significant trend of earlier onset of snow ablation during the 2001-2015 period was detected over the TP [262].By analyzing MODIS snow cover data from 2001 to 2011, other research [263] found about 34.14% (5.56%) of the TP area having a declining (significant declining) trend in snow duration, while 24.75% (3.9%) of the region showed increasing snow duration.To further enhance the accuracy of TP snow products, Chen et al. [264] integrated snow cover data from multiple sources to generate a gap-filled daily 5-km Tibetan Plateau snow cover extent record (TPSCE) from 1981-2016.As revealed by the TPSCE, the snow cover fraction increased in the northern interior TP river basins and upper reaches of the Yangtze, Mekong, and Brahmaputra River basins (Figure 7).Global and localized satellite snow products have been used for studying environmental changes in the TP, including snow impacts on the regional water cycle, ecosystems, and atmospheric circulation [261].
Based on the MODIS snow cover product (MOD10A2) and observations from 37 meteorological stations, a significant trend of earlier onset of snow ablation during the 2001-2015 period was detected over the TP [262].By analyzing MODIS snow cover data from 2001 to 2011, other research [263] found about 34.14% (5.56%) of the TP area having a declining (significant declining) trend in snow duration, while 24.75% (3.9%) of the region showed increasing snow duration.To further enhance the accuracy of TP snow products, Chen et al. [264] integrated snow cover data from multiple sources to generate a gap-filled daily 5-km Tibetan Plateau snow cover extent record (TPSCE) from 1981-2016.As revealed by the TPSCE, the snow cover fraction increased in the northern interior TP river basins and upper reaches of the Yangtze, Mekong, and Brahmaputra River basins (Figure 7).Beside the global FT products available from satellite microwave observations [160], Kou et al. [265] developed an enhanced resolution (0.05 degree) FT product over the TP by merging MODIS LST with AMSR-E Tb records.Li et al. [266] analyzed changes in the soil FT cycle over the TP using SSM/I data from 1988-2007.They identified a trend toward earlier onset of soil thaw by approximately 14 days/decade and later onset of soil freeze by approximately 10 days/decade.The observed changes in FT patterns over the TP are also closely related to regional climate warming [265].
There are more than 1000 lakes with an area greater than 1 km 2 on the TP.Generally, the total TP lake area is expanding, from 41,800 km 2 in 2005-2006 [267] to 46,600 km 2 in 2015 [268].The lake expansion can be directly measured through optical satellite remote sensing.Wang et al. [269] analyzed the trend of lake area changes during 1960-2000 by integrating aerial photos, satellite images (Landsat and CBERS-1 (China-Brazil Earth Resources Satellite)), and topography information.They found that most lakes in the central TP were expanding, while the lakes in the source regions of the Yellow River were shrinking.For lakes in the central TP, Wan et al. [270] analyzed Landsat TM/ETM+ and CBERS images, and found that lakes southeast of the Qiangtang area were expanding from 1975 to 2005.Another study [271] analyzed Landsat images from 1970 to 2010 and found a shrinking lake trend in the southwest TP contrasting with rapid lake expansion in the northeast TP.Yang et al. [272] investigated lake extent fluctuations in the Hindu Kush-Himalaya-Tibetan (HKHT) regions over the past 40 years using Landsat images obtained from the 1970s to 2014.They showed that the TP lake trends are distinct from region to region, with the most intensive lake shrinking observed in northeastern HKHT (HKHT Interior, Tarim, Yellow, Yangtze), while the most extensive expansion was observed in the western and southwestern HKHT (Amu Darya, Ganges Indus, and Brahmaputra), largely caused by the proliferation of small lakes in high-altitude regions from the 1970s to 1995.Lake water levels measured from LIDAR/Radar altimetry can also be used to quantify lake changes.Based on 10-year TOPEX/Poseidon altimetry data, the water level of La'nga lake in the western TP was found to decrease steadily from 1993 to 2001, while that of Ngangzi lake in the eastern TP decreased from 1993 to 1997 and then increased monotonically afterwards [273].Another study [274] used ICESat altimetry data to analyze water level changes from 2003-2009 for 74 TP lakes and identified an increasing lake level trend (~0.23 m yr −1 ) over 84% of the lakes represented.Consistent with the findings of [275], an average water level increase of 0.20 m yr −1 was detected from GLAS data over 154 TP lakes for the same period [276].By integrating multi-altimeter data from Envisat/RA-2, Cryosat-2/Siral, Jason-1/Poseidon-2, and Jason-2/Poseidon-3, Gao et al. [277] found that water levels increased by about 0.275 m yr −1 for over 82% of 51 TP lakes sampled, while major lake expansion and shrinking were identified over the northern and southern TP, respectively.
A general greening trend from the 1980s to 2010s was detected from satellite remote sensing over the TP [278].Seasonal analysis based on GIMMS (Global Inventory Monitoring and Modeling System) NDVI data further revealed that the largest NDVI increase occurred in autumn over 61% of the TP, while the smallest increase occurred in spring over 41% of the region [279].Vegetation changes in the northeast, southwest, mid-eastern, and southern TP regions were driven by three different factors, including changes in surface air temperature, water availability, and solar radiation, respectively.Anthropogenic disturbances may offset climate-driven vegetation greening and exacerbate vegetation browning, while ecosystem conservation efforts contributed to vegetation recovery in the TP Three-River Headwaters Region [280].
Based on the NDVI data from 1982 to 2014, an advancing start of growing season, delayed end of season and increasing length of growing season were identified for meadow areas in the eastern TP; while the opposite changes were found for the steppe and sparse herbaceous or sparse shrub areas in the northwest and western edges of the TP.The satellite-observed phenology changes were driven by a number of environmental factors including temperature [281][282][283], precipitation [284], sunshine duration [285], and snow cover [286]; and may be partially attributed to aerosol contamination in the satellite observations [287].

Limitations of Current Approaches
Despite great achievements in cold land remote sensing, comprehensive assessment of both long-term trends and relatively abrupt environmental changes requires quantifying hydrological and ecological variables with greater precision, including clearer delineation of different land components, better accuracy, data consistency, spatial resolution, and temporal sampling.
Ice-The spatial resolution of satellite gravimetry observations of glaciers and ice sheets are generally very coarse (around 100 km), which limits their ability to precisely detect and locate glacier mass change at finer scales.LIDAR-based measurements of glacier elevations are of high accuracy (decimeter level) but with relatively poor spatial coverage, which consequently requires data interpolation and extrapolation for applications over large areas.For monitoring glacier movement, feature tracking fails when low coherence occurs in fast moving glaciers or over prolonged time intervals between observations.Snow-Snow covered area is mainly derived using optical sensors onboard polar-orbiting or geostationary satellites.However, it is still challenging to distinguish between clouds and snow in satellite optical snow mapping.Long-term (>40-year) SWE data records have been generated for the globe through satellite passive microwave sensors, but the retrieval spatial resolution and accuracy has generally not met the requirements for regional climate, numerical weather prediction, and hydrological research.The retrieval uncertainties mainly come from four sources, including mixed pixel effects, terrain effects, large diversity in snow physical properties, and low microwave sensitivity to shallow snow.
Landscape Freeze/Thaw States-Most of the available global FT products represent aggregate landscape FT conditions that do not distinguish land components of the FT signal within a satellite footprint at several-10s km.Accurate estimates of FT metrics (e.g., spring thaw timing, frost days) pertaining to soil, snow, and vegetation components isolated from the integrated microwave FT signals with improved spatial resolution are required for better understanding of land and atmosphere interactions, including carbon, water, and energy exchanges.It is also needed to understand the scaling effects for heterogeneous terrain to bridge multi-resolution satellite observations [38].
Water Bodies-It is challenging to detect water within mixed pixels or under overlying vegetation for both optical-IR and microwave sensors.In addition, near real-time and fine-scale mapping of regional SW dynamics are urgently needed for monitoring flood hazards and arctic wetting and drying trends.
Vegetation-The short growing season and generally low vegetation productivity of Arctic-boreal regions can create challenges in terms of detecting changes in vegetation growth through use of remote sensing time series data.Other challenges in Arctic-boreal regions include saturated satellite signals over high-biomass vegetation, high occurrence of shallow water bodies, the presence of snow and ice cover, and spatially heterogenous snowmelt resulting in pixels mixed with snow and vegetation.Additionally, long-term and fine-scale delineation (1-200 m) of vegetation (e.g., shrub) patch distributions is lacking, but greatly needed in high-latitude ecosystem studies [288].
It is anticipated that SAR will play a greater role in cold land studies.The European Space Agency's Sentinel-1A/B satellites and the Canadian Space Agency's RADARSAT Constellation Mission (RCM) launched on 12 June 2019 [293] can provide multi-resolution (1-3 m to 100-m) C-band SAR and InSAR measurements every 1-6 days globally and at low cost for monitoring glacier movement, river ice conditions, water body dynamics, snow, and vegetation parameters [298].The next generation NISAR (National Aeronautics and Space Administration -Indian Space Research Organisation SAR) mission is expected to improve capabilities for dynamic mapping of global surface water at resolutions from 3-10 m [294] and provide new opportunities to derive SWE in complex terrain from singleand repeat-pass radar (L-band) interferometry [16].The surface water ocean topography (SWOT) mission has a projected launch in 2021 and will enable estimation of the changing volumes of global water bodies whose surface area exceeds 250 m by 250 m at sub-monthly time scales [178,291].The upcoming BIOMASS mission (launch 2021) offers opportunities for mapping vegetation biomass with enhanced P-band penetration capability [299], though National restrictions on P-band transmissions will eliminate BIOMASS coverage over North America [300].Another future direction for analyzing environmental changes in remote areas under harsh climatic conditions could also involve relatively low-cost multi-temporal remote sensing surveys.Unmanned aerial vehicle (UAV)-based remote sensing (e.g., thermal imaging, structure from motion photogrammetric techniques), given relatively low application barriers (e.g., low cost for repeat seasonal or yearly surveys), is expected to become a strategic tool for better monitoring and understanding of the geomorphologic dynamics of cold regions, and related impacts on ecosystems and infrastructures.
Improved retrieval algorithms will better leverage the remote sensing observations.For example, entropy-based multi-scale image matching and optical flow techniques potentially help overcome the limitations of conventional image matching approaches for monitoring fast moving glaciers [301,302].Backward reconstruction techniques using a temperature-index or energy-balance model provide another promising tool to estimate SWE through the melt season in mountainous regions and elsewhere [303,304].There is also great potential for improving snow depth and SWE retrievals using LIDAR measurements [305,306], microwave interferometry [147,307], and GNSS techniques [308], which largely avoid issues related to snow microstructure.Considering the complementarity of LIDAR, optical-IR, and microwave remote sensing, the use of data integration techniques such as deep-learning approaches [208] from a collection of multi-sensor observations may enable enhanced delineations of snow, water, soil, and vegetation elements, and accurate mapping of environmental variables at high spatial-temporal resolutions.

Conclusions
The rapid environmental changes in cold land regions have profound impacts on ecosystems, geomorphology, animal habitats, and human lives.Remote sensing is the most valuable and indispensable technique for large-scale monitoring of such changes, accurately quantifying both transient anomalies and longer-term trends, while providing observational benchmarks to test earth system model projections.Calibrated long-term observations from these sensors have produced relatively precise satellite data records documenting significant changes in landscape FT states, snow extent and depth, glacier mass and movement, water body dynamics, and lake ice [309] and vegetation phenology.Earlier onset of snow melt [310], soil thaw, and lake ice break-up, longer potential growing seasons, and expanding lakes were identified in both the northern high latitudes and TP; and continuous ice loss has been observed in both Greenland and Antarctica.Despite these achievements, pressing issues such as melting permafrost, disappearing glaciers, shrinking ice cover, and structural and functional changes in ecosystems require more timely and accurate interpretations from remote sensing observations.Multi-sensor data fusion approaches show promise in overcoming the drawbacks of single sensor observations, while strengthening capabilities for monitoring and detecting environmental changes in cold regions.Next generation satellite missions including SWOT, NISAR, and BIOMASS; emerging techniques such as micro-satellites and data mining; and coordinated research activities [311][312][313] will enable cold land mapping with unprecedented sampling frequency, spatial resolution, and accuracy.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/11/16/1952/s1,Table S1.Overview of satellite missions/sensors for cold land remote sensing and summarized in the review, Table S2.Vertical accuracy of DEMs from different sensors and methods over glacier surface.
Author Contributions: All authors contributed equally to the manuscript.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Figure 1 .
Figure 1.Components of snow emissions and scattering observed by space-borne microwave radiometer and radar observations.By neglecting atmosphere scattering and upward emission, the satellite signals mainly consist of contributions from the air-snow surface (as), snow pack (v), underlying soil (g), and snow-soil interactions (gv) .

Figure 1 .
Figure 1.Components of snow emissions and scattering observed by space-borne microwave radiometer and radar observations.By neglecting atmosphere scattering and upward emission, the satellite signals mainly consist of contributions from the air-snow surface (as), snow pack (v), underlying soil (g), and snow-soil interactions (gv).

Figure 3 .
Figure 3. Projected infrastructure hazard areas from degrading permafrost over the Northern Hemisphere by year 2050 [26].The figure was reproduced with permission from [26], which is licensed under a Creative Commons Attribution 4.0 International License.

Figure 3 .
Figure 3. Projected infrastructure hazard areas from degrading permafrost over the Northern Hemisphere by year 2050 [26].The figure was reproduced with permission from [26], which is licensed under a Creative Commons Attribution 4.0 International License.

3. 1 .
Northern High Latitudes The long-term (1979-2017) anomalies of the annual non-frozen season derived from the satellite microwave FT observational record show a lengthening non-frozen season trend over vegetated lands (excluding large water bodies and permanent ice/snow covered areas) in the high-northern latitudes (3.30 day decade −1 ; p-value <0.001), with similar trends over the Northern Hemisphere (3.98 day decade −1 ), Southern Hemisphere (3.61 day decade −1 ), and global domains (3.93 day decade −1 ) as shown in Figure 4. Remote Sens. 2019, 10, x FOR PEER REVIEW 13 of 36

Figure 4 .
Figure 4. Non-frozen season trend (1979-2017) for the Northern Hemisphere (NH), Southern Hemisphere (SH), high-northern latitudes (HNL), and global domain.The anomalies were calculated as annual differences from the long-term mean.Grid cells dominated by permanent snow/ice cover and large water bodies (open water fraction ≥ 20%) are excluded from the analysis.

Figure 4 .
Figure 4. Non-frozen season trend (1979-2017) for the Northern Hemisphere (NH), Southern Hemisphere (SH), high-northern latitudes (HNL), and global domain.The anomalies were calculated as annual differences from the long-term mean.Grid cells dominated by permanent snow/ice cover and large water bodies (open water fraction ≥ 20%) are excluded from the analysis.
Remote Sens. 2019, 10, x FOR PEER REVIEW 14 of 36 similar spatial patterns with the AMSR-E/2 VOD record [216], despite the different spatial resolutions and underlying physics of the observations (Figure 5a,b).Major greening and biomass growth trends are found in northern taiga and tundra regions, where both NDVI and VOD show significant increases (p-value < 0.05) (Figure 5c,d).However, declining NDVI and VOD trends are also widespread, indicating decreasing productivity.The declining trend areas largely occur in boreal forest but are generally less significant than positive trend areas.

Figure 5 .
Figure 5. Vegetation conditions and growth trend over the high northern latitudes from 2003 through 2017.(a) Mean summer MODIS NDVI; (b) Mean summer AMSR-E/2 VOD; (c) Trends in mean summer NDVI; (d) Trends in mean summer VOD.In c and d, positive values indicate an increasing trend and negative values indicate a decreasing trend; pixels with significant trend (p-value < 0.05) are shown with crosshatch.

Figure 5 .
Figure 5. Vegetation conditions and growth trend over the high northern latitudes from 2003 through 2017.(a) Mean summer MODIS NDVI; (b) Mean summer AMSR-E/2 VOD; (c) Trends in mean summer NDVI; (d) Trends in mean summer VOD.In c and d, positive values indicate an increasing trend and negative values indicate a decreasing trend; pixels with significant trend (p-value < 0.05) are shown with crosshatch.

Figure 6 .
Figure 6.Antarctic surface ice velocity map.The floating velocity field of the Drygalski Ice Tongue is enlarged in the close-up below.The figure was reproduced with permission from [60].

Figure 6 .
Figure 6.Antarctic surface ice velocity map.The floating velocity field of the Drygalski Ice Tongue is enlarged in the close-up below.The figure was reproduced with permission from [60].
[257] used ICESat to analyze glacier mass change in the Hindu Kush-Karakoram-Himalaya region during 2003-2008 and found a mass loss rate of −12.8 ± 3.5 Gt yr −1 in this region, which is faster than the rate previously estimated using GRACE [258].For the whole TP region, Neckel et al. [259] estimated an overall mass loss rate of −15.6 ± 10.1 Gt yr −1 using ICESat observations from 2003 to 2009.For the same period, Gardner et al. [260] estimated a total mass change of -29 ± 13 Gt yr −1 over High Mountain Asia by integrating GRACE and ICESat observations.

Figure 7 .
Figure 7. Climatology and changes of the TP snow cover fraction (%) calculated using MCD10A1-TP (a, b) and TPSCE (c, d) from 2001-2014.Black dots in (b) and (d) indicate the changes are statistically significant at the 95% level.The figure was reproduced with permission from [264].

Figure 7 .
Figure 7. Climatology and changes of the TP snow cover fraction (%) calculated using MCD10A1-TP (a, b) and TPSCE (c, d) from 2001-2014.Black dots in (b) and (d) indicate the changes are statistically significant at the 95% level.The figure was reproduced with permission from [264].