## Abstract

Modern geophysical data recorded during lava dome building eruptions indicate the presence of multiple connected magma storage regions. Most numerical models for lava dome eruptions assume a single magma chamber fed from below with a constant or prescribed time dependent influx rate and connected to the Earth surface through a conduit.

Here we present a development of the model of extrusive eruptions considering a system made of two magma chambers located in elastic rocks and connected by a dyke between each other. We use locations and volumes of the magma chambers inferred for Soufrière Hills Volcano, Montserrat from ground deformation studies and seismic tomography.

The model shows cyclic behaviour with a period that depends on the intensity of the influx rate, the volumes and shape of the two chambers and the degree of connectivity of the two reservoirs. For a weak connectivity the overpressure in the lower chamber stays nearly constant during the cycle and the influx of fresh magma into the shallow chamber is also nearly constant. For a strong connectivity between the chambers their overpressure increases or decreases during the cycle in a synchronous way.

Lava dome eruptions are typically characterized by large fluctuations in discharge rate showing cyclic patterns of seismicity, ground deformation and volcanic activity. Timescales associated with these fluctuations range from hours–days (short-term cycles) to years–decades (long-term cycles) with intermediate fluctuations of 5–7 weeks as, for example, observed on the Soufrière Hills Volcano (SHV), Montserrat. Concerning short-term cycles, cyclic patterns of seismicity, ground deformation and volcanic activity have been documented at Mount Pinatubo, the Philippines, in 1991 and SHV, Montserrat, British West Indies, in 1996–2010. At SHV, periodicity in seismicity and tilt ranged from approximately 3 to 30 h (Voight *et al.* 1999; Odbert & Wadge 2009).

Several models have been proposed to explain the observed cyclicity in the last decade. Concerning short-term cyclicity, in both the Denlinger & Hoblitt (1999) and Wylie *et al.* (1999) models, the lower part of the conduit acts like a capacitor that allows magma to be stored temporarily in order to be released during the intense phase of the eruption. Costa *et al.* (2012) generalized the Denlinger & Hoblitt (1999) model to a more realistic system made of a compressible magma flowing through an elastic dyke joined above to a shallower cylindrical conduit. Lensky *et al.* (2008) explained short-term cyclic behaviour at the SHV as being a result of gas diffusion into growing bubbles and filtration through the bubble network in a stagnated magma column before the critical overpressure is reached followed by magma motion, depressurization and stagnation of the plug at the top of the conduit.

Regarding long-term cycles, Costa *et al.* (2007*a*, *b*) considered an elastic dyke joined above to a cylinder and coupled below with a shallow (*c.* 5 km) magma chamber. They showed that degassing-induced crystallization within the conduit coupled with the wall-rock elasticity play a major control on processes for the long-term cyclicity during lava-dome-building eruptions. In particular, they showed that there is a regime where the period of pulsations is controlled by the elasticity of the dyke (*c.* weeks to months) and a regime where the period is controlled by the volume of the magma chamber (*c.* years). Intermediate regimes are possible. The geometry proposed by Costa *et al.* (2007*a*, *b*) is also consistent with inversion of strain data for the March 2004 SHV explosion that highlight the presence of a spherical pressure source and a shallow dyke (Linde *et al.* 2010; Hautmann *et al.* 2009).

Most of the existing conduit flow models for lava-dome-building eruptions assume a system consisting of a magma chamber (a large volume of magma located in elastic wall rocks) connected to the surface through a conduit. The magma chamber is assumed to be fed from below at a constant rate with new magma. In a few models (Dirksen *et al.* 2006; Maeda 2000), the influx rate into the magma chamber is specified as a function of time.

Modern geophysical data recorded during lava-dome-building eruptions indicate the presence of interconnected multiple magma storage regions. For the SHV system, different geometries have been proposed on the basis of both mineralogical and geodetic observations. Elsworth *et al.* (2008) proposed a model with two stacked magma reservoirs connected by vertical conduits (the shallow chamber was centred at 6 km and the deep chamber at 12 km below sea level). Surface efflux and global positioning system (GPS) station velocities were inverted to calculate crustal magma transfer. Results from the inversion also indicated that the eruptive episodes deplete the lower reservoir only, not the upper reservoir.

A vertically elongated chamber idealized as a prolate ellipsoid was proposed by Voight *et al.* (2010). This chamber was assumed stratified, with upper parts compressible owing to exsolved gas phases (magma sponge), and fed at the base by a deep influx estimated to be of the order of 2 m^{3} s^{−1}. The top of the ellipsoid was assumed to be 5 km below sea level and its centroid near 10 km.

Whether the SHV system is better represented by a single prolate chamber or a two-chamber model is still not resolved. However, Foroozan *et al.* (2010) showed that the deformations at GPS stations located near or far from the volcano are better explained by a two-chamber model.

In a more recent paper, Foroozan *et al.* (2011) simultaneously inverted geodetic data and measured magma efflux rates with the additional constraint of constant basal influx. In contrast to Elsworth *et al.* (2008), they used data from 6 to 13 stations (4 previously) and a compressible magma phase (incompressible previously). Their results suggest a basal supply of 1.2 m^{3} s^{−1} for the activity from 1995 to 2008. Centroids of the chambers were estimated to depths of about 5 and 19 km. In contrast to Elsworth *et al.* (2008), Foroozan *et al.* (2011) found that the data are characterized by synchronous deflation then inflation of both shallow and deep chambers.

In order to account for the large-strain effects and modulus softening that are unrepresented in the extrapolated one-dimensional (1D) seismic velocity data, Foroozan *et al.* (2011) decreased the shear modulus of the wall rock, *G*, by a factor of 10, estimated on the basis of a 1D seismic velocity profile (Shalev *et al.* 2010). Rather than this arbitrary approach, we will use, following Costa *et al.* (2009, 2011), an empirical relationship between static and dynamic values of the Young modulus, *E*, of rocks (Wang 2000):

However, the analysis of data for the period 2003–2010 from the borehole strainmeter network at SHV gives good constraints on where the pressure changes occur. These observations are also consistent with a system made of a deep reservoir coupled with a shallow chamber connected to the surface through a shallow dyke (Sacks *et al.* 2011).

Most numerical models of these eruptions assume a single magma chamber fed from below with a constant or prescribed time dependent influx rate and connected to the Earth's surface through a conduit. Here we present a development of the model of extrusive eruption by Costa *et al.* (2007*a*, *b*), considering a system made of two magma chambers located in elastic rocks and connected with a dyke between the two. The deep chamber is fed from below with a constant influx rate. For the sake of comparison with previous studies (e.g. Hautmann *et al.* 2009, 2010), we use locations and volumes of the magma chambers inferred for SHV, Montserrat, from ground deformation studies by Hautmann *et al.* (2010) and also the more recent estimations obtained by integrating seismic tomography with numerical models of the shallow magma chamber by Paulatto *et al.* (2012). The estimations of both chamber locations and volumes contain large uncertainties that need to be systematically investigated, but here we will focus on exploring the control of two magma chambers on the dynamics of lava dome extrusion only. In order to study the role of a dual-magma-chamber system, however, we explore a large range of chamber volumes and two different chamber shapes (spheroidal and prolate) for both deep and shallow magma chambers.

The model shows cyclic behaviour with a period that depends on the intensity of the influx rate and the chamber connectivity (described as the horizontal extent of the dyke connecting the two chambers). For a weak connectivity, the overpressure in the lower chamber stays nearly constant during the cycle and the influx of fresh magma into the shallow chamber is also nearly constant. For a strong connectivity between the chambers, their overpressures increase or decrease during the cycle in a synchronous way. Influx into the shallow chamber stays close to the extrusion rate at the surface.

## A model for magma transfer between multiple reservoirs and dykes

We consider a system like that proposed by Hautmann *et al.* (2010) and represented in Figure 3.1. The system is composed of a deep magma chamber fed from below and connected to a shallow magma chamber through a dyke. The shallow magma chamber is connected to the surface via a volcanic conduit comprising a dyke that is joined above to a cylinder, as described by Costa *et al.* (2007*a*, *b*).

The system of equations for the description of magma transport from the deep to the shallow chamber is:
(1)
Here μ denotes magma viscosity, Δ*P* = *p*−*p*_{e}(*z*) is a magmatic overpressure, *p* is the pressure inside the dyke, *p*_{e}(*z*)=ρ_{r}*gz*−σ_{t} denotes the lithostatic pressure (*g* is acceleration due to gravity) minus σ_{t}; that is, tensile stress along the dyke due to the presence of magma chamber and an extensional/compressional far-field stress. *S* is the cross-section area of the dyke, which is assumed to have an elliptical shape, *a* and *b* are the semi-axis of the ellipse. Wall rocks are assumed to be purely elastic with Young's modulus, *E*, and Poisson's ratio, ν. The Muskhelishvili solution (Muskhelishvili 1963) is used for calculation of the relationship between magmatic overpressure and the deformation of the dyke. Densities of magma, ρ, and wall rocks, ρ_{r}, in the deeper part of the magmatic system, are assumed to be constant, *z* is a vertical coordinate increasing upwards, *V* is the ascent velocity. In order to simplify the equation system, since typically *a*_{0}≫*b*_{0}, relative changes of *a* are negligible. Thus, we can assume *a* is constant and the initial value of *b* is equal to 0; that is, *b*_{0}=0. Moreover, at typical depths of several kilometres, we can reasonably neglect the presence of the free gas phase and assume incompressible magma of constant viscosity. For CO_{2}-rich magma, this assumption will be violated owing to deep exsolution of dissolved volatiles. However, for the case of SHV, even assuming the upper-bound estimation of CO_{2} content in the magma (600 ppm: Edmonds *et al.* 2014) is exsolved, magma density variation is negligible; that is, from 2520 to 2600 kg m^{−3} when pressure changes from 240 to 300 MPa (moreover, for such an CO_{2} content at 300 MPa, the magma bulk modulus is *c.* 3 GPa, only few times smaller than rock rigidity at that depth). Hence, we have:
(2)
where α=(1−ν)*a*_{0}/*E*. The above system can be rewritten as a non-linear parabolic equation:
(3)
where Ω=∂/∂*z*[*p*_{e}(*z*) + ρ*gz*]. In order to solve equation (3) we impose the following boundary conditions:
with *Q*=π*a*_{0}*bV* = −(π*a*_{0}*b*/4μ)(∂/∂*z*)[Δ*P*+*p*_{e}(*z*)+ρ*gz*] being the volumetric outflux from the deep magma chamber and *Q*_{in,d} being the influx into deep magma chamber from a deeper magmatic source. *Q*_{in,d} will be considered as a given constant parameter in this study, although it can be influenced by pressure variations in the magmatic source, local tectonic stress and thermal development of the feeding system.

Equation (3) is solved together with the system of equations that describes the flow from the shallow magma chamber through a dyke-shaped conduit that evolves to a cylinder in the vicinity of the surface. These equations are described in detail in Costa *et al.* (2007*a*, *b*) and are summarized in the Appendix.

## Parametric study

Here we present results of two separate studies. First, we analyse the effects of a far-field stress and a pressurized deep magma chamber on the steady-state solution of the flow equation in the upper conduit. Second, assuming a neutral far-field stress, we present results on the role of the connectivity between deep and shallow magma chambers during an extrusive eruption on the transient behaviour of the system.

### Effect of far-field stress

The presence of a deep reservoir affects the dynamics and the cyclicity of the shallow part of the system, but it also has an important effect in perturbing the stress field around it. A pressurized magma chamber under the effect of a far-field stress can modify the tensile stress along the axis of the conduit σ_{t} (see Costa *et al.* 2011). This stress is calculated using the general analytical solution by Gao (1996) obtained in the limit of a 2D plain geometry (for the limitations of this assumption, see Costa *et al.* 2011). We explore the effect of far-field stresses, σ_{ff}, typical of the arc volcanism environment, from compressional (−10 MPa) to extensional (10 MPa) (Roman *et al.* 2006; Hautmann *et al.* 2010), on the steady-state solution of the flow equation in the upper conduit.

Figure 3.2 reports a set of steady-state solutions for a fixed deep reservoir overpressure of 10 MPa and shows the control of the far-field stress field on the system sketched in Figure 3.1. A compressional far-field stress tends to move the steady-state solution towards larger pressures; whereas an extensional far-field stress moves the solution to lower pressures, allowing larger flow rates at low pressures. The conduit cross-sectional area is larger in the case of the extensional environment, thus the transition to the upper branch moves to higher values of discharge rate because ascent velocity for large σ_{ff} is smaller for the same discharge rate and crystals grow more efficiently (see fig. 12b in Melnik & Sparks 2005 for a more detailed explanation).

### Effect of coupling the two reservoirs

In this subsection we explore the role of the connectivity between deep and shallow magma chambers during an extrusive eruption. We consider only the effect due to the two pressurized magma chambers and will set σ_{ff}=0.

As an initial condition, we assume that both of the magma chambers have zero overpressure with respect to the lithostatic pressure, the flow rate is zero everywhere and the dyke between the two chambers has a negligible thickness. At time *t=*0, influx of magma into the deep magma chamber starts with a constant intensity, *Q*_{in,d}. Pressure in the deep chamber starts to grow, leading to magma rising and opening of the pre-existing dyke through the rocks between the two chambers. When the dyke reaches the shallow chamber, the influx of the fresh magma leads to pressure growth and initiates the process of the dome extrusion.

Figure 3.3 shows the evolution of the chamber overpressures normalized to their maximum values (upper plots), and the influx into the shallow magma chamber, *Q*_{in,s}, and discharge rate at the surface, *Q*_{out} (lower plots), for two cases: (i) weak connectivity between the chambers (when the value of *a*_{0} in equation 2 is small, i.e. 100 m); and (ii) strong connectivity (when the value of *a*_{0} in equation 2 is large, i.e. 600 m).

In the case of weak connectivity, the overpressure in the deep chamber increases to a high value (*c.* 70 MPa) and stays nearly constant, while the overpressure in the shallow chamber shows a pronounced cyclic activity. In this case, influx into the shallow chamber is almost constant, although *Q*_{out} varies.

In the case of strong connectivity, both chambers inflate and deflate simultaneously. Influx rate into the shallow chamber increases when the pressure in the shallow magma chamber decreases. *Q*_{out} follows changes in *Q*_{in,s}. For the intermediate values of *a*_{0}, the system stays in-between these two end-member cases. The case of strong connectivity is consistent with reconstructions of Foroozan *et al.* (2011). Behaviour for the weak connectivity is not supported by the recent analysis and interpretation of geodetic observations (Foroozan *et al.* 2011). In this case, high overpressure in the chamber suggests a horizontal extension of the dyke and a consequent improvement of connectivity between the two magma chambers.

In order to quantify the inflation/deflation phase shift between two chambers, we introduce a variable characterizing the value of offset as:
where *T* is the period of discharge rate variation. In the case when the two chambers work in phase, the value of the offset is equal to 0, while an offset of 1 corresponds to completely opposite directions of chamber evolution. Figure 3.4 shows the dependence of the offset on the intensity of deep magma chamber influx, *Q*_{in,d}. Curves are shown for the values of *a*_{0}=100 (short dashed line), 200, 300, 400, 500 and 600 m (solid line). The offset progressively decreases with an increase in the horizontal dyke extent, *a*_{0}, and decreases with an increase in *Q*_{in,d} because, for larger influx intensity, the pressure in the deep chamber reaches higher values (see Fig. 3.5) and the dyke cross-sectional area increases. This improves the connectivity between the chambers.

Figure 3.5 shows the maximum overpressure *P*_{d} that is reached in the deep chamber (a) and the amplitude of overpressure variation during the cycle Δ*P*_{d}; (b) as a function of *Q*_{in,d} for different *a*_{0}. Small values of *a*_{0} result in unrealistically high overpressures and small amplitude of pressure variations. However, in real situations if the overpressure increases above the rock strength, it is expected that a fracture will propagate at the tips of the dyke, thus increasing the dyke dimension and moving the system to lower, more realistic pressures. An increase in *Q*_{in,d} leads to an improvement in the connectivity between the chambers and, thus, increases the amplitude.

The time lag between the beginning of the influx into the deep chamber and the start of the eruption is shown in Figure 3.6 as a function of *Q*_{in,d} for different values of *a*_{0}. It varies from approximately 3 months to more than 1 year. Larger influx intensities and better chamber connectivity lead in shortening of this interval. Mattioli *et al.* (2010) showed that ground deformation and extrusion of the magma at the surface are closely correlated in time. This supports a hypothesis of good chamber connectivity for the SHV.

### Period of discharge rate variations

In the case of weak connectivity between the chambers, for fixed values of the shallow and deep magma chamber volumes (e.g. *V*_{ch,s}=4 km^{3} and *V*_{ch,d}=32 km^{3}), the period of pulsations (Fig. 3.7) reaches a minimum value for some intermediate value of *Q*_{in,d} corresponding to an unstable branch of the steady-state solution similar to the case of a single magma chamber (see fig. 7 in Melnik & Sparks 2005). When the connectivity is strong, the period monotonically increases with increase in the influx rate into the deep chamber.

In order to evaluate the role of magma chamber volumes on the period of discharge rate variation, we performed a series of runs varying the aspect ratios and volumes of both shallow (Fig. 3.8) and deep magma chambers (Fig. 3.9).

Informed by the volume estimates of Paulatto *et al.* (2012), we varied the volume of the shallow magma chamber to between 4 and 30 km^{3} for different values of the deep magma chamber (from 30 to 60 km^{3}), and fixing *Q*_{in,d}=2 m^{3} s^{−1}, a value that represents a realistic estimation for the deep magma influx (Voight *et al.* 2010; Odbert & Wadge 2009). Keeping fixed the volume of the deep chamber to *V*_{ch,d}=32 km^{3}, for *V*_{ch,s}≤20 km^{3} the period monotonically increases from 125 to 150 days, but above 20 km^{3} there is a bifurcation and the solution becomes multiperiodic. The case of *V*_{ch,s}=30 km^{3} results in a double period behaviour, with periods for small-amplitude oscillations of the order of 150 days, whereas the large peak of increase in discharge rate occurs every 400 days (Fig. 3.8). Similar double periodicity behaviour was discovered first in Costa *et al.* (2007*b*).

Figure 3.9 shows the dependence of the period and the shape of the pulsations as a function of the deep magma chamber volume. The volume of the shallow magma chamber is fixed at 15 km^{3} in accord with Paulatto *et al.* (2012). Two shapes of the shallow magma chamber (spherical and prolate, aspect ratio of 1 and 1.5, respectively) are studied. For small volumes of the deep magma chamber (*V*_{ch,d}<37 km^{3}), the period does not strongly depend upon the shallow chamber geometry, and pulsations have a unimodal shape pattern. For larger magma chambers, the pattern during the cycle consists of peaks in the discharge rate separated by intervals of periodic variation in the discharge rate with much lower amplitude that characterize the major part of the cycle.

We have discovered that a prolate shape of the shallow magma chamber leads to much longer periods of pulsations due to a stronger influence of the chamber pressure on the opening of the dyke that connects it to the surface. The opposite effect was observed for the aspect ratio of the deep magma chamber. A more spherical deep magma chamber tends to produce fluctuations with longer periodicity.

## Implications for SHV

Observations at SHV suggest long-term cyclic behaviour with a period of 2–3 years. This period, for example, can be obtained for a prolate shallow magma chamber with a volume of approximately 15 km^{3} and a deep magma chamber of about 50 km^{3}. Figure 3.10 represents the variation in extrusion rate with time for these parameters.

Although with these magma chamber values we can obtain reasonable periodicity, the simulated patterns of the flow rate variations, alternations between phases of high and low discharge rate, are not in agreement with observations for the SHV. Typical behaviour of the volcano involves repetitive periods of extrusion separated by repose periods with nearly no activity. Moreover, during the first three cycles, the periods of extrusion were longer (i.e. years) than the fourth and fifth cycles in which extrusion lasted for months, similar to the simulations showed in Figure 3.10. The behaviour of the first three cycles is better reproduced assuming a narrow shallow dyke, *d*_{s} (*a*_{0} of *c.* 150 m), as shown in Figure 3.11. A systematic parametric investigation also involving the shallow part of the system is an object of ongoing research.

It is also likely that magma rheology can play a key role in the behaviour of the extrusion pattern. For example, we would expect that for a Bingham rheology of the magma, the discharge rate between the two major pulses is zero until a critical chamber overpressure is reached (Melnik & Sparks 2005). Moreover, a more realistic description of the effective viscosity of magma during lava dome eruptions should account for the coupling with energy loss, viscous dissipation, and 2D and stick–slip effects (e.g. Costa & Macedonio 2005; Costa *et al.* 2007*c*, 2012). However, owing to the limitations of the current transient numerical code, we cannot investigate the effects of complex rheologies.

## Discussion and conclusion

The models show that the presence of two magma chambers connected by an elastic dyke significantly influences dynamics of magma extrusion in comparison with a case of a single magma chamber fed with a constant influx rate. The dynamics is strongly influenced by the amount of connectivity between the chambers measured as the length of the large semi-axis of the dyke, *a*_{0}. When *a*_{0} is small (weak connectivity) the pressure in the deep chamber reaches high values and the influx into the shallow chamber is close to a constant value. This case, in terms of surface extrusion, is similar to the case of a single magma chamber, but the geodetic signal in this case will be different. At the initial stage of the eruption, the deep source will inflate producing widespread deformation pattern. Later, the pressure in the deep reservoir will come to a high constant value and ground deformation will be controlled by changes in pressure in the shallow chamber.

In the case of strong connectivity (large *a*_{0}), the initial ground deformation pattern will be controlled by an inflation of the deep magma chamber. Later in the eruptions both pressure sources work in phase. For an intermediate *a*_{0}, there is a transition between two end-members.

It is clear that the assumption of a constant value of *a*_{0} greatly oversimplifies the physics. In the case of large overpressures, stress at the dyke tips will exceed the fracture toughness of the rocks and the dyke will expand horizontally. The maximum dyke extent will be controlled by horizontal chamber dimensions (Costa *et al.* 2011). When the deep magma chamber deflates, overpressure in the dyke will decrease and, as flow rate decreases, magma at the dyke tips can solidify, leading to a decrease in *a*_{0} (Kavanagh & Sparks 2011).

Another strong simplification of the model comes from assumed Newtonian rheology of the magma. Although the system produces periodic behaviour, the shape of discharge rate evolution during the cycle is different from that observed at the SHV. Measurements of the dome volume suggest relatively slow initiation of the eruption, with progressively increasing or nearly constant extrusion rate at later stages of the cycle followed by a rapid drop in discharge rate leading to stagnation of the eruption. Cycles produced by our model show a rapid increase in discharge rate, with more gradual reduction in eruption intensity. In the simulations with a shallow dyke of about 500 m, as in Costa *et al.* (2007*a*), discharge rate remains always larger than zero leading to continued extrusion of magma, whereas simulations using a narrower shallow dyke (*c.* 300 m) show better agreement with observations.

Finally, in order to understand whether simulations can give insights about the issue of a single versus a dual magma chamber system, we carried out some simulations assuming only the presence of the shallow chamber. Results show that for a shallow dyke of approximately 300 m, we can have 2–3 years periodicity assuming a chamber volume of about 30 km^{3}, implying that model results alone do not allow us to discriminate between a dual- and a single-chamber geometry.

Predictions of the model must be verified against measurements of ground deformation, dome volume, gas flux and other geophysical data that provide constraints on the possible scenarios of the eruption. Only a robust approach that is based on an integrated system of geophysical observations will be able to give strong constraints to the model, and to give some insight into the geometry and dynamics of the system.

## Acknowledgments

O. Melnik acknowledges the support from the ERC grant no. 228064 (VOLDIES) and the Russian Foundation for Basic Research (12–01-00465). A. Costa was supported by NERC research grant reference NE/H019928/1. We thank R. S. J. Sparks and G. Wadge for very useful comments on an earlier version of the manuscript.

## Appendix

### Governing equations in the shallow conduit (elliptical dyke to cylinder)

The ascent of magma is simulated by 1D transport equations. We assume that the conduit has an elliptical cross-section with area of *S*=*ab*, with *a* and *b* the major and minor semi-axes, respectively. Their vertical variation occurs at length-scales that are much larger than the dyke width. The set of cross-section averaged equations is presented below:
(A1)
(A2)
(A3)
(A4)
(A5)
Here *t* denotes time, *x* the vertical coordinate, ρ_{m}, ρ_{ph}, ρ_{mc}, ρ_{d} and ρ_{g} are the densities of melt, phenocrysts, microlites, dissolved gas and exsolved gas, respectively, and *V* and *V*_{g} are the velocities of magma and gas, respectively. *G*_{ph} and *G*_{mc} represent the mass transfer rate due to crystallization of phenocrysts and microlites, respectively, and *J* is the mass transfer rate due to gas exsolution. Equation (A1) represents the mass conservation for the melt phase, equations (A2) and (A3) are the conservation equations for microlites and phenocrysts, respectively, and equations (A4) and (A5) represent the conservation of the dissolved gas and of the exsolved gas, respectively.

Two momentum equations for the mixture as a whole (equation A6) and for the free gas phase (equation A7) describe the magma motion:
(A6)
(A7)
Here *p* is the pressure, ρ is the bulk density of magma, *g* is the acceleration due to gravity, *μ* is the magma viscosity, *k* is the magma permeability and *μ*_{g} is the gas viscosity.

The main source of the magma temperature change is the release of the latent heat of crystallization (equation A8).
(A8)
Here *C*_{m} is the bulk specific heat of magma, *T* is the bulk flow-averaged temperature and *L*_{*} is the latent heat of crystallization (viscous heating and heat flux to the host rocks are neglected).For details of equation of state and mass fluxes related to gas exsolution and crystallization see Costa *et al.* (2007*a*, *b*).

For parameterizations of magma permeability, *k*, and magma viscosity, *μ*, we use:
(A9)
(A10)
where *k* depends on the volume fraction of bubbles, α; μ is given by the product of (1) melt viscosity μ_{m} which depends on the water content, *c*, and temperature, *T*, and (2) a function describing viscosity variation due to crystal content, *β*, and a correction for the effect of bubble fraction, α, and capillary number, *Ca* (Costa *et al.* 2007*a*, *b*).

Variations of dyke semi-axis *a* and *b* with pressure are characterized by equation (1). In order to get a smooth transition from the dyke at depth to a cylindrical conduit, the value of *a*_{0} is parameterized as:
(A11)
Here *L*_{T} and *w*_{T} are the height and the vertical extent of the cylinder to ellipse transition zone, and constants *A*_{1} and *A*_{2} are calculated to satisfy conditions *a*_{0}(*L*)=*R* and *a*_{0}(0)=*a*_{0}, where *R* is the radius of the cylindrical part of the conduit and *a*_{0} is the length of major dyke semi-axis at the inlet of the dyke. The value of *b*_{0} is calculated in order to conserve the cross-section area of the unpressurized dyke, although it can also be specified independently.

- © The Geological Society of London 2014