Melt probabilities and surface temperature trends on the Greenland ice sheet using a Gaussian mixture model

. The Greenland ice sheet has experienced significant melt over the past six decades, with extreme melt events covering large areas of the ice sheet. Melt events are typically analysed using summary statistics, but the nature and characteristics of the events themselves are less frequently analysed. Our work examines melt events from a statistical perspective by modelling 19 years of Moderate Resolution Imaging Spectroradiometer (MODIS) ice surface temperature data using a Gaussian mixture model. We use a mixture model with separate model components for ice and meltwater temperatures at 1139 locations (cid:58)(cid:58)(cid:58)(cid:58) cells 5 spaced across the ice sheet. By considering the uncertainty of the ice surface temperature measurements, we use the two categories of model components to (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) define, (cid:58)(cid:58)(cid:58) for (cid:58)(cid:58)(cid:58)(cid:58) each (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) observation, (cid:58) a probability of melt for a given observation rather than using a (cid:58)(cid:58)(cid:58)(cid:58) which (cid:58)(cid:58) is (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) independent (cid:58)(cid:58) of (cid:58)(cid:58)(cid:58)(cid:58) any (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) pre-defined fixed melt threshold. This probability can then be used to estimate the expected number of melt events at a given location (cid:58)(cid:58)(cid:58) cell. Furthermore, the model can be used to estimate temperature quantiles at a given location (cid:58)(cid:58)(cid:58) cell, and analyse temperature and melt trends over time by fitting the model to subsets of time. Fitting the model to data 10 from 2001-2009 and 2010-2019 shows increases in melt probability (cid:58)(cid:58)(cid:58) and temperatures for significant portions of the ice sheet , as well as the yearly expected maximum temperatures.

2 Data and methods

MODIS IST data
We use MODIS IST data from MODIS/Terra Sea Ice Extent 5-Min L2 Swath 1km, Version 6 (MOD29) contained within a multilayer Greenland MODIS-based product (Hall and DiGirolamo, 2019) :::::::::::::: (Hall et al., 2018). MODIS records surface reflectance from 36 spectral bands of different wavelengths -including those used in IST -:::: near daily for the entire Earth. The :::: This dataset 60 spans the period 01/03/2000 to 31/12/2019 and has a spatial resolution of .78 km × .78 km. Here we discard the first 10 months of the data set, up to 01/01/2001, in order to work only with those years for which a full annual cycle is available. To reduce the computational burden of our model, we also subsample the data taking 1 in every 50 points :::: cells : in both x and y dimensions for a total of 1139 points :::: cells, roughly equally spaced across in latitude and longitude and thus covering a ::: the ::: full : range of glaciological and climatological settings across the ice sheet.
As a result of cloud masking, areas on the coast and in the north have a higher proportion of missing data than more central areas ( Figure 1). We also see that winter months have more missing data on average than summer months :::::: because ::: of ::::: cloud :::: cover, with a range of 65.1% of data available in December compared to 91.1% of data available in May. This is important to 75 bear in mind when interpreting the predictions made from the statistical models, as the IST distributions will be more heavily weighted towards warmer temperatures. This shouldn't affect our inference with regards to melt, however, as melt temperatures almost exclusively occur in the summer months which have a much lower proportion of missing data.

Modelling considerations
To create a statistical model that is parsimonious and applicable across sites with widely varying geophysical characteristics :: at 80 :: all :::: cells :::: over :: a :::: large :::: and :::::::::::::::::: geographically-varied :::::: region, we model the IST data using statistical methods that allow us to treat melting in a probabilistic manner. Exploratory data analysis shows that there is no clear point :::::: quantile : in the temperature distribution that can be attributed to :: as the onset of melting ( Figure 2). As a result, we model melting ice temperatures and non-melting ice temperatures separately and estimate the probability of melt occurring over a range of temperatures. This approach allows for some uncertainty in the observations from factors such as the precision of the dataset, which has a stated 85 uncertainty of ±1 • C). We hereby refer to temperatures associated with melting ice as "melt" temperatures and temperatures associated with non-melting ice as "ice" temperatures. A key feature of the dataset and a core modelling consideration is the soft upper limit at 0 • C. The melting point of the ice acts as a physical upper limit on ISTs, as once the ice exceeds this temperature it melts and may no longer form the surface of the ice sheet. Some sites have measurements above this limit, which arise due to meltwater sitting on top of the ice. However, 90 the ice under the water places a limit on these melt temperatures, hence the distribution of positive temperatures is truncated close to 0 • C. This soft upper limit of ISTs causes a significant peak in the distribution centred at approximately −0.5 • C, as any ISTs that would exceed 0 • C are truncated to small positive values close to 0 • C.

Model description
In order to accommodate spatial variability in the temperature distribution, we model IST using a truncated Gaussian mixture model in which components are assigned to model groups of temperatures that we assign to be either ice or melt. For n I ice components and a single melt component, let ϕ i be the weight associated with model component i such that for n c = n I + 1 total components, is the lower (upper) truncation point. Then the probability density function of ISTs x is: We set the upper and lower truncation points for the ice and melt components at values that bound each measurement type with relative certainty. For the ice components, a = −∞ as there is no hard lower limit on the temperature of ice (aside from 130 absolute zero), and b = 0 as, theoretically, ice temperature can't exceed 0 • C. This means that there is no limit on how low ice temperatures can go, but they can't exceed 0 • C. For the melt component, b = ∞ and a = −1.65, so that temperatures in the melt component can't go below −1.65 • C but are not upper truncated. We take a bound lower than zero here to account for uncertainty in the data and any potential impurities in the ice surface. −1.65 • C is the theoretical minimum temperature at which saline ice can melt (Hall et al., 2004), and thus should be a conservative estimate for this lower bound. Temperatures 135 between −1.65 • C and 0 • C can be modelled by either/both the ice and melt components as there is uncertainty as to whether they are associated with melting or non-melting ice. Appendix A). We used this method to obtain estimates of µ, σ, and ϕ for each model component at each location ::: cell.
We used Bayesian Information Criterion (BIC) to assess the most appropriate number of ice and melt components and found that three ice components and one melt component fit the data best. These components may be broadly interpreted as winter, autumn, spring, and the melt season for the three ice components and single melt component respectively.
When modelled with separate Gaussian components, the characteristics of the different modes of the data are much clearer 145 ( Figure 2). The melt component at each location ::: cell generally has a much lower variance than the ice components due to the soft upper limit of ISTs and the lower truncation point of the model, whereas the ice components have higher variances and more overlap between components. For the sites that experience melt regularly, a substantial proportion of the overall temperature distribution occurs in the overlap between true ice and true melt. A similar result is seen across sites located on or near the coasts, which further validates the decision to use a fixed melt threshold as the melt temperatures -and thereby the 150 melt process -appear to have consistent characteristics across locations :::: cells.

Defining melt
Using this model, the probability of melt occurring, which we denote by ρ(x), can be quantified as the ratio of the densities of the ice and melt components. For a given IST x, n I ice components melt component M , we have: .

155
Consequences of this definition are that for ISTs below −1.65 • C, the probability of melt is 0, for ISTs above 0 • C the probability of melt is 1, and between these values the melt probability depends on relative values of the melt and ice components' densities.  We then compare our modelled estimates to a simple threshold-based approach to defining melt, i.e. the average number of days per year with temperatures exceeding ::::: greater :::: than :: or ::::: equal :: to : −1 • C (Figure 4). The majority of the ice sheet -90.7% of locations ::: cells : from the expected melt from the model, 79.5% from a threshold of 170 −1 • C on the data -experiences some degree of melt on average each year, except for sites in the dry snow zone in the centre and north of the ice sheet :::::::::::: (Benson, 1960). Of the locations ::: cells : that experience melt, most (62.2% from the model, 57.3% from the data) sites on average see less than 2 days of melt per year, which makes up the rest of the dry snow zone and most of the percolation zone. The areas with the most melt are located around the coast and in the south and west as may be expected.
The main discrepancies between the two measures are at coastal locations ::: cells, particularly on the west and north coasts. Here, 175 the model estimates a larger amount of melt, with a maximum of 14 additional melt days at 1 specific location ::: cell : on the edge of the south east coast compared to the dataset. However, 89% of locations :::: cells have an absolute difference of less than 2 melt days, showing the broad agreement between the measures at central locations :::: cells.

Temperature quantiles
We now use the model fit to calculate quantiles of the ISTs at each location ::: cell ( Figure 5). This gives context to the overall the broader trends of high temperatures that aren't necessarily melt temperatures, as well as the 10% quantiles for temperatures that are as relatively low as the 90% quantiles are high. The estimated 10% and 90% quantiles broadly follow the same trends as elevation on the ice sheet. The 10% quantiles have a range from −53.84 • C in the centre of the ice sheet to a maximum of −15.75 • C at the south tip of the ice sheet. As would be expected, locations :::: cells at higher elevations have a lower 10% and 185 90% quantile. However, of more interest are the few (30/1139) locations ::: cells : located on the west, east and southern coasts that have a 90% quantile above 0 • C. At these locations ::: cells, we would expect at least (in some cases more than) 10% of observed temperatures to be above 0 • C and thereby melt temperatures.
We also calculate the 1-year return levels of each location ::: cell. This is the IST that is on average only exceeded once per year as estimated from each location :: cell's mixture model. The return levels range from a minimum of −7.08 • C in the centre 190 to a maximum of 7.24 • C on the west coast. Although, as previously discussed, ISTs should not be seen higher than 0 • C, these return levels reflect similar temperatures recorded by observations in the dataset and can be plausible temperatures when considering the effect that meltwater on the surface of the ice sheet has on the observations. The rarity of melt in certain central areas can be seen more clearly, as temperatures in many locations ::: cells : (519/1139) on average reach −1.65 • C less than once a year. The trends seen in the return levels also broadly agree with those seen in the quantiles, and are in reasonable agreement 195 with the elevations and distance to the coast of each location, with locations :::: cell, :::: with :::: cells at lower altitudes and closer to the coast generally experiencing more melt.

Decadal variability
To examine potential changes in melt over time, we fit mixture models at each location ::: cell for two separate decades: 2001 to 2009 and 2010 to 2019. Averaging over a decade helps to smooth some of the annual variability and thus highlight any 200 potential differences as a result of climatic change. To assess any changes in melt and high ISTs between the decades, we compare quantiles between the fitted models and the estimated melt probabilities at each location ::: cell : in each decade.

Temperature quantiles
Because some central areas of the ice sheet do not have many historical melt observations, we examine the 95% quantiles and yearly expected maximum temperatures, both of which give an indication of overall trends in high temperatures even if these 205 do not reach the level required for melt at some locations ::: cells. As previously, we take the estimated quantiles from each of the fitted decadal models for each location ::: cell.
Areas in the east show the largest increases -with the largest increase being 2.66 • C -however on the south west coast and particularly the north central area of the ice sheet there are also several locations ::: cells : that show a slight decrease in contrast to larger increases. More clearly than in the 95% quantiles, 1-year return levels at coastal locations ::: cells : do not increase as are already close to or above 0 • C. Because of the soft upper limit of the IST data, values already close to this limit can be partially constrained from further increases, so locations :::: cells that had a 1-year return level above 0 • C are less likely to show an increase than colder areas such as in the centre of the ice sheet. This makes the 1-year return level more informative for central locations ::: cells : than for coastal locations :::: cells. 220

Melt probability
We next compare the probability that each location ::: cell : experiences melt on any given day for each decade. Using the fitted models, we estimate the probability that each daily observation is a melt temperature, then take an average of all values within our defined decades. For the purposes of interpretation, we limit our discussions to the summer months (May to September, inclusive) when considering melt probabilities, due to the almost zero probability of observing melt outside of this period.

225
The two decades show very similar trends in their daily melt probabilities, particularly around the coast. However, decade 2 has more locations :::: cells with a non-zero probability of melt (1017) than decade 1 (853) -an increase of around 19.2% between the two decades -and 68.5% of locations ::: cells : saw an increased probability of melt between the decades. The location ::: cell with the single largest probability from either decade is from decade 2 (64.11,-49.93). This location ::: cell : has a probability of a melt temperature on any given summer day of 0.64 -equivalent to an expected 97.92 melt days per year.

230
Most of central Greenland has experienced minimal change in the probability of melt between the two decades ( Figure   7). This may be largely due to the probabilities being extremely small for these areas regardless of the time period chosen.
Coastal locations :::: cells : show clearer and larger cross-decadal variation in melt probability. South east and south west areas of the ice sheet were generally more likely to experience melt, in addition to some locations ::: cells : in the north east and north west areas that were less likely to experience melt, in the more recent decade. The largest increase is on the south east coast, where  Increases in ice melt in Greenland are of major concern due to the impact that it will have on sea levels etc (van den Broeke et al., 2016), however direct :::::: in-situ observations of ice melt are sparse, discontinuous ::::: spaced ::::::::: irregularly, and of coarse resolution.

255
We observe that melt is much more likely at coastal locations :::: cells and in the south of the ice sheet than in the centre, and that there is a non-trivial probability of melt occurring below −1 • C. The spatial melt trends are in keeping with previous work examining melt using surface mass balance data (van den Broeke et al., 2016) and satellite data (Mernild et al., 2011), including MODIS data (Nghiem et al., 2012). The fitted models also show a clear link between elevation and high ISTs similarly to previous studies linking temperature to elevation (Reeh, 1991), and the yearly expected maxima show the potential 260 for even central areas of the ice sheet to experience melt (Nghiem et al., 2012). Trends previously observed in the south (Mote, 2007) also appeared to have continued, as all locations :::: cells : examined south of 75.16 • N saw increases in high-temperature quantiles in the most recent decade.
Given the assumptions and intuition behind some of the modelling choices, this dataset ::: data ::: set : could alternatively be modelled using a Bayesian framework with prior distributions that reflect these assumptions. We would expect melt to have similar distributions at different locations ::: cells : even if there is less evidence of melt in some locations :::: cells than others. If this is the case, then a modelling framework could be established whereby the melt components of the model share information or 290 parameters, while the ice components are independent between locations ::: cells. This could be used to estimate melt probabilities even in locations ::: cells : where no melt temperatures have been observed, as melt components could still be estimated using information from other locations ::: cells : with more data resembling melt.
Fitting models to sub-decadal data sets would lead to insufficient data to fit the model; in particular, there would be many Let X ∼ T N (µ, σ 2 , a, b) where µ is the mean, σ is the standard deviation, and a (b) is the lower (upper) truncation point.
Furthermore, let α = a−µ σ and β = b−µ σ . Then X has probability density function: where f N and F N are the probability density function and the cumulative distribution function of a standard normal distribution respectively.

A2 Algorithm
Let (µ k , σ k , α k , β k ) denote the parameters for the kth truncated normal distribution. To initialise the algorithm, randomly sample without replacement three values of x ∈ X and set them as µ k for k = 1, 2, 3. We set µ 4 = 0 to ensure that one of the model components starts in the region of melt temperatures. Let σ k be the sample variance and the component weights ϕ k = 1/4 for k = 1, 2, 3, 4. For simplicity we refer to the truncated normal probability density function and cumulative distribution 315 function as f (x) and F (x) respectively. The EM algorithm consists of iterating between two stages, the expectation and maximisation steps, until convergence is obtained. For the expectation step, we set: whereγ ik is the estimated probability that observation i belongs to model component k.
For the maximisation step, let: We iterate between these two steps until the parameters converge to the final estimates (800 iterations was sufficient in this case). The algorithm is considered to have converged if the difference between parameters in each iteration is sufficiently small.

325
We found a difference of 10 −5 between iterations to be sufficient indication of convergence for all parameters.
Author contributions. DC and EE devised the model framework and carried out the data analysis. DC, EE, and AL worked on the interpretation of the model and results and wrote the manuscript.