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

. Surface melting is a major component of the Greenland ice sheet surface mass balance, affecting sea level rise 10 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-to date) and high temporal resolution (daily). However, a major limitation of PMW datasets has been the relatively coarse spatial resolution, being historically of the order of tens of kilometres. Here, we use a newly released PMW dataset (37 GHz, horizontal polarization) 15 made available through the 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 through data collected by Automatic Weather Stations (AWS) and the outputs of the MAR regional climate model. We found that sporadic melting is well captured using a dynamic algorithm based on the outputs of Microwave Emission Model of Layered Snowpack (MEMLS) model while a fixed threshold 20 of 245K is capable of detecting persistent melt. Our results indicate that, during the reference period 1979 – 2019 (1988 – 2019), surface melting over the ice sheet increased in terms of both duration, up to 4.5 (2.9) days 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 has up (2.5) days earlier and 7.0 (3.9) days 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 through a semi- 25 variogram approach. We found that the enhanced product is more sensitive to local scale processes, hence confirming the potential of this new enhanced product for monitoring surface melting over Greenland at a higher spatial resolution than the historical products and monitor its impact on sea level rise. This offers the opportunity to improve our understanding of the processes driving melting, to validate modelled melt extent at high resolution and potentially to assimilate this data in climate models.


Introduction
The Greenland ice sheet is the largest ice mass in the Northern Hemisphere with a glaciated surface area of about 1,800,000 km 2 , a thickness up to 3 km, and a stored water volume of about 2,900,000 km 3 , enough to rise the mean sea level by about 7.2 m (Aschwanden et al., 2019). In this regard, estimating mass losses from Greenland is crucial for better understanding climate system variability and the contribution of Greenland to current and future sea level rise. According to 35 data from the Gravity Recovery and Climate Experiment (GRACE) satellite mission, which records changes in Earth's gravitational field, Greenland lost mass at an average rate of 278±11 Gt y -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 was also accelerating at a rate of 21.9±1 Gt y -2 over the period 1992-2010 (Rignot et al., 2011) thus indicating that monitoring the Greenland ice sheet together with the Antarctic ice sheet is crucial to assess the impact of global warming on sea level rise and the global water 40 balance (Kargel et al., 2005;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. difference between accumulation (Frezzotti et al., 2007) and ablation, surface melting plays a crucial role, affecting direct loss through export of surface meltwater to the surrounding oceans and though the feedbacks between supraglacial, englacial and subglacial processes and their influence on ice dynamics (e.g., Fettweis et al., 2005Fettweis et al., , 2011Fettweis et al., , 45 2017van 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 et al., 2009;Tedesco, 2009;Fettweis et al., 2011). The PMW-based algorithms are based 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 50 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 (see for instance 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 55 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;Aschraft and Long, 2006;Macelloni et al., 2007;Tedesco et al., 2007;Kouki et al., 2019;Tedesco and Fettweis, 2020). 60 The capability of PMW sensors to collect useful data during both day and night and in 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 timeseries 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 on the order of tens of km). This 65 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 allow a better comparison between in-situ measured quantities and satellite-derived estimates, reducing uncertainties in the satellite products and allowing for potential 70 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 improve 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 our results of a study in which surface melting over Greenland is estimated making use of a recently released product developed within the framework of a NASA Making Earth System Data Records for Use in Research 75 Environments (MeASUREs) project (https://nsidc.org/data/nsidc-0630). The product contains daily maps of PMW Tb generated at an enhanced spatial resolution of a few kilometers (depending on 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;Knowles et al., 2006) were based on a "drop in the bucket" approach, in which the gridded value was 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 80 reconstruction algorithm adopted to build the Tb maps makes use of the so-called effective measurement response function (Long et al., 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 for meeting the 85 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 Section 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 Ka band (37 GHz), horizontal polarization in view of its 90 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 the Tb at Ka band are distributed at the highest spatial resolution of 3.125 km (Brodzik et al., 2018). The atmospheric effect on 37 GHz Tb is higher than in case 19 GHz at low values of Tb. However, when considering higher values of Tb, the difference of the atmospheric effect between 37 GHz and 19 GHz Tb decreases (Tedesco and Wang, 2009). We, 95 then, focus on assessing whether the noise introduced by the gridding algorithm might limit the application of the enhanced dataset to mapping surface meltwater. Then, we focus our attention on testing and assessing existing approaches to deriving melt from PMW data and propose an update on 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 (AWS) in terms of 100 melting timing and with the outputs of the regional climate model Modèle Atmosphérique Régional (MAR; Fettweis et al., 2017) in terms of melting timing and extent. Lastly, we focus on the analysis of melting patterns and trends over study period and investigate the information content in the enhanced resolution dataset through a semi-variogram analysis.

Enhanced resolution passive microwave data 105
We use Ka band (37 GHz), horizontal polarization Tb data produced within the framework of a NASA MeASUREs project and distributed at the spatial resolution of 3.125 km (Brodzik et al., 2018) 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, SSMI/-F11, SSM/I-F13 and the special sensor microwave imager/sounder SSMI/S-F17 because of its higher orbit stability (http://www.remss.com/support/crossing-times/). Currently, the product time series 110 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, together with other information. More information can be found at https://nsidc.org/data/nsidc-0630/versions/1. Historical gridding techniques for PMW spaceborne datasets (Armstrong et al., 1994, updated yearly;Knowles et al., 2000;Knowles et al., 2006) are relatively simplistic and were produced on grids (Brodzik and Knowles, 2002;Brodzik et al., 2012) 115 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 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;120 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 improving the spatial resolution of surface Tb fields (Farrar and Smith, 1992;Long and Daum, 1998). More information about the product can be found at https://nsidc.org/data/nsidc-0630. An example of Tb maps at 37 GHz, horizontal polarization, in the case of both the coarse and enhanced resolution products over 125 Greenland on 16 July 2001 is reported in Figure 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 (blue) and enhanced (red) PMW Tb (again at 37GHz, horizontal polarization) for the pixel containing the Swiss Camp station. From the figure we observe that the two timeseries are highly consistent with each other, with a mean difference of 0.895 K and standard deviation of 4.89 K, 130 indicating that the potential noise introduced by the enhancement process is not a major issue. Yet, differences do exist, as in the case of 3 April 2012 (DOY 93), when the enhanced product suggests the presence of melting while the coarse product does not ( Figure 2). This is likely due to the different spatial resolution 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.
To compare our results with precedent PMW surface melting products, we perform our calculations also to the widely 135 used PMW dataset by Mote (2014). The dataset consists in 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 in 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 simulating the Tb change for every grid cell, each year, through a microwave emission model (Mote, 2014). The dataset covers the temporal window 1979 -2012 and is available at the NSIDC website 140 (https://nsidc.org/data/NSIDC-0533/versions/1).

Greenland air temperature data
In order to assess the results obtained from PMW data, we use in-situ data collected by AWS 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. 145 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 17 selected stations reported in Table 1. We considered a validation period from 2000 (when all the considered AWS were in operation) to 2016 and used daily averaged values. More information about the GC-Net dataset can be found at http://cires1.colorado.edu/steffen/gcnet/. 150
MAR is a modular atmospheric model that uses the sigma-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 155 have been assessed over Greenland (e.g., Fettweis et al., 2005;Fettweis et al., 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. 160 In this study, we use the output from MAR version v3.11.2 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), 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 165 along the vertical profile, following Fettweis et al. (2007).

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, computed in different ways depending on the 170 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.
Aschraft and Long (2006) proposed a threshold based on dry (winter) and wet snow Tb as Tc= αM+(1-α) Twet where 175 M is the average of winter Tb (January and February), Twet 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 snowpack. Similarly, a method based on a fixed threshold (set to 245 K and derived from the outputs of electromagnetic model) above which melting is assumed to be occurring was proposed in . Several other studies have been detecting melting when Tb values exceed the mean winter value plus an additional value ∆Tb (M+∆Tb approaches) associated 180 with the insurgence of liquid water within the snowpack: Torinesi et al. (2003) proposed a value of ΔTb=Nσ with M and σ (standard deviation of the timeseries) varying in space (specific pixel) and time (specific year) but fixed 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 185 Microwave Emission Model of Layered Snowpack (MEMLS) electromagnetic model (Weisman and Matzler, 1999). In this case, an ensemble of outputs is generated by MEMLS model by varying the inputs (e.g., correlation length, LWC, density, etc.). 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: 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 values of LWC, potentially making the mapping of the wet snow areas inconsistent in terms of LWC values. For example: a snowpack with small grain size will require a relatively larger amount of LWC with respect to a snowpack with larger grain size for the Tb values to increase by 30 K. Or, from a complementary point of view, an increase of 30 K due to presence of 195 liquid water in the case of a snowpack with relatively coarse grains will correspond to a lower value of LWC than an increase occurring in a snowpack with 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 (R 2 =0.92). The Tb threshold value computed in this case can, therefore, be written as follows: 200 where (γ, ω) assume the values of 0.48 and 128 K.
Here, we focus five approaches: the M+∆Tb approach choosing ∆Tb equal to 30K and, to test the sensitivity to Zwally and Fiegles (1994), 35K and 40K (M+30, M+35 and M+40 from here on), the algorithm based on MEMLS in case of LWC=0.2% (referred simply as MEMLS from here on for brevity) and the 245 K fixed threshold (245K from here on). We selected M+∆Tb 205 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), proved in previous studies (Tedesco, 2009). We selected also the 245K to test a more conservative approach aimed to detect 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.

Inter-sensor calibration 210
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   2). Several approaches have been proposed in the literature to address this issue for the historical, coarser spatial resolution 215 gridded datasets. For example, Jezek et al. (1993) compared SMMR and SSM/I over the Antarctic ice sheet for K and Ka bands (~ 19 GHz and ~ 37 GHz) for both horizontal and vertical polarizations. Steffen et al. (1993) proposed an approach focusing over Greenland for the K-band; Abdalati et al. (1995) derived relations between SSM/I observations for the F08 and F11 platforms over Antarctica and Greenland for 19.35 GHz, 22.2 GHz and 37 GHz. Dai et al. (2015) intercalibrated SMMR, SSM/I (F08 and F13) and SSMI/S (F17) over snow covered pixels in China and SMMR, SSM/I and AMSR-E over the whole 220 Earth surface sampling hot and cold pixels.
Given the novelty of the Tb products used here and the absence of 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) 225 of the linear regression 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 i-th day we first compute mi, qi and the 230 coefficient of determination for the linear regression of Eq. (4) (R 2 i), and then we average them according to Eq. (5) and Eq. (6).
This choice assigns higher values to the weights obtained from pairs of data with higher correlation. In the second 235 method, we consider all values for all days when data from both platforms are available and then evaluate m and q through a linear regression fitting procedure based on least-square fitting. Using the estimated values of m and q, 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 matching 240 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: where hi is the absolute frequency of the i-th value of Tb of the two sensors (A and B) considered. The NSE is usually applied in calibration/validation procedures to assess the matching between measured and modelled quantities, as in Subsection 4.2. 245 After the application of linear relations found using Eq. (4) through (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 obtained as: where Di is the absolute difference between the two histograms A and B for the i-th value of Tb. Then, we sum the differences over the total number of pixels and compute the relative variation as follows: 250 where Doriginal and Dcorrected are, respectively, the summations of Di before and after the calibration. The relative variation d can range from -∞, indicating worsening in matching of the histograms, to 1, indicating a perfect matching of the histograms after the intercalibration.

Spatial autocorrelation: the variogram analysis 255
Variogram analysis is generally adopted in geostatistical analyses to evaluate 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 260 empirical variogram as where  is the semi-variance, N() is the number of pair measurements (i,j) spaced by distance  and xi and xj are the values of the i-th and j-th 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 265 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.

Inter-sensor calibration of enhanced resolution passive microwave data
At first, we show the results obtained for the inter-sensor calibration of the selected satellite constellation. The 270 overlapping periods for the different sensors are the following: SMMR and SSM/I-F08 overlap between 9 July 1987 and 20 August 1987 for a total of 22 days (one every two days as sensed by SMMR sensor); F08 and F11 overlap between 3 December 1991 and 18 December 1991 for a total of 16 days, F11 and F13 overlap between 3 May 1995 and 30 September 1995 for a total of 76 days; and F13 and F17 overlap for the period 1 March 2008 -10 December 2008 for a total of 71 days. In Figure 3 we show the scatter plots of the data used for the linear regression for Greenland for both evening and morning passes for the 275 SMMR and SSM/I-F08 sensors, reporting values of m, q and R 2 . 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, potentially having frozen snow (low values of Tb) early in the morning and late at night (SMMR ascending and SSMI/-F08 descending passes) and 280 presence of liquid water (high values of Tb) during the day. We report in Table 3 average values of the difference between pairs of Tb data and values of the NSE coefficient for the histograms of the same pairs. In Table 4 we report the values for slope and intercept obtained from the linear regression analysis of enhanced PMW Tb (37 GHz, H-pol.) over Greenland for SMMR vs. SSM/I-F08, F08 vs. F11, F11 vs. F13 and F13 vs. F17, together with the R 2 values and values of d computed according to Eq. (6). In the case of SSM/I and SSMI/S, R 2 values are higher, mostly around 0.98. In Figure 4 we also show 285 examples of histograms in the case of the SMMR and SSM/I F08 sensors. Large differences are obtained in the case of the SMMR and SSM/I-F08, for both evening and morning passes, likely because of the difference in overpass time and the presence/absence of melting in some of the scenes observed by one sensor but not present in the other (Table 3). On the other hand, in the case of 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 showing an average difference slightly larger than 1 K, consistent with previous results obtained by Abdalati et 290 al. (1995) in the case of the 25 km resolution data) together with NSE values extremely close to 1 (Table 3). Still in the case of SMMR and SSM/I-F08, the higher average difference (ranging between -3.4 K to -4.3 K) and the relatively lower NSE values (ranging between 0.89 and 0.96) show that these sensors show 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, SSM/I F13 -SSMI/S F17 datasets as, in this case, the linear correction worsened the agreement between the two sets of measurements. We 295 applied the correction coefficients obtained with the second method according to the higher relative improvement for the evening pass.

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 daily averaged from AWS as an index of surface melting (Braithwaite and Oelsen, 1989) and with the 300 liquid water content simulated by the regional climate model MAR. We first evaluate performances at local scale (at the specific locations of the selected AWS), comparing the number and the concomitance of melting days according to PMW and the ground truth reference. Then, according to the results obtained, we focus on MEMLS and 245K algorithms to evaluate at ice sheet scale the capability of the two approaches in describing the surface melt extent.

Assessment with AWS data 305
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 because of 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 of 0ºC, -1°C and -2°C, as in Tedesco (2009). We assessed the performance of the PMW-based algorithm by defining commission and omission errors. Commission error occurs when melting is detected by 310 PMW data but not by AWS data and omission error occurs when melting is detected by AWS data but not by PMW. The results of the error analysis are summarized in Table 5Table 5 (Table S2, Table S3 and Table S4) together with a map of the AWS network ( Figure S1). Table 5 also reports as more 315 general performance indicator the sum of the two errors, computed for each AWS case (C+O). We also report an average value of all the C+O (for both AWS and MAR assessments, presented in the next subsection) for each PMW algorithm (C+O Mean) as synthetic index of performance. Our results indicate that the 245K algorithm shows the lowest commission error (between 0.31% and 0.63%) and the highest omission error (5. 38%-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 320 error is achieved in the case of the M+30, M+35 and M+40 thresholds, particularly for the Humboldt and GITS stations (North-West Greenland), where the commission error is up to one order of magnitude larger than in the case of MEMLS and 245K algorithms for every ground-truth reference (e.g. from 0.70% for MEMLS and 0.09% for 245K to 5.63% for M+30, in case of air temperature equal to 0°). Moreover, we note in the case of the MEMLS algorithm the lowest omission error in 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 245K). The 325 coarse resolution dataset presents a commission error between 1.02% and 1.74% and an omission error between 4.06% and 7.12%, confirming the capability of historical data in detecting 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, respectively decreasing by about 1% and increasing by 3% when considering threshold values from 0°C to -2°C. 330 In order to better understand the sources of the relatively high values of the commission errors at some locations, we show in Figure 5 the timeseries of air temperature and Tb at 37 GHz H-pol. at three selected stations: a) Summit, b) Humboldt and c) Swiss Camp for the year 2005. The threshold values obtained with the different detection algorithms are also plotted as horizontal lines (black) as well as the 0°C air temperature threshold (magenta). We selected these three locations as they are representative of three different environmental and melting conditions. The timeseries recorded at Summit station (Figure 5a) 335 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 timeseries is typical of a location where melting is generally absent. The Tb timeseries collected in correspondence of Humboldt location (Figure 5b) shows a strong and sudden peak starting on 20 July, when the air temperature average is about -0.5 °C (detected by -1°C and -2°C air temperature thresholds). This event is detected by all algorithms. Nevertheless, the M+30 (and similar algorithms) 340 indicate the potential presence of melting also for the period preceding the July melting (between 17 June and 17 July). This melting is not confirmed by other algorithms or by the AWS analysis, suggesting that the threshold value used for these algorithms might be too low. Lastly, melting clearly occurs in the case of Swiss Camp (Figure 5c), characterized by the sharp and substantial increment of Tb beginning around mid-May. For this case, all algorithms detect melting, with the MEMLS providing the lowest threshold and the 245K fixed threshold being the most conservative. The computed rough estimation of 345 the average emissivity for the period 17 June -17 July (as Tb divided by the recorded air temperature) also suggests that melting is not occurring in the considered period in Humboldt case, presenting an average emissivity even lower than in Summit case. Figure 6 shows maps of surface melt extent obtained using the different approaches for July 13th, 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, while the 245K fixed-350 threshold approach shows, as expected, the most conservative estimates. As mentioned, the threshold algorithms for ΔTb (M+30, etc.) rely on a fixed ΔTb value, which could produce errors if there is a large seasonal range in Tb 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 takes into account the snow conditions before melting and, at the same time, follows a more consistent approach with respect to the 355 amount of LWC detected in the snowpack.

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 (MAR5cm) and the top 1m (MAR1m) of 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 to which we 360 assume melting is occurring to 0.2% for both depths. We selected two different depths for our analysis so 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 modelled LWC indicates that the more conservative 365 approaches (i.e., 245K) perform better when considering the case of MAR5cm. In fact, the 245K 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%, slightly lower than in case of MEMLS (7.56%). In the case of the top 1 m, all the algorithms present similar performances on average, with the best performance obtained again in the case of the 245K (C+O=5.53%). However, all the M+ΔTb algorithms present the same issue of larger commission error than MEMLS and 245K (e.g., from 0.99% for MEMLS and 0.26% for 245K to 4.62% for 370 M+30) in North-East Greenland (e.g., Humboldt and GITS stations, see Table S2 and Table S3 in the supplementary material).
This confirms the results we obtained from the comparison with AWS data, pointing out the overestimation of melting in some dry areas by M+ΔTb. For both the MAR1m and MAR5cm cases, for all the considered algorithms, we find a high commission error in the cases of 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 in Figure 7 further insights into the differences 375 between the PMW Tb and MAR outputs. Figure 7a and b show, respectively, the timeseries of LWC averaged over the first 5cm (MAR5cm) and 1m (MAR1m) obtained from MAR at the Swiss Camp site. In Figure 7c we report the Tb timeseries and the daily average air temperature (threshold values reported as horizontal lines). We first note an early melt event (labeled LWC=0.046% in Figure 7b for the 108th day of the year) detected by PMW MEMLS algorithm and at the AWS station but apparently undetected in MAR1m. A closer look to the time series shows that in fact MAR1m does estimate a LWC of 0.046% 380 on this day while MAR5cm 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 snowpack, as a consequence of the approximation of the electromagnetic outputs imposed by the linear fitting. We also note a melt event (labeled with LWC=0.5% in Figure 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 timeseries reveals a small peak, but the signal is not strong enough to 385 exceed any threshold. This corresponds to a rainfall event (simulated by MAR) suggesting that the sensitivity to liquid clouds of the 37 GHz channel could mask some melt events. Moreover, at the end of the melting season the Tb appears to be slightly lower than January/February average, possibly because of an increment in grain size after refreezing, leading to a lower emissivity.
The results discussed above (together with results from the comparison with AWS data) suggest that 245K is the most 390 conservative among the approaches we tested, providing the lowest (highest) commission (omission) error but being unable to detect sporadic melt events. On the contrary, the MEMLS and M+Tb algorithms can detect sporadic melt events and present lower omission error compared to 245K. However, M+Tb algorithm overestimate melting in some dry areas (North West of the ice sheet), suggesting melting when it is not actually occurring. Contrarily, MEMLS algorithm is not affected by the large commission error in dry areas, presents the lowest omission error in Swiss Camp area (together with M+30) and is still sensitive 395 to low levels of LWC. Considering the average error (C+O Mean in Table 5), the MEMLS algorithm shows the best performance (6.66%). In view of the presented analysis and the different sensitivity to surface and subsurface melting, in the following we focus on the 245K 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 Figure  400 8 we show the timeseries of melt extent integrated over the whole ice sheet for two selected years ( (a) 1983 and (b) 2005, selected randomly to present an example of SMMR and SSM/I cases) estimated according to MAR5cm and MAR1m together with the timeseries of the melt extent from the PMW data. The analogous figure for the coarse resolution dataset is reported in the supplementary material ( Figure S2). For each year, we compute the daily melt extent for the period 1 May to 15 September and use the Nash-Sutcliffe Efficiency (NSE) coefficient (Nash and Sutcliffe, 1970), described in Section 3, for a 405 quantitative analysis. Here, we remind the reader that NSE can assume values in the interval (-∞,1]. A perfect match between the two timeseries is achieved when the NSE value is 1. Values of NSE in the [0,1] interval indicate that the modelled variable is a better predictor of the measurements than the mean. If NSE is a negative number, the mean of the measured data describes the timeseries better than the modelled predictor. Here, we chose NSE=0.4 as efficiency threshold, considering that we compute melt extent at daily timescale and from datasets at two different resolutions (i.e., resulting in an intrinsic bias related 410 to the different pixel size). We compared the timeseries of melt extent (ME) obtained using the 245K algorithm with MAR1m (245K vs. MAR1m) and the MEMLS melt extent with MAR5cm (MEMLS vs. MAR5cm) due to the expected differences in sensitivity to detect persistent and sporadic melting between 245K and MEMLS, respectively. We compare the melt extent obtained from the coarse resolution dataset with MAR5cm, according to the similarity with MEMLS in terms commission/omission error. We report NSE coefficients computed for the 41-year (34-year in case of the coarse resolution 415 dataset) period in Table 6. At first, we notice that for the [1979][1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990][1991][1992]  obtained with 245K appears to better follow the temporal variability of melt extent from MAR during the melting season but still presenting a strong underestimation ). On the other hand, MEMLS-derived timeseries better matches the MAR-derived one, showing a largely satisfactory Nash-Sutcliffe Efficiency coefficient (0.782). In these two years the NSE computed for the coarse resolution case is negative. The magnitude of the errors is lower than in case of 245K algorithm, 435 indicating a weaker underestimation of the melt extent. This can be a consequence of the coarser resolution, lacking in capturing melting areas at the edges of the Greenland ice sheet.
In summary, 245K threshold, even if presenting acceptable results in terms of commission and omission error considering both AWS and MAR comparison, is too high to fully capture the melt extent everywhere over the ice sheet.
Contrarily, we found that MEMLS algorithm is suitable in capturing the evolution of melting over the Greenland ice sheet. 440 The comparison with 25 km historical surface melting dataset shows the same underestimation issue of 245K. Even showing lower errors than in case of 245K, we did not find acceptable values of NSE.

Surface melting trends
Here, we report results concerning trends of melt duration, length of the melting season and melt extent. We define melt duration (MD) as the total number of days when melting is detected. We compute trends of MD over the whole ice sheet 445 (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 timeseries; Tedesco et. al, 2007). Lastly, we define melt onset date (MOD) and melt end date (MED) as the first day when melting occurs 450 for two days in a row (MOD) and when melting does not occur for at least 2 days 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). (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% y -1 for MEMLS and 0.47% y -1 for 245K). The obtained trends computed using the 25 km dataset are equal to 1.31% y -1 .for the period 1979 -2012 and 0.91% y -1 for the period 1988 -2012, larger than 470 the trends computed for MEMLS (1.03% y -1 and 0.77% y -1 ) but smaller than the ones computed for 245K (1.50% y -1 and 1.23% y -1 ). In the case of MI (Figure 9c), we also found positive statistically significant trends of 9.166 x 10 5 km 2 d y -1 (MEMLS) and 5.862 x 10 5 km 2 d y -1 (245K) for the complete timeseries. When considering the reduced reference period, we found a 95% statistically significant trend of 5.726 x10 5 km 2 d y -1 only in case of MEMLS. The trends computed for the 25 km resolution dataset results equal to 0.999 x10 6 km 2 d y -1 (1979 -2012) and 1.019 x10 6 km 2 d y -1 (1988 -2012), smaller than 475 MEMLS (1.428 x10 6 km 2 d y -1 and 1.325 x10 6 km 2 d y -1 , respectively) and 245K for the period 1979-2012 (0.999 x10 6 km 2 d y -1 ) but larger than 245K for the period 1988 -2012 (1.019 x10 6 km 2 d y -1 ). Lastly, we report in Figure 9d the MOD and MED averaged spatially over the pixels with 95% significant trends in Figure 11. We found that average MOD (crosses in Figure10d) presents similar trends for 245K and MEMLS considering both the entire and shortened time series equal to -0.546 d y -1 and -0.273 d y -1 , respectively, in case of 245K and -0.404 d y -1 and -0.254 d y -1 in case of MEMLS. For the reference period 1979 480 -2012 (1988 -2012)  In Figure 10 we show the trends of MD, MOD and MED on a pixel-by-pixel basis for the complete time series . We found that the trend in MD exhibits the highest statistical significance (in terms of 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 MD in all pixels (Figure 10a and b), with higher values moving towards the coastline, maxima in the ablation zone of the Jakobshavn Isbrae (2.40 d y -1 ) for 245K and 2.66 d y -1 for MEMLS) and minima in high altitude areas. We averaged the 495 statistically significant trends, finding an average of 0.468 d y -1 for the 245K algorithm, and of 0.697 d y -1 in the case of the MEMLS. In case of MOD and MED, we found a lower number of statistically significant pixels. The statistically significant pixels exhibit a negative trend for MOD (Figure 10c and d) and a positive trend for MED (Figure 10e and f), with the melting season starting on average 0.694 d y -1 earlier and ending 0.680 d y -1 later according to 245K algorithm (0.360 d y -1 earlier and 0.909 d y -1 later for MEMLS). We point out that the average of the statistically significant trends is generally higher than the 500 trends computed at ice sheet scale since we computed the average over the statistically significant pixels only.

Spatial information content
In order to investigate the spatial information content of the enhanced resolution PMW data with respect to the coarser one, we also performed a variogram-based analysis of MD estimated from the two products when using either the 245K or the MEMLS algorithms. We point out that knowledge of scales is imperative for improving our understanding of the observed 505 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 MD estimated from the MEMLS and 245K algorithms for the months of May through August of 2012 when using either the enhanced or the coarse 510 resolution products. We also performed the same analysis applied to MD estimated using LWC modelled 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 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 Figure 12, where we show the empirical (blue crosses) and modelled (red line) semivariograms for Greenland MD computed applying the MEMLS and 245K algorithms to both 25 km and 3.125 km resolution 515 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 in case of MD obtained according to MAR1m and MAR5cm approaches (an analogous representation of Figure 12 is reported in the supplement in Figure S5). At first, we note that R 2 values of the fitting for the modeled variograms are consistently higher in the case of enhanced resolution data, suggesting that enhanced resolution data might be more suitable for a variogram-based analysis. For the coarse resolution data, we found R 2 values of comparable 520 magnitude with the enhanced resolution case only in May. When computing the spherical fitting of the empirical variograms of MD from MAR, we found for each case considered similar R 2 values (between 0.118 and 0.484). The values for the range in the case of the 3.125 and 25 km products are similar in May for the 245K algorithm (on the order of ~ 200 km) but they appear different in the case of the MEMLS algorithm, when the enhanced product shows a lower value of ~ 170 km against ~ 270 km in the case of the coarse product. This could be due to the fact that the MEMLS algorithm is more sensitive to sporadic 525 melting and when applied to the enhanced Tb dataset it allows to detect melting driven by processes whose scale cannot be captured by the coarser nature of the historical dataset. In case of MAR, the value of the range is lower in case of MAR5cm (187.70 km) than in case of MAR1m (199.17 km), suggesting again the affinity of MEMLS algorithm with melting strictly confined in the very first layer of snowpack. As the melting season progresses, the variograms of the coarse dataset shows either similar values for the range or a poor fit of the experimental variogram. Instead, in the case of the enhanced product, the 530 values of the range tend to decrease up to July and increase again in August. We found the same temporal variability in case of MAR1m, while in case of MAR5cm we found that the range increases until July and decreases in August. 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 nugget effect, we found larger values from the MAR outputs than in case of PMW. The decrease in the spatial 535 autocorrelation length in the case of 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, new snowfall, etc.) 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, 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 540 substantial melting occurs) is of the same order of magnitude of 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, together with the draining of the supraglacial lakes and rivers might justify the increase in the range values computed for this month. Our preliminary results, therefore, point to an increased information content of the enhanced spatial PMW product with respect to the historical, coarse one, offering the opportunity to better capture the spatial details of how surface melting evolved over the Greenland ice sheet 545 over the past ~ 40 years. Further analysis will help to shed light on the processes responsible for the recent acceleration of surface melting.

Conclusions and future work
We applied threshold-based melt detection algorithms to the 3.125 km resolution 37 GHz horizontal polarization PMW Tb to assess the skills of the enhanced PMW product to detect surface melting over the 1979-2019 period over the 550 Greenland Ice Sheet. As the product is composed of data acquired by different sensors onboard of different platforms, we first developed a cross calibration among all the sensors. Then, we compared surface melting detected from PMW data with that estimated from AWS air temperature data and the outputs of the regional climate model MAR. We found that the algorithm making use of a fixed threshold value on Tb values (245K) and the one based on the outputs of an electromagnetic model were the most suitable to detect persistent (245K) and sporadic (MEMLS) melting. Overall, we found that that MEMLS algorithm 555 showed the best performance (lowest commission and omission errors). We compared surface melting detected from PMW data with the one estimated by the MAR model when considering the two cases of integrating LWC over the top 5cm and 1m, respectively. We selected these two depths to study those conditions when melting occurs sporadically (5 cm) or persistently (1 m). We obtained good matching (i.e., NSE>0.4 or, at least, positive) in most of the years from 1992-2019 when comparing MEMLS derived melt extent with MAR liquid water content in the first 5 cm of snowpack. On the other hand, we found bad 560 matching in the period 1979-1992, possibly due to differences in sensor characteristics. In the case of melt extent retrieved by 245K, we found a strong underestimation of melt extent (largely negative values of NSE) from 1979 to 1987 likely because of the lower values of "wet" Tb in case of SMMR data, slightly improving from 1993 to 2019 but still negative. Accordingly, the results obtained applying MEMLS approach are more reliable than in case of 245K algorithm when considering the period 1979 -2019. When comparing with the coarse resolution dataset, we found that the melt extent timeseries derived from the 565 enhanced resolution data using MEMLS better agree with MAR simulations than the ones obtained using the 25 km resolution data.
After assessing the outputs of the PMW-based algorithms, we studied the MOD, MED, MMD and MMS for the period 1979 -2019. According to MEMLS algorithm, we found that the melting season has begun 0.404 (0.254) days earlier every year between 1979-2019 (1988-2019) and has ended 0.708 (0.396) days later every year between 1979-2019 (1988-570 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 MMD has increased every year by 0.451 d y -1 (0.291 d y -1 ) during the period 1979-2019 (1988-2019). We found differences in trends computed using the 25 km resolution with respect to the enhanced resolution product for the reference periods 1979 -2012 and 1988 -2012, possibly because of the different rationale behind the melt detection algorithms and the higher level of detailing of the enhanced resolution dataset. When we performed a spatial analysis of the trends for the 575 melt onset dates and duration, we found that the areas where the number of melting days has been increasing are mostly located in West Greenland. The MMS presents positive trends as well, with an increment of 0.69% (0.36%) every year respect to the Greenland ice sheet surface since 1979 (1988).
Finally, we explored the information content of the enhanced resolution dataset with respect to the one at 25 km and the MAR outputs through a semi-variogram approach. The results obtained showed a better fitting of the modelled spherical 580 function to the empirical semi-variogram in case of the 3.125 km and MAR maps of MD. Our analysis suggests that the enhanced resolution product is sensitive to local scale processes, with higher sensitivity in case of 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 using 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 585 and longwave radiation and complement this analysis with our previous work focusing on understanding the changes in atmospheric patterns that have been promoting enhanced melting in Greenland over the recent decades (Tedesco and Fettweis, 2020). Assessed the capability of this dataset and method in observing temporal trends, a further development can include a combination of the enhanced PMW product 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 590 Canadian Arctic Archipelago, the Himalayan Plateau and the Antarctic Peninsula, where the enhancement in spatial resolution can be fully exploited.       Table 4: Slope (m) and intercept (q) obtained from the linear regression analysis between the selected couples of satellites enhanced PMW Tb at 37 GHz, horizontal polarization over Greenland. The subscripts refer to the case when the coefficients are weighted by means of the R 2 (case 1, see Eq. (5) and Eq. (6)) or not (case 2). In the Table, we also report the values for the R 2 as well as the values of d computed according to Eq. (9).
Figure 4: Histograms of Tb before and after the application of the intercalibration relations, Greenland. Relations are applied for both evening (a) and morning (b) passes, reporting the histograms of the data and the distance (absolute value of the difference as in Eq. (8)) between the histograms for original data. The first column represents the uncorrected data, the second the results applying 820 the correction to SMMR data and the third the results applying the correction to the SSM/I data.