Introduction

Interactions between magma and ice occur on both Earth and Mars, and take a number of forms (Reference Smellie, Smellie and ChapmanSmellie and Chapman, 2002). Effusive or explosive volcanic eruptions can occur beneath (Reference W.H.Mathews, 1947) or adjacent to (Reference Larsen, Smellie and ChapmanLarsen, 2002) glaciers, with lava flows encroaching onto glacial ice (Reference EinarssonEinarsson, 1949; Reference KjartanssonKjartansson, 1951), and shallow igneous intrusions can occur into ice-rich crustal rocks or beneath glacial ice (Reference Smellie and ChapmanSmellie, 2002; Reference Loughlin, Smellie and ChapmanLoughlin, 2002). It is conventional to refer to the latter type of event as a subglacial eruption, but in an earlier analysis of volcano–ice interactions (Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002) we suggested that it can equally be regarded as the intrusion of a sill at the silicate– ice interface. This has implications for the pressure within the magma and its vesicularity and flow rate, all of which will differ from their values in an equivalent subaerial eruption. On the timescale of mafic dyke propagation, controlled by magma flow speeds typically of the order 0.1–1ms^–1, ice behaves as a brittle solid, and so in addition to subglacial sill formation, dyke penetration into glacial ice or ice-rich sediments is possible, with various consequences. In this paper we focus on re-assessing the bulk thermal effects of (i) lava flows erupted on top of glacial ice; (ii) sill intrusions beneath glacial ice evolving into subglacial lava flows and (iii) dyke intrusions into the interiors of glaciers. Most observations of the consequences of lavas advancing over ice or snow (Reference EinarssonEinarsson, 1949; Reference KjartanssonKjartansson, 1951) relate to basaltic-andesite ‘a’a lavas e.g. in Iceland or on Mt Etna, and most of the data on sub-glacial eruptions comes from Iceland where both mafic and evolved magmas are erupted (e.g. Reference Larsen, Smellie and ChapmanLarsen, 2002). Magma rheology controls the aspect ratios of both flows and intrusions and we show that, although the thickness of a magma body strongly influences the timescale of thermal interactions between magma and ice, it does not greatly influence the overall amount of heat exchange.

Supraglacial Lava Flows

If a lava flow advances over ice (Fig. 1), its base transfers heat to the ice and it also loses heat to the surrounding atmosphere, mainly through the top surface (the edges are neglected here). Let the temperature of the core of the flow be T _c (say 1450 K for a basaltic flow). Very quickly, relative to the duration of the emplacement of the flow and its subsequent cooling, the temperature of the top of the flow T _t either reaches the ambient temperature T _a (if a continuous cooled crust forms) or reaches some average value between T _c and T _a if there is a sheared crust on the flow surface, as occurs in a channelized flow. Provided that water can escape efficiently, the temperature of the base of the flow T _b rapidly becomes equal to, and then stays buffered at, the ice melting point, ∽273.15 K, and remains at this temperature as long as there is ice to be melted. The distribution of water beneath the flow will be determined by the extent and scale of irregularities in the lower flow surface and the pre-existing topography, and therefore the ease with which water can escape at the margins of the flow. The fact that explosions forming pseudo-craters can occasionally occur when lava advances over wet ground (Reference Greeley and FagentsGreeley and Fagents, 2001) implies that water may be trapped in contact with the base of a lava flow in ways that raise its temperature above the boiling point. For the moment, we assume that water can always escape fast enough so that this does not occur, and focus on the water production rate.

We assume that the time required for the emplacement of the flow is small compared with its subsequent cooling time: the links between lava effusion rate, flow thickness and flow length ensure that this is true for most lava flows (Reference Pinkerton and WilsonPinkerton and Wilson, 1994). The key issue is then that the temperatures T _t and T _b can be regarded as fixed, and the problem becomes that of determining the rate of heat flow out of the base of the flow into the ice as the core temperature T _c decreases after emplacement. This scenario is treated by Reference Carslaw and JaegerCarslaw and Jaeger (1959, Article 3.4, equation (1)) if the constant initial internal flow temperature T _c is substituted for their initial temperature distribution function f(x). Then the temperature T at any point a distance x below the surface of a flow of thickness X at any time t after the local arrival of the flow is given by

(1)

where κ is the thermal diffusivity of the lava, taken as 7×10^–7m² s^–1. For large values of t (more than ∽3×10⁴ s, i.e. 8 hours), the series expansion converges after only a few terms, but for shorter times progressively many more (tens of) terms must be included, a process made easy by spreadsheet implementation.

The conducted heat flux per unit area leaving the base of the flow is H = k dT/dx, where k is the thermal conductivity of the lava. For consistency we define k as equal to ₁ c where ρ ₁ is the bulk density of the lava and c the specific heat, averaging ∽900 J kg^–1 K^–1 over the lava cooling interval. Differentiating Equation (1) analytically to obtain the temperature gradient dT/dx involves an extensive amount of algebra, so this was derived numerically by evaluating the temperature at 98%, 99% and 99.9% of the depth of the flow and extrapolating to the gradient at the lower boundary. Given the heat flux H, the rate of decrease in the ice thickness R is given by

(2)

where ρ _i is the ice density, ∽917 kgm^–3, and L _i is the latent heat of fusion of the ice, ∽3.35×10⁵ J kg^–1.

Table 1 shows R as a function of flow thickness and time for representative conditions on Earth, with the flow surface temperature T _t equal to 270 K, consistent with the presence of stable ice before the eruption begins. For mafic compositions, the thinnest of the flows illustrated would be emplaced as pahoehoe lobes, but the other three examples would be ‘a’a flows. These commonly advance over a layer of rubble derived from the thermal and mechanical fragmentation and collapse of the outer skin of the flow front. The thickness of this basal layer λ is typically 0.1–0.3 m, and its presence means that the temperature in the bulk of the flow develops above a pre-existing cooled boundary layer of this thickness. This means that the interaction between the flow and the ice beneath it effectively starts after a time interval Ƭ approximately related to λ by the standard thermal boundary layer relationship λ ² ≈ ≈∽Ƭ. This effective delay is ∽4 hours if λ = 0.1 m and ∽1.5 days if λ = 0:3m. Thus there is never a time when the high heat flux represented by the first column of Table 1 is relevant, and the maximum ice melting rate probably only rarely exceeds 10 mms^–1 (even then only for a short period just after a flow encroaches onto ice).

Fig. 1. Illustration of the three scenarios considered in the text. Left: lava erupted subaerially advances over ice; centre: a dyke feeds a sill intruded at rock–ice interface and the sill develops into a subglacial lava flow; right: a dyke overshoots the rock–ice interface and propagates into the body of the glacier. Diagram based on figures 9 and 10 in Reference Wilson, Head, Smellie and ChapmanWilson and Head (2002), used with permission of the Geological Society of London.

Table 1. Ice melting rate R in mms–1 at base of lava flow as a function of flow thickness in metres and time after start of melting thickness

Table 1 shows that the ice melting rate is always a maximum very early in the melting process. If the melting rate is integrated with respect to the time, using a more detailed version of Table 1, the total amount of water produced can be evaluated as a function of time after the start of melting. As would be expected, it is found that the thicker the flow the longer it takes to melt a given fraction of the total amount of ice that can be melted. Thus, after one day, the absolute thicknesses of ice melted by the 1, 3, 10 and 30m thick flows are 4.95, 4.65, 4.0 and 1.8m respectively, corresponding to 99%, 31%, 8% and 1.2% of the total melting that is ultimately achievable. Conversely, 90% melting is achieved after about 20 hours, one week, 3 months and 3 years respectively. These lengthy delays for thicker flows may explain why various observers have reported a negligible amount of snow melting by advancing lava flows (Reference EinarssonEinarsson, 1949: 5m thick; Reference KjartanssonKjartansson, 1951: 10m thick). Nevertheless, the final total amounts of melting that are predicted by this simple analysis are impressive: the ratio of ice thickness melted/lava flow thickness is nearly constant and almost exactly equal to 5.

A simple comparison of the total heat content of the lava with the latent heat of fusion of the melted ice shows that this ratio is the inevitable consequence of all of the heat from the lower half of the flow being ultimately used to melt ice. There are two reasons why the ratio may be an over-estimate. The treatment so far assumes that no heat is used to warm the escaping water above its melting temperature, whereas there must be some water heating if the water is to escape a finite distance from under the flow. However, it would be necessary to heat all of the water to about 808C in order to halve the ice thickness melted/lava flow thickness ratio to ∽2.5. Field observations during future eruptions should be directed at measuring the temperature of released water.

A second issue is the value to be ascribed to the thermal conductivity of the rubble layer at the base of the flow, the significance of which in reducing the initial heat flow rate has already been addressed. The model represented by Equation (1) assumes uniform thermal properties for the lava flow, and the similarity of the temperature boundary conditions at the top and bottom means that nearly equal amounts of heat are transferred upward and downward from the central parts of the flow. If the thermal conductivity of the basal lava zone were less than that of the rest of the flow, then a disproportionate amount of heat would be conducted upward through the upper part of the flow rather than downward through the lower part, and thus ultimately lost to the atmosphere instead of being involved in ice melting. This scenario is clearly a possibility if either (a) the water layer on top of the ice does not completely flood the interfacial spaces, so that air, a very poor conductor of heat, is trapped at the interface, or (b) the water fills the available space and has a temperature in excess of ∽48C, the temperature of maximum density, in which case, being heated from above, it will not undergo spontaneous convection. An adjustment of this kind would further reduce the ice thickness melted/ lava flow thickness ratio, but it is difficult to see how the value could be reduced to less than ∽2.

A caveat about the above treatment is that it ignores the latent heat of crystallization L _m, liberated by the lava as it cools to the solidus. There will be a glassy rind contributing no latent heat on all flows, but this will not be a large fraction of the flow thickness even for the thinnest flows. An allowance for the latent heat, at least as regards the total thermal budget of the flow, can be made (e.g. Reference Carslaw and JaegerCarslaw and Jaeger, 1959) by increasing the initial lava temperature T _c by an amount ΔT such that the latent heat L _m (∽2×10⁵ J kg^–1) is approximated by the sensible heat c _mΔT, where c _m is the specific heat (∽1200 J kg^–1 K^–1 at magmatic temperatures), so that the required value of ΔT is ∽167 K. However, making this adjustment also artificially increases the heat transfer rate from the base of the flow. The error in completely neglecting the latent heat can be estimated by comparing the latent heat, ∽2×10⁵ J kg^–1, with the total amount of sensible heat liberated by the lava in cooling from ∽1450 K to the ambient (on Earth) ∽273 K, which is ∽(1450 – 273) K ×900 J kg^–1 K^–1 = 1.06 ×10⁶ J kg^–1; thus the error is at most ∽19% and will act in the sense of slightly underestimating the volume of water generated.

Subglacial Eruptions

We have argued above and elsewhere (Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002) that, at least in their early stages, subglacial eruptions should be regarded as the equivalent of sill intrusions at the glacier-bedrock interface (Fig. 1). To a first approximation, the model represented by Equation (1) can initially be applied to this geometry also, with the proviso that the vertical coordinate x is now measured upward from the pre-intrusion bedrock–ice interface. We note that Reference Tuffen, Pinkerton, McGarvie and GilbertTuffen and others (2002) examined heat loss from a subglacial intrusion into a growing hemispherical cavity in the ice; we prefer to treat planar geometry here for simpler and more direct comparison with our other two scenarios. Our model provides values for the initial heat transfer rate into the ice and then the growing water layer above the magma. The applicability breaks down progressively because the lower surface of the magma is not maintained at a constant temperature by efficient heat removal as is the case for a subaerial flow surface; instead the heat being conducted into the underlying rocks raises their temperature and soon greatly reduces the downward heat loss rate. This situation can again be modelled with an analytical solution from Carslaw and Jaeger (1959), in this case Article 2.4.iv. The model assumes that a layer of hot material at temperature T _c is emplaced onto an infinite half-space at temperature T _b, and its top surface is rapidly reduced to temperature T _t. The relevant values here are again T _c = 1450 K and T _t = 273.15 K, and the initial value of T _b is also 273.15 K. The model equation is:

(3)

where x is measured downward through the magma layer of thickness d and erf is the error function.

The extent to which this case differs from that of the subaerial lava overlying ice can be appreciated by comparing the heat loss rates at the lower magma boundaries in the two cases. Table 2 gives the heat fluxes, H = k dT/dx, as a function of time for the lower boundaries of a lava flow and a sill, both 3 m thick. At all times later than a few hours, the heat loss rate from the base of the subglacial flow is about half that from the base of a subaerial flow emplaced onto ice. Given that the subaerial flow lost approximately equal amounts of heat upward to the atmosphere and downward to the ice, this comparison suggests that a subglacial flow will lose about twice as much heat upward to the melting glacier as it loses to the rock substrate. In other words, about two thirds of its entire heat content will be used to melt glacial ice. The resulting ratio of melted ice thickness/magma layer thickness should be ∽6–7. This result is significantly different from that in earlier treatments (Reference Höskuldsson and SparksHöskuldsson and Sparks, 1997; Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002) where it was tacitly assumed that essentially all of the heat content of the magma was used to melt ice. In our treatment (Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002) we found a value of ∽14.5 for this ratio, but we note that the silicate specific heat that was used in the calculations (∽1200 J kg^–1 K^–1) did not make sufficient allowance for the temperature dependence of this quantity, and if we had used the value adopted here (∽900 J kg^–1 K^–1), the ratio would have been ∽10.5.

Table 2. Comparison of heat fluxes in Wm^–2 as a function of time through the lower boundaries of a subaerial lava flow emplaced onto ice and a subglacial flow emplaced onto old volcanic rocks underlying a glacier

Dyke Intrusions in to Glaciers

Our inference that ice can be treated as a brittle elastic solid on the timescale of mafic intrusions (Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002) implies that we should consider the heat transfer from a dyke intruded sub-vertically into ice (Fig. 1). As regards heat transfer from the dyke magma, this case is very similar to that of the subaerial flow emplaced on top of ice, because the boundary conditions using the simple model of Equation (1) are very similar: both faces of the magma body are maintained at a temperature close to 273.15 K. It is likely that very good thermal contact will be maintained for a while between the chilled dyke margin and the surrounding meltwater, thus eliminating concerns about reduced thermal conductivity at the boundary.

However, instabilities are likely to develop at the interface as the water layer becomes thicker and cooling cracks develop in the chilled margin; collapse of the disaggregating chilled rock will continually expose fresh hot magma to the water. The heat transfer rate will still be determined mainly by the temperature gradient in the outermost part of the magma, but this temperature gradient will propagate into the magma faster than if the dyke thickness were constant. Without a model of the fracturing mechanism dependent upon the chilling process, it is not possible to estimate by how much the heat transfer rate will be enhanced. But it is clear that, in this scenario, almost all of the heat content of the dyke will ultimately be transferred to the water. The ratio of ice thickness melted/lava flow thickness closely equal to 5 found earlier for a subaerial flow was the consequence of all of the heat from half of the flow being transferred to the ice. Thus in this case the ratio of total ice thickness melted/dyke thickness should approach close to double this value, i.e. 10. This is very close to the adjusted version (∽10.5, see previous section) of the value estimated in our earlier treatment (Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002) of subglacial intrusions, in which we also assumed that all of the heat from the intrusion was transferred to the water.

Discussion and Summary

The new values presented here for the efficiency of heat transfer in various magma–ice contact situations are improvements on earlier estimates. The maximum ratios of thickness of ice that can be melted to the thickness of magmatic heat source are likely to be ∽2–5 for subaerial lava flows encroaching onto glaciers, ∽6–7 for subglacial lava flows and ∽10 for dykes intruded into glacial ice. The ratio for lava advancing over ice is the least reliable as there are numerous uncertainties about the likely temperature and escape rate of the water produced. We stress that all of the models used here are approximations to the actual systems, and infer that there is a need for a more detailed, entirely numerical approach to these problems, which we plan to address in future work. Nevertheless, systematic measurements of volumes and release rates of water can place important constraints on the volumes and intrusion rates of magma involved in subglacial volcanic activity.