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

. Surface melting is a major component of the Greenland ice sheet (GrIS) surface mass balance, affecting sea level 10 rise through direct runoff and the modulation on 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 in view ofdue 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 passive microwave dataset (37 GHz, 15 horizontal polarization) 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 GrIS at an enhanced spatial resolution of 3.125 kKm. 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 20 a fixed threshold of 245K is suitable in describingcapable of detecting persistent melt. Our results indicate that, during the reference period 1979 – 2019 (1988 – 2019), surface melting over the GrIS 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 GrIS surface extent per decade, according to the MEMLS of dataset with respect to the one at 25 km and a better fitting spherical to the empirical semi-variogram in case of the 3.125 km and MAR melt duration maps. Our analysis suggests that the resolution product sensitive to local scale processes, 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 630 analysis of spatial scales to the atmospheric drivers of surface melting, such as incoming solar radiation, surface temperature 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 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 635 to investigate the evolution of the surface meltwater networks and the application of similar tools to other regions, such Arctic Archipelago, the Himalayan Plateau and the Antarctic Peninsula, where the enhancement in spatial resolution can be fully exploited.

rise through direct runoff and the modulation on 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 in view ofdue 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 passive microwave dataset (37 GHz,15 horizontal polarization) 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 GrIS at an enhanced spatial resolution of 3.125 kKm. 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 20 a fixed threshold of 245K is suitable in describingcapable of detecting persistent melt. Our results indicate that, during the reference period 1979 -2019 (1988 -2019), surface melting over the GrIS 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 GrIS surface extent per decade, according to the MEMLS algorithm. Furthermore, the melting season has started up to ~4 (2.5) days earlier and ended ~7 (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 25 through a semi-variogram approach. We found that the enhanced product is more sensitive to local scale processes, hence confirming the potential interest of this new enhanced product for studying 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. 30

Introduction
The Greenland ice sheet (GrIS) 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 m 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 (SLR). 35 According to 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 GrIS 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 GrIS together with the Antarctic ice sheet is crucial to assess the impact of global warming on SLR and the global water balance (Kargel et al., 2005;40 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 (SMB), 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., , 2017van den Broe ke 45 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;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;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 4 microwave data and propose an update on a recently proposed algorithm in which meltwater is detected when T b 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) and with the outputs of 100 the regional climate model Modèle Atmosphérique Régional (MAR; Fettweis et al., 2017). 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 (M) and evening (E) 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;115 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 the measurements within a given time falling into a specific pixel are averaged. In the reconstruction algorithm used for the enhanced Tbs, the so-called effective measurement response function (MRF), 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  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 interpolatin g and smoothing data to match the resolution between different channels (Robinson et al., 1992), and improving the spatial resolution of surface brightness temperature 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 125 the case of both the coarse and enhanced resolution products over 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 (69° 34' 06" N, 49° 18' 57" W as shown in Figure 1). From the figure we observe that the two 130 timeseries are highly consistent with each other, with a mean difference of 0.895 K and standard deviation of 4.89 K, 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 processs. 135 In order toTo compare our results with precedent PMW surface melting products, we perform our calculations also to the widely used PMW dataset by Mote (2014). The dataset consists in a daily melt maps product obtained from 37 GHz brightness temperatures 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 140 generated using a dynamically changing threshold (Mote, 2007) obtained simulating the brightness temperature 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 (https://nsidc.org/data/NSIDC-0533/versions/1). 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 (69° 34' 06" N, 49° 18' 57" W as shown in Figure 1). From the figure we observe that the two 145 timeseries are highly consistent with each other, with a mean difference of 0.895 K and standard deviation of 4.89 K, 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. 150

Greenland air/surface 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.
Specifically, we use data recorded by stations of the Greenland Climate Network (GC-Net; Steffen et al., 1996). The AWSs 155 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/.
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 (SISVAT) (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;Fettweis et al., 2017;Alexander et al., 2014). 165 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 SISVAT.
In this study, we use the output from MAR version v3.11.2 characterized by an enhanced computational efficiency and 170 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 meltwater extent obtained from PMW data, we average the LWC simulated by MAR along the vertical profile, following Fettweis et al. (2007). 175

Melt detection algorithms
Generally speaking, melt detection algorithms can be divided into threshold-based (T-B) and edge-detection (E-D) 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 algorithm. For example, Steffen et al. (1993) used the normalized gradient ratio GR=(Tb19H -Tb37H)/(Tb37H +Tb19H) to detect 180 wet pixels with a threshold value computed based on in-situ measurements. This method was later updated by  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= αTwinter+(1-α) Twet where Twinter is the average of winter Tb, and Twet fixed as 273 K and Tc indicates the threshold value (we keep the same 185 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 (Tedesco et al., 2007).
Several other studies have been detecting melting when Tb values exceed the mean winter value Twinter plus an additional value ∆T (M+∆T approaches, where M represents the January-February mean brightness temperature) associated with the insurgence of liquid water within the 190 snowpack: Torinesi et al. (2003) proposed a value of ΔT=Nσ with Twinter 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 ΔT=30 K. Tedesco (2009) proposed an alternative approach based on the outputs of the 195 Microwave Emission Model of Layered Snowpack (MEMLS) electromagnetic model (Weisman and Matzler, 1999). In this case, an ensemble of outputs is generated by MEMLS by varying the inputs (e.g., correlation length, LWC, density, etc.).
These outputs are, then, used to build a linear regression model for the ∆T 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 200 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 presen ce of liquid water in the case of a snowpack with relatively coarse grains will correspond to a lower value of LWC than an increase 205 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 two LWC values of 0.1 % and 0.2 %. The coefficients are φ = -0.2 and ω = 58 K (R 2 =0.91) in the case of LWC=0.1% and φ = -0.52 and ω = 128 K (R 2 =0.92) in case of LWC=0.2%. . The Tb threshold value computed in this case can, therefore, be written as follows: T c = T winter + φT winter + ω = (1 + φ)T winter + ω = γT winter + ω, (3) 210 where (γ, ω) assume the values of 0.80 and 58 K in the case of LWC=0.1% and 0.48 and 128 K in the case of LWC=0.2%.
Here, we focus five approaches: the M+∆T approach choosing ∆T equal to 30K and,, to test the sensitivity to Zwally and Fiegles (1994), 35K and 40K as sensitivity to Zwally and Fiegles (1994) (M+30, M+35 and M+40 from here on), the algorithm based on MEMLS in case of LWC=0.2% (MEMLS from here on) and the 245 K fixed threshold (245K from here on). We selected M+∆T and MEMLS (LWC=0.2%) due to their higher accuracy in detecting both sporadic and persistent melting with 215 respect to the other approaches presented above (i.e. Torinesi et al. (2003), Ashcraft and Long (2006) and MEMLS in case of LWC=0.1%) 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 (0.2 % LWC). As we explain below, this choice was driven by the performance of the different considered algorithms. Moreover, we found that the fixed-threshold algorithm is more sensitive 220 to persistent melting where the MEMLS-based one can detect sporadic melting. This allows us to analyze both melting conditions (sporadic vs. persistent) and analyze them within the long-term, large spatial scales that the PMW data can provide.
In view of the novel nature of the PMW dataset, 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 225 angle, altitude and Local-Time-Of-Day (LTOD) 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 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 ~ 37GHz) for both horizontal and vertical polarizations. Steffen et al. (1993) proposed an approach focusing over Greenland for the K-band;  derived relations between SSM/I 230 observations for the F08 and F11 platforms over Antarctica and Greenland for 19.35GHz, 22.2GHz 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 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 235 et al. (1998), we perform the intercalibration using only data collected over the Greenland ice sheet. The overlapping periods for the different sensors are the following: SMMR and SSM/I-F08 overlap between 07/09/9 July 1987 and 08/20/20 August 1987 for a total of 22 days (one every two days as sensed by SMMR sensor); F08 and F11 overlap between 12/03/ 3 December 1991 and 12/18/18 December 1991 for a total of 16 days, F11 and F13 overlap between 05/03/3 May1995 and 09/30/30 September 1995 for a total of 76 days; and F13 and F17 overlap for the period 03/01/1 March 2008 --12/10/10 December 240 2008 for a total of 71 days. 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 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 245 and intercept values from the regression of daily data. Considering n days, for every i-th day we first compute mi, qi and the 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 method, we consider all values for all days when data from both platforms are available and then evaluate m and q through a Formattato: Titolo 2 linear regression fitting procedure based on least-square fitting. As an example, 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 SMMR and SSM/I-F08 sensors, reporting values of m, q and R 2 . Using the estimated values of m and q, we then correct the values for one of the sensors by 255 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 brightness temperature values and evaluating the matching 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 brightness temperature 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. 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 Tbs obtained as: 265 where Di is the absolute difference between the two histograms A and B for the i-th value of brightness temperature. Then, we sum the differences over the total number of pixels and compute the relative variation as follows: where Doriginal and Dcorrected are, respectively, the summations of Di before and after the calibration. The relative variation d can 270 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
Variogram analysis is generally adopted in geostatistical analyses to evaluate autocorrelation of spatial data 275 (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 We point out that the overlap between SMMR and SSM/I-F08 data occurs in the months of July and 280 August. During these months, the differences between acquisition times might lead to biases and errors associated with snow conditions (e.g., wet vs. dry). 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 Tbs (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 285 and SSMI/S, R2 values are higher, mostly around 0.98. In Figure 4 we also show 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. 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 290 larger than 1 K, consistent with previous results obtained by  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 295 worsened the agreement between the two sets of measurements.
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 300 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 305 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 May1995 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.
As an example, inIn Figure 3 we show the scatter plots of the data used for the linear regression for Greenland for 310 both evening and morning passes for the 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 315 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 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  (Table 3). On the other hand, in the case of the SSM/I and SSMI/S sensors, the average difference 325 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  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 330 -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 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 335 situ air/surface temperature daily averaged from AWS as an index of surface melting (Braithwaite and Oelsen, 1989) daily averaged from AWS and with the 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. 340

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 because of radiative ha formattato: Pedice ha formattato: Pedice ha formattato: Apice 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 345 PMW-based algorithm by defining commission and omission errors. Commission error occurs when melting is detected by 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 5  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 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-355 truth reference (e.g. from 0.70% for MEMLS and 0.09% for 245K to 5.63% for M+30, in case of Tair=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 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 360 terms of C+O when applying the M+40 and MEMLS algorithms. The sensitivity to the air/surface 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.
In order to better understand the sources of the relatively high values of the commission errors at some locations, we show in Figure 5  plotted as horizontal lines (black) as well as the 0°C air/surface 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) shows the sensitivity of brightness temperature to physical temperature and its seasonal variations. In this case, the surface/air temperature remains below 0°C throughout year and the Tb signal does not exceed any threshold value 370 (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/surface 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) 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 375 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 being the most sensitiveproviding the lowest threshold and the 245K fixed threshold being the most conservative. The computed rough estimation of the average emissivity for the period 17 June -17 July (as Tb divided by the recorded surface temperature) also suggests that melting is not occurring in the considered 380 period in Humboldt case, presenting an average emissivity even lower than in Swiss CampSummit 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+30K and M+35K algorithms suggest melting up to high elevations, within the dry snow zone, where it likely did not occur. The M+40K and MEMLS algorithms show similar results, while the 245K fixed-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 385 could produce errors if there is a large seasonal range in brightness temperatures 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 brightness temperature 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 amount of LWC detected in the snowpack. 390

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 assume melting is occurring to 0.2% for both depths. We selected two different depths for our analysis so we could study two 395 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 approaches (i.e., 245K) perform better when considering the LWC averaged on the first 1m 5 cm of snowpack. In fact, the 400 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 5cm1 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 M+30) in North-East Greenland (e.g., Humboldt and GITS stations). This confirms 405 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 between the PMW brightness temperature and MAR outputs. Figure 7a and b show, respectively, the timeseries of LWC 410 averaged over the first 5cm (LWC5cm) and 1m (LWC1m) obtained from MAR at the Swiss Camp site. In Figure 7c we report the Tb timeseries and the daily average surface 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% on this day while MAR5cm a LWC of 0.6%. This suggests that in some cases (before the main melt 415 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 (Tair>-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 exceed any threshold. This corresponds to a rainfall event (simulated by MAR) suggesting that the sensitivity to liquid clouds 420 of the 37 GHz channel could mask some melt events. Moreover, at the end of the melting season the brightness temperature 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 conservative among the approaches we tested, providing the lowest (highest) commission (omission) error but being unable to 425 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 GrIS), 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 to low levels of LWC. Considering the overall average error (C+O Mean in Table 5), the MEMLS algorithm shows the best 430 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 (ME) with that estimated from MAR outputs. In coefficient (Nash and Sutcliffe, 1970), described in Section 3, for a quantitative analysis. We remind hereHere, we remind the reader that NSE can assume values in the interval (-∞,1]. A perfect match between the two timeseries is achieved when the 440 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 to the different pixel size). We compared ha formattato: Pedice 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 ME 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 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]  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, indicating a weaker underestimation of the ME. This can be a consequence of the coarser resolution, lacking in capturing melting areas at the edges 470 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.
The comparison with 25 km historical surface melting dataset shows the same underestimation issue of 245K. Even showing 475 lower errors than in case of 245K, we did not found 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 (mean melt duration, MMD, averaged over the GrIS area) and at a pixel by pixel scale. We also study the maximum melting 480 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 GrIS 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 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 485 computed according to the coarse resolution dataset with reference to the time period of data availability ( 1979 -2012). suggest that the MMS has been increasing by 0.69% y -1 in the case of MEMLS and 0.94% y -1 in the case of the 245K algorithm for the period 1979-2019 (percentage with respect to the whole GrIS 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 are on average 3.4% lower than the ones computed using the 25 km dataset are equal to 1.31% y -1 .for the period 1979 505 -2012 and 0.91% y -1 for the period 1988 -2012, larger than 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 ha formattato: Apice ha formattato: Apice ha formattato: Apice ha formattato: Apice case of MEMLS. The increase of the trends passing from the 25 km to the 3.125 km resolution dataset is on average equal to 510 10.8%As in case 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 MEMLS (1.428 x10 6 km 2 d y -1 and 1.325 x10 6 km 2 d y -1 , respectively) and 245K1979-2012. (0.999 x10 6 km 2 d y -1 ) but larger than 245K1979-2012 (1.019 x10 6 km 2 d y -1 ). Lastly, we report in Figure 9d the melt onset date and melt end date averaged spatially over the pixels with 95% significant trends in Figure 11. We found that average In Figure 10 we show the trends of melt duration (MD), melt onset date (MOD) and melt end date (MED) on a pixelby-pixel basis for the complete time series . We found that the trend in melt duration exhibits the highest statistical significance (in terms of number of statistically significant pixels), being the most stable and reliable trend among the pixel-530 by-pixel parameters analyzed. We found mostly positive trends in melt duration 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 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 535 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 245 K 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 trends computed at ice sheet scale since we computed the average over the statistically significant pixels only.

Spatial information content 540
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 melt duration estimated from the two products when using either the ha formattato: Pedice ha formattato: Pedice ha formattato: Pedice 245K or the MEMLS algorithms. 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 545 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 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 550 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 melt duration estimated from the MEMLS and 245K algorithms for the months of May through August of 2012 when using either the enhanced or the coarse resolution products. We plan to extend this analysis to other variables and periods as a part of our future plans. We also performed the same analysis applied to melt duration estimated using LWC modelled with MAR, according to the same rationale described in the previous sections. 555 Here, we compute the melt duration 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) semi-variograms for Greenland melt duration computed applying the MEMLS and 245K algorithms to both 25 km 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 in case of melt 560 duration obtained according to LWC1m and LWC5cm approaches. 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 magnitude with the enhanced resolution case only in May. When computing the spherical fitting of the empirical variograms of melt duration from MAR, we found for each case considered similar R 2 values (between 0.118 and 0.484). The values for the range 565 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 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 results is lower in case of 570 LWC5cm (187.70 km) than in case of LWC1m (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 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 till until July and decreases in August. 575 Moreover, a proper fitting of the experimental variograms is achieved for all cases for the enhanced resolution PMW and the MAR derived melt duration. 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 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, 580 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 meltwat er 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 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 585 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 over the past ~ 40 years. Further analysis will help to shed light on the processes responsible for the recent acceleration of surface melting. 590

Conclusions and future work
We applied threshold-based melt detection algorithms to the 3.125 km resolution 37 GHz horizontal polarization passive microwave brightness temperatures to assess the skills of the enhanced PMW product to detect surface melting over the 1979-2019 period over the 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 595 detected from PMW data with that estimated from AWS surface/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 (245 K) 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 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 600 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). Following Fettweis et al., (2007), we set up to 0.2 % the minimum value for the LWC for a pixel in MAR to be flagged as "melting". We obtained good matching (i.e., NSE>0.4 or, at least, positive) in most of the years from 19791992-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 matching in the period 1979-1992, possibly due to 605 differences in sensor characteristics. In the case of melt extent retrieved by 245K, we found a strong underestimation of melt extent (largelyy negative values of NSE) from 1979 to 1987 likely because of the lower values of "wet" brightness temperature in case of SMMR data, slightly improving from 1993 to 2019 but still negative. Accordingly, the results obtained applying comparing with the coarse resolution dataset, we found that the melt extent timeseries derived from the enhanced resolution 610 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 melt onset date, melt end date, the duration of the melting season and the melting areal extent for the period 1979 -2019. According to MEMLS algorithm, we found that Finally, we explored the information content of the enhanced resolution dataset with respect to the one at 25 km and 625 the MAR outputs through a semi-variogram approach. The results obtained showed a better fitting of the modelled spherical function to the empirical semi-variogram in case of the 3.125 km and MAR melt duration maps. 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 pav es the way for using this dataset in statistically downscaling model outputs. In this regard, as a future work, we plan to extend the 630 analysis of spatial scales to the atmospheric drivers of surface melting, such as incoming solar radiation, surface temperature 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 635 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 in spatial resolution can be fully exploited.       Table 4: Slope (m) and intercept (q) obtained from the linear regression analysis between SMMR and SSM/I-the selected couples of satellites F08 enhanced PMW brightness temperatures 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).