A continuum model for meltwater flow through compacting snow

Meltwater is produced on the surface of glaciers and ice sheets when the seasonal energy forcing warms the snow to its melting temperature. This meltwater percolates into the snow and subsequently runs off laterally in streams, is stored as liquid water, or refreezes, thus warming the subsurface through the release of latent heat. We present a continuum model for the percolation process that includes heat conduction, meltwater percolation and refreezing, as well as mechanical compaction. The model is forced by surface mass and energy balances, and the percolation process is described using Darcy’s law, allowing for both partially and fully saturated pore space. Water is allowed to run off from the surface if the snow is fully saturated. The model outputs include the temperature, density, and water-content profiles and the surface runoff and water storage. We compare the propagation of freezing fronts that occur in the model to observations from the Greenland Ice Sheet. We show that the model applies to both accumulation and ablation areas and allows for a transition between the two as the surface energy forcing varies. The largest average firn temperatures occur at intermediate values of the surface forcing when perennial water storage is predicted.


Introduction
Meltwater percolation into surface snow and firn plays an important role in determining the impact of climate forcing on glacier and ice-sheet mass balance.Percolated meltwater may refreeze, run off, or be stored as liquid water.Since meltwater that runs off from the surface ultimately contributes to sea-level rise, and can influence ice dynamics if it is routed to the ocean via the ice-sheet bed, understanding the proportion of meltwater that runs off is important in assessing the health of glaciers and ice sheets under atmospheric warming (Harper et al., 2012;Enderlin et al., 2014;Forster et al., 2014;Koenig et al., 2014;Machguth et al., 2016).The balance between runoff, refreezing, and storage is controlled by the mechanics and thermodynamics of the porous snow.These processes also underlie the rate of compaction of firn into ice and therefore control the average temperature and accumulation rate that provide surface boundary conditions to numerical ice-sheet models (which typically do not include the compacting firn layer explicitly).
Liquid water that is produced at the surface holds a substantial quantity of latent heat.If the meltwater percolates into the snow and refreezes, it releases the latent heat to warm the snow.Humphrey et al. (2012) observe that the snow at 10 m depth in Greenland is often more than 10 • C warmer than the mean annual air temperature because of the refreezing of meltwater.If, however, this water runs off through supraglacial streams or drains to the bed through moulins, the latent heat is carried away and subsequent cooling of the surface in the winter means that the remaining snow is relatively cold.Since the capacity to store and/or refreeze meltwater is tied to the porosity of the snow, which is in turn linked to the amount of storage and refreezing that have occurred in previous years, it is of interest to know how the partitioning of meltwater between runoff, refreezing, and storage, as well as the firn temperature and density profiles, depends on climatic forcing (air temperature and radiative forcing as well as accumulation).This question is of interest even under steady climate conditions (i.e., seasonally periodic, without any yearon-year trend), and this forms the focus of our study.A further question of current interest is how the firn responds transiently to year-on-year increases in melting (Harper et al., 2012;Koenig et al., 2014), but we consider the steady problem a prerequisite to understanding such transient response.
Our approach in this paper is to construct a continuum model for meltwater percolating through porous snow, along similar lines to Gray (1996).This contrasts cell-based numerical models that are often applied to the Greenland Ice Sheet, such as the Firn Densification Model (IMAU-FDM) that is incorporated in the regional climate model RACMO and is described by Ligtenberg et al. (2011) and Kuipers Munneke et al. (2014, 2015).That model includes mechanical compaction and a "tipping-bucket" hydrology scheme, where the firn is divided into distinct layers and water fills each layer up to the irreducible water content and then trickles instantaneously into the lower layers.Runoff occurs when the water reaches an impermeable layer and the water is removed (representing the lateral flow that occurs in reality).Steger et al. (2017a, b) used a similar tipping-bucket method in the SNOWPACK model (Bartelt and Lehning, 2002) and compared the results to the IMAU-FDM.Wever et al. (2014) compared the tipping-bucket and Richards equation formulations within SNOWPACK to field observations and found that using Richards' equation provided a better fit.The Richards equation approach, in which water flow is driven by gravity and capillary pressure, is similar to the model we adopt in this study.This has also been used in a number of more theoretical models for the percolation of meltwater through snow (Colbeck, 1972(Colbeck, , 1974(Colbeck, , 1976)).Gray andMorland (1994, 1995) and Gray (1996) provide detailed descriptions of this approach in the context of mixture theory.
We now summarize an outline of the paper.In Sect.2, we construct our continuum model for the firn layer and describe its conversion to an enthalpy formulation that facilitates the numerical solution method.In Sect.3, we analyze two test problems that involve the propagation of a refreezing front moving into cold snow and a saturation front filling pore space.These act as benchmarks for the numerics and elucidate some of the generic dynamics that occur within the model.In Sect.4, we impose a more realistic surface energy forcing, corresponding to a periodic seasonal cycle, to examine the effect of climate variables on the fate of the meltwater and the resulting thermal structure of the snow.

Percolation through porous ice
Here we describe our model for the flow of meltwater through porous, compacting snow.We keep track of the flow of water, mechanical compaction, and the melt-refreezing of water into the snow.A volume fraction 1−φ is solid ice while the void space φ is composed of water and air.We define the saturation S as the fraction of the void space that is filled by water (see the schematic in Fig. 1).Conservation of mass for ice, water, and air is expressed as The water partially saturates the snow near the surface (S < 1) whereas, at depth, all of the air is replaced by water and the snow is fully saturated (S = 1).Panel (b) shows an ablation area where the there is fully saturated porous snow in a thin layer near the surface and the underlying ice is solid, advecting into the domain from upstream.Ice grains make contact in the third dimension (into the page) and similarly many of the air and water pockets are connected in the third dimension.

∂(Sφρ
where the subscripts i, w, and a indicate ice, water, and air, respectively.The densities ρ i , ρ w , and ρ a are constants.The velocities of the ice, water, and air are given by u i , u w , and u a .The variable density of the snow is The rate at which ice melts and turns into meltwater internally is given by m and is therefore a source in Eq. ( 2) and a sink in Eq. ( 1).This term is always negative, i.e., refreezing, and in fact is zero except on refreezing interfaces.We assume that the air density is negligible and henceforth neglect Eq. (3).
The flow of water is governed by Darcy's law, i.e., where p w is the water pressure, k(φ) is the permeability, k r (S) is the relatively permeability, and µ is the viscosity of the water.For the permeability we use a simplified Carman- where d p is a typical grain size (Gray, 1996).Table 1 provides the parameter values we use later.
We must now distinguish between partially saturated (S < 1) and fully saturated (S = 1) flow.When the snow is partially saturated, capillary forces drive flow along liquid bridges connecting ice crystals (Bear, 1972).Thus, we relate the water pressure to the capillary pressure p c by p w = p a − p c , where p a is the air pressure (taken as zero).Both the capillary pressure and the relative permeability are prescribed functions of the saturation S. We take where γ is surface tension, and we choose the exponents α and β such that β = α + 1, which avoids a singularity in k r (S)p c (S) at S = 0 (Gray, 1996).
If the snow is fully saturated, water pressure p w is constrained by mass conservation.Combining Eqs. ( 1), (2), and (4) gives which is an elliptic problem for p w .Boundary conditions for this equation are provided by the constraint that p w must be continuous across the interfaces between partially and fully saturated regions and the constraint of no flow across impermeable boundaries (e.g., ice lenses).

Compaction
One of the difficult aspects of modeling firn in a percolation zone is that both mechanical compaction and refreezing combine to control the changes in snow density.There are various empirical parametrizations of dry compaction that can be used; these typically relate the rate of change of density, or equivalently porosity, to quantities such as depth, accumulation rate, temperature, and grain size.In our context, these can be expressed using the material derivative where C is a parametrization of the rate of compaction (Arthern et al., 2010).The appropriateness of such models for snow containing meltwater is uncertain.The parametrizations represent the rearrangement and growth of snow crystals and the accompanying closure of air voids as functions of temperature and accumulation rate, and these processes may be modified by the presence of liquid between crystals.In the absence of a more developed theory for wet compaction, we take the approach of using these dry parametrizations but modify the material density derivative to include the rate of refreezing that is calculated from the thermodynamics.Therefore, we have where the specific compaction rate we choose is the Herron and Langway (1980) model, which is written as C = cφ.The coefficient c in units of yr −1 is given by where a is the accumulation rate in meters of water equivalent per year and T is the absolute temperature.This is an empirical parametrization, and the two forms reflect a change in dominant compaction processes at a certain snow density.
Other parametrizations for compaction that could easily be incorporated in this framework are discussed by Zwally and Li (2002), Reeh (2008), and Morris and Wingham (2014).
We have chosen to use the Herron and Langway model here for simplicity; from the experiments we have conducted, different formulations do not appear to qualitatively change our results.
Combining with Eq. ( 1), we note that Eq. ( 9) is equivalent to where w i is the vertical component of the ice velocity.

Temperature
We assume that ice and water are at the same temperature and therefore any region containing meltwater (S > 0) is at the melting point T m .In regions without water, we solve the temperature evolution equation, where the heat capacity is c p and the thermal conductivity is K = (1 − φ)K.The latent heat term −L m operates on interfaces of refreezing, where it is singular and causes discontinuities to occur in the temperate gradient.

Surface boundary conditions
Here we write boundary conditions on the surface z s (t), which we assume is locally flat, and we write w i and w w as the vertical velocities.The kinematic conditions are where żs is the velocity of the surface, M is the rate of melting, a is the accumulation rate, R is the rainfall rate, and r is www.the-cryosphere.net/11/2799/2017/The Cryosphere, 11, 2799-2813, 2017 runoff, all expressed in units of water equivalent per year.The compaction Eq. ( 9) requires a boundary condition, φ = φ 0 , when the accumulation rate is greater than the rate of melting (i.e., w i − żs < 0), where (1 − φ 0 )ρ i is the bulk density of freshly deposited snow.The energy balance on the surface provides a boundary condition for the temperature equation when the surface temperature is less than T m , and determines the rate of melting M when T = T m .These conditions are combined as along with the conditions M = 0 when T < T m , and M ≥ 0 when T = T m .The forcing energy flux Q(t) includes the combined effects of radiative, turbulent, and sensible heat fluxes.We assume that this is prescribed in order to provide a simple parametrization of the climate forcing.However, it can be related to more specific components of the energy balance as described in Appendix A. The heat transfer coefficient h represents a combination of radiative and turbulent heat transfer.We expect Q to have a typical magnitude on the order of Q 0 = 200 W m −2 with a comparable seasonal amplitude and take h = 14.8 W m −2 K −1 as a representative constant (Cuffey and Paterson, 2010;van den Broeke et al., 2011).

Numerical method
Our complete model is given by ice and water conservation (1) and ( 2), Darcy's law (4), compaction (9), and temperature evolution (12), subject to the boundary conditions ( 13)-( 15).The model is forced by a prescribed energy flux Q, accumulation a, and precipitation R, and it predicts the temperature, porosity, and saturation profiles as well as the surface melt rate, runoff, refreezing, and storage of liquid water.
In this section, we rewrite the equations in a form that we use for our numerical solutions.There are two steps: first, we combine the equations as conservation equations for total water (ice and liquid water) and enthalpy (sensible and latent heat).Using this approach, commonly referred to as the enthalpy method, we can avoid tracking the phase change interfaces and can solve for their location using inequalities (Hutter, 1982;Aschwanden et al., 2012;Hewitt and Schoof, 2017).The second step is to change variables into a frame that moves with the ice surface.At this stage we also simplify the model to write it in one vertical dimension, and we make the Boussinesq approximation to ignore density differences so that ρ i = ρ w = ρ.
We define the total water as W, which is the sum of liquid and solid fractions, i.e., and we define the enthalpy as the sum of sensible and latent heat as The inverse relationships that relate the enthalpy H and total water W to the temperature T , saturation S, and porosity φ are and We define the depth below the ice surface Z, and the relative downward ice velocity w i , as We combine the conservation Eqs. ( 1) and ( 2) with Darcy's law (4) and temperature evolution (12) as where the downward water flux is Combining ice conservation (1) with compaction (9) gives and water pressure is given by The surface boundary conditions (at Z = 0), re-expressed in terms of W and H, are given as In these conditions, the runoff r is assumed to be zero unless the snow at the surface reaches full saturation, in which case Eqs. ( 25) and (26) determine r.At the bottom of our domain, we assume that the conductive heat flux and pressure gradients vanish to replicate effective matching conditions to the deep interior of the ice sheet.On internal interfaces between fully saturated and partially saturated regions we apply p w = −γ /d p to ensure pressure continuity.
We discretize the conserved fluxes in space using a finite volume method implemented in MATLAB (see the Supplement for code).In this construction, the value of each variable is constant in each cell center and the velocities and fluxes are evaluated at cell edges, thereby transferring fluxes of each variable from one cell to another.We evolve Eqs. ( 20)-( 21) in time using explicit forward Euler time stepping, which involves evaluating the fluxes on the cell edges using the quantities from the previous timestep.For advection, we use an upwind scheme where the value of the variable advected depends on the velocity direction.For edges between partially saturated cells, we evaluate the water fluxes using the capillary pressure for p w on the adjacent cells.For edges between fully saturated cells, we solve Eq. ( 20) with S = 1 as an elliptic equation for p w on the saturated cells, which we then use to evaluate the water fluxes.In order to allow cells to switch from fully to partially saturated, we compute the fluxes using both of these methods on edges between fully and partially saturated cells and choose that which gives the largest flux away from the saturated region.

Test problems
In this section we consider two test problems that demonstrate the model behavior and validate the numerical method.The two problems that we consider here are designed to explore the boundaries between frozen and unfrozen snow (refreezing interfaces) as well as the boundaries between partially and fully saturated snow (saturated interfaces).Both problems ignore mechanical compaction.We start by describing the propagation of rainwater into dry snow.This is similar to the problem studied by Colbeck (1972), Gray (1996), and Durey (2014) and has an approximate analytical solution that provides a useful test case for the enthalpy method.We also compare the results of the analytical solution for the propagation of the meltwater front to temperature data from Humphrey et al. (2012).Secondly, to test the propagation of saturated fronts, we consider an isothermal problem in which the porosity profile is prescribed to decrease with depth.We again investigate how rainwater propagates into the snow, with saturation increasing as the front propa- In panel (a), the snow is initially cold with T = T ∞ and dry, with uniform porosity φ 0 .The rainwater percolates through the snow, refreezes at the interface Z f (t), and releases latent heat that warms the snow.The refreezing decreases the porosity in the upper region so that φ + < φ 0 .In panel (b), the snow is temperate, T = T m , with a porosity profile that decays exponentially with depth.After the snow fully saturates two saturation fronts emerge with Z l propagating downward and Z u upward.
gates down.At a certain point the snow fully saturates and a saturated front propagates up toward the snow surface.

Rainfall into cold snow
We consider the infiltration of rain into cold, dry snow as a test problem.We start with a patch of dry snow (S = 0) with constant porosity (φ = φ 0 ) and temperature (T = T ∞ < T m ).We assume no accumulation and ignore compaction so that the ice is stationary.Furthermore, the porosity is large enough that the snow never fully saturates.At time t = 0 a fixed flux of rain R with a temperature T = T m is applied at the surface Z = 0 and a wetting front at Z = Z f moves down at velocity Żf (we show a schematic in Fig. 2 and the numerical solutions in Fig. 3).Since the capillary pressure gradients are small and the flow is largely driven by gravity, the wetting front behaves as a smoothed shock front.Some of the water at the shock front refreezes, warming the snow ahead.
As shown in more detail in Appendix C, the behavior of this shock can be understood by ignoring the diffusive capillary pressure term.This approximation relegates Eqs. ( 1) and ( 2) to hyperbolic partial differential equations for the porosity and saturation as well as simplifying the temperature Eq. ( 12) www.the-cryosphere.net/11/2799/2017/The Cryosphere, 11, 2799-2813, 2017 so that where k r (S) = dk r /dS, and with initial and boundary conditions Equations ( 28)-( 30) have corresponding jump conditions across the shock which incorporate the refreezing rate −m i at that front.These are Żf where + refers to the region above the front (Z < Z f ).Using these jump conditions and the solutions to Eqs. ( 28)-( 30) subject to the boundary conditions ( 31)-( 34), we find an approximate expression for the front velocity: Note that if T ∞ = T m , i.e., isothermal snow, the front simply propagates at the speed of the draining rainwater R/(φ + S + ) (Bear, 1972).The effect of refreezing due to T ∞ < T m is to slow the front and to cause a decrease in the porosity as the front passes by an amount . This is a mechanism by which ice lenses can form: if the pre-existing porosity is small enough, the porosity above the front can decrease to zero and the pores freeze shut.In this case the front stops propagating, and a saturated region forms above the lens in a similar way to that described in Sect.3.3.We also determine approximate analytical solutions for temperature and saturation, which are compared to the numerical solutions in Fig. 3.The agreement between the numerical and approximate solutions is very good.The approximate temperature profile ahead of the refreezing front is given by

Data comparison
The refreezing and release of latent heat as a front of meltwater moves through a firn layer allows the percolation of meltwater to be observed in englacial temperature data.Harper et al. (2012) and Humphrey et al. (2012) collected temperature data in the accumulation zone on the western flank of the Greenland Ice Sheet and inferred the movement of meltwater by warming of the snow due to the release of latent heat.They set up a vertical string of thermistors to determine the temperature profile in the upper 10 m of the ice sheet.Data from one vertical string between the dates of 5 and 25 July 2007 (days 185-203) are shown in Fig. 4. From these data it is clear that the ice at depth progressively warmed, likely due to the refreezing of liquid meltwater.Over the 12 days between day 185 and day 197, the warming front propagated about a meter, while over the course of the next 6 days from day 197 to day 203 the meltwater penetrated two additional meters, showing a 4-fold increase in front velocity.Humphrey et al. (2012) infer that the warming spike on day 199 is due to an influx of meltwater from lateral sources.A minimum temperature is observed at around 5 m depth and the temperature recorded on the lower thermistors is warmer, which could be due to prior warming by meltwater pulses or a manifestation of the seasonal thermal wave.We now compare these data to the approximate solution for the temperature field ahead of a refreezing front, as given in Eq. ( 39).We fit front speed Żf for the days 185-197 and a larger front speed for days 197-203.The increase in the front speed is likely due to an increase in surface melt.We set the melting temperature T m = 0 • C, fit a constant farfield temperature T ∞ , and use the heat diffusivity for ice K/(ρc p ) = 1.1 × 10 −6 m 2 s −1 (Table 1).In light of the simplified analysis, the fit between Eq. (39) and the Humphrey et al. (2012) data is quite good.

Isothermal saturation fronts
We now consider the propagation of rainwater into isothermal, temperate snow of decreasing porosity such that fully saturated fronts develop.The porosity decreases exponentially with depth as where Z 0 is a constant.We continue to ignore compaction and accumulation, and since the snow is isothermal the porosity is therefore constant in time.Initially, the rain partially saturates the snow and a wetting front moves downward, as shown in Fig. 5a.Then, at a certain depth, the maximum saturation reaches unity and two saturation fronts emerge, one that propagates up and the other down, as shown in Fig. 5b and c.
In Appendix D, we derive the locations of the upper Z u and lower Z l fronts by neglecting flow driven by gradients in capillary pressure.This analysis results in two differential The Cryosphere, 11, 2799-2813, 2017 www.the-cryosphere.net/11/2799/2017/Equations ( 28)-( 30) have corresponding jump conditions across the shock which incorporate the refreezing rate m I at that front.These are Żf [ ] (1 where + refers to the region above the front (Z < Z f ).Using these jump conditions and the solutions to ( 28)-( 30) subject to the boundary conditions ( 31)-( 34), we find an approximate expression for the front velocity, Note that if T 1 = T m , i.e. isothermal snow, the front simply propagates at the speed of the draining rain water R/( + S + ) (Bear, 1972).The effect of refreezing due to T 1 < T m is to slow the front and to cause a decrease in the porosity as the front 10 passes, by an amount 0 + = (1 0 )c p (T m T 1 )/L.This is a mechanism by which ice lenses can form: if the pre-existing porosity is small enough, the porosity above the front can decrease to zero and the pores freeze shut.In this case the front stops propagating, and a saturated region forms above the lens in a similar way to that described in section 3.3.We also determine approximate analytical solutions for temperature and saturation, which are compared to the numerical solutions in figure 3.
The agreement between the numerical and approximate solutions is very good.The approximate temperature profile ahead of 15 the refreezing front is given by  2012) show the propagation of refreezing fronts in Greenland firn.We overlay the approximate temperature solution for the temperature ahead of a refreezing front (black lines, Eq. 39).The speed of the front varies over the 18-day record: dashed lines use the initial speed and the dotted line uses the final speed.The far-field temperature is assumed to be constant in the model whereas the data show a local minimum in temperature at around 5 m, which could be due to prior freezing fronts or the seasonal wave.
equations for the evolution of upper front Z u and lower front Z l : subject to the initial conditions where Z 1 and t 1 are the location and time at which full saturation initiates.We solve these coupled, nonlinear ordinary differential equations (ODEs) using a numerical integrator in MATLAB and compare these semi-analytical solutions to the full numerical solutions in Fig. 5 (the dashed black lines).The slight differences are due to neglecting the gradient in capillary pressure in our approximate solutions.

Results
We now examine the solutions to the full model with prescribed seasonal energy forcing, which we parametrize as a sinusoid, using the annual mean as a control parameter.In principle, we could also incorporate diurnal periodicity, but we choose to ignore it because we expect diurnal variability to affect only a small surface layer (∼ 1 m depth) and we are interested in the full firn column (∼ tens of meters of depth).
For cold ice, the variation of surface energy flux leads to a seasonal temperature wave (Cuffey and Paterson, 2010) and a dry-compaction density profile.This solution breaks down if the surface temperature reaches the melting point during summer, at which point the surface snow melts and the meltwater can percolate through the snow and refreeze, thereby warming the snow through the release of latent heat.Even with a small amount of melting, the resulting temperature profiles become very different from the thermal wave.
We apply an oscillating surface forcing in Eq. ( 26) of the form where Q is the annual mean surface forcing, and we take the amplitude Q 0 = 200 W m −2 and period t 0 = 1 year.For simplicity, we assume a constant accumulation rate and ignore rainfall.
We run a suite of numerical simulations varying the accumulation rate and annual mean surface forcing, each time allowing the dynamics to reach an annual periodic state (typwww.the-cryosphere.net/11/2799/2017/The Cryosphere, 11, 2799-2813, 2017

Isothermal saturation fronts
We now consider the propagation of rain water into isothermal, temperate snow of decreasing porosity such that fully saturated fronts develop.The porosity decreases exponentially with depth as where Z 0 is a constant.We continue to ignore compaction and accumulation and since the snow is isothermal, the porosity is 5 therefore constant in time.Initially, the rain partially saturates the snow and a wetting front moves downward, as shown in figure 5(a).Then, at a certain depth, the maximum saturation reaches unity and two saturation fronts emerge, one that propagates up and the other down, as shown in figures 5(b) and 5(c).
In Appendix D, we derive the locations of the upper Z u and lower Z l fronts by neglecting flow driven by gradients in capillary pressure.This analysis results in two differential equations for the evolution of upper front Z u and lower front Z l , subject to the initial conditions where Z 1 and t 1 are the location and time at which full saturation initiates.We solve these coupled, nonlinear ODEs using a numerical integrator in MATLAB, and compare these semi-analytical solutions to the full numerical solutions in figure 5 (the 15 dashed black lines).The slight differences are due to neglecting the gradient in capillary pressure in our approximate solutions.

12
Nondimensional quantities Nondimensional quantities Nondimensional quantities Figure 5. Evolution of fully saturated fronts at three instances in time, showing (red), water flux (cyan), and water pressure (magenta).The porosity (green) decreases exponentially with depth over length scale Z 0 = /2, where is the characteristic length scale defined in Appendix B and given in Table 1.Panel (a) shows the position of the front before the firn fully saturates.Panels (b, c) show the bidirectional motion of the fully saturated fronts.Dashed black lines show semi-analytical solutions from solving Eq. ( 41).
ically this takes around 10 years).Four representative spacetime diagrams these simulations are shown in Fig. 6.
Each case shows a different value of Q with same accumulation (1.7 m water equivalent per year) and porosity of fresh snow φ 0 = 0.64.While the ice surface moves up and down during the simulation, we plot the quantities as a function of depth below the surface Z = z s (t)−z and plot ice streamlines to show the relative motion of the ice.In Fig. 7 we show how the mean temperature at the bottom of the domain T ∞ and the mean surface temperature T s change as the mean surface forcing varies, for three different values of accumulation rate.Each point in this figure corresponds to an annual average of a periodic simulation such as those in Fig. 6 (and which are labeled in Fig. 7b).
The four simulations in Fig. 6 represent the spectrum of possible surface types on glaciers and ice sheets, encompassing both accumulation and ablation regions.If we interpret increasing Q as a parametrization of slow climate warming, we might expect a location that is initially an accumulation area to transition through each of these states.Figure 6I is an accumulation area where there is no melting at any point during the year.The ice streamlines show that the ice advects downward as more snow accumulates on the surface.The snow compaction is visible from a convergence of the streamlines with time.The temperature variation with depth in this case is just the thermal wave and the variations in surface temperature are only felt around K/(ρc p ω) ∼ 6 m into the snow.
Increasing Q above −Q 0 leads to melting during summer.Figure 6II shows an accumulation area where the temperature and porosity profiles are significantly affected by the meltwater that drains into the snow during the summer.Here there is water below 10 m throughout the year fed by percolation each summer.This is a perennial aquifer, as found in a number of field observations (Forster et al., 2014;Koenig et al., 2014).
Figure 6III shows a region which is an accumulation area but with more melting than in Fig. 6II.Interestingly, this situation no longer has a perennial aquifer and all of the meltwater that is produced refreezes.Although still a percolation zone it is different in character than the region shown in Fig. 6II.The porosity decreases more rapidly with depth in this case so that despite more water being produced on the surface during the summer, this larger quantity of water is not able to percolate as far into the snow.As a consequence, it is not so well insulated from the cold surface during the winter and all of the water refreezes.This greater quantity of refreezing is in turn responsible for the more rapid decrease in porosity with depth that prevents the liquid water percolating as deep as it does in Fig. 6II (more refreezing means the pore space is filled in more effectively with ice).In contrast, the reason a perennial aquifer is sustained in Fig. 6II is because the water penetrates sufficiently far that it is insulated from the cold surface (Kuipers Munneke et al., 2014).
Above a critical Q there is too much melting for the firn to accommodate and runoff begins (this occurs at a value of Q intermediate between Fig. 6III and IV and is clearest to see in Fig. 7b).The transition from an accumulation area to an ablation area occurs when runoff exceeds the accumulation.Figure 6IV shows an ablation area where the surface meltwater is only able to enter a few meters into the snow and reaches the impermeable barrier of the glacial ice surface.During the course of the summer all of the snow is melted as well as some of the glacial ice.The streamlines show net upward motion in this case indicating that there is net ablation over the course of the year.
In Fig. 7 we also calculate the total quantity of surface melt and the partitioning of the melt between runoff, liquid storage in the ice, and refreezing in the firn.Runoff and melt are calculated from the model output, liquid storage is taken to be the total water flux passing out of the bottom of the domain (the domain represents only the surface firn layer, so this represents water that is stored within the upper part of the  ice sheet), and the amount of refreezing is computed as the residual.As shown in Fig. 7b and c, the maximum storage is 0.56 and 1.5 m w.e.yr −1 for accumulation rates of 1.7 and 3.4 m w.e.−1 , respectively.
For Q < −Q 0 no melting occurs and the domain top and bottom temperatures are identical.However, as soon as the annual mean surface forcing increases above −Q 0 , the domain top and bottom temperatures diverge due to the re-lease of latent heat which warms the snow.Depending on the accumulation rate, the average bottom firn temperature can reach the melting point, corresponding to a perennial firn aquifer.This does not occur for smaller accumulation rates, i.e., Fig. 7a, but does for larger accumulation rates, i.e., Fig. 7b and c.Additionally, all three panels show that when Q increases further the bottom firn temperature decreases again.This corresponds to the second type accumulation area shown in Fig. 6III, in which water only penetrates part of the way into the domain before refreezing.When Q is large enough such that the region has become an ablation zone, the bottom temperature (now the temperature of incoming glacial ice) is almost the same as the surface temperature.The largest bottom temperatures occur at intermediate values of surface forcing, considerably lower than the value required to transition to an ablation region.
The thermal structure and water content of the lower firn are strongly tied to the amount of meltwater produced, which in this model is tied directly to the annual mean surface forcing.In a warming world, one can imagine a particular location transitioning from an accumulation to ablation region.Our results in Fig. 7 show that storage and refreezing can accommodate much of the melt that occurs when the warming is not too large.Once the forcing is sufficient for runoff to start, the amount of refreezing decreases slightly so that an increasingly large fraction of the melt runs off.Most of this runoff is presumably routed to the glacier bed and then the ocean.As well as a form of mass loss, the timing and quantity of meltwater delivery to the bed will determine the style of subglacial drainage system that develops and the subsequent ice dynamics (Zwally et al., 2002;Schoof, 2010;Tedstone et al., 2015).

Conclusions
We have described a continuum model for the evolution of firn hydrology, compaction, and thermodynamics.The model is capable of determining the evolution of the firn including the temperature, porosity, and water content.The model differs from other models of firn hydrology in its treatment of the percolation of water, for which we use Darcy's law and a parametrization of capillary pressure.Our treatment for runoff also differs in that we assume that water runs off when the surface layer of snow is fully saturated rather than assuming runoff at depth when the percolating water first reaches an impermeable ice layer.
The model applies to both accumulation and ablation areas.Given the forcing (energy flux and accumulation rate), the model selects which of these applies to any particular region.One of the useful outputs of the model is an indication of how the firn may change as function of climate warming, as revealed by moving from left to right in Fig. 7.In agreement with Kuipers Munneke et al. ( 2014) and Steger et al. (2017a) we find that perennial firn aquifers occur when there Table 2. Typical numerical values for the surface energy balance (Cuffey and Paterson, 2010;van den Broeke et al., 2011). is sufficiently high accumulation and sufficient melting occurs.
In the future, we hope to extend this work beyond the one-dimensional solutions presented here.In principle the model applies to fully three-dimensional geometries, when the slope of the saturated surface (the "water table" in the firn) will allow meltwater to flow laterally as well as vertically.The data from Humphrey et al. (2012) suggest the occurrence of "piping events", where meltwater forms a vertical channel and breaks through to depths where the snow is much colder.These events could be captured in a twodimensional framework, and it is possible that a theory allowing the solid ice and liquid water to have different temperatures may help explain these features.On a larger scale, the horizontal scales of the ice sheet are much larger than the depth of the firn, so a reduced, vertically integrated version of this theory may also be useful.
The use of Darcy's law requires an estimate for the permeability and the relative permeability.The comparison of our model behavior with the data from Humphrey et al. (2012) in Fig. 4 is encouraging and suggests that these parameters could be determined with detailed measurements of surface melt and snow temperatures.Here we have interpreted the porosity and the permeability as grain-scale properties.An alternative interpretation that might be appropriate on larger scales would treat these as averages over fractures, pipes, and ice lenses to give a macroscopic effective porosity and permeability.
Although we have focused on idealized, periodic simulations, the model can be forced by real climatological data or coupled to a regional atmospheric model.The model could also be coupled to an ice-sheet model, using the deep firn temperature T ∞ as the surface boundary condition for the ice sheet.
Data availability.The data associated with this paper are contained in Humphrey et al. (2012) or can be produced from the code attached in the Supplement.which states that the mass −m i that freezes from the liquid phase enters the solid phase and that the latent heat from refreezing warms the dry ice below.We can simplify these equations since θ = 0 in the upper portion (+), φ = φ 0 and S = 0 in the lower portion (−), and the ice velocity w i is zero, so q − φ + S + Żf = −m i , (C11) After a short initial transient, the solution approximates a traveling wave in which the upper region 0 < Z < Z f has θ = 0, φ = φ + (to be determined shortly), and q = R. Since B 1, this means Uk(φ + )k r (S + ) ≈ R, which determines the constant S + in the upper region (there is a narrow boundary layer behind the front, in which S + changes rapidly but q − φ + S + Żf is constant; see below).
We next solve for the temperature evolution in the lower region.Assuming that the freezing front moves quickly, i.e., | Żf | 1 (this is appropriate since U is large), we can move into a translating frame Z = Z − Z f and neglect the time dependence so that is constant (set by the far-field temperature), and hence θ ≈ θ ∞ 1 − e −Pe Żf Z .This is the approximate solution given dimensionally in Eq. ( 39).From the temperature field we can determine the melt rate using Eq.(C13) as which is negative, corresponding to freezing, since θ ∞ < 0, and Eq.(C12) therefore determines the porosity jump: Finally, the jump condition for water conservation, Eq. (C11), determines the speed of the front as ) This result corroborates the front velocity derived by Colbeck (1972), Gray (1996), andDurey (2014).
To capture the smoothing of the front due to capillary pressure, we can examine the narrow boundary layer behind the front.The relevant scale for this region is of order 1/B, so we write Z − Z f = Ẑ/B and determine the leading-order quasistatic approximation where ψ = φ Żf U k(φ) , which can be integrated to give

−φ Żf
which is similar to the result derived by Gray (1996).

Appendix D: Saturation fronts
Here we calculate the motion of the fully saturated fronts for isothermal conditions with fixed porosity φ = φ 0 e −Z/Z 0 , as in Sect.3.3.We again make use of dimensionless variables.
In the time before full saturation initiates, and in the limit B 1, conservation of water at the wetting front Z f (t) is given, as in the Appendix C with θ ∞ = 0, by where φ f (t) = φ 0 e −Z f /Z 0 is the porosity at the front and S f is the saturation.Using permeability k(φ) = φ 3 and relative permeability k r (S) = S 2 , we can calculate the saturation induced by the rainfall as Thus, the initial evolution equation for the front before full saturation is ) which can be integrated to give The Cryosphere, 11, 2799-2813, 2017 www.the-cryosphere.net/11/2799/2017/

Figure 1 .
Figure1.The three components of meltwater-infiltrated snow: air, water, and ice.Panel (a) shows water infiltrating an accumulation area where the snow density increases with depth and snow advects down.The water partially saturates the snow near the surface (S < 1) whereas, at depth, all of the air is replaced by water and the snow is fully saturated (S = 1).Panel (b) shows an ablation area where the there is fully saturated porous snow in a thin layer near the surface and the underlying ice is solid, advecting into the domain from upstream.Ice grains make contact in the third dimension (into the page) and similarly many of the air and water pockets are connected in the third dimension.

Figure 2 .
Figure2.Schematic of the test problems considered in (a) Sect.3.1 and (b) Sect.3.3.In both panels, rain falls at a rate R on the surface of the snow.White shading indicates dry snow (S = 0), grey indicates partially saturated snow (0 < S < 1), and dark shading indicates fully saturated snow (S = 1).In panel (a), the snow is initially cold with T = T ∞ and dry, with uniform porosity φ 0 .The rainwater percolates through the snow, refreezes at the interface Z f (t), and releases latent heat that warms the snow.The refreezing decreases the porosity in the upper region so that φ + < φ 0 .In panel (b), the snow is temperate, T = T m , with a porosity profile that decays exponentially with depth.After the snow fully saturates two saturation fronts emerge with Z l propagating downward and Z u upward.

Figure 3 .
Figure 3. Evolution of a refreezing front at three instances of time, partitioned between the three components of the enthalpy.The green, red, and yellow colors show the porosity, saturation, and temperature profiles, respectively.The dashed lines show the approximate analytical solutions described in Appendix C. The temperature is made nondimensional by T = Tm + (T1 Tm) T and the parameters are 0 = 0.4, and R = 0.54, along with other values in table 1. 5

Figure 3 .Figure 4 .
Figure 3. Evolution of a refreezing front at three instances of time, partitioned between the three components of the enthalpy.The green, red, and yellow colors show the porosity, saturation, and temperature profiles, respectively.The dashed lines show the approximate analytical solutions described in Appendix C. The temperature is made nondimensional by T = T m + (T ∞ − T m ) T and the parameters φ 0 = 0.4 and R = 0.54, along with other values in Table1.

Figure 5 .
Figure 5. Evolution of fully saturated fronts at three instances in time, showing saturation (red), water flux (cyan), and water pressure (magenta).The porosity (green) decreases exponentially with depth over length scale Z0 = `/2, where `is the characteristic length scale defined in Appendix B and given in table 1. Panel (a) shows the position of the front before the firn fully saturates.Panels (b) and (c) show the bidirectional motion of the fully saturated fronts.Dashed black lines show semi-analytical solutions from solving equation (41).

Figure 6 .Figure 6 .
Figure 6.Space-time diagrams showing the evolution of porosity (top), saturation S (middle), and temperature T (bottom) as a function of time for the accumulation rate a = 1.7 mwe yr 1 and four values of forcing: (I) cold accumulation zone where the mean forcing is Q = Q0 .(II) accumulation area with mean forcing Q = 0.707Q0.In this case, a clear perennial aquifer develops.(III) accumulation area with larger forcing Q = 0.575Q0.(IV) ablation zone with mean forcing Q = 0.146Q0.In all simulations the porosity of the falling snow is 0 = 0.64 and the black lines show ice streamlines.

Figure 7 .
Figure 7. Average meltwater partition (right) and annual mean temperature at the ice surface T s and bottom of the domain T ∞ (left) as a function of the annual mean surface forcing, with accumulation increasing from left to right: (a) a = 0.68 m w.e.yr −1 , (b) a = 1.7 m w.e.yr −1 , and (c) a = 3.4 m w.e.yr −1 .For Q > −Q 0 melting occurs at the surface and meltwater percolation warms the bottom of the domain.Dashed lines in panels (a, b) mark the transition from an accumulation to ablation zone and the roman numerals in panel (b) correspond to the solutions in Fig. 6.

Table 1 .
Table of physical values, derived scales, and nondimensional parameter values (defined in Appendix B).