Generally, snowmelt detection methods involve the problem of classifying melt signals (X) for dry snow and wet snow. The melt signal is some quantification of dry snow and wet snow. For example, in the XPGR method, the XPGR value is used as melt signal to indicate wet snow or dry snow, and the critical value (Liu et al. 2005) is used as melt signal in wavelet transformation-based edge detection approach. Let m represent the melt detection, with m=1 indicating melt and m=0 indicating non-melt: 1 m={0,if,X<T1,if,X≥T

The classification precision (determined by the threshold T) of dry snow and wet snow will directly influence the accuracy of the calculated results of ice-sheet snowmelt distribution and is therefore formulated as a determination of the optimal threshold of dry snow and wet snow. If the melt signal X of dry snow and wet snow is considered as a random variable, the optimal threshold can be determined by using statistical inference based on the statistical characteristics of the melt signal of dry snow and wet snow. Let h(X _l ), X _l =0,1,…,L−1 be the histogram of melt signal samples of dry snow and wet snow, where L stands for the number of the maximum melt signal. The histogramh(X _l ) can be considered an approximation of the actual probability density function P(Xl)=P(Xl∣ωdry)P(ωdry)+P(Xl∣ωwet)P(ωwet) of the mixture population describing the wet snow and dry snow. In the Kittler-Illingworth (KI) threshold selection criterion (Kittler & Illingworth 1986), the selection of an appropriate decision threshold T∈{0,1,…,L−1} is based on the optimization of the criterion function: 2 J(T)=∑Xl=0L-1h(Xl)c(Xl,T)

The function c(X _l ,T) over the histogram h(X _l ) is a cost function, which measures the cost of classifying snow by comparing their melt signal with the threshold T. The optimal threshold that minimizes the classification error is the one that minimizes the following cost function: 3 T*=argminJ(T)T=0,1,...,L-1

Determining cost function depends on the probability density function (pdf) of melt signal of dry snow and wet snow. Because different detection methods use different melt signals, the pdfs of melt signals of different methods is likely to be different. For this reason, we considered a model that should be adaptable to the pdfs of different melt signals; that is, it should be capable of spanning a large variety of statistical behaviours. Moreover, for easy implementation during regular operation, it should not require the estimation of an excessively large number of parameters. Among the possible models, the GG distribution is a particularly attractive candidate. The analytical expression of the GG distribution considered in our approach for modelling the two class-conditional pdfs is given by (Sharifi & Leon-Garcia 1995; Niehsen 1999): 4 p(Xl∣ωi)=aie-[bi∣Xi-mi∣]βi,i=wet,dry,

where the positive constants a _i and b _i are given by: 5 ai=biβi2Γ(1βi),bi=1σiΓ(3βi)Γ(1βi)

The terms m _i , σi2 and β _i are the mean, the variance and the shape parameters of the distribution, respectively, and Γ(·) is the well-known Gamma function. We can see from Eqn. 5 that the GG model requires the estimation of only three parameters, which are m _i , σi2 and β _i . The shape parameter β _i , in the GG model, tunes the decay rate of the density function. When β _i =2, the GG density function is actually the Gaussian density function, and β _i =1, it is the Laplacian density function. Moreover, in the limit case β _i →0 or β _i →∞, it approaches an impulsive function or uniform distribution, respectively. Thus, the GG model can adapt to a large class of statistical distributions.

Under the GG distribution assumption, it can be proved that the cost function is optimized as follows (see the Supplementary File): 6 J(T)=∑Xl=0Th(Xl)[bdry(T)∣Xl-mdry(T)∣]βdry(T)+∑Xl=T+1L-1h(Xl)[bwet(T)∣Xl-mwet(T)∣]βwet(T)+H(Ω,T)-[Pdry(T)lnadry(T)+Pwet(T)lnawet(T)],

where: 7 (Pdry(T)=∑Xl=0Th(Xl)mdry(T)=1Pdry(T)∑Xl=0TXlh(Xl)σdry2(T)=1Pdry(T)∑Xl=0T[Xl-mdry(T)]2h(Xl)Pwet(T)=1-Pdry(T),mwet(T)=1Pwet(T)∑Xl=T+1L-1Xlh(Xl)σwet2(T)=1Pwet(T)∑Xl=T+1L-1[Xl-mwet(T)]2h(Xl)

P _dry (T) and P _wet (T) are the prior probabilities; m _dry (T) and m _wet (T) are the means; σdry2(T) and σwet2(T) are the variances; β _dry (T) and β _wet (T) are the shape parameters of the distribution. These parameters are all associated with the dry snow and wet snow samples, respectively, for a given value of the threshold T. H(Ω,T) stands for the entropy associated with the binary set of classes Ω={ω _dry , ω _wet1}. For the steps to estimate a _i , b _i , β _i , please refer to Sharifi & Leon-Garcia (1995).

As shown above, the advantages of the GG model method are: it requires less manual input parameters (the parameters are computed automatically); the classification result is unique (namely, the threshold is unique); and it is capable of spanning a large variety of statistical behaviours.

The melt signal adaptation method can be applicable to many different melt signals in many different snowmelt detection methods, for example the XPGR value in XPGR method, the normalized gradient ratio (GR) value in GR method GR=(Tb19H-Tb37H)/(Tb19H+Tb37H) , the critical value in the ice-sheet snowmelt detection method based on wavelet transformation, the melt signal X in α-based T _b melt detection method X=(Tb-Tbdry)/(Twet+Tbdry) , and so on.

Here, we only apply the method to the existing XPGR and wavelet transformation-based ice-sheet snowmelt detection methods. Further work will focus on the GR method, the α-based method and other methods.

The XPGR method proposed by Abdalati & Steffen (1995) mainly utilized the different response of two frequencies (19 GHz and 37 GHz) and two polarizations (horizontal and vertical) on ice sheets to detect snowmelt. The advantage of this method is that the response differences of frequency and polarization to emissivity and moisture contents were used which clearly showed snowmelt information and that the model was stable. The XPGR equation is: 8 XPGR=Tb19H-Tb37VTb19H+Tb37V,

where Tb19H is the microwave brightness temperature at the 19-GHz horizontally (19H) polarized channel and Tb37V is the microwave brightness temperature at the 37-GHz vertically (37V) polarized channel.

Fig. 2 shows the XPGR long time series change derived from SMMR and SSM/I data sets from 1978 to 2010 at the Wilkins Ice Shelf area in Antarctica, which experienced melting every year. By sampling XPGR values both in the summer season (November, January and February) and other seasons, we analysed the statistical characteristics of the XPGR values of wet and dry pixels. Fig. 3 shows the histogram of XPGR values. The distribution characteristics of XPGR values of dry snow generally have a low mean value whereas for wet snow the mean values are generally higher. As we can see in Fig. 3, the XPGR values are divided into two dominant categories with a valley between them: to the left are low values that indicate the distribution characteristics of dry snow; to the right are higher values representing wet snow. To model a bimodal histogram of XPGR values, we can consider it as a mixture of two GG models. Then, we computationally obtain an optimal threshold value for separating wet pixels from dry pixels by applying the melt signal adaptation method. The XPGR optimal threshold calculated from the method is −0.0239 (see Figs. 3 and 4). Fig. 4 shows that the optimal XPGR value corresponds to sharply changed regions, namely, it corresponds to time points when snow changed from dry to wet or from wet to frozen. We propose a modified XPGR detection method, based on this. Fig. 5 is a flowchart of the modified XPGR detection method.

Ice-sheet snowmelt detection using a wavelet transformation-based method is based on the strong and significant edges in the brightness temperature (Tb) time series curve that signify snow melting and refreezing events (Liu et al. 2005). The basic steps are as follows. First, the time series brightness temperature values of each pixel are decomposed into multi-scale components through a wavelet transform. Second, the local extrema of the wavelet transform modulus (modulus extrema), which indicate the strength of the edges caused by sharp transitions and the time when the sharp variations occur, are tracked and analysed across scales. Third, for each pixel, a critical value as a melt signal to differentiate the sharp transition induced by melting and refreezing processes from other transitions caused by non-melt processes is determined by variance analysis. Then, based on critical values of samples of wet pixels and dry pixels, a bimodal Gaussian (BG) model fitting technique is used to statistically determine an optimal edge threshold to classify wet snow and dry snow. Finally, based on the principle of spatial autocorrelation, a spatial neighbourhood operator is used to detect and correct possible errors brought about by the strong noise of data pixels. This method can determine not only whether an area experienced melt, but also when an area experienced melt by detecting and tracking strong and significant edges in the brightness temperature time series curve.

Liu et al. (2005) used the BG model fitting to determine the optimal threshold. As stated in an earlier section in this article, the GG model is also adaptable to the BG model of melt signal for the wavelet-based method. It is therefore possible to apply the melt signal adaptation method to the existing ice-sheet snowmelt detection method based on wavelet transformation. Through selection of the samples of wet and dry snow in Antarctica based on the metrics of Liu et al. (2005), sampling both wet pixels that experienced melting and dry pixels that remained dry, we use the GG models to determine the thresholds of dry snow and wet snow based on the Antarctic SMMR 18-GHz horizontal polarization channel data and the SSM/I and SSMIS 19-GHz horizontal channel data from 1978 to 2010.

Fig. 6 shows the optimal threshold—10.30—derived from the melt signal adaptation method. Using this, we propose a modified ice-sheet snowmelt detection method based on wavelet transformation (Fig. 7).

The difference between the modified XPGR detection method and the XPGR method lies in the method in which the threshold is determined. The XPGR method uses in situ measurements to determine the threshold. The modified XPGR method uses the GG model to automatically seek an optimal threshold. Therefore, the modified XPGR method is solely dependent on the XPGR value. To analyse and compare the results of the method proposed by Abdalati & Steffen (1997), SSM/I data from 1988 to 2000 at the same experimental area (Swiss Camp, Greenland, 69°34′N, 49°17′W, 1150 m a.s.l.) is used.

Fig. 8 shows the XPGR threshold determination respectively based on Abdalati & Steffen's (1997) method and the modified method in the Greenland study area. It is evident that the difference between the XPGR threshold obtained by Abdalati & Steffen (1997) and the unsupervised XPGR threshold is very small.

Fig. 9 shows ice-sheet snowmelt distribution drawn from the two methods on 2 August 1996. The white indicates melt area, and grey indicates the non-melt area or bare rock area. The results of the two methods are almost the same, with only small differences seen. Near the AWS at Swiss Camp, the modified XPGR method did not detect melting on snow surface. Meanwhile, the maximum temperature measured by this station was −2.59°C (temperature measured at 1 m above the snow surface) on 2 August 1996, which indicates that it is highly possible that no melting occurred around this area. This result supports the accuracy of the proposed method.

Analysis and comparison of the modified and original wavelet transformation detection method results are based on the Antarctic SSM/I data from 1978 to 2010. Fig. 7 is a flowchart of the ice-sheet snowmelt detection method based on wavelet transformation (Liu et al. 2005). The difference between the modified wavelet transformation detection method and the original method is that the original uses the BG model to determine optimal threshold, whereas the modified wavelet transformation detection method uses the GG model-based method to seek an optimal threshold. In Fig. 6, we can see that the optimal threshold value for the optimal threshold derived from the BG method l is 10.28 whereas that for the modified method is 10.30.

Fig. 10 shows the ice-sheet snowmelt distributions obtained from the modified wavelet transformation detection method and the original method during 2007–08. The two results are similar, but the result obtained from the BG distribution model shows that there is some inland melting in Antarctica, which is unlikely. In contrast, the GG model automatic threshold segmentation indicates no such melting. Fig. 11 shows melt onset, melt end and melt duration distribution obtained from the modified wavelet transformation detection method in Antarctica during 2007–08.

This comparison shows that the optimal threshold based on the GG model is almost the same as those based on the BG model. However, the melt signal adaptation method increases threshold determination efficiency, has higher accuracy and is independent of estimation and initial value.

In addition to comparing the results of the modified and original methods, we have performed quantitative validations on the effectiveness and precision of the modified ice-sheet snowmelt detection method.

We use in situ data to validate the daily snowmelt detection results of each method. In general, surface melting occurrence spatially corresponds to the spatial pattern of temperature. Therefore, we validate the snowmelt detection results derived from two methods with the corresponding surface air temperature acquired by AWSs. An AWS cannot cover the same area observed by large-scale remote sensing and can therefore not completely reflect large-scale climate. However, to minimize the influence of the environment, we choose AWSs in those areas where the spatial gradients are relatively simple because the temperature data collected on these stations can reflect average temperature of large areas to some degree. Based on the Radarsat Antarctic Mapping Project Digital Elevation Model, we use the 1-km data to calculate the minimum height, maximum height and standard deviation of height in the 625 km² area around the AWS. We use those parameters as a standard of evaluation of spatial complexity. We choose the AWS in the area where the standard deviation of height and difference between the extreme values of height and elevation of AWS are smallest.

Because of the limited temperature records and strong fluctuations in surface melt on the Ross Ice Shelf, on the Ronne-Filchner Ice Shelf and in Marie Byrd Land AWSs were not chosen in those regions. The AWSs on Butler Island, at Bonaparte Point, on Larsen Shelf and at Amery G3 were chosen (Table 2).

The daily mean air temperature data together with the corresponding daily detection results of each modified snowmelt detection method were used to quantitatively validate the snowmelt detection results of each method. For processing the AWS data, we chose the daily mean air temperature data from the summer months (November–March) in Antarctica. The total number of days N (i.e., the total number of daily mean air temperature data) chosen in each AWS can be seen in Table 2. The daily mean air temperature data are made up of the mean of three daily maximum air temperatures; a daily mean air temperature data above 0°C is regarded as snowmelt occurrence. A melt day is the day with daily mean air temperature data above 0°C and a non-melt day is the day whose daily mean air temperature data is below 0°C. Once the corresponding daily detection results have been obtained by using the snowmelt detection method, the following defined quantities (Celik 2010) are computed for comparing the daily detection results against the corresponding daily mean air temperature data:

True positives (TP): the number of melt days that were correctly detected, and P _TP =TP/N.

Based on these four quantities, three assessment indices are proposed to evaluate the snowmelt detection result of the proposed method. The three assessment indices are as follows: Correct Detection Rate (CDR), Priori True Positives Rate (Priori TPR) and Posterior True Positives Rate (Posterior TPR). CDR is the ratio between the number of correctly detection results and the total number of detection results (or the total number of temperature data); Priori TPR is the ratio between the number of correctly detected “melt” results and the number of temperature data above 0°C; Posterior TPR is the ratio between the number of correctly detected “melt” results and the number of detected “melt” results: 9 CDR=(TP+TN)/N 10 PorioriTPR=TP/(TP+FN) 11 PosteriorTPR=TP/(TP+FP)

The statistical analysis results are shown in Table 3. For CDR, two methods have relatively high probabilities in several areas around the AWSs. The percentages of CDR exceed 65%, except for the XPGR methods in the area around the Bonaparte Point AWS. The CDR indicated that the improvements on two snowmelt detections are effective.

However, we were primarily concerned with Priori TPR and Posterior TPR. Although the comparative statistics with in situ data were collected at the Bonaparte Point AWS and the Larsen Ice Shelf AWS, the result values of both methods show relatively high probabilities, except the XPGR methods in the area around the Bonaparte Point AWS. Moreover, the Posterior TPR of both methods are over 90%. For the Amery G3 AWS, the Priori TPR of both methods are over 65%; however, the Posterior TPR show relatively low probabilities. For the Butler Island AWS, both methods have low Priori TPR and low Posterior TPR. This is because the Butler Island AWS is very close to the seaside, so it does not reflect the average temperature of a large adjacent inland area.

Using temperature data to validate the snowmelt detection results of the modified XPGR detection method and the modified wavelet transformation-based method shows that they have relatively high precision.