A model of the weathering crust and microbial activity on an ice-sheet surface

. Shortwave radiation penetrating beneath an ice-sheet surface can cause internal melting and the formation of a near-surface porous layer known as the weathering crust, a dynamic hydrological system that provides home to impurities and microbial life. We develop a mathematical model, incorporating thermodynamics and population dynamics, for the evolution of such layers. The model accounts for conservation of mass and energy, for internal and surface-absorbed radiation, and for logistic growth of a microbial species mediated by nutrients that are sourced from the melting ice. It also accounts for potential melt–albedo and microbe–albedo feedbacks, through the dependence of the absorption coef-ﬁcient on the porosity or microbial concentration. We investigate one-dimensional steadily melting solutions of the model, which give rise to predictions for the weathering crust depth, water content, melt rate, and microbial abundance, depending on a number of parameters. In particular, we examine how these quantities depend on the forcing energy ﬂuxes, ﬁnding that the relative amounts of shortwave (surface-penetrating) radiation and other heat ﬂuxes are particularly important in determining the structure of the weathering crust. The results explain why weathering crusts form and disappear under different forcing conditions and suggest a range of possible changes in behaviour in response to climate change.


Introduction
The weathering crust on the surface of the south-western Greenland ice sheet affects the transport and storage of meltwater and microbes, as well as the albedo of the ice sheet (Muller and Keeler, 1969;Cooper et al., 2018;Irvine-Fynn et al., 2021).The weathering crust is a layer, around a metre thick, of very porous ice at the surface of the ice sheet, with impermeable ice below (Cooper et al., 2018) (see Fig. 1).The porous structure is formed by shortwave radiation penetrating below the surface, causing "internal" melting.The weathering crust responds to changes in weather conditions.Clear, sunny conditions favour formation of the weathering crust, whereas weathering crust removal typically occurs during periods of cloud cover, warm air temperatures, and strong winds (Muller and Keeler, 1969;Schuster, 2001).Most of the thickness of the weathering crust is saturated with meltwater, forming what is known as the "weathering crust aquifer" (e.g.Christner et al., 2018).
The weathering crust is home to many microorganisms, including dark-coloured glacier ice algae at the surface (following the naming convention recommended by Cook et al., 2019, andHotaling et al., 2021) and algae and bacteria in the weathering crust aquifer (Irvine-Fynn et al., 2021).The growth and transport of these microbes affects the cycling of nutrients such as carbon (e.g.Stibal et al., 2012) and nitrogen (e.g.Telling et al., 2012).For example, microbial activity affects whether an ice-sheet surface is a source or a sink of CO 2 , and microbes and nutrients can be transported by the flow of meltwater, impacting downstream habitats such as oceans and soils (Stibal et al., 2012).
The weathering crust can transiently store meltwater and affect how the water is transported across the ice-sheet surface (Cooper et al., 2018).Moreover, the presence of the weathering crust can affect the amount of meltwater produced, since it modulates how much shortwave radiation is reflected or absorbed by the ice.For example, the weathering crust reflects more radiation than the glazed, bare ice left behind when the crust is removed (Tedstone et al., 2020).
Published by Copernicus Publications on behalf of the European Geosciences Union.(b) Schematic of a weathering crust -a near-surface layer of porous ice (above the dashed line), below which the ice is impermeable (though the "solid" ice shown here may contain small, closed-off air bubbles).The porosity of the weathering crust decreases with depth.Much of it is saturated with meltwater, with a thin unsaturated layer at the surface.The saturated layer is known as the "weathering crust aquifer".Schematic based on Cooper et al. (2018); Moure et al. (2022).
In addition to changes to the ice-sheet surface, biological factors -including glacier ice algae (Benning et al., 2014) contribute to the darkening of the Greenland ice sheet (Hotaling et al., 2021), lowering the albedo and increasing melting.Work has been done to assess and quantify the contribution of glacier ice algae to albedo and surface runoff (e.g.Cook et al., 2017Cook et al., , 2020;;Tedstone et al., 2020;Chevrollier et al., 2022).For example, Cook et al. (2020) estimated that 10 %-13 % of the total runoff from south-western Greenland in summer 2017 could be attributed to algae growth.However, it is difficult to isolate the biological impact from the impact due to changes to the weathering crust structure, as well as other factors including non-biological impurities (Tedesco et al., 2016;Cook et al., 2017;Tedstone et al., 2020).The potential positive feedback loop associated with glacier ice algae and melting, whereby increased melting enables more algal growth, demonstrates the intrinsic link between the structure of the weathering crust and the microbial activity within it (Tedstone et al., 2020;McCutcheon et al., 2021).All of this emphasises the need to understand the microbial community and the ice structure as one.
In this paper, we develop a continuum model for the evolution of the weathering crust, to which we couple a model for microbes and nutrients in the weathering crust aquifer.We reduce our model to one dimension and focus on steadily melting solutions in which constant radiation and atmospheric conditions cause the surface to lower at a constant rate.We investigate how the structure of the weathering crust, the amount of melting, and the amount and distribution of microbes are affected by the forcing conditions.In Sect. 2 we describe the model in the absence of any microbes, and we discuss solutions for the weathering crust structure.In Sect. 3 we extend the model to include microbes and nutrients and discuss solutions for the microbial and nutrient distributions as well as the influence of radiative forcings.In Sect.4, we use our model to investigate feedbacks between microbial growth, albedo, and melting.Finally, in Sect. 5 we discuss our results and conclude.
2 Weathering crust model

Model
Firstly, we present a model for the evolution of the structure of the weathering crust without any microbes or nutrients present.We use a continuum model, in contrast to the layer model used by Schuster (2001) -the only existing attempt to model the weathering crust that we are aware of.We model the weathering crust as porous ice with upper surface z = h(t).The pores can contain both water and air.The ice takes up a volume fraction 1 − φ, where φ is the porosity.The saturation S describes the proportion of the pore space taken up by water, so φS is the volume fraction of water and φ(1 − S) is the volume fraction of air.The porosity and saturation can evolve by the ice melting internally at a rate m int (with negative m int corresponding to water freezing in the pores) and by the water flowing through the porous ice.The internal melting is driven by shortwave (solar) radiation penetrating below the ice surface.We now present the governing equations.

Mass conservation
Mass conservation equations for the ice, water, and air, respectively, are where ρ i , ρ w , and ρ a are the densities of ice, water, and air, respectively, and u i , u w , and u a are their velocities.Since air is much less dense than ice and water (ρ a ρ i , ρ w ), we can neglect Eq. (3).A relationship between the ice and water velocities could come from Darcy's law for flow through a porous medium, and u i could be prescribed to account for the upward advection of ice in the ablation zone.However, for our purposes in this study, we will later set both u i and u w to zero for the specific one-dimensional solution we consider.

Energy equation
We also need to model the temperature T , which we assume to be homogeneous across the ice-water-air mixture.We assume that if water is present (φ > 0, S > 0) then the region is at the melting temperature T m .If there is no water present (S = 0), then the ice can be below the melting temperature T ≤ T m .Conservation of energy leads to the following equation: The terms on the left-hand side represent the rate of change and advection of heat energy.The terms on the right-hand side correspond to heat conduction according to Fourier's law, the heat source provided by internal shortwave radiation F (discussed further below), and the energy available for melting the ice.Here, c i and c w are the specific heat capacities and k i and k w are the thermal conductivities for ice and water, respectively, and L is the latent heat of fusion of water.
Consistent with assuming negligible density of air, we ignore both the heat content and the conductive heat flux in the air.
We have assumed a simple volume averaging for the effective thermal conductivity of the combined ice-water mixture.
The energy equation captures two cases.Firstly, when the system is below the melting temperature (T < T m ), there cannot be any internal melting (m int = 0) and there is no water present (S = 0).The energy Eq. ( 4) determines the temperature.Alternatively, when the system is at the melting temperature (T = T m ), melting is possible (m int ≥ 0) and water can be present (S ≥ 0).Since we know that T = T m , the energy Eq. ( 4) instead determines the internal melt rate m int .

Radiation model
The internal melting of the ice is caused by some of the shortwave radiation from the sun penetrating and being absorbed below the surface of the ice, acting as an internal heat source.We model the internal shortwave radiation flux F following the two-stream approximation (e.g.Perovich, 1990;Liston et al., 1999;Taylor andFeltham, 2004, 2005).We assume that the net radiation flux can be decomposed into upward F + and downward F − travelling components.The net radiation flux is then F = (F + −F − ) k, where k is the upward-pointing unit vector, with F ± satisfying where α and r are absorption and scattering coefficients, respectively, with α, r > 0. We assume that α and r represent bulk properties of the ice-water mixture and, for now, we also assume that α and r are constant.These equations tell us that, as we move through the ice, F + is decreased by some of this radiation being absorbed and reflected but increased by some of F − being reflected (and thus changing from a downward to an upward flux).A similar argument applies for F − .

Kinematic conditions
The kinematic conditions for the ice and the water at the surface z = h(t) are, respectively, where u i • k and u w • k are the upward ice and water velocities (which, in the case of water, we expect to be negative), ḣ = dh/dt is the velocity of the surface, m surf is the rate of surface melting, and r off is the rate of runoff.We assume that any water at the surface drains down into the ice below, provided that the ice is unsaturated.If the surface becomes saturated, we assume that any excess water immediately runs off.In reality, this must be assumed to happen via horizontal flow to other parts of the ice surface and ultimately into the ocean.The runoff term r off is therefore only non-zero when the surface is saturated.

Surface energy balance
The energy balance at the surface z = h(t) is (see, for example, Hoffman et al., 2008;van den Broeke et al., 2008van den Broeke et al., , 2011) where the terms in order from left to right represent conduction of heat from the porous ice, incoming shortwave radiation absorbed at the surface (Q si is the incoming shortwave (solar) radiation, a is the albedo, χ is the proportion of the total shortwave radiation absorbed by the ice which is absorbed at the surface, discussed further below), emitted longwave radiation given by the Stefan-Boltzmann law for a black body https://doi.org/10.5194/tc-17-1967-2023 The Cryosphere, 17, 1967Cryosphere, 17, -1987Cryosphere, 17, , 2023 ( is the surface emissivity, a correction for the fact that the ice is not a true black body, and σ is the Stefan-Boltzmann constant), incoming longwave radiation (for example, from reflection off clouds), turbulent heat transfer (C is a turbulent heat transfer coefficient, and T a is the reference air temperature 2 m above the ice surface), and latent heat flux available for surface melting at rate m surf .The turbulent heat flux term takes into account the fact that, if the air is warmer than the surface, winds can cause turbulent eddies that transfer heat from the air to the surface.For simplicity, we are neglecting latent heat transfer due to evaporation and ablation since van den Broeke et al. (2008Broeke et al. ( , 2011) ) showed this contribution to be small compared to the other terms, at least in the ablation zone of the western Greenland ice sheet.
Often, it is assumed that shortwave radiation incident on an ice-sheet surface is either absorbed at the surface or reflected straight back out, with no radiation being transmitted below the surface.This assumption results in a (1 − a)Q si term in the surface energy balance.However, a key part of weathering crust formation is the penetration of shortwave radiation below the ice surface.Hence, based on methods used by Hoffman et al. (2014) and Law et al. (2020), we introduce the parameter χ < 1 -the proportion of the total amount of shortwave radiation absorbed by the ice which is absorbed at the surface -to allow some of the shortwave radiation to penetrate into the subsurface ice, where we model it explicitly as the internal shortwave radiation flux F (see discussion around Eqs. 5 and 6).In particular, this means that the downward internal radiation flux F − at the surface Note that in this model the albedo a is not an independent parameter.It is a result of the internal radiation calculation, dependent on the scattering and absorption coefficients α and r (see Sect. 2.3).
By linearising about the melting temperature, the surface energy balance can be re-written as where ) is a combination of longwave radiation and turbulent heat flux, which we will prescribe as a function of time, and υ = C + 4 σ T 3 m is an effective heat transfer coefficient, combining turbulence and radiation.

Model summary
The full weathering crust model is given by the mass conservation Eqs. ( 1) and (2) for ice and water, the energy Eq. ( 4), and the two-stream Eqs. ( 5) and ( 6) for the internal shortwave radiation, along with a prescription of the ice and water velocities u i and u w and suitable boundary conditions.As well as the kinematic conditions in Eqs. ( 7) and (8), the surface energy balance Eq. ( 10), and the condition for downward-travelling radiation F − at the surface, we also require a condition on the upward-travelling radiation F + and the temperature T in the far-field (discussed below).The incoming shortwave radiation Q si and the "other" surface radiation Q 0 (a combination of longwave radiation and turbulent heat fluxes) are prescribed functions of time, representing the climate forcing.In addition to the temperature T , porosity φ, and saturation S, additional outputs from the model are the surface and internal melt rates, m surf and m int , and surface runoff, r off .The default parameter values used in the model are shown in Table 1.

Steadily melting states
We now make some simplifications to our model.Firstly, we make the approximation that the ice and water have the same properties, introducing a shared density ρ = ρ i = ρ w , thermal conductivity k = k i = k w , and specific heat capacity c = c i = c w .Further, we reduce our model to one dimension considering a vertical column of ice.We also assume that the porous ice is saturated (S = 1) and that there is no flow of ice or water (u i = u w = 0).Some of these simplifications could quite easily be relaxed but implementing them allows for a more transparent interpretation of the model solutions.
We now seek what we will refer to as "steadily melting" solutions.These are solutions in which the surface is lowering at a constant rate (− ḣ), with the subsurface ice melting at exactly a rate m int to maintain a steady porosity profile below the surface.Therefore, relative to the surface, the solution looks the same at each snapshot in time, as if the column of porous ice is simply being translated downwards, as shown in Fig. 2a.However, this effect is produced by internal and surface melting, not by the ice moving (recall that the ice velocity is zero, by assumption; in reality the downward melting of the ice surface may be counteracted by the upward advection of ice from the interior of the ice sheet).The steadily melting states correspond to an equilibrium between internal melting and surface lowering, which Schuster (2001) observed can occur on cloudy days when the crust is transitioning between growth and decay.It is not expected that this equilibrium lasts long; the weathering crust is very changeable, growing and decaying over a few hours as weather conditions change (Schuster, 2001).Therefore, it would be preferable to have a time-dependent model that can capture this variation.However, for simplicity and to gain some initial insight, we chose to use steady forcings.
The steadily melting setup suggests a natural change of coordinates to where Z is the depth below the surface z = h(t) (see Fig. 2).Note that Z points downwards.In the Z frame, the melting surface of the ice is fixed at Z = 0, and the ice below gets advected upwards towards the surface at rate − ḣ, as shown   The Cryosphere, 17, 1967Cryosphere, 17, -1987Cryosphere, 17, , 2023 in Fig. 2b.Making this transformation and using the simplifications previously discussed, we can turn our general weathering crust equations into equations for the steadily melting solution.The two-stream Eqs. ( 5) and ( 6) for the internal radiation Both the ice and water mass conservation Eqs. ( 1) and ( 2) reduce to and the energy Eq. ( 4) becomes The advection terms ḣd/dZ come from the time derivatives under the transformation of z → Z.As for boundary conditions, the kinematic conditions ( 7) and ( 8) tell us that and that the runoff rate is The linearised surface energy balance (10) becomes where we have used the fact that T = T m at the surface in order for surface melting to be possible.In addition to these, we also prescribe a cold temperature T ∞ < T m at depth, so Finally, the boundary conditions for the two-stream Eqs. ( 12) and ( 13) become and since we assume there is no source of shortwave radiation at depth.We also have which serves to determine the albedo a (which is defined as the overall ratio of the reflected and incident shortwave radiation).Note that in some cases it may be more natural to prescribe the albedo a as a known (i.e.measured) quantity.
In this model, the albedo is related to χ , α, and r, so independently prescribing a would mean that one of these other parameters must be prescribed consistently (see Eq. 24 below).
There are three forcings in our model, which we prescribe.These are the incoming shortwave radiation Q si , the combined longwave radiation and turbulent heat fluxes Q 0 , and the temperature of the ice at depth T ∞ .The outputs of the model are the internal shortwave radiation F , porosity φ and temperature T , the internal melt rate m int , surface melt rate m surf , runoff r off , and rate of change of surface height − ḣ.We choose values of the forcings such that the surface is melting (m surf > 0) and the deep ice is cold.Therefore, as shown in Fig. 2, we expect there to be a temperate region 0 < Z < Z m near the surface where the ice is porous (φ > 0) and at the melting temperature (T = T m ).This region is the weathering crust.In Z > Z m , the ice is solid (φ = 0) and cold (T < T m ).We refer to Z m as the weathering crust thickness or the melting depth, which must also be determined by the model.Note that having impermeable, solid ice at Z = Z m , combined with our assumption that the ice is saturated with water, justifies our assumption of no water flow (u w = 0), at least in the vertical direction.Furthermore, Irvine-Fynn et al. ( 2021) have observed lateral flow through the weathering crust to be slow (less than about 0.13 m d −1 ), which somewhat justifies neglecting horizontal water velocities and using a one-dimensional model.

Results
In the case where we have constant absorption coefficient α, scattering coefficient r, and albedo a, we can find the steadily melting solution analytically.Firstly, the two-stream Eqs. ( 12) and ( 13) can be solved with boundary conditions ( 20) and ( 21) to find that the internal radiation F = F − − F + decreases exponentially with depth, where the extinction coefficient κ is related to the absorption and scattering coefficients α and r by κ = √ α 2 + 2αr and s = (κ − α)/(κ + α) (e.g.Taylor andFeltham, 2004, 2005).The extra boundary condition ( 22) for F + tells us that This is the relationship telling us that we can only prescribe three of α, r, a, and χ (or equivalently, three of κ, s, a, and χ).Using the condition (24), we can eliminate s from F + to obtain This is the familiar Beer's law, previously used to model radiation attenuation in sea ice (e.g.Maykut and Untersteiner, 1971) and glaciers (e.g.Kelliher et al., 1996).We chose to use the two-stream approximation to model the internal radiation rather than assuming Beer's law straight away so that our model can later be extended to non-constant absorption and scattering coefficients α and r.We also note that Cooper et al. (2021) recently made observations of the spectral attenuation of light through bare ice on the Greenland ice sheet, observing that attenuation is greater than through pure ice due to the presence of impurities.However, since our model considers a single band of wavelengths of radiation, we choose to follow Taylor and Feltham (2005) in using κ = 1.5 m −1 .This is towards the lower end of the range of values observed by Cooper et al. (2021) (∼ 0.9-8 m −1 for 350-750 nm), but it is at least in the right ballpark.
We can now analytically solve for the steadily melting temperature T , porosity φ, and internal melt rate m int .As mentioned previously, we seek solutions with a temperate region 0 < Z < Z m above a cold region Z > Z m .The boundary conditions between the two regions are where the final condition comes from continuity of conductive heat flux.Also, since T = T m everywhere in 0 < Z < Z m , there is no conductive heat flux at the surface, so the surface energy balance Eq. ( 18) directly tells us the surface melt rate: In the cold region, we integrate the temperature Eq. ( 4) and apply the far-field condition (19), along with the fact that a uniform far-field temperature means that the heat flux tends to zero, to obtain for Z > Z m .Evaluating this at Z = Z m and using the boundary conditions (26) gives us an expression which can be rearranged to tell us Z m in terms of ḣ: The expression ( 28) can be integrated once more using the boundary condition ( 26) for T to find that the temperature is In the temperate region, the energy Eq. ( 15) determines the internal melt rate m int from the internal radiation F as follows: In turn, φ can be calculated from m int using the mass conservation Eq. ( 14), giving (32) Finally, the kinematic condition ( 16) determines the rate of change of surface height ḣ as follows: which can be substituted into the expression (29) for the melting depth to give As is automatic from Eq. ( 17), the surface runoff in this situation is equal to the rate of surface lowering, r off = − ḣ.Example solutions for temperature, porosity, internal radiation, and internal melt rate are shown in Fig. 3 for two different values of the deep-ice temperature T ∞ .Unsurprisingly, a colder deep-ice temperature (blue lines) leads to a smaller porous region, with the dotted lines indicating Z = Z m .
For our assumption that the surface is melting to be valid, we need the surface melt rate m surf to be positive.Using the expression ( 27) for m surf , this gives a restriction on the possible values of the prescribed radiations Q si and Q 0 , Furthermore, the solution found here, with a porous ice region occupying 0 < Z < Z m , is only valid if the calculated Z m is positive, and hence the argument of the logarithm in Eq. ( 34) must be greater than 1.Physically, this can be interpreted as requiring that the raising of the internal temperature dominates the surface lowering, as we will now explain.Firstly, as more of the internal ice is raised to the melting temperature, the size of the porous region (and hence Z m ) increases.This is represented by the term in the expression (34) for Z m .The numerator is the amount of radiation absorbed internally, i.e. the energy available to heat and melt the internal ice.The denominator represents the energy needed to raise the ice temperature from T ∞ at depth to the melting temperature, i.e. the energy required to raise the internal ice to the melting temperature.Therefore, the fraction (36) can be thought of as a measure of raising the internal ice to the melting temperature.
On the other hand, as the surface lowers, the size of the porous region (and hence Z m ) decreases, since the porous ice at the surface gets removed as meltwater.Surface lowering is represented by the term  in the expression (34) for Z m .The numerator is the total amount of radiation absorbed by the ice, accounting for both internal and surface absorption.This is the total energy available for melting the ice (both at the surface and internally).
The denominator represents the energy required to melt cold ice; the temperature must first be raised from T ∞ to T m , then latent heat is needed to effect the phase change.Indeed, it is this fraction that also gives the surface lowering rate − ḣ in Eq. ( 33).It is different from the surface melt rate m surf , since melting occurs in two stages: first the ice can melt internally, lowering the porosity, then the remaining ice melts at the surface.
As previously mentioned, we expect Z m > 0 when raising the internal temperature dominates surface lowering.Using the above interpretations of the fractions (36) and (37), this happens when and Eq. ( 34) then gives Z m > 0.
The conditions ( 35) and ( 38) tell us the range of parameter values for which we can have steadily melting solutions with a non-trivial porous region, i.e. for which a weathering crust exists.Figure 4a illustrates this graphically, showing the range of values of shortwave radiation Q si and combined longwave radiation and turbulent heat fluxes Q 0 for which the steadily melting states with a non-trivial weathering crust are possible (blue to yellow region).For a fixed amount of shortwave radiation Q si , increasing the other surface radiation Q 0 increases the amount of surface lowering but not the internal heating.This is because Q 0 is completely absorbed at the surface.If Q 0 is too large, surface lowering is so large that ice gets advected to the surface so quickly that it does not have the chance to get heated up to the melting temperature internally first.Hence there is no internal melting.If Q 0 is too small, there is not enough energy to melt the surface.
The weathering crust thicknesses Z m found from our model broadly agree with the ice core measurements taken by Cooper et al. (2018) in the south-western Greenland ice sheet ablation zone.They found the weathering crust in that region to be at least 1 m and up to about 2 m thick, in line with our results in Fig. 4b.Moreover, Fig. 4b and c show how the weathering crust thickness Z m depends on the incoming shortwave radiation Q si and the combined longwave radiation and turbulent heat fluxes Q 0 .There is a qualitative difference between the behaviour of Z m with Q si , depending on the sign of Q 0 .We see that Z m increases with Q si when Q 0 > 0 but decreases with Q si when Q 0 < 0. The reason for this relates back to our interpretation of Z m being determined by a balance between internal heating and surface lowering.If Q 0 < 0, doubling Q si doubles the internal heating term (36) and more than doubles the surface lowering term (37).Therefore, as Q si increases, surface lowering becomes more significant compared to internal heating, so the porous region (and hence Z m ) decreases.On the other hand, if Q 0 > 0, surface lowering becomes less significant as Q si increases, so the porous region (and hence Z m ) increases.In both cases, a weathering crust can only exist if the shortwave radiation Q si is large enough.Physically, positive Q 0 corresponds to conditions of warm air temperature, strong wind, and cloud cover, when the ice surface would melt even in the absence of solar radiation.Such conditions favour removal of the weathering crust (Muller and Keeler, 1969;Schuster, 2001), as we observe in Fig. 4 when Q 0 becomes too large.Negative Q 0 corresponds to conditions when the surface would freeze in the absence of solar radiation.In the ablation zone of the western Greenland ice sheet, we would typically expect the net longwave radiation to be negative and the turbulent heat fluxes to be positive, with the combination giving Q 0 as slightly negative (van den Broeke et al., 2008Broeke et al., , 2011).An exception is in regions near the margin where ice-free tundra causes positive turbulent heat fluxes to dominate (van den Broeke et al., The other parameter values are those shown in Table 1.2008,2011).However, the weather conditions on the Greenland ice sheet are very changeable (e.g.Fausto et al., 2021), rapidly changing between clear, sunny conditions that favour weathering crust formation and cloudy, windy, and warm weather that favours removal of the crust (Muller and Keeler, 1969;Schuster, 2001).Therefore, we choose to investigate both positive and negative Q 0 to explore the different behaviours exhibited during these changeable conditions.The sensitivity of the weathering crust thickness Z m to the values of Q si and Q 0 , demonstrated by our model might partly explain why the weathering crust tends to form and be removed by changing weather conditions.
3 Microbes and nutrients model

Model
We now extend our model to include microbes and nutrients.Our model focusses on the microbes within the weathering crust aquifer and does not account for the ice algae found right on the ice-sheet surface (see Sect. 5 for further discussion).In order to keep things simple, we model a single "representative" nutrient, the concentration of which influences the microbial growth.In reality, there may be multiple sources of such nutrients located within the modelled space.The microbial growth will also depend on the shortwave radiation flux F .We also model a single representative type of phototrophic microbe (a microbe that obtains energy from solar radiation), rather than separately considering different primary and secondary producers (see discussion in Sect.5).We assume that both the microbes and the nutrients can exist in both the ice and the water, and when the ice melts they are transferred from the ice to the water.Microbes in the water can consume nutrients from the water but not from the ice, and microbes in the ice cannot consume any nutrients (see Fig. 5).Therefore, only microbes in the water are able to proliferate.We model the microbial abundance of microbes in the water A w and ice A i (cells mL −1 ) separately, and we do the same for the nutrient concentrations C w and C i (µmol L −1 ).
To model the water microbial abundance A w , we adapt a simple logistic growth model (Murray, 1993), https://doi.org/10.5194/tc-17-1967-2023 The Cryosphere, 17, 1967Cryosphere, 17, -1987Cryosphere, 17, , 2023 T. Woods and I. J. Hewitt: A model of the weathering crust and microbial activity on an ice-sheet surface On the left-hand side, we have the time derivative and an advection term, with the microbes being advected at the water velocity u w .On the right-hand side we have a growth term, a death term (with a death rate d w ), and source and sink terms due to melting and freezing, respectively.The growth term, of the form gφA w (1 − A w /A max ), is the standard form for logistic growth, with A max being the maximum microbial abundance.We base the form of the logistic growth rate g on that used by Jørgensen and Bendoricchio (2001) for algal growth, in which g depends on the intensity of photosynthetically active radiation, the volume-averaged nutrient concentration φC w and the temperature, since these are limiting factors for growth of photosynthetic microbes.(We can ignore the temperature dependence, since we assume that the water is always at the melting temperature.)We follow Jørgensen and Bendoricchio (2001) in using a Michaelis-Menten form of the growth term, where α PAR = 0.56 is the proportion of the total incident solar radiation that is photosynthetically active (i.e. that can be used for photosynthesis).Photosynthesis, like weathering crust formation, is driven by shortwave radiation.Here, the broadband radiation F represents all wavelengths of shortwave radiation.The factor α PAR accounts for the fact that the two processes in fact use slightly different (overlapping) parts of the spectrum.k PAR and k C are constants used to represent the fact that, above a certain limit, increasing the radiation or nutrient concentration further will no longer lead to increased algal growth.β A is also a constant, which represents the growth rate if there is abundant radiation and nutrients (in which case the two bracketed factors are essentially 1).We refer to β A as the growth rate constant.The ice microbial abundance equation is simpler.We assume that the microbes in the ice cannot grow but can die (with a death rate d i ) and are fixed within the ice (so they are advected at the ice velocity u i ).Hence A i satisfies We expect that d i ≥ d w , meaning that the death rate of microbes in the ice is greater than the death rate in the water.Nutrients are transported through the weathering crust aquifer, often from nutrient-rich cryoconite holes (Cook et al., 2016).To model the nutrient transport in the water, we use similar ideas to those used to model bioreactors (e.g.Shakeel et al., 2013).The nutrients in the water can be advected and consumed by microbes.The melting of ice provides a source of nutrients into the water (and a sink from the ice), and water freezing provides a corresponding sink.The nutrients in the ice are advected with the ice but cannot be consumed because microbes cannot access the nutrients in the ice.The equations describing this are for the nutrients in the water and for the nutrients in the ice.The uptake term has a similar form to the growth rate in Eq. ( 39) since we expect the consumption of nutrients to be related to the growth rate of the microbes.The constant β C represents a maximum nutrient uptake rate and controls the degree to which microbes influence the nutrient concentration.
A further assumption we make is that microbes and nutrients get instantaneously removed in surface runoff.This would need to be modified if we extended our model to include surface ice algae, which are able to persist on the ice surface, even in the presence of surface runoff (see Sect. 5).

Steadily melting equations
We now return to the steadily melting states considered in Sect.2.2.Under these same assumptions and making the same change of coordinates (11), Eqs. ( 39) and (42) for the microbial abundance and nutrient concentration in the water become As an additional simplification, we make the assumption that the microbes in the water and in the ice do not die, so d w = d i = 0.In this particular situation, Eqs. ( 41) and ( 43) for A i and C i simply tell us that the microbial abundance and nutrient concentration in the ice are constant.Since we prescribe the microbial abundance and nutrient concentration in the deep ice to be A ∞ and C ∞ , respectively, this leads to Referring back to Fig. 2, water is only present in the porous region 0 < Z < Z m , so it is in this region that we solve Eqs. ( 44) and ( 45) for A w and C w .For continuity at the interface Z = Z m between the porous ice and the solid ice, we prescribe For simplicity, we assume for now that absorption of shortwave radiation is independent of microbial abundance (with this assumption being removed in Sect.4).As a result, the microbe and nutrient problem decouples from the weathering crust problem.The analytical solution to the weathering crust problem found in Sect.2.3 can first be used to calculate φ, m int , ḣ, and Z m .Then, using these results, Eqs. ( 44) and ( 45) can be solved (numerically) with the boundary conditions ( 46) and ( 47) for A w and C w .We do this using the boundary value problem solver bvp4c in MATLAB.Our code is made available at https://doi.org/10.5281/zenodo.7199159(Woods, 2022).Before looking at some solutions, we discuss what we expect to be typical values of the microbial abundance and nutrient concentrations.Motivated by Christner et al. ( 2018), we choose A max = 10 4 cells mL −1 and A ∞ = 10 2 cells mL −1 as sensible values for the maximum microbial abundance and the microbial abundance deep in the ice.The nutrient concentration deep in the ice C ∞ gives the maximum value of C w in the steadily melting model.To estimate a value for C ∞ , we note that Holland et al. ( 2019) measured concentrations of dissolved organic and inorganic nitrogen, phosphorous, and carbon near the surface of the Greenland ice sheet dark zone varying in orders of magnitude from about 0.1 µmol L −1 to over 100 µmol L −1 .The nutrient concentration of interest will be that of the most limiting nutrient, which could potentially be any of these.Therefore, as a ballpark figure, we take C ∞ = 1 µmol L −1 .In fact, the actual numerical values of nutrient concentration are not particularly important to the model, and we will show the nutrient concentration solutions normalised by C ∞ .

Effect of microbial growth rate
Figure 6 shows an example solution of how the volumeaveraged microbial abundance φA w and nutrient concentration φC w vary with depth in the weathering crust shown in Fig. 3.The solutions are shown for two different values of microbial growth rate: β A = 20 d −1 (dark blue line) and β A = 5 d −1 (grey line).The latter value is loosely based on observed doubling times of algal blooms on the Greenland ice sheet (Stibal et al., 2017;Williamson et al., 2018).We note that these observations are of lateral spreading of algal blooms on the ice-sheet surface, not microbial growth within the weathering crust.Therefore, our approximate values of β A are crude but still serve to demonstrate the key behaviour it is possible to capture with our model.Even in the porous region (above the dotted line showing where Z = Z m ), there is an inner region where the microbial abundance A w and nutrient concentration C w are still approximately their deep-ice values, and hence φA w /A max ≈ φA ∞ /A max (= 0.01φ) and φC w /C ∞ ≈ φ.Above this, A w increases towards the surface and C w decreases.This is because there is more shortwave radiation nearer the surface, enabling more microbial growth.In turn, more microbes mean that more nutrients get consumed.We can see that, as we would expect, there are more microbes when the growth rate is higher (dark blue line).The turning point in the nutrient profiles is the point above which the nutrient concentration C w changes significantly from the deep-ice value C ∞ , causing the bulk nutrient concentration φC w to move off the exponential porosity profile.

Effect of nutrient parameters
We now investigate how the distribution of microbes is affected by parameters related to the uptake of the nutrient.In particular, we look at the parameters β C and k C .Firstly, β C is the maximum nutrient uptake rate, with a larger value of β C meaning that the microbes consume more nutrients to achieve the same amount of growth.Figure 7a shows the distribution of microbial abundance A w and nutrient concentration C w with depth for two values of β C -one 100 times larger (light blue) than the other (yellow).Consistent with our interpretation of β C , the nutrient concentration C w is lower when β C is larger.Correspondingly, the microbial abundance A w is lower because the nutrient has been used up more quickly by the microbes, so less nutrient is available for further microbial growth.
microbes.The larger k C is, the higher the volume-averaged nutrient concentration φC w must be before the microbial growth rate becomes independent of φC w .Equivalently, the smaller k C is, the lower the volume-averaged nutrient concentration must get before the nutrient becomes a limiting factor.Figure 7b shows the distribution of microbial abundance A w and nutrient concentration C w with depth for two values of k C -one 10 times larger (dark blue) than the other (green).The effect of decreasing k C is to increase the size of the region in which there is a significant number of microbes.This is because the porosity decreases with depth, which in turn means the volume-averaged nutrient concentration φC w decreases with depth in the region near Z = Z m where C w ≈ C ∞ .Therefore, as k C becomes smaller, the depth at which the nutrient becomes a limiting factor (controlled by the size of φC w relative to k C ) increases, so the microbes spread deeper into the ice.Note that decreasing β C does not affect the depth to significant microbial growth because the size of β C does not change the conditions required for microbial growth -it only affects the rate at which nutrients are consumed.

Total microbial abundance
As we did for the weathering crust, we can examine how the steadily melting solution varies with radiative forcing, assuming fixed microbe and nutrient parameters.Firstly, we show solutions for different Q 0 and Q si in Fig. 8. Secondly, we look at the dependence of the total amount of microbes in the water A tot,w , defined as This has units of cells m −2 .As discussed at the end of Sect.2.3, changing weather conditions (represented by the radiative forcings Q si and Q 0 in our model) lead to the formation and removal of the weathering crust.Therefore, we choose to investigate the response of microbes to changes in both Q si and Q 0 , including positive and negative Q 0 , to capture the range of conditions experienced by microbes in the weathering crust.For a weathering crust to form, the shortwave radiation Q si must be large enough.Positive Q 0 corresponds to conditions of high turbulent heat flux (warm, cloudy, windy), which promote removal of the weathering crust.Negative Q 0 means that the surface would freeze in the absence of shortwave radiation.The observed behaviour is quite different in these two cases.
Figure 9a shows that when Q 0 > 0 and for small values of Q si -below the critical value required to produce internal melting -there are no water microbes because there is no water.Then, as Q si increases, A tot,w increases due to the increase in internal radiation leading to an increase in microbial growth and an increase in the size of the weathering crust Z m (i.e. the region in which microbes can grow).As shown in Fig. 9c, if Q si is increased further (to less physical values), the total amount of microbes A tot,w starts to decrease again.This is because the internal radiation reaches the level above which it is no longer a limiting factor for the microbes.Increasing the shortwave radiation further does not increase the microbial growth but does continue to increase the surface lowering and runoff, meaning that the removal of microbes in runoff begins to dominate microbial growth, causing a reduction in A tot,w .
However, when Q 0 < 0, the total number of microbes A tot,w decreases monotonically with increasing Q si (compared to increasing when Q 0 > 0).The reason for this relates back to Fig. 4. As previously discussed, when Q 0 < 0, the melting depth Z m decreases as Q si increases because sur-  shows contours of A tot,w in the region for which steadily melting states with both internal and surface melting are possible.The contour spacing is 0.02 × 10 10 cells m −2 .The regions for which there is no surface melting and no internal melting are also shown (grey, see Fig. 4).The other parameter values are those shown in Table 1. https://doi.org/10.5194/tc-17-1967-2023 The Cryosphere, 17, 1967Cryosphere, 17, -1987Cryosphere, 17, , 2023 T. Woods and I. J. Hewitt: A model of the weathering crust and microbial activity on an ice-sheet surface face lowering increases at a faster rate than internal melting.Therefore, the region in which microbes can grow gets smaller, with the microbe distribution shifting upwards, as shown in Fig. 8a. Figure 9b also shows how A tot,w depends on the other surface radiation Q 0 , which is a combination of net longwave radiation and turbulent heat fluxes, for a fixed value of shortwave radiation Q si .We see that A tot,w is a decreasing function of Q 0 .This is because increasing Q 0 increases the surface melting without affecting the internal melting or microbial growth.Hence, the surface runoff increases without microbial growth increasing to compensate.When Q 0 gets large enough, there is no time for any microbial growth before the microbes get removed in runoff, so the microbial abundance remains close to the small background value A ∞ .
The combined effect of Q si and Q 0 on A tot,w is shown in the contour plot in Fig. 9f.In addition to the behaviour captured by Fig. 9a-c, this shows how increasing Q 0 to more positive values causes the maximum of the blue A tot,w against Q si curve in Fig. 9c to decrease and move to a higher value of Q si .As discussed, the total number of microbes A tot,w is determined by a balance between the microbial growth in the interior and removal of microbes at the surface in runoff, with the position of the maximum being controlled by the point past which increasing Q si further no longer increases the microbial growth significantly.When Q 0 is more positive, this point is delayed to larger Q si because the size of the region in which microbial growth is possible, Z m , takes longer to tend to its maximum value (see Fig. 4).Hence, for more positive Q 0 , Qsi must get larger before A tot,w starts decreasing with Q si .
Importantly, the key result from Fig. 9f is that the total number of microbes in the water is largest when Q 0 is negative and Q si is just large enough to allow surface melting, which are the conditions under which weathering crust growth is most favoured.

Residence times
Another quantity of interest is the residence time of the microbes -that is, how long the microbes remain in the weathering crust before being removed in runoff.This is of interest because it affects the carbon fluxes from the ice sheetin particular whether or not the ice sheet is a net source or sink of CO 2 .Cook et al. (2012) showed the Greenland ice sheet to be in a stable state of autotrophy (net removal of CO 2 from the atmosphere) due to the presence of cryoconite and surface algae, with increased warming leading to increased biomass and carbon fixation.However, they warned that too much warming could lead to the removal of microbes in runoff, reducing the removal of CO 2 from the atmosphere.
In our steadily melting model, we approximate the average residence time by where is the microbial runoff rate at the surface.Equation (49) for t res is the ratio between the amount of microbes in the water of the weathering crust and the rate at which the microbes are removed from the weathering crust.Therefore, we can think of t res as being the average amount of time it would take to remove all the microbes from the weathering crust water.We can also interpret this as the average amount of time each microbe spends in the weathering crust, i.e. the residence time.Figure 9d shows the dependence of t res on the incoming shortwave radiation Q si .This is shown for positive (dark blue line) and negative (light blue line) values of Q 0 .
As with the total number of microbes A tot,w , the observed behaviour is quite different depending on the sign of Q 0 .When Q 0 > 0, the global behaviour (ignoring the local minimum in Fig. 9d which we will discuss shortly) is that as Q si increases, t res initially increases then decreases again, for the same reasons as those just discussed around A tot,w : as Q si initially increases, the microbial growth rate increases so there are more microbes in the water, and these have to wait longer until they can be removed.Once Q si is large enough to no longer be a limiting factor, the removal of microbes at the surface starts to dominate.This agrees with the previously discussed behaviour expected by Cook et al. (2012): increased warming leads to increased biomass (and residence times) until the microbes get removed by runoff.The local minimum observed in Fig. 9d occurs because the microbial abundance in the ice is small, A i /A max = A ∞ /A max = 0.01.As Q si increases, the microbial abundance A w initially does not increase very much from A ∞ (seen from the flat section in Fig. 9a for small Q si ).When significant microbial growth does begin, growth first occurs at the surface, with the surface microbial abundance becoming quite large before the region of significant growth spreads downwards.When nearsurface growth dominates, the microbial runoff rate q (which depends on the surface value of the microbial abundance) grows much more rapidly with Q si than the total number of microbes A tot,w because the microbial abundance is only growing significantly near the surface.Hence, the residence time t res decreases.Once significant microbial growth begins spreading into the interior, A tot,w grows more quickly than q, so t res increases again.
However, when Q 0 < 0, t res decreases monotonically with increasing Q si (as does A tot,w ).The residence time decreases for the same reason that it decreases for large Q si when Q 0 > 0 -increasing Q si increases the surface lowering but does not increase the growth rate because radiation is high enough that it is no longer limiting.The non-monotonic behaviour for smaller Q si is missed because Q si cannot take such small values when Q 0 < 0, otherwise surface melting, and hence steadily melting states, would not be possible (see Fig. 4 and surrounding discussion in Sect.2.3).
Figure 9e also shows that the approximate residence time t res is a decreasing function of Q 0 for a fixed value of Q si .This is for the same reasons (already discussed) that A tot,w is a decreasing function of Q 0 for a fixed value of Q si .

Summary
Our steadily melting state model has shown that microbial abundance in the weathering crust generally decreases with depth, as expected by Cook et al. (2016) and observed by Christner et al. (2018).This is because the radiation -which the microbes need to grow -decreases with depth.We have seen that significant microbial growth often only occurs in a top portion of the weathering crust, with the size of this region being controlled by microbial growth and nutrient parameters.
Furthermore, we have seen that the radiative forcings Q si and Q 0 have a significant impact on the total number of microbes in the weathering crust and their residence time.This impact is via the balance between internal radiation (which favours growth of microbes and the weathering crust) and surface radiation (which increases surface melting and removal of microbes via runoff).In general, increasing internal radiation at a greater rate than surface radiation leads to an increase in the number of microbes and their residence time.In line with this, we observe different behaviour depending on the size of the shortwave radiation Q si and the sign of the combined longwave radiation and turbulent heat fluxes Q 0 , suggesting that the behaviour of the microbial population varies during the changeable conditions that lead to the formation and removal of the weathering crust.Recall that Q si must be large enough for weathering crust formation to occur and that Q 0 > 0 favours removal (Muller and Keeler, 1969).

Feedbacks between microbial growth, albedo, and melting
The weathering crust and the microbes within it give rise to some potential albedo feedbacks that can lead to amplified reductions in albedo and increases in melting.Here we will discuss the melt-albedo feedback (e.g.Box et al., 2012) and the microbe-albedo feedback (e.g.Tedstone et al., 2020;Mc-Cutcheon et al., 2021) and explore the ability of our model to capture and quantify these feedbacks.We first briefly outline the potential feedbacks.
As the surface of an ice sheet melts, the surface transitions from snow covered to bare ice to having meltwater on the surface.The change in surface type from snow-covered to water-covered surfaces leads to a reduction in the albedo.This in turn means that more radiation is absorbed, so more melting occurs, reducing the albedo further.The result is a positive feedback: the melt-albedo feedback (e.g.Box et al., 2012).
The surface of ice sheets can also be made darker by biological factors, especially glacier ice algae and cryoconite (Hotaling et al., 2021), with glacier ice algae thought to be a key contributor to the darkening of the Greenland ice sheet (Benning et al., 2014).There is thought to be a potential positive feedback between ice algae and meltwater, whereby the presence of algae on the ice-sheet surface lowers the albedo, increasing melting, which releases nutrients that had been frozen in the ice, enabling more algal growth (e.g.Tedstone et al., 2020;McCutcheon et al., 2021).The darkening effect is particularly strong for surface ice algae because the algae have a special dark pigment to protect them from high levels of UV and photosynthetically active radiation (Williamson et al., 2019(Williamson et al., , 2020)).However, a similar, but weaker, feedback can be expected for microbes in the weathering crust.Furthermore, the framework we provide here for studying microbe-albedo feedback in the weathering crust can be easily adapted to ice algae.

Model
To introduce the melt-albedo and microbe-albedo feedbacks into our steadily melting model, we make the absorption coefficient α a linearly increasing function of either the porosity φ or the bulk microbial abundance A bulk = φA w +(1−φ)A i .This captures the fact that more shortwave radiation is absorbed when there is more water (compared to ice) or when there are more microbes.The albedo a is then calculated as an output of the model from the additional boundary condition ( 22), as discussed in Sect. 2. The absorption coefficient functions we chose are to capture the melt-albedo feedback and to capture the microbe-albedo feedback, where we have used that A i = A ∞ in our steadily melting setup.Here, α 0 is the absorption coefficient of pure ice containing the small background abundance of microbes A ∞ .The specific forms of Eqs. ( 51) and ( 52) have been chosen as a simple way to include the feedbacks in our model.The value 3 we have chosen in Eqs. ( 51) and ( 52) is arbitrary and using other positive values leads to qualitatively the same results.We are not aware of any literature detailing the dependence of the absorption coefficient on porosity and microbial abundance.

Results
By making the absorption coefficient a function of the porosity φ and/or the microbial abundance A w , the steadily melting weathering crust problem can no longer be solved analytically.Furthermore, a dependence on A w fully couples the microbe and nutrient problem to the weathering crust https://doi.org/10.5194/tc-17-1967-2023 The Cryosphere, 17, 1967Cryosphere, 17, -1987Cryosphere, 17, , 2023 1.
problem, meaning that we can no longer solve the weathering crust problem on its own first.Therefore, we must solve the full problem numerically.It is still possible to find the solution in Z > Z m analytically, whereas the solution in Z < Z m is found numerically using the boundary value problem solver bvp4c in MATLAB.Our code is made available at https://doi.org/10.5281/zenodo.7199159(Woods, 2022).
As expected, making the absorption of shortwave radiation increase with water content or microbial abundance causes a reduction in albedo in our model (Fig. 11c).We find that including the meltwater coupling (Eq.51) reduces the albedo to 0.47 compared to 0.6 in the uncoupled case with α = α 0 .Similarly, the microbe coupling (Eq.52) reduces the albedo of the steadily melting solution to 0.51 from 0.6.
Figure 10 shows a comparison between the porosity, internal radiation, and microbial abundance distributions in the coupled and uncoupled cases.The effects of including meltwater coupling and microbe coupling are similar.In both cases, more radiation is absorbed by the ice than in the uncoupled case (seen in the reduction in albedo).However, the internal radiation does not increase everywhere -it is higher than the uncoupled case near the surface but lower at depth.This is because there is more water and more mi-crobial biomass near the surface, so absorption of radiation is greatest near the surface, which in turn means that there is less radiation left to penetrate deeper into the ice.
The overall effect of including the coupling is to decrease the size of the weathering crust (decrease Z m ), lower the porosity, and reduce the total amount of microbes (Fig. 10).This might seem counter-intuitive; the feedbacks previously discussed would suggest that making the absorption depend on water content or microbial abundance would lead to an increase in water content and microbial abundance since more water and microbes leads to greater absorption, resulting in even more melting and microbial growth.However, we are observing the previously discussed removal of microbes in runoff proposed by Cook et al. (2012).The increase in absorption has led to an increase in surface melting, lowering, and runoff, with enhanced removal of microbes in runoff dominating any increase in microbial growth.Therefore, including the coupling leads to an upward translation of the microbial abundance distribution, as we saw in Sect.3.3.2.
Figure 11 reproduces the plots of the total number of microbes in the water A tot,w and residence time t res as a function of Q si from Fig. 9; this time it shows the results for meltwater coupling, microbe coupling, and the uncoupled case  1. from Fig. 9 (grey).The variation with Q si for Q 0 < 0 (top row) and Q 0 > 0 (middle row) is shown, as well as the variation with Q 0 (bottom row).These show that including the meltwater or microbe coupling reduces the total number of microbes and reduces their residence time, with the main effect when Q 0 > 0 being to reduce the second local maximum.This is because the region of significant algal growth is smaller with coupling (Fig. 10), so the increase in t res from the local minimum to the second local maximum, which corresponds to the spread of the significant microbial growth region deeper into the ice (discussed around Fig. 9), is reduced.
Note that Fig. 11 shows that the total number of microbes, residence time and albedo are all lower with meltwater cou-pling than microbe coupling, but this should not be read into too much.The forms of the absorption coefficient functions ( 51) and ( 52) are illustrative only and could easily be modified such that the microbe coupling has a larger effect (e.g. by replacing 3 in Eq. ( 52) with a larger value).The main point to take away from this investigation is that including coupling between meltwater, microbes, and radiation absorption can lead to a reduction in albedo but also a reduction in the size of the weathering crust and the total microbial abundance within the weathering crust aquifer.Due to excluding surface ice algae from our model, we are unable to draw any conclusions about the total microbial abundance of the whole https://doi.org/10.5194/tc-17-1967-2023 The Cryosphere, 17, 1967Cryosphere, 17, -1987Cryosphere, 17, , 2023 system, comprised of both microbes in the weathering crust aquifer and algae on the ice surface.

Discussion and conclusion
We have presented a mathematical model for the weathering crust and the microbes within it, as well as a nutrient for the microbes.We have focussed on one-dimensional steadily melting solutions, where the internal and surface melt rates balance to produce a weathering crust of a steady size.We have assumed that all surface meltwater runs off, removing any microbes or nutrients that were in the water.In reality, the lateral flow of meltwater across the ice-sheet surface and through the weathering crust is important for the transport of nutrients and microbes across the ice (e.g.Stibal et al., 2012), but our one-dimensional model is unable to capture this.
We have made a number of simplifications in order to produce a reduced problem from which we can gain informative insights.One simplification of reality is our treatment of the microbes and their nutrients.We have assumed that there is a single type of phototrophic microbe and a single representative nutrient for the microbe.In reality, the weathering crust contains a range of different microbes: primary producers (e.g.algae and filamentous cyanobacteria), which use solar radiation as their energy source, and heterotrophs, which use autochthonous organic carbon produced by the primary producers as their energy source (Stibal et al., 2012).Our model could be thought of as describing a single representative type of primary producer, which requires both solar radiation and nutrients for growth.Alternatively, it could be interpreted as modelling (crudely) the combined influence of all microbes.Clearly there are more complex interactions that could be added at the expense of introducing further parameters and potentially losing interpretability.
Similarly, the reality is that microbes require multiple nutrients to survive and grow, namely carbon (C), nitrogen (N), phosphorous (P), and a variety of micronutrients.The representative nutrient in our model can be interpreted as whichever of the three macronutrients (C, N, or P) is limiting.In the case of glacier ice algae and cryoconite microbial communities, P is usually the limiting nutrient (McCutcheon et al., 2021;Stibal and Tranter, 2007;Bagshaw et al., 2013).
Our model has given insight into the porosity and microbe distributions in the weathering crust.Our model shows that the porosity decreases non-linearly with depth, in agreement with previous observations (Schuster, 2001;Cooper et al., 2018) and modelling (Schuster, 2001).However, the surface values of porosity in our model are larger than observed by Schuster (2001) and Cooper et al. (2018) (∼ 0.85 compared to ∼ 0.65).This could be due to our assumption that the entire crust is saturated, neglecting the observed unsaturated surface layer.The presence of near-surface water in our model leads to increased absorption and melting compared to if the near-surface was unsaturated.Furthermore, we would expect the crust to disintegrate once the porosity gets above a certain limit, but this is not included in our model.We have also found the microbial abundance in the weathering crust generally decreases with depth.This agrees with the expectation of Cook et al. (2016) and the observations of Christner et al. (2018).Since radiation decreases with depth, it is expected that there will be more primary producers (with protective pigments) near the surface and more heterotrophs (which are less adapted to light) deeper down (Cook et al., 2016).The decrease in radiation with depth also leads to a decrease in porosity and melting with depth.
Furthermore, we have observed that the climate forcing (represented by the shortwave radiation Q si and combined longwave radiation and turbulent heat fluxes Q 0 in our model) has a clear impact on the size of the weathering crust and the residence time and abundance of microbes in the weathering crust.Our model shows that it is the balance between internally absorbed radiation and surface forcings that controls these quantities.The relationship is complex, but, in general, increasing the surface forcings relative to internal radiation (for example, increasing the turbulent heat fluxes and longwave radiation Q 0 ) causes the size of the weathering crust, the residence time, and the total number of microbes in the crust to all decrease.This is due to enhanced surface melting and runoff causing removal of the weathering crust and the microbes within it.This is in line with the observation that warm, windy, overcast conditions (high turbulent heat fluxes) favour weathering crust removal (Muller and Keeler, 1969), as well as the suggestion that excessive warming could lead to a removal of microbes in runoff, threatening the Greenland ice sheet's state of autotrophy (net removal of CO 2 ) (Cook et al., 2012).
Our model also allows for the investigation of positive feedbacks.We have explored melt-albedo and microbealbedo feedbacks by making the absorption of shortwave radiation a function of the water and microbe content of the weathering crust, respectively.These feedbacks lead to enhanced melting of the ice-sheet surface.Our model suggests that both feedbacks act to reduce the weathering crust size and number of microbes in the weathering crust aquifer, as well as their residence times.This is due to enhanced melting leading to an increased runoff rate, so more microbes get removed from the surface.Hence, these feedbacks are selflimiting -the reduction in albedo cannot keep going unlimited because the increased melting will lead to the removal of the microbes and meltwater that cause the albedo reduction.
However, we caution that we have focussed on microbes that live within the weathering crust aquifer (transported by meltwater), and have ignored the glacier ice algae that live on the surface of the weathering crust (affixed to the ice).An area for future work is to extend our model to include a population of surface algae that can persist in the presence of meltwater runoff, necessitating a removal of our assumption that all microbes get washed away in runoff.This would enable more to be said about the total system micro-bial abundance of the weathering crust, including that stored on the surface.Moreover, the surface algae provide a major contribution to the albedo via the algal blooms that occur on the south-western Greenland ice sheet.Extending our modified surface ice algae model to higher dimensions would enable an investigation of the spatio-temporal dynamics of algal blooms and their impact on albedo.
In a similar vein, our model could be used to investigate the temporal dynamics of the weathering crust across a full melt season at locations in the field, such as in the "dark zone" on the south-western Greenland ice sheet.This would allow us to better understand the changeable behaviour of the weathering crust in response to variable weather conditions.Another benefit would be to provide a physical approach to developing a bare-ice albedo scheme.The bare-ice albedo schemes used in major climate models are often very simple and do not compare well with in situ observations (e.g.Fettweis et al., 2017), so alternative approaches are desirable.
In conclusion, we have shown that a mathematical model can be used to capture and explore key characteristics and behaviours of the weathering crust and the microbes within it.In particular, even with a number of simplifications, it shows that the response of the weathering crust and its microbial community to changes in atmospheric forcing is complex, and can lead to both increases and decreases in microbial abundance.We acknowledge that the model has crudely approximated many aspects of the dynamics but hope that this study can provide a framework into which more detailed biogeochemistry and physics can be incorporated in the future.

Figure 1 .
Figure 1.(a) Photo of a weathering crust (credit: Eva Doting).(b) Schematic of a weathering crust -a near-surface layer of porous ice (above the dashed line), below which the ice is impermeable (though the "solid" ice shown here may contain small, closed-off air bubbles).The porosity of the weathering crust decreases with depth.Much of it is saturated with meltwater, with a thin unsaturated layer at the surface.The saturated layer is known as the "weathering crust aquifer".Schematic based on Cooper et al. (2018); Moure et al. (2022).

Figure 2 .
Figure 2. Schematic of the steadily melting setup.Panel (a) shows how the surface of the ice z = h(t) lowers over time in the steadily melting state, with the ice melting in such a way that it looks like the ice column is simply being translated downwards.Panel (b) shows this column of ice in the Z = h(t) − z frame.The forcings Q si , Q 0 , and T ∞ and the surface melt rate m surf are also shown.

Figure 3 .
Figure 3. Example solution for T ∞ = −10 • C (dark blue) and T ∞ = −1 • C (grey).The dotted lines show the melting depth Z m , with Z m ≈ 1.8 m when T ∞ = −10 • C and Z m ≈ 3.3 m when T ∞ = −1 • C. The rate of change in surface height is ḣ ≈ −1.6 cm d −1 when T ∞ = −10 • C and ḣ ≈ −1.7 cm d −1 when T ∞ = −1 • C. The surface melt rate is m surf ≈ 0.25 cm d −1 in both cases.The other parameter values are those shown in Table 1.The shaded region on the left panel corresponds to the smaller region shown in the other panels.

Figure 4 .
Figure 4. (a, b) The range of Q si and Q 0 values for which steadily melting states with both internal and surface melting are possible, as well as the regions for which there is no surface melting and no internal melting (grey).Contours of (a) rate of surface lowering − ḣ and (b) weathering crust thickness Z m are shown, ranging from − ḣ = 0 cm d −1 to − ḣ = 20 cm d −1 and from Z m = 0 m to Z m = 1.9 m with 1 cm d −1 and 0.1 m spacing between contours, respectively.(c) The dependence of Z m on the shortwave radiationQ si when Q 0 = 20 W m −2 , Q 0 = −20 W m −2 , and Q 0 = 0 W m −2 .The other parameter values are those shown in Table1.

Figure 5 .
Figure 5. Schematic illustration of how microbes (big purple dots) and nutrients (small pink dots) are included in the model.Both are transferred from ice (grey) to water (blue) as the ice melts (red arrows), and the microbes in the water then consume nutrients in the water (bright blue arrow).

Figure 6 .
Figure 6.Volume-averaged microbial abundance φA w and nutrient concentration φC w solutions scaled with the maximum concentrations A max and C ∞ .Solutions correspond to the weathering crust solutions (dark blue lines) in Fig. 3.The solutions are shown for β A = 20 d −1 (dark blue line) and β A = 5 d −1 (grey line).The dotted line shows Z = Z m .The other parameter values are those shown in Table1.

Figure 7 .
Figure 7. Volume-averaged microbial abundance φA w and nutrient concentration φC w solutions scaled with the maximum concentrations A max and C ∞ .Panel (a) shows the solutions for β C /C ∞ = 10 −3 mL cell −1 d −1 (light blue) and β C /C ∞ = 10 −5 mL cell −1 d −1 (yellow).Panel (b) shows the solutions for k C = 2C ∞ (dark blue) and k C = 0.2C ∞ (green).These solutions correspond to the weathering crust solutions (dark blue lines) in Fig. 3.The dotted lines show the melting depth Z m ≈ 1.78 m.The other parameter values are those shown in Table1.

Figure 8 .
Figure 8.The effect of Q si and Q 0 on volume-averaged microbial abundance φA w and nutrient concentration φC w solutions scaled with the maximum concentrations A max and C ∞ .Panel (a) shows the solutions for Q 0 = −20 W m −2 (negative) in the case of Q si = 150 W m −2 (light blue) and Q si = 300 W m −2 (dark blue).Panel (b) shows the solutions for Q 0 = 20 W m −2 (positive) in the case of Q si = 50 W m −2 (green) and Q si = 300 W m −2 (dark blue).The dotted lines show the weathering crust thickness Z m .The other parameter values are those shown in Table1.

Figure 9 .
Figure 9. Dependence of (a, b, c, f) total amount of microbes in the water A tot,w and (d, e) approximate residence time of the microbes t res on the incoming shortwave radiation Q si and other surface radiation Q 0 .Panels (a), (c), and (d) show the solutions for Q 0 = 20 W m −2 (dark blue line), Q 0 = −20 W m −2 (light blue line), and Q 0 = 0 W m −2 (grey dashed line), respectively.Panels (b) and (e) show solutions with Q si = 200 W m −2 .Panel (c) is an extension of panel (a) to larger Q si values.Panel (f)shows contours of A tot,w in the region for which steadily melting states with both internal and surface melting are possible.The contour spacing is 0.02 × 10 10 cells m −2 .The regions for which there is no surface melting and no internal melting are also shown (grey, see Fig.4).The other parameter values are those shown in Table1.

Figure 10 .
Figure 10.Comparison of steadily melting solutions with (dark and light blue) and without (grey) albedo feedbacks.Panel (a) shows meltwater coupling, where α = α 0 (1 + 3φ), so that the absorption increases with water content.The albedo in the coupled case is 0.47, compared to 0.6 in the uncoupled case.Panel (b) shows microbe coupling, where α = α 0 (1+3φ(A w −A ∞ )), so that the absorption increases with the bulk microbial abundance.The albedo in the coupled case is 0.51.The dotted lines show the melting depth Z m .The parameter values are those shown in Table1.

Figure 11 .
Figure 11.Dependence of (a, d, g) total amount of microbes in the water A tot,w , (b, e, h) approximate residence time of the microbes t res , and (c, f, i) albedo a on the incoming shortwave radiation Q si (a-f) and other surface radiation Q 0 (e-g), in the uncoupled (grey), meltwater coupling (dark blue), and microbe coupling (light blue) cases.The value of Q 0 is −20 W m −2 (negative) in (a)-(c) and 20 W m −2 (positive) in (d)-(f).The value of Q si is 200 W m −2 in (g)-(i).The other parameter values are those shown in Table1.

Table 1 .
Default values of the parameters in the model.Unless otherwise stated, these are the parameter values used.