Articles | Volume 15, issue 6
The Cryosphere, 15, 2623–2646, 2021
The Cryosphere, 15, 2623–2646, 2021

Research article 14 Jun 2021

Research article | 14 Jun 2021

Surface melting over the Greenland ice sheet derived from enhanced resolution passive microwave brightness temperatures (1979–2019)

Surface melting over the Greenland ice sheet derived from enhanced resolution passive microwave brightness temperatures (1979–2019)
Paolo Colosio1, Marco Tedesco2,3, Roberto Ranzi1, and Xavier Fettweis4 Paolo Colosio et al.
  • 1DICATAM, Università degli Studi di Brescia, Brescia 25121, Italy
  • 2Lamont-Doherty Earth Observatory, Columbia University, New York 10964, USA
  • 3NASA GISS, New York 10025, USA
  • 4Department of Geography, University of Liège, Liège 4000, Belgium

Correspondence: Paolo Colosio (


Surface melting is a major component of the Greenland ice sheet surface mass balance, and it affects sea level rise through direct runoff and the modulation of ice dynamics and hydrological processes, supraglacially, englacially and subglacially. Passive microwave (PMW) brightness temperature observations are of paramount importance in studying the spatial and temporal evolution of surface melting due to their long temporal coverage (1979–present) and high temporal resolution (daily). However, a major limitation of PMW datasets has been the relatively coarse spatial resolution, which has historically been of the order of tens of kilometers. Here, we use a newly released PMW dataset (37 GHz, horizontal polarization) made available through a NASA “Making Earth System Data Records for Use in Research Environments” (MeASUREs) program to study the spatiotemporal evolution of surface melting over the Greenland ice sheet at an enhanced spatial resolution of 3.125 km. We assess the outputs of different detection algorithms using data collected by automatic weather stations (AWSs) and the outputs of the Modèle Atmosphérique Régional (MAR) regional climate model. We found that sporadic melting is well captured using a dynamic algorithm based on the outputs of the Microwave Emission Model of Layered Snowpack (MEMLS), whereas a fixed threshold of 245 K is capable of detecting persistent melt. Our results indicate that, during the reference period from 1979 to 2019 (from 1988 to 2019), surface melting over the ice sheet increased in terms of both duration, up to 4.5 (2.9) d per decade, and extension, up to 6.9 % (3.6 %) of the entire ice sheet surface extent per decade, according to the MEMLS algorithm. Furthermore, the melting season started up to 4.0 (2.5) d earlier and ended 7.0 (3.9) d later per decade. We also explored the information content of the enhanced-resolution dataset with respect to the one at 25 km and MAR outputs using a semi-variogram approach. We found that the enhanced product is more sensitive to local-scale processes, thereby confirming the potential of this new enhanced product for monitoring surface melting over Greenland at a higher spatial resolution than the historical products and for monitoring its impact on sea level rise. This offers the opportunity to improve our understanding of the processes driving melting, to validate modeled melt extent at high resolution and, potentially, to assimilate these data in climate models.

1 Introduction

The Greenland ice sheet is the largest ice mass in the Northern Hemisphere with a glaciated surface area of about 1 800 000 km2, a thickness of up to 3 km and a stored water volume of about 2 900 000 km3, which is enough to raise the mean sea level by about 7.2 m (Aschwanden et al., 2019). In this regard, estimating mass losses from Greenland is crucial in order to better understand climate system variability and the contribution of Greenland to current and future sea level rise (Mouginot et al., 2019). According to data from the Gravity Recovery and Climate Experiment (GRACE) satellite mission, which records changes in the Earth's gravitational field, Greenland lost mass at an average rate of 278 ± 11 Gt yr−1 between 2002 and 2016 (IPCC, 2019), contributing to a sea level rise of 7.9 mm per decade. The contribution of the Greenland ice sheet to sea level rise also accelerated at a rate of 21.9 ± 1 Gt yr−2 over the period from 1992 to 2010 (Rignot et al., 2011), indicating that monitoring the Greenland ice sheet and the Antarctic ice sheet is crucial to assess the impact of global warming on sea level rise and the global water balance (Kargel et al., 2005, 2014; Le Meur et al., 2018).

Mass can be lost through surface (e.g., runoff) and dynamic (e.g., calving) processes with total mass loss roughly split in half between the two (Flowers, 2018). Among the processes influencing the surface mass balance, i.e., the difference between accumulation (Frezzotti et al., 2007) and ablation, surface melting plays a crucial role, affecting direct loss through the export of surface meltwater to the surrounding oceans and though feedbacks between supraglacial, englacial, and subglacial processes and their influence on ice dynamics (e.g., Fettweis et al., 2005, 2011, 2017; van den Broeke et al., 2016; Alexander et al., 2016).

Passive microwave (PMW) brightness temperatures (Tb) are a crucial tool for studying the evolution of surface melting over the Greenland and Antarctica ice sheets (i.e., Jezek et al., 1993; Steffen et al., 1993; Abdalati et al., 1995; Tedesco, 2007; Tedesco et al., 2009; Tedesco, 2009, 2014; Fettweis et al., 2011). PMW-based algorithms are founded on the fact that the emission of a layer of dry snow in the microwave region is dominated by volume scattering (e.g., Macelloni et al., 2001); as snow melts, the presence of liquid water within the snowpack increases the imaginary part of the electromagnetic permittivity by several orders of magnitude with respect to dry snow conditions, with the ultimate effect of considerably increasing Tb (Ulaby et al., 1986; Hallikainen et al., 1987), as shown by in situ measurement campaigns (e.g., Cagnati et al., 2004). Because of the large difference between dry and wet snow emissivity, even relatively small amounts of liquid water have a dramatic effect on the Tb values (e.g., Tedesco, 2009), making PMW data extremely suitable for mapping the extent and duration of melting at large spatial scales and high temporal resolution (in view of their insensitivity to atmospheric conditions at the low frequencies of the microwave spectrum). Consequently, PMW data have been widely adopted in melt detection studies, and different remote sensing techniques have been proposed in the literature (e.g., Steffen et al., 1993; Abdalati and Steffen, 1995; Joshi et al., 2001; Liu et al., 2005; Ashcraft and Long, 2006; Macelloni et al., 2007; Tedesco et al., 2007; Kouki et al., 2019; Tedesco and Fettweis, 2020).

The capability of PMW sensors to collect useful data during both day and night and under all weather conditions allows surface melt mapping at a high temporal resolution (at least twice a day over most of the Earth). PMW Tb records are also among the longest available remote sensing continuous time series and an irreplaceable tool in climatological and hydrological studies, complementing in situ long-term observations where they are absent or too coarse. The trade-off associated with the high temporal resolution of PMW data is the relative coarse spatial resolution (historically of the order of tens of kilometers). This can represent a limiting factor when studying surface melting as a substantial portion of meltwater production and runoff occurs along the margins of the ice sheet, with some of these areas being relatively narrow (of the order of a few tens of kilometers or smaller, depending on the geographic position and time of the year). The use of a product with a finer spatial resolution would allow a more effective mapping of surface melting and would also permit a better comparison between quantities measured in situ and satellite-derived estimates, reducing uncertainties in the satellite products and allowing for potential improvements to retrieval algorithms. Lastly, finer spatial resolution tools could be helpful, should they be proven effective, in improving mapping of meltwater over ice shelves in Antarctica and furthering our understanding of the processes leading to ice shelf collapse or disintegration (e.g., van den Broeke, 2005; Tedesco, 2009).

In this paper, we report the results of a study in which surface melting over Greenland is estimated by making use of a recently released product developed within the framework of a NASA “Making Earth System Data Records for Use in Research Environments” (MeASUREs) project (, last access: 26 May 2021). The product contains daily maps of PMW Tb generated at an enhanced spatial resolution of a few kilometers (depending on the frequency, as explained below) between 1979 and 2019. The historical gridding techniques for PMW sensors (Armstrong et al., 1994, updated yearly; Knowles et al., 2000, 2006) are based on a “drop in the bucket” approach, in which the gridded value is obtained by averaging the Tb data falling within the area defined by a specific pixel. In the case of the enhanced-spatial-resolution product, the reconstruction algorithm adopted to build the Tb maps makes use of the so-called effective measurement response function (Long and Daum, 1998), determined by the antenna gain pattern, which is unique for each sensor and sensor channel. This pattern is used in conjunction with the scan geometry and the integration period, allowing for “weighting” of measurements within a certain area. The approach used to generate the enhanced resolution product, a radiometer version of the scatterometer image reconstruction algorithm, also addresses another issue in the historical PMW dataset, which is the need to meet the requirements of modern Earth system data records or climate data records, most notably in the areas of inter-sensor calibration and consistent processing methods. More details are reported in Sect. 2.1.

We divide the results of our study into two main parts. In the first part, we report the results of the cross-calibration of different PMW sensors over the Greenland ice sheet to assure a consistent and calibrated Tb time series. Specifically, we use the newly developed spatially enhanced PMW product at the Ka band (37 GHz, horizontal polarization) in view of its sensitivity to the presence of liquid water within the snowpack (Ulaby et al., 1986; Macelloni et al., 2005). We prefer this frequency to the  19 GHz, generally used in the literature as it is less sensitive to liquid water clouds (Fettweis et al., 2011; Mote, 2007), because Tb values at the Ka band are distributed at the highest spatial resolution of 3.125 km (Brodzik et al., 2016). The atmospheric effect on the 37 GHz Tb is higher than that at 19 GHz at low Tb values. However, when considering higher Tb values, the difference in the atmospheric effect between 37 GHz and 19 GHz Tb decreases (Tedesco and Wang, 2006). We then focus on assessing whether the noise introduced by the gridding algorithm might limit the application of the enhanced dataset to mapping surface meltwater. Following this, we focus our attention on testing and assessing existing approaches to deriving melt from PMW data and put forward an update to a recently proposed algorithm in which meltwater is detected when Tb exceeds a threshold computed using the outputs of an electromagnetic model (Tedesco, 2009). We compare results from these algorithms with estimates of surface melting obtained from data collected by automatic weather stations (AWSs) in terms of melt timing and with the outputs of the Modèle Atmosphérique Régional (MAR; Fettweis et al., 2017) regional climate model in terms of melt timing and extent. Lastly, we focus on the analysis of melting patterns and trends over the study period and investigate the information content in the enhanced-resolution dataset using a semi-variogram analysis.

2 Datasets and methods

2.1 Enhanced-resolution passive microwave data

We use Ka band (37 GHz) horizontally polarized Tb data produced within the framework of a NASA MeASUREs project and distributed at the spatial resolution of 3.125 km (Brodzik et al., 2016) over the Northern Hemisphere. Specifically, we use data collected by the Scanning Multichannel Microwave Radiometer (SMMR), SMMR-Nimbus 7; the special sensor microwave/imager (SSM/I), SSM/I-F08, SSM/I-F11 and SSM/I-F13; and the special sensor microwave imager/sounder, SSMI/S-F17 (Table 1), because of its higher orbit stability (, last access: 26 May 2021). Currently, the product time series begins in 1979 and ends in 2019. Data are provided twice a day, as morning and evening passes. Beginning and ending acquisition times for the morning and evening passes are contained within the product's metadata, along with other information. More information can be found at (last access: 26 May 2021).

Table 1Characteristics of the PMW sensors used for this work. LST denotes local solar time.

Download Print Version | Download XLSX

Historical gridding techniques for PMW spaceborne datasets (Armstrong et al., 1994, updated yearly; Knowles et al., 2000, 2006) are relatively simplistic and were produced on grids (Brodzik and Knowles, 2002; Brodzik et al., 2012) that are not easily accommodated in modern software packages. Specifically, the coarse-resolution gridding methodology is based on a simple “drop-in-the-bucket” average, i.e., all of the measurements within a given time falling into a specific pixel are averaged. In the reconstruction algorithm used for the enhanced Tb, the so-called effective measurement response function, determined by the antenna gain pattern and being unique for each sensor and sensor channel, is used in conjunction with the scan geometry and the integration period. The gridding approach uses the Backus–Gilbert technique (Backus and Gilbert, 1967, 1968), a general method for inverting integral equations, which has been applied for solving sampled signal reconstruction problems (Caccin et al., 1992; Stogryn, 1978; Poe, 1990), for spatially interpolating and smoothing data to match the resolution between different channels (Robinson et al., 1992), and for improving the spatial resolution of surface Tb fields (Farrar and Smith, 1992; Long and Daum, 1998). An example of Tb maps at 37 GHz, horizontal polarization, in the case of both the coarse- and enhanced-resolution products over Greenland on 16 July 2001 is reported in Fig. 1. The higher detail captured by the enhanced spatial resolution is clearly visible, especially along the ice sheet edges, where melting generally occurs at the beginning of the season and lasts for the remaining part of the summer. Figure 2 shows an example of time series of both coarse and enhanced PMW Tb (again at 37 GHz, horizontal polarization) for the pixel containing the Swiss Camp station. From the figure, we observe that the two time series are highly consistent with each other, with a mean difference of 0.895 K and a standard deviation of 4.89 K, indicating that the potential noise introduced by the enhancement process is not a major issue. However, differences do exist, as in the case of 3 April 2012 (day of the year 93), when the enhanced product suggests the presence of melting while the coarse product does not. This is likely due to the different spatial resolutions between the two products, as we discuss in the following sections, and shows the added value of using the 37 GHz frequency in detecting small-scale features of the melting process.

Figure 1Maps of PMW Tb at 37 GHz, horizontal polarization, acquired over Greenland on 16 July 2001 with the (a) coarse (25 km) and (b) enhanced (3.125 km) resolution products. Panels (c) and (d) refer to the area highlighted in the square in panels (a) and (b).

To compare our results with precedent PMW surface melting products, we also perform our calculations with the commonly used PMW dataset by Mote (2014). The dataset consists of a daily melt maps product obtained from 37 GHz Tb from SMMR, SSM/I and SSMI/S available on a 25 km grid. This comparison enables us to see the advantages of using the enhanced-resolution product with respect to a coarser-resolution surface melt product. The dataset is generated using a dynamically changing threshold (Mote, 2007) obtained by simulating the Tb change for every grid cell, each year, using a microwave emission model (Mote, 2014). The dataset covers the temporal window from 1979 to 2012 and is available at the National Snow and Ice Data Center (NSIDC) website (, last access: 26 May 2021).

2.2 Greenland air temperature data

In order to assess the results obtained from PMW data, we use in situ data collected by AWSs distributed over the Greenland ice sheet. In the absence of direct observations of melting, we use air temperature (3 m above the surface) to extrapolate instances when liquid water is present, following the procedure adopted by Tedesco (2009) for Antarctica. Specifically, we use data recorded by stations of the Greenland Climate Network (GC-Net; Steffen et al., 1996). The AWSs provide continuous measurements of air temperature, wind speed, wind direction, humidity, pressure and other parameters. We focus on air temperature data collected every hour by the 17 selected stations reported in Table 2. We considered a validation period from 2000 (when all of the considered AWSs were in operation) to 2016 and used daily averaged values. More information about the GC-Net dataset can be found at (last access: 26 May 2021).

Table 2Locations of the automatic weather stations of the Greenland Climate Network (GC-Net) sites used to validate the results in this study.

Download Print Version | Download XLSX

2.3 The MAR model

We assess the enhanced PMW-based surface melt maps with the outputs of the Modèle Atmosphérique Régional (MAR, e.g., Alexander et al., 2014; Fettweis et al., 2013, 2017; Tedesco et al., 2013) regional climate model. MAR is a modular atmospheric model that uses the σ-vertical coordinate to simulate airflow over complex terrain and the soil ice snow vegetation atmosphere transfer scheme (e.g., De Ridder and Gallée, 1998) as the surface model. MAR outputs have been assessed over Greenland (e.g., Fettweis et al., 2005, 2017; Alexander et al., 2014). The snow model in MAR, which is based on the Crocus model of Brun et al. (1992), calculates albedo for snow and ice as a function of snow grain properties, which in turn depend on energy and mass fluxes within the snowpack. Lateral and lower boundary conditions of the atmosphere are prescribed from reanalysis datasets. Sea surface temperature and sea ice cover are prescribed using the same reanalysis data. The atmospheric model within MAR interacts dynamically with the surface model.

In this study, we use the output from MAR version v3.11.2, which is characterized by an enhanced computational efficiency and improved snow model parameters (Fettweis et al., 2017; Delhasse et al., 2020). The model is forced at the boundaries using ERA5 reanalysis (Hersbach et al., 2020), which is the newest generation of global atmospheric reanalysis data that superseded ERA-Interim (Dee et al., 2011), and output is produced at a horizontal spatial resolution of 6 km. In order to compare output from MAR with estimates of melt extent obtained from PMW data, we average the liquid water content (LWC) simulated by MAR along the first 5 cm and 1 m from the surface of the vertical profile of the snowpack, following Fettweis et al. (2007).

2.4 Melt detection algorithms

Generally speaking, melt detection algorithms can be divided into threshold-based and edge-detection algorithms (e.g., Liu et al., 2005; Joshi et al., 2001; Steiner and Tedesco, 2014). Here, we focus on threshold-based algorithms, detecting melting when Tb values (or their combination) exceed a defined threshold, which is computed in different ways depending on the algorithm. For example, Steffen et al. (1993) used the normalized gradient ratio GR =(Tb19H-Tb37H)/(Tb37H+Tb19H) to detect wet pixels with a threshold value computed based on in situ measurements. This method was later updated by Abdalati and Steffen (1995) who introduced the cross-polarized gradient ratio XPGR =(Tb19H-Tb37V)/(Tb37H+Tb19V), where the Ka-band component of the algorithm was switched from horizontally to vertically polarized.

Ashcraft and Long (2006) proposed a threshold based on dry (winter) and wet snow Tb as Tc=αM+(1-α)Twet, where M is the average of winter Tb (January and February), Twet is fixed as 273 K and Tc indicates the threshold value (we keep the same notation in the following). The mixing coefficient α=0.47 was derived considering LWC = 1 % in the first 4.7 cm of the snowpack. Similarly, a method based on a fixed threshold (set to 245 K and derived from the outputs of an electromagnetic model) above which melting is assumed to be occurring was proposed in Tedesco et al. (2007). Several other studies have detected melting when the Tb values exceed the mean winter value plus an additional value ΔTb (MTb approaches) associated with the insurgence of liquid water within the snowpack:

(1) T c = M + Δ T b .

Torinesi et al. (2003) proposed a value of ΔTb=Nσ with M and σ (standard deviation of the time series) varying in space (specific pixel) and time (specific year) but a fixed value of N=3 from the analysis of weather station temperature data. Zwally and Fiegles (1994) used a fixed value of ΔTb=30 K. Tedesco (2009) proposed an alternative approach based on the outputs of the Microwave Emission Model of Layered Snowpack (MEMLS) electromagnetic model (Wiesmann and Mätzler, 1999). In this case, an ensemble of outputs is generated by the MEMLS model by varying the inputs (e.g., correlation length, LWC and snow density). These outputs are then used to build a linear regression model for the ΔTb that is a function of the winter Tb value as follows:

(2) Δ T b = φ M + ω ,

with the values of the coefficients obtained from the linear regression. This is done to account for the increment related to the presence of LWC within the snowpack as a function of the snow properties: a fixed increase would correspond to different LWC values, potentially making the mapping of the wet snow areas inconsistent in terms of LWC values. For example, a snowpack with a small grain size would require a relatively larger amount of LWC compared with a snowpack with a larger grain size for the Tb values to increase by 30 K. From a complementary point of view, an increase of 30 K due to the presence of liquid water in the case of a snowpack with relatively coarse grains would correspond to a lower LWC value than an increase occurring in a snowpack with a smaller grain size. In summary, the adoption of this approach provides consistency in terms of the minimum LWC that is detected by the algorithm. Building on Tedesco (2009), we considered the LWC value of 0.2 %. The coefficients are φ=-0.52 and ω=128 K (R2=0.92). The Tb threshold value computed in this case can, therefore, be written as follows:

(3) T c = M + φ M + ω = ( 1 + φ ) M + ω = γ M + ω ,

where γ and ω assume the values of 0.48 and 128 K, respectively.

Here, we focus five approaches: the MTb approach with ΔTb equal to 30 K and, to test the sensitivity to Zwally and Fiegles (1994), 35 and 40 K (M+30, M+35 and M+40 from here on); the algorithm based on MEMLS for LWC = 0.2 % (referred to simply as MEMLS from here on for brevity); and the 245 K fixed threshold (245 K from here on). We selected MTb and MEMLS due to their higher accuracy in detecting both sporadic and persistent melting with respect to the other approaches presented above, i.e., Torinesi et al. (2003) and Ashcraft and Long (2006), as proven in previous studies (Tedesco, 2009). We also selected the 245 K approach to test a more conservative technique aimed at detecting persistent melting only. In the following sections, we report the results of two algorithms, namely the one using a fixed threshold of 245 K and the one based on MEMLS.

2.5 Inter-sensor calibration

In view of the novelty of this PMW dataset introduced by the enhancement in spatial resolution thanks to the improvement of the gridding technique, we first focus on the cross-calibration of the data acquired by the different sensors. This initial processing step aims to account for biases and differences associated with swath width, view angle, altitude and local time of day as well as the specific intrinsic differences associated with each sensor on the different platforms (Table 1). Several approaches have been proposed in the literature to address this issue for the historical, coarser-spatial-resolution gridded datasets (e.g. Cavalieri et al., 2012). For example, Jezek et al. (1993) compared SMMR and SSM/I over the Antarctic ice sheet for the K and Ka bands ( 19 and  37 GHz) for both horizontal and vertical polarizations. Steffen et al. (1993) proposed an approach that was focused over Greenland for the K band. Abdalati et al. (1995) derived relations between SSM/I observations for the SSM/I-F08 and SSM/I-F11 platforms over Antarctica and Greenland for 19.35, 22.2 and 37 GHz. Dai et al. (2015) intercalibrated SMMR, SSM/I (SSM/I-F08 and SSM/I-F13) and SSMI/S-F17 over snow-covered pixels in China and SMMR, SSM/I and AMSR-E (Advanced Microwave Scanning Radiometer on NASA's Earth Observing System) over the whole Earth surface sampling hot and cold pixels.

Given the novelty of the Tb products used here and the absence of the specific intercalibration of data collected from different platforms for this product, we developed an ad hoc intercalibration for the enhanced PMW dataset. Following Stroeve et al. (1998), we perform the intercalibration using only data collected over the Greenland ice sheet. We perform a linear regression between the data acquired by two sensors over the Greenland ice sheet and calculate the slope (m) and intercept (q) of the linear regression:

(4) y = m x + q .

In Eq. (4), x and y represent the Tb values from coincident data from the two overlapping satellite products. We consider two approaches to compute the m and q values in Eq. (4). In the first method, we compute the weighted average of the daily slope and intercept values from the regression of daily data. Considering n days, for every ith day we first compute mi, qi and the coefficient of determination for the linear regression of Eq. (4) (Ri2); we then average them according to Eqs. (5) and (6):


This choice assigns higher values to the weights obtained from pairs of data with higher correlation. In the second method, we consider all values for all days when data from both platforms are available and then evaluate m and q using a linear regression fitting procedure based on least squares fitting. Using the estimated m and q values, we then correct the values for one of the sensors by applying Eq. (4) to the Tb values of one sensor (x, e.g., SMMR) to obtain new corrected Tb values (y).

We perform an additional comparison using the average difference between the Tb values and evaluating the match between histograms of the overlapping data (Dai et al., 2015) by means of the Nash–Sutcliffe efficiency (NSE) coefficient (Nash and Sutcliffe, 1970), defined as

(7) NSE = 1 - h i ( T b A ) - h i ( T b B ) 2 h i ( T b B ) - h ¯ ( T b B ) 2 .

Here, hi is the absolute frequency of the ith value of Tb of the two sensors (A and B) considered. The NSE is usually applied in calibration and validation procedures to assess the match between measured and modeled quantities, as in Sect. 4.2. After the application of linear relations found using Eqs. (4)–(6), in order to quantitatively assess the impacts of the intercalibration on Tb values, we computed the absolute difference between the values of the histograms of the Tb as follows:

(8) D i = h i ( T b A ) - h i ( T b B ) ,

where Di is the absolute difference between the two histograms A and B for the ith value of Tb. We then summed the differences over the total number of pixels and computed the relative variation as follows:

(9) d = D original - D corrected D original ,

where Doriginal and Dcorrected are the summations of Di before and after the calibration, respectively. The relative variation in d can range from −∞, indicating a worse match of the histograms, to 1, indicating a perfect match of the histograms after the intercalibration.

2.6 Spatial autocorrelation: the variogram analysis

Variogram analysis is generally adopted in geostatistical analyses to evaluate the autocorrelation of spatial data (Delhomme, 1978; Edward et al., 1989) with variograms being characterized by three parameters: the sill, the range and the nugget effect. The sill is the variance value at which the variogram becomes flat. The range is the distance at which the variogram reaches the sill. Beyond this value, the data are no longer autocorrelated. The nugget effect is the variance value at null distance, theoretically zero and resulting from measurement errors or highly localized variability. We computed the empirical variogram as

(10) γ ( δ ) = 1 2 N ( δ ) i , j N ( δ ) ( x i - x j ) 2 ,

where γ is the semi-variance, N(δ) is the number of pair measurements (i,j) spaced by distance δ, and xi and xj are the respective values of the ith and jth measured variable. Generally, the semi-variance γ increases as the distance δ increases according to the principle that close events are more likely to be correlated than distant events. The experimental variogram is the graphical representation of the semi-variance γ as a function of the distance δ. Finally, the experimental variogram is fitted with a function (here, we use a spherical function) to compute the sill, the range and the nugget effect.

3 Results and discussion

3.1 Inter-sensor calibration of enhanced-resolution passive microwave data

First, we show the results obtained for the inter-sensor calibration of the selected satellite constellation. The overlapping periods for the different sensors are as follows: SMMR and SSM/I-F08 overlap between 9 July and 20 August 1987 for a total of 22 d (1 d every 2 d as sensed by the SMMR sensor); SSM/I-F08 and SSM/I-F11 overlap between 3 and 18 December 1991 for a total of 16 d; SSM/I-F11 and SSM/I-F13 overlap between 3 May and 30 September 1995 for a total of 76 d; and SSM/I-F13 and SSM/S-F17 overlap for the period from 1 March to 10 December 2008 for a total of 71 d. In Fig. 3, we show the scatterplots of the data used for the linear regression for Greenland for both evening and morning passes for the SMMR and SSM/I-F08 sensors, reporting m, q and R2 values. We point out that the overlap between SMMR and SSM/I-F08 data occurs in the months of July and August. During these months, the differences between acquisition times (Table 2) might lead to biases and errors associated with snow conditions (e.g., wet vs. dry). Specifically, we expect larger errors at the beginning of the melting season when snow undergoes thawing–refreezing cycles during the day, with frozen snow (low Tb values) potentially being present early in the morning and late at night (SMMR ascending and SSM/I-F08 descending passes) and liquid water (high Tb values) being present during the day. In Table 3, we report the average values of the difference between pairs of Tb data and the NSE coefficient values for the histograms of the same pairs. In Table 4, we report the values for the slope and intercept obtained from the linear regression analysis of enhanced PMW Tb (37 GHz, horizontal polarization) over Greenland for SMMR vs. SSM/I-F08, SSM/I-F08 vs. SSM/I-F11, SSM/I-F11 vs. SSM/I-F13 and SSM/I-F13 vs. SSM/S-F17, along with the R2 values and d values computed according to Eq. (6). In the case of SSM/I and SSMI/S, R2 values are higher, mostly around 0.98. In Fig. 4, we also show examples of histograms for the SMMR and SSM/I-F08 sensors. Large differences are obtained for these sensors for both the evening and morning passes, likely because of the difference in overpass time and the presence or absence of melting in some of the scenes observed by one sensor but not present in the other (Table 3). Conversely, for the SSM/I and SSMI/S sensors, the average difference is close to 0 K (with the exception of the F-08 and F-11 satellites, which show an average difference slightly larger than 1 K, consistent with previous results obtained by Abdalati et al. (1995) in the case of the 25 km resolution data), and the NSE values are extremely close to 1 (Table 3). For of SMMR and SSM/I-F08, the higher average difference (ranging between 3.4 and 4.3 K) and the relatively lower NSE values (ranging between 0.89 and 0.96) show that these sensors have the largest bias. Lastly, we only applied the correction to SMMR, and we did not apply the linear regression to the SSM/I-F08–SSM/I-F11, SSM/I-F11–SSM/I-F13 and SSM/I-F13–SSMI/S-F17 datasets as, in this case, the linear correction worsened the agreement between the two sets of measurements. We applied the correction coefficients obtained with the second method according to the higher relative improvement for the evening pass.

Table 3Average enhanced-resolution Tb differences at 37 GHz, horizontal polarization, for the different PMW sensors and NSE coefficients computed for the histograms of Tb.

Download Print Version | Download XLSX

Table 4Slope (m) and intercept (q) obtained from the linear regression analysis between the selected satellite pairs' enhanced PMW Tb at 37 GHz, horizontal polarization, over Greenland. The subscripts refer to the case where the coefficients are weighted by means of the R2 (case 1, see Eqs. 5 and 6) or cases where they are not weighted by means of the R2 (case 2). We also report the values for the R2 and the d values computed according to Eq. (9).

Download Print Version | Download XLSX

Figure 2(a) Time series of Tb at 37 GHz, horizontal polarization, for the year 2012 for the pixel containing the Swiss Camp site in the case of the coarse (blue) and enhanced (red) products. Threshold values, shown as horizontal lines, are obtained from two approaches considered in this study: 245 K and MEMLS. (b) The difference between the 3.125 km and the 25 km Tb time series for the same pixel (mean of 0.895 K and standard deviation of 4.89 K).


Figure 3Density scatterplots of SMMR and SSM/I-F08 Tb data sensed during the overlap period (9 July–20 August 1987) of the two sensors over the Greenland ice sheet for (a) morning and (b) evening passes. The solid black lines show the linear fitting, and the dashed black lines show the 1:1 line. The color palette indicates the relative frequency.


Figure 4Histograms of Tb before and after the application of the intercalibration relations (for Greenland). Relations are applied for both evening (a) and morning (b) passes, and the histograms of the data and the distance (absolute value of the difference as in Eq. 8) between the histograms for original data are reported. The left column represents the uncorrected data, the central column represents the results applying the correction to SMMR data and the right column represents the results applying the correction to the SSM/I data.


3.2 Assessment of melt detection algorithms

In order to assess the capability of the selected algorithms, we compare the outputs obtained by PMW data with in situ air temperature averaged daily from AWSs as an index of surface melting (Braithwaite and Oelsen, 1989) and with the liquid water content simulated by the MAR regional climate model. We first evaluate performance at the local scale (at the specific locations of the selected AWSs), comparing the number and the concomitance of melting days according to PMW data and the ground-truth reference. According to the results obtained, we then focus on the MEMLS and 245 K algorithms to evaluate the capability of the two approaches with respect to describing the surface melt extent at the ice sheet scale.

3.2.1 Assessment with AWS data

Historically, the presence of liquid water within the snowpack using data from AWS has been estimated when recorded air temperature exceeds a certain threshold during the day. Because melting can also occur due to radiative forcing (i.e., solar radiation) and the air temperature does not necessarily represent the snow surface temperature, we tested three threshold values for air temperature, 0, 1 and 2 C, as in Tedesco (2009). We assessed the performance of the PMW-based algorithm by defining commission and omission errors. Commission errors occur when melting is detected by PMW data but not by AWS data, and omission errors occur when melting is detected by AWS data but not by PMW data. The results of the error analysis are summarized in Table 5 for the different algorithms and for the Mote (2014) dataset for the different threshold values on the AWS air temperature values. In Table 5, commission and omission error values are reported as an average over all stations. Specific results for each AWS location are reported in the Supplement (Tables S2, S3 and S4) along with a map of the AWS network (Fig. S1). As a more general performance indicator, Table 5 also reports the sum of the two commission and omission errors, computed for each AWS case (C+O). Moreover, we show an average value of all of the C+O (for both AWS and MAR assessments, presented in the next subsection) for each PMW algorithm (C+O Mean) as a synthetic index of performance. Our results indicate that the 245 K algorithm shows the lowest commission error (between 0.31 % and 0.63 %) and the highest omission error (between 5.38 % and 9.19 %). This is consistent with this algorithm being the most conservative among those considered (i.e., the algorithm is not sensitive to sporadic melting). In contrast, a higher commission error is achieved for the M+30, M+35 and M+40 thresholds, particularly for the Humboldt and GITS stations (northwestern Greenland), where the commission error is up to 1 order of magnitude larger than the MEMLS and 245 K algorithms for every ground-truth reference (e.g., from 0.70 % for MEMLS and 0.09 % for 245 K to 5.63 % for M+30, for an air temperature of 0 C). Moreover, for the MEMLS algorithm, we note the lowest omission error for the Swiss Camp, JAR-1 and JAR-2 sites (6.9 % for MEMLS and 8.4 % for M+30, 10.0 % for M+35, 12.4 % for M+40 and 17.4 % for 245 K). The coarse-resolution dataset presents a commission error of between 1.02 % and 1.74 % and an omission error of between 4.06 % and 7.12 %, confirming the capability of historical data to detect surface melting over the Greenland ice sheet. However, the enhanced-resolution dataset presents better results in terms of C+O when applying the M+40 and MEMLS algorithms. The sensitivity to the air temperature threshold is low, with commission and omission error decreasing by about 1 % and increasing by 3 %, respectively, when considering threshold values from 0 to 2 C.

Table 5Performance of the PMW melt detection algorithms studied with AWS and MAR data. Five thresholds are used to detect melt at a 3.125 km resolution, and the Mote (2014) melting dataset is used as a 25 km resolution comparison. For each case, three thresholds (0, 1 and 2 C) are applied to the AWS data and two approaches (MAR1 m and MAR5 cm) are applied to the MAR-simulated LWC to detect melt. The performance of the respective PMW melting products is computed in terms of commission and omission errors averaged for all of the AWS sites considered. C+O refers to the total error considering both commission and omission. The average of the C+O of each melting dataset (C+O Mean) is reported as a synthetic index of performance.

Download Print Version | Download XLSX

In order to better understand the sources of the relatively high values of the commission errors at some locations, we show the time series of air temperature and Tb at 37 GHz, horizontal polarization, at three selected stations – (a) Summit, (b) Humboldt and (c) Swiss Camp – for the year 2005 (Fig. 5). The threshold values obtained with the different detection algorithms are also plotted as horizontal lines (black) along with the 0 C air temperature threshold (magenta). We selected these three locations as they are representative of three different environmental and melting conditions. The time series recorded at Summit station (Fig. 5a) shows the sensitivity of Tb to physical temperature and its seasonal variations. In this case, the air temperature remains below 0 C throughout year and the Tb signal does not exceed any threshold value (horizontal lines). This time series is typical of a location where melting is generally absent. The Tb time series collected for the Humboldt location (Fig. 5b) shows a strong and sudden peak starting on 20 July, when the air temperature average is about 0.5 C (detected by 1 and 2 C air temperature thresholds). This event is detected by all algorithms. Nevertheless, M+30 (and similar algorithms) also indicates the potential presence of melting for the period preceding the July melting (between 17 June and 17 July). This melting is not confirmed by the other algorithms or by the AWS analysis, suggesting that the threshold value used for these algorithms might be too low. Lastly, melting clearly occurs for the Swiss Camp site (Fig. 5c), which is characterized by a sharp and substantial increase in Tb beginning around mid-May. For this case, all algorithms detect melting, with MEMLS providing the lowest threshold and the 245 K fixed threshold being the most conservative. The computed rough estimation of the average emissivity for the period from 17 June to 17 July (Tb divided by the recorded air temperature) also suggests that melting is not occurring in the period considered for the Humboldt case, presenting an average emissivity that is even lower than in the Summit case. Figure 6 shows maps of surface melt extent obtained using the different approaches for 13 July 2008. Consistent with the results discussed above, the M+30 and M+35 algorithms suggest melting up to high elevations, within the dry snow zone, where it likely did not occur. The M+40 and MEMLS algorithms show similar results, whereas the 245 K fixed-threshold approach shows the most conservative estimates, as expected. As mentioned, the threshold algorithms for ΔTb (M+30 etc.) rely on a fixed ΔTb value, which could produce errors if a large seasonal range in Tb existed due to temperature variability. In contrast, the MEMLS algorithm is based on the linear regression of the ΔTb as function of different combinations of dry snow conditions (LWC =0, i.e., different winter Tb means). This provides an appropriate threshold value that considers the snow conditions before melting and, at the same time, follows a more consistent approach with respect to the amount of LWC detected in the snowpack.

Figure 5Time series of (blue) enhanced-resolution Tb 37 GHz, horizontal polarization, and (red) air temperature at the (a) Summit, (b) Humboldt and (c) Swiss Camp stations for the year 2005. Threshold values obtained with the different detection algorithms are reported as horizontal black lines (solid line, MT; dashed line, 245 K; and dot-dashed line, MEMLS), and the 0 C threshold is reported as a magenta solid line. The 30 d window between 17 June and 17 July is shown in the shaded orange area and reports the average estimated emissivity (ε) values.


3.2.2 Assessment with MAR outputs

For the comparison between PMW-based and MAR outputs, we averaged the vertical profiles of LWC computed by MAR to the top 5 cm (MAR5 cm) and the top 1 m (MAR1 m) of the snowpack following Fettweis et al. (2007). In order to be consistent with the minimum LWC to which the MEMLS algorithm is sensitive, we set the threshold on the LWC values at which we assume melting is occurring to 0.2 % for both depths. We selected two different depths for our analysis so that we could study two types of melting events: (1) sporadic surface melting, affecting the first few centimeters of the snowpack, and (2) persistent subsurface melting, affecting the snowpack from the surface up to around the first meter. For consistency with the AWS analysis, we report the results averaged over those MAR pixels containing the AWS stations discussed in the section above in Table 5. The comparison between the results obtained from the PMW and modeled LWC indicates that the more conservative approaches (i.e., 245 K) perform better when considering the case of MAR5 cm. In fact, the 245 K threshold shows the lowest overall error for this case (C+O = 6.06 %). The coarse-resolution dataset shows a C+O error equal to 7.39 %, which is slightly lower than for MEMLS (7.56 %). For the top 1 m, all of the algorithms present similar performance on average, with the best performance again obtained for 245 K (C+O = 5.53 %). However, all of the MTb algorithms present the same issue of larger commission errors compared with MEMLS and 245 K (e.g., from 0.99 % for MEMLS and 0.26 % for 245 K to 4.62 % for M+30) in northeastern Greenland (e.g., Humboldt and GITS stations, see Tables S2 and S3 in the Supplement). This confirms the results that we obtained from the comparison with AWS data, pointing out the overestimation of melting in some dry areas by MTb. For both the MAR1 m and MAR5 cm cases, for all of the algorithms considered, we find a high commission error for the JAR-1, JAR-2 and Swiss Camp sites (between 10 % and 22 %).

In order to better understand the origins of these errors, we show further insights into the differences between the PMW Tb and MAR outputs in Fig. 7. Figure 7a and b show the time series of LWC averaged over the first 5 cm (MAR5 cm) and 1 m (MAR1 m) obtained from MAR at the Swiss Camp site, respectively. In Fig. 7c, we report the Tb time series and the daily average air temperature (threshold values reported as horizontal lines). We first note an early melt event (labeled LWC = 0.046 % in Fig. 7b for the 108th day of the year) detected by the PMW MEMLS algorithm and at the AWS station but apparently undetected in MAR1 m. A closer look at the time series shows that MAR1 m does in fact estimate a LWC of 0.046 % on this day, whereas MAR5 cm estimates a LWC of 0.6 %. This suggests that in some cases (before the main melt season) the MEMLS algorithm is actually sensitive to the LWC in the first 5 cm of the snowpack, as a consequence of the approximation of the electromagnetic outputs imposed by the linear fitting. We also note a melt event (labeled as LWC = 0.5 % in Fig. 7b) at the end of the melting season detected by both AWS data (air temperature larger than 1 C) and MAR (both in the first 5 cm and 1 m of snow) but not by any PMW algorithm. The Tb time series reveals a small peak, but the signal is not strong enough to exceed any threshold. This corresponds to a rainfall event (simulated by MAR), suggesting that the sensitivity of the 37 GHz channel to liquid clouds could mask some melt events. Moreover, the Tb appears to be slightly lower than January–February average at the end of the melting season, possibly due to an increase in grain size after refreezing, leading to a lower emissivity.

Figure 6Melting maps obtained using the (a) M+30, (b) M+35, (c) M+40, (d) 245 K and (e) MEMLS algorithms over the Greenland ice sheet on 13 July 2008. An example of an area presenting the false detection problem is shown in the red circle.

Figure 7LWC from MAR averaged in the first 5 cm (a) and the first 1 m (b) of the snowpack. (c) Time series of the 37 GHz horizontally polarized Tb (3.125 km, blue), air temperature from AWS (black) and 245 K (dashed magenta line), MT (solid magenta lines) and MEMLS (dot-dashed magenta line) thresholds for the Swiss Camp site in the year 2001. Melting days according to MAR are marked as vertical light blue lines in panels (a) and (b).


The results discussed above (and the results from the comparison with AWS data) suggest that 245 K is the most conservative among the approaches tested in this study, providing the lowest (highest) commission (omission) errors, although being unable to detect sporadic melt events. On the contrary, the MEMLS and MTb algorithms can detect sporadic melt events and present lower omission error compared with 245 K. However, the MTb algorithms overestimate melting in some dry areas (the northwest of the ice sheet), suggesting melting when it is not actually occurring. On the contrary, the MEMLS algorithm is not affected by the large commission error in dry areas, presents the lowest omission error in the Swiss Camp area (along with M+30) and is still sensitive to low LWC levels. Considering the average error (C+O Mean in Table 5), the MEMLS algorithm shows the best performance (6.66 %). In view of the analysis presented and the different sensitivity to surface and subsurface melting, in the following we focus on the 245 K and MEMLS algorithms to study the extent of persistent and sporadic surface melting, respectively.

As a further analysis, we compared the PMW-retrieved melt extent with that estimated from MAR outputs. In Fig. 8 we show the time series of the melt extent integrated over the whole ice sheet for 2 selected years – (a) 1983 and (b) 2005, selected randomly to present an example of SMMR and SSM/I cases – estimated according to MAR5 cm and MAR1 m as well as the time series of the melt extent from the PMW data. The analogous figure for the coarse-resolution dataset is reported in the Supplement (Fig. S2). For each year, we compute the daily melt extent for the period from 1 May to 15 September and use the Nash–Sutcliffe efficiency (NSE) coefficient (Nash and Sutcliffe, 1970), described in Sect. 3, for a quantitative analysis. Here, we remind the reader that the NSE coefficient can assume values in the (-,1] interval. A perfect match between the two time series is achieved when the NSE coefficient is 1. NSE coefficient values in the [0,1] interval indicate that the modeled variable is a better predictor of the measurements than the mean. If the NSE coefficient is a negative number, the mean of the measured data describes the time series better than the modeled predictor. Here, we chose NSE = 0.4 as the efficiency threshold, considering that we compute melt extent at a daily timescale and from datasets at two different resolutions (i.e., resulting in an intrinsic bias related to the different pixel size). We compared the time series of the melt extent obtained using the 245 K algorithm with MAR1 m (245 K vs. MAR1 m) and the MEMLS melt extent with MAR5 cm (MEMLS vs. MAR5 cm) due to the expected differences in the sensitivity to detect persistent and sporadic melting between 245 K and MEMLS, respectively. We compare the melt extent obtained from the coarse-resolution dataset with MAR5 cm, according to the similarity with MEMLS in terms commission and omission error. We report NSE coefficients computed for the 41-year (34-year for the coarse-resolution dataset) period in Table 6. At first, we notice that the comparison between 245 K and MAR1 m produces large negative NSE coefficients for the 1979–1992 period, indicating an unsatisfactory match between PMW- and MAR-derived melt extents. The comparison between MEMLS and MAR5 cm also presents negative NSE coefficient values (unsatisfactory results). Similarly, the coarse-resolution dataset shows negative NSE coefficient values, larger than MEMLS but smaller than 245 K. Between 1987 and 1992, we found larger but still negative NSE coefficients presenting smaller absolute values. Between 1993 and 2019, we found negative NSE coefficient values for 245 K and positive values for MEMLS for every year, indicating satisfactory results only for the latter algorithm. However, the coarse-resolution dataset only presents positive (but not satisfactory) NSE coefficient values in 2003 and 2012 (0.111 and 0.208, respectively). The time series in Fig. 8a reveals a strong underestimation of the 245 K-derived melt extent relative to MAR1 m (the cause of the low NSE value, NSE =-151.596) and shows a slightly better match for MEMLS (NSE =-0.540). This result suggests that, for SMMR data, Tb values cannot always reach the 245 K threshold, even if the snowpack is saturated with liquid water and surface melting is developed, possibly due to a persistent bias after the intercalibration of the dataset. As a consequence, the 245 K threshold might be too high in the first part of the dataset, resulting in an underestimation of the melt extent. On the contrary, the MEMLS threshold, which is generally lower, can better capture the spatiotemporal evolution of surface melting, even if the melt extent is still underestimated. A possible consequence of the melt extent being underestimated in the first part of the time series is a slightly overestimated long-term trend. To address this possible implication, we compute long-term trends considering both the 1979–2019 and 1987–2019 reference periods in the next section. In Fig. 8b, the time series obtained with 245 K appears to better follow the temporal variability of the melt extent from MAR during the melting season, although it still presents a strong underestimation (NSE =-5.250). On the other hand, the MEMLS-derived time series better matches the MAR-derived time series, showing a largely satisfactory NSE coefficient (0.782). In these 2 years the NSE coefficient computed for the coarse-resolution case is negative. The magnitude of the errors is lower than for the 245 K algorithm, indicating a weaker underestimation of the melt extent. This can be a consequence of the coarser resolution, which causes issues with respect to capturing melting areas at the edges of the Greenland ice sheet.

Table 6Nash–Sutcliffe efficiency (NSE) coefficients computed for the comparison of the retrieved melt extent using the 245 K and MEMLS algorithms applied to the enhanced-resolution PMW Tb and MAR liquid water content outputs averaged in the first 1 m and the first 5 cm of the snowpack. The Nash–Sutcliffe efficiency coefficients for the comparison of the coarse-resolution dataset (Mote, 2014) are computed considering MAR5 cm.

Download Print Version | Download XLSX

Figure 8Melt extent estimation from PMW 37 GHz horizontally polarized Tb (red) and the MAR (blue) regional climate model. Time series were obtained using the 245 K algorithm and the LWC average in the first 1 m of the snowpack (left), and the MEMLS algorithm and the LWC average in the first 5 cm of the snowpack (right), for the years (a) 1983 and (b) 2005.


In summary, the 245 K threshold (even if presenting acceptable results in terms of commission and omission error considering both the AWS and MAR comparison) is too high to fully capture the melt extent everywhere over the ice sheet. Contrarily, we found that the MEMLS algorithm is suitable for capturing the evolution of melting over the Greenland ice sheet. The comparison with the 25 km historical surface melting dataset shows the same underestimation issue for 245 K. Although the 245 K case showed lower errors, we did not find acceptable NSE coefficient values.

3.3 Surface melting trends

Here, we report results concerning the trends in the melt duration, length of the melting season and melt extent. We define melt duration (MD) as the total number of days on which melting is detected. We compute trends in the MD over the whole ice sheet (mean melt duration, MMD, averaged over the total ice sheet area) and at a pixel-by-pixel scale. We also study the maximum melting surface (MMS, maximum extent of melting area, i.e., the sum of the pixel areas in which melting has been detected at least once over a period, expressed as a fraction of the total ice sheet area) and the cumulative melting surface or melting index (MI, the sum of the melting pixel days multiplied by the area subjected to melting, i.e., the integral of the MD time series; Tedesco et al., 2007). Lastly, we define melt onset date (MOD) and melt end date (MED) as the first day when melting occurs for 2 d in a row and when melting does not occur for at least 2 d in a row. We report the comparison with the trends computed according to the coarse-resolution dataset with reference to the time period of data availability (1979–2012). Figures related to this analysis can be found in the Supplement. In Fig. 9, we report the time series of the annual MMD, MMS and MI values for the 1979–2019 and 1988–2019 reference periods. We decided to look at two different reference periods in view of the fact that SMMR data are collected every other day and that the SMMR and SSM/I sensors are fundamentally different from one another (although this is not true for the remaining SSM/I sensors). We show the results obtained by applying the MEMLS algorithm (the one that presented the best performance in all of the cases considered) and the 245 K threshold (because it presents good performance in the omission and commission error analysis, even if it has the limitation of strongly underestimating the melt extent from the comparison with MAR outputs). For the MMD (Fig. 9a), we obtain a positive statistically significant (p value <0.05) trend from both the 245 K and MEMLS algorithms (except for the 1988–2019 period for 245 K): 0.249 d yr−1 (0.108 d yr−1) for the 245 K algorithm for the 1980–2019 (1988–2019) period and 0.451 d yr−1 (0.291 d yr−1) for the MEMLS algorithm. The trends computed using the coarse-resolution dataset results (Fig. S3a in the Supplement) equal 0.587 d yr−1 for the 1979–2012 period and 0.595 d yr−1 for the 1988–2012 period, which are smaller than MEMLS (0.704 and 0.671 d yr−1, respectively) but larger than 245 K (0.457 and 0.418 d yr−1, respectively). For the MMS (Fig. 9b), both the 245 K and MEMLS algorithms also indicate statistically significant positive trends (p value < 0.05 for every case and p value < 0.1 for MEMLS for the 1988–2019 period). The computed trends suggest that the MMS increased by 0.69 % yr−1 in the case of MEMLS and by 0.94 % yr−1 in the case of the 245 K algorithm for the 1979–2019 period (percentage with respect to the whole ice sheet surface area). For the 1988–2019 period, we also found that the trends are statistically significant but smaller in value (0.36 % yr−1 for MEMLS and 0.47 % yr−1 for 245 K). The obtained trends computed using the 25 km dataset (Fig. S3c in the Supplement) are equal to 1.31 % yr−1 for the 1979–2012 period and 0.91 % yr−1 for the 1988–2012 period, which are larger than the trends computed for MEMLS (1.03 % yr−1 and 0.77 % yr−1) but smaller than those computed for 245 K (1.50 % yr−1 and 1.23 % yr−1). In the case of MI (Fig. 9c), we also found positive statistically significant trends of 9.166 × 105 km2 d yr−1 (MEMLS) and 5.862 × 105 km2 d yr−1 (245 K) for the complete time series. When considering the reduced reference period, we found a 95 % statistically significant trend of 5.726 × 105 km2 d yr−1 only in the case of MEMLS. The trends computed for the 25 km resolution dataset results (Fig. S3b in the Supplement) equal 0.999 × 106 km2 d yr−1 (1979–2012) and 1.019 × 106 km2 d yr−1 (1988–2012), which are smaller than those for MEMLS (1.428 × 106 and 1.325 × 106 km2 d yr−1, respectively) and 245 K for the 1979–2012 period (0.999 × 106 km2 d yr−1) but larger than 245 K for the 1988–2012 period (1.019 × 106 km2 d yr−1). Lastly, in Fig. 9d, we report the MOD and MED averaged spatially over the pixels with 95 % significant trends. We found that the average MOD (crosses in Fig. 9d) presents similar trends for 245 K and MEMLS considering both the entire and shortened time series: 0.546 and 0.273 d yr−1 for 245 K, respectively, and 0.404 and 0.254 d yr−1 for MEMLS, respectively. For the 1979–2012 (1988–2012) reference period, the trend in the MOD computed from the 25 km resolution dataset (Fig. S4a in the Supplement) was equal to 0.585 d yr−1 (0.562 d yr−1), which is equal to and larger in absolute value than MEMLS between 1979 and 2012 (0.585 d yr−1) and MEMLS between 1988 and 2012 (0.494 d yr−1), respectively, but smaller than 245 K between 1979 and 2012 (0.801 d yr−1) and 245 K between 1988 and 2012 (0.568 d yr−1). On the contrary, we found larger differences when considering the reduced time series for the average MED (in red, Fig. 9), with results equal to 0.687 d yr−1 for 245 K (1979–2019), 0.708 d yr−1 for MEMLS (1979–2019) and 0.396 d yr−1 for MEMLS (1988–2019). The 245 K algorithm does not present a statistically significant trend over the 1988–2019 period. This difference suggests that the 245 K algorithm may have stronger limitations with respect to capturing the last portion of the melting season for SMMR data, thereby confirming the problems observed with melt detection using this source of data vs. MAR1 m simulations. The trends computed from the coarse-resolution dataset (Fig. S4b in the Supplement) are equal to 0.850 for the 1979–2012 period and 0.716 for the 1988–2012 period. For the 1979–2012 (1988–2012) period, we found a delay of 0.937 d yr−1 (0.621 d yr−1) for MEMLS and 1.046 d yr−1 (0.521 d yr−1) for 245 K.

Figure 9Time series of annual (a) mean melt duration (MMD), (b) maximum melting surface fraction (MMS, expressed as fraction of the surface area of the ice sheet), (c) melt index (MI), and (d) melt onset date (MOD) and melt end date (MED). Regression lines were computed for the 1979–2019 (solid line) and 1988–2019 (dot-dashed line) periods. The MMD is averaged over all of the Greenland ice sheet pixels. Red (blue) lines refer to the 245 K (MEMLS) algorithm; in panel (d), squares (crosses) refer to MED (MOD).


In Fig. 10 we show the trends in the MD, MOD and MED on a pixel-by-pixel basis for the complete time series (1979–2019). We found that the trend in the MD exhibits the highest statistical significance (in terms of the number of statistically significant pixels), being the most stable and reliable trend among the pixel-by-pixel parameters analyzed. We found mostly positive trends in the MD in all pixels (Fig. 10a, b), with higher values moving towards the coastline, maxima in the ablation zone of the Jakobshavn Glacier (2.40 d yr−1 for 245 K and 2.66 d yr−1 for MEMLS) and minima in high-altitude areas. We averaged the statistically significant trends, finding an average of 0.468 d yr−1 for the 245 K algorithm and an average of 0.697 d yr−1 for MEMLS. For MOD and MED, we found a lower number of statistically significant pixels. The statistically significant pixels exhibit a negative trend for MOD (Fig. 10c, d) and a positive trend for MED (Fig. 10e, f), with the melting season starting on average 0.694 d yr−1 earlier and ending 0.680 d yr−1 later according to the 245 K algorithm (0.360 d yr−1 earlier and 0.909 d yr−1 later for MEMLS). We point out that the average of the statistically significant trends is generally higher than the trends computed at the ice sheet scale as we computed the average over the statistically significant pixels only.

Figure 10Maps of 95 % significant trends (1979–2019) obtained with the 245 K (a, c, e) and MEMLS (b, d, f) algorithms for melt duration (MD; panels a and b), melt onset date (MOD; panels c and d) and melt end date (MED; panels e and f). MOD and MED are defined as the first and last 2 melting days in a row.

3.4 Spatial information content

In order to investigate the spatial information content of the enhanced-resolution PMW data with respect to the coarser-resolution dataset, we also performed a variogram-based analysis of the MD estimated from the two products when using either the 245 K or the MEMLS algorithms. We point out that knowledge of scales is imperative for improving our understanding of the observed changes because processes and related relationships change with scale. Moreover, quantifying the variability of processes across scales is a critical step, ultimately leading to proper observation and modeling scale resolution. In this regard, the relationship between processes, observation and modeling scales controls the ability of a tool to detect and describe the constituent processes. Here, we show our preliminary results of a variogram-based analysis applied to the MD estimated from the MEMLS and 245 K algorithms for the months of May through August of 2012 when using either the enhanced or the coarse-resolution products. We also performed the same analysis applied to the MD estimated using LWC modeled with MAR, according to the same rationale described in the previous sections. Here, we compute the MD for each month of the melting season at the pixel scale as the number of days of the month (May, June, July or August) detected as melting for the specific pixel. The results of our analysis are summarized in Fig. 11, where we show the empirical (blue crosses) and modeled (red line) semi-variograms for the Greenland MD computed by applying the MEMLS and 245 K algorithms to both 25 and 3.125 km resolution data for the months of May through August of 2012, and in Table 7, where we report the parameters of the spherical fitting of the empirical semi-variogram for the MD obtained according to the MAR1 m and MAR5 cm approaches (an analogous representation of Fig. 11 is reported in the Supplement Fig. S5). At first, we note that the R2 values of the fitting for the modeled variograms are consistently higher for the enhanced-resolution data, suggesting that enhanced-resolution data might be more suitable for a variogram-based analysis. For the coarse-resolution data, we only found R2 values of comparable magnitude with the enhanced-resolution case in May. When computing the spherical fitting of the empirical variograms of the MD from MAR, we found considerably similar R2 values (between 0.118 and 0.484) for each case. The values of the range for the 3.125 and 25 km products are similar in May for the 245 K algorithm (of the order of  200 km), but they appear to be different for the MEMLS algorithm – the enhanced product shows a lower value of  170 km in comparison with  270 km for the coarse product. This could be due to the fact that the MEMLS algorithm is more sensitive to sporadic melting, and, when applied to the enhanced Tb dataset, it allows for the detection of melting driven by processes whose scale cannot be captured by the coarser nature of the historical dataset. For MAR, the value of the range is lower in the case of MAR5 cm (187.70 km) than for MAR1 m (199.17 km), again suggesting the affinity of MEMLS algorithm with melting strictly confined to the very first layer of the snowpack. As the melting season progresses, the variograms of the coarse dataset show either similar values for the range or a poor fit of the experimental variogram. In contrast, for the enhanced product, the values of the range tend to decrease up to July and increase again in August. We found the same temporal variability for MAR1 m, whereas we found that the range increases until July and decreases in August for MAR5 cm (see Fig. S5 in the Supplement). Moreover, a proper fitting of the experimental variograms is achieved for all cases for the enhanced-resolution PMW and the MAR-derived MD. This suggests that the 25 km spatial resolution might be too coarse to capture the spatial autocorrelation of melting processes. In terms of the nugget effect, we found larger values from the MAR outputs than for the PMW data. The decrease in the range for the enhanced product may be a consequence of the local processes that drive melting as the melting season progresses (e.g., impact of bare ice exposure, cryoconite holes and new snowfall) and of a more developed network of surface meltwater, the presence of supraglacial lakes and, in general, the fact that the processes driving surface meltwater distribution (e.g., albedo and temperature) promote a stronger spatial dependency of meltwater production at smaller spatial scales. This is even more important when considering that the width of regions such as the bare ice area (where substantial melting occurs) is of the same order of magnitude as the resolution of the coarse PMW dataset. In August, the start of freezing of the surface runoff system and the covering of bare ice, cryoconite holes, and the draining of the supraglacial lakes and rivers might justify the increase in the range values computed for this month. Therefore, our preliminary results point to an increased information content of the enhanced-spatial-resolution PMW product with respect to the historical coarse-resolution product, offering the opportunity to better capture the spatial details of how surface melting has evolved over the Greenland ice sheet over the past  40 years. Further analysis will help to shed light on the processes responsible for the recent acceleration of surface melting.

Table 7Parameters of the spherical function fitted to the empirical semi-variogram for the maps of melt duration (MD) obtained by cumulating the LWC simulated by MAR over the first 1 m and 5 cm of the snowpack.

Download Print Version | Download XLSX

Figure 11Empirical (blue crosses) and modeled (red line) semi-variograms for the Greenland melt duration (MD) computed by applying the (a, b) MEMLS and (c, d) 245 K algorithms to both the (a, c) 25 km and (b, d) 3.125 km resolution data for each month of the melting season (May, June, July and August). The range (r), sill (s), nugget (n) and R2 values are reported.


4 Conclusions and future work

We applied threshold-based melt detection algorithms to the 3.125 km resolution 37 GHz horizontally polarized PMW Tb to assess the skill of the PMW enhanced-resolution data to detect surface melting for the 1979–2019 period over the Greenland ice sheet. As the product is composed of data acquired by different sensors on board different platforms, we first developed a cross-calibration among all of the sensors. We then compared surface melting detected from PMW enhanced-resolution data with that estimated from AWS air temperature data and the outputs of the MAR regional climate model. We found that the algorithm making use of a fixed threshold value for Tb values (245 K) and the algorithm based on the outputs of an electromagnetic model were the most suitable for detecting persistent (245 K) and sporadic (MEMLS) melting. Overall, we found that the MEMLS algorithm showed the best performance (lowest commission and omission errors). We compared surface melting detected from PMW enhanced-resolution data with that estimated by the MAR model when considering the two cases of integrating LWC over the top 5 cm and 1 m, respectively. We selected these two depths to study conditions in which melting occurs sporadically (5 cm) or persistently (1 m). We obtained a good match (i.e., NSE > 0.4 or, at least, positive) in most of the years from 1992 to 2019 when comparing the MEMLS-derived melt extent with MAR LWC in the first 5 cm of the snowpack. In contrast, we found a bad match between the two in the 1979–1992 period, possibly due to differences in sensor characteristics. For the melt extent retrieved by 245 K, we found a strong underestimation of the melt extent (largely negative NSE coefficient values) from 1979 to 1987 that was likely due to the lower values of “wet” Tb for SMMR data; this slightly improved from 1993 to 2019 but was still negative. Accordingly, the results obtained by applying the MEMLS approach are more reliable than those for the 245 K algorithm when considering the 1979–2019 period. In the comparison with the PMW coarse-resolution dataset (25 km), we found that the melt extent time series derived from the PMW enhanced-resolution data using MEMLS showed better agreement with the MAR simulations than those obtained using the 25 km resolution data.

After assessing the outputs of the PMW-based algorithms, we studied the melt onset date, melt end date, mean melt duration and maximum melting surface for the 1979–2019 period. According to the MEMLS algorithm, we found that the melting season began 0.404 (0.254) d earlier every year between 1979 and 2019 (1988–2019) and ended 0.708 (0.396) d later every year between 1979 and 2019 (1988–2019). These values are averaged over the whole ice sheet, and the trends are statistically significant at a 95 % level (p value < 0.05). The mean melt duration increased every year by 0.451 d yr−1 (0.291 d yr−1) during the 1979–2019 (1988–2019) period. We found differences in trends computed using the PMW coarse-resolution data with respect to the PMW enhanced-resolution data for the 1979–2012 and 1988–2012 reference periods, possibly because of the different rationale behind the melt detection algorithms and the higher level of detail in the PMW enhanced-resolution dataset. When we performed a spatial analysis of the trends for the melt onset dates and duration, we found that the areas where the number of melting days has been increasing are mostly located in western Greenland. The maximum melting surface also presents positive trends, with an increment of 0.69 % (0.36 %) every year with respect to the Greenland ice sheet surface since 1979 (1988).

Finally, we explored the information content of the PMW enhanced-resolution dataset with respect to the one at 25 km and the MAR outputs using a semi-variogram approach. The results obtained showed a better fitting of the modeled spherical function to the empirical semi-variogram for the enhanced-resolution data and MAR maps of melt duration. Our analysis suggests that the enhanced-resolution product is sensitive to local-scale processes, with higher sensitivity in the case of the MEMLS algorithm. This offers the opportunity to improve our understanding of the spatial scale of the processes driving melting and potentially paves the way for the use of this dataset in statistically downscaling model outputs. In this regard, as a future work, we plan to extend the analysis of spatial scales to the atmospheric drivers of surface melting, such as incoming solar radiation, surface temperature and longwave radiation, and to complement this analysis with our previous work, which is focused on understanding the changes in atmospheric patterns that have been promoting enhanced melting in Greenland over the recent decades (Tedesco and Fettweis, 2020). As we have assessed the capacity of this dataset and method to observe temporal trends, further development could include a combination of the PMW enhanced-resolution dataset with higher-resolution satellite data (optical sensors or lower frequencies) in order to investigate the evolution of the surface meltwater networks and the application of similar tools to other regions, such as the Canadian Arctic Archipelago, the Himalayan Plateau and the Antarctic Peninsula, where the enhancement of the spatial resolution could be fully exploited.

Code availability

The melting maps were obtained by implementing the algorithms described in this paper and are available from Paolo Colosio upon request. The MAR code is available at (last access: 10 June 2021, MAR model, 2021).

Data availability

Enhanced-resolution passive microwave data are available on the NSIDC portal at (Brodzik et al., 2016). Air temperature data from the Greenland Climate Network (GC-Net) can be found at (last access: 9 June 2021, Steffen Research Group, 2021). MAR outputs used in this study are available at (last access: 10 June 2021).


The supplement related to this article is available online at:

Author contributions

PC, MT, RR and XF conceived the work. PC and MT selected the remote sensing methodology. XF ran the MAR simulations and provided the outputs. PC performed the data analysis. All the authors contributed to paper preparation and refinements.

Competing interests

The authors declare that they have no conflict of interest.


Funding for this paper was provided by the US National Science Foundation (grant nos. ANS 1713072, PLR- 1603331, 1604058 and OPP19-01603), the National Aeronautics and Space Administration (grant nos. 80NSSC18K0814, 80NSSC17K0351 and NNX17AH04G), the Heising-Simons Foundation (grant no. HS 2019–1160) and the Italian Agency for Development Cooperation (AICS; grant no. D24I20000430005). The first author is also grateful for the support provided by the University of Brescia and the Honors Center of Italian Universities (H2CU).

Financial support

This research has been supported by the US National Science Foundation (grant nos. ANS 1713072, PLR-1603331, 1604058 and OPP19-01603); the National Aeronautics and Space Administration, Goddard Space Flight Center (grant nos. 80NSSC18K0814, 80NSSC17K0351 and NNX17AH04G); the Heising-Simons Foundation (grant no. HS 2019–1160); the Italian Agency for Development Cooperation (grant no. D24I20000430005); and the Cariplo Foundation partly funded this research through the FLORIMAP project (grant no. 2017/0708).

Review statement

This paper was edited by Kenichi Matsuoka and reviewed by Masashi Niwano and one anonymous referee.


Abdalati, W. and Steffen, K.: Passive microwave-derived snow melt regions on the Greenland ice sheet, Geophys. Res. Lett., 22, 787–790,, 1995. 

Abdalati, W., Steffen, K., Otto, C., and Jezek, K. C.: Comparison of brightness temperatures from SSMI instruments on the DMSP F8 and FII satellites for Antarctica and the Greenland ice sheet, Int. J. Remote Sens., 16, 1223–1229,, 1995. 

Alexander, P. M., Tedesco, M., Fettweis, X., van de Wal, R. S. W., Smeets, C. J. P. P., and van den Broeke, M. R.: Assessing spatio-temporal variability and trends in modelled and measured Greenland Ice Sheet albedo (2000–2013), The Cryosphere, 8, 2293–2312,, 2014. 

Alexander, P. M., Tedesco, M., Schlegel, N.-J., Luthcke, S. B., Fettweis, X., and Larour, E.: Greenland Ice Sheet seasonal and spatial mass variability from model simulations and GRACE (2003–2012), The Cryosphere, 10, 1259–1277,, 2016. 

Armstrong, R., Knowles, K., Brodzik, M., and Hardman, M. A.: DMSP SSM/I-SSMIS Pathfinder Daily EASE-Grid Brightness Temperatures, Version 2, NASA DAAC at the National Snow and Ice Data Center, Boulder, Colorado, USA, available at: (last access: 26 May 2021), 1994. 

Aschwanden, A., Fahnestock, M. A., Truffer, M., Brinkerhoff, D. J., Hock, R., Khroulev, C., Mottram, R., and Khan, S. A.: Contribution of the Greenland Ice Sheet to sea level over the next millennium, Science Advances, 5, eaav9396,, 2019. 

Ashcraft, I. S. and Long, D. G.: Comparison of methods for melt detection over Greenland using active and passive microwave measurements, Int. J. Remote Sens., 27, 2469–2488,, 2006. 

Backus, G. E. and Gilbert, J. F.: Numerical applications of a formalism for geophysical inverse problems, Geophys. J. Roy. Astr. S., 13, 247–276,, 1967. 

Backus, G. E. and Gilbert, J. F.: The Resolving Power of Gross Earth Data, Geophys. J. Int., 16, 169–205,, 1968. 

Braithwaite, R. J. and Oelsen, O. B.: Calculation of glacier ablation from air temperature, West Greenland, in: Glacier Fluctuations in Cimatic Change, edited by: Oerlemens, J., Kluwer Academic Publishers, Dordrecht, The Netherlands, 219–233, 1989. 

Brodzik, M. J. and Knowles, K. W.: EASE-Grid: a versatile set of equal-area projections and grids, in: Discrete Global Grids, edited by: Goodchild, M., National Center for Geographic Information & Analysis, Santa Barbara, California, USA, 2002. 

Brodzik, M. J., Billingsley, B., Haran, T., Raup, B., and Savoie, M. H.: EASE-Grid 2.0: Incremental but Significant Improvements for Earth-Gridded Data Sets, ISPRS Int. J. Geo-Inf., 1, 32–45,, 2012. 

Brodzik, M. J., Long, D. G., Hardman, M. A., Paget, A., and Armstrong, R. L.: MEaSUREs Calibrated Enhanced-Resolution Passive Microwave Daily EASE-Grid 2.0 Brightness Temperature ESDR, Version 1, National Snow and Ice Data Center, Boulder, Colorado, USA [data set],, 2016. 

Brun, E., David, P., Sudul, M., and Brunot, G.: A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting, J. Glaciol., 38, 13–22, 1992. 

Caccin, B., Roberti, C., Russo, P., and Smaldone, L. A.: The Backus-Gilbert inversion method and the processing of sampled data, IEEE T. Signal Proces., 40, 2823–2825,, 1992. 

Cagnati, A., Crepaz, A., Macelloni, G., Pampaloni, P., Ranzi, R., Tedesco, M., and Valt, M.: Study of the snow melt–freeze cycle using multi-sensor data and snow modelling, J. Glaciol., 50, 419–426,, 2004. 

Cavalieri, D. J., Parkinson, C. L., DiGirolamo, N., and Ivanoff, A.: Intersensor Calibration Between F13 SSMI and F17 SSMIS for Global Sea Ice Data Records, IEEE Geosci. Remote S., 9, 233–236,, 2012. 

Dai, L., Che, T., and Ding, Y.: Inter-Calibrating SMMR, SSM/I and SSMI/S Data to Improve the Consistency of Snow-Depth Products in China, Remote Sens.-Basel, 7, 7212–7230,, 2015. 

De Ridder, K. and Gallée, H.: Land Surface-Induced Regional Climate Change in Southern Israel, J. Appl. Meteorol., 37, 1470–1485, 1998. 

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597,, 2011. 

Delhasse, A., Kittel, C., Amory, C., Hofer, S., van As, D., S. Fausto, R., and Fettweis, X.: Brief communication: Evaluation of the near-surface climate in ERA5 over the Greenland Ice Sheet, The Cryosphere, 14, 957–965,, 2020. 

Delhomme, J.: Kriging in the hydrosciences, Adv. Water Resour., 1, 251–266,, 1978. 

Edward, H., Isaaks, E. H., and Srivastava, R. M.: An introduction to applied geostatistics, Oxford University Press, Oxford, UK, 561 pp., 1989. 

Farrar, M. R. and Smith, E. A.: Spatial resolution enhancement of terrestrial features using deconvolved SSM/I microwave brightness temperatures, IEEE T. Geosci. Remote, 30, 349–355,, 1992. 

Fettweis, X., Gallée, H., Lefebre, F., and van Ypersele, J.-P.: Greenland surface mass balance simulated by a regional climate model and comparison with satellite-derived data in 1990–1991, Clim. Dynam., 24, 623–640,, 2005. 

Fettweis, X., van Ypersele, J.-P., Gallee, H., Lefebre, F., and Lefebvre, W.: The 1979–2005 Greenland ice sheet melt extent from passive microwave data using an improved version of the melt retrieval XPGR algorithm, Geophys. Res. Lett., 34, L05502,, 2007. 

Fettweis, X., Tedesco, M., van den Broeke, M., and Ettema, J.: Melting trends over the Greenland ice sheet (1958–2009) from spaceborne microwave data and regional climate models, The Cryosphere, 5, 359–375,, 2011. 

Fettweis, X., Franco, B., Tedesco, M., van Angelen, J. H., Lenaerts, J. T. M., van den Broeke, M. R., and Gallée, H.: Estimating the Greenland ice sheet surface mass balance contribution to future sea level rise using the regional atmospheric climate model MAR, The Cryosphere, 7, 469–489,, 2013. 

Fettweis, X., Box, J. E., Agosta, C., Amory, C., Kittel, C., Lang, C., van As, D., Machguth, H., and Gallée, H.: Reconstructions of the 1900–2015 Greenland ice sheet surface mass balance using the regional climate MAR model, The Cryosphere, 11, 1015–1033,, 2017. 

Flowers, G. E.: Hydrology and the future of the Greenland Ice Sheet, Nat. Commun., 9, 2729,, 2018. 

Frezzotti, M., Urbini, S., Proposito, M., Scarchilli, C., and Gandolfi, S.: Spatial and temporal variability of surface mass balance near Talos Dome, East Antarctica, J. Geophys. Res.-Earth, 112, F02032,, 2007. 

Hallikainen, M. T., Ulaby, F. T., and Van Deventer, T. E.: Extinction Behavior of Dry Snow in the 18-to 90-GHz Range, in: IEEE Transactions on Geoscience and Remote Sensing, GE-25, 737–745,, 1987. 

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J. N.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049,, 2020. 

IPCC: IPCC Special Report on the Ocean and Cryosphere in a Changing Climate, edited by: Pörtner, H.-O., Roberts, D. C., Masson-Delmotte, V., Zhai, P., Tignor, M., Poloczanska, E., Mintenbeck, K., Alegría, A., Nicolai, M., Okem, A., Petzold, J., Rama, B., and Weyer, N. M., 755 pp., Geneva, 2019. 

Jezek, K. C., Merry, C. J., and Cavalieri, D. J.: Comparison of SMMR and SSM/I passive microwave data collected over Antarctica, Ann. Glaciol., 17, 131–136,, 1993. 

Joshi, M., Merry, C. J., Jezek, K. C., and Bolzan, J. F.: An edge detection technique to estimate melt duration, season and melt extent on the Greenland Ice Sheet using Passive Microwave Data, Geophys. Res. Lett., 28, 3497–3500,, 2001. 

Kargel, J. S., Abrams, M. J., Bishop, M. P., Bush, A., Hamilton, G., Jiskoot, H., Kääb, A., Kieffer, H. H., Lee, E. M., Paul, F., Rau, F., Raup, B., Shroder, J. F., Soltesz, D., Stainforth, D., Stearns, L., and Wessels, R.: Multispectral imaging contributions to global land ice measurements from space, Remote Sens. Environ., 99, 187–219,, 2005. 

Kargel, J. S., Leonard, G. J., Bishop, M. P., Kääb, A., and Raup, B. H. (Eds.): Global land ice measurements from space, Springer-Praxis, Berlin, Heidelberg, 876 pp., 2014. 

Knowles, K., Njoku, G., Armstrong, R., and Brodzik, M.: Nimbus-7 SMMR Pathfinder Daily EASE-Grid Brightness Temperatures, version 1, NASA National Snow Ice Data Center Distributed Active Archive Center, Boulder, Colorado, USA,, 2000. 

Knowles, K., Savoie, M., Armstrong R., and Brodzik, M.: AMSR-E/Aqua Daily EASE-Grid Brightness Temperatures, NASA DAAC at the National Snow and Ice Data Center, Boulder, Colorado, USA,, 2006. 

Kouki, K., Anttila, K., Manninen, T., Luojus, K., Wang, L., and Riihelä, A.: Intercomparison of Snow Melt Onset Date Estimates From Optical and Microwave Satellite Instruments Over the Northern Hemisphere for the Period 1982–2015, J. Geophys. Res.-Atmos., 124, 11205–11219,, 2019. 

Le Meur, E., Magand, O., Arnaud, L., Fily, M., Frezzotti, M., Cavitte, M., Mulvaney, R., and Urbini, S.: Spatial and temporal distributions of surface mass balance between Concordia and Vostok stations, Antarctica, from combined radar and ice core data: first results and detailed error analysis, The Cryosphere, 12, 1831–1850,, 2018. 

Liu, H., Wang, L., and Jezek, K. C.: Wavelet-transform based edge detection approach to derivation of snowmelt onset, end and duration from satellite passive microwave measurements, Int. J. Remote Sens., 26, 4639–4660,, 2005. 

Long, D. G. and Daum, D. L.: Spatial resolution enhancement of SSM/I data, IEEE T. Geosci. Remote, 36, 407–417,, 1998. 

Macelloni, G., Paloscia, S., Pampaloni, P., and Tedesco, M.: Microwave emission from dry snow: a comparison of experimental and model results, IEEE T. Geosci. Remote, 39, 2649–2656,, 2001. 

Macelloni, G., Paloscia, S., Pampaloni, P., Brogioni, M., Ranzi, R., and Crepaz, A.: Monitoring of melting refreezing cycles of snow with microwave radiometers: The Microwave Alpine Snow Melting Experiment (MASMEx 2002–2003), IEEE T. Geosci. Remote, 43, 2431–2442,, 2005. 

Macelloni, G., Brogioni, M., Pampaloni, P., and Cagnati, A.: Multifrequency microwave emission fromthe dome-c area on the east antarctic plateau: temporal and spatial variability, IEEE T. Geosci. Remote, 45, 2029–2039,, 2007. 

MAR model: MAR, available at:, last access: 10 June 2021. 

Mote, T. L.: Greenland surface melt trends 1973–2007: evidence of a large increase in 2007, Geophys. Res. Lett., 34, L22507,, 2007. 

Mote, T. L.: MEaSUREs Greenland Surface Melt Daily 25 km EASE-Grid 2.0, Version 1, NASA National Snow and Ice Data Center Distributed Active Archive Center, Boulder, Colorado, USA,, 2014. 

Mouginot, J., Rignot, E., Bjørk, A. A., van den Broeke, M., Millan, R., Morlighem, M., Noël, B., Scheuchl, B., and Wood, M.: Forty-six years of Greenland Ice Sheet mass balance from 1972 to 2018, P. Natl. Acad. Sci. USA, 116, 9239–9244,, 2019. 

Nash, J. E. and Sutcliffe, J. V.: River flow forecasting through conceptual models, Part 1 – a discussion of principles, J. Hydrol., 10, 282–290, 1970. 

Poe, G. A.: Optimum interpolation of imaging microwave radiometer data, IEEE T. Geosci. Remote, 28, 800–810,, 1990. 

Rignot, E., Velicogna, I., van den Broeke, M. R., Monaghan, A., and Lenaerts, J. T. M.: Acceleration of the contribution of the Greenland and Antarctic ice sheets to sea level rise, Geophys. Res. Lett., 38, L05503,, 2011. 

Robinson, W. D., Kummerow, C., and Olson, W. S.: A technique for enhancing and matching the resolution of microwave measurements from the SSM/I instrument, IEEE T. Geosci. Remote, 30, 419–429,, 1992. 

Steffen, K., Abdalati, W., and Stroeve, J.: Climate sensitivity studies of the Greenland ice sheet using satellite AVHRR, SMMR, SSM/I and in situ data, Meteorol. Atmos. Phys., 51, 239–258,, 1993. 

Steffen, K., Box, J. E., and Abdalati, W.: Greenland Climate Network: GC-Net, in: CRREL 96–27 Special report on glaciers, ice sheets and volcanoes, trib. to Meier, M., edited by: Colbeck, S. C., US Army Cold Regions Research and Engineering Laboratory (CRREL), Hanover, New Hampshire, USA, 98–103, 1996. 

Steffen Research Group: Greenland Climate Network (GC-Net), available at:, last access: 9 June 2021. 

Steiner, N. and Tedesco, M.: A wavelet melt detection algorithm applied to enhanced-resolution scatterometer data over Antarctica (2000–2009), The Cryosphere, 8, 25–40,, 2014. 

Stogryn, A.: Estimates of brightness temperatures from scanning radiometer data, IEEE T. Antenn. Propag., 26, 720–726,, 1978. 

Stroeve, J., Maslanik, J., and Xiaoming, L.: An Intercomparison of DMSP F11- and F13-Derived Sea Ice Products, Remote Sens. Environ., 64, 132–152,, 1998. 

Tedesco, M.: Snowmelt detection over the Greenland ice sheet from SSM/I brightness temperature daily variations, Geophys. Res. Lett., 34, L02504,, 2007. 

Tedesco, M.: Assessment and development of snowmelt retrieval algorithms over Antarctica from K-band spaceborne brightness temperature (1979–2008), Remote Sens. Environ., 113, 979–997,, 2009. 

Tedesco, M.: Remote sensing of the cryosphere, John Wiley & Sons, UK, 2014. 

Tedesco, M. and Fettweis, X.: Unprecedented atmospheric conditions (1948–2019) drive the 2019 exceptional melting season over the Greenland ice sheet, The Cryosphere, 14, 1209–1223,, 2020. 

Tedesco, M. and Wang, J. R.: Atmospheric correction of AMSR-E brightness temperatures for dry snow cover mapping, IEEE Geosci. Remote S., 3, 320–324,, 2006. 

Tedesco, M., Abdalati, W., and Zwally, H. J.: Persistent surface snowmelt over Antarctica (1987–2006) from 19.35 GHz brightness temperatures, Geophys. Res. Lett., 34, L18504,, 2007.  

Tedesco, M., Brodzik, M., Armstrong, R., Savoie, M., and Ramage, J.: Pan arctic terrestrial snowmelt trends (1979–2008) from spaceborne passive microwave data and correlation with the Arctic Oscillation, Geophys. Res. Lett., 36, L21402,, 2009. 

Tedesco, M., Fettweis, X., Mote, T., Wahr, J., Alexander, P., Box, J. E., and Wouters, B.: Evidence and analysis of 2012 Greenland records from spaceborne observations, a regional climate model and reanalysis data, The Cryosphere, 7, 615–630,, 2013. 

Torinesi, O., Fily, M., and Genthon, C.: Variability and Trends of the Summer Melt Period of Antarctic Ice Margins since 1980 from Microwave Sensors, J. Climate, 16, 1047–1060,<1047:vatots>;2, 2003. 

Ulaby, F. T., Moore, R. K., and Fung, A. K.: Microwave remote sensing: Active and passive, Volume 3 – From theory to applications, Norwood, MA, 1986. 

van den Broeke, M.: Strong surface melting preceded collapse of Antarctic Peninsula ice shelf, Geophys. Res. Lett., 32, L12815,, 2005. 

van den Broeke, M. R., Enderlin, E. M., Howat, I. M., Kuipers Munneke, P., Noël, B. P. Y., van de Berg, W. J., van Meijgaard, E., and Wouters, B.: On the recent contribution of the Greenland ice sheet to sea level change, The Cryosphere, 10, 1933–1946,, 2016. 

Wiesmann, A. and Mätzler, C.: Microwave Emission Model of Layered Snowpacks, Remote Sens. Environ., 70, 307–316,, 1999. 

Zwally, H. J. and Fiegles, S.: Extent and duration of Antarctic surface melting, J. Glaciol., 40, 463–475,, 1994. 

Short summary
We use a new satellite dataset to study the spatiotemporal evolution of surface melting over Greenland at an enhanced resolution of 3.125 km. Using meteorological data and the MAR model, we observe that a dynamic algorithm can best detect surface melting. We found that the melting season is elongating, the melt extent is increasing and that high-resolution data better describe the spatiotemporal evolution of the melting season, which is crucial to improve estimates of sea level rise.