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 coefficient 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 fluxes, finding that the relative amounts of shortwave (surface-penetrating) radiation and other heat fluxes 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.
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, and Hotaling 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 CO2, 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).
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., 2017, 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.
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 mint (with negative mint 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.
2.1.1 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 ui, uw, and ua are their velocities. Since air is much less dense than ice and water (), 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 ui 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 ui and uw to zero for the specific one-dimensional solution we consider.
2.1.2 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 () then the region is at the melting temperature Tm. If there is no water present (S=0), then the ice can be below the melting temperature T≤Tm.
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, ci and cw are the specific heat capacities and ki and kw are the thermal conductivities for ice and water, respectively, and ℒ 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<Tm), there cannot be any internal melting (mint=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=Tm), melting is possible (mint≥0) and water can be present (S≥0). Since we know that T=Tm, the energy Eq. (4) instead determines the internal melt rate mint.
2.1.3 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 and Feltham, 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 , where is the upward-pointing unit vector, with F± satisfying
where α and r are absorption and scattering coefficients, respectively, with . 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−.
2.1.4 Kinematic conditions
The kinematic conditions for the ice and the water at the surface z=h(t) are, respectively,
where and are the upward ice and water velocities (which, in the case of water, we expect to be negative), is the velocity of the surface, msurf is the rate of surface melting, and roff 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 roff is therefore only non-zero when the surface is saturated.
2.1.5 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., 2008, 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 (Qsi 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 (ϵ 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 Ta is the reference air temperature 2 m above the ice surface), and latent heat flux available for surface melting at rate msurf.
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. (2008, 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)Qsi 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 z=h(t) is . The upward flux F+ is aQsi. 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 is an effective heat transfer coefficient, combining turbulence and radiation.
2.1.6 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 ui and uw 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 Qsi and the “other” surface radiation Q0 (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, msurf and mint, and surface runoff, roff. The default parameter values used in the model are shown in Table 1.
a Motivated by Hoffman et al. (2014). b Taylor and Feltham (2005). c Calculated from κ, a, and χ using Eq. (24) and definitions of κ and s. d Jørgensen and Bendoricchio (2001). e Motivated by Christner et al. (2018). f Inferred from van den Broeke et al. (2011). g Motivated by Holland et al. (2019). Many of the ice properties are taken from Cuffey and Paterson (2010).
2.2 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 , thermal conductivity , and specific heat capacity . 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 (). 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 mint 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 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 become
Both the ice and water mass conservation Eqs. (1) and (2) reduce to
and the energy Eq. (4) becomes
The advection terms 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=Tm at the surface in order for surface melting to be possible. In addition to these, we also prescribe a cold temperature T∞<Tm at depth, so
Finally, the boundary conditions for the two-stream Eqs. (12) and (13) become
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 Qsi, the combined longwave radiation and turbulent heat fluxes Q0, 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 mint, surface melt rate msurf, runoff roff, and rate of change of surface height . We choose values of the forcings such that the surface is melting (msurf>0) and the deep ice is cold. Therefore, as shown in Fig. 2, we expect there to be a temperate region near the surface where the ice is porous (ϕ>0) and at the melting temperature (T=Tm). This region is the weathering crust. In Z>Zm, the ice is solid (ϕ=0) and cold (T<Tm). We refer to Zm 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=Zm, combined with our assumption that the ice is saturated with water, justifies our assumption of no water flow (uw=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.
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 decreases exponentially with depth,
where the extinction coefficient κ is related to the absorption and scattering coefficients α and r by and (e.g. Taylor and Feltham, 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 mint. As mentioned previously, we seek solutions with a temperate region above a cold region Z>Zm. The boundary conditions between the two regions are
where the final condition comes from continuity of conductive heat flux. Also, since T=Tm everywhere in , 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>Zm. Evaluating this at Z=Zm and using the boundary conditions (26) gives us an expression which can be rearranged to tell us Zm 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 mint from the internal radiation F as follows:
In turn, ϕ can be calculated from mint using the mass conservation Eq. (14), giving
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, . 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=Zm.
For our assumption that the surface is melting to be valid, we need the surface melt rate msurf to be positive. Using the expression (27) for msurf, this gives a restriction on the possible values of the prescribed radiations Qsi and Q0,
Furthermore, the solution found here, with a porous ice region occupying , is only valid if the calculated Zm 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 Zm) increases. This is represented by the term
in the expression (34) for Zm. 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 Zm) decreases, since the porous ice at the surface gets removed as meltwater. Surface lowering is represented by the term
in the expression (34) for Zm. 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 Tm, 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 msurf, 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 Zm>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 Zm>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 Qsi and combined longwave radiation and turbulent heat fluxes Q0 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 Qsi, increasing the other surface radiation Q0 increases the amount of surface lowering but not the internal heating. This is because Q0 is completely absorbed at the surface. If Q0 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 Q0 is too small, there is not enough energy to melt the surface.
The weathering crust thicknesses Zm 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 Zm depends on the incoming shortwave radiation Qsi and the combined longwave radiation and turbulent heat fluxes Q0. There is a qualitative difference between the behaviour of Zm with Qsi, depending on the sign of Q0. We see that Zm increases with Qsi when Q0>0 but decreases with Qsi when Q0<0. The reason for this relates back to our interpretation of Zm being determined by a balance between internal heating and surface lowering. If Q0<0, doubling Qsi doubles the internal heating term (36) and more than doubles the surface lowering term (37). Therefore, as Qsi increases, surface lowering becomes more significant compared to internal heating, so the porous region (and hence Zm) decreases. On the other hand, if Q0>0, surface lowering becomes less significant as Qsi increases, so the porous region (and hence Zm) increases. In both cases, a weathering crust can only exist if the shortwave radiation Qsi is large enough. Physically, positive Q0 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 Q0 becomes too large. Negative Q0 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 Q0 as slightly negative (van den Broeke et al., 2008, 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., 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 Q0 to explore the different behaviours exhibited during these changeable conditions. The sensitivity of the weathering crust thickness Zm to the values of Qsi and Q0, demonstrated by our model might partly explain why the weathering crust tends to form and be removed by changing weather conditions.
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 Aw and ice Ai (cells mL−1) separately, and we do the same for the nutrient concentrations Cw and Ci (µmol L−1).
To model the water microbial abundance Aw, we adapt a simple logistic growth model (Murray, 1993),
On the left-hand side, we have the time derivative and an advection term, with the microbes being advected at the water velocity uw. On the right-hand side we have a growth term, a death term (with a death rate dw), and source and sink terms due to melting and freezing, respectively. The growth term, of the form , is the standard form for logistic growth, with Amax 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 ϕCw 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. kPAR and kC 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 di) and are fixed within the ice (so they are advected at the ice velocity ui). Hence Ai satisfies
We expect that di≥dw, 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).
Equations (39), (41), (42), and (43) for the microbial abundances Aw and Ai and nutrient concentrations Cw and Ci are coupled to Eqs. (1), (2), (4), (5), and (6) for the weathering crust evolution since they depend on both porosity ϕ and radiation F.
3.2 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 . In this particular situation, Eqs. (41) and (43) for Ai and Ci 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 Ai=A∞ and Ci=C∞ everywhere.
Referring back to Fig. 2, water is only present in the porous region , so it is in this region that we solve Eqs. (44) and (45) for Aw and Cw. For continuity at the interface Z=Zm 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 ϕ, mint, , and Zm. Then, using these results, Eqs. (44) and (45) can be solved (numerically) with the boundary conditions (46) and (47) for Aw and Cw. 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 cells mL−1 and A∞=102 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 Cw 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∞.
3.3.1 Effect of microbial growth rate
Figure 6 shows an example solution of how the volume-averaged microbial abundance ϕAw and nutrient concentration ϕCw 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=Zm), there is an inner region where the microbial abundance Aw and nutrient concentration Cw are still approximately their deep-ice values, and hence (=0.01ϕ) and . Above this, Aw increases towards the surface and Cw 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 Cw changes significantly from the deep-ice value C∞, causing the bulk nutrient concentration ϕCw to move off the exponential porosity profile.
3.3.2 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 kC. 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 Aw and nutrient concentration Cw 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 Cw is lower when βC is larger. Correspondingly, the microbial abundance Aw is lower because the nutrient has been used up more quickly by the microbes, so less nutrient is available for further microbial growth.
The second nutrient parameter kC is a measure of the point at which the nutrient stops being a limiting factor for the microbes. The larger kC is, the higher the volume-averaged nutrient concentration ϕCw must be before the microbial growth rate becomes independent of ϕCw. Equivalently, the smaller kC 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 Aw and nutrient concentration Cw with depth for two values of kC – one 10 times larger (dark blue) than the other (green). The effect of decreasing kC 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 ϕCw decreases with depth in the region near Z=Zm where Cw≈C∞. Therefore, as kC becomes smaller, the depth at which the nutrient becomes a limiting factor (controlled by the size of ϕCw relative to kC) 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.
3.3.3 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 Q0 and Qsi in Fig. 8. Secondly, we look at the dependence of the total amount of microbes in the water Atot,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 Qsi and Q0 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 Qsi and Q0, including positive and negative Q0, to capture the range of conditions experienced by microbes in the weathering crust. For a weathering crust to form, the shortwave radiation Qsi must be large enough. Positive Q0 corresponds to conditions of high turbulent heat flux (warm, cloudy, windy), which promote removal of the weathering crust. Negative Q0 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 Q0>0 and for small values of Qsi – below the critical value required to produce internal melting – there are no water microbes because there is no water. Then, as Qsi increases, Atot,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 Zm (i.e. the region in which microbes can grow). As shown in Fig. 9c, if Qsi is increased further (to less physical values), the total amount of microbes Atot,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 Atot,w.
However, when Q0<0, the total number of microbes Atot,w decreases monotonically with increasing Qsi (compared to increasing when Q0>0). The reason for this relates back to Fig. 4. As previously discussed, when Q0<0, the melting depth Zm decreases as Qsi increases because surface 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 Atot,w depends on the other surface radiation Q0, which is a combination of net longwave radiation and turbulent heat fluxes, for a fixed value of shortwave radiation Qsi. We see that Atot,w is a decreasing function of Q0. This is because increasing Q0 increases the surface melting without affecting the internal melting or microbial growth. Hence, the surface runoff increases without microbial growth increasing to compensate. When Q0 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 Qsi and Q0 on Atot,w is shown in the contour plot in Fig. 9f. In addition to the behaviour captured by Fig. 9a–c, this shows how increasing Q0 to more positive values causes the maximum of the blue Atot,w against Qsi curve in Fig. 9c to decrease and move to a higher value of Qsi. As discussed, the total number of microbes Atot,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 Qsi further no longer increases the microbial growth significantly. When Q0 is more positive, this point is delayed to larger Qsi because the size of the region in which microbial growth is possible, Zm, takes longer to tend to its maximum value (see Fig. 4). Hence, for more positive Q0, Qsi must get larger before Atot,w starts decreasing with Qsi.
Importantly, the key result from Fig. 9f is that the total number of microbes in the water is largest when Q0 is negative and Qsi is just large enough to allow surface melting, which are the conditions under which weathering crust growth is most favoured.
3.3.4 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 sheet – in particular whether or not the ice sheet is a net source or sink of CO2. Cook et al. (2012) showed the Greenland ice sheet to be in a stable state of autotrophy (net removal of CO2 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 CO2 from the atmosphere.
In our steadily melting model, we approximate the average residence time by
is the microbial runoff rate at the surface. Equation (49) for tres 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 tres 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 tres on the incoming shortwave radiation Qsi. This is shown for positive (dark blue line) and negative (light blue line) values of Q0. As with the total number of microbes Atot,w, the observed behaviour is quite different depending on the sign of Q0.
When Q0>0, the global behaviour (ignoring the local minimum in Fig. 9d which we will discuss shortly) is that as Qsi increases, tres initially increases then decreases again, for the same reasons as those just discussed around Atot,w: as Qsi 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 Qsi 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, . As Qsi increases, the microbial abundance Aw initially does not increase very much from A∞ (seen from the flat section in Fig. 9a for small Qsi). 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 near-surface growth dominates, the microbial runoff rate q (which depends on the surface value of the microbial abundance) grows much more rapidly with Qsi than the total number of microbes Atot,w because the microbial abundance is only growing significantly near the surface. Hence, the residence time tres decreases. Once significant microbial growth begins spreading into the interior, Atot,w grows more quickly than q, so tres increases again.
However, when Q0<0, tres decreases monotonically with increasing Qsi (as does Atot,w). The residence time decreases for the same reason that it decreases for large Qsi when Q0>0 – increasing Qsi 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 Qsi is missed because Qsi cannot take such small values when Q0<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 tres is a decreasing function of Q0 for a fixed value of Qsi. This is for the same reasons (already discussed) that Atot,w is a decreasing function of Q0 for a fixed value of Qsi.
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 Qsi and Q0 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 Qsi and the sign of the combined longwave radiation and turbulent heat fluxes Q0, 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 Qsi must be large enough for weathering crust formation to occur and that Q0>0 favours removal (Muller and Keeler, 1969).
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; McCutcheon 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, 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.
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 . 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 Ai=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.
By making the absorption coefficient a function of the porosity ϕ and/or the microbial abundance Aw, the steadily melting weathering crust problem can no longer be solved analytically. Furthermore, a dependence on Aw fully couples the microbe and nutrient problem to the weathering crust 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>Zm analytically, whereas the solution in Z<Zm 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 microbial 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 Zm), 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 Atot,w and residence time tres as a function of Qsi from Fig. 9; this time it shows the results for meltwater coupling, microbe coupling, and the uncoupled case from Fig. 9 (grey). The variation with Qsi for Q0<0 (top row) and Q0>0 (middle row) is shown, as well as the variation with Q0 (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 Q0>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 tres 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 coupling 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 system, comprised of both microbes in the weathering crust aquifer and algae on the ice surface.
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 Qsi and combined longwave radiation and turbulent heat fluxes Q0 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 Q0) 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 CO2) (Cook et al., 2012).
Our model also allows for the investigation of positive feedbacks. We have explored melt–albedo and microbe–albedo 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 self-limiting – 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 microbial 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.
The MATLAB code used to produce the results and figures in this article is made available at https://doi.org/10.5281/zenodo.7199159 (Woods, 2022).
No data sets were used in this article.
TW and IJH developed the model. TW wrote the code and produced the results. TW wrote the paper with input from IJH.
The contact author has declared that neither of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We thank Liz Bagshaw and Chris Williamson for useful discussions about microbes, as well as reviewers Andrew Tedstone and Sammie Buzzard for their very helpful comments. We also thank Eva Doting for the weathering crust photo in Fig. 1 and Ian Stevens for providing us with more photos that did not make the final manuscript.
Tilly Woods was supported by the UK Engineering and Physical Sciences Research Council through a doctoral studentship. Ian J. Hewitt was supported by the UK Natural Environment Research Council (grant no. NE/S011390/1).
This paper was edited by Elizabeth Bagshaw and reviewed by Sammie Buzzard and Andrew Tedstone.
Bagshaw, E. A., Tranter, M., Fountain, A. G., Welch, K., Basagic, H. J., and Lyons, W. B.: Do cryoconite holes have the potential to be significant sources of C, N, and P to downstream depauperate ecosystems of Taylor Valley, Antarctica?, Arct. Antarct. Alp. Res., 45, 440–454, https://doi.org/10.1657/1938-4246-45.4.440, 2013. a
Benning, L. G., Anesio, A. M., Lutz, S., and Tranter, M.: Biological impact on Greenland's albedo, Nat. Geosci., 7, 691, https://doi.org/10.1038/ngeo2260, 2014. a, b
Box, J. E., Fettweis, X., Stroeve, J. C., Tedesco, M., Hall, D. K., and Steffen, K.: Greenland ice sheet albedo feedback: thermodynamics and atmospheric drivers, The Cryosphere, 6, 821–839, https://doi.org/10.5194/tc-6-821-2012, 2012. a, b
Chevrollier, L.-A., Cook, J. M., Halbach, L., Jakobsen, H., Benning, L. G., Anesio, A. M., and Tranter, M.: Light absorption and albedo reduction by pigmented microalgae on snow and ice, J. Glaciol., 69, 333–341, https://doi.org/10.1017/jog.2022.64, 2022. a
Christner, B. C., Lavender, H. F., Davis, C. L., Oliver, E. E., Neuhaus, S. U., Myers, K. F., Hagedorn, B., Tulaczyk, S. M., Doran, P. T., and Stone, W. C.: Microbial processes in the weathering crust aquifer of a temperate glacier, The Cryosphere, 12, 3653–3669, https://doi.org/10.5194/tc-12-3653-2018, 2018. a, b, c, d, e
Cook, J., Flanner, M., Williamson, C., and McKenzie Skiles, S.: Bio-optical Properties of Terrestrial Snow and Ice, in: Springer Series in Light Scattering, edited by: Kokhanovsky, A., Springer, Cham, https://doi.org/10.1007/978-3-030-20587-4_3, 2019. a
Cook, J. M., Hodson, A. J., Anesio, A. M., Hanna, E., Yallop, M., Stibal, M., Telling, J., and Huybrechts, P.: An improved estimate of microbially mediated carbon fluxes from the Greenland ice sheet, J. Glaciol., 58, 1098–1108, https://doi.org/10.3189/2012JoG12J001, 2012. a, b, c, d
Cook, J. M., Hodson, A. J., and Irvine-Fynn, T. D.: Supraglacial weathering crust dynamics inferred from cryoconite hole hydrology, Hydrol. Process., 30, 433–446, https://doi.org/10.1002/hyp.10602, 2016. a, b, c, d
Cook, J. M., Hodson, A. J., Gardner, A. S., Flanner, M., Tedstone, A. J., Williamson, C., Irvine-Fynn, T. D. L., Nilsson, J., Bryant, R., and Tranter, M.: Quantifying bioalbedo: a new physically based model and discussion of empirical methods for characterising biological influence on ice and snow albedo, The Cryosphere, 11, 2611–2632, https://doi.org/10.5194/tc-11-2611-2017, 2017. a, b
Cook, J. M., Tedstone, A. J., Williamson, C., McCutcheon, J., Hodson, A. J., Dayal, A., Skiles, M., Hofer, S., Bryant, R., McAree, O., McGonigle, A., Ryan, J., Anesio, A. M., Irvine-Fynn, T. D. L., Hubbard, A., Hanna, E., Flanner, M., Mayanna, S., Benning, L. G., van As, D., Yallop, M., McQuaid, J. B., Gribbin, T., and Tranter, M.: Glacier algae accelerate melt rates on the south-western Greenland Ice Sheet, The Cryosphere, 14, 309–330, https://doi.org/10.5194/tc-14-309-2020, 2020. a, b
Cooper, M. G., Smith, L. C., Rennermalm, A. K., Miège, C., Pitcher, L. H., Ryan, J. C., Yang, K., and Cooley, S. W.: Meltwater storage in low-density near-surface bare ice in the Greenland ice sheet ablation zone, The Cryosphere, 12, 955–970, https://doi.org/10.5194/tc-12-955-2018, 2018. a, b, c, d, e, f, g
Cooper, M. G., Smith, L. C., Rennermalm, A. K., Tedesco, M., Muthyala, R., Leidman, S. Z., Moustafa, S. E., and Fayne, J. V.: Spectral attenuation coefficients from measurements of light transmission in bare ice on the Greenland Ice Sheet, The Cryosphere, 15, 1931–1953, https://doi.org/10.5194/tc-15-1931-2021, 2021. a, b
Cuffey, K. M. and Paterson, W. S. B. (Eds.): The Physics of Glaciers, Elsevier, 4th edn., ISBN 9780123694614, 2010. a
Fausto, R. S., van As, D., Mankoff, K. D., Vandecrux, B., Citterio, M., Ahlstrøm, A. P., Andersen, S. B., Colgan, W., Karlsson, N. B., Kjeldsen, K. K., Korsgaard, N. J., Larsen, S. H., Nielsen, S., Pedersen, A. Ø., Shields, C. L., Solgaard, A. M., and Box, J. E.: Programme for Monitoring of the Greenland Ice Sheet (PROMICE) automatic weather station data, Earth Syst. Sci. Data, 13, 3819–3845, https://doi.org/10.5194/essd-13-3819-2021, 2021. a
Fettweis, X., Box, J. E., Agosta, C., Amory, C., Kittel, C., Lang, C., van As, D., Machguth, H., and Gallée, H.: Reconstructions of the 1900–2015 Greenland ice sheet surface mass balance using the regional climate MAR model, The Cryosphere, 11, 1015–1033, https://doi.org/10.5194/tc-11-1015-2017, 2017. a
Hoffman, M. J., Fountain, A. G., and Liston, G. E.: Surface energy balance and melt thresholds over 11 years at Taylor Glacier, Antarctica, J. Geophys. Res.-Earth Surf., 113, F04014, https://doi.org/10.1029/2008JF001029, 2008. a
Hoffman, M. J., Fountain, A. G., and Liston, G. E.: Near-surface internal melting: A substantial mass loss on Antarctic Dry Valley glaciers, J. Glaciol., 60, 361–374, https://doi.org/10.3189/2014JoG13J095, 2014. a, b
Holland, A. T., Williamson, C. J., Sgouridis, F., Tedstone, A. J., McCutcheon, J., Cook, J. M., Poniecka, E., Yallop, M. L., Tranter, M., Anesio, A. M., and The Black & Bloom Group: Dissolved organic nutrients dominate melting surface ice of the Dark Zone (Greenland Ice Sheet), Biogeosciences, 16, 3283–3296, https://doi.org/10.5194/bg-16-3283-2019, 2019. a, b
Hotaling, S., Lutz, S., Dial, R. J., Anesio, A. M., Benning, L. G., Fountain, A. G., Kelley, J. L., McCutcheon, J., Skiles, S. M. K., Takeuchi, N., and Hamilton, T. L.: Biological albedo reduction on ice sheets, glaciers, and snowfields, Earth-Sci. Rev., 220, 103728, https://doi.org/10.1016/j.earscirev.2021.103728, 2021. a, b, c
Irvine-Fynn, T. D., Edwards, A., Stevens, I. T., Mitchell, A. C., Bunting, P., Box, J. E., Cameron, K. A., Cook, J. M., Naegeli, K., Rassner, S. M., Ryan, J. C., Stibal, M., Williamson, C. J., and Hubbard, A.: Storage and export of microbial biomass across the western Greenland Ice Sheet, Nat. Commun., 12, 3960, https://doi.org/10.1038/s41467-021-24040-9, 2021. a, b, c
Jørgensen, S. E. and Bendoricchio, G.: Fundamentals of Ecological Modelling, Elsevier, Amsterdam, New York, 3rd edn., vol. 21, ISBN 9780080440156, 2001. a, b, c
Kelliher, F. M., Owens, I. F., Sturman, A. P., Byers, J. N., Hunt, J. E., and Mcseveny, T. M.: Radiation and ablation on the névé of Franz Josef Glacier, J. Hydrol., 35, 131–150, 1996. a
Law, R., Arnold, N., Benedek, C., Tedesco, M., Banwell, A., and Willis, I.: Over-winter persistence of supraglacial lakes on the Greenland Ice Sheet: Results and insights from a new model, J. Glaciol., 66, 362–372, https://doi.org/10.1017/jog.2020.7, 2020. a
Liston, G. E., Winther, J.-G., Bruland, O., Elvehøy, H., and Sand, K.: Below-surface ice melt on the coastal Antarctic ice sheet, J. Glaciol., 45, 273–285, https://doi.org/10.3189/S0022143000001775, 1999. a
Maykut, G. A. and Untersteiner, N.: Some results from a time- dependent thermodynamic model of sea ice, J. Geophys. Res., 76, 1550–1575, https://doi.org/10.1029/jc076i006p01550, 1971. a
McCutcheon, J., Lutz, S., Williamson, C., Cook, J. M., Tedstone, A. J., Vanderstraeten, A., Wilson, S. A., Stockdale, A., Bonneville, S., Anesio, A. M., Yallop, M. L., McQuaid, J. B., Tranter, M., and Benning, L. G.: Mineral phosphorus drives glacier algal blooms on the Greenland Ice Sheet, Nat. Commun., 12, 570, https://doi.org/10.1038/s41467-020-20627-w, 2021. a, b, c, d
Moure, A., Jones, N., Pawlak, J., Meyer, C., and Fu, X.: A thermodynamic nonequilibrium model for preferential infiltration and refreezing of melt in snow, ESS Open Archive, https://doi.org/10.1002/essoar.10512834.1, 2022. a
Muller, F. and Keeler, C. M.: Errors in short-term ablation measurements on melting ice surfaces, J. Glaciol., 8, 91–105, https://doi.org/10.3189/S0022143000020785, 1969. a, b, c, d, e, f
Murray, J. D. (Ed.): Mathematical Biology, Springer-Verlag, Berlin, Germany, Heidelberg, Germany, second, corrected edn., ISBN 9783662085424, 1993. a
Perovich, D. K.: Theoretical estimates of light reflection and transmission by spatially complex and temporally varying sea ice covers, J. Geophys. Res., 95, 9557–9567, https://doi.org/10.1029/jc095ic06p09557, 1990. a
Schuster, C. J.: Weathering crust processes on melting glacier ice (Alberta, Canada), PhD thesis, https://scholars.wlu.ca/etd/489/ (last access: 16 December 2022), 2001. a, b, c, d, e, f, g, h, i
Shakeel, M., Matthews, P. C., Graham, R. S., and Waters, S. L.: A continuum model of cell proliferation and nutrient transport in a perfusion bioreactor, Math. Med. Biol., 30, 21–44, https://doi.org/10.1093/imammb/dqr022, 2013. a
Stibal, M. and Tranter, M.: Laboratory investigation of inorganic carbon uptake by cryoconite debris from Werenskioldbreen, Svalbard, J. Geophys. Res.-Biogeo., 112, G04S33, https://doi.org/10.1029/2007JG000429, 2007. a
Stibal, M., Šabacká, M., and Žárský, J.: Biological processes on glacier and ice sheet surfaces, Nat. Geosci., 5, 771–774, https://doi.org/10.1038/ngeo1611, 2012. a, b, c, d
Stibal, M., Box, J. E., Cameron, K. A., Langen, P. L., Yallop, M. L., Mottram, R. H., Khan, A. L., Molotch, N. P., Chrismas, N. A., Calì Quaglia, F., Remias, D., Smeets, C. J., van den Broeke, M. R., Ryan, J. C., Hubbard, A., Tranter, M., van As, D., and Ahlstrøm, A. P.: Algae Drive Enhanced Darkening of Bare Ice on the Greenland Ice Sheet, Geophys. Res. Lett., 44, 463–11, https://doi.org/10.1002/2017GL075958, 2017. a
Taylor, P. D. and Feltham, D. L.: A model of melt pond evolution on sea ice, J. Geophys. Res.-Oceans, 109, 1–19, https://doi.org/10.1029/2004JC002361, 2004. a, b
Taylor, P. D. and Feltham, D. L.: Multiple stationary solutions of an irradiated slab, J. Cryst. Growth, 276, 688–697, https://doi.org/10.1016/j.jcrysgro.2004.11.417, 2005. a, b, c, d
Tedesco, M., Doherty, S., Fettweis, X., Alexander, P., Jeyaratnam, J., and Stroeve, J.: The darkening of the Greenland ice sheet: trends, drivers, and projections (1981–2100), The Cryosphere, 10, 477–496, https://doi.org/10.5194/tc-10-477-2016, 2016. a
Tedstone, A. J., Cook, J. M., Williamson, C. J., Hofer, S., McCutcheon, J., Irvine-Fynn, T., Gribbin, T., and Tranter, M.: Algal growth and weathering crust state drive variability in western Greenland Ice Sheet ice albedo, The Cryosphere, 14, 521–538, https://doi.org/10.5194/tc-14-521-2020, 2020. a, b, c, d, e, f
Telling, J., Stibal, M., Anesio, A. M., Tranter, M., Nias, I., Cook, J., Bellas, C., Lis, G., Wadham, J. L., Sole, A., Nienow, P., and Hodson, A.: Microbial nitrogen cycling on the Greenland Ice Sheet, Biogeosciences, 9, 2431–2442, https://doi.org/10.5194/bg-9-2431-2012, 2012. a
van den Broeke, M., Smeets, P., Ettema, J., van der Veen, C., van de Wal, R., and Oerlemans, J.: Partitioning of melt energy and meltwater fluxes in the ablation zone of the west Greenland ice sheet, The Cryosphere, 2, 179–189, https://doi.org/10.5194/tc-2-179-2008, 2008. a, b, c, d
van den Broeke, M. R., Smeets, C. J. P. P., and van de Wal, R. S. W.: The seasonal cycle and interannual variability of surface energy balance and melt in the ablation zone of the west Greenland ice sheet, The Cryosphere, 5, 377–390, https://doi.org/10.5194/tc-5-377-2011, 2011. a, b, c, d, e
Williamson, C. J., Anesio, A. M., Cook, J., Tedstone, A., Poniecka, E., Holland, A., Fagan, D., Tranter, M., and Yallop, M. L.: Ice algal bloom development on the surface of the Greenland Ice Sheet, FEMS Microbiology Ecology, 94, 3, https://doi.org/10.1093/femsec/fiy025, 2018. a
Williamson, C. J., Cameron, K. A., Cook, J. M., Zarsky, J. D., Stibal, M., and Edwards, A.: Glacier algae: A dark past and a darker future, Front. Microbiol., 10, 519, https://doi.org/10.3389/fmicb.2019.00524, 2019. a
Williamson, C. J., Cook, J., Tedstone, A., Yallop, M., Mccutcheon, J., Poniecka, E., Campbell, D., Irvine-Fynn, T., Mcquaid, J., Tranter, M., Perkins, R., and Anesio, A.: Algal photophysiology drives darkening and melt of the Greenland Ice Sheet, P. Natl. Acad. Sci. USA, 117, 5694–5705, https://doi.org/10.1073/pnas.1918412117, 2020. a
Woods, T.: Weathering crust and microbial activity code, Zenodo [code], https://doi.org/10.5281/zenodo.7199159, 2022. a, b, c