Formation and persistence of glaciovolcanic voids explored with analytical and numerical models

Abstract

One fifth of Earth's volcanoes are covered by snow or ice and many have active geothermal systems that interact with the overlying ice. These glaciovolcanic interactions can melt voids into glaciers, and are subject to controls exerted by ice dynamics and geothermal heat output. Glaciovolcanic voids have been observed to form prior to volcanic eruptions, which raised concerns when such features were discovered within Job Glacier on Qw̓elqw̓elústen (Mount Meager Volcanic Complex), British Columbia, Canada. In this study we model the formation, evolution, and steady-state morphology of glaciovolcanic voids using analytical and numerical models. Analytical steady-state void geometries show cave height limited to one quarter of the ice thickness, while numerical model results suggest the void height h scales with ice thickness H and geothermal heat flux $\dot {Q}$ as $h/H = a H^b \dot {Q}^c$, with exponents b = −n/2 and c = 1/2 where n is the creep exponent. Applying this scaling to the glaciovolcanic voids within Job Glacier suggests the potential for total geothermal heat flux in excess of 10 MW. Our results show that relative changes in ice thickness are more influential in glaciovolcanic void formation and evolution than relative changes in geothermal heat flux.

1. Introduction

A significant portion of Earth's volcanoes are currently directly overlain by, or in near proximity to, glaciers or perennial or seasonal snow/firn (e.g., Curtis and Kyle, 2017; Edwards and others, 2020). Ice-clad volcanoes pose increased risk to societies and infrastructure compared to many of their ice-free counterparts through processes such as jökulhlaups (e.g., Björnsson, 2010), lahars (e.g., Pierson and others, 1990) and heightened explosivity of eruptions due to magma–ice interaction (e.g., Sigmundsson and others, 2010). Heat fluxes from volcanoes are capable of significantly influencing glacier mass balance and dynamics on regional (e.g., Rogozhina and others, 2016; Barr and others, 2018; Burton-Johnson and others, 2020; Jóhannesson and others, 2020) and local (e.g., Björnsson, 1975; Barr and others, 2018) scales, sometimes forming distinct features such as ice cauldrons, subglacial lakes and caverns (e.g., Barr and others, 2018). Variations in geothermal heat output, and consequent morphological changes of these glaciovolcanic features have been observed as precursory activity before subglacial eruptions (Bleick and others, 2013; Reynolds and others, 2018).

In 2016, localised subglacial geothermal activity was observed to have caused bed-to-surface melt-through of Job Glacier on Qw̓elqw̓elústen or the Mount Meager Volcanic Complex, British Columbia, Canada (e.g., Williams-Jones, 2016; Roberti and others, 2018) (Fig. 1). Three distinct features were documented at the time: two where complete melt-through had occurred and active venting of geothermal gases was observed (labelled 1a and 2a in Fig. 1), and one localised area of surface subsidence suspected to be a result of subglacial geothermal activity (labelled 3a in Fig. 1) (see Supp. Mat. for observational details). These glaciovolcanic features are all believed to be formed by fumaroles, i.e., vents that emit geothermal or volcanic gases. The discovery raised the question whether the geothermal features had simply been hidden beneath the glacier prior to 2016, or if the appearance of these features could be interpreted as a manifestation of increased volcanic activity within the system. Roberti and others (2018) hypothesised that glacier thinning was the cause, rather than an increase in volcanic activity.

Figure 1. A map of Qw̓elqw̓elústen or the Mount Meager Volcanic Complex in British Columbia, Canada. The red star indicates the site of the glaciovolcanic activity within Job Glacier. Shaded glacier areas are from 1985 to 1986 (GeoBC, 2010), and white glacier areas are from 2004–2006 (RGI Consortium, 2017). Dark contours are m.a.s.l., and coordinates are given in UTM Zone 10 North. Photographs show the evolution of the glaciovolcanic voids within Job Glacier, seen from helicopters looking south. The glaciovolcanic voids are enumerated based on their inferred subglacial geothermal source. Labels are further subdivided with letters, in order of appearance, to represent voids that share the same inferred geothermal source. Labels in parentheses indicate inactive or extinct voids. Photo credits: G. Roberti (2016), G. Williams-Jones (2019), and T. Unnsteinsson (2020, 2021).

This problem of distinguishing between possible causes, however, demonstrates the current lack of established relationships between glaciological and volcanological conditions that allow for the formation and persistence of these glaciovolcanic features. Previous analytical work has focused on the formation of voids during subglacial eruptions (Höskuldsson and Sparks, 1997; Tuffen and others, 2002; Tuffen, 2007), while numerical studies have either been highly case specific (Curtis, 2016) or focused on water-filled voids (e.g., Einarsson and others, 2017). In this paper we aim to model the formation and evolution of water-drained glaciovolcanic voids in a holistic manner using analytical and numerical models. We start by combining existing analytical models of symmetrical void closure rates with the induced geothermal melt rates, to derive steady-state void geometries. We also employ a full-Stokes numerical ice-flow model to further investigate steady-state void geometries and melt-through occurrence, with a more detailed treatment of ice dynamics. Insight from both the analytical and numerical models is then applied to Job Glacier to estimate the subglacial geothermal heat fluxes.

1.1 Glaciovolcanic caves and chimneys

One mechanism for the transfer of geothermal heat is the release of volatiles from volcanic and hydrothermal systems through various processes, such as magma degassing, boiling of water, and evaporation (e.g., Fischer and Chiodini, 2015). When volatile release in a subnivean or subglacial environment is sufficiently intense that ice dynamics and accumulation can not fully compensate for the melted volume, and if the meltwater produced by this glaciovolcanic interaction is effectively drained away, then voids may form in the overlying snow, firn or ice (e.g., Barr and others, 2018). We define glaciovocanic voids as subnivean and sub- or englacial macroscopic voids that are not saturated with water and are formed, kept open by or accumulate volcanic or geothermal gases. We use the term void instead of terms like cavity to distinguish it from the terminology used for subglacial hydrology. The size and morphology of these glaciovolcanic voids is determined by the rate and magnitude of the heat output from the geothermal system, the manner and efficiency of heat transfer and the dynamics and rheological properties of the overlying snow, firn or ice, along with meteorological forcings (Zimbelman and others, 2000; Curtis and Kyle, 2011).

1.1.1 Classification of glaciovolcanic voids

A spectrum of glaciovolcanic voids exists that depends on the geothermal activity and the magnitude of effects it has on the overlying ice. Here we describe several distinct examples of glaciovolcanic voids based on the degree of glaciovolcanic interaction and the resulting void morphologies. Passive glaciovolcanic voids are voids in ice or snow that are neither formed nor maintained by geothermal activity, though still contain volcanic/geothermal gases. These voids can pose a significant risk to humans (e.g., Hill, 2000; Hansell and Oppenheimer, 2004; Williams-Jones and Rymer, 2015), but are not investigated further here.

Voids that form, or are kept open, primarily as a result of geothermal/volcanic gas influx are here termed active. These voids can form either by direct melting or by exploitation of existing weaknesses or pathways within the ice, where growth or maintenance of the voids is a direct result of the influx of gases. Note that evaporation of geothermal waters at the ice–bed interface is also considered a vapour influx. We further categorise active glaciovolcanic voids into two sub-categories: caves, which are (near) horizontal or ground-hugging, and chimneys, which are (near) vertical (Fig. 2).

Figure 2. (a) Schematic illustration of a glaciovolcanic cave (horizontal) and chimney (vertical) within a glacier of uniform thickness (blue arrow indicates flow direction), above bedrock that hosts a hydrothermal system (multicoloured arrows). (b) A glaciovolcanic chimney within Job Glacier in the Mount Meager Volcanic Complex, British Columbia, Canada (2b in 2020 photo in Fig. 1). The chimney formed as a result of subglacial fumarolic activity and is ~15 m wide at the glacier surface. (c) A glaciovolcanic cave, roughly 10 m in diameter, beneath Kverkjökull, Iceland, formed due to presence of geothermal waters in the subglacially originating stream Volga. Photos (b) and (c) were taken by T. Unnsteinsson.

Glaciovolcanic caves form when the supplied geothermal heat flux is sufficient to melt a void into the overlying snow, firn or ice, and compensate for any influx of ice, but not enough to completely melt through. Due to the presence of permeable debris at the bed (or possibly fractured bedrock) in volcanic settings (Brugman and Meier, 1981; Walder and others, 2005), hydrothermal or volcanic gases from the source may be able to percolate along the ice–bed interface where the associated heat slowly melts an elliptical or semi-cylindrical passageway, i.e., a cave. The term glaciovolcanic cave is adopted here as it is more descriptive and inclusive of features variously referred to as glacier/ice/firn cave (Kiver and Mumma, 1971; Lyon and Giggenbach, 1973; Kiver and Steele, 1975; Wardell and others, 2003; Cousins and others, 2013; Oddsson, 2016; Pflitsch and others, 2017), volcanic ice-cave (Zimbelman and others, 2000), geothermal firn/ice cave (Giggenbach, 1976; Anderson and others, 1998; Anderson and Vining, 1999), fumarole/fumarolic ice cave (Curtis and Kyle, 2011; Ilanko and others, 2019; Florea and others, 2021) and subglacial melt cavity (Barr and others, 2018). The term glaciovolcanic cave has now also been adopted in recent literature (e.g., Sobolewski and others, 2022).

If the melt opening due to influx of geothermal heat is larger than ice closure, vertical growth of a void is sustained. Convection of the geothermal gases and decreasing closure rates with increasing void height promote vertical void expansion until melt-through occurs at the glacier surface, forming a vertical shaft through the ice. Smellie (2002) refers to a feature of this kind, formed by fumarole activity, as a chimney to conform with the terminology used by Gudmundsson and others (1997). However, Gudmundsson and others (1997) and Jakobsson and Gudmundsson (2008) use the word chimney to describe the opening in the glacier once the surface of an ice cauldron, above a subglacial eruption, collapses. Other terms that have been used to describe this phenomenon and the associated features are breach, vent, perforation or steaming crevasse (Eichelberger and others, 1976), ice pits (Malone and Frank, 1975; Friedman and Frank, 1980), melt pit/hole or ice crater (Barr and others, 2018) and ice perforation (Friedman and Frank, 1980). We use the term glaciovolcanic chimneys to describe these features, to be inclusive of different glaciovolcanic processes (e.g., subglacial volcanic eruptions) that lead to their formation.

1.1.2 Conceptual model of the glaciovolcanic features within Job Glacier

The glaciovolcanic voids within Job Glacier consist of both caves and chimneys (see Supplementary Material for observational details). To explain the formation and evolution of the glaciovolcanic voids within Job Glacier, general physical relationships between void geometries and the relevant glaciovolcanic parameters are required. That is, the effects of these parameters on creep closure and melt opening must be explored. We thus model glaciovolcanic voids generally, rather than attempting to represent the peculiarities of any individual void, and then use the model results to constrain the parameters applicable to naturally occurring voids. Our analytical and numerical models apply to both glaciovolcanic caves and chimneys, the latter of which can be open at the glacier surface as observed at Job Glacier. To simplify the analysis, we assume that: (i) heat transfer is efficient, with most of the geothermal energy going directly into melting and negligible amounts escaping; (ii) all melt water is efficiently drained away subglacially; (iii) energy-depleted geothermal gases do not accumulate within voids but escape without otherwise affecting the system; and (iv) glaciovolcanic voids are independent from each other.

2. Methods

Given the observationally based conceptual model of the glaciovolcanic features within Job Glacier as described above, we carry out several steps to investigate these features. We use an analytical model to derive: (i) the maximum height of glaciovolcanic caves; (ii) the total heat flux associated with idealised glaciovolcanic chimneys; (iii) a physics-based conceptual model of chimney evolution. To further investigate the influence of background flow fields we use a numerical model to: (iv) simulate the formation of glaciovolcanic voids; and (v) relate steady-state void geometries to glaciological and volcanological parameters. The results of these steps are then used to infer the geothermal heat fluxes responsible for the glaciovolcanic voids within Job Glacier.

2.1 Analytical modelling

To model glaciovolcanic voids there are two processes that need to be considered: the closure of the void as a result of glacier dynamics, and the geothermally induced melt rate. Glacier ice can flow as a result of three mechanisms: basal slip, bed deformation and ice creep (e.g., Cuffey and Paterson, 2010). The first two mechanism are important for macro-scale ice-flow, but for this analysis we consider only ice creep, i.e., the internal deformation of the ice mass, and assume that any movement at the base is entirely modulated by ice creep. Ice creep can be described by a constitutive relation that relates the deformation of the ice to the applied stresses. For isotropic ice, this relationship is described as a power law:

(1)$$\dot{\epsilon} = A \tau^n,\; $$

where $\dot {\epsilon }$ is the effective strain rate, τ is the effective stress, A is the creep parameter, which depends on the physical properties of the ice, and n is the creep exponent (Glen, 1955). We adopt commonly used values of the creep parameter and exponent for temperate ice: A = 2.4 × 10⁻²⁴ Pa⁻³ s⁻¹ and n = 3, respectively (Cuffey and Paterson, 2010).

The radial creep closure of radially symmetric subglacial voids, within an infinitely large and uniform ice layer, is

(2)$$-{1\over r} {{\rm d} r\over {\rm d} t} = k A \left({\,p_{\rm i}\over n}\right)^n,\; $$

where r is the radial distance in both cylindrical and spherical coordinates, p _i = ρ _igH is the ice-overburden pressure, $\rho _{\rm i} = {917}{\rm \, kg\, m^{-3} }$ the ice density, $g = {9.81}{\rm \, ms^{-2} }$ the gravitational acceleration, H the ice thickness, and k = 1 or k = 2⁻ⁿ 3^(n−1)/2 for cylindrical or spherical void geometry, respectively (Nye, 1953). Assuming a uniform specific heat flux $\dot {q}$ over the void surface and ice at the pressure-melting point, the radial melt rate is

(3)$${{\rm d} r\over {\rm d} t} = {\dot{q}\over \rho_{\rm i} L},\; $$

where L is the latent heat of fusion for ice. Balancing creep closure Eqn (2) and geothermal melt-opening Eqn (3) gives

(4)$$r = {\dot{q}\over k A \rho_{\rm i} L} \left({n\over p_{\rm i}}\right)^n,\; $$

which describes the steady-state radial geometry of a void as a function of the ice-overburden pressure and specific heat flux. Equation (4) is analogous to the formulations used by Tuffen and others (2002), Tuffen (2007) and Höskuldsson and Sparks (1997) to model subglacial volcanic eruptions.

2.1.1 Glaciovolcanic caves

An analytical solution to an arbitrary, nonradially symmetric void geometry does not exist. Despite the theoretical limitations associated with Eqn (2), the equation is still regarded to hold merit for geometries that strictly violate the assumptions used for its derivation (e.g., Nye, 1976). It has been used to model nonradially symmetric hydrological passageways and cavities (e.g., Walder, 1986; Hooke and others, 1990; Creyts and Schoof, 2009), as well as large (r = H/4) symmetrical voids formed due to subglacial volcanic eruptions (Tuffen, 2007). Similarly to Hooke and others (1990), we employ Eqn (2) to describe the creep closure of nonradially symmetric glaciovolcanic caves by defining an effective void radius that varies along the void surface. Instead of void radius r in Eqn (2), we introduce normal radius r _n (Fig. 3). The normal radius is the radial distance between the cave surface and the z-axis, perpendicular to the cave surface, and is defined as $r_n = \mathbf{r} \;\cdot$ $\hat{\mathbf{n}}$, where $\mathbf{r}$ is the point on the cave surface, and $\hat{\mathbf{n}}$ is the unit normal vector to the cave surface (Fig. 3). We approximate the normal stress as the hydrostatic stress, thus replacing ρ _igH in Eqn (2) with ρ _ig (H − z). Based on these assumptions, and Eqn (2), the creep closure at any point along the surface can then be written

(5)$$-{1\over r_n} {{\rm d} r_n\over {\rm d} t} = k A \left({\rho_{\rm i} g ( H-z) \over n}\right)^n.$$

Note that in the limit r/H → 0, Eqn (5) reduces to Eqn (2) of Nye (1953). Similarly, the melt rate normal to the cave surface is

(6)$${{\rm d} r_n\over {\rm d} t} = {\dot{q}\over \rho_{\rm i} L},\; $$

and the balance between creep closure and melt opening becomes

(7)$$r_n = {\dot{q}\over k A \rho_{{\rm i}} L} \left({n\over \rho_{\rm i} g \left(H - z\right)}\right)^n.$$

The geometries given by Eqn (7), and the associated closure rate Eqn (5), were tested with a numerical ice-flow model and show improved results compared to radially symmetric geometries Eqn (4) with uniform closure rates Eqn (2) (Fig. S1).

Figure 3. Schematic diagram of an arbitrary noncylindrical cave geometry. At any point on the cave surface, described by $\mathbf{r}$, the surface normal at that point, $\hat {\mathbf{n}}$, defines a radial distance between the surface and the vertical z-axis that is normal to the cave surface, given by ${\bf r}_n$. The point of intersection with the ${\cal Z}$-axis is O _n and the azimuth angles of $\mathbf{r}$ and ${\bf r}_n$ are θ and θ _n, respectively.

2.1.2 Glaciovolcanic chimneys

We model chimneys as radially symmetric about a vertical central axis, with the possibility of a depth-dependent radius, leading to variable cross-sectional area. The radial closure of englacial tunnels Eqn (2), and variations thereof, has previously been used to model vertical hydraulic conduits (Röthlisberger, 1972) and moulins (e.g., Catania and Neumann, 2010; Andrews and others, 2022). We thus apply the formulation to glaciovolcanic chimneys, which are geometrically similar vertical features. Allowing a variable cross-sectional area strictly violates the assumptions used to derive the equation for the contraction of cylindrical voids Eqn (2), as nonradial strain rates would be induced along the axis of chimneys. However, realistic chimney geometries should not vary greatly in radius over short distances, and they should have vanishing strain rates at the glacier surface, just as described by Eqn (2). To account for the variable cross-sectional area, we rewrite Eqn (4) as

(8)$$r( z) = {\dot{q}( z) \over k A \rho_{\rm i} L} \left({n\over \rho_{\rm i} g ( H-z) }\right)^n.$$

The total heat flux required to keep a chimney in steady state is therefore

(9)$$\dot{Q} = 2 \pi \rho_{{\rm i}} L A {\int}_{0}^H r( z) ^2 \left({\rho_{{\rm i}} g ( H-z) \over n} \right)^n dz.$$

2.1.3 Glaciovolcanic chimney advection

The above formulations neglect ice dynamics aside from those induced by the chimneys themselves. Such chimneys would only evolve according to variations in the geothermal input and radially symmetric creep closure. To determine the effect of glacier dynamics on chimneys, background flow fields must be considered. Incorporating background flow fields into the analytical chimney models is nontrivial, and is thus not attempted here. Rather, we propose a physics-based conceptual model of the evolution of a chimney. For simplicity let us consider a glacier of constant thickness H, and infinite width and length resting on an inclined bed. The coordinate system is aligned with the glacier: ${\cal X}$ is along the bed in the direction of flow, ${\cal Y}$ is along the bed but perpendicular to ${\cal X}$, and ${\cal Z}$ is normal to the bed. The glacier is subjected to a bed-parallel velocity field, which is chosen to be represented by the shallow-ice approximation:

(10)$$v_{\cal X}( {\cal Z}) = {2 A\over n + 1} \tau_b^n H \left(1 - \left(1 - {{\cal Z}\over H} \right)^{n + 1} \right),\; $$

where τ _b is the basal shear stress, assumed equal to the driving stress (Cuffey and Paterson, 2010). An incipient cylindrical chimney of radius R is present at the centre of the coordinate system, and we assume it is aligned with the ${\cal Z}$-axis to simplify the analysis. This chimney now begins to passively advect downstream, i.e., we assume the chimney does not induce creep closure nor affect the flow field, in a manner determined by the velocity profile Eqn (10), which here is set constant in time (Fig. S2a). The chimney is elongated over time (Fig. S2b), and its surface area is hence simultaneously enlarged Eqn. S3. However, as the height and the ${\cal X}{\cal Y}$-parallel cross-sectional area of the chimney remain constant, Cavalieri's principle dictates that the chimney volume stays constant as well. If the geothermal input stays fixed, then the constant volume of geothermal gases must dissipate the heat over the growing surface area. The specific melt rate will thus decrease, allowing creep closure to constrict the chimney. Further, as the chimney is advected down glacier, its lower effective cross-sectional area, perpendicular to the centreline of the chimney, decreases (Fig. S2c). These two mechanisms can close the chimney, or facilitate blockages of incoming debris, ice, or snow. Any sliding along the bed, as may be expected for temperate mountain glaciers, would serve to accelerate the advection of chimneys away from their sources.

2.2 Numerical modelling

To incorporate background flow fields into our models, a numerical model is needed as solving for the more complex stress regimes becomes impractical, if even possible, using analytical methods. Glaciovolcanic caves and chimneys are often observed within mountain glaciers that straddle the slopes of volcanoes (e.g., Curtis and Kyle, 2017), and have flow regimes that depend on the topography of the landscape. Jarosch and Gudmundsson (2007) showed that ignoring the flow of ice into a geothermally active region may significantly underestimate predictions of the basal heat flux. Modelling what effect the downslope flow of ice has on the formation and evolution of glaciovolcanic caves and chimneys will therefore give more realistic estimates of the relationships between volcanological and glaciological parameters.

2.2.1 Model setup and parameters

To further investigate the interactions between geothermally induced ice melt and glacier dynamics, we perform a set of forward numerical model simulations using synthetic glacier geometries. Our model employs two modules that run sequentially to capture the transient evolution of glaciovolcanic voids (Fig. S3).

1. 1. The ice-flow module employs Elmer/Ice, a three-dimensional (3-D) finite-element ice-flow model (Gagliardini and others, 2013), to solve the Stokes equations (e.g., Pattyn and others, 2008). This module takes in a mesh of the model domain, runs a diagnostic Elmer/Ice simulation to initialise variables, and then continuously deforms the model domain, i.e., the free surfaces associated with the void and glacier, by running a fixed number of prognostic Elmer/Ice simulations. A diagnostic simulation is required before each prognostic simulation to initialise the variables for the input mesh, as extrapolating variables after remeshing yielded poor results. The default setup runs four iterations using 6-hour time steps during each prognostic simulation, as this combination yielded the most stable simulations for all parameter ensembles during model testing. Note that the melt is not applied to the void within this module.

2. 2. The melt and remeshing module is written in Python and instantaneously applies the geothermally induced melt associated with the modelled time interval to the void, and then remeshes the model domain. This module takes in the deformed mesh and extracts the free surfaces of the void and glacier. Both surfaces are remeshed using the Python application programming interface of the open source 3-D finite-element mesh generator GMSH (Geuzaine and Remacle, 2009), and the melt is applied to the void surface using the Python package PyVista (Sullivan and Kaszynski, 2019). Finally, the void and glacier surfaces are used to reassemble the model domain and the volume is remeshed. Note that deformation due to ice dynamics is not applied within this module.

During model development, we attempted to simulate ice flow and geothermal melt entirely within Elmer/Ice, instead of treating these processes sequentially in the two modules described above, but this strategy proved unsuccessful for three main reasons that follow (Unnsteinsson, 2022). (i) Instabilities quickly arise and can distort the void surface into a irregular nonphysical geometry. This problem arises from the relatively large deformation and melt rates that can warp mesh elements to the point that they overlap. A solution would require both reduced time steps, on the order of minutes, and very fine, sub-metre, mesh resolutions in most cases, along with frequent remeshing to guarantee the integrity of the mesh. Such a solution would be immensely computationally inefficient, and cannot be applied to features of large spatial and temporal scales. (ii) Elmer/Ice does not currently have a solver that can create and remesh an unstructured 3-D mesh, thus requiring an external solution. (iii) The way surface normal vectors are computed in numerical models results in nodes that define the void–bed interface to be lifted off the bed when melt is applied. Hence the void–bed boundary is transferred to the distal end of an affected cell, resulting in additional synthetic void growth on the scale of the mesh resolution. The solution is either the same as for the first problem, which suffers from the same drawbacks, or to manipulate melt implementation at the boundary externally.

For the reasons above, we devised the two-step procedure, which allowed us to model glaciovolcanic voids for an array of realistic physical parameters, while maintaining computational feasibility. The model was run for a total of 1232 parameter ensembles, with four input parameters that were varied between simulations. The input parameters are: glacier thicknesses ${50}\, {\rm m } \leq H \leq {200}\, {\rm m }$, chosen to be representative of mountain glaciers where glaciovolcanic voids have been observed; bed/surface slopes ${0}^{\rm o} \leq \alpha \leq \, {15}^{\rm o}$; total geothermal heat fluxes ${0.5}\, {\rm MW } \leq \dot {Q} \leq {10}\, {\rm MW }$, inspired by measured values for fumaroles and hot springs (Gaudin and others, 2016); and four different heat transfer mechanisms (see below). Simulations run until they either crash, or are terminated when (i) the void width becomes 75 % of the domain width; (ii) the void height drops below 2 m; or (iii) time or iteration number limits are reached. The terminations due to void width and height limitations are imposed to prevent interference with boundary conditions and the disappearance of a free surface, respectively. Note that the model will crash when a void reaches the surface, as Elmer/Ice does not handle the merging of free surfaces. While we allow these simulations to run until they crash, we define melt-through to occur when the void height exceeds 90 % of the glacier thickness. With this step-wise approach, 48 hours of run time, in serial with 4000 MB allocated memory on a computer cluster, delivered up to 800 days of simulated time, depending on input parameters.

2.2.2 Model domain

With the aim of investigating a broad glaciological and volcanological parameter space, laterally limited domains of uniform ice thickness are used (Fig. 4). Compared to modelling fully realistic glacier geometries, this approach offers several advantages: (i) simplicity in changing glaciological parameters, i.e., glacier thickness and the bed/surface slope, hence velocity field; (ii) reduction of necessary computational resources due to considering only the area likely to be influenced by a glaciovolcanic void, and (iii) simplification of the interpretation of simulation results. The initiation of void formation cannot be modelled numerically, thus requiring an incipient void to be prescribed. The domain mesh is created using the Python application programming interface for GMSH, and all geometrical entities are automatically labelled so as to be compatible with Elmer/Ice (Unnsteinsson, 2022).

Figure 4. The numerical model domain (a) cross-sectional view and (b) bottom view. The glacier thickness H and bed/surface slope α are given, and the horizontal extent of the domain is set to be double the glacier thickness in both directions, with vertical boundaries on all sides. Note the mesh refinement around the proto-void (enlarged here for better visualisation).

2.2.3 Initial and boundary conditions

Both the glacier surface and the surface of the void are modelled as stress-free, i.e.,

(11)$$\sigma \cdot \mathbf{n}_{{\rm S}} = 0,\; $$

where $\mathbf{n}_{{\rm S}}$ is the normal vector to the ice surface. The void- and glacier surface are both subject to zero surface mass balance within the Elmer/Ice module, as the geothermally induced melt is implemented within the Python module. The four vertical sides of the model domain can be split into two groups that we will term passive boundaries (PBs) and active boundaries (ABs). Passive boundaries are those parallel to glacier flow on sloped model domains and all vertical boundaries on horizontal model domains. All passive boundaries have a prescribed no-flux condition:

(12)$$\mathbf{v} \cdot \mathbf{n}_{\rm PB} = 0,\; $$

where $\mathbf{n}_{{\rm PB}}$ is the normal vector to the passive boundary. While it would be possible to implement stress conditions here to simulate, for example, the drag of valley walls, in our case, this would only serve to alter glacier velocities which are already controlled by the slope. Active boundaries (ABs) are those that are perpendicular to glacier flow on sloped model domains, i.e., the upstream and downstream boundaries. For active boundaries we prescribe a vertical profile of slope-parallel velocity according to the shallow-ice approximation Eqn. (10). We chose the shallow-ice approximation boundary condition to have a standardised flux boundary condition across all model domains. Periodic boundary conditions could also have been used.

This study is restricted to fully water-drained glaciovolcanic voids, the observed existence of which suggests that efficient subglacial drainage is possible in these environments such that water pressure cannot be maintained. Observations of the crater glacier of Mount St. Helens (WA, USA), which hosts a glaciovolcanic cave network (e.g., Anderson and others, 1998), suggest that permeable bed geology inhibits basal water accumulation and thus the water pressures that would enhance sliding (Walder and others, 2005). It is therefore reasonable to assume that fast sliding is likely not a significant contributor to ice flow where glaciovolcanic voids form. Nonetheless, a temperate glacier flowing down the flank of a volcano would still be expected to demonstrate basal slip, at least in the form of regelation and enhanced creep (e.g, Weertman, 1957), with the possible addition of bed deformation (e.g., Iverson and others, 1998). Implementing sliding/friction laws into the model would, however, increase the number of model parameters that could be varied, increase the number of simulations needed and consequently complicate the analysis. Investigating the impact of sliding on the formation and evolution of glaciovolcanic voids is thus reserved for future work. We therefore impose a no-slip boundary condition at the bed, i.e., $\mathbf{v}( {\it z}_{\rm bed}) = \rm 0$.

For the Elmer/Ice module, zero velocity and pressure fields are used as initial conditions for the diagnostic and prognostic simulations. During these simulations, the Stokes equations are solved to satisfy the boundary conditions above. Thus, the void is only subjected to ice deformation during the prognostic Elmer/Ice simulations of each module iteration. The geothermally induced melt is only applied to the void after the prognostic step, via the Python module.

2.2.4 Melt implementation

The manner in which heat fluxes are distributed within the void is of vital importance, as the expansion/shrinkage of different sectors of a void is dictated by the heat transfer to those sectors. Gas flow within glaciovolcanic voids is complex (e.g., Florea and others, 2021), and accounting for the various mass- and heat-transfer mechanisms among the gas, glacier ice, and associated melt water would render such physically detailed modelling approaches impractical. Assumptions about the mass flux, vapour and chemical concentrations, and the initial temperature would have to be made to estimate the effect of processes such as vapour condensation, radiation and liquid film dynamics (Woodcock and others, 2015, 2016). Modelling gas flow within glaciovolcanic voids is computationally intensive in itself (Curtis, 2016), so accounting for these processes would demand more computational resources and further complicate analysis of results. We therefore rely on simple mathematical approximations of the heat transport within the voids instead of resolving the gas flow and associated processes. The total heat flux available for melting is specified and remains constant, while the specific heat fluxes at the void surface are computed based on the approximations below. The four different approximations of heat transport within voids, each meant to represent a physically plausible scenario, are:

1. M1 The case where low temperature geothermal sources are only capable of elevating the ambient temperature within a void, but do not drive vigorous vertical convection. The total heat flux, $\dot {Q}$, is thus distributed uniformly along the entire void surface, ${\cal S}$, with specific heat flux $\dot {q} = \dot {Q} / {\cal S}$.

2. M2 A geothermal gas plume driven by either natural or forced convection. The plume is assumed to be cylindrical, with radius r _plume, and the total heat flux is evenly distributed over the plume surface. Outside the plume, the specific heat flux decays radially as dictated by thermal radiation, $\dot {q} \propto \dot {Q} / r$. Should any portion of the void surface be within the plume, r < r _plume, the specific heat flux acting on that portion is prescribed to be the same as at the plume surface.

3. M3 A geothermal gas plume driven by either natural or forced convection, like M2, but the plume conforms to void geometry. This is done by setting the centre line of the plume to follow the centre line of the void, analogous to fluid flow in a tube. As for M2, the specific heat flux decays radially from the plume, and is kept constant within the plume.

4. M4 A geothermal point source with inefficient heat transport, with the specific heat flux away from the source set as thermal radiation from a point source: $\dot {q} \propto \dot {Q} / r^2$. As with M2 and M3, the specific heat flux is kept constant within a given radial distance.

The heat transport modes can be classified as “vertical modes” and “nonvertical modes”, depending on the presence or absence of vertical preference, respectively. M1 and M4 are nonvertical as their formulations essentially describe them as point sources, which distribute specific heat fluxes based on either surface area or total radial distance. For the vertical modes, M2 and M3, the heat sources are assumed to be (sub-)vertical and specific heat fluxes are distributed uniformly along their (sub-)vertical axes, but decay radially away from the axes.

For a given total heat flux and heat transport mode, the specific heat flux at each node of the void surface can be computed. The melt rate normal to the ice surface at each node is subsequently computed according to Eqn (6). Any node of the mesh, P_t, at time t, is then translated to node

(13)$$\eqalign{{\rm P}_{t + \Delta t}& = {\rm P}_{t} + {{\rm d}r_n\over {\rm d}t} \Delta t \hat{\mathbf{n}}_{\rm P} \cr & = {\rm P}_{t} + {\dot{q}( {\rm P}_{t}) \Delta t\over \rho_{{\rm i}} L} \hat{\mathbf{n}}_{\rm P},\; }$$

at time t + Δt, with $\dot {q}( {\rm P}_{t})$ the specific heat flux at node P_t and $\hat {\mathbf{n}}_{\rm P}$ the unit normal vector to the void surface at node P_t.

3. Results and discussion

3.1 Analytical models

3.1.1 Glaciovolcanic caves

By incorporating depth dependency into the creep closure formulation Eqn (7), we find that glaciovolcanic caves have a vertical preferential growth. The maximum physical value of the normal radius r _n in Eqn (7) can be found graphically for example, and is r _n = z = H/4, i.e., the maximum height of caves is one quarter of the ice thickness (Fig. 5). This limit appears to match that of observed void geometries during subglacial volcanic eruptions before brittle deformation of the overlying ice takes over (Tuffen, 2007). However, it is unlikely that our model accurately captures such a scenario to claim it explains the observations. For a given glacier thickness there must exist a corresponding maximum geothermal heat flux, $\dot {q}_{{\rm max}}$. If the maximum heat flux is surpassed, Eqn (7) has no physical solution for the normal radius r _n. Melt opening can therefore not be entirely counteracted by creep closure and voids will grow unstably, meaning that no steady state can exist under those conditions. The geothermal heat fluxes or glacier thicknesses that lead to these nonphysical solutions of Eqn (7) must therefore result in either complete melt-through, i.e., formation of a chimney, or the closure of the cave. Theoretical steady-state geometries for glaciovolcanic caves can therefore only exist for geothermal heat fluxes less than or equal to a critical value, and can never exceed a height of one quarter of the ice thickness. Equation (7) only yields physical results when the normal radius has a vertical component greater than or equal to zero. For the lower part of the cave geometry, we use a uniform radius equal to the radius where the normal angle becomes zero (Fig. 5a). The difference between this new cave geometry (solid lines in Fig. 5b) and Eqn (7) is always less than that of assuming a constant radius equal to the maximum height of the void, i.e., Nye's tunnel/sphere radius (dashed lines in Fig. 5b).

Figure 5. Analytically derived steady-state cave geometries for five fractions of the theoretical maximum void height: z _max = H/4, and the horizontal difference between the geometries and Eqn (7). (a) The new steady-state cave geometries given by a combination of Eqs. (7) and (4) (solid line) and radially symmetric geometries Eqn (4) (dashed line). (b) The horizontal difference between the geometries in (a) and Eqn (7). All axes are normalised to glacier thickness.

The departure of the geometry from a perfect cylinder or sphere violates the assumptions on which Eqn (2) is based. It is therefore clear that Eqn (7) does not perfectly describe glaciovolcanic cave geometries and the parameter space that allows for their existence. There are, however, two reasons that justify use of Eqn (7). (i) By computing the horizontal difference between Eqn (7) and our new derived cave geometries—a combination of Eqs. (7) and (4)—as well as the difference between Eqn (7) and the radially symmetric geometries Eqn (4), we can assess the quality of the new geometries (Fig. 5). In all cases our model performs better than the traditional radially symmetric geometries derived by Nye (1953), i.e., Eqn (4). However, the portions of the new geometries that are assumed radially symmetric naturally have a larger horizontal difference than those directly computed from Eqn (7). Our result that caves can not reach heights above one quarter of a glacier thickness may also support use of Eqn (2), as it may be seen as a weak form of r < <H. (ii) Our approach was deliberately designed to use, and expand on, simple existing analytical models Eqn (2) as has been commonly done in the literature for both symmetrical and asymmetrical voids (e.g., Walder, 1986; Hooke and others, 1990; Höskuldsson and Sparks, 1997; Tuffen and others, 2002; Tuffen, 2007).

3.1.2 Glaciovolcanic chimneys

We model glaciovolcanic chimneys similarly to glaciovolcanic caves, except the void axis is oriented perpendicular to the bed. Equation (9) permits the use of either known or estimated glaciovolcanic chimney geometries to determine the minimum subglacial geothermal heat fluxes needed to sustain chimneys in steady state. If unknown, the ice thickness could be assessed from global glacier thickness estimates such as those of Farinotti and others (2019) or Millan and others (2022) and the creep parameter from Table 3.4 in Cuffey and Paterson (2010), based on the glacier temperature. Estimating subglacial volcanic heat fluxes using Eqn (9) could prove useful for volcano monitoring efforts and hazard assessments.

As a chimney is advected downstream due to glacier flow, it becomes elongated. Thus, even if the radius of the chimney R does not change, its surface area grows over time with its length, according to Eqn. S3. As the chimney surface grows due to advection with glacier flow, the chimney requires a higher input of geothermal energy to maintain its geometry, based on Eqn (9). This need for increasing geothermal heat fluxes over time would mean that a chimney with a fixed heat flux budget will be advected into a state of thermal disequilibrium where the chimney must reduce its radius to reach a steady state (Fig. S2). Further, the effective cross-sectional area of chimneys, i.e., the cross-sectional area normal to a chimney centreline, shrinks as chimneys are deformed. These constrictions may lead to the complete closure of the chimney, or facilitate closure by other processes such as ice collapse.

The gas flow within the chimney is either dominated by forced or natural convection and can be described with the Boussinesq approximation (e.g., Bird and others, 2006). Since open nonwater filled voids within a glacier are considered here, the pressure within voids must approach atmospheric pressure. The pressure gradient between the geothermal source at the bed and chimney opening is therefore small and decreases even further the longer the chimney becomes. Under conditions where the pressure gradient becomes insignificant, buoyancy dominates the gas flow. Thus, if a chimney is sufficiently elongated, or becomes otherwise constricted, vertical growth of a new void is initiated. Such a process would lead to a cycle where: (i) a chimney forms above a geothermal source; (ii) the chimney is advected into an unstable thermal state; (iii) buoyancy-driven flow within the chimney becomes dominant and initiates the growth of a new chimney above the geothermal source resulting in the abandonment of the original chimney. There is some evidence for such a cycle having occurred within Job Glacier, where the advection of a chimney, and the formation of a new one directly upstream, have been observed (chimneys 2a and 2b in Fig. 1).

3.2 Numerical model

A total of 1232 numerical simulations were carried out to model the formation and evolution of glaciovolcanic voids. Four approximations of heat transport were used to melt voids into synthetic glacier geometries, over a range of glacier thicknesses, bed/surface slopes and total geothermal heat fluxes (Figs. 6, S4–S7). Intuition and analytical calculations based on the contraction of radially symmetric glacier voids, derived by Nye (1953) and expanded above, affirm that: (i) for a given glacier thickness, the void size should increase with increasing geothermal heat flux; and (ii) for a given geothermal heat flux, void sizes should decrease with increasing glacier thickness. Our numerical model results conform with this statement for all parameter ensembles, independent of heat transport modes and bed slope. Several key features can be seen in the timeseries of the model results: (i) Due to the basal melt parameterisation, void heights generally reach a steady state before either void areas or volumes (Figs. 6, S4–S7). (ii) The uniform specific heat flux (M1) generally produces a maximum void height, before flattening the void due to melting at the base and consequent roof lowering (Fig. S4). (iii) The radially controlled heat transport modes (M2, M3 and M4) reach and maintain a steady-state void height until the end of most simulations. (iv) The specific heat fluxes associated with heat transport modes M1, M2 and M3 tend towards a limit over time. The limit value increases with surface/bed slope in increments that grow with glacier thickness. (v) The void height in some simulations shows step-like behaviour. A step-like increase in height may be an indication of a positive feedback when a cauldron forms above a void, reducing the ice thickness which further slows the creep closure. A step-like decrease in height could be due to dynamical adjustment to widening at the base.

Figure 6. Example of the evolution of the model domain and velocity components x (1a–d) and z (2a–d), and stress components xx (3a–d) and xz (4a–d), for ice thickness $H = {100}{\rm m }$, bed/surface slope $\alpha = {10}^{\rm o}$, heat flux $\dot {Q} = {5.0}\, {\rm MW }$ and heat transport mode M2. The velocities are normalised to the maximum velocity given by the shallow-ice approximation Eqn (10), $v_{\rm SIA} = v_{{\cal X}}( {\cal Z} = H)$, and the stresses are normalised to the maximum ice overburden pressure p _i. The variables shown are a combination of the background fields and the fields induced by the presence of a void. Examples of other heat transport modes and additional physical variables can be seen in Figs. S4, S5, S6 and S7. Note that the figures are only meant to show equally spaced snapshots through the entire timeseries, not necessarily the evolution toward a steady state.

Akin to the results of Jarosch and Gudmundsson (2007), the effect of ice influx is apparent in our model results. Void heights, surface areas and volumes are limited by increasing surface slopes, here equal to bed slopes, for all heat transport modes. The degree of void suppression induced by bed slopes varies based on the glacier thickness and total geothermal heat flux. Thicker glaciers are generally more effective in limiting void sizes with increasing bed slope, due to higher gravitational driving stresses and therefore strain rates, compared to those of their thinner counterparts.

3.2.1 Steady-state voids

The height, area and volume of the void are useful to determine when voids reach a steady state. Note that here “void height” refers to vertical void height. To find steady states, the timeseries are first smoothed with a 11-day boxcar filter, and then differentiated. We define steady state as $\pm 1\%$ change in the time derivative of the variable in question (coloured points in Figs. 7 and S8).

Figure 7. Steady-state void heights of numerical model simulations, and the physical analytical approximation, Eqs. (14) and (S14), of the results using the theoretical values of the exponents. Numerical model simulation results are shown as filled circles r-rlying contours of the analytical approximations. X's indicate simulations that did not reach steady state, or otherwise failed. White circles are the nonsteady-state simulation results where void height exceeded 90 % of ice thickness. Results of all void-height approximations can be seen in Fig. S9.

To link the numerical model results to the glaciological and volcanological input variables, two analytical approximations were constructed (see Void-height approximation in Supp. Mat.). Both approximations give void height, h, as a function of glacier thickness, H, and geothermal heat flux, $\dot {Q}$, in the form of a general approximation

(14)$$h/H = a H^b \dot{Q}^c,\; $$

where the exponents have the same theoretical values, b = −n/2 and c = 1/2 Eqn. (S8). The difference between the two approximations lies in how the slope-dependence is incorporated. The first approximation equates the vertical melt rate to the creep closure and then uses a linear fit to the sine of the bed slope α. This yields Eqn. (S9), where coefficient a from the general approximation Eqn (14) becomes $a = a_{\rm L1} ( 1 + a_{\rm L2} \sin ( \alpha ) )$, and is here termed as the “linear” approximation. The second approximation computes the balance between the vertical melt rate, creep closure, and downward component of the ice velocity, and is termed the “physical” approximation. In this balance, the ice velocity introduces a dependence on bed slope by yielding $a = a_{\rm P1} ( 1 + a_{\rm P2} \sin ^{n + 1}( \alpha ) ) ^{-1/2}$. The latter approximation has universal theoretical values for coefficients a _P1 and a _P2, but for the former approximation, a _L1 only has a theoretical value for the case with zero bed slope. The constants a _L/_P1 and a _L/_P2 were fitted for both approximations (Supp. Mat.), using either the best fit for the exponents b and c or their theoretical values (Fig. S11 and Table S1).

The difference between the quality of the approximations is minimal, but the physical approximation slightly outperforms the linear one. The mean and median errors of the approximations, compared to simulation results, are <31 % and follow similar trends between heat transport modes and different bed slopes. For approximations using the best-fit values of constants b and c, the error increases with bed slope, often rapidly once a bed slope of 15° is reached. However, if the theoretical values of b and c are used, then the mean and median errors become more uniform over the range of tested bed slopes. These values are elevated compared to those of approximations using fitted b and c for lower bed slopes. The maximum mean and median errors are comparable between theoretical and best-fit values of b and c. We recommend using the physical approximation with the theoretical exponents b = −n/2 and c = 1/2, for the most consistent approximation results across the range of physical parameters tested here, yielding an uncertainty of ±30 % (Unnsteinsson, 2022).

3.2.2 Void aspect ratios

The aspect ratio of the resulting glaciovolcanic voids (Fig. 8) encodes important information about the manner of heat transport within the void. The height-to-mean-width ratio clearly distinguishes between vertical (green and yellow circles in Fig. 8) and nonvertical modes (blue and red circles in Fig. 8) of heat transport, with the former producing voids with aspect ratios >1 that are more likely to melt through the glacier. Nonvertical modes dissipate a significant portion of the heat budget to lateral expansion of the void, whereas vertical modes concentrate heat to melt upwards. Vertical modes are the only ones able to fully melt through a glacier for the parameters tested here. Based on our results, we hypothesise that vertical heat transport modes are required to form chimneys such as those observed within Job Glacier on the Mount Meager Volcanic Complex, British Columbia. There, hydrothermal plumes can be seen rising from the chimneys, indicating that buoyant convection is driving gas flow. The visible supraglacial gas plumes further demonstrate that the heat fluxes required to form a chimney exceed those required to sustain one. Using Eqn (14), which was derived for glaciovolcanic voids prior to melt-through, to estimate heat flux from an open chimney would lead to minimum estimates of the escaping heat flux. However, other methods would be required to better constrain the actual heat escape through open chimneys.

Figure 8. Aspect ratios of steady-state voids, using (a) void base width, or (b) an approximated average void width $\bar {w}_{\rm void} \approx \bar {w}_{\rm base} / 2$. The heights and widths are normalised to the ice thickness H. Symbol transparency scales with heat flux, with more opacity indicating higher heat fluxes. Marker size scales with the glacier thickness.

Steady-state aspect ratios (Fig. 8) show a distinct geometrical difference in the results between the vertical and nonvertical heat transport modes. Nonvertical heat transport modes (blue and red circles Fig. 8) yield voids with a lower height-to-base-width ratio than that of a hemispherical cave (1:2), but a height-to-mean-width ratio between one and a hemispherical cave (1:π/2). Vertical heat transport modes (green and yellow circles Fig. 8) give a height-to-base-width ratio between that of one and a hemispherical cave, but a height-to-mean-width ratio of nearly 2:1. The height-to-mean-width ratio has a tighter clustering of results, and is hence better suited for estimating heat transport modes based on void morphology.

Nonvertical modes may still lead to perforation of the glacier surface if glacier thinning is maintained as a result of elevated subglacial melt. Thinning may either facilitate a melt-through for a given heat flux Eqn (14), or eventually lead to the formation of a nunatak within the glacier, as in chimney 1a within Job Glacier (Supplementary Material). Either process could be further expedited by roof lowering, or collapse, the first stages of which may be seen in Fig. S4, although the physics required for complete roof lowering or collapse are not included in our model. Events associated with roof lowering and subsequent collapse are documented to have occurred on Mt. Baker (WA, USA) in 1975, when a sudden increase in geothermal activity formed new glaciovolcanic voids (Eichelberger and others, 1976). These roof collapses took place above voids likely formed from both vertical and nonvertical modes of heat transport, driven by fumaroles and geothermal waters, respectively.

The maximum specific heat flux reached in our model simulations is 32 kW m⁻², an order of magnitude lower than values associated with subglacial eruptions as determined by Woodcock and others (2015, 2016) in their analysis of heat transfer within vapour-dominated glaciovolcanic voids. They find that during subglacial eruptions most of a glaciovolcanic void surface is likely never subjected to specific heat fluxes higher than ~100 kW m⁻², but could reach as high as 25 MW m⁻². The lower maximum specific heat flux reported here are to be expected, as lower total heat fluxes are used in this study, and demonstrate that the approximations of the heat transfer within the voids produce realistic specific heat fluxes.

3.2.3 Model shortcomings

Several model simulations neither reach a steady state nor indicate chimney formation (black X's in Fig. 7). The reasons behind simulation failure can be split into three categories: (i) physical limitations of the model domain have been reached; (ii) failure during Elmer/Ice simulations, e.g., by not converging to a solution; and (iii) failure during remeshing. The code has inbuilt subroutines that rerun failed iterations with smaller timesteps, but some parameter combinations still yield meshes that are unusable.

Our numerical model has several relevant shortcomings (Unnsteinsson, 2022). First, the transient simulations must be done step-wise with an Elmer/Ice module that accounts for the ice deformation and a Python module that applies the melt and remeshes the model domain. A more efficient and satisfying solution would be to have the transient free-surface evolution of the void simulated entirely within a single module.

The model setup described here does not account for sliding at the glacier bed. The no-sliding condition is implemented by setting tangential ice velocities to zero as a boundary condition at the bed. This means that there is no closure rate at the bed that is able to counteract the geothermally induced melting, and hence void closures are only a result of roof lowering. This often causes the voids to expand laterally at the bed and may result in lower void heights due to the increase in surface area. It is therefore likely that heat fluxes are overestimated if computed from the void-height approximation.

The model is currently not capable of accounting for brittle roof collapse as voids approach the glacier surface. Future efforts could implement damage mechanics in the numerical model (e.g., Krug and others, 2014), informed by visco-elastic parameters from observational analysis (Ultee and others, 2020), to account for brittle roof collapse. Another possibility is to use an analytical model of the elastic deformation of the ice roof to determine radial stresses, similar to Ultee and others (2020), and assume a roof collapse once the crevasse depth, based on Nye (1955), at the boundary equals the roof thickness. Elmer/Ice does not accommodate the merging of two free surfaces, so some innovation would be required to deal with this limitation.

4. Implications

The glaciovolcanic voids at Job Glacier have shown remarkable variability and dynamic activity since their discovery in 2016. They have opened and closed, grown and shrunk throughout the period, with some remaining fixed in space while others were advected with the glacier ice. The most interesting example is that of chimneys 2a and 2b, which likely underwent a cyclical transformation (Fig. 1). Our hypothesis is that chimney 2a was advected into an unstable thermal state where it could not sustain itself, initiating the formation of chimney 2b (Fig. S2). This cycle resembles that described by Catania and Neumann (2010) for moulins in the Greenland Ice Sheet, whereby moulins are progressively advected away from their stationary sources (supraglacial lakes in this case) and then reform periodically in their original locations. While the source of the void is surface water instead of basal heat, as in the case of Job Glacier, the cyclic void formation due to void advection relative to a stationary source is comparable. Such a cycle could continue indefinitely if the supply of ice were not diminished, as is expected at Job Glacier over the coming decades. Alternatively, chimneys 2b and 4a could act like chimney 1a and create new nunataks within the glacier. While the fate of these chimneys remains to be seen, estimates of the heat fluxes associated with the glaciovolcanic voids can be made using our analytical and numerical model results.

Given the measured physical parameters, our analytical formulations can be used to estimate the geothermal heat fluxes associated with the glaciovolcanic chimneys within Job Glacier. Equation (9) gives the total subglacial geothermal heat flux for a chimney with a known radial geometry r(z). If the radius is uniform with depth, r(z) = R, the heat flux becomes

(15)$$\dot{Q} = 2 \pi \rho_i L A \left({\rho_i g\over n}\right)^n {H^{n + 1}\over n + 1} R^2.$$

The radii of the chimneys of Job Glacier are on the order of 5 m to 10 m and the ice thickness near them was 60 m to 90 m in 2018 (Warwick, 2020). If the radii are assumed uniform with depth, the heat fluxes associated with the chimneys range from 0.01 MW to 0.2 MW. These heat flux estimates come with significant caveats: background flow fields and the escape of gases are not accounted for, while the real void geometries are more complex. Heat flux estimates made from Eqn (15) are expected to be relatively low, as Eqn (15) only includes heat fluxes required to maintain a perfectly vertical chimney in steady state. Observations indicate that the heat sources beneath Job Glacier are fumaroles (Supp. Mat.), but Eqn (15) yields heat fluxes more akin to hot springs or degassing (Gaudin and others, 2016) which seems unlikely. We therefore consider these heat flux estimates as conservative lower bounds, with actual heat fluxes likely being significantly higher.

The void-height approximations based on the numerical model results according to Eqs. (14), (S9) and (S14) can also provide important insights into the mechanisms behind the formation of glaciovolcanic voids within Job Glacier. Since these glaciovolcanic voids have melted through the glacier, we can assert that heat fluxes are driven upwards due to forced or natural convection. That is, coefficients for heat transport mode M2 are the most appropriate to analyse the heat fluxes associated with the voids (Table S1). Using glacier thicknesses of 60 m to 90 m and a bed and surface slope of 15° around the chimneys yields total heat fluxes of 3 MW to 8 MW. These heat flux estimates are subject to uncertainties due to model limitations, such as the lack of void closure at the bed, the possible escape of heat, and disregard of roof collapse.

A discrepancy between the analytical and numerical model results was to be expected, but we regard the latter to have more validity based on its less extensive shortcomings. As a conservative estimate, we expect the heat fluxes associated with the smaller glaciovolcanic voids within Job Glacier to be $\lesssim$2 MW, and those associated with the larger voids to be 5 ± 3 MW. Heat fluxes of this magnitude would yield a total heat flux for the subglacial fumarole field beneath Job Glacier well in excess of 10 MW. A total heat flux of this magnitude, if it were to be uniformly distributed over the 3.2 km² (RGI Consortium, 2017) area of Job Glacier, would yield a uniform basal melt rate of ~30 cm a⁻¹. By comparison, the regional specific geothermal heat flux, <100 mW m⁻² (Grasby and others, 2009) would yield an estimated basal melt rate of $\lesssim$1 cm a⁻¹.

Changes in glacier thickness are more likely to drive the melt-through of voids, rather than increases in heat flux. The effect of glacier thickness changes are proportional to the exponent b, and the effect of changes in geothermal heat flux are proportional to the exponent c (Supp. Mat.). The theoretical values for exponents b and c are −1.5 and 0.5, respectively. Thus any change in glacier thickness will affect the normalised void height to greater extent than changes in heat flux. This aligns with observations at Job Glacier. A glaciovolcanic chimney may have been present in 1947 (Roberti and others, 2021) after a period of glacier thinning, but is not visible a few years later when glacier mass balance in the region is in equilibrium or positive (Moore and others, 2009). The chimneys appear again in 2015/2016 after a sustained period of glacier thinning (Roberti and others, 2018), while multiple new voids become visible after substantial melt during a heat wave in 2021 (WMO, 2022; White and others, 2023). While the correlation between glacier thinning and glaciovolcanic void emergence as observed at the surface of Job Glacier has a mechanistic explanation, we cannot unequivocally determine causation without more comprehensive data to independently constrain the heat fluxes. The continued thinning of glaciers worldwide could increase the probability of glaciovolcanic voids fully melting through any overlying ice. The formation of glaciovolcanic voids has been observed as a precursory signal to volcanic eruptions (Bleick and others, 2013; Reynolds and others, 2018). Their appearance could therefore be misattributed to volcanic activity, rather than glacier thinning, hence potentially raising false alarms of volcanic unrest.

5. Conclusion

In this study we used a combination of analytical and numerical models to investigate the formation, evolution and steady-state existence of water-drained voids within glaciers that arise from glaciovolcanic interactions. We derived steady-state glaciovolcanic cave geometries, and derived an expression from which the total heat flux associated with a known glaciovolcanic chimney geometry can be inferred. The analytical model suggests glaciovolcanic cave height is limited to one quarter of the ice thickness before chimney formation is initiated.

We use a numerical model to ascertain the relative importance of glaciological and volcanological parameters by running simulations over a range of ice thicknesses, bed/surface slopes, total geothermal heat fluxes and heat transport modes. Results confirm that the presence of background flow fields within glaciers suppresses the formation of glaciovolcanic voids. By combining our analytical models and numerical model results we devised an analytical approximation that links glaciovolcanic void height, h, to glacier thickness, H, and total geothermal heat flux, $\dot {Q}$, through a power law: $h/H = a H^b \dot {Q}^c$.

Using the analytical models, and the void-height approximation of the numerical model results, we evaluated the power of the subglacial geothermal heat sources beneath Job Glacier. Conservative estimates give a heat flux of up to 8 MW for each fumarole, and a total power well in excess of 10 MW for the subglacial geothermal area. This total heat flux results in modest mass loss for Job Glacier, but is capable of maintaining the current glaciovolcanic features.

The void-height approximation demonstrates that a reduction in glacier thickness is three times as effective as an increase in heat flux in driving glaciovolcanic void expansion. We thus expect that as glaciers worldwide continue to lose mass, more glaciovolcanic voids will form and melt through glaciers. Caution must be exercised in attributing the appearance of new glaciovolcanic voids to volcanic unrest, as they may instead be manifestations of global warming.