Introduction

The initial phases of subglacial eruptions occur within cavities melted into the glacier base by magmatic heat (e.g. Reference Höskuldsson and SparksHöskuldsson and Sparks, 1997; Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002; Reference Guðmundsson, White, Smellie and ClagueGuðmundsson, 2003). The eruptive behaviour involves interaction between rising magma and meltwater, which may lead to the formation of pillow lavas, quench-hyaloclastite or phreatomagmatic tephra, depending upon the mechanism of magma–water interaction (e.g. Reference Moore, Hickson and CalkMoore and others, 1995; Reference Smellie and HoleSmellie and Hole, 1997; Reference Tuffen, Gilbert and McGarvieTuffen and others, 2001).

Studies of subaqueous volcanism have shown that the mechanism of magma–water interaction depends upon the confining pressure (e.g. Reference Moore and SchillingMoore and Schilling, 1973), magma–water ratio (Reference WohletzWohletz, 1983), magma volatile content (Reference McBirneyMcBirney, 1963; Reference Dixon, Stolper and HollowayDixon and others, 1995; Reference Wright, Gamble and ShaneWright and others, 2003) and the ability of magma and water to mix (e.g. Reference Zimanowski, Buttner and LorenzZimanowski and others, 1997). Similarly, the mechanisms of subglacial eruptions are controlled by the cavity conditions such as confining pressure and abundance of meltwater (Reference Tuffen, Gilbert and McGarvieTuffen, 2001; Reference Kelman, Russel, Hickson, Smellie and ChapmanKelman and others, 2002; Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002; Reference Höskuldsson, Sparks and CarrollHöskuldsson and others, 2006). These conditions may change very rapidly, even in non-volcanic areas (e.g. Reference Kavanaugh and ClarkeKavanaugh and Clarke, 2000) and are determined by patterns of melting, ice deformation and meltwater drainage (e.g. Reference HookeHooke, 1984; Reference BjörnssonBjörnsson, 1988). Since the melting rate during subglacial eruptions is itself determined by the eruption mechanisms (e.g. Reference Guðmundsson, White, Smellie and ClagueGuðmundsson, 2003), the eruption mechanism will in turn influence cavity conditions, and the volcano–ice system is coupled (Reference Tuffen, Gilbert and McGarvieTuffen and others, 2001). The models presented here constitute a first attempt to simulate the behaviour of this coupled system.

Existing Models of Subglacial Eruption Mechanisms

Many useful insights have come from recent models of subglacial eruption mechanisms. Reference Höskuldsson and SparksHöskuldsson and Sparks (1997) considered pressure changes during effusive eruptions within isolated cavities, from which meltwater cannot escape (closed system). Their calculations showed that the system volume will inevitably increase during effusion of rhyolitic magma, as the volume of ice melted is insufficient to accommodate the erupted magma. Volume increases were also predicted to occur during effusive basaltic eruptions, unless energy transfer was highly efficient (>80%) and the magma was vesicle-poor. The results were used to argue that meltwater may accumulate at the vent during some basaltic eruptions, but cannot do so during any rhyolitic eruption as an increase in volume will lead to increased pressure that will drive meltwater away from the vent area.

Reference Tuffen, Pinkerton, McGarvie and GilbertTuffen and others (2002a) modelled localized melting of a glacier base during an intrusive rhyolite eruption. The competing effects of melting and inward ice deformation on the evolving size of the cavity were modelled using Nye’s law (Reference NyeNye, 1953). Reference Tuffen, Pinkerton, McGarvie and GilbertTuffen and others (2002a) assumed that melting occurs evenly over the roof of a hemispherical cavity and that heat loss from magma was the rate-limiting step in the transfer of thermal energy from magma to ice, rather than the melting process itself. The validity of the latter assumption is brought into question by observations of the Gjálp eruption, as heated meltwater may have left the vent area (Reference Guðmundsson, Sigmundsson, Björnsson and HögnadóttirGuðmundsson and others, 2004).

Reference Wilson, Head, Smellie and ChapmanWilson and Head (2002) provided a novel approach to the mechanisms of subglacial basaltic eruptions, suggesting that dykes may propagate into the glacier base, melting narrow, blade-like cavities, and that the early stages of eruptions may be regarded as the intrusion of sills at the ice-bedrock interface. The dykes would collapse once the supporting ice was melted away.

The Influence of Confining Pressure on Subglacial Eruption Mechanisms

Subglacial basaltic tuyas typically consist of basal pillow lavas overlain by increasingly vesicular hydroclastic tephra (e.g. Reference Smellie, Sigurdsson, Houghton, Rymer, Stix and McNuttSmellie, 2000), indicating that the mechanism of magma–water interaction changes during edifice growth as the confining pressure decreases. In general, pillow lavas generally form beneath ice ≥500m thick and hydroclastic tephras beneath thinner ice (e.g. Reference Moore, Hickson and CalkMoore and others, 1995), but there are exceptions to this rule, as other factors such as the magma volatile content, effusion rate and magma–water ratio are also important (e.g. Reference McBirneyMcBirney, 1963; Reference Wright, Gamble and ShaneWright and others, 2003; Reference Griffiths and FinkGriffiths and Fink, 1992; Reference WohletzWohletz, 1983). Explosive activity may also occur under thicker ice if subglacial pressures are less than glaciostatic (Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002; Reference Schopka, Gudmundsson and TuffenSchopka and others, 2006).

In contrast, there is no obvious relationship between stratigraphic position and facies type for the subglacial phase of rhyolite tuya-building eruptions. Eruptions appear to be either dominantly explosive, generating phreatomagmatic tephra, or intrusive, with the emplacement of quench hyaloclastite (Reference Tuffen, Gilbert and McGarvieTuffen and others, 2001; Reference Tuffen, McGarvie, Gilbert, Pinkerton, Smellie and ChapmanTuffen and others, 2002b; Reference StevensonStevenson, 2005). It thus appears that factors other than confining pressure are of predominant importance during rhyolitic eruptions, most probably the effusion rate and magma volatile content.

Cavity Pressure, Meltwater Drainage and Ice Deformation

Reference Höskuldsson and SparksHöskuldsson and Sparks (1997) showed that cavity pressure in a closed system may be less than glaciostatic (with the difference here defined as the underpressure) if magma–ice energy transfer is highly efficient. Although such efficiency may seldom occur in real systems (Reference Guðmundsson, White, Smellie and ClagueGuðmundsson, 2003), underpressures may nonetheless develop if hydraulic connectivity with low-pressure cavities is established and meltwater drainage occurs (e.g. Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002; Reference Schopka, Gudmundsson and TuffenSchopka and others, 2006). Indeed, there is increasing evidence that meltwater can drain subglacially at the onset of most subglacial eruptions (e.g. Reference Tuffen, Gilbert and McGarvieTuffen and others, 2001; Reference Guðmundsson, Sigmundsson, Björnsson and HögnadóttirGuðmundsson and others, 2004).

Cavity underpressures depend on aspects of the local hydrology such as tunnel length, meltwater flux and temperature (e.g. Reference HookeHooke, 1984). Underpressure leads to inward roof deformation, which will be ductile if cavity radius≪ice thickness. The deformation rate for cylindrical or spherical cavities is given by Nye’s law (Reference NyeNye, 1953). Roof closure leads to depression of the ice surface above, forming an ice cauldron. As the cavity size increases, the ice eventually fails brittly, typically once the cauldron is ≥50m deep (as observed in Iceland, e.g. Reference Guðmundsson, Sigmundsson and BjörnssonGuðmundsson and others, 1997). The depth of brittle failure then increases as the cavity grows until the cavity roof fails and a subaerial eruption commences (Reference Guðmundsson, Sigmundsson and BjörnssonGuðmundsson and others, 1997, 2004).

Evidence from the 1996 Gjá Lp Eruption and Models of Guðmundsson

Observations and models of the 1996 Gjálp eruption, as described in Reference Guðmundsson, White, Smellie and ClagueGuðmundsson (2003) and Reference Guðmundsson, Sigmundsson and BjörnssonGuðmundsson and others (1997, Reference Guðmundsson, Pálsson, Björnsson, Högnadóttir, Smellie and Chapman2002, Reference Guðmundsson, Sigmundsson, Björnsson and Högnadóttir2004) have provided important new constraints on many aspects of magma–ice interaction. The following information is useful for the models presented here: (1) Deformation of the 450–650m thick ice above the eruption site was initially ductile. (2) Although meltwater drained from the eruption site, a column of meltwater accumulated above the edifice. (3) The well-drained subglacial cavity above the southern part of the fissure had an underpressure of about 2MPa. (4) The heat transfer efficiency between magma and ice at the eruption site f _i was between 0.55 and 0.66. (5) Less than 10% of the erupted products were washed away by meltwater. (6) The heat flux during the eruption was 100 times that during slow cooling of the Gjálp edifice.

Outline of the Model

This paper employs a simple model of the volcano–ice system to explore the effects of various parameters on eruption mechanisms. It considers rhyolitic eruptions as point sources at the base of horizontal ice sheets with constant effusion rates. It is assumed that the feeder dyke does not overshoot the bedrock–ice interface (cf. Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002). Magmatic heating causes melting of ice, forming a hemispherical cavity at the glacier base (Fig. 1). The cavity roof closes by ductile deformation, approximated by Nye’s law, at a rate determined by the cavity size and underpressure. The cavity size is determined by the relative rates of enlargement by melting and closure by ice deformation. Meanwhile, a volcanic edifice accumulates within the cavity (Fig. 1). Although meltwater is assumed to escape from the cavity, any gap above the edifice will be filled with water and steam (Fig. 1a). The model predicts the evolving sizes of the cavity and the edifice together with the depth of water/steam above the edifice and the depth of the ice cauldron (Fig. 2). A range of ice thicknesses, effusion rates and magma vesicularities is employed. Tables 1 and 2 list the variables and constants used.

Fig. 1. Representations of the (a) explosive and (b) intrusive styles of subglacial rhyolite eruptions modelled in this paper. There is a substantial gap above the edifice in (a), filled by meltwater or steam allowing explosive magma–water interaction to occur. In (b) the edifice completely fills the cavity and rising magma intrudes within the edifice. Modified from Reference Tuffen, Gilbert and McGarvieTuffen (2001).

Fig. 2. Model results for (a) large and (b) small effusion rate eruptions beneath ice 450m thick, indicating the changing elevation of the edifice radius, cavity radius and, in (a) only, the ice surface.

Table 1. Variables and constants used in the models. Data from Reference Höskuldsson and SparksHöskuldsson and Sparks, 1997

Table 2. Properties of ice and rhyolite magma (Reference Höskuldsson and SparksHöskuldsson and Sparks, 1997)

Heat transfer from magma to ice

The thermal energy released per unit volume of rhyolitic magma during cooling and quenching E is given by

(1)

where ϴ denotes magma vesicularity, ρ _l the unvesiculated magma density, c _l magma heat capacity and ΔT _l the temperature change of the magma (Reference Höskuldsson and SparksHöskuldsson and Sparks, 1997). This assumes that (1) the deposits are glassy, which is appropriate for the vast majority of subglacial rhyolites studied, so that the latent heat of crystallisation can be ignored, and (2) the heat capacity of gases within vesicles is negligible.

Mechanisms and efficiency of heat transfer

Only heat released rapidly from the magma is considered (within seconds or minutes of its emplacement); later edifice cooling is ignored. The heat transfer efficiency f _i is the fraction of heat transferred instantaneously to melt the ice

(Reference Guðmundsson, Sigmundsson, Björnsson and HögnadóttirGuðmundsson and others, 2004). Therefore, the heat flux for ice melting, Q _m, is equal to

(2)

where Ql_(t) is the instantaneous effusion rate and E is the total thermal energy of the magma as before. The heat transfer efficiencies for two possible eruption mechanisms are considered. During (a) phreatomagmatic fragmentation, heat loss from magma is rapid due to the large surface area for heat exchange (Reference Zimanowski, Buttner and LorenzZimanowski and others, 1997). Heat loss models and observations from Gjálp suggest that f _i may range between 0.55 and 0.66 (Reference Guðmundsson, White, Smellie and ClagueGuðmundsson, 2003), but for simplicity and because these models are only approximate, a value of 1 will be used for these models. During (b) intrusive eruptions, rising magma is quenched by meltwater within the slowly cooling edifice (Reference Tuffen, Pinkerton, McGarvie and GilbertTuffen and others, 2002a). The heat transfer rate is much smaller than during phreatomagmatic fragmentation and is similar to that from a newly erupted edifice. As the heat flux from the Gjálp edifice immediately after the eruption was approximately 5–10% of the heat flux during it (Reference Guðmundsson, Sigmundsson, Björnsson and HögnadóttirGuðmundsson and others, 2004), f _i is estimated to be 0.1.

Melting of ice

Magmatic heat may cause melting of ice and also heating and vaporization of the meltwater. Although meltwater temperatures may be 208C or more (Reference Guðmundsson, Sigmundsson and BjörnssonGuðmundsson and others, 1997), for simplicity meltwater will be assumed to be at 08C. The fraction of heat converted to kinetic energy during phreatomagmatic explosions is likely to be small (Reference WohletzWohletz, 1983) and will be ignored. The volume of ice melted per unit volume of magma v _i is thus given by

(3)

where ρ _i and L _i are the density and heat of fusion of ice, respectively.

It is assumed that melting occurs evenly over the roof of a hemispherical cavity of radius R. The melt-back rate r0_mis thus given by

(4)

where v _i is the volume of ice melted per unit volume of magma and Q _l is the effusion rate.

Ice deformation

According to Nye’s law, the roof of a hemispherical cavity within an infinite ice sheet will deform at a rate r0_d given by

(5)

where ΔP is the cavity underpressure, n is Glen’s flow law constant (3) and B is an ice deformation parameter, taken as 5.3×10⁷ Pa s^1/3 (e.g. Reference HookeHooke, 1984).

Nye’s law applies when the cavity radius is small compared with the ice thickness (Reference NyeNye, 1953) but breaks down at large cavity sizes. The simulations have been run until the ice roof is less than 50m thick, when it is assumed to fail. The results will therefore only be approximate at large cavity sizes, once Nye’s law begins to break down.

Ice will flow horizontally into the ice cauldron, as its surface seeks to regain an equilibrium shape (Reference Aðalgeirsdóttir, Guðmundsson and BjörnssonAðalgeirsdóttir and others, 2000). At Gjálp, deformation speeds were approximately 10^–5–10^–6ms^–1, meaning that substantial changes in the surface profile would require years. To illustrate this, at a reasonable range of effusion rates of 1–100m³ s^–1 for silicic dome-building eruptions and 10³– 10⁶m³ s^–1 for plinian eruptions (Reference Pyle, Sigurdsson, Houghton, Rymer, Stix and McNuttPyle, 2000), we find that a 500m high hemispherical edifice would be constructed in months to years during a dome-type eruption and hours to days in plinian eruptions. This suggests that horizontal closure beneath ice thicknesses of the order of 500 m may only be significant during prolonged eruptions with small effusion rates. Nonetheless, two end member scenarios are envisaged in the models: one in which no horizontal ice deformation occurs, the other where closure is so rapid that the surface remains horizontal at all times (here termed ‘perfect horizontal deformation’).

Cavity pressure and roof closure

For the purposes of this study, cavities are assumed to be at atmospheric pressure due to meltwater drainage and the establishment of a hydraulic connection with low-pressure cavities. Although this may be unrealistic for some eruptions, it avoids the difficulty of estimating cavity pressures, which depend on so many aspects of the local hydrology and the melting process. A justification for this assumption is that many rhyolitic eruptions occur at central volcanoes with powerful geothermal systems (e.g. Torfajökull, Krafla and Kerlingarfjöll). There is therefore likely to be a high geothermal heat flux in the vicinity of the vent directly preceding the eruption, which would favour the development of low-pressure cavities, as found in the firn-filled crater of Mount Rainier where there is a strong geothermal field (Reference Kiver and SteeleKiver and Steele, 1975). The rates of ice deformation predicted by the model are therefore maximum rates. Models exploring the effects of different cavity pressures on the mechanisms of subglacial eruptions will be presented elsewhere.

The glaciostatic pressure at the cavity roof P _gðtÞ at time t is given by

(6)

where h _{r(t )} is the ice roof thickness at time t, expressed as

(7)

where h ₀ is the initial ice thickness and r0_mis the melt-back rate as defined in Equation (4). The pressure driving deformation of the cavity roof, P _{(t )} is given by

(8)

where P _c is the cavity pressure, assumed to be 0.1 MPa in these simulations. Note that when there is no horizontal deformation, the roof will progressively thin during the eruption, reducing the weight of the ice roof and leading to slower rates of closure.

Cavity size balance

The cavity size is determined by the relative rates of enlargement by melting and closure by deformation. The rate of cavity enlargement dR/dt is thus given by the relation:

(9)

Substituting Equations (4) and (5) into Equation (9) yields

(10)

where Q _l is the volume flux of magma (m³ s ^–1). Equations (3) and (8) can therefore be combined with Equation (10) to give

(11)

which provides the changing cavity size with time. Solutions to Equation (11) are approximated numerically using a spreadsheet program. Values of melt-back rate, ice deformation rate and cavity radius are calculated incrementally for a series of time steps such that t _nþ1 = 1:01t _n. Simulations are continued until the roof is less than 50m thick, which typically takes 106–10⁷ s. Roof failure will then lead to a subaerial eruption. The cavity radius initially increases rapidly (dR/dT ∽ 10^–2–10^–3ms^–1), because the melting occurs over a small area and roof closure is slow. As the cavity size increases, the melt-back rate reduces and roof closure accelerates, although substantial roof thinning will reduce the closure rate.

Edifice volume and space in cavity

The edifice volume V _e is given by

(12)

where S is a packing coefficient, taken to be 1 for intrusive lava and 2 for pyroclastic deposits. It is assumed that a negligible amount of material is removed by meltwater, justified by the estimate that only 10% of the Gjálp products were washed away (Reference Guðmundsson, Sigmundsson, Björnsson and HögnadóttirGuðmundsson and others, 2004). The edifice radius is thus given by

(13)

Two scenarios are envisaged, one in which the edifice is considerably smaller than the cavity (A), and one in which the cavity is completely filled by the edifice (B).

A. Edifice radius < cavity radius

There will be a significant water column above the growing edifice, meaning that rising magma can interact with abundant meltwater. This allows potentially explosive magma–water interaction to occur (Figs 1a and 2a). The development of convection cells in the water may assist transfer of heat from fragmenting magma to ice, leading to rapid melting.

B. Edifice fills cavity

If the edifice fills the cavity, rising magma cannot interact with a column of water (Figs 1b and 2b). Instead, it is likely to intrude and be quenched within the poorly consolidated, waterlogged edifice. This may generate facies similar to pumiceous peperites (Reference Hunns and McPhieHunns and McPhie, 1999; Reference Tuffen, Gilbert and McGarvieTuffen and others, 2001), where the explosivity of the magma–water interaction is suppressed by the surrounding fragmental material, despite the magma being volatile-rich. A small proportion of the magma may reach the glacier base, where it may generate small ice caves through localized melting (Reference Tuffen, Pinkerton, McGarvie and GilbertTuffen and others, 2002a) or intrude as a sill at the bedrock– ice interface (Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002).

Results of the Models

The model calculates the radii of the growing edifice and evolving subglacial cavity, together with the ice surface elevation. Figure 2 shows typical graphs produced using the model, for two values of effusion rate beneath ice 450m thick. Deformation is insignificant when the effusion rate is large. Melting is more rapid (Fig. 2a) than roof closure and a considerable water- or steam-filled gap develops above the edifice until roof failure, when the cauldron is ∽70m deep. The smaller effusion rate leads to slower melting, and ice deformation becomes more important (Fig. 2b). With time, the edifice begins to fill the cavity (see arrow on Fig. 2b), leading to an intrusive style of eruption. Once cavities become filled, it is assumed that eruptions will remain intrusive as the heat flux to the glacier is greatly reduced. The evolving ratio of edifice radius to cavity radius for a range of effusion rates is displayed in Figure 3a, illustrating how cavities become filled at small effusion rates but a considerable meltwater/steam-filled gap develops at large effusion rates. The effects of magma vesicularity are shown in Figure 3b. If all other parameters are unchanged, increasing the magma vesicularity can lead to the cavity becoming filled at intermediate effusion rates. This is due to the smaller amount of heat per unit volume carried by highly vesicular magma.

Fig. 3. Ratio of edifice radius to cavity radius during eruptions. Once the ratio becomes equal to 1, the eruption becomes intrusive and is assumed to remain so. (a) Results for eruptions with a range of effusion rates (Q _l values indicated) beneath ice 450m thick. X indicates when the roof is <50m thick and fails, so that eruptions become subaerial. (b) Results for variably vesicular magma under ice 500m thick (% vesicularity indicated). When the vesicularity is small, there is sufficient melting to maintain space above the edifice and explosive eruptions occur. Greater vesicularities lead to filling of the cavity and the transition to intrusive eruptions.

Will eruptions be explosive or intrusive?

Solutions were found for the critical effusion rate Q _cr for a given ice thickness and with either no or perfect horizontal ice deformation. At rates less than Q _cr, the cavity becomes completely filled by the edifice, which is likely to lead to an intrusive eruption mechanism. Above Q _cr, the edifice never completely fills the cavity, and the eruption has the capacity to remain explosive until it becomes subaerial. This process was repeated for a range of ice thicknesses, and the results plotted on Figure 4. Intrusive eruptions are favoured by thick ice and small values of Q _l. In addition to the main part of the modelling, it has also been shown that high meltwater temperatures and low cavity pressures favour intrusive eruptions (Reference Tuffen, Gilbert and McGarvieTuffen, 2001), although these parameters are themselves dependent on the melting rate and eruption mechanisms and so are not strictly independent.

Fig. 4. The critical effusion rate Q _cr plotted as a function of ice thickness. The two curves represent the two end-member scenarios: no and perfect horizontal ice deformation. The shaded fields illustrate plausible effusion rates and ice thicknesses for subglacial rhyolite eruptions at Torfajökull: B-Bláhnúkur, D-Dalakvísl, R-SE Rauðufossafjöll. Modified from Reference Tuffen, Gilbert and McGarvieTuffen (2001).

Mixed explosive-intrusive eruption: Dalakvísl

A subglacial eruption at Dalakvísl, Rauðufossafjöll, was intermediate in style, including both explosive and intrusive phases (Tuffen and others, in press) (Fig. 5). The explosive phase initially involved fragmentation of >70% vesicular magma within a subglacial cavity, but this gave way to foam collapse and the rising magma began to be trapped and quenched within the vesicular deposits. The relatively high dissolved water content (>0.5 wt%) measured using Fourier Transform Infrared (FTIR) spectroscopy on samples of the intruded lava indicates that this transition was not due to the exhaustion of magmatic volatiles (Tuffen and others, in press). Instead, it appears that the poorly-consolidated tephra had filled the cavity and was intruded by partly-fragmented foam, which no longer had the space to fragment or opportunity to interact explosively with meltwater. The highly vesicular nature of the magma may have encouraged the cavity to fill, due to the small amount of thermal energy carried per unit volume of the magma (Fig. 3b).

Fig. 5. An illustration of the mechanisms of the subglacial rhyolite eruption at Dalakvísl, Torfajökull (after Tuffen and others, 2007). In one location, an initially explosive eruption within a subglacial cavity became intrusive, apparently when the cavity was filled with tephra. Elsewhere, intrusive lava lobes were emplaced within waterlogged hyaloclastite breccias, some quenching at the glacier base. This may reflect changes in the effusion rate during the eruption. Meltwater drained away subglacially during the eruption.

Discussion

Our model suggests that the effusion rate and ice thickness are important parameters influencing the mechanisms of subglacial rhyolite eruptions. The effusion rates required for intrusive and explosive eruptions at Torfajökull are consistent with independent estimates. Viewed broadly, the explosive eruptions are equivalent to plinian eruptions and the intrusive eruptions to effusive eruptions, although the elevated volatile contents in intrusive facies (Tuffen and others, in press) suggest that these contrasting styles are principally due to contrasting cavity conditions rather than different degrees of degassing. This contrasts with basaltic eruptions, where cavity pressures appear to be more important (e.g. Reference Guðmundsson, Sigmundsson, Björnsson and HögnadóttirGuðmundsson and others, 2004).

An explanation may be related to the initial volatile content of the magma. Whereas modest (∽5MPa) changes in cavity pressure may significantly affect the vesiculation and ascent of water-poor basaltic magma (Reference Höskuldsson, Sparks and CarrollHöskuldsson and others, 2006), rhyolitic magma may be so water-saturated that the bulk of vesiculation occurs beneath the glacier base, reducing the impact of such pressure changes on its rate of ascent. This highlights the need to understand the relationship between volatile degassing, effusion rate and subglacial eruption mechanisms more fully, and this will be investigated in future work.

Alternatively, the lack of palagonitisation in subglacial rhyolitic deposits may explain the difference, meaning that they remain loose and waterlogged even during sustained eruptions. This makes it more difficult for rising magma to reach the glacier base: instead it is trapped and quenched within the edifice. In contrast, basaltic hyaloclastite rapidly palagonitises (Reference JakobssonJakobsson, 1978), possibly fast enough to cause cementation during eruptions and facilitate propagation of dykes through to the glacier base.

Regarding glacier hydrology and cavity pressures, it is important to remember that the models developed here are simplistic, and represent only a first step towards the formulation of a general model of subglacial rhyolite eruptions. In addition to magma volatile contents, other issues to further consider include cavity morphologies (Reference Höskuldsson and SparksHöskuldsson and Sparks, 1997; Reference Wilson, Head, Smellie and ChapmanWilson and Head, 2002) and how cavity pressures are influenced by the glacier and bedrock topography, melting rate and the existing subglacial cavity system (e.g. Reference BjörnssonBjörnsson, 1988; Reference Guðmundsson, Sigmundsson, Björnsson and HögnadóttirGuðmundsson and others, 2004; Reference Höskuldsson, Sparks and CarrollHöskuldsson and others, 2006).

Finally, our models predict that ice cauldron depths may depend on the effusion rate, with the deepest cauldrons formed in small effusion rate eruptions when there is most time for ice surface deformation. Although these models only poorly simulate the formation of ice cauldrons, as brittle processes are clearly important, they do raise the intriguing possibility that apparent palaeo-ice thicknesses, as recorded by the ash-lava transition at tuyas, may depend upon the effusion rate, making palaeoclimatic reconstructions more difficult. This highlights the need for more attempts to quantitatively model ice deformation and melting during subglacial eruptions, and to improve our understanding of past climate as well as the hazards from explosive eruptions and jökulhlaups.

Conclusions

The results of simple models of subglacial rhyolite eruptions suggest that their mechanisms are influenced by the effusion rate and ice thickness. At small effusion rates, subglacial cavities may become completely filled with volcanic debris as roof closure is important, favouring an intrusive style of eruption. At larger eruption rates, there is always space in the cavity for an explosive eruption to occur, due to the faster rate of melting. The models are used to estimate the effusion rate during some Quaternary eruptions at Torfajökull and the values found are consistent with independent estimates. This has been a first attempt to model the feedbacks between eruption mechanisms, ice melting and ice deformation during subglacial eruptions, and illustrates the challenge the complexly coupled volcano–ice system poses for future modellers.