Model architecture

GlacierLake was written in MATLAB R2017b using the equations presented below, with source code and equations adapted from Benedek (Reference Benedek2014) (who adapted source code from Lüthje and others, Reference Lüthje, Feltham, Taylor and Worster2006), Buzzard and others (Reference Buzzard, Feltham and Flocco2018), Saloranta and Andersen (Reference Saloranta and Andersen2007) and Essery (Reference Essery2015). GlacierLake represents a single point in x, y space modelled as a $1\, {\rm m}^{2}$ column as any scaling factors cancel. The model is divided into five stages (Fig. 1 and Table 1). Movement through GlacierLake, and processes at each stage, are detailed in Table 2 and occur once a given threshold has been reached. The model is also allowed to move backwards through stages (for example stage 5 to 3 in Fig. 1) to mimic year-on-year evolution. The model was developed using weather data from the Program for Monitoring of the GrIS (PROMICE, vanAs and others, Reference van As and Fausto2010). The AWS UPE-U, at an elevation of 980 m a.s.l. and $72.89^{\circ }{\rm N}$, $53.58^{\circ }{\rm W}$, was used for development, as the record for this station was the most comprehensive of any PROMICE station in the ablation zone of the GrIS.

Fig. 1. Schematic diagram of model stages with explanations in Table 1. Width of each stage schematically represents expected residence time. Segment shapes represent schematic evolution of phase throughout stage. The initial section of stage 4 is snow free to show that lid can function with or without snow cover. Height of each column indicates the overall depth of the model domain.

Table 1. Description of model stages

Table 2. Progression through model stages with brief descriptions of processes occurring at each stage

Trackers monitor events such as initialisation of snow layer to ensure repetition does not occur. Italicised text indicates setup and sub-functions, normal text indicates conditional statements.

GlacierLake comprises two modules. The main module contains cells of constant size throughout the model run (default 0.1 m for upper 15 m and, 1 m below). These cells are in one of three states: ice, water or ‘slush’, where ice and water coexist in a ratio dependent upon the total enthalpy of the cell. The density of water was used for slush, ice and water to avoid a change in cell depth with state as the difference is small (11%) and therefore assumed to be insignificant (Lüthje and others, Reference Lüthje, Feltham, Taylor and Worster2006; Benedek, Reference Benedek2014; Buzzard and others, Reference Buzzard, Feltham and Flocco2018). A second module deals with snow on the surface when present and incorporates the snow physics of Essery (Reference Essery2015) (see the Snow cover section and Fig. 2) to avoid computationally expensive resizing at each timestep as with Buzzard and others (Reference Buzzard, Feltham and Flocco2018). If the domain depth varies due to meltwater input, a temporary lake insert comprised cells of standard size is used within the main cell column (Fig. 2). A bucket filling method is used for each new cell to keep cell sizes uniform, so a new cell is only added when additional meltwater inflow exceeds 10 cm (as default). If the meltwater inflow input ceases, the lake insert is combined with the main column in a one-off operation to reduce complexity in subsequent stages.

Fig. 2. Construction of separate snow module (left) and the lake insert contained within the main module (right) during meltwater inflow. The lake insert size is increased discretely whereas the snow module size is increased continuously. The lake insert is combined with the main grid after meltwater inflow is complete. The grid size for snow/ice/slush cells remains constant in the upper section of the model but can increase in deeper cells (not shown in this figure).

Surface energy balance

Surface energy flux in ${\rm W}\, {\rm m}^{-2}$ is calculated after Buzzard and others (Reference Buzzard, Feltham and Flocco2018) who follow Ebert and Curry (Reference Ebert and Curry1993), as

(1)$$F_{{\rm total}}=\varepsilon F_{{LW}_{{\rm in}}}+\lpar 1-\alpha\rpar F_{{SW}}- \varepsilon \sigma T^{4}+F_{{\rm sens}}+F_{{\rm lat}} {\comma\; }$$

where ɛ is the surface emissivity, $F_{{LW}_{{\rm in}}}$ is incoming longwave radiation (${\rm W}\, {\rm m}^{-2}$), α is albedo, $F_{{SW}}$ is incoming shortwave radiation (${\rm W}\, {\rm m}^{-2}$, separated into water-penetrating and surface radiation when a lake is present using eqn (12)), σ is the Stefan–Boltzmann constant (${\rm W}\, {\rm m}^{-2}\, {\rm K}^{-4}$), T is temperature (K), $F_{{\rm sen}}$ is sensible heat flux (${\rm W}\, {\rm m}^{-2}$) and $F_{{\rm lat}}$ is latent heat flux (${\rm W}\, {\rm m}^{-2}$). Sensible heat flux is calculated as

(2)$$F_{{\rm sens}} = \rho_{a}C_{\,p}^{{\rm air}}C_{T}v\lpar T_{a}-T_{s}\rpar {\comma\; }$$

where $\rho _{a}$ is the density of dry air ($1.275\, {\rm kg}\, {\rm m}^{-3}$), $C_{p}^{{\rm air}}$ is the specific heat capacity of dry air (${\rm J}\, {\rm kg}^{-1}\, {\rm K}^{-1}$), $C_{T}$ is a function of atmospheric stability described in Ebert and Curry (Reference Ebert and Curry1993) (eqns (4) and (5)), v is the wind speed (${\rm m}\, {\rm s}^{-1}$) and $T_{a}$ and $T_{s}$ are air and surface temperature respectively. Latent heat flux is calculated as

(3)$$F_{{\rm lat}} = \rho_{a}L_{\,f}C_{T}v\lpar q_{a}-q_{s}\rpar {\comma\; }$$

where $L_{f}$ is the latent heat of fusion of ice (${\rm J}\, {\rm kg}^{-1}$), and $q_{a}$ and $q_{s}$ are the air and surface specific humidities.

$C_{T}$ is calculated as

(4)$$\left\{\matrix{ C_{T} = C_{T_{s}}\left(1-{2b'{Ri}\over 1+c\vert {{Ri}}\vert ^{1/2}}\right)\hfill & {\rm if}\ {Ri} \lt 0\hfill \cr C_{T} = C_{T_{s}}\lpar 1+b'{Ri}\rpar ^{-2} \hfill & {\rm if}\ {Ri} \geq 0 {\comma\; } \hfill }\right.$$

where $C_{T_{s}} = 1.3\times 10^{-3}$, $b' = 20$ and c = 50.986 are constants and Ri is the bulk Richardson number,

(5)$${Ri} = {g\lpar T_{a}-T_{s}\rpar \Delta z\over T_{a}v_{a}^{2}} {\comma\; }$$

where $\Delta z$ is equal to 10 m following Ebert and Curry (Reference Ebert and Curry1993).

Humidity is provided as relative (%) by PROMICE weather stations and converted to specific humidity by first obtaining the saturation vapour pressure, $e_{s}\lpar T\rpar$ at temperature T ($^{\circ }{\rm C}$), following Tetens (Reference Tetens1930) as

(6)$$e_{s}\lpar T\rpar =0.611 \times 10^{{7.5T}/\lpar {T+237.3}\rpar } {.}$$

Vapour pressure is obtained using the definition of relative humidity (RH) as

(7)$$e = {RH}e_{s} {.}$$

The mixing ratio of water vapour, w, is calculated after American Meteorological Society (2012) as

(8)$$w = {eR_{d}\over R_{v}\lpar p-e\rpar } {\comma\; }$$

where p is the pressure (Pa), $R_{d}$ is the specific gas constant for dry air (${\rm J}\, {\rm kg}^{-1}\, {\rm K}^{-1}$) and $R_{v}$ is the specific gas constant for water vapour (${\rm J}\, {\rm kg}^{-1}\, {\rm K}^{-1}$). Lastly, the specific humidity, q, can be obtained after American Meteorological Society (2012) using

(9)$$q={w\over w+1} {.}$$

Lake albedo, α, was calculated following Lüthje and others (Reference Lüthje, Feltham, Taylor and Worster2006) based on a two-stream approximation from Taylor and Feltham (Reference Taylor and Feltham2004) and adapted for the generally lower albedo for glacier ice than sea ice as

(10)$$\alpha={9702+1000e^{3.6z_{l}}\over -539+20000e^{3.6z_{l}}} {\comma\; }$$

where $z_{l}$ is the depth of the lake.

Shortwave propagation

The Beer–Lambert law was used to calculate the transfer of shortwave radiation through water and ice as

(11)$$F_{i}=F_{b}e^{-\tau z_{i}}-F_{b}e^{-\tau z_{i+1}} {\comma\; }$$

where $F_{i}$ is the flux at cell i (${\rm W}\, {\rm m}^{-2}$), τ is the shortwave extinction coefficient (${\rm m}^{-1}$) and $z_{i}$ is the depth of cell i. $F_{b}$, the shortwave radiation entering the water column, is calculated from total incoming shortwave radiation as

(12)$$F_{b}=I_{0}\lpar 1-\alpha\rpar F_{{sw}} {\comma\; }$$

where $I_{0}$ is the proportion of shortwave radiation absorbed at the lake surface. The net radiation for the surface cell is then reduced to $\lpar 1- I_{0}\rpar \lpar 1-\alpha \rpar F_{{sw}}$.

Lake convection and indexing

The lake is modelled as turbulently convecting at ${\geq }0.1\, {\rm m}$, following Buzzard and others (Reference Buzzard, Feltham and Flocco2018), using the four thirds rule,

(13)$$F_{c}\lpar T^{\ast }\rpar ={\rm sign}\lpar \bar{T}-T^{\ast }\rpar \lpar \rho C\rpar _{l}J\vert \bar{T}-T^{\ast }\vert ^{{4}/{3}} {\comma\; }$$

where $F_{c}\lpar T^{\ast }\rpar$ is the energy flux at the lake boundary of temperature $T^{\ast }$, $\bar {T}$ is the average lake temperature, $\lpar \rho C\rpar _{l}$ is the volumetric specific heat capacity of water and J is the turbulent heat flux factor ($1.907 \times 10^{-5}\, {\rm ms}^{1}\, {\rm K}^{-{1}/{3}}$). Turbulent heat flux is applied to the first slush or ice cell immediately surrounding the lake.

To correctly apply turbulent mixing, correctly indexing slush, ice and water cells was necessary. This was implemented by progressing upwards from the base of the domain and recording the first and last instances of water between ice sections. This prevented two lakes being separated due to any transient appearance of one slush cell in the middle of the lake and encourages slush and ice encroachment from the upper and lower bounds of the lake.

Lake water is stratified based on temperature following the approach of Saloranta and Andersen (Reference Saloranta and Andersen2007) by using the UNESCO International Equation of State of Seawater (Loucks and van Beek, Reference Loucks and van Beek2005). This equation is applicable here as it deals with both saline content and temperature. Once lake input is complete in stage 3 the entire lake insert is combined with the main model domain in a one-off resizing to simplify model code in other stages.

Snow cover

The snow layer is adapted to be coupled to ice, rather than soil as developed by Essery (Reference Essery2015). Further details, including changing snow conductivity, compaction of snow over time, water content calculation and the initialisation of snow conductivity, can be found within Essery (Reference Essery2015). Realistic values for snow density ($200 {-} 500 \, {\rm kg}\, {\rm m}^{-3}$) are used for the conduction subfunction. Calculations for snow are completed at the beginning of each time step, to allow snow properties to be used subsequently in the conduction subfunction. When the snow layer becomes liquid, the snow layer is eliminated. Following results of sensitivity tests that found the influence of this water to be minor in comparison with meltwater flow into the lake, this water is not introduced into the lake. If refreezing occurs after partial melting of the snow layer, the density and thermal conductivity are updated to account for the proportion of the snow layer that is then solid ice based on the relative proportion of snow, ice and water.

Enthalpy and heat diffusion

By default, the model is initialised as an ice column with temperature linearly interpolated between a surface and basal temperature specified by the user. Temperature is calculated from the cell's enthalpy as below, or in reverse if necessary as

(14)$$\left\{\matrix{ T_{i}={E_{i}\over \rho_{{\rm water}}C_{{\rm ice}}\, {\rm d}z} \hfill & {\rm if}\ E_{i}\leq E_{{\rm ice}_{ \rm i}}\hfill \cr T_{i}=T_{\rm melt} \hfill & {\rm if}\ E_{\rm ice_{\rm i}}\lt E_{i}\lt E_{\rm water_{i}}\hfill \cr T_{i}={E_{i}\over \rho_{\rm water}C_{\rm ice}\, {\rm d}z}-{1\over C_{\rm water}} \lpar C_{\rm ice}T_{\rm melt}+L_{\,f}-C_{\rm water}T_{\rm melt}\rpar \hfill & {\rm if}\ E_{i}\geq E_{\rm water_{i}} {\comma\; } \hfill }\right.$$

where $T_{i}$ is the temperature at each grid cell (K), $E_{i}$ is the enthalpy at each grid cell (${\rm J}\, {\rm m}^{-2}$), $\rho _{\rm water}$ is the density of water (${\rm kg}\, {\rm m}^{-3}$), $C_{\rm ice}$ is the specific heat capacity of ice (${\rm J}\lpar {\rm kg}\, {\rm K}\rpar ^{-1}$), $C_{\rm water}$ is the specific heat capacity of water (${\rm J}\lpar {\rm kg}\, {\rm K}\rpar ^{-1}$) and ${\rm d}z$ is the thickness of the cell. $E_{\rm ice_{i}}$ and $E_{\rm water_{i}}$ are the enthalpy thresholds for ice and water respectively (${\rm J}\, {\rm m}^{-2}$). The model space uses unit metre squares in the x and y direction.

The proportion of each cell that is water, $\lambda _{i}$, is calculated as

(15)$$\left\{\matrix{ \lambda_{i}=0 \hfill & {\rm if}\ E_{i}\leq E_{\rm ice_{i}}\hfill \cr \lambda_{i}={E_{i}-\rho_{\rm water}C_{\rm ice}\, {\rm d}zT_{\rm melt}L_{\,f}\over \rho_{\rm water}\, {\rm d}z} \hfill & {\rm if}\ E_{\rm ice_{i}}\lt E_{i}\lt E_{\rm water_{ i}}\hfill \cr \lambda_{i}=1 \hfill & {\rm if}\ E_{i}\geq E_{\rm water_{i}} {.} \hfill }\right.$$

Heat diffusion is calculated using a backward-time, centred-space approach to enable much longer time steps (tested up to 2 h) whilst maintaining numerical stability when compared to the use of forward-time, centred-space after Lüthje and others (Reference Lüthje, Feltham, Taylor and Worster2006) and Alexiades and Solomon (Reference Alexiades and Solomon1993). The enthalpy change is calculated from temperature difference at each step, which remains valid unless a cell changes from ice to water or vice versa in one time step (an eventuality which depends on the time step and grid cell thickness and unlikely to occur whilst the model remains stable due to the high latent heat of melting/freezing). The one-dimensional heat diffusion equation,

(16)$${\partial T\over \partial t} = K{\partial^2 T\over \partial z^2} {\comma\; }$$

(here excluding surface energy flux which is incorporated earlier in each timestep) is discretised to allow for non-uniform layer spacing as a result of meltwater input, snowfall and larger deep cells as

(17)$${T_{i}^{n+1}-T_{i}^{n}\over \Delta t}=K_{i+{1}/{2}}^{n+1} {T_{i+1}^{n+1}-T_{i}^{n}\over \lpar \Delta z_{i+{1}/{2}}^{n+1}\rpar ^2}-K_{i-{1}/{2}}^{n+1} {T_{i}^{n+1}-T_{i-1}^{n}\over \lpar \Delta z_{i-{1}/{2}}^{n+1}\rpar ^2} {\comma\; }$$

where n is the time step, i is the cell index, $\Delta x$ is the distance between the midpoints of two adjacent cells (m) and $\Delta t$ is the time step. The thermal diffusivity of each cell is calculated as

(18)$$K_{i}^{n}=\lambda_{i}^{n}K_{\rm water}+\lpar 1-\lambda_{i}^{n}\rpar K_{\rm ice}\comma\; $$

and the intermediate values between cells are obtained as depth weighted averages. For convenience, eqn (17) is simplified with the use of the $S_{i}^{n+1}$ term,

(19)$$S_{i}^{n+1}=K_{i}^{n+1}{\Delta t\over {\lpar \Delta z_{i}^{n+1}}\rpar ^2} {\comma\; }$$

to give

(20)$$T_{i}^{n}=-S_{i+{1}/{2}}^{n+1}T_{i+1}^{n+1}+\lpar 1 + S_{i+{1}/{2}}^{n+1} + S_{i-{1}/{2}}^{n+1}\rpar T_{i }^{n+1}-S_{i-{1}/{2}}^{n+1}T_{i-1}^{n+1} {.}$$

A zero-flux boundary condition was specified for the base of the domain with the ice temperature of the bottom cell initially set at $-5^{\circ }{\rm C}$. The default size of the upper cells is 0.1 m, enlarged to 1 m for the lowermost 30 cells to allow a large spatial domain without significantly increasing run-time. Equation (20) can then be entered into a system of simultaneous equations using a tridiagonal matrix generalised to a given number of model layers and solved as a matrix equation at each time step. The temperature of the surface cell is updated at each time step according to eqn (1) prior to heat diffusion.

Sensitivity testing

Testing was conducted to determine the sensitivity of important GlacierLake outputs to parameter uncertainty. Normalised sensitivity coefficients for seven input parameters and six important output measures (Table 3) were calculated following Loucks and van Beek (Reference Loucks and van Beek2005) as

(21)$${\left \vert Q\lpar P_{0}+\Delta P\rpar -Q\lpar P_{0}-\Delta P\rpar \right \vert \over 2 \Delta P}\times {P_{0}\over Q\lpar P_{0}\rpar } {\comma\; }$$

where $Q\lpar P_{0} \pm \Delta P\rpar$ is the output value Q when forced with a parameter of value $P_{0} \pm \Delta P$, $P_{0}$ is the initial parameter value and $\Delta P$ is variation away from $P_{0}$.

Table 3. Parameters, parameter values and outputs used in sensitivity testing

Temperature units are ${}^{\circ }{\rm C}$, depth and thickness in m and density in ${\rm kg}\, {\rm m}^{-3}$.

Figure 3 shows sensitivity coefficients. Most uncertainty is confined to albedo and the $I_{0}$ term (proportion of shortwave radiation absorbed at the lake surface) with the model output being fairly insensitive (sensitivity coefficient ≤0.2 excluding lake depth sensitivity to meltwater input) to other parameters. Many parameters are hard to obtain precisely and are subject to temporal and spatial variability, so this low sensitivity to parameter uncertainty promotes confidence in GlacierLake's results. The $I_{0}$ term merits further consideration and is covered in the discussion.

Fig. 3. Inter-comparison of sensitivity coefficients. Ordered by greatest sum uncertainty in columns and then rows. Abbreviations shown in Table 3. Note exponential scale for colour bar.

Tuning

The only published in-situ measurements of lake bottom ablation and lake depth come from Tedesco and others (Reference Tedesco2012); data which were also used by Banwell and others (Reference Banwell, Arnold, Willis, Tedesco and Ahlstrøm2012a). The record from Lake Ponting (69.589 N, −49.783 E, 962 m, data from 15th–19th June 2011) was compared to a GlacierLake run from 24th September 2009–1st August 2010, which was driven with GC-Net AWS data as processed by Banwell and others (Reference Banwell2012b), with incoming long wave radiation calculated at each time step within GlacierLake using equations from Banwell and others (Reference Banwell2012b). For precipitation input, RACMO2.3p2 was used (Noël and others, Reference Noël2018). The change in the $I_{0}$ value resulted in a large variation in the model output. An $I_{0}$ value of 0.8, slightly greater than the value of 0.6 used by Lüthje and others (Reference Lüthje, Feltham, Taylor and Worster2006), Tedesco and others (Reference Tedesco2012), Benedek (Reference Benedek2014) and Buzzard and others (Reference Buzzard, Feltham and Flocco2018) produced the output seen in Fig. 4. Values set as default following tuning can be found within the main function of the supplied code. GlacierLake is not intended to predict the evolution of specific, individual lakes with absolute accuracy, but this test shows that it is well suited for modelling broad trends and for understanding the behaviour of supraglacial lakes in Greenland. Large-scale comparison to remotely sensed lid thickness is beyond the scope of this study.

Fig. 4. Comparison of modelled lake bottom ablation to measured lake bottom ablation data from Tedesco and others (Reference Tedesco2012). Dashed light grey is modelled basal ice ablation excluding lake formation, black is measured lake bottom ablation. Jumps in modelled ablation around day 175 result from threshold behaviour with cell enthalpy near the slush-water boundary.

Air temperature and precipitation variation

Precipitation and temperature vary across the spatial domain of lake occurrence in Greenland (Bromwich and others, Reference Bromwich, Chen, Bai, Cassano and Li2001; Gledhill and Williamson, Reference Gledhill and Williamson2017) and exert strong controls over ice lid formation within GlacierLake. While arbitrary, the following temperature and precipitation ranges were chosen to observe the behaviour of GlacierLake within, and slightly beyond, conditions where supraglacial lakes may reasonably be expected to form across the GrIS. Testing the relative importance of precipitation and temperature on lid formation enables a range of expected lid thicknesses at the beginning of the melt season to be obtained, as well as to find if any threshold behaviour is apparent. Weather data from PROMICE UPE U, with 2 m total of meltwater input (over 5 days using Lake Ponting data from Tedesco and others, Reference Tedesco2012) was used to form a lake with conditions varied upon lid formation to simulate the effect of precipitation and temperature differences on lid formation from a common starting point. Precipitation was multiplied by a factor ranging from 0 to 3 resulting in start of melt season snow depths between 0 and 3.45 m. Temperature from the day of lid formation to the end of the model run was added to or subtracted to with a range from $+10$ to $-20^{\circ }{\rm C}$. Average temperature over this period with default settings is $-12.4^{\circ }{\rm C}$ giving a range from $-2.4$ to $-32.4^{\circ }{\rm C}$. The maximum lid thickness before lid breakup (or at 298 days of lid coverage if no lid breakup occurs) was recorded for a range of 25 inputs for precipitation and temperature giving 625 values in total.

Figure 5 shows that a change in snow cover between 0 and 3.9 m over winter while air temperature is constant results in a greater range of ice lid thickness than a change in over-winter temperature between $+10$ and $-20^{\circ }{\rm C}$ (0.8–2.5 m against 1.4–2.3 m). Lid thickness varies smoothly with a parameter change with no threshold behaviour observed. Figure 5 also shows that even under the most extreme conditions tested (no snow cover and average winter temperatures of $-32.4^{\circ }{\rm C}$) lid thickness does not exceed 3.3 m. Koenig and others (Reference Koenig2015) find an average ice lid thickness of 1.4 m, with 0.65 m of snow cover present using observations obtained through Operation IceBridge. GlacierLake predicts lid thickness ranging from 1.7–2.7 m depending on temperature when snow cover at the start of the melt season reaches 0.65 m Koenig and others (Reference Koenig2015). This translates to thinner lids (1.5–2.5 m) during April and May when IceBridge flights are performed but suggests GlacierLake may slightly over predict lid thickness.

Fig. 5. Change in the maximum lid thickness with temperature and precipitation variation. Red vertical and horizontal lines indicate default UPE U conditions. White is where lid formation malfunctioned within GlacierLake. White contours show 0.5 m maximum thickness intervals.

Advection and basal freeze-up

Although GlacierLake does not have inbuilt allowance for advection, its effect was simulated by resetting ice temperature to its initial value on the day of lid formation, effectively removing the impact of warming earlier in the model run (Fig. 6). As surface energy flux is effectively isolated from the base upon lid formation this represents the earliest time in the model run where near-surface temperatures are not atmospherically controlled (e.g. cold wave in Fig. 6, left). A lower temperature limit of $-15^{\circ }{\rm C}$ was chosen as limited near-surface ice temperature profiles of land terminating sectors of the GrIS show this to be at the lower end of expected values (Meierbachtol and others, Reference Meierbachtol, Harper and Humphrey2013). Figure 7 shows that including advection increases basal freeze-up by at most 0.1 m compared to runs with no simulated advection. Total basal freeze-up with an initial ice temperature of $-15^{\circ }{\rm C}$ is limited to 0.5 or 0.6 m with simulated advection. Altering the initial and advecting ice temperature made no difference to the magnitude of lid freezing (constant at 1.6 m throughout the test).

Fig. 6. Simulated advection. Panels a and b: no alteration to ice temperature initially set at $-10^{\circ }{\rm C}$. Panels c and d: a constant temperature of $-10^{\circ }{\rm C}$ was set initially, and repeated upon lid formation. Basal freeze-up was 0.4 and 0.5 m respectively. White space in panels b and d is free space before meltwater inflow input. Black line in panels a and c is snow depth, dashed blue line is water depth in snow layer.

Fig. 7. Total basal freeze-up whilst lid is present (stages 4–5) with and without simulated advection. Whilst the trend is upwards, small decreases in basal freeze-up are observed with and without advection included. This is a result of threshold behaviour between slush, water and ice within cells where the temperature is at, or within $0.05^{\circ }{\rm C}$ of $0^{\circ }{\rm C}$. For example a colder initial ice temperature may result in an increase of three slush cells above the ice-slush boundary but a decrease of one ice cell. Refer to Fig. 6 for test setup.

Temperature change at depth

The assumption that the warming influence of supraglacial lakes on underlying ice is not significant (Poinar and others, Reference Poinar, Joughin, Lenaerts and Van Den Broeke2017) has not been tested in the literature. We examine this by extending the 1 m deep cells at the bottom of the model to a total of 50 m (Fig. 8) making a 60 m domain in total. The initial temperature profile is set uniformly to $-10^{\circ }{\rm C}$ with weather data from PROMICE station UPE-U (Fausto and van As, Reference Fausto and van As2019). Figure 8 shows the lake exerts a strong influence (up to $8^{\circ }{\rm C}$) to around 10 m below the lake, however the warming effect beyond 15 m below the lake is limited at less than $1^{\circ }{\rm C}$.

Fig. 8. Ice temperature at depth using AWS data from the PROMICE UPE-U weather station, commencing 1st January 2010. Note the change in depth scale and jump from absolute temperature to temperature change from day 0 between panels b and c. Although the overall model domain extends to 60 m, only the upper 30 m are shown here. Black line in panel a is snow depth, blue dashed line is water depth.

GlacierLake behaviour at UPE-U AWS

To test steady-state behaviour, and the general behaviour of GlacierLake, it was run with PROMICE UPE-U data for 1000 days from 1st of January, 2010, repeating forcing data year on year. Figure 9 shows that although the thermal impact of a lake on underlying ice is limited beyond 15 m it is sufficient to prompt faster lake formation the following melt season. As annual melt exceeds annual refreezing, no steady state arises.

Fig. 9. Example model output including a short duration meltwater inflow input and snow cover, with data from PROMICE UPE-U AWS from 1st January 2010 repeating year-on-year for 1000 days. Panel a shows snow cover where black is snow depth and dashed blue is water depth. Panel b showsa lake temperature profile. White space in panel b is free space before meltwater inflow input.