# Modelling the evolution of an ice sheet’s weathering crust

Tilly Woods<EMAIL_ADDRESS>Mathematical Institute, University of
Oxford, Oxford, UK  Ian J. Hewitt  Mathematical Institute, University of
Oxford, Oxford, UK

###### Abstract

The weathering crust is a layer of porous ice that can form at the surface of
an ice sheet. It grows and decays in response changing weather and climate
conditions, affecting the albedo, the melt rate, and the transport of
meltwater across the surface. To understand this behaviour, we seek time-
dependent solutions to a continuum, thermodynamic model for the porosity,
temperature and thickness of the weathering crust, and the internal and
surface melt rates. We find solutions using a numerical enthalpy method,
presented in this study. We use idealised ‘switching’ and sinusoidal forcings
to explore the different dynamics exhibited during growth and decay, the
timescales involved, and the impact of diurnal vs. annual variations. The
results demonstrate qualitative agreement with observations, and provide
insight into the relative importance of different surface heat fluxes during
the growth and decay of the crust. The model therefore provides a useful tool
for exploring the response of the weathering crust to climate change.

Keywords: ice sheet; phase change; enthalpy method; climate.

##  I Introduction

The weathering crust is a layer of porous ice that forms at the surface of ice
sheets and glaciers [1, 2]. The structure is formed by shortwave radiation
from the sun penetrating below the surface of the ice, where it acts as an
internal heat source, causing internal melting [3, 4]. Weathering crusts occur
on bare (snow-free) ice in the ablation zone [5], where there is a net loss of
ice mass at the surface, for example on the south-western Greenland ice sheet
[e.g. 6]. The weathering crust is on the order of 1 m thick, with impermeable
ice below, and is mostly saturated with meltwater [2, 7, 8].

The weathering crust is of interest for multiple reasons. Firstly, it affects
the transport of meltwater across the surface of the ice sheet [2, 8], into
nearby surface streams. This in turn impacts the ice sheet mass loss and
consequent effect on sea level [9]. For example, surface meltwater can be
transiently stored or refrozen in the weathering crust [10, 2], and it can be
routed to the bed of the ice sheet via moulins [10], where it affects the ice
sheet motion [11, 12] before being expelled downstream (e.g. into the ocean)
through a subglacial drainage network [e.g. 13]. Moreover, microbial biomass
is transported in the surface meltwater [7, 8], affecting the provision of
nutrients to downstream environments [14] and the carbon fluxes between the
ice sheet and the atmosphere [15].

Secondly, the appearance and structure of the weathering crust affect how much
radiation is reflected or absorbed by the ice [16], which in turn affects how
much the ice melts. For example, a weathering crust reflects more radiation
than the glazed, bare ice that is left behind when the weathering crust is
removed [16]. The reflectivity of a surface is quantified by the albedo – the
proportion of outgoing to incoming shortwave radiation. Darker surfaces have a
lower albedo and absorb more radiation, leading to increased melt. In addition
to the ice structure, microbes (e.g. algae) and impurities within the
weathering crust interact with the ice structure [16, 17] and act to lower the
albedo [18, 19, 20].

The weathering crust is a naturally transient structure. It only forms in the
summer months when there is enough shortwave radiation to cause internal
(subsurface) melting [e.g. 7]. On a shorter timescale, the weathering crust
responds to rapidly changing weather conditions, growing and decaying on the
order of hours and days [e.g. 21]. Clear, sunny conditions favour formation of
the crust, whereas warm, windy, cloudy conditions favour removal [1, 21].

In order to understand the behaviour of the weathering crust and what impact
it has on melting and the transport of meltwater, we need to understand its
transient response to changing weather and climate conditions. Existing
studies have primarily focussed on understanding the weathering crust through
observations from the field [e.g. 1, 2, 7]. Until now, there have been limited
attempts to model the weathering crust [21, 22, 23, 24], and these have tended
to focus on specific aspects. Schuster [21], Hoffman  _et al._ [22, 23], Woods
and Hewitt [24] all present one-dimensional (vertical) models forced by the
surface energy balance components (e.g. shortwave and longwave radiation and
turbulent heat fluxes). Schuster [21] presented a physically-based 6-layer
model for the density of the weathering crust, with density changes driven by
internal melting resulting from internal shortwave radiation. Hoffman  _et
al._ [22, 23] modelled the temperature within the ice using the heat equation
with a heat source provided by internal shortwave radiation, gaining
predictions of surface lowering and internal melt extent over time. The more
complex continuum model we presented in Woods and Hewitt [24] models the
surface lowering, internal and surface melt, weathering crust porosity and
thickness, and internal ice temperature simultaneously, based on mass
conservation, energy conservation, and internal shortwave radiation. That
study focussed only on steady solutions. However, studying the full time-
dependent problem is crucial for extending our understanding of the physical
processes governing weathering crust evolution, and exploring the response of
the weathering crust, melting and runoff to climate change.

In the present study, we extend our work from Woods and Hewitt [24] to seek
time-dependent solutions of the weathering crust model presented in that
study. We begin in Section II by summarising the model and presenting the
model equations relevant for this work. We then present a numerical enthalpy
method for solving the time-dependent equations and summarise relevant
background in Sections III and IV. In Section V we consider idealised
‘switching’ experiments to gain insight into the physical processes
controlling the growth and decay of the weathering crust. In Section VI we use
sinusoidal forcings to explore behaviour closer to the reality. Finally, we
discuss our findings, avenues for future work, and conclude in Section VII.

##  II Time-dependent model

Figure 1: The setup for the weathering crust model in (a) the general two-
dimensional case and (b) the one-dimensional case. In the simplest case, the
domain consists of two regions: a porous region (the weathering crust) above a
cold, solid ice region. Variables are defined within the main text. In the
one-dimensional domain, we transform coordinates to $Z=h(t)-z$, meaning that
the ice surface is fixed and the cold ice below gets advected towards to
surface at rate $-\dot{h}$. Any surface meltwater is removed instantaneously.

In Woods and Hewitt [24], we presented a fully three-dimensional continuum
model for the evolution of the weathering crust, accounting for mass
conservation, energy conservation and internal shortwave radiation. We model
the porous weathering crust and the solid ice below it, as shown in Figure 1.
Here, we present a simplified version of that model under the following
assumptions. We restrict ourselves to one (vertical) dimension, and assume
that the weathering crust is completely saturated with meltwater. Given that
there is always a layer of impermeable ice beneath the weathering crust, this
results in no flow of meltwater.

The volume fraction of the water-filled pores is given by the porosity $\phi$,
and the volume fraction of ice is $1-\phi$. The temperature $T$ is assumed to
be homogeneous at the pore scale, and we assume local thermal equilibrium,
meaning that if both ice and the water are present (i.e. in the porous
weathering crust) then both phases are at the melting temperature
$T_{\mathrm{m}}$. In regions of pure ice ($\phi=0$), we must have $T\leq
T_{\mathrm{m}}$, and in regions of pure water ($\phi=1$), $T\geq
T_{\mathrm{m}}$. Furthermore, we assume that the ice and water have shared
properties, such as density $\rho$, thermal conductivity $k$ and specific heat
capacity $c$.

The ice surface (at the interface with the air) is at $z=h(t)$, which evolves
due to surface melting (see kinematic condition and surface energy balance
below). This melting and subsequent surface lowering produce a thin layer of
pure meltwater above the surface, which we do not model explicitly. It is
assumed that all such surface meltwater immediately runs off laterally into
nearby surface streams. As done in [24], we transform coordinates to
$Z=h(t)-z$, the depth below the ice surface. The equations are presented using
the $Z$ coordinate, where the ice surface is now at $Z=0$. The model setup is
shown in Figure 1.

###  II.1 Mass conservation

The mass conservation equations for the water and ice reduce to

$\rho\Bigg{(}\frac{\partial\phi}{\partial
t}+\dot{h}\frac{\partial\phi}{\partial Z}\Bigg{)}=m_{\mathrm{int}},$ (1)

where $m_{\mathrm{int}}$ is the internal melt rate, and $\dot{h}$ ($\leq 0$)
is the rate of change of surface height. The second term on the left-hand side
accounts for advection of ice towards the surface. In the physical $z$-frame,
the ice is stationary and melting causes the surface $z=h(t)$ to lower at rate
$-\dot{h}$. In the $Z$-frame, in which the ice surface is fixed at $Z=0$, this
manifests itself as the ice advecting backwards towards $Z=0$.

###  II.2 Energy equation

The energy equation is

$\rho c\Bigg{(}\frac{\partial T}{\partial t}+\dot{h}\frac{\partial T}{\partial
Z}\Bigg{)}=k\frac{\partial^{2}T}{\partial Z^{2}}-\frac{\partial F}{\partial
Z}-\mathcal{L}m_{\mathrm{int}},$ (2)

where $F$ is the net downwelling internal shortwave radiation (discussed
below), and $\mathcal{L}$ is the latent heat of fusion for water. The second
term on the left-hand side again accounts for advection of ice towards the top
of the domain $Z=0$. The terms of the right-hand side represent the conduction
of heat following Fourier’s law, the internal heat source provided by the
internal shortwave radiation $F$, and the energy sink associated with melting
the ice. (This can also be a source, when $m_{\mathrm{int}}<0$, which
corresponds to refreezing).

###  II.3 Internal shortwave radiation

Shortwave radiation from the sun - in particular its penetration below the ice
surface - is a key ingredient for forming a weathering crust. Some shortwave
radiation penetrates below the surface and is absorbed internally, acting as
an internal heat source and producing the porous structure of the weathering
crust. Some is absorbed at the surface, and some is reflected back into the
atmosphere. In our model, we partition the absorbed shortwave radiation
between the surface and the subsurface using a parameter $\chi$ [23, 25, 24].

The total amount of shortwave radiation absorbed is $(1-a)Q_{\mathrm{si}}$,
where $Q_{\mathrm{si}}$ is the incoming shortwave radiation and $a$ is the
albedo (the ratio of outgoing to incoming shortwave radiation). Of that,
$\chi(1-a)Q_{\mathrm{si}}$ is absorbed at the surface and
$(1-\chi)(1-a)Q_{\mathrm{si}}$ is absorbed internally, distributed over depth,
as determined by the flux term $F$, which we now describe.

Both downwards and upwards internal radiation can be modelled explicitly using
a two-stream approximation, which accounts for the absorption and scattering
within the ice [e.g. 26, 27, 28, 29, 24]. If the absorption and scattering
coefficients are assumed uniform, as we do here, they can be solved for
analytically, and the result is that the _net_ downwards radiation $F$ is
given by

$F=(1-\chi)(1-a)Q_{\mathrm{si}}e^{-\kappa Z},$ (3)

where $\kappa=\sqrt{\alpha^{2}+2\alpha r}$ is the extinction coefficient,
related to the absorption coefficient $\alpha$ and the scattering coefficient
$r$. In this study, we also assume $\alpha$ and $r$, and consequently
$\kappa$, to be constant in time. The internal shortwave radiation, and the
internal heat source it provides, decay exponentially with depth (this is
sometimes referred to as Beer’s law). The explicit modelling of the upwelling
radiation coming back out of the ice surface also determines a relationship
between the absorption and scattering coefficients and the albedo [24]; that
is

$a=\frac{(1-\chi)(\kappa-\alpha)}{(1-\chi)\kappa+(1+\chi)\alpha}.$ (4)

###  II.4 Kinematic condition

The kinematic condition for the ice at the ice surface $Z=0$ is

$-\rho(1-\phi)\dot{h}=\rho m_{\mathrm{surf}}\quad\mathrm{at}\quad Z=0,$ (5)

where $m_{\mathrm{surf}}$ is the surface melt rate. We assume that any excess
water at the ice surface runs off instantaneously (rather than allowing it to
pool on the surface); combining the rate of surface melting
$m_{\mathrm{surf}}$ with the rate $-\phi\dot{h}$ at which already-liquid water
is released from the lowering surface, this means that the runoff occurs at
rate $r_{\mathrm{off}}=-\dot{h}$.

###  II.5 Surface energy balance

The surface energy balance is written (see Woods and Hewitt [24]) as

$k\frac{\partial T}{\partial
Z}+\chi(1-a)Q_{\mathrm{si}}+Q_{0}-\upsilon(T-T_{\mathrm{m}})=\rho\mathcal{L}m_{\mathrm{surf}}\quad\mathrm{at}\quad
Z=0,$ (6)

where $Q_{\mathrm{si}}$ is the incoming shortwave radiation, and
$Q_{0}-\upsilon(T-T_{\mathrm{m}})$ is a linear approximation to the net
incoming longwave radiation and turbulent heat flux. It is possible to
calculate $Q_{\mathrm{si}}$, $Q_{0}$ and $\upsilon$ from observed weather data
(see Woods and Hewitt [24]). In this study, however, we prescribe
$Q_{\mathrm{si}}$ and $Q_{0}$ as idealised, time-dependent forcings, and take
$\upsilon$ to be constant.

The way that the surface energy balance is interpreted depends on whether or
not the surface is melting. When the surface is melting, the surface energy
balance determines the surface melt rate $m_{\mathrm{surf}}>0$, and the
boundary condition for the temperature is $T=T_{\mathrm{m}}$ at $Z=0$. When
the surface is not melting, $m_{\mathrm{surf}}=0$ and the surface energy
balance provides the boundary condition for the temperature $T\leq
T_{\mathrm{m}}$ at $Z=0$.

###  II.6 Far field boundary conditions

We prescribe that the ice is cold at depth,

$T\to T_{\infty}\quad\mathrm{as}\quad Z\to\infty,$ (7)

where $T_{\infty}<T_{\mathrm{m}}$ is the deep temperature of the ice.
Consistent with previous assumptions, this also means that it is solid ice,
$\phi=0$.

##  III Enthalpy method

We are interested in solutions of the time-dependent model presented in
Section II. To find these, we reframe the problem in terms of enthalpy, and
solve it numerically.

There are three different types of regions we could have in the model domain:
(1) cold regions of only ice ($T<T_{\mathrm{m}}$, $\phi=0$), (2) temperate
regions with a mixture of ice and water ($T=T_{\mathrm{m}}$, $0<\phi<1$), and
(3) warm regions of only water ($T\geq T_{\mathrm{m}}$, $\phi=1$). We know the
porosity in the only-water and only-ice regions, and we know the temperature
in the ice-and-water regions. Hence, in each of the three regions, we only
need to solve for one of $\phi$ and $T$: we solve the mass conservation
equation (1) for $\phi$ in ice-and-water regions, and solve the energy
equation (2) for $T$ in the only-ice and only-water regions. To implement this
directly would require tracking the free-boundaries between neighbouring
regions and solving different equations in different parts of the domain.

The enthalpy method [30, 31, 32, 33, 34] is an alternative approach which we
choose to use here. It removes the need to track the free boundaries between
the three different types of regions, and allows us to solve a single equation
for the enthalpy $H$ in the whole domain. This also has the advantage of
accommodating situations in which there are multiple free boundaries, and in
which these boundaries come and go, automatically.

The enthalpy is defined as

$H=\rho c(T-T_{\mathrm{m}})+\rho\mathcal{L}\phi,$ (8)

i.e. the sum of the sensible and latent heat. The porosity $\phi$ and
temperature $T$ can be recovered from the enthalpy $H$ as follows:

$\begin{split}&\phi=\begin{cases}0,&H\leq 0,\\\
\frac{H}{\rho\mathcal{L}},&0<H<\rho\mathcal{L},\\\
1,&H\geq\rho\mathcal{L},\end{cases}\\\
&T=\begin{cases}T_{\mathrm{m}}+\frac{H}{\rho c},&H\leq 0,\\\
T_{\mathrm{m}},&0<H<\rho\mathcal{L},\\\ T_{\mathrm{m}}+\frac{H}{\rho
c}-\frac{\mathcal{L}}{c},&H\geq\rho\mathcal{L}.\end{cases}\end{split}$ (9)

The three different cases correspond to there being only ice ($H\leq 0$), both
ice and water ($0<H<\rho\mathcal{L}$), and only water
($H\geq\rho\mathcal{L}$). The first holds in the cold, solid ice beneath the
weathering crust, and the second holds in the weathering crust. In the
simplest weathering crust setup (as shown in Figure 1), we have only these two
regions. In general, we could have a combination of only-ice, only-water and
ice-and-water layers (such as the examples shown in Figure 2). The enthalpy
method means that a more complicated situation like this should be no more
difficult to solve than the simple two layer case.

By combining the mass conservation equation (1) and the energy equation (2) to
eliminate the internal melt rate $m_{\mathrm{int}}$, we get a single equation
for the enthalpy,

$\frac{\partial H}{\partial t}+\dot{h}\frac{\partial H}{\partial
Z}=k\frac{\partial^{2}T}{\partial Z^{2}}-\frac{\partial F}{\partial Z},$ (10)

which is solved everywhere in the domain. Once we have solved the enthalpy
equation (10) for $H$, we can recover $T$ and $\phi$ using (9). The boundaries
between ice-only and ice-and-water region are identified as where $H=0$, and
boundaries of water-only regions are identified by $H=\rho\mathcal{L}$.

Details of the numerical implementation of the enthalpy method are found in
Appendix A. This includes a description of how we decide which of the two
different cases (melt or no-melt) of the surface energy balance (6) to use in
an implicit timestepping setting, and how to deal with the diffusion term in
the enthalpy equation (10).

Figure 2: Possible transient configurations of the weathering crust model.
Configurations (a) and (b) are encountered in this study. Configurations (c)
and (d) are other, more general, possibilities.

##  IV Solution preliminaries

In Sections V and VI below, we use solutions to the time-dependent model
presented in Section II to understand the physical processes involved in the
growth and decay of the weathering crust. Before doing this, we discuss how
the forcings in the model relate to the real world, and we formulate the
steady-state solutions.

###  IV.1 Connecting model and observations

It is qualitatively observed that clear, sunny conditions favour formation of
the weathering crust, whereas warm, windy, cloudy conditions favour removal
[1, 21, 16]. These conditions can be related to the forcings $Q_{\mathrm{si}}$
and $Q_{0}$ used in our model. $Q_{\mathrm{si}}$ is the shortwave radiation
(direct radiation from the sun), which is high in sunny conditions and low in
cloudy or dark conditions. It is always positive. $Q_{0}$ is a combination of
longwave radiation and turbulent heat fluxes [24], and can be positive or
negative. Turbulent heat fluxes account for the fact that, if the air is
warmer than the surface, winds can cause turbulent eddies which transfer heat
from the air to the surface [5]. The turbulent heat flux is greater when the
wind is stronger and the air temperature is warmer. Therefore, $Q_{0}$ is
larger in warm, windy conditions. Increased cloud cover also leads to an
increase in $Q_{0}$, since longwave radiation emitted by the Earth gets
reflected back down to the surface in the presence of clouds. Hence, based on
observations [1, 21], we expect that large $Q_{\mathrm{si}}$ favours
weathering crust formation, and large, positive $Q_{0}$ favours weathering
crust removal.

Note that negative $Q_{0}$ corresponds to conditions, such as cold air
temperatures, in which the surface would freeze in the absence of shortwave
radiation. The heat flux $Q_{0}$ is completely absorbed at the surface,
whereas the absorption of $Q_{\mathrm{si}}$ is split between the surface and
the subsurface (with the proportion given by the parameter $\chi$). Therefore,
a positive value of $Q_{0}$ will result in surface melting, and sufficiently
negative value of $Q_{0}$ (compared to $Q_{\mathrm{si}}$) will result in
surface freezing.

###  IV.2 Steadily melting solutions

Figure 3: The types of steady ($\dot{h}=0$) and steadily melting ($\dot{h}<0$)
solutions that are possible for different values of $Q_{\mathrm{si}}$ and
$Q_{0}$. The plot shows the regions of the $Q_{\mathrm{si}}$-$Q_{0}$ plane for
which there are: (A) steady (no melt) solutions with pure ice, (B) steady (no
melt) solutions with a subsurface pond, (C) steadily melting solutions with a
porous weathering crust, and (D) steadily melting solutions with no crust.
Contours of weathering crust thickness $Z_{\mathrm{m}}$ are shown in region C.
Typical temperature profiles in each region are shown, with the porosity in
the weathering crust also shown in region C. In regions A and B, the far-field
temperature $T_{\infty}$ is determined as part of the solution. In regions C
and D, $T_{\infty}=-10^{\mathrm{o}}$C is prescribed. All other parameter
values are as listed in Appendix B. The coloured arrows in the main plot show
the change in forcing used in Figures 4, 5, 10 (magenta), Figure 6 (cyan) and
Figure 9 (yellow). Adapted from Woods and Hewitt [24].

These ideas about the conditions under which a weathering crust forms or not
were demonstrated by Woods and Hewitt [24] using steadily melting solutions to
the one-dimensional weathering crust problem. These are solutions to the model
in Section II with constant forcings $Q_{\mathrm{si}}$ and $Q_{0}$ such that
the surface lowers at a constant rate $\dot{h}<0$, the weathering crust
occupies a region between the surface and a constant depth $Z_{\mathrm{m}}$,
and the temperature and porosity profiles in the $Z$-frame are constant in
time. The ice melts in such a way that it appears the column of ice is simply
being translated downwards over time (see Figure 2 in Woods and Hewitt [24]).
The bottom of the weathering crust is located at
$z=z_{\mathrm{m}}=h-Z_{\mathrm{m}}$. Figure 3 [modified from 24] shows the
region of $Q_{0}$-$Q_{\mathrm{si}}$ parameter space for which steadily melting
solutions with a weathering crust exist (region C) and the corresponding
weathering crust depth. It is also possible to have steadily melting states
with no weathering crust, just solid ice (region D). Equations for the
boundaries of region C are presented by Woods and Hewitt [24].

###  IV.3 Steady solutions (non-melting)

In regions A and B of Figure 3, there is not enough energy to melt the
surface, so a steadily melting state is not possible. Instead, it is possible
to have steady states with no surface lowering ($\dot{h}=0$). Depending on the
values of $Q_{\mathrm{si}}$ and $Q_{0}$, these steady states consist of
either: pure ice (region A), or a layer of surface ice below which there is
pure water (region B). The latter is not physically realistic, but is included
here for completeness.

The temperature profile in the non-melting steady states is found by solving
the energy equation (2) in the steady state. This gives

$T=T_{\mathrm{m}}+\frac{1}{\upsilon}\big{(}(1-a)Q_{\mathrm{si}}+Q_{0}\big{)}+\frac{1}{\kappa
k}(1-\chi)(1-a)Q_{\mathrm{si}}\big{(}1-e^{-\kappa Z}\big{)}.$ (11)

In the steadily melting solutions, ice is being advected towards the surface
in the $Z$-frame as a result of the physical surface $z=h(t)$ lowering at rate
$-\dot{h}>0$. This allows us to prescribe the temperature $T_{\infty}$ at
depth. In contrast, there is no advection in the non-melting steady states, so
it is not possible to prescribe the temperature at depth, but only to
prescribe that $\frac{\partial T}{\partial Z}$ tends to zero as $Z\to\infty$.
The far-field temperature $T_{\infty}$ is instead determined as part of the
solution. The sign of $T_{\infty}-T_{\mathrm{m}}$ determines whether we have a
pure-ice solution ($T_{\infty}\leq T_{\mathrm{m}}$, region (b)) or subsurface
pond solutions ($T_{\infty}\geq T_{\mathrm{m}}$), with the boundary between
the two regions given by the line

$Q_{0}+\bigg{(}\frac{\upsilon}{k\kappa}(1-\chi)+1\bigg{)}(1-a)Q_{\mathrm{si}}=0.$
(12)

###  IV.4 Time-dependent solutions

The key result from above is that it is the balance between shortwave
radiation $Q_{\mathrm{si}}$ and longwave radiation and turbulent heat fluxes
$Q_{0}$ that controls the growth and decay of the weathering crust. However,
real ice sheets rarely have constant forcing. In the following sections, we
use the time-dependent model to explore the dynamic growth and decay of the
weathering crust when conditions transition between those favouring growth and
those favouring decay.

In the time-dependent solutions we are seeking, we remove the assumption made
in the steadily melting solutions that the surface of the ice is melting. It
is possible to have a solid layer of ice above the top of the weathering crust
(see Figure 2(b)). In this case, we call the location of the top of the
weathering crust (the porous region) $z=z_{\mathrm{f}}=h-Z_{\mathrm{f}}$,
where $Z_{\mathrm{f}}$ is the depth of the top of the weathering crust below
the surface. The weathering crust thickness is defined as
$Z_{\mathrm{t}}=Z_{\mathrm{m}}-Z_{\mathrm{f}}$. Schematics of possible
transient configurations are shown in Figure 2.

##  V ‘Switching’ forcing - growth and decay of the weathering crust

###  V.1 Setup

To gain some initial understanding of the growth and decay of the weathering
crust, we consider ‘switching’ experiments, in which we start with a steadily
melting state (with constant forcings $Q_{\mathrm{si}}^{\mathrm{old}}$,
$Q_{0}^{\mathrm{old}}$), then instantaneously switch to a new set of constant
forcings ($Q_{\mathrm{si}}^{\mathrm{new}}$, $Q_{0}^{\mathrm{new}}$), and
investigate how the weathering crust responds. Such changes in forcing are
demonstrated in Figures 4, 5, 6.

Based on Figure 3, we expect that a steadily melting weathering crust can be
removed by (i) melting from the surface if we switch to forcings in region D
(low $Q_{\mathrm{si}}$, positive $Q_{0}$), and by (ii) freezing from the
surface if we switch to forcings in regions A or B (negative enough $Q_{0}$).

To grow a weathering crust, we can do the reverse: switch to forcings for
which there is a steadily melting crust solution.

###  V.2 Growth of a weathering crust

Figure 4: Growth of a weathering crust from a steadily melting, no-crust
solution. The second and third plots in panel (a) show the evolution of the
temperature $T$ and porosity $\phi$ profiles over times. The bottom of the
weathering crust, $z=z_{\mathrm{m}}$, is shown as a black line. Contours of
$T$ ($1^{\mathrm{o}}$C spacing) and $\phi$ (0.05 spacing) are shown in grey.
The corresponding forcings $Q_{\mathrm{si}}$, $Q_{0}$, and the resulting
surface lowering $-\dot{h}$, surface melt rate $m_{\mathrm{surf}}$ and the
weathering crust thickness $Z_{\mathrm{t}}$ are shown in the first and fourth
plots. Snapshots are shown of (b) temperature and (c) porosity at 10 day
intervals, from day 0 to day 70. The forcings were switched at time $t=0$.

We first consider the growth of a weathering crust from a no-crust steadily
melting state, as shown in Figure 4. We start with no shortwave radiation
$Q_{\mathrm{si}}$, then turn it on to $Q_{\mathrm{si}}=100$ W m-2 whilst
keeping $Q_{0}=50$ W m-2 fixed. The snapshots of the temperature and porosity
profiles in panels (b) and (c) show that the weathering crust develops to
close to its final thickness quite quickly (over a few days), but with very
low porosity initially. The porosity in the crust then continues to increase,
as does the temperature in the top $\sim 6$ m of the ice, with the crust
thickness changing very little. The resulting porosity and temperature
profiles are always roughly exponential, due to the exponential decay of the
internal shortwave radiation $F$ (3), which acts as the internal heat source.

###  V.3 Removal via surface melting

Figure 5: Removal of a steadily melting weathering crust by melting from the
surface. The second and third plots in panel (a) show the evolution of the
temperature $T$ and porosity $\phi$ profiles over times. The bottom of the
weathering crust, $z=z_{\mathrm{m}}$, is shown as a black line. Contours of
$T$ ($1^{\mathrm{o}}$C spacing) and $\phi$ (0.05 spacing) are shown in grey.
The corresponding forcings $Q_{\mathrm{si}}$, $Q_{0}$, and the resulting
surface lowering $-\dot{h}$, surface melt rate $m_{\mathrm{surf}}$ and the
weathering crust thickness $Z_{\mathrm{t}}$ are shown in the first and fourth
plots. Snapshots are shown of temperature (b) and porosity (c) at 10 day
intervals, from day 0 to day 40. The forcings were switched at time $t=0$.

We now consider removing the weathering crust by melting from the surface.
Figure 5 shows the reverse experiment to that shown in Figure 4. This time, we
are removing a steadily melting crust by turning off the shortwave radiation
$Q_{\mathrm{si}}$ \- an exact reversal of the forcings just considered. We see
that removal happens through a combination of melting from the surface and
freezing from below. Also note, from the porosity contour plot, that the
porosity at each point in space remains unchanged until that point is either
frozen or melted. That is, the lines of constant $\phi$ lie horizontal in the
physical $z$ frame. The only way the porosity can change in the bulk of the
domain is through melting due to the internal heat source provided by the
internal radiation $F$. Therefore, when $Q_{\mathrm{si}}$ (and hence $F$) is
turned off, the internal melting ceases. Freezing also changes the porosity,
but this can only happen at interfaces between porous and solid ice (e.g. at
$Z=Z_{\mathrm{m}}$), and results in a jump in porosity (see Figure 5(c)). The
jump in porosity satisfies the Stefan condition (in the $z$ frame)

$-k\frac{\partial T}{\partial
z}=-\rho\mathcal{L}\phi\dot{z}_{\mathrm{m}}\quad\mathrm{at}\quad
z=z_{\mathrm{m}},$ (13)

where $z_{\mathrm{m}}=h-Z_{\mathrm{m}}$ is the bottom of the weathering crust,
$\phi$ is the porosity value on the crust side of the interface (with $\phi=0$
on the other side), and the temperature derivative is that on the no-crust
side (being equal to zero on the crust side). The Stefan condition (13) is
satisfied automatically by the solution of the enthalpy equation when
discretised in a conservative way.

The slight jump in the surface lowering rate $-\dot{h}$ that can be observed
at the instant the crust is removed (around $t=30$ d) can be explained by a
combination of (1) the freezing from below causing a discontinuity of $\phi$
at the freezing front at the bottom of the crust, and (2) the surface lowering
rate $-\dot{h}$ depending on the surface porosity $\phi$ (with
$m_{\mathrm{surf}}$ being constant). When the final thickness of crust is
removed, the surface porosity jumps from non-zero to zero, but
$m_{\mathrm{surf}}$ remains the same, resulting in a jump in $\dot{h}$.
Physically, this is because the same amount of energy is now being used to
melt solid ice rather than porous ice, so a smaller decrease in surface height
is achieved for the same amount of energy (see discussion).

###  V.4 Removal via surface freezing

Figure 6: Removal of a steadily melting weathering crust by melting from the
surface. The second and third plots in panel (a) show the evolution of the
temperature $T$ and porosity $\phi$ profiles over times. The bottom of the
weathering crust, $z=z_{\mathrm{m}}$, is shown as a black line. Contours of
$T$ ($1^{\mathrm{o}}$C spacing) and $\phi$ (0.05 spacing) are shown in grey.
The corresponding forcings $Q_{\mathrm{si}}$, $Q_{0}$, and the resulting
surface lowering $-\dot{h}$, surface melt rate $m_{\mathrm{surf}}$ and the
weathering crust thickness $Z_{\mathrm{t}}$ are shown in the first and fourth
plots. Snapshots are shown of temperature (b) and porosity (c) at 5 day
intervals, from day 0 to day 15. The forcings were switched at time $t=0$.

The final case we consider is the switching experiment in which we remove a
steadily melting weathering crust by freezing from the surface, as well as
from below. We start with the same steadily melting state as used in the
removal-by-melting test, and again switch to $Q_{\mathrm{si}}=0$ W m-2, but
this time the new value of $Q_{0}$ takes negative, instead of positive, values
- we are moving into the left-hand region of Figure 3 in which there is no
surface melting. The solution is shown in Figure 6. In this case, the upper
boundary of the weathering crust does not coincide with the ice surface
$z=h(t)$. Instead, there is a layer of solid ice above the weathering crust,
and a freezing front at $z=z_{\mathrm{f}}(t)$. As in Figure 6, when
$Q_{\mathrm{si}}=0$ W m-2, there is no internal melting, so the porosity at
each point in the crust remains unchanged until it jumps to zero when it is
passed by one of the freezing fronts (the top $z=z_{\mathrm{f}}$, or the
bottom $z=z_{\mathrm{m}}$). As a result, the porosity profile simply gets
‘truncated’ over time, as shown particularly in Figure 6(c).

Despite the internal heat source being removed (by turning off the shortwave
radiation $Q_{\mathrm{si}}$), the temperature deeper in the ice increases over
time (see Figure 6(b)). This is because the more dominant effect is that the
surface lowering has stopped, meaning that there is no longer any advection of
cold ice at depth upwards towards the surface. The freezing of the water in
the weathering crust releases heat, which conducts downwards and initially
warms the ice lower down. Meanwhile the cooling of the surface due to negative
$Q_{0}$ causes the temperature near the surface to decrease.

###  V.5 Dependence on the value of $Q_{0}^{\mathrm{new}}$

Figure 7: Removal of the weathering crust by ((a), (b)) melting from the
surface and ((c), (d)) freezing from the surface, for different values of
$Q_{0}$. Initially, there was a steadily melting crust with
$Q_{\mathrm{si}}=25$ W m-2 and $Q_{0}=50$ W m-2 ((a), (c), (e), (f)), and
$Q_{\mathrm{si}}=100$ W m-2 and $Q_{0}=50$ W m-2 ((b), (d), (g), (h)). In both
cases, at $t=0$, the forcings were switched to
$Q_{\mathrm{si}}^{\mathrm{new}}=0$ W m-2 and values of
$|Q_{0}^{\mathrm{new}}|$ from 0 to 100 W m-2, in steps of 10 W m-2 (positive
for melting, negative for freezing). Solid lines show the top ($z=h(t)$ when
melting, $z=z_{\mathrm{f}}$ when freezing) and bottom ($z=z_{\mathrm{m}}$) of
the weathering crust for different values of $Q_{0}$. The analytical
approximation (15) for $z=h(t)$ for small $\phi$ is also shown (dotted).
Panels (e) and (g) show how the time $t_{\mathrm{rem}}$ taken to remove the
crust depends on $Q_{0}^{\mathrm{new}}$. Panels (f) and (h) shown the porosity
profiles of the initial steadily melting solutions, taking values of 0.05 and
0.61 at the surface, respectively.

We also explore how the _rate_ of the removal of a weathering crust depends on
the change in the values of $Q_{\mathrm{si}}$ and $Q_{0}$. Firstly, we
investigate the impact of $Q_{0}^{\mathrm{new}}$. We explore what happens when
we turn off the shortwave radiation (set $Q_{\mathrm{si}}^{\mathrm{new}}=0$ W
m-2) and vary the new value of $Q_{0}^{\mathrm{new}}$ (Figure 7). We consider
both melting ($Q_{0}^{\mathrm{new}}>0$) and freezing
($Q_{0}^{\mathrm{new}}<0$) from the surface. (Figures 5 and 6 have already
shown examples of such solutions for the specific cases
$Q_{0}^{\mathrm{new}}=50$ W m-2 and $Q_{0}^{\mathrm{new}}=-50$ W m-2.)

Broadly, we see that the time taken to remove the weathering crust by melting
is shorter when $Q_{0}^{\mathrm{new}}$ is more positive (‘+’ in Figure 7 (e),
(g)), and the time taken to remove the weathering crust by freezing is shorter
when $Q_{0}^{\mathrm{new}}$ is more negative (‘o’ in Figure 7 (e), (g)).

Comparing the melting and freezing cases, we see that the freezing from the
surface removes the weathering crust much more quickly than melting from the
surface when the porosities are small ($\leq 0.05$, left-hand side of Figure
7). Smaller porosities mean that much more energy is required to melt the ice
fraction (0.95) than freeze the water fraction (0.05). On the contrary, when
the porosities are larger (0.61 at the surface initially, right-hand side of
Figure 7), the removal times are more similar. This time, melting removes the
crust more quickly than freezing, since there is less ice to melt than water
to freeze.

The fact that it takes less time to remove the weathering crust by freezing
from the surface when $Q_{0}^{\mathrm{new}}$ is more negative can be explained
by the surface energy balance (6) (with $m_{\mathrm{surf}}=0$ in this case).
To maintain the surface energy balance, a decrease in $Q_{0}$ must be
compensated for by either decreasing the surface temperature $T$ or by
increasing the heat flux $k\frac{\partial T}{\partial Z}$ towards the surface
from the ice below. This acts to cool the internal ice.

The reason that the time taken to remove the crust by melting is shorter when
$Q_{0}^{\mathrm{new}}$ is more positive is that increasing $Q_{0}$ increases
the amount of surface melting. The surface melt rate given by the surface
energy balance (6) (when $Q_{\mathrm{si}}=0$, as here) is

$m_{\mathrm{surf}}=\frac{Q_{0}}{\rho\mathcal{L}}.$ (14)

This in turn leads to an increase in the rate of surface lowering $-\dot{h}$,
according to the kinematic condition (5). For small surface porosity, the
surface lowering rate can be approximated analytically from (5) as

$-\dot{h}=\frac{m_{\mathrm{surf}}}{1-\phi}\approx m_{\mathrm{surf}},$ (15)

giving an approximately constant rate of surface lowering. This is observed in
Figure 7(a), where there is good agreement between the analytical (dotted) and
numerical (solid) solutions for $z=h$. In this case, the small $\phi$
assumption is valid, with $\phi\approx 0.05$ at the surface in the initial
steadily melting state. However, the agreement is less good in Figure 7(b)
where $\phi\approx 0.61$ at the surface, making the small $\phi$ assumption
less valid.

We also observe that the evolution of the bottom of the weathering crust
$z=z_{\mathrm{m}}$ is independent of the value of $Q_{0}$. We can see this
mathematically by writing down the problem for the temperature in the solid
ice below the weathering crust in the absence of shortwave radiation. This is

$\rho c\frac{\partial T}{\partial t}=k\frac{\partial^{2}T}{\partial z^{2}}$
(16)

in $z<z_{\mathrm{m}}(t)$, with three boundary conditions: the far-field
condition (7) as $z\to-\infty$, and the Stefan condition (13), together with
continuity of temperature $T=T_{\mathrm{m}}$, at $z=z_{\mathrm{m}}$.

We have one more boundary condition than required to solve the heat equation,
enabling us to also solve for the moving boundary $z=z_{\mathrm{m}}(t)$. We
note that this problem is independent of the boundary conditions at $z=h$,
hence $z=z_{\mathrm{m}}$ can be solved for independent of $Q_{0}$ (although,
with the given initial temperature and porosity profiles, it is unfortunately
not straightforward to solve analytically).

###  V.6 Dependence on the value of $Q_{\mathrm{si}}^{\mathrm{new}}$

Figure 8: Removal of the weathering crust by melting from the surface for
different values of $Q_{\mathrm{si}}$. Initially, there was a steadily melting
crust with $Q_{\mathrm{si}}=25$ W m-2 and $Q_{0}=50$ W m-2. The forcings were
switched at $t=0$ to $Q_{0}=50$ W m-2 (unchanged) and various values of
$Q_{\mathrm{si}}$, from 0 to 14 W m-2, in steps of 2 W m-2. Panel (a) shows
the evolution of the top (solid, $z=h$) and bottom (dashed,
$z=z_{\mathrm{m}}$) of the weathering crust for different values of
$Q_{\mathrm{si}}^{\mathrm{new}}$. Panel (b) show how the time
$t_{\mathrm{rem}}$ taken to remove the crust depends on the new value of
$Q_{\mathrm{si}}$. Panel (c) shows the porosity profile of the initial
steadily melting solutions, taking a value of 0.05 at the surface.

We now explore how the removal of the weathering crust is affected by the new
value of $Q_{\mathrm{si}}$, which affects both the internal melting (via the
internal radiation $F$ (3)) and surface melting (via the surface energy
balance (6)). We do this by starting with the steadily melting state with
$Q_{\mathrm{si}}^{\mathrm{old}}=25$ W m-2 and $Q_{0}^{\mathrm{old}}=50$ W m-2,
then switching $Q_{\mathrm{si}}$ to a new, lower value, whilst keeping $Q_{0}$
fixed. Figure 8(a) shows that the effect of this is to change the rate of
freezing at $z=z_{\mathrm{m}}$, whereas the rate of surface lowering changes
very little. The surface lowering is approximately constant in this example
because all the values of $Q_{\mathrm{si}}^{\mathrm{new}}$ are small (ranging
from 0 to 14 W m-2), so the surface melt rate (14) is dominated by the $Q_{0}$
term. This is amplified by $Q_{\mathrm{si}}$ being multiplied by
$\chi(1-a)\approx 0.2<1$, representing the fact that not all the shortwave
radiation is absorbed at the surface: some is reflected ($aQ_{\mathrm{si}}$)
and some is absorbed internally ($(1-\chi)(1-a)Q_{\mathrm{si}}$). If
$Q_{0}^{\mathrm{new}}$ was closer to 0 and $Q_{\mathrm{si}}^{\mathrm{new}}$
was larger, the change in surface lowering with
$Q_{\mathrm{si}}^{\mathrm{new}}$ would be more noticeable.

The main effect of varying $Q_{\mathrm{si}}^{\mathrm{new}}$ is to change the
internal heat source, generated by the internal shortwave radiation $F$, given
by (3). The smaller $Q_{\mathrm{si}}^{\mathrm{new}}$ is, the greater is the
reduction in the internal heat source (as $Q_{\mathrm{si}}$ switches from
$Q_{\mathrm{si}}^{\mathrm{old}}$ to $Q_{\mathrm{si}}^{\mathrm{new}}$ in (3)),
and hence the crust gets frozen from below more quickly. There is a maximum
value of $Q_{\mathrm{si}}^{\mathrm{new}}$ for which the crust gets removed
(12.7 W m-2 here), given by the boundary between regions C and D in the
contour plot in Figure 3 (with the equation for this line given by Woods and
Hewitt [24]). Any value larger than this results in a new steadily melting
crust of a different thickness being formed, determined by Figure 3. As
$Q_{\mathrm{si}}^{\mathrm{new}}$ tends towards its maximum value for removal,
the removal time tends to infinity (Figure 8(b)).

##  VI Periodic forcing

Although the switching experiments in Section V provide insight into the basic
behaviours of weathering crust growth and removal, such sudden forcing
transitions are not very realistic. In reality, we do not expect the forcings
to remain constant for sustained periods of time. As previously discussed, the
weathering crust and the weather conditions that force it are very variable,
changing on the order of hours and days [e.g. 35], with the weathering crust
growing and decaying over short timescales as a result [21]. Therefore, a more
relevant test is sinusoidal variations of forcings, representing, for example,
diurnal or annual variations.

In all of the experiments considered in this section, a sinusoidal forcing is
run until a roughly periodic state is reached. Note that the forcings here are
not necessarily realistic diurnal or annual variations (since these are not
sinusoidal in amplitude), but serve to demonstrate the kind of periodic
behaviour that is possible from our model. As noted by Schuster [21], the
initial conditions of the weathering crust (e.g. how porous it is) will likely
have a strong effect on the amount and type (surface or internal) of melt
resulting from given forcing conditions. Periodic solutions provide a way to
examine the behaviour of the weathering crust whilst eliminating the impact of
the initial conditions.

###  VI.1 Periodic growth and complete removal

Figure 9: Periodic weathering crust generated by sinusoidal forcing, with
complete removal each period by surface melting. The second and third plots of
panel (a) show the evolution of the temperature $T$ and porosity $\phi$
profiles over times. The bottom of the weathering crust, $z=z_{\mathrm{m}}$,
is shown by a black line. Contours of $T$ ($0.2^{\mathrm{o}}$C spacing) and
$\phi$ ($0.5\times 10^{-3}$ spacing) are shown in grey. The corresponding
forcings $Q_{\mathrm{si}}$, $Q_{0}$, and the resulting surface lowering
$-\dot{h}$, surface melt rate $m_{\mathrm{surf}}$ and weathering crust
thickness $Z_{\mathrm{t}}$ are shown in the first and fourth plots. Snapshots
are shown of temperature (b) and porosity (c) at 1 hour intervals over the
course of 1 day.

Figure 9 shows an example where the weathering crust is periodically
completely removed (by melting from the surface and freezing from below), then
regrown, as a result of daily variations in $Q_{\mathrm{si}}$, with $Q_{0}$
remaining constant. In this example, the porosity in the weathering crust is
very small ($\sim 10^{-3}$), due to the limited time for it to grow. A key
point to note is the non-symmetric behaviour in the growth and removal of the
crust, with the crust growing at a faster rate than it decays. We see that
when the crust is growing (blue to green snapshots in panel (c)), the porosity
smoothly increases from zero at the bottom of the crust, with the depth of the
crust and the porosity at each depth increasing smoothly (as in Figure 4). In
contrast, when the crust is removed, the freezing from below results in a jump
in the porosity at the bottom of crust, according to the Stefan condition (13)
(as in Figures 5 and 6). Compared to the switching examples in Figures 5 and 6
in which the $Q_{\mathrm{si}}=0$, here $Q_{\mathrm{si}}>0$, so there continues
to be internal melting whilst freezing from the bottom of the crust is
occurring. We see this in the fact that the lines of constant porosity slope
downwards in the $t$-$z$ plane here, compared to being horizontal in Figures 5
and 6. Also note the lag in the crust thickness. The crust thickness continues
to increase for a short while after $Q_{\mathrm{si}}$ has begun to decrease,
and $Q_{\mathrm{si}}$ almost reaches its maximum before the crust begins to
re-form. This is explained by the fact that the weathering crust does not
respond instantaneously to the changes in $Q_{\mathrm{si}}$ and $Q_{0}$, as
demonstrated in the switching experiments in Section V. It takes time for heat
to conduct through the domain.

###  VI.2 Diurnal vs. annual forcing

Figure 10: Periodic weathering crust generated by periodic forcing, for the
same forcings but over different timescales (a) diurnal and (b) annual. The
simulations have been run to reach a roughly periodic state: 100 days for the
daily forcings, 10 years for the yearly. The second and third plots in panels
(a) and (b) show the evolution of the temperature $T$ and porosity $\phi$
profiles over times. The bottom of the weathering crust, $z=z_{\mathrm{m}}$,
is shown by a black line. Contours of $T$ ($0.5^{\mathrm{o}}$C spacing) and
$\phi$ ($0.02$ spacing) are shown in grey. The corresponding forcings
$Q_{\mathrm{si}}$, $Q_{0}$, and the resulting surface lowering $-\dot{h}$,
surface melt rate $m_{\mathrm{surf}}$ and weathering crust thickness
$Z_{\mathrm{t}}$ are shown in the first and fourth plots. Panel (c) compares
the time mean solutions to the steadily melting solution for the mean
forcings.

We now explore the impact of the period of the forcing. Figure 10 shows two
periodic solutions forced by the same periodic forcings ($Q_{\mathrm{si}}$
varying between 0 and 100 W m-2, $Q_{0}$ constant at 50 W m-2), but with
different periods: 1 day and 1 year. We see that with the 1 day period, the
crust thickness undergoes small variations ($\sim 0.1$ m), with the crust
remaining present and of approximately constant thickness throughout. With the
1 year period, the crust gets completely removed (and regrown) periodically.
This is in line with the switching experiment in Figure 5, in which it took
$\sim 1$ month to remove the steadily melting weathering crust corresponding
to $Q_{\mathrm{si}}=100$ W m-2, $Q_{0}=50$ W m-2 by turning off the shortwave
radiation (i.e. the switching version of the periodic forcing in Figure 10).

We can also compare the periodic solutions to steadily melting solutions for
the same mean forcing (Figure 10(c)). We see that the mean $\phi$, $T$ and
$Z_{\mathrm{m}}$ of the periodic solution are close to, but not equal to, the
steadily melting solution of the mean forcing ($Q_{\mathrm{si}}=50$ W m-2,
$Q_{0}=50$ W m-2). Firstly, the mean porosity is non-zero as far down as the
greatest depth to which the porosity is non-zero at any time during the
period. This is deeper than the mean crust depth in the steadily melting
porosity profile. The depth to which non-zero porosity is reached is largest
for the annual forcing since there is more time for heat to conduct deeper
into the ice. However, the annual mean crust thickness is lower than both the
steadily melting crust thickness and the mean thickness for the diurnal period
because the crust is completely removed for periods of time. The crust would
continue to thin further during some of this time if there was more ice
available to melt (i.e. if the crust was thicker in the first place). Since
this is not possible, it leads to a proportionately thinner mean crust
thickness. On the other hand, the mean crust thickness for the diurnal period
is quite close to the steadily melting crust thickness. This is due to the
diurnal variations being faster than the timescale over which the crust
significantly evolves, so the oscillations in forcing are effectively averaged
out. The mean crust thickness for the diurnal variations is slightly larger
that the steadily melting crust thickness, and we attribute this to the
asymmetric growth and decay of the weathering crust. As seen in Figure
9(a),(c), growth occurs more quickly than removal, favouring a crust that is
thicker on average than if the constant mean of the forcing was applied.

###  VI.3 Superimposed diurnal and annual forcing

As well as considering diurnal and annual variations in isolation, it is
informative to consider diurnal variations superimposed onto annual
variations. This is representative of the increased solar radiation during the
day compared to at night, and during the summer compared to winter. When
modelling the evolution of a system that responds to the solar cycle, such as
the weathering crust, we would like to know how important the short-term
diurnal variations in forcing are. When considering longer timescales (many
years, e.g. to investigate the response to changing climate conditions), it
would be preferable to be able to ignore the diurnal variations, enabling us
to take longer timesteps in numerical computations. Here, we explore the
impact of these superimposed diurnal variations on the growth and decay of the
weathering crust.

To do this, we take the shortwave radiation to be

$Q_{\mathrm{si}}(t)={Q}_{\mathrm{si}}^{\mathrm{year}}(t)(1+\sin(2\pi
t/t_{\mathrm{day}})),$ (17)

where ${Q}_{\mathrm{si}}^{\mathrm{year}}(t)$ is the annually-varying diurnal-
mean,

${Q}_{\mathrm{si}}^{\mathrm{year}}(t)=\bar{Q}_{\mathrm{si}}+A\sin(2\pi
t/t_{\mathrm{year}}),$ (18)

with mean $\bar{Q}_{\mathrm{si}}$ and amplitude $A$. Here, $t_{\mathrm{day}}$
is the period of 1 day, $t_{\mathrm{year}}$ is the period of 1 year. This form
of $Q_{\mathrm{si}}$ takes a minimum of 0 in each daily cycle (at night). We
take ${Q}_{\mathrm{si}}^{\mathrm{year}}(t)$ to be the same as the shortwave
radiation forcing used for the annual variations in Figure 10, with
$A=\bar{Q}_{\mathrm{si}}=50$ W m-2. The result is shown in Figure 11.

Figure 11: Periodic weathering crust generated by diurnal sinusoidal
variations superimposed onto annual sinusoidal variations (panel (a)). The
diurnal variations are such that $Q_{\mathrm{si}}$ goes to zero each night.
The simulations have been run for 22 years to reach a periodic state. Panel
(b) shows the _difference_ between the solution in panel (a) and the solution
with only annual variations shown in Figure 10 (‘$\Delta=$ diurnal-and-annual
– annual-only’). The second and third plots in panels (a) and (b) show the
evolution of the temperature $T$ and porosity $\phi$ profiles over time. The
black dashed line shows the bottom of the weathering crust,
$z=z_{\mathrm{m}}$. The corresponding forcings $Q_{\mathrm{si}}$, $Q_{0}$, and
the resulting surface lowering $-\dot{h}$, surface melt rate
$m_{\mathrm{surf}}$ and weathering crust thickness $Z_{\mathrm{t}}$ are shown
in the first and fourth plots. Panel (c) shows weathering crust thickness
$Z_{\mathrm{t}}$, the surface lowering $-h$, cumulative surface melt
$M_{\mathrm{surf}}$, and cumulative internal melt $M_{\mathrm{int}}$ for the
diurnal-superimposed-onto-annual forcings (panel (a)) and the purely annual
forcings (Figure 10(b)).

The main effect of superimposing diurnal variations onto the annual variations
is to create diurnal oscillations in the porosity, temperature, weathering
crust thickness, internal and surface melt rates, and surface lowering rate.
The overall behaviour of the solution is similar to that with the annual
variations only.

The oscillations in temperature are of order $0.01^{\mathrm{o}}$C, and the
oscillations in porosity are of order $0.001$ (Figure 11(b)). The oscillations
in weathering crust thickness (Figure 11(b)) are once again of order 0.1 m,
consistent with the diurnal variations in Figure 10. As already seen in Figure
10(a),(c), the daily mean of these oscillations is slightly larger than the
crust thickness variation resulting from the annual-only forcing (Figure
11(c)), attributed to the asymmetric growth and decay of the weathering crust.

We see that the cumulative surface lowering $-h$ ($=\int-\dot{h}\
\mathrm{d}t$), cumulative surface melt $M_{\mathrm{surf}}$ ($=\int
m_{\mathrm{surf}}\ \mathrm{d}t$) and cumulative internal melt
$M_{\mathrm{int}}$ ($=\frac{1}{\rho}\int\int\max(m_{\mathrm{int}},0)\
\mathrm{d}z\ \mathrm{d}t$) over 1 year are largely unaffected by the diurnal
variations (Figure 11(c)). The surface lowering has contributions from surface
melting and internal melting. Internal melting lowers the porosity, and the
resulting non-zero surface porosity means that more surface lowering occurs
than would be accounted for by surface melting alone. Any already-liquid water
at the surface (produced by internal melt) adds to the surface melt and runs
off as the surface lowers. In general, not all the internal meltwater runs off
at the surface; some of it can instead refreeze. In this example, however, the
amount that refreezes is small, so almost all the internal meltwater is
removed at the surface. Hence the surface lowering is approximately equal to
the sum of the internal and surface melting.

Nevertheless, slightly more internal melting does occur in the diurnal+annual
solution than the annual-only solution, due to the diurnal refreezing and then
re-melting at the bottom of the weathering crust. This is demonstrated by the
difference in porosity between the solutions with and without diurnal
variations (Figure 11(b)). The amount of additional internal melting produced
as a result of this diurnal freeze-melt cycle is small, since the porosity at
the bottom of the weathering crust is close to zero. For larger variations in
$Q_{0}$ as well as $Q_{\mathrm{si}}$, it can be larger.

##  VII Conclusion and Discussion

We have explored time-dependent solutions of the one-dimensional continuum
model for the weathering crust we presented in Woods and Hewitt [24]. We have
used an enthalpy method to numerically solve the model equations, removing the
need to track the boundaries between the porous, pure-ice and pure-water
regions of the domain. Our choice of idealised ‘switching’ and sinusoidal
forcings has enabled us to gain insight into the physical processes
controlling the growth and decay of the weathering crust in response to
changing weather conditions (represented in our model by the shortwave
radiation $Q_{\mathrm{si}}$ and the longwave radiation and turbulent heat
fluxes $Q_{0}$). We find that our model produces broadly realistic behaviours
when compared qualitatively to observations (see below), suggesting that it
could be used as a tool to make predictions about the response of the
weathering crust to future weather and climate conditions. There are, however,
improvements and extensions that could be made to the model, discussed below.

###  VII.1 Growth and decay behaviour and timescales

In line with observations [21, 1], we have seen that it is the balance between
shortwave radiation ($Q_{\mathrm{si}}$) and longwave radiation and turbulent
heat fluxes ($Q_{0}$) that controls the growth and decay of the weathering
crust, with clear, sunny conditions (high $Q_{\mathrm{si}}$, low $Q_{0}$)
favouring growth, and warm, windy, cloudy conditions favouring decay (low
$Q_{\mathrm{si}}$, high $Q_{0}$). Furthermore, we have observed in our model
different mechanisms of growth and decay. The weathering crust grows via
internal melting (Figure 4), resulting from some shortwave radiation
penetrating beneath the surface. Decay/removal of the crust can occur by (1)
melting from the surface (Figure 5), or (2) freezing from the surface (Figure
6), combined with (3) freezing from below (Figures 5 and 6).

Moreover, there is an asymmetry in the growth and decay behaviour, made clear
by our study of sinusoidal forcings. Growth maintains a smooth porosity
profile ($\phi\to 0$ smoothly at the bottom of the weathering crust), whereas
decay by freezing (from above or below) results in a jump in porosity as phase
change occurs. In a periodic setting, growth occurs more quickly than decay.
This is true across a number of different cases: (1) when the weathering crust
is periodically completely removed by surface melting (Figure 9), (2) when it
oscillates in thickness but is never removed (Figure 10(a)), or (3) when it is
periodically completely removed by surface freezing (not shown here). The
asymmetry in timescales means that a weathering crust produced by a sinusoidal
forcing about a mean has a thicker weathering crust on average than if the
forcing was held constant at the mean value, unless the crust is completely
removed during the cycle (Figure 10).

The timescale for removal in the specific examples studied here is typically
on the order of several days. The removal time decreases with increasing
magnitude of $Q_{0}$ (with $Q_{0}>0$ corresponding to removal by surface
melting, and $Q_{0}<0$ corresponding to removal by surface freezing). The
porosity, particularly at the surface, also affects the rate of removal of the
weathering crust: a more porous crust is easier to remove by melting, and
harder to remove by freezing.

###  VII.2 Actual vs. potential surface lowering

In Section V.3, we observed that the surface porosity $\phi_{\mathrm{surf}}$
affects the rate of surface lowering, with a jump in $\phi$ at the surface
resulting in a jump in the rate $-\dot{h}$ of surface lowering. In Section
VI.3, we saw that the surface lowering is less than the surface melting due to
the non-zero surface porosity (see kinematic condition (5)). This ties into
ideas discussed by Schuster [21]. They studied the different stages of
weathering crust growth and decay in terms of the ratio $R$ between the actual
(observed) surface lowering and the potential surface lowering (based on the
total energy available to melt). In our case,

$R=\frac{\mathrm{actual\ surface\ lowering}}{\mathrm{potential\ surface\
lowering}}=\frac{-\dot{h}}{\big{(}(1-a)Q_{\mathrm{si}}+Q_{0}\big{)}/\rho\mathcal{L}},$
(19)

where the potential surface lowering comes from the total melt energy: the sum
of the energy absorbed at the surface, $\chi(1-a)Q_{\mathrm{si}}+Q_{0}$, and
the energy absorbed internally, $(1-\chi)(1-a)Q_{\mathrm{si}}$. It is expected
that $R<1$ during crust growth because some of the melt energy is used for
internal melting, which is not visible at the surface. During crust decay
(removal via melting), it is expected that $R>1$ because the porous ice
requires less energy to remove than solid ice [21]. Phrasing this differently,
the observed surface lowering would suggest more (less) melting is occuring
during decay (growth) than is truly the case [1]. This can lead to errors in
surface ablation (mass loss) predictions based on observed surface lowering,
although these errors resulting from short-term crust variations are expected
to cancel out over longer timescales [1].

Using the kinematic condition (5) and the expression (14) for the surface melt
rate $m_{\mathrm{surf}}$, we can rewrite $R$ as

$R=\Bigg{(}\frac{\chi(1-a)Q_{\mathrm{si}}+Q_{0}}{(1-a)Q_{\mathrm{si}}+Q_{0}}\Bigg{)}\Bigg{(}\frac{1}{1-\phi_{\mathrm{surf}}}\Bigg{)},$
(20)

where the first fraction is the proportion of the total energy available for
melt that is used to melt the surface. This shows that the key factors causing
a discrepancy between the observed and expected surface lowering are (1) the
fact that some of the shortwave radiation penetrates and is absorbed below the
surface, determined by the parameter $\chi$, and (2) the high surface porosity
of the weathering crust.

###  VII.3 Future/ongoing work

####  VII.3.1 Including drainage of internal meltwater

One key process missing from the current model is lateral flow of meltwater
_within_ the weathering crust, and the possibility for the crust to become
partially saturated (have air as well as water in the pore space), as occurs
in reality [2, 8]. We assume that surface meltwater runs off instantaneously,
but any internal meltwater can only stay in place and potentially refreeze
later on. The modelling work of Hoffman  _et al._ [23] demonstrated that
lateral drainage of internal meltwater in the weathering crust is vital for
correctly reproducing observed surface lowering, porosity and temperature in
an Antarctic setting. They found that, without lateral drainage, surface
lowering and porosity are underestimated. The drainage of meltwater internally
enables the weathering crust to become unsaturated and hence refreeze with
non-zero porosity. Less energy is then required to melt the porous ice,
compared to solid ice, during the next day/year. This results in increased
surface lowering (see previous discussion of the effect of porosity on surface
lowering).

In future work, we will explore this effect by including a parametrisation for
lateral drainage in our one-dimensional model. Including a lateral drainage
parametrsiation also requires including a model for gravity-driven vertical
drainage, to fill the ‘gaps’ left by the internally drained meltwater. The
vertical and lateral drainage can both be modelled using Darcy’s law for flow
through a porous medium.

####  VII.3.2 Accumulation

During the winter, snowfall will accumulate on the frozen surface on which the
weathering crust will form in the summer. The layer of snow can insulate the
ice below during the winter, and some energy that would otherwise be used to
melt the ice is required to melt the snow at the start of the summer. Our
model currently neglects such accumulation, but incorporating this is a target
for future work.

####  VII.3.3 Real-world forcing

The forcings used in this study were idealised ‘switching’ and sinusoidal
forcings. Although informed by the conditions under which we expect growth and
removal of the weathering crust (see Figure 3), these forcings do not exactly
replicate reality. In ongoing work, we are extending this study by forcing the
model with observed weather data from the Greenland Ice Sheet, and comparing
to more specific observations of surface lowering, porosity and temperature.

## Acknowledgments

Tilly Woods thanks Brian Wetton for valuable discussions about the enthalpy
method.

## Funding

This work was supported by the UK Engineering and Physical Sciences Research
Council [doctoral studentship to TW].

##  Appendix A Details of numerical enthalpy method

We solve the enthalpy equation (10) using a semi-implicit finite volume scheme
in MATLAB. The flux $F$ is evaluated on cell edges, and the enthalpy $H$ is
evaluated on cell centres. We include $\dot{h}$ and $F$ explicitly, $H$ in the
advection term implicitly, and $T$ in the diffusion term semi-implictly, as
follows

$\frac{H^{m+1}-H^{m}}{\Delta t}=-\dot{h}^{m}\frac{\partial H^{m+1}}{\partial
Z}+k\frac{\partial^{2}\big{(}T_{H}(H^{m+1})\big{)}}{\partial
Z^{2}}-\frac{\partial F^{m}}{\partial Z},$ (21)

where the superscripts $m$ and $m+1$ represent solutions at the current and
next timestep. $T_{H}$ is a function that we use to express $T$ in terms of
$H$, based on the transformation (9).

We want to include $H$ in the diffusion term implicitly, but how $T$ is
expressed in terms of $H$ depends on whether there is only ice, both ice and
water, or only water (the three cases in (9)), which we do not know until we
know the value of $H$. Similarly, the implementation of the surface boundary
condition is different depending on whether or not there is surface melting,
which, again, we do not know a priori. We tackle these issues by using an
‘enthalpy flag’ and a ‘melt flag’, and iterating within each timestep to
determine the appropriate values of these.

###  A.1 Melt flag

At each timestep, the surface could be in one of two ‘cases’, which we label
with a flag $f_{\mathrm{melt}}$: the melt case ($f_{\mathrm{melt}}=1$) or the
no-melt case ($f_{\mathrm{melt}}=0$), with different surface boundary
conditions being applied in each case. In the melt case, the surface is at the
melting temperature, $T=T_{\mathrm{m}}$ at $Z=0$ is our boundary condition,
and the surface energy balance (6) tells us the surface melt rate
$m_{\mathrm{surf}}>0$. In the no-melt case, there is no surface melting
($m_{\mathrm{surf}}=0$), $T\leq T_{\mathrm{m}}$ at the surface, and the
surface energy balance (6) is the boundary condition.

###  A.2 Enthalpy flag

Similarly, at each timestep, each point in the domain can be in one of the
three states in (9). We label these using a flag $\mathbf{f}_{H}$, referred to
as the enthalpy flag:

$\mathbf{f}_{H}=\begin{cases}1,&\mathrm{if\ there\ is\ only\
ice}\quad(\phi=0,\ T\leq T_{\mathrm{m}},\ H\leq 0),\\\ 2,&\mathrm{if\ there\
is\ both\ ice\ and\ water}\quad(0<\phi<1,\ T=T_{\mathrm{m}},\
0<H<\rho\mathcal{L}),\\\ 3,&\mathrm{if\ there\ is\ only\ water}\quad(\phi=1,\
T\geq T_{\mathrm{m}},\ H\geq\rho\mathcal{L}).\end{cases}$ (22)

Note that $\mathbf{f}_{H}$ is a vector with its size equal to the number of
grid cells.

###  A.3 Iterations within timesteps

To ensure we use the correct form of the surface boundary condition and the
correct conversion from $T$ to $H$, we iterate the melt flag
$f_{\mathrm{melt}}$ and enthalpy flag $\mathbf{f}_{H}$ as follows. At the
start of each timestep, we assume the values of the flags to be the same as
the previous timestep (e.g. if the surface was melting in the previous
timestep then we assume it is in the new timestep too). Then within each
timestep, we carry out 3 iterations. In each iteration, we

  1. 1.

Solve for $H^{m+1}$ using the values of $f_{\mathrm{melt}}$ and
$\mathbf{f}_{H}$ from the previous iteration. If $f_{\mathrm{melt}}=1$,
$T=T_{\mathrm{m}}$ is used as the surface boundary condition. If
$f_{\mathrm{melt}}=0$, the surface energy balance (6) is used with
$m_{\mathrm{surf}}=0$. The value of $\mathbf{f}_{H}$ tells us how to express
$T$ in the enthalpy equation (10) in terms of $H$, following (9).

  2. 2.

Update the melt flag if the values of $m_{\mathrm{surf}}$ and
$T_{\mathrm{surf}}$ (surface temperature) in the solution contradict the
assumed melt flag.

     * •

Switch to $f_{\mathrm{melt}}=0$ if we assumed $f_{\mathrm{melt}}=1$ and
calculated $m_{\mathrm{surf}}\leq 0$ (i.e. if we assumed there was surface
melting but calculated none).

     * •

Switch to $f_{\mathrm{melt}}=1$ if we assumed $f_{\mathrm{melt}}=0$ and
calculated $T_{\mathrm{surf}}\geq T_{\mathrm{m}}$ (i.e. if we assumed there
was no melting but calculated that the surface was above the melting
temperature).

  3. 3.

Update the enthalpy flag $\mathbf{f}_{H}$ if $T$ and $\phi$ in the solution
contradict the assumed enthalpy flag.

     * •

If we assumed only ice ($\mathbf{f}_{H}=1$) but found $T>T_{\mathrm{m}}$
($H>0$), set $\mathbf{f}_{H}=2$ (try assuming there is both ice and water).

     * •

If we assumed there is both ice and water ($\mathbf{f}_{H}=2$) but found
$\phi<0$ ($H<0$), set $\mathbf{f}_{H}=1$ (try assuming there is only ice).

     * •

If we assumed there is both ice and water ($\mathbf{f}_{H}=2$) but found
$\phi>1$ ($H>\rho\mathcal{L}$), set $\mathbf{f}_{H}=3$ (try assuming there is
only water).

     * •

If we assumed only water ($\mathbf{f}_{H}=3$) but found $T<T_{\mathrm{m}}$
($H<0$), set $\mathbf{f}_{H}=2$ (try assuming there is both ice and water).

  4. 4.

If the assumed melt flag and the enthalpy flag at each point in the domain
were ‘correct’ (i.e. the values from the solution are equal to what we
assumed), stop iterating.

Note that the inequalities in step 2 are weak. This is to prevent the surface
from getting ‘stuck’ in a no-melt state or vice versa. Having a solution with
$m_{\mathrm{surf}}=0$ and $T_{\mathrm{surf}}=T_{\mathrm{m}}$ is possible in
both states, so having equality is not a contradiction. However, if
$m_{\mathrm{surf}}=0$ is achieved having previously been in a state with
$m_{\mathrm{surf}}>0$, it suggests that the surface will likely want to switch
to being frozen in the next timestep, and not allowing this can cause issues.

##  Appendix B Parameter values

The parameter values used throughout this study are listed in Table 1.

Quantity | Symbol | Default value | Units
---|---|---|---
Density of ice | $\rho$ | 910 |  kg m-3
Specific heat capacity of ice at 0oC | $c$ | 2097 |  J kg-1 K-1
Thermal conductivity | $k$ | 2.1 |  kg m-3 K-1
Melting temperature | $T_{\mathrm{m}}$ | 273.15 | K
Latent heat of melting | $\mathcal{L}$ | 334000 |  m2 s-2
Albedo | $a$ | 0.6 |
Proportion of radiation absorbed at the surfacea | $\chi$ | 0.36 |
Extinction coefficient of internal radiationb | $\kappa$ | 1.5 |  m-1
Absorption coefficientc | $\alpha$ | 0.2637 |  m-1
Scattering coefficientc | $r$ | 4.1338 |  m-1
Incoming shortwave radiation flux | $Q_{\mathrm{si}}$ | various |  W m-2
Turbulent heat and longwave radiation flux | $Q_{0}$ | various |  W m-2
Effective heat transfer coefficientd | $\upsilon$ | 14.8 |  W m-2 K-1
Temperature at depth | $T_{\infty}$ | -10 |  oC
Table 1: Default values of the parameters in the model. Unless otherwise
stated, these are the parameter values used. aMotivated by Hoffman  _et al._
[23]. bTaylor and Feltham [29]. cCalculated from $a$, $\kappa$ and $\chi$ as
in Woods and Hewitt [24]. dCalculated as $\nu=C+4\epsilon\sigma
T_{\mathrm{m}}^{3}$ [24], where $C=10.3$ W m-2 K-1 is a heat transfer
coefficient [5], $\epsilon=0.97$ is the surface emissivity [31],
$\sigma=5.7\times 10^{-8}$ W m-2 K-4 is the Stefan-Boltzmann constant, and
$T_{\mathrm{m}}=273.15$ K is the melting temperature of ice [5]. Many of the
ice properties are taken from Cuffey and Paterson [5].

## References

  * Muller and Keeler [1969] F. Muller and C. M. Keeler, Errors in short-term ablation measurements on melting ice surfaces, Journal of Glaciology 8, 91 (1969).
  * Cooper  _et al._ [2018] M. G. Cooper, L. C. Smith, A. K. Rennermalm, C. Mige, L. H. Pitcher, J. C. Ryan, K. Yang, and S. W. Cooley, Meltwater storage in low-density near-surface bare ice in the Greenland ice sheet ablation zone, Cryosphere 12, 955 (2018).
  * Brandt and Warren [1993] R. E. Brandt and S. G. Warren, Solar-heating rates and temperature profiles in Antarctic snow and ice, Journal of Glaciology 39, 99 (1993).
  * Cooper  _et al._ [2021] M. G. Cooper, L. C. Smith, A. K. Rennermalm, M. Tedesco, R. Muthyala, S. Z. Leidman, S. E. Moustafa, and J. V. Fayne, Spectral attenuation coefficients from measurements of light transmission in bare ice on the Greenland Ice Sheet, Cryosphere 15, 1931 (2021).
  * Cuffey and Paterson [2010] K. M. Cuffey and W. S. B. Paterson,  _The Physics of Glaciers_ , 4th ed. (Elsevier, Burlington, MA, 2010).
  * van den Broeke  _et al._ [2008] M. van den Broeke, P. Smeets, J. Ettema, C. van der Veen, R. van de Wal, and J. Oerlemans, The Cryosphere Partitioning of melt energy and meltwater fluxes in the ablation zone of the west Greenland ice sheet, The Cryosphere 2, 179 (2008).
  * Cook  _et al._ [2016] J. M. Cook, A. J. Hodson, and T. D. Irvine-Fynn, Supraglacial weathering crust dynamics inferred from cryoconite hole hydrology, Hydrological Processes 30, 433 (2016).
  * Irvine-Fynn  _et al._ [2021] T. D. Irvine-Fynn, A. Edwards, I. T. Stevens, A. C. Mitchell, P. Bunting, J. E. Box, K. A. Cameron, J. M. Cook, K. Naegeli, S. M. Rassner, J. C. Ryan, M. Stibal, C. J. Williamson, and A. Hubbard, Storage and export of microbial biomass across the western Greenland Ice Sheet, Nature Communications 12, 3960 (2021).
  * Andersen  _et al._ [2015] M. L. Andersen, L. Stenseng, H. Skourup, W. Colgan, S. A. Khan, S. S. Kristensen, S. B. Andersen, J. E. Box, A. P. Ahlstrøm, X. Fettweis, and R. Forsberg, Basin-scale partitioning of Greenland ice sheet mass balance components (2007–2011), Earth and Planetary Science Letters 409, 89 (2015).
  * Smith  _et al._ [2017] L. C. Smith, K. Yang, L. H. Pitcher, B. T. Overstreet, V. W. Chu, . K. Rennermalm, J. C. Ryan, M. G. Cooper, C. J. Gleason, M. Tedesco, J. Jeyaratnam, D. Van As, M. R. Van Den Broeke, W. J. Van De Berg, B. Noël, P. L. Langen, R. I. Cullather, B. Zhao, M. J. Willis, A. Hubbard, J. E. Box, B. A. Jenner, and A. E. Behar, Direct measurements of meltwater runoff on the Greenland ice sheet surface, Proceedings of the National Academy of Sciences of the United States of America 114, E10622 (2017).
  * Bell [2008] R. Bell, The role of subglacial water in ice-sheet mass balance., Nature Geoscience 1, 297 (2008).
  * van de Wal  _et al._ [2008] R. S. W. van de Wal, W. Boot, M. R. van den Broeke, C. J. P. P. Smeets, C. H. Reijmer, J. J. A. Donker, and J. Oerlemans, Large and Rapid Melt-Induced Velocity Changes in the Ablation Zone of the Greenland Ice Sheet, Science 321, 111 (2008).
  * Chandler  _et al._ [2013] D. M. Chandler, J. L. Wadham, G. P. Lis, T. Cowton, A. Sole, I. Bartholomew, J. Telling, P. Nienow, E. B. Bagshaw, D. Mair, S. Vinen, and A. Hubbard, Evolution of the subglacial drainage system beneath the Greenland Ice Sheet revealed by tracers, Nature Geoscience 6, 195 (2013).
  * Stibal  _et al._ [2012] M. Stibal, M. Šabacká, and J. Žárský, Biological processes on glacier and ice sheet surfaces, Nature Geoscience 5, 771 (2012).
  * Cook  _et al._ [2012] J. M. Cook, A. J. Hodson, A. M. Anesio, E. Hanna, M. Yallop, M. Stibal, J. Telling, and P. Huybrechts, An improved estimate of microbially mediated carbon fluxes from the Greenland ice sheet, Journal of Glaciology 58, 1098 (2012).
  * Tedstone  _et al._ [2020] A. J. Tedstone, J. M. Cook, C. J. Williamson, S. Hofer, J. McCutcheon, T. Irvine-Fynn, T. Gribbin, and M. Tranter, Algal growth and weathering crust state drive variability in western Greenland Ice Sheet ice albedo, Cryosphere 14, 521 (2020).
  * McCutcheon  _et al._ [2021] J. McCutcheon, S. Lutz, C. Williamson, J. M. Cook, A. J. Tedstone, A. Vanderstraeten, S. A. Wilson, A. Stockdale, S. Bonneville, A. M. Anesio, M. L. Yallop, J. B. McQuaid, M. Tranter, and L. G. Benning, Mineral phosphorus drives glacier algal blooms on the Greenland Ice Sheet, Nature Communications 12, 570 (2021).
  * Benning  _et al._ [2014] L. G. Benning, A. M. Anesio, S. Lutz, and M. Tranter, Biological impact on Greenland’s albedo, Nature Geoscience 7, 691 (2014).
  * Tedesco  _et al._ [2016] M. Tedesco, S. Doherty, X. Fettweis, P. Alexander, J. Jeyaratnam, and J. Stroeve, The darkening of the Greenland ice sheet: Trends, drivers, and projections (1981-2100), Cryosphere 10, 477 (2016).
  * Hotaling  _et al._ [2021] S. Hotaling, S. Lutz, R. J. Dial, A. M. Anesio, L. G. Benning, A. G. Fountain, J. L. Kelley, J. McCutcheon, S. M. K. Skiles, N. Takeuchi, and T. L. Hamilton, Biological albedo reduction on ice sheets, glaciers, and snowfields, Earth-Science Reviews 220, 103728 (2021).
  * Schuster [2001] C. J. Schuster,  _Weathering crust processes on melting glacier ice (Alberta, Canada)_ , Ph.D. thesis (2001).
  * Hoffman  _et al._ [2008] M. J. Hoffman, A. G. Fountain, and G. E. Liston, Surface energy balance and melt thresholds over 11 years at Taylor Glacier, Antarctica, Journal of Geophysical Research: Earth Surface 113, F04014 (2008).
  * Hoffman  _et al._ [2014] M. J. Hoffman, A. G. Fountain, and G. E. Liston, Near-surface internal melting: A substantial mass loss on Antarctic Dry Valley glaciers, Journal of Glaciology 60, 361 (2014).
  * Woods and Hewitt [2023] T. Woods and I. J. Hewitt, A model of the weathering crust and microbial activity on an ice-sheet surface, Cryosphere 17, 1967 (2023).
  * Law  _et al._ [2020] R. Law, N. Arnold, C. Benedek, M. Tedesco, A. Banwell, and I. Willis, Over-winter persistence of supraglacial lakes on the Greenland Ice Sheet: Results and insights from a new model, Journal of Glaciology 66, 362 (2020).
  * Perovich [1990] D. K. Perovich, Theoretical estimates of light reflection and transmission by spatially complex and temporally varying sea ice covers, Journal of Geophysical Research 95, 9557 (1990).
  * Liston  _et al._ [1999] G. E. Liston, J.-G. Winther, O. Bruland, H. Elvehøy, and K. Sand, Below-surface ice melt on the coastal Antarctic ice sheet, Journal of Glaciology 45, 273 (1999).
  * Taylor and Feltham [2004] P. D. Taylor and D. L. Feltham, A model of melt pond evolution on sea ice, Journal of Geophysical Research C: Oceans 109, 1 (2004).
  * Taylor and Feltham [2005] P. D. Taylor and D. L. Feltham, Multiple stationary solutions of an irradiated slab, Journal of Crystal Growth 276, 688 (2005).
  * Alexiades and Solomon [1993] V. Alexiades and A. D. Solomon,  _Mathematical Modeling of Melting and Freezing Processes_ (Hemisphere Publishing Corporation, Washington; London, 1993).
  * Meyer and Hewitt [2017] C. R. Meyer and I. J. Hewitt, A continuum model for meltwater flow through compacting snow, Cryosphere 11, 2799 (2017).
  * Voller and Cross [1981] V. Voller and M. Cross, Accurate solutions of moving boundary problems using the enthalpy method, International Journal of Heat and Mass Transfer 24, 545 (1981).
  * Aschwanden  _et al._ [2012] A. Aschwanden, E. Bueler, C. Khroulev, and H. Blatter, An enthalpy formulation for glaciers and ice sheets, Journal of Glaciology 58, 441 (2012).
  * Hewitt and Schoof [2017] I. J. Hewitt and C. Schoof, Models for polythermal ice sheets and glaciers, Cryosphere 11, 541 (2017).
  * Fausto  _et al._ [2021] R. S. Fausto, D. Van As, K. D. Mankoff, B. Vandecrux, M. Citterio, A. P. Ahlstrøm, S. B. Andersen, W. Colgan, N. B. Karlsson, K. K. Kjeldsen, N. J. Korsgaard, S. H. Larsen, S. Nielsen, A. Pedersen, C. L. Shields, A. M. Solgaard, and J. E. Box, Programme for Monitoring of the Greenland Ice Sheet (PROMICE) automatic weather station data, Earth System Science Data 13, 3819 (2021).