# Ascent of volatile-rich felsic magma in dykes: a numerical model applied to deep-sourced porphyry intrusions

R. Abdullin, O. Melnik, A. Rust, J. Blundy, E. Lgotina, S. Golovin

► To cite this version:

R. Abdullin, O. Melnik, A. Rust, J. Blundy, E. Lgotina, et al. Ascent of volatile-rich felsic magma in dykes: a numerical model applied to deep-sourced porphyry intrusions. Geophysical Journal International, 2024, 236, pp.1863-1876. 10.1093/gji/ggae027. insu-04604401

HAL Id: insu-04604401

https://insu.hal.science/insu-04604401

Submitted on 7 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Ascent of volatile-rich felsic magma in dykes: a numerical model applied to deep-sourced porphyry intrusions

R. Abdullin, O. Melnik, A. Rust, J. Blundy, E. Lgotina, S. Golovin

► To cite this version:

R. Abdullin, O. Melnik, A. Rust, J. Blundy, E. Lgotina, et al. Ascent of volatile-rich felsic magma in dykes: a numerical model applied to deep-sourced porphyry intrusions. Geophysical Journal International, 2024, 236, pp.1863-1876. 10.1093/gji/ggae027. insu-04604401

HAL Id: insu-04604401

https://insu.hal.science/insu-04604401v1

Submitted on 7 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ascent of volatile-rich felsic magma in dykes: a numerical model applied to deep-sourced porphyry intrusions

R. Abdullin<sup>1,2</sup> O. Melnik<sup>3,4</sup> A. Rust<sup>5</sup> J. Blundy<sup>6</sup> E. Lgotina<sup>1</sup> and S. Golovin<sup>1,2</sup>

Accepted 2024 January 17. Received 2023 December 22; in original form 2023 February 24

## SUMMARY

Dyke propagation is a mechanism for more rapid ascent of felsic magmas through the crust than is possible via diapirs or percolative flow. As it ascends, the magma undergoes complex physical and chemical transformations induced by decompression and cooling. These processes dramatically change the magma density and viscosity, which in turn affect magma ascent rate and the depth at which the dyke arrests. We present a mathematical model of dyke propagation for silicic magmas taking into account the presence of multiple volatile species ($H_2O$ and $CO_2$), bubble growth, heat advection and loss, crystallization and latent heat release. We consider conditions for dykes associated with porphyry ore deposits, which may represent an end-member in rapid ascent of felsic magmas from depth. In particular, we simulate the propagation of dykes launched from a deep (900 MPa), volatile-saturated magma source, testing the effects of the magma $H_2O/CO_2$ content, temperature and mass on its ascent rate and final emplacement depth. The model predicts short ascent times (hours to days), with a large increase in viscosity at shallow depth, leading to stagnation and solidification of the dyke. Higher initial water content, higher temperature and larger mass of the magma in the dyke promote faster propagation and shallower arrest. Volatile loss from ascending magma remains limited until the stagnation depth, providing a potential mechanism for transfer of deep volatiles to hypabyssal blind intrusions associated with porphyry ore deposits. Our findings are applicable to the problem of silicic magma ascent through the crust more generally.

**Key words:** Dyke propagation model; Volatile exsolution; Magma viscosity; Buoyancy; Porphyry deposits.

## 1 INTRODUCTION

Buoyant felsic magmas can ascend through the crust via dykes, diapirs and percolative flow. Although all three mechanisms can play critical roles in the ascent of magma from mid- or lower-crustal sources to the shallow emplacement depths (km) of most felsic plutons, subvolcanic reservoirs and hypabyssal intrusions associated with porphyry ore deposits, only dykes can enable magma transfer to be very rapid (e.g. <10 yr; e.g. Brandon et al. 1996; Petford et al. 2000). Our understanding of the ascent rates of felsic dykes through the crust is based largely on theoretical models, with an emphasis on understanding the conditions (e.g. dyke width, magma properties) for which dykes can survive thermally and so transit through the crust rather than solidify close to source (e.g. Petford et al. 1993, 1994; Rubin 1995). There are some constraints from the geological record, including ascent rates from deep sources (>600 MPa) to the

final level of emplacement (e.g. 200 MPa) based on the kinetics of epidote resorption in felsic intrusions (e.g. up to 14 km yr⁻<fnref n="1" /> for a case study of rhyodacite dykes; Brandon et al. 1996). However, most studies have focused on shallower regimes, and the best constraints on ascent rates of felsic magmas through the upper crust are from magmas that reached the surface in eruptions (Rutherford 2008; Browne & Szramek 2015). In historical eruptions, the ascent rate can be estimated from the duration and volume of erupted magma; other constraints come from the kinetics of crystallization and vesiculation, and timescales of diffusion through melts and crystals. For example, crystallization kinetics indicate that rhyolite magma of the 2008 eruption of Chaiten volcano ascended at about 1 m s⁻<fnref n="1" /> from >5 km depth to near surface in ~4 hr (Castro & Dingwell 2009). Such high ascent rates in the shallow crust will be aided by vesiculation of volatile-saturated magma. Volatiles are expected to have less effect on felsic magma ascent in the lower- and mid-crust because

<footnote n="1"><sup>1</sup>Novosibirsk State University, Novosibirsk-90, 2 Pirogova Str. 630090, Russia</footnote>

<footnote n="2"><sup>2</sup>Lavrentyev Institute of Hydrodynamics SB RAS, 15, Ac. Lavrentieva ave., Novosibirsk, 630090, Russia</footnote>

<footnote n="3"><sup>3</sup>Institut des Sciences de la Terre, University Grenobles Alpes, CS 40700 38058 GRENOBLE Cedex 9, France E-mail: oemelnik@gmail.com</footnote>

<footnote n="4"><sup>4</sup>Institute of Mechanics, Moscow State University, 1 Michurinskii prospekt, 119192, Moscow, Russia</footnote>

<footnote n="5"><sup>5</sup>School of Earth Sciences, University of Bristol, Queen's Road, Bristol, BS8 1RJ Bristol, UK</footnote>

<footnote n="6"><sup>6</sup>Department of Earth Sciences, University of Oxford, South Park Road, Oxford, OX1 3AN Oxford, UK</footnote>
# Ascent of volatile-rich felsic magma in dykes: a numerical model applied to deep-sourced porphyry intrusions

R. Abdullin<sup>1,2</sup> O. Melnik<sup>3,4</sup> A. Rust<sup>5</sup> J. Blundy<sup>6</sup> E. Lgotina<sup>1</sup> and S. Golovin<sup>1,2</sup>

Accepted 2024 January 17. Received 2023 December 22; in original form 2023 February 24

## SUMMARY

Dyke propagation is a mechanism for more rapid ascent of felsic magmas through the crust than is possible via diapirs or percolative flow. As it ascends, the magma undergoes complex physical and chemical transformations induced by decompression and cooling. These processes dramatically change the magma density and viscosity, which in turn affect magma ascent rate and the depth at which the dyke arrests. We present a mathematical model of dyke propagation for silicic magmas taking into account the presence of multiple volatile species ($H_2O$ and $CO_2$), bubble growth, heat advection and loss, crystallization and latent heat release. We consider conditions for dykes associated with porphyry ore deposits, which may represent an end-member in rapid ascent of felsic magmas from depth. In particular, we simulate the propagation of dykes launched from a deep (900 MPa), volatile-saturated magma source, testing the effects of the magma $H_2O/CO_2$ content, temperature and mass on its ascent rate and final emplacement depth. The model predicts short ascent times (hours to days), with a large increase in viscosity at shallow depth, leading to stagnation and solidification of the dyke. Higher initial water content, higher temperature and larger mass of the magma in the dyke promote faster propagation and shallower arrest. Volatile loss from ascending magma remains limited until the stagnation depth, providing a potential mechanism for transfer of deep volatiles to hypabyssal blind intrusions associated with porphyry ore deposits. Our findings are applicable to the problem of silicic magma ascent through the crust more generally.

**Key words:** Dyke propagation model; Volatile exsolution; Magma viscosity; Buoyancy; Porphyry deposits.

## 1 INTRODUCTION

Buoyant felsic magmas can ascend through the crust via dykes, diapirs and percolative flow. Although all three mechanisms can play critical roles in the ascent of magma from mid- or lower-crustal sources to the shallow emplacement depths (km) of most felsic plutons, subvolcanic reservoirs and hypabyssal intrusions associated with porphyry ore deposits, only dykes can enable magma transfer to be very rapid (e.g. <10 yr; e.g. Brandon et al. 1996; Petford et al. 2000). Our understanding of the ascent rates of felsic dykes through the crust is based largely on theoretical models, with an emphasis on understanding the conditions (e.g. dyke width, magma properties) for which dykes can survive thermally and so transit through the crust rather than solidify close to source (e.g. Petford et al. 1993, 1994; Rubin 1995). There are some constraints from the geological record, including ascent rates from deep sources (>600 MPa) to the

final level of emplacement (e.g. 200 MPa) based on the kinetics of epidote resorption in felsic intrusions (e.g. up to 14 km yr⁻<fnref n="1" /> for a case study of rhyodacite dykes; Brandon et al. 1996). However, most studies have focused on shallower regimes, and the best constraints on ascent rates of felsic magmas through the upper crust are from magmas that reached the surface in eruptions (Rutherford 2008; Browne & Szramek 2015). In historical eruptions, the ascent rate can be estimated from the duration and volume of erupted magma; other constraints come from the kinetics of crystallization and vesiculation, and timescales of diffusion through melts and crystals. For example, crystallization kinetics indicate that rhyolite magma of the 2008 eruption of Chaiten volcano ascended at about 1 m s⁻<fnref n="1" /> from >5 km depth to near surface in ~4 hr (Castro & Dingwell 2009). Such high ascent rates in the shallow crust will be aided by vesiculation of volatile-saturated magma. Volatiles are expected to have less effect on felsic magma ascent in the lower- and mid-crust because

<footnote n="1"><sup>1</sup>Novosibirsk State University, Novosibirsk-90, 2 Pirogova Str. 630090, Russia</footnote>

<footnote n="2"><sup>2</sup>Lavrentyev Institute of Hydrodynamics SB RAS, 15, Ac. Lavrentieva ave., Novosibirsk, 630090, Russia</footnote>

<footnote n="3"><sup>3</sup>Institut des Sciences de la Terre, University Grenobles Alpes, CS 40700 38058 GRENOBLE Cedex 9, France E-mail: oemelnik@gmail.com</footnote>

<footnote n="4"><sup>4</sup>Institute of Mechanics, Moscow State University, 1 Michurinskii prospekt, 119192, Moscow, Russia</footnote>

<footnote n="5"><sup>5</sup>School of Earth Sciences, University of Bristol, Queen's Road, Bristol, BS8 1RJ Bristol, UK</footnote>

<footnote n="6"><sup>6</sup>Department of Earth Sciences, University of Oxford, South Park Road, Oxford, OX1 3AN Oxford, UK</footnote>
volatile solubilities are greater at higher pressure, and even if the volatile content is sufficient for saturation, the resulting bubbles of supercritical fluid will be much denser and less compressible than magmatic vapours exsolved at lower pressures in the shallow crust. However, high dissolved water content can also exercise important controls on the ascent of deep felsic dykes because of the reduced viscosity, increased buoyancy and lower liquidus temperature (and so deeper onset of decompression-induced crystallization) of melts with high dissolved water contents.

Here we develop a numerical model for felsic dyke ascent, which we use to explore the role of high volatile ($H_2O \pm CO_2$) contents of a lower crustal, deep (900 MPa) felsic magma source on the ascent rate and final emplacement depths of dykes. Although the model has broad applications to magma ascent through the crust, the conditions we test here are motivated by the generation of porphyry Cu ($\pm Mo \pm Au$) deposits, which may represent an end-member in fast ascent of volatile-rich, deep-sourced felsic magmas. The ore metals ultimately derive from fluids exsolved from magma, yet mass balance requires the metals (and carrier fluids) be sourced from a substantially greater volume of magma than the hypabyssal intrusions that are spatially and temporally associated with mineralization (Sillitoe 2010). Porphyry intrusions typically comprise intermediate to felsic composition arc magmas with adakite-like geochemical signatures (e.g., high Sr/Y) attributed to high water contents at depth that suppressed plagioclase saturation, but encouraged amphibole ($\pm$ garnet) crystallization (Richards 2011; Loucks 2014; Rezeau & Jagoutz 2020). Near-Moho magma source depths and >10 wt per cent dissolved $H_2O$ are considered plausible from a geochemistry perspective (Chiaradia & Caricchi 2017; Loucks 2021). Inclusions of melt that were trapped in crystals as they grew can provide more direct constraints on the volatile contents of the parental melts of porphyry deposits, although such data are relatively scarce in the literature. A recent study of $H_2O$ and $CO_2$ dissolved in re-homogenized zircon-hosted melt inclusions obtained volatile saturation pressures corresponding to depths >20 km in all three porphyry deposits sampled, with a maximum saturation depth of ~30 km for Grasberg (Indonesia), with ~5.5 wt per cent $H_2O$ and a remarkably high $CO_2$ content of ~4000 ppm (Butters 2022). Further evidence for high $H_2O$ and $CO_2$ contents of magmas associated with porphyry deposits comes from quartz-hosted glassy melt inclusions from the fallout of an explosive volcanic eruption of similar age, composition and mineralogy to the Pine Grove (USA) Mo-porphyry system (Lowenstein 1994). The inclusions show a range of saturation pressures consistent with magma degassing as it ascended from ~16 km depth; inclusions with saturation pressures >400 MPa have >7 wt per cent $H_2O$ ($\le 7.99$ per cent) and >800 ppm $CO_2$ ($\le 972$ ppm). Although melt inclusion data relevant to porphyry deposits are limited, they do provide strong evidence that the sources are deep and wet, as suggested by geochemical data (e.g. high rock Sr/Y) and that porphyry-forming magmas have high $CO_2$ contents. It has been argued on petrological grounds (Blundy et al. 2010) that arc magmas generally are more $CO_2$-rich than has been previously recognized.

The significance of the high volatile contents of magmas associated with porphyry ore deposits has largely focused on chemical consequences; here, we consider the physical implications, along with the effect of the temperature and volume of magma ascending in dykes. In particular, we hypothesize that high volatile contents can enable magma ascent that is fast relative to both magma solidification and magmatic fluid loss by permeable flow into host rocks. In this way, a volatile-rich, deep-sourced magma can reach the shallow crust with most of its fluid (and metal) cargo intact. However,

to form an ore deposit, the fast ascent cannot be a runaway process. If magma accelerates to the surface and erupts explosively, it will release most of its fluids (and metals) to the atmosphere, as occurs in $CO_2$-charged ultramafic kimberlite eruptions (Sparks 2013), and felsic Plinian eruptions (e.g. Pinatubo in 1991, Chaiten in 2008). On this basis, it has been proposed that porphyry-forming systems represent failed giant ignimbrite eruptions (Chiaradia & Caricchi 2022). Instead of reaching the surface, the rapidly ascending magma must decelerate so that it is emplaced in the shallow crust to form a hypabyssal intrusion; porphyry ore deposit-related dykes often appear to be ‘blind’ meaning that the magma did not reach the surface and erupt (Sillitoe 2010). The deceleration of magma to form an intrusion is most readily explained by magma rheology changes due to crystallization induced by cooling and/or degassing during decompression, a process termed ‘viscous death’ by Annen et al. (2006). At the point of magma deceleration and stalling, entrained fluids exsolved during decompression (or first boiling) will be released, along with additional dissolved fluids produced due to crystallization (or second boiling). Explosive fluid release features, such as veining and brecciation, are commonly associated with blind porphyry dykes in mineralized systems (Sillitoe 2010).

Exsolution of volatiles at depth will tend to aid rapid dyke ascent through increased buoyancy; however, exsolving $H_2O$ can also impede magma ascent because dehydration increases melt viscosity and drives crystallization through its effect on liquidus temperature. The more $CO_2$ is present, the deeper the magma becomes volatile-saturated, and the greater the amount of exsolved fluid for a given pressure and bulk $H_2O$ content. Although dissolved $CO_2$ has little effect on melt properties and phase relations, it does reduce the activity of $H_2O$ in the melt, which in turn affects the liquidus temperature. This arises because of the differing solubilities of $H_2O$ and $CO_2$. At the onset of volatile saturation of a melt containing dissolved $H_2O$ and $CO_2$, preferential loss of $CO_2$ to the fluid raises $H_2O$ activity in the residual melt, lowering the liquidus temperature. This delays the onset of degassing-induced crystallization, which is a key factor in slowing magma ascent, a phenomenon called ‘viscous death’ (Annen et al. 2006). Instead, degassing takes place over a much greater depth range unaccompanied by crystallization. This effect, called ‘liquidus retreat’ by Blundy et al. (2010), has important consequences for the physics of magma ascent. Note that a similar process occurs during decompression of $CO_2$-free, $H_2O$-undersaturated magmas but no volatile phase is exsolved prior to $H_2O$-saturation. As the presence of bubbles impacts the density (and buoyancy) of the magma, there are important differences between the ascent of volatile-saturated and volatile-undersaturated magmas having the same initial $H_2O$ content.

To quantify the effects of $H_2O$ and $CO_2$, temperature and volume on magma ascent rate and emplacement depth, we present a model for dyke ascent in an elastic medium that includes volatile exsolution, crystallization and heat exchange with host rocks. The magma is modelled as a Newtonian fluid with a viscosity that evolves with changes in dissolved $H_2O$ content, temperature and crystal fraction. Crystallization follows equilibrium; if crystallization kinetics were included, it would make the ascent even more rapid (due to lower crystallinity and hence viscosity) and cause more abrupt deceleration of the magma due to a delayed burst of non-equilibrium crystallization. We model the end-member scenario where all the exsolved volatiles travel with the magma (i.e. closed-system degassing), which maximizes buoyancy. The limitations and implications of this and other simplifications are considered in the discussion.
volatile solubilities are greater at higher pressure, and even if the volatile content is sufficient for saturation, the resulting bubbles of supercritical fluid will be much denser and less compressible than magmatic vapours exsolved at lower pressures in the shallow crust. However, high dissolved water content can also exercise important controls on the ascent of deep felsic dykes because of the reduced viscosity, increased buoyancy and lower liquidus temperature (and so deeper onset of decompression-induced crystallization) of melts with high dissolved water contents.

Here we develop a numerical model for felsic dyke ascent, which we use to explore the role of high volatile ($H_2O \pm CO_2$) contents of a lower crustal, deep (900 MPa) felsic magma source on the ascent rate and final emplacement depths of dykes. Although the model has broad applications to magma ascent through the crust, the conditions we test here are motivated by the generation of porphyry Cu ($\pm Mo \pm Au$) deposits, which may represent an end-member in fast ascent of volatile-rich, deep-sourced felsic magmas. The ore metals ultimately derive from fluids exsolved from magma, yet mass balance requires the metals (and carrier fluids) be sourced from a substantially greater volume of magma than the hypabyssal intrusions that are spatially and temporally associated with mineralization (Sillitoe 2010). Porphyry intrusions typically comprise intermediate to felsic composition arc magmas with adakite-like geochemical signatures (e.g., high Sr/Y) attributed to high water contents at depth that suppressed plagioclase saturation, but encouraged amphibole ($\pm$ garnet) crystallization (Richards 2011; Loucks 2014; Rezeau & Jagoutz 2020). Near-Moho magma source depths and >10 wt per cent dissolved $H_2O$ are considered plausible from a geochemistry perspective (Chiaradia & Caricchi 2017; Loucks 2021). Inclusions of melt that were trapped in crystals as they grew can provide more direct constraints on the volatile contents of the parental melts of porphyry deposits, although such data are relatively scarce in the literature. A recent study of $H_2O$ and $CO_2$ dissolved in re-homogenized zircon-hosted melt inclusions obtained volatile saturation pressures corresponding to depths >20 km in all three porphyry deposits sampled, with a maximum saturation depth of ~30 km for Grasberg (Indonesia), with ~5.5 wt per cent $H_2O$ and a remarkably high $CO_2$ content of ~4000 ppm (Butters 2022). Further evidence for high $H_2O$ and $CO_2$ contents of magmas associated with porphyry deposits comes from quartz-hosted glassy melt inclusions from the fallout of an explosive volcanic eruption of similar age, composition and mineralogy to the Pine Grove (USA) Mo-porphyry system (Lowenstein 1994). The inclusions show a range of saturation pressures consistent with magma degassing as it ascended from ~16 km depth; inclusions with saturation pressures >400 MPa have >7 wt per cent $H_2O$ ($\le 7.99$ per cent) and >800 ppm $CO_2$ ($\le 972$ ppm). Although melt inclusion data relevant to porphyry deposits are limited, they do provide strong evidence that the sources are deep and wet, as suggested by geochemical data (e.g. high rock Sr/Y) and that porphyry-forming magmas have high $CO_2$ contents. It has been argued on petrological grounds (Blundy et al. 2010) that arc magmas generally are more $CO_2$-rich than has been previously recognized.

The significance of the high volatile contents of magmas associated with porphyry ore deposits has largely focused on chemical consequences; here, we consider the physical implications, along with the effect of the temperature and volume of magma ascending in dykes. In particular, we hypothesize that high volatile contents can enable magma ascent that is fast relative to both magma solidification and magmatic fluid loss by permeable flow into host rocks. In this way, a volatile-rich, deep-sourced magma can reach the shallow crust with most of its fluid (and metal) cargo intact. However,

to form an ore deposit, the fast ascent cannot be a runaway process. If magma accelerates to the surface and erupts explosively, it will release most of its fluids (and metals) to the atmosphere, as occurs in $CO_2$-charged ultramafic kimberlite eruptions (Sparks 2013), and felsic Plinian eruptions (e.g. Pinatubo in 1991, Chaiten in 2008). On this basis, it has been proposed that porphyry-forming systems represent failed giant ignimbrite eruptions (Chiaradia & Caricchi 2022). Instead of reaching the surface, the rapidly ascending magma must decelerate so that it is emplaced in the shallow crust to form a hypabyssal intrusion; porphyry ore deposit-related dykes often appear to be ‘blind’ meaning that the magma did not reach the surface and erupt (Sillitoe 2010). The deceleration of magma to form an intrusion is most readily explained by magma rheology changes due to crystallization induced by cooling and/or degassing during decompression, a process termed ‘viscous death’ by Annen et al. (2006). At the point of magma deceleration and stalling, entrained fluids exsolved during decompression (or first boiling) will be released, along with additional dissolved fluids produced due to crystallization (or second boiling). Explosive fluid release features, such as veining and brecciation, are commonly associated with blind porphyry dykes in mineralized systems (Sillitoe 2010).

Exsolution of volatiles at depth will tend to aid rapid dyke ascent through increased buoyancy; however, exsolving $H_2O$ can also impede magma ascent because dehydration increases melt viscosity and drives crystallization through its effect on liquidus temperature. The more $CO_2$ is present, the deeper the magma becomes volatile-saturated, and the greater the amount of exsolved fluid for a given pressure and bulk $H_2O$ content. Although dissolved $CO_2$ has little effect on melt properties and phase relations, it does reduce the activity of $H_2O$ in the melt, which in turn affects the liquidus temperature. This arises because of the differing solubilities of $H_2O$ and $CO_2$. At the onset of volatile saturation of a melt containing dissolved $H_2O$ and $CO_2$, preferential loss of $CO_2$ to the fluid raises $H_2O$ activity in the residual melt, lowering the liquidus temperature. This delays the onset of degassing-induced crystallization, which is a key factor in slowing magma ascent, a phenomenon called ‘viscous death’ (Annen et al. 2006). Instead, degassing takes place over a much greater depth range unaccompanied by crystallization. This effect, called ‘liquidus retreat’ by Blundy et al. (2010), has important consequences for the physics of magma ascent. Note that a similar process occurs during decompression of $CO_2$-free, $H_2O$-undersaturated magmas but no volatile phase is exsolved prior to $H_2O$-saturation. As the presence of bubbles impacts the density (and buoyancy) of the magma, there are important differences between the ascent of volatile-saturated and volatile-undersaturated magmas having the same initial $H_2O$ content.

To quantify the effects of $H_2O$ and $CO_2$, temperature and volume on magma ascent rate and emplacement depth, we present a model for dyke ascent in an elastic medium that includes volatile exsolution, crystallization and heat exchange with host rocks. The magma is modelled as a Newtonian fluid with a viscosity that evolves with changes in dissolved $H_2O$ content, temperature and crystal fraction. Crystallization follows equilibrium; if crystallization kinetics were included, it would make the ascent even more rapid (due to lower crystallinity and hence viscosity) and cause more abrupt deceleration of the magma due to a delayed burst of non-equilibrium crystallization. We model the end-member scenario where all the exsolved volatiles travel with the magma (i.e. closed-system degassing), which maximizes buoyancy. The limitations and implications of this and other simplifications are considered in the discussion.
## 1.1 Dyke propagation models

We start by reviewing dyke propagation models to put the new modelling into context. Dyke propagation modelling involves the coupling of several processes:

(i) Mechanical deformation of host rocks induced by magma pressure on the fracture surfaces.

(ii) Transport of magma within the fracture, taking into account the evolution of physical and rheological properties of the magma due to pressure and temperature changes.

(iii) Heat exchange between hot magma and cold surrounding rocks.

(iv) Fracture tip propagation in host rocks.

The dynamics of dyke propagation involves some of the most difficult problems in mathematical physics (Rivalta et al. 2015) whose solution requires computationally efficient and stable algorithms. Currently, the formulation requires many simplifications.

Early modelling of dyke propagation focused on semi-analytical solutions based on assumptions of a constant overpressure and elliptical shape for fracture opening (Gudmundsson 1983; Lister & Kerr 1991), or considered a toughness-dominated fracture with a linear pressure gradient and no fluid flow (Weertman 1971). These approaches provide rough estimates for dyke geometry, but do not take into account the heterogeneity of the host rock or changes of the physical properties of the magma during ascent.

Magma properties strongly depend on temperature and pressure, thus it is necessary to consider thermal effects during magma ascent and account for non-linear pressure distribution within the dyke. A large body of work investigates magma transport in dykes and conduits during different types of volcanic eruptions after the magma reaches the surface. Most relevant studies come from models of lava dome-building eruptions related to the slow ascent of high-viscosity magma (Costa et al. 2007a). During water-saturated magma ascent, as pressure decreases the viscosity of the magma starts to increase due to water exsolution. This also leads to degassing-induced crystallization because the liquidus temperature of the magma increases progressively with reduced dissolved water (Blundy et al. 2010). There are several major controls on magma temperature evolution: release of the latent heat of crystallization, shear heating, and heat loss to the surrounding rocks (Rubin 1995; Costa & Macedonio 2005; Costa et al. 2007b). The models cited above do not consider the initial stages of eruptions during which magma propagates to the surface. Maimon et al. (2012) consider two-phase magma ascent in a dyke with account for bubble growth due to CO₂ exsolution from kimberlitic magma at large depth. Simulation results show different propagation regimes depending on host rock fracture toughness, with either fracture-controlled or magma-controlled propagation. Formation of separate gas-filled cavity could accelerate propagation in viscosity dominated regime.

Tracking propagation of the fracture front is necessary to understand whether magma ascent will lead to an eruption or will be arrested at some depth and solidify to form an intrusion. It was shown previously that the presence of a stress barrier caused by tectonic stresses or load of a volcanic edifice can lead to a deviation of the dyke from vertical or completely stop its propagation (Dahm 2000; Rivalta et al. 2015). At the level of neutral buoyancy (the depth at which the magma density is equal to the crust density, Walker 1989) the fracture begins to grow in the horizontal direction until the stress intensity factor becomes equal to the material fracture toughness (Chen et al. 2011; Zia & Lecampion 2020).

The mathematical formulation of the problem of magma ascent in a dyke is closely related to hydraulic fracturing (HF) models used in the petroleum industry (Adachi et al. 2007). The simulation of the HF process originates from the 1950s, and many numerical methods have since been developed. The boundary element method (BEM) is well suited to the problem of dyke propagation, since it links the fluid overpressure inside the dyke with its opening using an integral equation. Together with linear elastic fracture mechanics (LEFM), mass conservation inside the dyke and lubrication theory, it allows formulation of a closed system of equations for the propagation of a hydraulic fracture. To construct accurate and efficient numerical methods near-tip asymptotics are used (Peirce & Detournay 2008; Detournay 2016; Dontsov & Peirce 2017; Zia & Lecampion 2020).

In this paper, we present a model of buoyant magma ascent in a dyke taking into account the multiphase nature of the magma, heat transfer, volatile exsolution, and melt crystallization. The model links advanced physics from conduit flow models with a classical HF approach. The numerical algorithm is based on the implicit level set method (Peirce & Detournay 2008), modified to take into account the compressibility of the magma (mass balance instead of volume balance). To consider changes of physical and rheological properties of magma the energy conservation law is incorporated. The effects of exsolution of mixed volatiles (H₂O and CO₂), melt crystallization, viscous heating, heat loss to wall rocks and the release of latent heat of crystallization are all taken into account.

## 2 PROBLEM FORMULATION

This section presents the mathematical formulation of the model for dyke ascent in an elastic medium taking into account the presence of the two principal magmatic volatile species (H₂O and CO₂), crystallization and heat exchange with host rocks. To decouple magma propagation dynamics from magma generation processes, we assume release of a fixed mass of magma from a point source located at depth $x_{ch}$ with a specified discharge rate $Q_{ch}$ for a fixed duration $t_{ch}$. If the timescale of the magma release is significantly shorter than the dyke propagation time, its dynamics are controlled mainly by the mass of the injected magma. The point source condition mimics the accumulation of a melt layer at the top of a crystallizing magma reservoir and its sudden release as soon as the critical buoyancy is reached. Melt can accumulate due to partial crystallization of a sill at the top of a magma reservoir or by percolative ascent of melt through a crystal-rich mush. The latter situation mimics that described by the dynamic hot zone models of Jackson et al. (2018).

The model consists of laws for the conservation of mass, momentum, and energy for the magma inside the dyke coupled by boundary conditions with elastic deformation of host rocks and fracture propagation criteria. The magma is assumed to be a Newtonian fluid with a viscosity that depends on composition, water content, temperature and crystal content. We assume equilibrium crystallization of the magma at all stages of its ascent. Magma flow is assumed to be laminar (as a result of the low Reynolds number).

### 2.1 Geomechanics

We assume that the dyke is a vertical plane fracture growing along the axes Ox as shown on Fig. 1. Ascent of the magma is initiated by influx from the melt layer and is mainly controlled by the buoyancy force and viscous resistance. Magma pressure p inside the dyke can be represented as:

$$p = p_{\text{lith}} + p_{e}, \qquad (1)$$
## 1.1 Dyke propagation models

We start by reviewing dyke propagation models to put the new modelling into context. Dyke propagation modelling involves the coupling of several processes:

(i) Mechanical deformation of host rocks induced by magma pressure on the fracture surfaces.

(ii) Transport of magma within the fracture, taking into account the evolution of physical and rheological properties of the magma due to pressure and temperature changes.

(iii) Heat exchange between hot magma and cold surrounding rocks.

(iv) Fracture tip propagation in host rocks.

The dynamics of dyke propagation involves some of the most difficult problems in mathematical physics (Rivalta et al. 2015) whose solution requires computationally efficient and stable algorithms. Currently, the formulation requires many simplifications.

Early modelling of dyke propagation focused on semi-analytical solutions based on assumptions of a constant overpressure and elliptical shape for fracture opening (Gudmundsson 1983; Lister & Kerr 1991), or considered a toughness-dominated fracture with a linear pressure gradient and no fluid flow (Weertman 1971). These approaches provide rough estimates for dyke geometry, but do not take into account the heterogeneity of the host rock or changes of the physical properties of the magma during ascent.

Magma properties strongly depend on temperature and pressure, thus it is necessary to consider thermal effects during magma ascent and account for non-linear pressure distribution within the dyke. A large body of work investigates magma transport in dykes and conduits during different types of volcanic eruptions after the magma reaches the surface. Most relevant studies come from models of lava dome-building eruptions related to the slow ascent of high-viscosity magma (Costa et al. 2007a). During water-saturated magma ascent, as pressure decreases the viscosity of the magma starts to increase due to water exsolution. This also leads to degassing-induced crystallization because the liquidus temperature of the magma increases progressively with reduced dissolved water (Blundy et al. 2010). There are several major controls on magma temperature evolution: release of the latent heat of crystallization, shear heating, and heat loss to the surrounding rocks (Rubin 1995; Costa & Macedonio 2005; Costa et al. 2007b). The models cited above do not consider the initial stages of eruptions during which magma propagates to the surface. Maimon et al. (2012) consider two-phase magma ascent in a dyke with account for bubble growth due to CO₂ exsolution from kimberlitic magma at large depth. Simulation results show different propagation regimes depending on host rock fracture toughness, with either fracture-controlled or magma-controlled propagation. Formation of separate gas-filled cavity could accelerate propagation in viscosity dominated regime.

Tracking propagation of the fracture front is necessary to understand whether magma ascent will lead to an eruption or will be arrested at some depth and solidify to form an intrusion. It was shown previously that the presence of a stress barrier caused by tectonic stresses or load of a volcanic edifice can lead to a deviation of the dyke from vertical or completely stop its propagation (Dahm 2000; Rivalta et al. 2015). At the level of neutral buoyancy (the depth at which the magma density is equal to the crust density, Walker 1989) the fracture begins to grow in the horizontal direction until the stress intensity factor becomes equal to the material fracture toughness (Chen et al. 2011; Zia & Lecampion 2020).

The mathematical formulation of the problem of magma ascent in a dyke is closely related to hydraulic fracturing (HF) models used in the petroleum industry (Adachi et al. 2007). The simulation of the HF process originates from the 1950s, and many numerical methods have since been developed. The boundary element method (BEM) is well suited to the problem of dyke propagation, since it links the fluid overpressure inside the dyke with its opening using an integral equation. Together with linear elastic fracture mechanics (LEFM), mass conservation inside the dyke and lubrication theory, it allows formulation of a closed system of equations for the propagation of a hydraulic fracture. To construct accurate and efficient numerical methods near-tip asymptotics are used (Peirce & Detournay 2008; Detournay 2016; Dontsov & Peirce 2017; Zia & Lecampion 2020).

In this paper, we present a model of buoyant magma ascent in a dyke taking into account the multiphase nature of the magma, heat transfer, volatile exsolution, and melt crystallization. The model links advanced physics from conduit flow models with a classical HF approach. The numerical algorithm is based on the implicit level set method (Peirce & Detournay 2008), modified to take into account the compressibility of the magma (mass balance instead of volume balance). To consider changes of physical and rheological properties of magma the energy conservation law is incorporated. The effects of exsolution of mixed volatiles (H₂O and CO₂), melt crystallization, viscous heating, heat loss to wall rocks and the release of latent heat of crystallization are all taken into account.

## 2 PROBLEM FORMULATION

This section presents the mathematical formulation of the model for dyke ascent in an elastic medium taking into account the presence of the two principal magmatic volatile species (H₂O and CO₂), crystallization and heat exchange with host rocks. To decouple magma propagation dynamics from magma generation processes, we assume release of a fixed mass of magma from a point source located at depth $x_{ch}$ with a specified discharge rate $Q_{ch}$ for a fixed duration $t_{ch}$. If the timescale of the magma release is significantly shorter than the dyke propagation time, its dynamics are controlled mainly by the mass of the injected magma. The point source condition mimics the accumulation of a melt layer at the top of a crystallizing magma reservoir and its sudden release as soon as the critical buoyancy is reached. Melt can accumulate due to partial crystallization of a sill at the top of a magma reservoir or by percolative ascent of melt through a crystal-rich mush. The latter situation mimics that described by the dynamic hot zone models of Jackson et al. (2018).

The model consists of laws for the conservation of mass, momentum, and energy for the magma inside the dyke coupled by boundary conditions with elastic deformation of host rocks and fracture propagation criteria. The magma is assumed to be a Newtonian fluid with a viscosity that depends on composition, water content, temperature and crystal content. We assume equilibrium crystallization of the magma at all stages of its ascent. Magma flow is assumed to be laminar (as a result of the low Reynolds number).

### 2.1 Geomechanics

We assume that the dyke is a vertical plane fracture growing along the axes Ox as shown on Fig. 1. Ascent of the magma is initiated by influx from the melt layer and is mainly controlled by the buoyancy force and viscous resistance. Magma pressure p inside the dyke can be represented as:

$$p = p_{\text{lith}} + p_e, \qquad (1)$$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic view of a dyke-host rock system. Magma in the dyke is characterized by its viscosity $\mu$, density $\rho$, pressure $p$ and temperature $T$. Rocks are linearly elastic with constant modulus $E$ and $\nu$. Dyke opening $w$ depends on the difference between the magma and lithostatic ($p_{\text{lith}}$) pressures.</figcaption></figure>

where $p_{\text{lith}}$ is the lithostatic pressure given by

$$p_{\text{lith}}(x) = p_{\text{atm}} + g \int_{x}^{0} \rho_{r}(x')dx'. \quad (2)$$

Here $g$ is the gravity acceleration, $\rho_r(x)$ is the density of the surrounding rocks at depth $x$, $p_{\text{atm}}$ is the pressure at the surface of Earth. We use eq. (2) from Tomiya et al. (2010) for a characteristic crust density

$$\rho_r(x) = \rho_1 - r_1 x - (\rho_1 - \rho_0)e^{x/r_2}, \quad (3)$$

where $\rho_1$, $\rho_2$, $r_1$ and $r_2$ are the fitting parameters.

The overpressure $p_e$ related to the fracture opening with width $w$ in an elastic medium can be written as an integral equation (see Muskhelishvili & Radok 1963) of the form:

$$p_e(x, t) = -\frac{E}{4\pi(1 - \nu^2)} \int_{l_b(t)}^{l_b(t)} \frac{\partial w(s, t)}{\partial s} \frac{ds}{s - x}, \quad (4)$$

where $E$ is the Young modulus of the rock, $\nu$ is the Poisson coefficient, $l_b$ and $l_t$ are the depths of fracture fronts at the bottom and top, respectively.

## 2.2 Mass conservation law

Magma is a multiphase medium and its bulk density can be represented by a sum of the bulk densities of individual components

$$\rho = \rho_m + \rho_c + \rho_d + \rho_g, \quad (5)$$

where $\rho_m$, $\rho_c$, $\rho_d$ and $\rho_g$ are the bulk densities of the melt, crystal, dissolved and exsolved gas, respectively. Bulk densities are related to pure component densities as

$$\rho_g = \alpha \rho_g^0, \quad (6)$$

$$\rho_c = (1 - \alpha)\beta\rho_c^0, \quad (7)$$

$$\rho_d = (1 - \alpha)(1 - \beta)c\rho_m^0, \quad (8)$$

$$\rho_m = (1 - \alpha)(1 - \beta)(1 - c)\rho_m^0, \quad (9)$$

Here, $\alpha$ is the volume concentration of the exsolved gas (bubbles), $\beta$ is the volume concentration of the crystal phase in the (bubble-free) melt and $c$ is the mass concentration of the dissolved gas in the melt, $\rho_g^0$ and $\rho_c^0$ are the densities of pure gas and pure crystal phases, respectively. $\rho_m^0$ denotes the mean density of the pure melt phase, which depends on its dissolved volatile content (Lange 2018). For example, dissolving 3 wt per cent water decreases the melt density by approximately 250 kg m$^{-3}$ due to a large partial molar volume of the water in rhyolite (see fig. 5.2 in Lesher & Spera 2015).

To describe mixed H$_2$O–CO$_2$ fluid saturation, we use the Magmatic Sat model of Ghiorso & Gualda (2015) implemented in VESIcal software (Iacovino et al. 2021). For a silicic magma composition typical for porphyritic dykes (Gustafson & Hunt 1975)—at 850 °C the relationship between pressure and the amount of dissolved volatiles is shown in Fig. 2. The presence of a relatively small amount of CO$_2$ dramatically reduces water solubility at high pressures (Mangan et al. 2021). For example, at 900 MPa the magma can dissolve more than 17 wt per cent H$_2$O if no CO$_2$ is present but only half that amount if there is 0.55 wt per cent dissolved CO$_2$. Fig. 2 shows three closed-system degassing paths that start with a volatile-saturated magma at 900 MPa, with each of three curves corresponding to a different initial H$_2$O–CO$_2$ ratio. At high pressures magma degassing occurs mostly by CO$_2$ exsolution, with bubbles containing 70–90 wt per cent CO$_2$. At lower pressures, rapid degassing of water occurs and gas bubbles become progressively more water-rich.

Mass conservation laws for condensed and free gas phases are described by the following equations

$$\frac{\partial (\rho_m + \rho_c) w}{\partial t} + \frac{\partial (\rho_m + \rho_c) vw}{\partial x} = 0, \quad (10)$$

$$\frac{\partial (\rho_d + \rho_g) w}{\partial t} + \frac{\partial (\rho_d v + \rho_g v_g) w}{\partial x} = 0. \quad (11)$$

Here $v_g$ is the free gas velocity that is controlled by magma permeability ($k_m$) because bubble rise in highly viscous magma is negligibly slow (Melnik & Sparks 1999). For the rest of the paper, we assume that $v_g = \nu$ during dyke ascent; the validity of this assumption will be assessed in Section 5. In this case, the mass conservation law for the magma in the dyke becomes:

$$\frac{\partial \rho w}{\partial t} + \frac{\partial \rho q}{\partial x} = \rho_{ch} Q_{ch} \delta(x - x_{ch}), \quad (12)$$

where $Q_{ch}$ is the volumetric flow rate of the magma from the source, and $\delta(x)$ represents the Dirac delta function of the coordinate. The volumetric flow rate is assumed to be constant during the injection period $t_{ch}$ and zero thereafter. The total mass of magma inside the dyke is simply the product of the flow rate $Q_{ch}$, magma density at the source, and injection duration.

The volume flux of the magma in the fracture, $q = vw$, can be expressed by means of lubrication theory

$$q = -\frac{w^3}{12\mu} \left( \frac{\partial p}{\partial x} + \rho g \right). \quad (13)$$

The viscosity $\mu$ of the magma depends mainly on the magma temperature $T$, melt composition $X$ and the crystal content $\beta$ (Costa et al. 2007a). The effect of dissolved CO$_2$ on reducing melt viscosity
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic view of a dyke-host rock system. Magma in the dyke is characterized by its viscosity $\mu$, density $\rho$, pressure $p$ and temperature $T$. Rocks are linearly elastic with constant modulus $E$ and $\nu$. Dyke opening $w$ depends on the difference between the magma and lithostatic ($p_{\text{lith}}$) pressures.</figcaption></figure>

where $p_{\text{lith}}$ is the lithostatic pressure given by

$$p_{\text{lith}}(x) = p_{\text{atm}} + g \int_{x}^{0} \rho_{r}(x')dx'. \quad (2)$$

Here $g$ is the gravity acceleration, $\rho_r(x)$ is the density of the surrounding rocks at depth $x$, $p_{\text{atm}}$ is the pressure at the surface of Earth. We use eq. (2) from Tomiya et al. (2010) for a characteristic crust density

$$\rho_r(x) = \rho_1 - r_1 x - (\rho_1 - \rho_0)e^{x/r_2}, \quad (3)$$

where $\rho_1$, $\rho_2$, $r_1$ and $r_2$ are the fitting parameters.

The overpressure $p_e$ related to the fracture opening with width $w$ in an elastic medium can be written as an integral equation (see Muskhelishvili & Radok 1963) of the form:

$$p_e(x, t) = -\frac{E}{4\pi(1-\nu^2)} \int_{l_b(t)}^{l_b(t)} \frac{\partial w(s, t)}{\partial s} \frac{ds}{s-x}, \quad (4)$$

where $E$ is the Young modulus of the rock, $\nu$ is the Poisson coefficient, $l_b$ and $l_t$ are the depths of fracture fronts at the bottom and top, respectively.

## 2.2 Mass conservation law

Magma is a multiphase medium and its bulk density can be represented by a sum of the bulk densities of individual components

$$\rho = \rho_m + \rho_c + \rho_d + \rho_g, \quad (5)$$

where $\rho_m$, $\rho_c$, $\rho_d$ and $\rho_g$ are the bulk densities of the melt, crystal, dissolved and exsolved gas, respectively. Bulk densities are related to pure component densities as

$$\rho_g = \alpha \rho_g^0, \quad (6)$$

$$\rho_c = (1 - \alpha)\beta\rho_c^0, \quad (7)$$

$$\rho_d = (1 - \alpha)(1 - \beta)c\rho_m^0, \quad (8)$$

$$\rho_m = (1 - \alpha)(1 - \beta)(1 - c)\rho_m^0, \quad (9)$$

Here, $\alpha$ is the volume concentration of the exsolved gas (bubbles), $\beta$ is the volume concentration of the crystal phase in the (bubble-free) melt and $c$ is the mass concentration of the dissolved gas in the melt, $\rho_g^0$ and $\rho_c^0$ are the densities of pure gas and pure crystal phases, respectively. $\rho_m^0$ denotes the mean density of the pure melt phase, which depends on its dissolved volatile content (Lange 2018). For example, dissolving 3 wt per cent water decreases the melt density by approximately 250 kg m⁻³ due to a large partial molar volume of the water in rhyolite (see fig. 5.2 in Lesher & Spera 2015).

To describe mixed H₂O–CO₂ fluid saturation, we use the Magmatic Sat model of Ghiorso & Gualda (2015) implemented in VESIcal software (Iacovino et al. 2021). For a silicic magma composition typical for porphyritic dykes (Gustafson & Hunt 1975)—at 850 °C the relationship between pressure and the amount of dissolved volatiles is shown in Fig. 2. The presence of a relatively small amount of CO₂ dramatically reduces water solubility at high pressures (Mangan et al. 2021). For example, at 900 MPa the magma can dissolve more than 17 wt per cent H₂O if no CO₂ is present but only half that amount if there is 0.55 wt per cent dissolved CO₂. Fig. 2 shows three closed-system degassing paths that start with a volatile-saturated magma at 900 MPa, with each of three curves corresponding to a different initial H₂O–CO₂ ratio. At high pressures magma degassing occurs mostly by CO₂ exsolution, with bubbles containing 70–90 wt per cent CO₂. At lower pressures, rapid degassing of water occurs and gas bubbles become progressively more water-rich.

Mass conservation laws for condensed and free gas phases are described by the following equations

$$\frac{\partial (\rho_m + \rho_c) w}{\partial t} + \frac{\partial (\rho_m + \rho_c) vw}{\partial x} = 0, \quad (10)$$

$$\frac{\partial (\rho_d + \rho_g) w}{\partial t} + \frac{\partial (\rho_d v + \rho_g v_g) w}{\partial x} = 0. \quad (11)$$

Here $v_g$ is the free gas velocity that is controlled by magma permeability ($k_m$) because bubble rise in highly viscous magma is negligibly slow (Melnik & Sparks 1999). For the rest of the paper, we assume that $v_g = \nu$ during dyke ascent; the validity of this assumption will be assessed in Section 5. In this case, the mass conservation law for the magma in the dyke becomes:

$$\frac{\partial \rho w}{\partial t} + \frac{\partial \rho q}{\partial x} = \rho_{ch} Q_{ch} \delta(x - x_{ch}), \quad (12)$$

where $Q_{ch}$ is the volumetric flow rate of the magma from the source, and $\delta(x)$ represents the Dirac delta function of the coordinate. The volumetric flow rate is assumed to be constant during the injection period $t_{ch}$ and zero thereafter. The total mass of magma inside the dyke is simply the product of the flow rate $Q_{ch}$, magma density at the source, and injection duration.

The volume flux of the magma in the fracture, $q = vw$, can be expressed by means of lubrication theory

$$q = -\frac{w^3}{12\mu} \left( \frac{\partial p}{\partial x} + \rho g \right). \quad (13)$$

The viscosity $\mu$ of the magma depends mainly on the magma temperature $T$, melt composition $X$ and the crystal content $\beta$ (Costa et al. 2007a). The effect of dissolved CO₂ on reducing melt viscosity
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Volatile solubility in magma. Dashed lines are isobars from the MagmaSat model with associated pressure marked in MPa. The coloured curves are three closed-system degassing paths starting volatile-saturated at 900 MPa with a mixture of CO₂ and H₂O but with different initial dissolved H₂O contents (see the legend).</figcaption></figure>

is not considered, but in silica-rich melts is likely less than approximately 30 per cent for CO₂ concentrations below 0.5 wt per cent (Brearley & Montana 1989). Bubbles can increase or decrease magma viscosity (Mader et al. 2013; Kolzenburg et al. 2022) but the effect is within one order of magnitude for the vesicularities (α) produced by the model scenarios and is neglected here for model simplicity (see Supporting Information). If, in future, the model were applied to scenarios where the magma erupted rather than being intruded at several km depth, then the effect of bubbles on rheology could become significant. In this study we assume that magma viscosity can be formulated as a product of the melt viscosity and the relative viscosity due to crystals as:

$$ \mu = \mu_m(T, X) \cdot \theta(\beta), \quad (14) $$

where $\mu_m(T, X)$ is the viscosity of the melt phase computed by the model of Giordano et al. (2008). $\theta(\beta)$ is the relative viscosity due to the presence of crystals defined as (Costa 2005)

$$ \theta(\beta) = \phi(\beta) \frac{1 + \left(\frac{\beta}{\beta_*}\right)^{\delta}}{\left(1 - \varepsilon \operatorname{erf}\left\{\frac{\sqrt{\pi\beta}}{2\varepsilon\beta_*}\left[1 + \frac{\beta}{\beta_*}\right]^\gamma\right\}\right)^{2.5\beta_*}} \quad (15) $$

Here, $\phi(\beta)$ is the correction due to the change in the composition of the residual melt parametrized in Dirksen et al. (2006)

$$ \phi(\beta) = \exp[a_1(\beta - \beta_{ch}) + a_2(\beta - \beta_{ch})^2], \quad (16) $$

where $a_1$ and $a_2$ are empirical coefficients and $\beta_{ch}$ is the initial crystal content in the magma source. Parametrization was done for Shiveluch volcano, Kamchatka, for a range of SiO₂ contents from 69 to 77 wt per cent, but should be equally applicable to a porphyry melt composition (see Table 1). A more elaborate approach could incorporate an external thermodynamic package (e.g. MELTS; Gualda et al. 2012) but that would significantly complicate the model and slow down calculations.

<table><caption>**Table 1.** Magma composition (in wt. per cent) for Pinatubo (1991) dacite magma and L-porphyry from El Salvador deposit, Chile (Gustafson & Hunt 1975).</caption><thead><tr><th>Oxide</th><th>Pinatubo dacite, 1991</th><th>L-porphyry</th></tr></thead><tbody><tr><td>SiO<sub>2</sub></td><td>64.6</td><td>65.6</td></tr><tr><td>Al<sub>2</sub>O<sub>3</sub></td><td>16.5</td><td>16.3</td></tr><tr><td>TiO<sub>2</sub></td><td>0.53</td><td>0.71</td></tr><tr><td>FeO<sub>tot</sub></td><td>4.37</td><td>4.14</td></tr><tr><td>MgO</td><td>2.39</td><td>1.6</td></tr><tr><td>CaO</td><td>5.23</td><td>4.34</td></tr><tr><td>Na<sub>2</sub>O</td><td>4.49</td><td>4.89</td></tr><tr><td>K<sub>2</sub>O</td><td>1.54</td><td>1.83</td></tr></tbody></table>

In the case of closed-system degassing (no gas escapes from the ascending magma), the density and viscosity of the magma can be expressed as a function of pressure and temperature only.

To describe the phase relationships of volatile-bearing silicic magma as a function of pressure *p*, temperature *T*, and volatile composition, we have parametrized the extensive volatile-saturated experimental dataset for the 1991 Pinatubo dacite magma (Rutherford & Devine 1996; Scaillet & Evans 1999; Hammer & Rutherford 2002, 2003). Although other compositions could be selected for study, the abundance of experimental data for Pinatubo over a wide p-T-X<sub>H₂O</sub> range is exceptional. To adapt the parametrization to more evolved magma compositions, such as those of dykes associated with some porphyry copper deposits (e.g. Gustafson & Hunt 1975), can be done by simply reducing the liquidus temperature, while maintaining the same solidus temperature and the same functional form of crystal content variation with temperature between liquidus and solidus.

For low-X<sub>H₂O</sub> magmas (shown as numbers on the curves in Fig. 3) solidus and liquidus temperature increase with pressure, while for high-X<sub>H₂O</sub> magmas liquidus temperature decreases sharply with
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Volatile solubility in magma. Dashed lines are isobars from the MagmaSat model with associated pressure marked in MPa. The coloured curves are three closed-system degassing paths starting volatile-saturated at 900 MPa with a mixture of CO₂ and H₂O but with different initial dissolved H₂O contents (see the legend).</figcaption></figure>

is not considered, but in silica-rich melts is likely less than approximately 30 per cent for CO₂ concentrations below 0.5 wt per cent (Brearley & Montana 1989). Bubbles can increase or decrease magma viscosity (Mader et al. 2013; Kolzenburg et al. 2022) but the effect is within one order of magnitude for the vesicularities (α) produced by the model scenarios and is neglected here for model simplicity (see Supporting Information). If, in future, the model were applied to scenarios where the magma erupted rather than being intruded at several km depth, then the effect of bubbles on rheology could become significant. In this study we assume that magma viscosity can be formulated as a product of the melt viscosity and the relative viscosity due to crystals as:

$$ \mu = \mu_m(T, X) \cdot \theta(\beta), \quad (14) $$

where $\mu_m(T, X)$ is the viscosity of the melt phase computed by the model of Giordano et al. (2008). $\theta(\beta)$ is the relative viscosity due to the presence of crystals defined as (Costa 2005)

$$ \theta(\beta) = \phi(\beta) \frac{1 + \left(\frac{\beta}{\beta_*}\right)^{\delta}}{\left(1 - \varepsilon \operatorname{erf}\left\{\frac{\sqrt{\pi\beta}}{2\varepsilon\beta_*}\left[1 + \frac{\beta}{\beta_*}\right]^\gamma\right\}\right)^{2.5\beta_*}} \quad (15) $$

Here, $\phi(\beta)$ is the correction due to the change in the composition of the residual melt parametrized in Dirksen et al. (2006)

$$ \phi(\beta) = \exp[a_1(\beta - \beta_{ch}) + a_2(\beta - \beta_{ch})^2], \quad (16) $$

where $a_1$ and $a_2$ are empirical coefficients and $\beta_{ch}$ is the initial crystal content in the magma source. Parametrization was done for Shiveluch volcano, Kamchatka, for a range of SiO₂ contents from 69 to 77 wt per cent, but should be equally applicable to a porphyry melt composition (see Table 1). A more elaborate approach could incorporate an external thermodynamic package (e.g. MELTS; Gualda et al. 2012) but that would significantly complicate the model and slow down calculations.

<table><caption>**Table 1.** Magma composition (in wt. per cent) for Pinatubo (1991) dacite magma and L-porphyry from El Salvador deposit, Chile (Gustafson & Hunt 1975).</caption><thead><tr><th>Oxide</th><th>Pinatubo dacite, 1991</th><th>L-porphyry</th></tr></thead><tbody><tr><td>SiO<sub>2</sub></td><td>64.6</td><td>65.6</td></tr><tr><td>Al<sub>2</sub>O<sub>3</sub></td><td>16.5</td><td>16.3</td></tr><tr><td>TiO<sub>2</sub></td><td>0.53</td><td>0.71</td></tr><tr><td>FeO<sub>tot</sub></td><td>4.37</td><td>4.14</td></tr><tr><td>MgO</td><td>2.39</td><td>1.6</td></tr><tr><td>CaO</td><td>5.23</td><td>4.34</td></tr><tr><td>Na<sub>2</sub>O</td><td>4.49</td><td>4.89</td></tr><tr><td>K<sub>2</sub>O</td><td>1.54</td><td>1.83</td></tr></tbody></table>

In the case of closed-system degassing (no gas escapes from the ascending magma), the density and viscosity of the magma can be expressed as a function of pressure and temperature only.

To describe the phase relationships of volatile-bearing silicic magma as a function of pressure *p*, temperature *T*, and volatile composition, we have parametrized the extensive volatile-saturated experimental dataset for the 1991 Pinatubo dacite magma (Rutherford & Devine 1996; Scaillet & Evans 1999; Hammer & Rutherford 2002, 2003). Although other compositions could be selected for study, the abundance of experimental data for Pinatubo over a wide p-T-X<sub>H₂O</sub> range is exceptional. To adapt the parametrization to more evolved magma compositions, such as those of dykes associated with some porphyry copper deposits (e.g. Gustafson & Hunt 1975), can be done by simply reducing the liquidus temperature, while maintaining the same solidus temperature and the same functional form of crystal content variation with temperature between liquidus and solidus.

For low-X<sub>H₂O</sub> magmas (shown as numbers on the curves in Fig. 3) solidus and liquidus temperature increase with pressure, while for high-X<sub>H₂O</sub> magmas liquidus temperature decreases sharply with
<figure><img src="image_4.png" /><figcaption>**Figure 3.** In the left image, solidus $T_s$ (dashed line; eq. 18) and liquidus $T_l$ (solid lines; eq. 17) temperatures as a function of pressure and molar H₂O fraction in the coexisting fluid, from 1 = pure H₂O to 0 = pure CO₂. In the right image, crystal content as a function of reduced temperature (eq. 19). See text for data sources and Table 1 for Pinatubo magma composition.</figcaption></figure>

increasing pressure before becoming almost isothermal. This behaviour reflects the pressure dependence of H₂O solubility. CO₂ has little effect on the crystallizing assemblage, but does have the effect of diluting the activity of H₂O in any coexisting fluid and driving earlier volatile exsolution. To a reasonable approximation, the phase relationships of H₂O-undersaturated magmas can be modelled as volatile-saturated magmas coexisting with a mixed H₂O-CO₂ volatile phase. This is the approach that we take with the Pinatubo experimental data.

From a total of 24 experimental constraints over the pressure range 50–960 MPa and XH₂O from 1 to 0.2 the Pinatubo dacite liquidus temperature is parametrized as a function of pressure and water content. The melt composition at the solidus of Pinatubo dacite (and other relatively evolved magmas) approximates that of a near-eutectic granitic melt. Thus, for the solidus we parametrize experimental determinations in the haplogranite system from Holtz et al. (1992), Holtz & Johannes (1994) and Holtz et al. (2001). The resulting expressions are:

$$T_L = a_L + b_L p - X_{H_2O} \left( c_L + f_L p - \frac{d_L}{p+e_L} \right) + \Delta T_L, \quad (17)$$

$$T_S = a_S + b_S p - X_{H_2O} \left( c_S + f_S p - \frac{d_S}{p+e_S} \right). \quad (18)$$

Here $p$ is the pressure in kbar, $X_{H_2O} = C_{H_2O}/(C_{H_2O} + C_{CO_2})$ is the molar fraction of water in the dissolved gas, and $a_L = 1205.7$, $b_L = 6.0$, $c_L = 285.7$, $d_L = 200.0$, $e_L = 0.7$, $f_L = 11.0$, $a_S = 854.1$, $b_S = 6.0$, $c_S = 224.1$, $d_S = 80.0$, $e_S = 0.36$ and $f_S = 6.0$ are the fitting parameters. Where required, a modification to Pinatubo liquidus temperature $\Delta T_L = 1287.6 - 20.15C_{SiO_2}$ °C can be introduced, to account for higher silica content of some porphyry magmas. The effect of varying magma silica content is explored in Supporting Information.

Crystal weight fraction from a total of 86 Pinatubo experiments over the pressure range 5–980 MPa and temperature range 780–900 °C is parametrized as a function of the reduced temperature $T' = (T - T_S)/(T_L - T_S)$ in the form:

$$\beta = \left[ 1 + \exp \left( a_F + b_F T' + c_F T'^2 + d_F T'^3 \right) \right]^{-1} \quad (19)$$

where $a_F = -4.974$, $b_F = 28.623$, $c_F = -52.708$ and $d_F = 34.816$. Crystal fractions were taken from the original experimental studies, where reported, or obtained from mass balance of the analysed

run product phases. Eq. (19) reproduces crystal fractions with an absolute average deviation of 0.065 over the entire P–T range.

The density and specific internal energy of H₂O–CO₂ mixture are calculated with the help of NIST REFPROP database (Lemmon et al. 2018).

## 2.3 Energy conservation law

To calculate the variation of magma temperature $T$, we consider the energy conservation law for magma flow inside the dyke:

$$\frac{\partial}{\partial t}(\rho we) + \frac{\partial}{\partial x}(\rho eq) = 2q_n - p \left( \frac{\partial w}{\partial t} + \frac{\partial q}{\partial x} \right) + \Phi + Q_c. \quad (20)$$

Here, $e(p, T)$ is the internal energy of the magma and defined as a function of temperature $T$ and pressure $p$

$$e(p, T) = \frac{e_g(p, T)\rho_g + C_m(\rho_m + \rho_d)T + C_c\rho_c T}{\rho}, \quad (21)$$

where $e_g(p, T)$ is the specific internal energy of the exsolved gas.

The first term on the right-hand side of eq. (20) is the heat loss due to conduction between the magma and surrounding rocks, which is given by

$$q_n = -k_r \frac{\partial T_r}{\partial y}, \quad (22)$$

where $k_r$ is the rock thermal conductivity, $T_r$ is the temperature of the host rocks, $y$ is a coordinate, perpendicular to the margin of the dyke, that can be assumed horizontal due to large ratio between the length and the width of the dyke. In order to calculate the temperature distribution around the ascending dyke we solve the 1-D-heat conduction equation

$$\rho_r C_r \frac{\partial T_r}{\partial t} - \frac{\partial}{\partial y} \left( k_r \frac{\partial T_r}{\partial y} \right) = 0, \quad (23)$$

in a narrow layer ($0 < y < L_y$, $L_y = \omega\sqrt{\frac{k}{\rho_c C_c}t_{max}}$), where $t_{max}$ is the dyke ascent time, and the value of $\omega$ is chosen to be large enough to ovoid the influence of the outer boundary condition on the temperature distribution near the dyke boundary.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** In the left image, solidus $T_s$ (dashed line; eq. 18) and liquidus $T_l$ (solid lines; eq. 17) temperatures as a function of pressure and molar H₂O fraction in the coexisting fluid, from 1 = pure H₂O to 0 = pure CO₂. In the right image, crystal content as a function of reduced temperature (eq. 19). See text for data sources and Table 1 for Pinatubo magma composition.</figcaption></figure>

increasing pressure before becoming almost isothermal. This behaviour reflects the pressure dependence of H₂O solubility. CO₂ has little effect on the crystallizing assemblage, but does have the effect of diluting the activity of H₂O in any coexisting fluid and driving earlier volatile exsolution. To a reasonable approximation, the phase relationships of H₂O-undersaturated magmas can be modelled as volatile-saturated magmas coexisting with a mixed H₂O-CO₂ volatile phase. This is the approach that we take with the Pinatubo experimental data.

From a total of 24 experimental constraints over the pressure range 50–960 MPa and XH₂O from 1 to 0.2 the Pinatubo dacite liquidus temperature is parametrized as a function of pressure and water content. The melt composition at the solidus of Pinatubo dacite (and other relatively evolved magmas) approximates that of a near-eutectic granitic melt. Thus, for the solidus we parametrize experimental determinations in the haplogranite system from Holtz et al. (1992), Holtz & Johannes (1994) and Holtz et al. (2001). The resulting expressions are:

$$T_L = a_L + b_L p - X_{H_2O} \left( c_L + f_L p - \frac{d_L}{p+e_L} \right) + \Delta T_L, \quad (17)$$

$$T_S = a_S + b_S p - X_{H_2O} \left( c_S + f_S p - \frac{d_S}{p+e_S} \right). \quad (18)$$

Here $p$ is the pressure in kbar, $X_{H_2O} = C_{H_2O}/(C_{H_2O} + C_{CO_2})$ is the molar fraction of water in the dissolved gas, and $a_L = 1205.7$, $b_L = 6.0$, $c_L = 285.7$, $d_L = 200.0$, $e_L = 0.7$, $f_L = 11.0$, $a_S = 854.1$, $b_S = 6.0$, $c_S = 224.1$, $d_S = 80.0$, $e_S = 0.36$ and $f_S = 6.0$ are the fitting parameters. Where required, a modification to Pinatubo liquidus temperature $\Delta T_L = 1287.6 - 20.15C_{SiO_2}$ °C can be introduced, to account for higher silica content of some porphyry magmas. The effect of varying magma silica content is explored in Supporting Information.

Crystal weight fraction from a total of 86 Pinatubo experiments over the pressure range 5–980 MPa and temperature range 780–900 °C is parametrized as a function of the reduced temperature $T' = (T - T_S)/(T_L - T_S)$ in the form:

$$\beta = \left[ 1 + \exp \left( a_F + b_F T' + c_F T'^2 + d_F T'^3 \right) \right]^{-1} \quad (19)$$

where $a_F = -4.974$, $b_F = 28.623$, $c_F = -52.708$ and $d_F = 34.816$. Crystal fractions were taken from the original experimental studies, where reported, or obtained from mass balance of the analysed

run product phases. Eq. (19) reproduces crystal fractions with an absolute average deviation of 0.065 over the entire P–T range.

The density and specific internal energy of H₂O–CO₂ mixture are calculated with the help of NIST REFPROP database (Lemmon et al. 2018).

## 2.3 Energy conservation law

To calculate the variation of magma temperature $T$, we consider the energy conservation law for magma flow inside the dyke:

$$\frac{\partial}{\partial t}(\rho we) + \frac{\partial}{\partial x}(\rho eq) = 2q_n - p \left( \frac{\partial w}{\partial t} + \frac{\partial q}{\partial x} \right) + \Phi + Q_c. \quad (20)$$

Here, $e(p, T)$ is the internal energy of the magma and defined as a function of temperature $T$ and pressure $p$

$$e(p, T) = \frac{e_g(p, T)\rho_g + C_m(\rho_m + \rho_d)T + C_c\rho_c T}{\rho}, \quad (21)$$

where $e_g(p, T)$ is the specific internal energy of the exsolved gas.

The first term on the right-hand side of eq. (20) is the heat loss due to conduction between the magma and surrounding rocks, which is given by

$$q_n = -k_r \frac{\partial T_r}{\partial y}, \quad (22)$$

where $k_r$ is the rock thermal conductivity, $T_r$ is the temperature of the host rocks, $y$ is a coordinate, perpendicular to the margin of the dyke, that can be assumed horizontal due to large ratio between the length and the width of the dyke. In order to calculate the temperature distribution around the ascending dyke we solve the 1-D-heat conduction equation

$$\rho_r C_r \frac{\partial T_r}{\partial t} - \frac{\partial}{\partial y} \left( k_r \frac{\partial T_r}{\partial y} \right) = 0, \quad (23)$$

in a narrow layer ($0 < y < L_y$, $L_y = \omega\sqrt{\frac{k}{\rho_c C_c}t_{max}}$), where $t_{max}$ is the dyke ascent time, and the value of $\omega$ is chosen to be large enough to ovoid the influence of the outer boundary condition on the temperature distribution near the dyke boundary.
Thermal conductivity of the rocks $k_r$ is obtained by the approximation of the measurements from Wen et al. (2015)

$$k_r = \frac{k_0}{1 + B(T - T_0)}, \quad (24)$$

where $k_0$, $B$ and $T_0$ are the fitting parameters (Table 2). Eq. (23) is solved by means of an implicit scheme on a non-uniform mesh to avoid time-step restrictions. It is subject to two boundary and one initial conditions

$$ \begin{aligned} T_r|_{y=0} &= T, \\ T_r|_{y=L_y} &= T_\infty, \\ T_r|_{t=0} &= T_\infty. \end{aligned} \quad (25) $$

We assume that the initial temperature distribution inside the rock is governed by a steady-state solution of the heat conduction equation taking account of radioactive heating:

$$\frac{\partial}{\partial x} \left( k_r \frac{\partial T_r}{\partial x} \right) = A_0 e^{-x/x_*} + A_1. \quad (26)$$

We modified the approach from Chapman (1986) where exponential decay of radioactive heating is assumed in the upper crust ($x < x_*$) and a constant value is postulated beneath. Rock temperature at the surface is assumed to be 0 °C, while at the depth of the magma chamber it is equal to the temperature of the magma. As an extra parameter, we specify the temperature gradient at the surface and modify parameters $A_0$ and $A_1$ in order to keep the temperature distribution monotonic with depth.

The second term on the right-hand side of eq. (20) is the work of the pressure forces associated with the change in density, which is averaged across the dyke width. Shear heating $\Phi$ can be calculated analytically for the parabolic velocity profile and is written as (Rubin 1995)

$$\Phi = \int_{-\frac{w}{2}}^{\frac{w}{2}} \mu \Phi dy = \frac{w^3}{12\mu} \left( \frac{\partial p}{\partial x} + \rho g \right)^2 . \quad (27)$$

Latent heat release due to decompression-induced crystallization can cause magma temperature increase by several tens of degrees during ascent (Blundy et al. 2006). $Q_c$ in eq. (20) is a dyke width averaged release of the latent heat

$$Q_c = \int_{-\frac{w}{2}}^{\frac{w}{2}} L_* \frac{d\rho_c}{dt} dy = L_* \frac{d(\rho_c w)}{dt}. \quad (28)$$

## 2.4 Propagation criteria and boundary conditions

The governing equations are complemented by the boundary conditions at the top $l_b(t)$ and bottom $l_b(t)$ dyke fronts. Due to the assumption that the fracture propagation is determined by linear elastic fracture mechanics (LEFM), in which case the mode I stress intensity factor is equal to the fracture toughness, dyke growth can be described by

$$\lim_{s \to 0} \frac{w(s)}{s^{1/2}} = \frac{K'}{E'}, \quad (29)$$

where $K' = \sqrt{\frac{32}{\pi} K_{lc}}$, $E' = \frac{E}{1 - v^2}$, and $s$ is the distance to the dyke front (Rice et al. 1968). In addition, zero flux conditions at the fracture tips are enforced

$$q(l_t, t) = q(l_b, t) = 0. \quad (30)$$

This means that magma reaches the front of the propagating dyke and no gas cavity exists at the front. This assumption, discussed later, is justified if there is no significant gas flux relative to the magma.

## 3 NUMERICAL SCHEME

In this section, we briefly describe the numerical algorithm for solving the governing equations for dyke ascent on a regular Eulerian mesh (see Supporting Information for details). For this, we use the variant of the implicit level set algorithm (Peirce & Detournay 2008). Each time-step iteration in this algorithm includes two sub-steps: first, we fix the dyke fronts and solve the combined system of eqs (4), (12) and (20). For this, the dyke width $w$ and magma parameters are discretized by a piecewise constant approximation on the mesh elements

$$f(x, t) \approx \sum_{i=1}^{N} f_i(t) H_i(x), \quad (31)$$

$$H_i(x) = \begin{cases} 1, & \text{if } x \in A_i, \\ 0, & \text{otherwise,} \end{cases} \quad (32)$$

where $A_i = [x_{i-1/2}, x_{i+1/2}]$ is a mesh element, $x_i$ and $x_{i \pm 1/2}$ are the centre and boundaries of the element $A_i$, respectively. The numerical solution of this system of equations involves the displacement discontinuous method (DDM, Crouch et al. 1983) for discretizing the elasticity eq. (4) and finite volume method for discretizing mass (eq. 12) and energy (eq. 20) conservation laws. Note that these equations are non-linear partial differential equations; therefore, an iterative process is necessary to find the solution. First, we find dyke width $w$ and magma pressure $p$ by fixing other parameters in the elasticity equation and the mass conservation law. Then, we calculate the magma temperature from the energy conservation law. Finally, knowing the magma pressure $p$ and temperature $T$, we can update magma density $\rho$ and viscosity $\mu$ until convergence is achieved.

In the second substep of the algorithm, we update the position of the front given the new estimate of the dyke width $w$. To obtain accurate solutions on a relatively coarse mesh, the propagation condition (eq. 29) can be replaced by an asymptotic solution $w_a$ (e.g. Peirce & Detournay 2008; Dontsov & Peirce 2017) that has a much larger region of validity

$$w(s) = w_a(s). \quad (33)$$

After convergence of the second sub-step, the code advances to the next step in time. The algorithm was validated against known analytical solutions. The numerical solver algorithm is depicted in Fig. 4. Details can be found in the Supporting Information.

## 4 SIMULATION RESULTS

We first present the results of simulations for the basic set of parameters referenced in Table 2 (Figs 5 and 6, Movie S1) and then analyse the influence of governing parameters on dyke propagation. The dyke propagates by buoyancy, and after release of a fixed mass of magma from the source, the dyke starts to close at depth and magma accumulates at the front resulting in an increase in elastic overpressure. For the reference case scenario (Table 2), the width of the dyke increases towards its top and reaches a maximum of 7.5 m and then rapidly decreases towards the front, where the width is zero (Fig. 5a). Elastic overpressure remains nearly constant in most of the dyke, which keeps it open, and increases near the tip
Thermal conductivity of the rocks $k_r$ is obtained by the approximation of the measurements from Wen et al. (2015)

$$k_r = \frac{k_0}{1 + B(T - T_0)}, \quad (24)$$

where $k_0$, $B$ and $T_0$ are the fitting parameters (Table 2). Eq. (23) is solved by means of an implicit scheme on a non-uniform mesh to avoid time-step restrictions. It is subject to two boundary and one initial conditions

$$\begin{align*}
T_r|_{y=0} &= T, \\
T_r|_{y=L_y} &= T_\infty, \\
T_r|_{t=0} &= T_\infty.
\end{align*} \quad (25)$$

We assume that the initial temperature distribution inside the rock is governed by a steady-state solution of the heat conduction equation taking account of radioactive heating:

$$\frac{\partial}{\partial x} \left( k_r \frac{\partial T_r}{\partial x} \right) = A_0 e^{-x/x_*} + A_1. \quad (26)$$

We modified the approach from Chapman (1986) where exponential decay of radioactive heating is assumed in the upper crust ($x < x_*$) and a constant value is postulated beneath. Rock temperature at the surface is assumed to be 0 °C, while at the depth of the magma chamber it is equal to the temperature of the magma. As an extra parameter, we specify the temperature gradient at the surface and modify parameters $A_0$ and $A_1$ in order to keep the temperature distribution monotonic with depth.

The second term on the right-hand side of eq. (20) is the work of the pressure forces associated with the change in density, which is averaged across the dyke width. Shear heating $\Phi$ can be calculated analytically for the parabolic velocity profile and is written as (Rubin 1995)

$$\Phi = \int_{-\frac{w}{2}}^{\frac{w}{2}} \mu \Phi dy = \frac{w^3}{12\mu} \left( \frac{\partial p}{\partial x} + \rho g \right)^2 . \quad (27)$$

Latent heat release due to decompression-induced crystallization can cause magma temperature increase by several tens of degrees during ascent (Blundy et al. 2006). $Q_c$ in eq. (20) is a dyke width averaged release of the latent heat

$$Q_c = \int_{-\frac{w}{2}}^{\frac{w}{2}} L_* \frac{d\rho_c}{dt} dy = L_* \frac{d(\rho_c w)}{dt} . \quad (28)$$

## 2.4 Propagation criteria and boundary conditions

The governing equations are complemented by the boundary conditions at the top $l_b(t)$ and bottom $l_b(t)$ dyke fronts. Due to the assumption that the fracture propagation is determined by linear elastic fracture mechanics (LEFM), in which case the mode I stress intensity factor is equal to the fracture toughness, dyke growth can be described by

$$\lim_{s \to 0} \frac{w(s)}{s^{1/2}} = \frac{K'}{E'}, \quad (29)$$

where $K' = \sqrt{\frac{32}{\pi} K_{lc}}$, $E' = \frac{E}{1 - v^2}$, and $s$ is the distance to the dyke front (Rice et al. 1968). In addition, zero flux conditions at the fracture tips are enforced

$$q(l_t, t) = q(l_b, t) = 0. \quad (30)$$

This means that magma reaches the front of the propagating dyke and no gas cavity exists at the front. This assumption, discussed later, is justified if there is no significant gas flux relative to the magma.

## 3 NUMERICAL SCHEME

In this section, we briefly describe the numerical algorithm for solving the governing equations for dyke ascent on a regular Eulerian mesh (see Supporting Information for details). For this, we use the variant of the implicit level set algorithm (Peirce & Detournay 2008). Each time-step iteration in this algorithm includes two sub-steps: first, we fix the dyke fronts and solve the combined system of eqs (4), (12) and (20). For this, the dyke width $w$ and magma parameters are discretized by a piecewise constant approximation on the mesh elements

$$f(x, t) \approx \sum_{i=1}^{N} f_i(t) H_i(x), \quad (31)$$

$$H_i(x) = \begin{cases} 1, & \text{if } x \in A_i, \\ 0, & \text{otherwise,} \end{cases} \quad (32)$$

where $A_i = [x_{i-1/2}, x_{i+1/2}]$ is a mesh element, $x_i$ and $x_i \pm 1/2$ are the centre and boundaries of the element $A_i$, respectively. The numerical solution of this system of equations involves the displacement discontinuous method (DDM, Crouch et al. 1983) for discretizing the elasticity eq. (4) and finite volume method for discretizing mass (eq. 12) and energy (eq. 20) conservation laws. Note that these equations are non-linear partial differential equations; therefore, an iterative process is necessary to find the solution. First, we find dyke width $w$ and magma pressure $p$ by fixing other parameters in the elasticity equation and the mass conservation law. Then, we calculate the magma temperature from the energy conservation law. Finally, knowing the magma pressure $p$ and temperature $T$, we can update magma density $\rho$ and viscosity $\mu$ until convergence is achieved.

In the second substep of the algorithm, we update the position of the front given the new estimate of the dyke width $w$. To obtain accurate solutions on a relatively coarse mesh, the propagation condition (eq. 29) can be replaced by an asymptotic solution $w_a$ (e.g. Peirce & Detournay 2008; Dontsov & Peirce 2017) that has a much larger region of validity

$$w(s) = w_a(s). \quad (33)$$

After convergence of the second sub-step, the code advances to the next step in time. The algorithm was validated against known analytical solutions. The numerical solver algorithm is depicted in Fig. 4. Details can be found in the Supporting Information.

## 4 SIMULATION RESULTS

We first present the results of simulations for the basic set of parameters referenced in Table 2 (Figs 5 and 6, Movie S1) and then analyse the influence of governing parameters on dyke propagation. The dyke propagates by buoyancy, and after release of a fixed mass of magma from the source, the dyke starts to close at depth and magma accumulates at the front resulting in an increase in elastic overpressure. For the reference case scenario (Table 2), the width of the dyke increases towards its top and reaches a maximum of 7.5 m and then rapidly decreases towards the front, where the width is zero (Fig. 5a). Elastic overpressure remains nearly constant in most of the dyke, which keeps it open, and increases near the tip
<table><caption>**Table 2.** Parameters used in the simulations.</caption><thead><tr><th>Notation</th><th>Description</th><th>Base value</th><th>Range</th></tr></thead><tbody><tr><td>C<sub>c</sub></td><td>Crystal specific heat [J (kg·K)<sup>-1</sup>]</td><td>1200</td><td>—</td></tr><tr><td>C<sub>m</sub></td><td>Melt specific heat [J (kg·K)<sup>-1</sup>]</td><td>1200</td><td>—</td></tr><tr><td>C<sub>r</sub></td><td>Host rock specific heat [J (kg·K)<sup>-1</sup>]</td><td>1200</td><td>—</td></tr><tr><td>c<sub>H2O</sub></td><td>Concentration of dissolved water in chamber</td><td>6.23 wt%</td><td>3–11 wt%</td></tr><tr><td>E</td><td>Young's modulus of the rocks (GPa)</td><td>15</td><td>—</td></tr><tr><td>K<sub>Ic</sub></td><td>Fracture toughness (MPa·m<sup>1/2</sup>)</td><td>1</td><td>—</td></tr><tr><td>L<sub>*</sub></td><td>Latent heat of crystallization (J kg<sup>-1</sup>)</td><td>350000</td><td>—</td></tr><tr><td>M</td><td>Total magma mass per meter of dyke breadth (M kg m<sup>-1</sup>)</td><td>100</td><td>50–150</td></tr><tr><td>T<sub>ch</sub></td><td>Temperature in the magma source (°C)</td><td>850</td><td>850–900</td></tr><tr><td>x<sub>ch</sub></td><td>Depth of magma source (km)</td><td>30</td><td>—</td></tr><tr><td>ρ<sub>c</sub></td><td>Density of the pure crystal phase (kg m<sup>-3</sup>)</td><td>2700</td><td>—</td></tr><tr><td>ρ<sub>r</sub></td><td>Host rock density (kg m<sup>-3</sup>)</td><td>2700</td><td>—</td></tr><tr><td>ν</td><td>Poisson's ratio of the rocks</td><td>0.25</td><td>—</td></tr><tr><td colspan="4"><strong>Rheological model parameters</strong></td></tr><tr><td>β<sub>*</sub>, γ, δ</td><td>Viscosity change with crystal fraction</td><td>0.67, 3.99, 16.94</td><td>—</td></tr><tr><td>a<sub>1</sub>, a<sub>2</sub></td><td>Viscosity change with melt composition</td><td>4.33, 10.48</td><td>—</td></tr><tr><td colspan="4"><strong>Fitting parameters in eq. (3)</strong></td></tr><tr><td>ρ<sub>0</sub>, ρ<sub>1</sub></td><td></td><td>2092.9, 2732.3</td><td>—</td></tr><tr><td>r<sub>1</sub>, r<sub>2</sub></td><td></td><td>1.918 × 10<sup>-3</sup>, 1724.7</td><td>—</td></tr><tr><td colspan="4"><strong>Fitting parameters in eq. (24)</strong></td></tr><tr><td>k<sub>0</sub>, B, T<sub>0</sub></td><td></td><td>1.9, 2 × 10<sup>-3</sup>, 273.15</td><td>—</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Numerical solver algorithm.</figcaption></figure>

where the buoyancy force increases due to a decrease in the magma density (Fig. 5b). Magmatic temperature decreases progressively from the source region due to heat loss to the host rocks. Near the tip of the dyke, there is a slight increase in temperature due to the larger dyke width that prevents rapid cooling and the shorter time that the magma interacts with cold rocks (Fig. 5c). Magma viscosity increases by 10¹² times from the source to the tip of the dyke due to both reduced temperature and decompression-induced crystallization. The volume fraction of crystals increases up to nearly 80 per cent in the tip region, leading to a dramatic increase in

the relative viscosity given by eq. (15). The low density of silica-rich magmas (e.g. compared to basalt) leads to high buoyancy of the magma at depth. Due to crystallization, the magma density increases in the middle part of the dyke but decreases again at low pressures due to vesiculation (Fig. 5b). The large buoyancy of the magma results in a high propagation speed (Fig. 6) reaching 1.9 m s⁻¹ during the initial propagation period. The increase in magma viscosity due to crystallization leads to deceleration of the dyke despite the buoyancy increase. The magma cools down, crystallizes (driving further volatile exsolution), and becomes extremely viscous. The propagation speed of the dyke decreases dramatically, leading to the final arrest at about 8.7 km depth for the parameter values of the reference case scenario (Table 2).

Our parametric study examines how the variation of the water content $c_{H2O}$, the mass of intruded magma $Q_{ch} \cdot t_{ch}$, starting magma temperature $T_{ch}$ and the rate of injection affect the velocity of the dyke tip propagation and the final position of the top of the dyke. The starting water content of the magma has a strong influence on the propagation speed and the arrest depth (Fig. 7a). We compare the three degassing paths shown in Fig. 2, all starting from the same initial pressure and volatile-saturated, but with different $H_2O-CO_2$ ratios. The amount of $CO_2$ in all cases does not exceed 0.5 wt per cent. A relatively water-poor magma ($C_{H2O} = 3.85$ wt per cent) propagates less quickly and becomes arrested in the deep crust due to its much higher starting melt viscosity. A higher water content leads to faster propagation and more shallow arrest. For each simulation, a fixed mass of magma is released from the source at a constant rate for a finite time. A greater magma mass creates a larger buoyancy force that leads to faster dyke propagation and a greater dyke thickness, suppressing cooling (Fig. 7b). As a result, the magma is arrested at a shallower depth and ascends to its final destination faster. Fig. 7(c) shows the influence of the initial magma temperature. During the initial phase of ascent, propagation velocity is only slightly higher for higher starting temperature because the propagation is mainly influenced by the injection of the magma. At later stages, hotter magmas cool more slowly and become arrested at a more shallow level. The intensity of magma injection controls the
<table><caption>**Table 2.** Parameters used in the simulations.</caption><thead><tr><th>Notation</th><th>Description</th><th>Base value</th><th>Range</th></tr></thead><tbody><tr><td>C<sub>c</sub></td><td>Crystal specific heat [J (kg·K)<sup>-1</sup>]</td><td>1200</td><td>—</td></tr><tr><td>C<sub>m</sub></td><td>Melt specific heat [J (kg·K)<sup>-1</sup>]</td><td>1200</td><td>—</td></tr><tr><td>C<sub>r</sub></td><td>Host rock specific heat [J (kg·K)<sup>-1</sup>]</td><td>1200</td><td>—</td></tr><tr><td>c<sub>H2O</sub></td><td>Concentration of dissolved water in chamber</td><td>6.23 wt%</td><td>3–11 wt%</td></tr><tr><td>E</td><td>Young's modulus of the rocks (GPa)</td><td>15</td><td>—</td></tr><tr><td>K<sub>Ic</sub></td><td>Fracture toughness (MPa·m<sup>1/2</sup>)</td><td>1</td><td>—</td></tr><tr><td>L<sub>*</sub></td><td>Latent heat of crystallization (J kg<sup>-1</sup>)</td><td>350000</td><td>—</td></tr><tr><td>M</td><td>Total magma mass per meter of dyke breadth (M kg m<sup>-1</sup>)</td><td>100</td><td>50–150</td></tr><tr><td>T<sub>ch</sub></td><td>Temperature in the magma source (°C)</td><td>850</td><td>850–900</td></tr><tr><td>x<sub>ch</sub></td><td>Depth of magma source (km)</td><td>30</td><td>—</td></tr><tr><td>ρ<sub>c</sub></td><td>Density of the pure crystal phase (kg m<sup>-3</sup>)</td><td>2700</td><td>—</td></tr><tr><td>ρ<sub>r</sub></td><td>Host rock density (kg m<sup>-3</sup>)</td><td>2700</td><td>—</td></tr><tr><td>ν</td><td>Poisson's ratio of the rocks</td><td>0.25</td><td>—</td></tr><tr><td colspan="4"><strong>Rheological model parameters</strong></td></tr><tr><td>β<sub>*</sub>, γ, δ</td><td>Viscosity change with crystal fraction</td><td>0.67, 3.99, 16.94</td><td>—</td></tr><tr><td>a<sub>1</sub>, a<sub>2</sub></td><td>Viscosity change with melt composition</td><td>4.33, 10.48</td><td>—</td></tr><tr><td colspan="4"><strong>Fitting parameters in eq. (3)</strong></td></tr><tr><td>ρ<sub>0</sub>, ρ<sub>1</sub></td><td></td><td>2092.9, 2732.3</td><td>—</td></tr><tr><td>r<sub>1</sub>, r<sub>2</sub></td><td></td><td>1.918 × 10<sup>-3</sup>, 1724.7</td><td>—</td></tr><tr><td colspan="4"><strong>Fitting parameters in eq. (24)</strong></td></tr><tr><td>k<sub>0</sub>, B, T<sub>0</sub></td><td></td><td>1.9, 2 × 10<sup>-3</sup>, 273.15</td><td>—</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Numerical solver algorithm.</figcaption></figure>

where the buoyancy force increases due to a decrease in the magma density (Fig. 5b). Magmatic temperature decreases progressively from the source region due to heat loss to the host rocks. Near the tip of the dyke, there is a slight increase in temperature due to the larger dyke width that prevents rapid cooling and the shorter time that the magma interacts with cold rocks (Fig. 5c). Magma viscosity increases by 10¹² times from the source to the tip of the dyke due to both reduced temperature and decompression-induced crystallization. The volume fraction of crystals increases up to nearly 80 per cent in the tip region, leading to a dramatic increase in

the relative viscosity given by eq. (15). The low density of silica-rich magmas (e.g. compared to basalt) leads to high buoyancy of the magma at depth. Due to crystallization, the magma density increases in the middle part of the dyke but decreases again at low pressures due to vesiculation (Fig. 5b). The large buoyancy of the magma results in a high propagation speed (Fig. 6) reaching 1.9 m s⁻¹ during the initial propagation period. The increase in magma viscosity due to crystallization leads to deceleration of the dyke despite the buoyancy increase. The magma cools down, crystallizes (driving further volatile exsolution), and becomes extremely viscous. The propagation speed of the dyke decreases dramatically, leading to the final arrest at about 8.7 km depth for the parameter values of the reference case scenario (Table 2).

Our parametric study examines how the variation of the water content $c_{H2O}$, the mass of intruded magma $Q_{ch} \cdot t_{ch}$, starting magma temperature $T_{ch}$ and the rate of injection affect the velocity of the dyke tip propagation and the final position of the top of the dyke. The starting water content of the magma has a strong influence on the propagation speed and the arrest depth (Fig. 7a). We compare the three degassing paths shown in Fig. 2, all starting from the same initial pressure and volatile-saturated, but with different $H_2O-CO_2$ ratios. The amount of $CO_2$ in all cases does not exceed 0.5 wt per cent. A relatively water-poor magma ($C_{H2O} = 3.85$ wt per cent) propagates less quickly and becomes arrested in the deep crust due to its much higher starting melt viscosity. A higher water content leads to faster propagation and more shallow arrest. For each simulation, a fixed mass of magma is released from the source at a constant rate for a finite time. A greater magma mass creates a larger buoyancy force that leads to faster dyke propagation and a greater dyke thickness, suppressing cooling (Fig. 7b). As a result, the magma is arrested at a shallower depth and ascends to its final destination faster. Fig. 7(c) shows the influence of the initial magma temperature. During the initial phase of ascent, propagation velocity is only slightly higher for higher starting temperature because the propagation is mainly influenced by the injection of the magma. At later stages, hotter magmas cool more slowly and become arrested at a more shallow level. The intensity of magma injection controls the
Time = 8.3h

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Reference case scenario after *t* = 8.3 hr from the beginning of magma injection. At this time the dyke velocity decreases to 10⁻³ m s⁻¹ and further propagation becomes impossible. (a) Dyke width and elastic overpressure, (b) magma density and buoyancy, (c) distributions of the temperature and viscosity and (d) crystal and bubble content.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Evolution of the dyke tip velocity and position with time for the reference case scenario.</figcaption></figure>

initial velocity of the dyke tip propagation but not the arrest depth (Fig. 7d) if the same mass of magma is injected and the duration of injection period is shorter than the total ascent duration.

We have shown that the mass of injected magma and its water content play a major role in dyke ascent dynamics. Fig. 8 shows in colour the depth of magma arrest for different values of injected magma mass (*M*) and initial water content (*c*<sub>H₂O</sub>). Black dashed lines show the ascent time to the stagnation point. For small *M*, the arrest depth and the propagation time are strongly influenced by *c*<sub>H₂O</sub> (through the magma viscosity and earlier buoyancy increase by vesiculation at larger *c*<sub>H₂O</sub>). At larger *M* the buoyancy of the magma

remains close to the initial melt buoyancy because we assume constant properties of injected material. An increase in *M* decreases the traveltime of the magma.

## 5 DISCUSSION

Our simulations reveal several general features of the ascent of volatile-saturated felsic magmas from the deep crust. Notably, the buoyancy of wet rhyolitic melts is nearly an order of magnitude higher than of basaltic magmas, which enables rapid ascent despite the relatively high viscosity of felsic magmas. In particular, simulations show that the ascent velocity of water-rich felsic magma can exceed 1–2 m s⁻¹ in agreement with eq. (13). Calculated dyke widths during propagation are greater than the critical values required for dyke propagation estimated by Petford *et al.* (1994) using a simpler model that includes cooling by conduction but without degassing-induced changes in viscosity or buoyancy during ascent. For the initially high magma buoyancy (600 kg m⁻³, see Fig. 5) caused by the low density of water-rich silicic melts, and magma viscosity of 10⁴ Pa·s the critical dyke width estimated by Petford *et al.* (1994) is 2 m. The buoyancy can decrease slightly due to crystallization during magma ascent, but at shallow depth the buoyancy increases again due to significant volatile exsolution. In our simulations, the magma is initially volatile-saturated but lacks an exsolved fluid phase; if there were also bubbles present at source, then the calculated ascent velocities would be even greater, as would the total mass flux of volatiles from the source to the shallow crust. Near the stagnation depth dramatic increase in viscosity significantly increases the critical width for dyke propagation estimated
Time = 8.3h

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Reference case scenario after *t* = 8.3 hr from the beginning of magma injection. At this time the dyke velocity decreases to 10⁻³ m s⁻¹ and further propagation becomes impossible. (a) Dyke width and elastic overpressure, (b) magma density and buoyancy, (c) distributions of the temperature and viscosity and (d) crystal and bubble content.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Evolution of the dyke tip velocity and position with time for the reference case scenario.</figcaption></figure>

initial velocity of the dyke tip propagation but not the arrest depth (Fig. 7d) if the same mass of magma is injected and the duration of injection period is shorter than the total ascent duration.

We have shown that the mass of injected magma and its water content play a major role in dyke ascent dynamics. Fig. 8 shows in colour the depth of magma arrest for different values of injected magma mass (*M*) and initial water content (*c*<sub>H₂O</sub>). Black dashed lines show the ascent time to the stagnation point. For small *M*, the arrest depth and the propagation time are strongly influenced by *c*<sub>H₂O</sub> (through the magma viscosity and earlier buoyancy increase by vesiculation at larger *c*<sub>H₂O</sub>). At larger *M* the buoyancy of the magma

remains close to the initial melt buoyancy because we assume constant properties of injected material. An increase in *M* decreases the traveltime of the magma.

## 5 DISCUSSION

Our simulations reveal several general features of the ascent of volatile-saturated felsic magmas from the deep crust. Notably, the buoyancy of wet rhyolitic melts is nearly an order of magnitude higher than of basaltic magmas, which enables rapid ascent despite the relatively high viscosity of felsic magmas. In particular, simulations show that the ascent velocity of water-rich felsic magma can exceed 1–2 m s⁻¹ in agreement with eq. (13). Calculated dyke widths during propagation are greater than the critical values required for dyke propagation estimated by Petford *et al.* (1994) using a simpler model that includes cooling by conduction but without degassing-induced changes in viscosity or buoyancy during ascent. For the initially high magma buoyancy (600 kg m⁻³, see Fig. 5) caused by the low density of water-rich silicic melts, and magma viscosity of 10⁴ Pa·s the critical dyke width estimated by Petford *et al.* (1994) is 2 m. The buoyancy can decrease slightly due to crystallization during magma ascent, but at shallow depth the buoyancy increases again due to significant volatile exsolution. In our simulations, the magma is initially volatile-saturated but lacks an exsolved fluid phase; if there were also bubbles present at source, then the calculated ascent velocities would be even greater, as would the total mass flux of volatiles from the source to the shallow crust. Near the stagnation depth dramatic increase in viscosity significantly increases the critical width for dyke propagation estimated
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Evolution of the ascending dyke tip position (solid lines) and velocity (dashed lines) as a function of time for different water contents $c_{H_2O}$ of magma (a), total masses of magma $M$ (b), initial temperatures of magma in chamber $T_{ch}$ (c) and activity time of chamber $t_{inj}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Dependence of dyke stagnation depth (colours) on the total masses of magma $M$ and concentration of dissolved water of saturated magma at source. Dashed lines indicate ascent times from the source to the stagnation depth.</figcaption></figure>

by Petford et al. (1994). For a magma viscosity of $10^8$ Pa·s Petford et al. (1994) predict a critical dyke width of >15 m which is larger than that produced in our simulations; further dyke propagation, as predicted by their simplified model, becomes impossible. Our simulations show rapid decrease in the propagation velocity prior to final stagnation of the dyke front.

Other models and experiments that do not account for the variable viscosity of the magma predict acceleration of the dyke in the vicinity of the free surface (Rivalta & Dahm 2006) due to modification of the stress field at shallow depth. Increase in the propagation velocity is observed for the fracture toughness-dominated regime and can cause a twofold to threefold increase in dyke propagation
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Evolution of the ascending dyke tip position (solid lines) and velocity (dashed lines) as a function of time for different water contents $c_{H_2O}$ of magma (a), total masses of magma $M$ (b), initial temperatures of magma in chamber $T_{ch}$ (c) and activity time of chamber $t_{inj}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Dependence of dyke stagnation depth (colours) on the total masses of magma $M$ and concentration of dissolved water of saturated magma at source. Dashed lines indicate ascent times from the source to the stagnation depth.</figcaption></figure>

by Petford et al. (1994). For a magma viscosity of $10^8$ Pa·s Petford et al. (1994) predict a critical dyke width of >15 m which is larger than that produced in our simulations; further dyke propagation, as predicted by their simplified model, becomes impossible. Our simulations show rapid decrease in the propagation velocity prior to final stagnation of the dyke front.

Other models and experiments that do not account for the variable viscosity of the magma predict acceleration of the dyke in the vicinity of the free surface (Rivalta & Dahm 2006) due to modification of the stress field at shallow depth. Increase in the propagation velocity is observed for the fracture toughness-dominated regime and can cause a twofold to threefold increase in dyke propagation
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The ratio of the velocities of gas bubbles and the magma front $v_g/v$, depending on the distance from stagnation point (black line). Blue and red lines in the inset are the magma permeability $\lambda$ and volume concentration of exsolved gas $\alpha$.</figcaption></figure>

speed near the surface. In such a case, eruption at the surface is inevitable once magma reaches a certain depth.

Our simulation parameters mean that the dyke regime is viscous-dominated and rock fracture toughness has no significant influence on dyke propagation dynamics in the range between 1 and 1000 MPa $m^{1/2}$ (see Fig. 5 in Supporting Information). Rivalta et al. (2015), based on interpretation of field data by various authors, propose that the fracture toughness can reach several GPa $m^{1/2}$ for long dykes. The maximum value of 4 GPa $m^{1/2}$ was estimated by Olson (2003) assuming that the aspect ratio of the dyke can be estimated from a solution from Pollard & Segall (1987) that assumes constant pressure distribution along the dyke. This approach was criticized by Gill et al. (2022) who argue that constant pressure along large dykes is not valid if the viscosity of the magma is not negligible. Therefore, extreme values of fracture toughness are not physically possible.

Magma viscosity progressively increases due to cooling, degassing and crystallization. This leads to significant dyke deceleration and, for all modelled scenarios, final arrest in the crust without eruption to the surface, that is the dykes are blind. Thus, blind dykes can originate without the presence of rheological barriers or external loading or decrease in the host rock densities (reaching the level of neutral buoyancy) but as a consequence of rheological stiffening of the magma.

Although the model presented in this manuscript has several advantages compared to those published previously, it still has several significant limitations that should be overcome in future. First, our model does not account for external tectonic stress. Homogeneous compressive or extensional regional stresses can affect magma propagation speed significantly. However, these effects can be easily accounted for by changing effective rock density, and thus the magma buoyancy. The presence of non-homogeneous field stresses or layers with different elastic properties will also affect dyke propagation direction (Rubin 1995). This effect cannot be easily incorporated into the current model, which allows only vertical dyke propagation.

Second, our model does not account for escape of exsolved gas from the ascending magma. Based on the results of the simulations, we show that vertical gas escape through the magma can only become significant shortly prior to dyke arrest and will affect only the vicinity of the dyke tip where magma porosity (bubble volume fraction, see Fig. 9) is large. Gas velocity relative to magma can be

estimated from Darcy's law:

$$V_g = V - \frac{\lambda}{\mu_g} \left( \frac{dp}{dx} - \rho_g g \right), \quad (34)$$

where $\rho_g$ and $\mu_g$ are the density and the viscosity of the gas, respectively. Magma permeability $\lambda$ is a function of the bubble content $\lambda = \lambda_0\alpha^n$. During most of its ascent, magma remains relatively impermeable and only close to the stagnation depth does the volume fraction of bubbles increase significantly due to the effect of second-boiling (crystallization). That opens pathways for gas escape as the magma velocity approaches zero. For all simulations, the distance prior to dyke arrest depth at which vertical gas escape becomes significant remains within a few meters. At this stage of magma ascent, significant fluid release may occur, leading to fracturing and brecciation of the country rocks. Thus, our simulations predict blind, stalled dykes associated with extensive fluid-fracturing features, characteristics that are widely observed around mineralizing dykes in porphyry copper deposits. Horizontal gas escape from the magma into wall-rocks will be controlled mostly by the permeability of the rocks through which the dyke propagates. At shallow depth this process can become significant, but at depths of a few km in a hot ductile environment will play a minor role in dyke ascent dynamics. Our results are in contradiction with Maimon et al. (2012) where the possibility of a gas cavity formation at the top of the dyke is shown for low viscosity kimberlitic magmas where gas can separate from the melt at large depth.

Third, account of cross-dyke viscosity variation is necessary due to its strong temperature dependence. 2-D models for magma flow in fixed conduits (Costa et al. 2007b) show that the velocity profile can be far from parabolic with non-monotonic temperature distribution and formation of shear layers. No dyke propagation model yet accounts for non-parabolic velocity profile.

Forth, 3D effects (Zia & Lecampion 2020; Davis et al. 2023) can strongly influence the horizontal extent of the dyke (assumed infinite in our model) and lead to flow rate instability at the tip of the dyke with formation of finger structures (Touvet et al. 2011). 3D simulations using Pyfrac code (Davis et al. 2023) show that after the initial period, the asymptotic dyke-tip propagation approximates a 2D solution (Lister & Kerr 1991) and the horizontal extent of the dyke becomes more or less constant. Tip instabilities will not affect averaged tip velocity which is mainly controlled by the competition between magma buoyancy and viscous resistance.

Finally, our model does not account for the kinetics of crystallization. At high ascent velocities, crystal content can be far from the equilibrium predicted by eq. (19) (see Fig. 3). Delayed crystallization causes reduced magma viscosity (i.e. a delay in the increase in viscosity), greater propagation velocity, less efficient cooling and, thus leads to a more shallow dyke arrest depth. The textural characteristics of dykes will be profoundly influenced by crystallization kinetics, notably in the development of phenocrysts when crystal growth dominates over nucleation and microlites, when nucleation dominates over growth. Thus, the textures of dyke rocks have the potential to reflect magma ascent dynamics, provided that quantitative account is taken of crystallization kinetics. The diagnostic porphyritic texture of mineralizing dykes in porphyry systems may be significant in this regard.

# 6 CONCLUSION

A model of felsic magma ascent in dykes is developed, accounting for the presence of multiple volatile species ($H_2O$ and $CO_2$), heat exchange with host rocks and magma crystallization. The model could
<figure><img src="image_11.png" /><figcaption>Figure 9. The ratio of the velocities of gas bubbles and the magma front $v_g/v$, depending on the distance from stagnation point (black line). Blue and red lines in the inset are the magma permeability $\lambda$ and volume concentration of exsolved gas $\alpha$.</figcaption></figure>

speed near the surface. In such a case, eruption at the surface is inevitable once magma reaches a certain depth.

Our simulation parameters mean that the dyke regime is viscous-dominated and rock fracture toughness has no significant influence on dyke propagation dynamics in the range between 1 and 1000 MPa $m^{1/2}$ (see Fig. 5 in Supporting Information). Rivalta et al. (2015), based on interpretation of field data by various authors, propose that the fracture toughness can reach several GPa $m^{1/2}$ for long dykes. The maximum value of 4 GPa $m^{1/2}$ was estimated by Olson (2003) assuming that the aspect ratio of the dyke can be estimated from a solution from Pollard & Segall (1987) that assumes constant pressure distribution along the dyke. This approach was criticized by Gill et al. (2022) who argue that constant pressure along large dykes is not valid if the viscosity of the magma is not negligible. Therefore, extreme values of fracture toughness are not physically possible.

Magma viscosity progressively increases due to cooling, degassing and crystallization. This leads to significant dyke deceleration and, for all modelled scenarios, final arrest in the crust without eruption to the surface, that is the dykes are blind. Thus, blind dykes can originate without the presence of rheological barriers or external loading or decrease in the host rock densities (reaching the level of neutral buoyancy) but as a consequence of rheological stiffening of the magma.

Although the model presented in this manuscript has several advantages compared to those published previously, it still has several significant limitations that should be overcome in future. First, our model does not account for external tectonic stress. Homogeneous compressive or extensional regional stresses can affect magma propagation speed significantly. However, these effects can be easily accounted for by changing effective rock density, and thus the magma buoyancy. The presence of non-homogeneous field stresses or layers with different elastic properties will also affect dyke propagation direction (Rubin 1995). This effect cannot be easily incorporated into the current model, which allows only vertical dyke propagation.

Second, our model does not account for escape of exsolved gas from the ascending magma. Based on the results of the simulations, we show that vertical gas escape through the magma can only become significant shortly prior to dyke arrest and will affect only the vicinity of the dyke tip where magma porosity (bubble volume fraction, see Fig. 9) is large. Gas velocity relative to magma can be

estimated from Darcy's law:

$$V_g = V - \frac{\lambda}{\mu_g} \left( \frac{dp}{dx} - \rho_g g \right), \quad (34)$$

where $\rho_g$ and $\mu_g$ are the density and the viscosity of the gas, respectively. Magma permeability $\lambda$ is a function of the bubble content $\lambda = \lambda_0\alpha^n$. During most of its ascent, magma remains relatively impermeable and only close to the stagnation depth does the volume fraction of bubbles increase significantly due to the effect of second-boiling (crystallization). That opens pathways for gas escape as the magma velocity approaches zero. For all simulations, the distance prior to dyke arrest depth at which vertical gas escape becomes significant remains within a few meters. At this stage of magma ascent, significant fluid release may occur, leading to fracturing and brecciation of the country rocks. Thus, our simulations predict blind, stalled dykes associated with extensive fluid-fracturing features, characteristics that are widely observed around mineralizing dykes in porphyry copper deposits. Horizontal gas escape from the magma into wall-rocks will be controlled mostly by the permeability of the rocks through which the dyke propagates. At shallow depth this process can become significant, but at depths of a few km in a hot ductile environment will play a minor role in dyke ascent dynamics. Our results are in contradiction with Maimon et al. (2012) where the possibility of a gas cavity formation at the top of the dyke is shown for low viscosity kimberlitic magmas where gas can separate from the melt at large depth.

Third, account of cross-dyke viscosity variation is necessary due to its strong temperature dependence. 2-D models for magma flow in fixed conduits (Costa et al. 2007b) show that the velocity profile can be far from parabolic with non-monotonic temperature distribution and formation of shear layers. No dyke propagation model yet accounts for non-parabolic velocity profile.

Forth, 3D effects (Zia & Lecampion 2020; Davis et al. 2023) can strongly influence the horizontal extent of the dyke (assumed infinite in our model) and lead to flow rate instability at the tip of the dyke with formation of finger structures (Touvet et al. 2011). 3D simulations using Pyfrac code (Davis et al. 2023) show that after the initial period, the asymptotic dyke-tip propagation approximates a 2D solution (Lister & Kerr 1991) and the horizontal extent of the dyke becomes more or less constant. Tip instabilities will not affect averaged tip velocity which is mainly controlled by the competition between magma buoyancy and viscous resistance.

Finally, our model does not account for the kinetics of crystallization. At high ascent velocities, crystal content can be far from the equilibrium predicted by eq. (19) (see Fig. 3). Delayed crystallization causes reduced magma viscosity (i.e. a delay in the increase in viscosity), greater propagation velocity, less efficient cooling and, thus leads to a more shallow dyke arrest depth. The textural characteristics of dykes will be profoundly influenced by crystallization kinetics, notably in the development of phenocrysts when crystal growth dominates over nucleation and microlites, when nucleation dominates over growth. Thus, the textures of dyke rocks have the potential to reflect magma ascent dynamics, provided that quantitative account is taken of crystallization kinetics. The diagnostic porphyritic texture of mineralizing dykes in porphyry systems may be significant in this regard.

## 6 CONCLUSION

A model of felsic magma ascent in dykes is developed, accounting for the presence of multiple volatile species ($H_2O$ and $CO_2$), heat exchange with host rocks and magma crystallization. The model could
be applied to a range of scenarios relevant to volcanology, petrology and ore genesis. To illustrate the potential of the model, we have focussed on dyke ascent from a very deep (~30 km) volatile-rich source, consistent with geochemical data from intrusions associated with porphyry ore deposits. The simulations show that the competition of high magma buoyancy and high viscosity leads to rapid ascent from depth, followed by deceleration at shallow depth due to volatile exsolution and degassing-induced crystallization. The decrease in ascent velocity dramatically increases heat loss and leads to the dyke arrest and complete solidification. Until the stagnation point, volatile loss from the magma is negligible because the dyke propagation is fast relative to percolative flow of the exsolved volatile phase out of the magma; thus, volatile-rich felsic dykes can efficiently transfer volatiles from deep hot-zone reservoirs to the shallow crust. In that case, dyke arrest and extensive fluid release tend to occur at a similar depth in the crust. This behaviour is consistent with the association of blind felsic dykes with brecciated rocks and mineral veins in porphyry ore systems; the concentration of metal-bearing fluids from a much larger mass of magma than the local dyke intrusions might be promoted by enhanced permeability structures associated with dyke emplacement.

The simulations presented demonstrate that volatile-rich felsic magma can ascend by dyke propagation from 30 to 10 km depth in just hours, orders of magnitude faster than feasible by diapirs or percolative melt flow. This provides a mechanism for rapid transport of magma and volatiles from the deep crust, as proposed on petrological grounds (Richards 2011; Loucks 2014; Rezeau & Jagoutz 2020). What our model does not address is the preconditions that would be required for sufficient mass of buoyant volatile-rich felsic magmas to accumulate at depth and then for the dyke to be triggered. Future developments of the felsic dyke ascent model itself should include crystallization kinetics. Non-equilibrium crystallization could lead to even more rapid magma ascent and subsequent deceleration because delayed crystallization delays increases in magma viscosity. Kinetics also affect the distribution of sizes and shapes of crystals such that the crystal textures of dykes may contain important quantitative information about magma ascent dynamics. Similarly, the emplacement depths of porphyry dykes in mineral systems may contain information about magmatic water contents and magma source depths.

Although our study is focussed on the potential for dykes in porphyry ore-forming systems to involve magmas and volatiles derived from the deep crust, it also has implications for the 1991 dacite eruption of Mount Pinatubo on which the phase relations were parametrized. On the basis of their phase relationships, it has been proposed that Pinatubo dacites originated in the lower crust or uppermost mantle at pressures of 950 MPa before ascending rapidly to the surface and erupting (Prouteau & Scaillet 2003). The dacites contain spinel peridotite xenoliths sourced from mantle rocks at or close to the Moho beneath Luzon (Yoshikawa et al. 2016). The xenoliths show limited evidence of interaction with their host dacite magma, indicative of rapid ascent through the crust. Finally, earthquake hypocentres following the eruption extend down to depths of almost 30 km (Möri & Lecampion 2022) suggesting magma mobility at these depths. Thus, our models indicate that not only can some porphyry-forming magmas be sourced and ascend from lower crustal depths, but so too can arc dacites. What controls why some such magmas stagnate within the shallow crust to form ore deposits, while others erupt at the surface, is likely a result of the complex interplay of the magmatic parameters and processes discussed here. On the basis of our simulations, factors that favor eruption over stalling include higher initial water content, higher

temperature and larger mass of the magma in the dyke. Additional factors could include an exsolved volatile phase at source, lower crystallinity, a hotter geotherm, favourable tectonic stresses and changes in conduit geometry at shallow levels.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

### suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## ACKNOWLEDGMENTS

We thank BHP for financial support. JB acknowledges support from a Royal Society Research Professorship (RP\R1\201048), OM acknowledges support from College de France, PAUSE program (C7H-PUB23A59).

## DATA AVAILABILITY

The MATLAB code and test case scenarios are available from http://github.com/rustamNSU/dikes_article.2023. See README.md for details.

## REFERENCES

* Adachi, J., Siebrits, E., Peirce, A. & Desroches, J., 2007. Computer simulation of hydraulic fractures, *Int. J. Rock Mech. Min. Sci.*, **44**(5), 739–757.

* Annen, C., Blundy, J. & Sparks, R., 2006. The genesis of intermediate and silicic magmas in deep crustal hot zones, *J. Petrol.*, **47**(3), 505–539.

* Blundy, J., Cashman, K. & Humphreys, M., 2006. Magma heating by decompression-driven crystallization beneath andesite volcanoes, *Nature*, **443**(7107), 76–80.

* Blundy, J., Cashman, K. V., Rust, A. & Witham, F., 2010. A case for CO₂-rich arc magmas, *Earth planet. Sci. Lett.*, **290**(3–4), 289–301.

* Brandon, A. D., Creaser, R. A. & Chacko, T., 1996. Constraints on rates of granitic magma transport from epidote dissolution kinetics, *Science*, **271**(5257), 1845–1848.

* Brearley, M. & Montana, A., 1989. The effect of CO₂ on the viscosity of silicate liquids at high pressure, *Geochim. Cosmochim. Acta*, **53**(10), 2609–2616.

* Browne, B. & Szramek, L., 2015. Rates of magma ascent and storage, in *The Encyclopedia of Volcanoes*, pp. 203–214, Elsevier.

* Butters, D., 2022. Microanalysis of zircon-hosted melt inclusions from mineralised and barren magmatic systems, *PhD thesis*, University of Bristol.

* Castro, J.M. & Dingwell, D.B., 2009. Rapid ascent of rhyolitic magma at Chaitén Volcano, Chile, *Nature*, **461**(7265), 780–783.

* Chapman, D., 1986. Thermal gradients in the continental crust, *Geol. Soc. Lond., Spec. Publ.*, **24**(1), 63–70.

* Chen, Z., Jin, Z.-H. & Johnson, S., 2011. Transient dike propagation and arrest near the level of neutral buoyancy, *J. Volc. Geotherm. Res.*, **203**(1), 81–86.

* Chiaradia, M. & Caricchi, L., 2017. Stochastic modelling of deep magmatic controls on porphyry copper deposit endowment, *Sci. Rep.*, **7**(1), doi:10.1038/srep44523.

* Chiaradia, M. & Caricchi, L., 2022. Supergiant porphyry copper deposits are failed large eruptions, *Commun. Earth Environ.*, **3**(1), doi:10.1038/s43247-022-00440-7.

* Costa, A., 2005. Viscosity of high crystal content melts: dependence on solid fraction, *Geophys. Res. Lett.*, **32**(22), doi:10.1029/2005GL024303.
be applied to a range of scenarios relevant to volcanology, petrology and ore genesis. To illustrate the potential of the model, we have focussed on dyke ascent from a very deep (~30 km) volatile-rich source, consistent with geochemical data from intrusions associated with porphyry ore deposits. The simulations show that the competition of high magma buoyancy and high viscosity leads to rapid ascent from depth, followed by deceleration at shallow depth due to volatile exsolution and degassing-induced crystallization. The decrease in ascent velocity dramatically increases heat loss and leads to the dyke arrest and complete solidification. Until the stagnation point, volatile loss from the magma is negligible because the dyke propagation is fast relative to percolative flow of the exsolved volatile phase out of the magma; thus, volatile-rich felsic dykes can efficiently transfer volatiles from deep hot-zone reservoirs to the shallow crust. In that case, dyke arrest and extensive fluid release tend to occur at a similar depth in the crust. This behaviour is consistent with the association of blind felsic dykes with brecciated rocks and mineral veins in porphyry ore systems; the concentration of metal-bearing fluids from a much larger mass of magma than the local dyke intrusions might be promoted by enhanced permeability structures associated with dyke emplacement.

The simulations presented demonstrate that volatile-rich felsic magma can ascend by dyke propagation from 30 to 10 km depth in just hours, orders of magnitude faster than feasible by diapirs or percolative melt flow. This provides a mechanism for rapid transport of magma and volatiles from the deep crust, as proposed on petrological grounds (Richards 2011; Loucks 2014; Rezeau & Jagoutz 2020). What our model does not address is the preconditions that would be required for sufficient mass of buoyant volatile-rich felsic magmas to accumulate at depth and then for the dyke to be triggered. Future developments of the felsic dyke ascent model itself should include crystallization kinetics. Non-equilibrium crystallization could lead to even more rapid magma ascent and subsequent deceleration because delayed crystallization delays increases in magma viscosity. Kinetics also affect the distribution of sizes and shapes of crystals such that the crystal textures of dykes may contain important quantitative information about magma ascent dynamics. Similarly, the emplacement depths of porphyry dykes in mineral systems may contain information about magmatic water contents and magma source depths.

Although our study is focussed on the potential for dykes in porphyry ore-forming systems to involve magmas and volatiles derived from the deep crust, it also has implications for the 1991 dacite eruption of Mount Pinatubo on which the phase relations were parametrized. On the basis of their phase relationships, it has been proposed that Pinatubo dacites originated in the lower crust or uppermost mantle at pressures of 950 MPa before ascending rapidly to the surface and erupting (Prouteau & Scaillet 2003). The dacites contain spinel peridotite xenoliths sourced from mantle rocks at or close to the Moho beneath Luzon (Yoshikawa et al. 2016). The xenoliths show limited evidence of interaction with their host dacite magma, indicative of rapid ascent through the crust. Finally, earthquake hypocentres following the eruption extend down to depths of almost 30 km (Möri & Lecampion 2022) suggesting magma mobility at these depths. Thus, our models indicate that not only can some porphyry-forming magmas be sourced and ascend from lower crustal depths, but so too can arc dacites. What controls why some such magmas stagnate within the shallow crust to form ore deposits, while others erupt at the surface, is likely a result of the complex interplay of the magmatic parameters and processes discussed here. On the basis of our simulations, factors that favor eruption over stalling include higher initial water content, higher

temperature and larger mass of the magma in the dyke. Additional factors could include an exsolved volatile phase at source, lower crystallinity, a hotter geotherm, favourable tectonic stresses and changes in conduit geometry at shallow levels.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

### suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## ACKNOWLEDGMENTS

We thank BHP for financial support. JB acknowledges support from a Royal Society Research Professorship (RP\R1\201048), OM acknowledges support from College de France, PAUSE program (C7H-PUB23A59).

## DATA AVAILABILITY

The MATLAB code and test case scenarios are available from http://github.com/rustamNSU/dikes_article.2023. See README.md for details.

## REFERENCES

* Adachi, J., Siebrits, E., Peirce, A. & Desroches, J., 2007. Computer simulation of hydraulic fractures, *Int. J. Rock Mech. Min. Sci.*, **44**(5), 739–757.

* Annen, C., Blundy, J. & Sparks, R., 2006. The genesis of intermediate and silicic magmas in deep crustal hot zones, *J. Petrol.*, **47**(3), 505–539.

* Blundy, J., Cashman, K. & Humphreys, M., 2006. Magma heating by decompression-driven crystallization beneath andesite volcanoes, *Nature*, **443**(7107), 76–80.

* Blundy, J., Cashman, K. V., Rust, A. & Witham, F., 2010. A case for CO₂-rich arc magmas, *Earth planet. Sci. Lett.*, **290**(3–4), 289–301.

* Brandon, A. D., Creaser, R. A. & Chacko, T., 1996. Constraints on rates of granitic magma transport from epidote dissolution kinetics, *Science*, **271**(5257), 1845–1848.

* Brearley, M. & Montana, A., 1989. The effect of CO₂ on the viscosity of silicate liquids at high pressure, *Geochim. Cosmochim. Acta*, **53**(10), 2609–2616.

* Browne, B. & Szramek, L., 2015. Rates of magma ascent and storage, in *The Encyclopedia of Volcanoes*, pp. 203–214, Elsevier.

* Butters, D., 2022. Microanalysis of zircon-hosted melt inclusions from mineralised and barren magmatic systems, *PhD thesis*, University of Bristol.

* Castro, J.M. & Dingwell, D.B., 2009. Rapid ascent of rhyolitic magma at Chaitén Volcano, Chile, *Nature*, **461**(7265), 780–783.

* Chapman, D., 1986. Thermal gradients in the continental crust, *Geol. Soc., Lond., Spec. Publ.*, **24**(1), 63–70.

* Chen, Z., Jin, Z.-H. & Johnson, S., 2011. Transient dike propagation and arrest near the level of neutral buoyancy, *J. Volc. Geotherm. Res.*, **203**(1), 81–86.

* Chiaradia, M. & Caricchi, L., 2017. Stochastic modelling of deep magmatic controls on porphyry copper deposit endowment, *Sci. Rep.*, **7**(1), doi:10.1038/srep44523.

* Chiaradia, M. & Caricchi, L., 2022. Supergiant porphyry copper deposits are failed large eruptions, *Commun. Earth Environ.*, **3**(1), doi:10.1038/s43247-022-00440-7.

* Costa, A., 2005. Viscosity of high crystal content melts: dependence on solid fraction, *Geophys. Res. Lett.*, **32**(22), doi:10.1029/2005GL024303.
Costa, A. & Macedonio, G., 2005. Viscous heating effects in fluids with temperature-dependent viscosity: triggering of secondary flows, *J. Fluid Mech.*, **540**, 21–38.

Costa, A., Melnik, O. & Sparks, R., 2007a. Controls of conduit geometry and Wallrock elasticity on lava dome eruptions, *Earth planet. Sci. Lett.*, **260**(1), 137–151.

Costa, A., Melnik, O. & Vedeneeva, E., 2007b. Thermal effects during magma ascent in conduits, *J. geophys. Res.*, **112**(B12), doi:10.1029/2007JB004985.

Crouch, S.L., Starfield, A.M. & Rizzo, F., 1983. Boundary element methods in solid mechanics, *J. Appl. Mech.*, **50**(3), 704–705.

Dahm, T., 2000. Numerical simulations of the propagation path and the arrest of fluid-filled fractures in the Earth, *Geophys. J. Int.*, **141**(3), 623–638.

Davis, T., Rivalta, E., Smittarello, D. & Katz, R.F., 2023. Ascent rates of 3-D fractures driven by a finite batch of buoyant fluid, *J. Fluid Mech.*, **954**, doi:10.1017/jfm.2022.986.

Detournay, E., 2016. Mechanics of hydraulic fractures, *Ann. Rev. Fluid Mech.*, **48**(1), 311–339.

Dirksen, O., Humphreys, M., Pletchov, P., Melnik, O., Demyanchuk, Y., Sparks, R. & Mahony, S., 2006. The 2001–2004 dome-forming eruption of Shiveluch Volcano, Kamchatka: observation, petrological investigation and numerical modelling, *J. Volc. Geotherm. Res.*, **155**(3), 201–226.

Dontsov, E. & Peirce, A., 2017. A multiscale implicit level set algorithm (ILSA) to model hydraulic fracture propagation incorporating combined viscous, toughness, and leak-off asymptotics, *Comput. Methods Appl. Mech. Eng.*, **313**, 53–84.

Ghiorso, M.S. & Gualda, G. A.R., 2015. An H₂O–CO₂ mixed fluid saturation model compatible with rhyolite-melts, *Contrib. Mineral. Petrol.*, **169**(6), d oi:10.1007/s00410-015-1141-8.

Gill, S. P.A., Walker, R.J., McCaffrey, K. J.W. & Greenfield, C., 2022. No unique scaling law for igneous dikes, *J. geophys. Res.*, **127**(9), e2022JB024120, doi:10.1029/2022JB024120.

Giordano, D., Russell, J.K. & Dingwell, D.B., 2008. Viscosity of magmatic liquids: a model, *Earth planet. Sci. Lett.*, **271**(1), 123–134.

Gualda, G.A., Ghiorso, M.S., Lemons, R.V. & Carley, T.L., 2012. Rhyolite-melts: a modified calibration of melts optimized for silica-rich, fluid-bearing magmatic systems, *J. Petrol.*, **53**(5), 875–890.

Gudmundsson, A., 1983. Form and dimensions of dykes in eastern Iceland, *Tectonophysics*, **95**(3), 295–307.

Gustafson, L.B. & Hunt, J.P., 1975. The porphyry copper deposit at El Salvador, Chile, *Econ. Geol.*, **70**(5), 857–912.

Hammer, J.E. & Rutherford, M.J., 2002. An experimental study of the kinetics of decompression-induced crystallization in silicic melt, *J. geophys. Res.*, **107**(B1), ECV 8–1–ECV 8-24.

Hammer, J.E. & Rutherford, M.J., 2003. Petrologic indicators of preeruptión magma dynamics, *Geology*, **31**(1), 79–82.

Holtz, F. & Johannes, W., 1994. Maximum and minimum water contents of granitic melts: implications for chemical and physical properties of ascending magmas, *Lithos*, **32**, 149–159.

Holtz, F., Pichavant, M., Barbey, P. & Johannes, W., 1992. Effects of H₂O on liquidus phase relations in the haplogranite system at 2 and 5 kbar, *Am. Mineralog.*, **77**(11-12), 1223–1241.

Holtz, F., Johannes, W., Tamic, N. & Behrens, H., 2001. Maximum and minimum water contents of granitic melts generated in the crust: a reevaluation and implications, *Lithos*, **56**, 1–14.

Iacovino, K., Matthews, S., Wieser, P.E., Moore, G.M. & Bégué, F., 2021. Vesical. Part I: an open-source thermodynamic model engine for mixed volatile (H₂O-CO₂) solubility in silicate melts, *Earth Space Sci.*, **8**(11), e2020EA001584, doi:10.1029/2020EA001584.

Jackson, M., Blundy, J. & Sparks, R., 2018. Chemical differentiation, cold storage and remobilization of magma in the Earth's crust, *Nature*, **564**(7736), 405–409.

Kolzenburg, S., Chevrel, M.O. & Dingwell, D.B., 2022. Magma/suspension rheology, *Rev. Mineral. Geochem.*, **87**(1), 639–720.

Lange, R.A., 2018. The effect of H₂O CO₂ and f on the density and viscosity of silicate melts, in *Volatiles in Magmas*, pp. 331–370, De Gruyter.

Lemmon, E.W., Bell, I., Huber, M.L. & McLinden, M.O., 2018. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 10.0, National Institute of Standards and Technology.

Lesher, C.E. & Spera, F.J., 2015. Thermodynamic and transport properties of silicate melts and magma, in *The Encyclopedia of Volcanoes*, pp. 113–141, Elsevier.

Lister, J.R. & Kerr, R.C., 1991. Fluid-mechanical models of crack propagation and their application to magma transport in dykes, *J. geophys. Res.*, **96**(B6), 10 049–10 077.

Loucks, R., 2014. Distinctive composition of copper-ore-forming arcmagmas, *Aust. J. Earth Sci.*, **61**(1), 5–16.

Loucks, R.R., 2021. Deep entrapment of buoyant magmas by orogenic tectonic stress: its role in producing continental crust, adakites, and porphyry copper deposits, *Earth-Sci. Rev.*, **220**, doi:10.1016/j.earscirev.2021.103744.

Lowenstern, J.B., 1994. Dissolved volatile concentrations in an ore-forming magma, *Geology*, **22**(10), 893–896.

Mader, H., Llewellin, E. & Mueller, S., 2013. The rheology of two-phase magmas: a review and analysis, *J. Volc. Geotherm. Res.*, **257**, 135–158.

Maimon, O., Lyakhovsky, V., Melnik, O. & Navon, O., 2012. The propagation of a dyke driven by gas-saturated magma, *Geophys. J. Int.*, **189**(2), 956–966.

Mangan, M.T., Sisson, T.W., Hankins, W.B., Shimizu, N. & Vennemann, T., 2021. Constraints on deep, CO₂-rich degassing at arc volcanoes from solubility experiments on hydrous basaltic andesite of Pavlof Volcano, Alaska Peninsula, at 300 to 1200 MPa, *Am. Mineral.*, **106**(5), 762–773.

Melnik, O. & Sparks, R.S., 1999. Nonlinear dynamics of lava dome extrusion, *Nature*, **402**(6757), 37–41.

Möri, A. & Lecampion, B., 2022. Three-dimensional buoyant hydraulic fracture growth: constant release from a point source, *J. Fluid Mech.*, **950**, doi:10.1017/jfm.2022.800.

Muskelishvili, N. & Radok, J., 1963. *Some Basic Problems of the Mathematical Theory of Elasticity: Fundamental Equations, Plane Theory of Elasticity, Torsion, and Bending*, P. Noordhoff.

Olson, J.E., 2003. Sublinear scaling of fracture aperture versus length: an exception or the rule?, *J. geophys. Res.*, **108**(B9), doi:10.1029/2001JB000419.

Peirce, A. & Detournay, E., 2008. An implicit level set method for modeling hydraulically driven fractures, *Comput. Methods Appl. Mech. Eng.*, **197**(33), 2858–2885.

Petford, N., Kerr, R.C. & Lister, J.R., 1993. Dike transport of granitoid magmas, *Geology*, **21**(9), 845–848.

Petford, N., Lister, J.R. & Kerr, R.C., 1994. The ascent of felsic magmas in dykes, *Lithos*, **32**(1–2), 161–168.

Petford, N., Cruden, A., McCaffrey, K. & Vigneresse, J.-L., 2000. Granite magma formation, transport and emplacement in the earth's crust, *Nature*, **408**(6813), 669–673.

Pollard, D.D. & Segall, P., 1987. 8 - theoretical displacements and stresses near fractures in rock: with applications to faults, joints, veins, dikes, and solution surfaces, in *Fracture Mechanics of Rock*, pp. 277–349, ed. Atkinson, B.K., Academic Press.

Prouteau, G. & Scaillet, B., 2003. Experimental constraints on the origin of the 1991 Pinatubo Dacite, *J. Petrol.*, **44**(12), 2203–2241.

Rezeau, H. & Jagoutz, O., 2020. The importance of H₂O in arc magmas for the formation of porphyry Cu deposits, *Ore Geol. Rev.*, **126**, doi:10.1016/j.oregeorev.2020.103744.

Rice, J.R. et al., 1968. Mathematical analysis in the mechanics of fracture, *Fracture*, **2**, 191–311.

Richards, J.P., 2011. High sr/y arc magmas and porphyry Cu±Mo±Au deposits: just add water, *Econ. Geol.*, **106**(7), 1075–1081.

Rivalta, E. & Dahm, T., 2006. Acceleration of buoyancy-driven fractures and magmatic dikes beneath the free surface, *Geophys. J. Int.*, **166**(3), 1424–1439.

Rivalta, E., Taisne, B., Bunger, A. & Katz, R., 2015. A review of mechanical models of dike propagation: schools of thought, results and future directions, *Tectonophysics*, **638**, 1–42.
Costa, A. & Macedonio, G., 2005. Viscous heating effects in fluids with temperature-dependent viscosity: triggering of secondary flows, *J. Fluid Mech.*, **540**, 21–38.

Costa, A., Melnik, O. & Sparks, R., 2007a. Controls of conduit geometry and Wallrock elasticity on lava dome eruptions, *Earth planet. Sci. Lett.*, **260**(1), 137–151.

Costa, A., Melnik, O. & Vedeneeva, E., 2007b. Thermal effects during magma ascent in conduits, *J. geophys. Res.*, **112**(B12), doi:10.1029/2007JB004985.

Crouch, S.L., Starfield, A.M. & Rizzo, F., 1983. Boundary element methods in solid mechanics, *J. Appl. Mech.*, **50**(3), 704–705.

Dahm, T., 2000. Numerical simulations of the propagation path and the arrest of fluid-filled fractures in the Earth, *Geophys. J. Int.*, **141**(3), 623–638.

Davis, T., Rivalta, E., Smittarello, D. & Katz, R.F., 2023. Ascent rates of 3-D fractures driven by a finite batch of buoyant fluid, *J. Fluid Mech.*, **954**, doi:10.1017/jfm.2022.986.

Detournay, E., 2016. Mechanics of hydraulic fractures, *Ann. Rev. Fluid Mech.*, **48**(1), 311–339.

Dirksen, O., Humphreys, M., Pletchov, P., Melnik, O., Demyanchuk, Y., Sparks, R. & Mahony, S., 2006. The 2001–2004 dome-forming eruption of Shiveluch Volcano, Kamchatka: observation, petrological investigation and numerical modelling, *J. Volc. Geotherm. Res.*, **155**(3), 201–226.

Dontsov, E. & Peirce, A., 2017. A multiscale implicit level set algorithm (ILSA) to model hydraulic fracture propagation incorporating combined viscous, toughness, and leak-off asymptotics, *Comput. Methods Appl. Mech. Eng.*, **313**, 53–84.

Ghiorso, M.S. & Gualda, G. A.R., 2015. An H₂O–CO₂ mixed fluid saturation model compatible with rhyolite-melts, *Contrib. Mineral. Petrol.*, **169**(6), d oi:10.1007/s00410-015-1141-8.

Gill, S. P.A., Walker, R.J., McCaffrey, K. J.W. & Greenfield, C., 2022. No unique scaling law for igneous dikes, *J. geophys. Res.*, **127**(9), e2022JB024120, doi:10.1029/2022JB024120.

Giordano, D., Russell, J.K. & Dingwell, D.B., 2008. Viscosity of magmatic liquids: a model, *Earth planet. Sci. Lett.*, **271**(1), 123–134.

Gualda, G.A., Ghiorso, M.S., Lemons, R.V. & Carley, T.L., 2012. Rhyolite-melts: a modified calibration of melts optimized for silica-rich, fluid-bearing magmatic systems, *J. Petrol.*, **53**(5), 875–890.

Gudmundsson, A., 1983. Form and dimensions of dykes in eastern Iceland, *Tectonophysics*, **95**(3), 295–307.

Gustafson, L.B. & Hunt, J.P., 1975. The porphyry copper deposit at El Salvador, Chile, *Econ. Geol.*, **70**(5), 857–912.

Hammer, J.E. & Rutherford, M.J., 2002. An experimental study of the kinetics of decompression-induced crystallization in silicic melt, *J. geophys. Res.*, **107**(B1), ECV 8–1–ECV 8-24.

Hammer, J.E. & Rutherford, M.J., 2003. Petrologic indicators of preeruptión magma dynamics, *Geology*, **31**(1), 79–82.

Holtz, F. & Johannes, W., 1994. Maximum and minimum water contents of granitic melts: implications for chemical and physical properties of ascending magmas, *Lithos*, **32**, 149–159.

Holtz, F., Pichavant, M., Barbey, P. & Johannes, W., 1992. Effects of H₂O on liquidus phase relations in the haplogranite system at 2 and 5 kbar, *Am. Mineralog.*, **77**(11-12), 1223–1241.

Holtz, F., Johannes, W., Tamic, N. & Behrens, H., 2001. Maximum and minimum water contents of granitic melts generated in the crust: a reevaluation and implications, *Lithos*, **56**, 1–14.

Iacovino, K., Matthews, S., Wieser, P.E., Moore, G.M. & Bégué, F., 2021. Vesical. Part I: an open-source thermodynamic model engine for mixed volatile (H₂O-CO₂) solubility in silicate melts, *Earth Space Sci.*, **8**(11), e2020EA001584, doi:10.1029/2020EA001584.

Jackson, M., Blundy, J. & Sparks, R., 2018. Chemical differentiation, cold storage and remobilization of magma in the Earth’s crust, *Nature*, **564**(7736), 405–409.

Kolzenburg, S., Chevrel, M.O. & Dingwell, D.B., 2022. Magma/suspension rheology, *Rev. Mineral. Geochem.*, **87**(1), 639–720.

Lange, R.A., 2018. The effect of H₂O CO₂ and f on the density and viscosity of silicate melts, in *Volatiles in Magmas*, pp. 331–370, De Gruyter.

Lemmon, E.W., Bell, I., Huber, M.L. & McLinden, M.O., 2018. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 10.0, National Institute of Standards and Technology.

Lesher, C.E. & Spera, F.J., 2015. Thermodynamic and transport properties of silicate melts and magma, in *The Encyclopedia of Volcanoes*, pp. 113–141, Elsevier.

Lister, J.R. & Kerr, R.C., 1991. Fluid-mechanical models of crack propagation and their application to magma transport in dykes, *J. geophys. Res.*, **96**(B6), 10 049–10 077.

Loucks, R., 2014. Distinctive composition of copper-ore-forming arcmagmas, *Aust. J. Earth Sci.*, **61**(1), 5–16.

Loucks, R.R., 2021. Deep entrapment of buoyant magmas by orogenic tectonic stress: its role in producing continental crust, adakites, and porphyry copper deposits, *Earth-Sci. Rev.*, **220**, doi:10.1016/j.earscirev.2021.103744.

Lowenstern, J.B., 1994. Dissolved volatile concentrations in an ore-forming magma, *Geology*, **22**(10), 893–896.

Mader, H., Llewellin, E. & Mueller, S., 2013. The rheology of two-phase magmas: a review and analysis, *J. Volc. Geotherm. Res.*, **257**, 135–158.

Maimon, O., Lyakhovsky, V., Melnik, O. & Navon, O., 2012. The propagation of a dyke driven by gas-saturated magma, *Geophys. J. Int.*, **189**(2), 956–966.

Mangan, M.T., Sisson, T.W., Hankins, W.B., Shimizu, N. & Vennemann, T., 2021. Constraints on deep, CO₂-rich degassing at arc volcanoes from solubility experiments on hydrous basaltic andesite of Pavlof Volcano, Alaska Peninsula, at 300 to 1200 MPa, *Am. Mineral.*, **106**(5), 762–773.

Melnik, O. & Sparks, R.S., 1999. Nonlinear dynamics of lava dome extrusion, *Nature*, **402**(6757), 37–41.

Möri, A. & Lecampion, B., 2022. Three-dimensional buoyant hydraulic fracture growth: constant release from a point source, *J. Fluid Mech.*, **950**, doi:10.1017/jfm.2022.800.

Muskelishvili, N. & Radok, J., 1963. *Some Basic Problems of the Mathematical Theory of Elasticity: Fundamental Equations, Plane Theory of Elasticity, Torsion, and Bending*, P. Noordhoff.

Olson, J.E., 2003. Sublinear scaling of fracture aperture versus length: an exception or the rule?, *J. geophys. Res.*, **108**(B9), doi:10.1029/2001JB000419.

Peirce, A. & Detournay, E., 2008. An implicit level set method for modeling hydraulically driven fractures, *Comput. Methods Appl. Mech. Eng.*, **197**(33), 2858–2885.

Petford, N., Kerr, R.C. & Lister, J.R., 1993. Dike transport of granitoid magmas, *Geology*, **21**(9), 845–848.

Petford, N., Lister, J.R. & Kerr, R.C., 1994. The ascent of felsic magmas in dykes, *Lithos*, **32**(1–2), 161–168.

Petford, N., Cruden, A., McCaffrey, K. & Vigneresse, J.-L., 2000. Granite magma formation, transport and emplacement in the earth’s crust, *Nature*, **408**(6813), 669–673.

Pollard, D.D. & Segall, P., 1987. 8 - theoretical displacements and stresses near fractures in rock: with applications to faults, joints, veins, dikes, and solution surfaces, in *Fracture Mechanics of Rock*, pp. 277–349, ed. Atkinson, B.K., Academic Press.

Prouteau, G. & Scaillet, B., 2003. Experimental constraints on the origin of the 1991 Pinatubo Dacite, *J. Petrol.*, **44**(12), 2203–2241.

Rezeau, H. & Jagoutz, O., 2020. The importance of H₂O in arc magmas for the formation of porphyry Cu deposits, *Ore Geol. Rev.*, **126**, doi:10.1016/j.oregeorev.2020.103744.

Rice, J.R. et al., 1968. Mathematical analysis in the mechanics of fracture, *Fracture*, **2**, 191–311.

Richards, J.P., 2011. High sr/y arc magmas and porphyry Cu±Mo±Au deposits: just add water, *Econ. Geol.*, **106**(7), 1075–1081.

Rivalta, E. & Dahm, T., 2006. Acceleration of buoyancy-driven fractures and magmatic dikes beneath the free surface, *Geophys. J. Int.*, **166**(3), 1424–1439.

Rivalta, E., Taisne, B., Bunger, A. & Katz, R., 2015. A review of mechanical models of dike propagation: schools of thought, results and future directions, *Tectonophysics*, **638**, 1–42.
Rubin, A.M., 1995. Propagation of magma-filled cracks, *Ann. Rev. Earth planet. Sci.*, **23**(1), 287–336.

Rutherford, M.J., 2008. Magma ascent rates, *Rev. Mineral. Geochem.*, **69**(1), 241–271.

Rutherford, M.J. & Devine, J.D., 1996. Preeruption pressure-temperature conditions and volatiles in the 1991 dacitic magma of Mount Pinatubo, *Fire and Mud: Eruptions and Lahars of Mount Pinatubo, Philippines*, pp. 751–766, eds Newhall, Christopher G. & Punongbayan, Andraymundo S., Univ. Washington Press.

Scaillet, B. & Evans, B.W., 1999. The 15 June 1991 eruption of Mount Pinatubo. I. Phase equilibria and pre-eruption P–T–fO$_2$–fH$_2$O conditions of the dacite magma, *J. Petrol.*, **40**(3), 381–411.

Sillitoe, R.H., 2010. Porphyry copper systems, *Econ. Geol.*, **105**(1), 3–41.

Sparks, R., 2013. Kimberlite volcanism, *Ann. Rev. Earth planet. Sci.*, **41**, 497–528.

Tomiya, A., Takahashi, E., Furukawa, N. & Suzuki, T., 2010. Depth and evolution of a silicic magma chamber: melting experiments on a low-k rhyolite from usu volcano, japan, *J. Petrol.*, **51**(6), 1333–1354.

Touvet, T., Balmforth, N.J., Craster, R.V. & Sutherland, B.R., 2011. Fingering instability in buoyancy-driven fluid-filled cracks, *J. Fluid Mech.*, **672**, 60–77.

Walker, G.P.L., 1989. Gravitational (density) controls on volcanism, magma chambers and intrusions, *Aust. J. Earth Sci.*, **36**(2), 149–165.

Weertman, J., 1971. Theory of water-filled crevasses in glaciers applied to vertical magma transport beneath oceanic ridges, *J. geophys. Res.*, **76**(5), 1171–1183.

Wen, H.,., Xiao, Y. & Deng, J., 2015. Temperature dependence of thermal conductivity, diffusion and specific heat capacity for coal and rocks from coalfield, *Thermochim. Acta*, **619**, 41–47.

Yoshikawa, M. et al., 2016. Aqueous fluids and sedimentary melts as agents for mantle wedge metasomatism, as inferred from peridotite xenoliths at pinatubo and iraya volcanoes, luzon arc, *Philippine Islands*, **262**, 355–368.

Zia, H. & Lecampion, B., 2020. Pyfrac: a planar 3D hydraulic fracture simulator, *Comput. Physics Commun.*, **255**, doi:10.1016/j.cpc.2020.107368.
Rubin, A.M., 1995. Propagation of magma-filled cracks, *Ann. Rev. Earth planet. Sci.*, **23**(1), 287–336.

Rutherford, M.J., 2008. Magma ascent rates, *Rev. Mineral. Geochem.*, **69**(1), 241–271.

Rutherford, M.J. & Devine, J.D., 1996. Preeruption pressure-temperature conditions and volatiles in the 1991 dacitic magma of Mount Pinatubo, *Fire and Mud: Eruptions and Lahars of Mount Pinatubo, Philippines*, pp. 751–766, eds Newhall, Christopher G. & Punongbayan, Andraymundo S., Univ. Washington Press.

Scaillet, B. & Evans, B.W., 1999. The 15 June 1991 eruption of Mount Pinatubo. I. Phase equilibria and pre-eruption P–T–fO$_2$–fH$_2$O conditions of the dacite magma, *J. Petrol.*, **40**(3), 381–411.

Sillitoe, R.H., 2010. Porphyry copper systems, *Econ. Geol.*, **105**(1), 3–41.

Sparks, R., 2013. Kimberlite volcanism, *Ann. Rev. Earth planet. Sci.*, **41**, 497–528.

Tomiya, A., Takahashi, E., Furukawa, N. & Suzuki, T., 2010. Depth and evolution of a silicic magma chamber: melting experiments on a low-k rhyolite from usu volcano, japan, *J. Petrol.*, **51**(6), 1333–1354.

Touvet, T., Balmforth, N.J., Craster, R.V. & Sutherland, B.R., 2011. Fingering instability in buoyancy-driven fluid-filled cracks, *J. Fluid Mech.*, **672**, 60–77.

Walker, G.P.L., 1989. Gravitational (density) controls on volcanism, magma chambers and intrusions, *Aust. J. Earth Sci.*, **36**(2), 149–165.

Weertman, J., 1971. Theory of water-filled crevasses in glaciers applied to vertical magma transport beneath oceanic ridges, *J. geophys. Res.*, **76**(5), 1171–1183.

Wen, H.,., Xiao, Y. & Deng, J., 2015. Temperature dependence of thermal conductivity, diffusion and specific heat capacity for coal and rocks from coalfield, *Thermochim. Acta*, **619**, 41–47.

Yoshikawa, M. et al., 2016. Aqueous fluids and sedimentary melts as agents for mantle wedge metasomatism, as inferred from peridotite xenoliths at pinatubo and iraya volcanoes, luzon arc, *Philippines, Lithos*, **262**, 355–368.

Zia, H. & Lecampion, B., 2020. Pyfrac: a planar 3D hydraulic fracture simulator, *Comput. Physics Commun.*, **255**, doi:10.1016/j.cpc.2020.107368.