2.1 Linearized response of a fluid parcel

The response to small-amplitude perturbations about an equilibrium or steady-state solution can be studied using linearized versions of the governing equations. We thus write fluid variables velocity, pressure, density, and gas fraction as the sum of a background value, denoted with a bar, and a perturbation:

(2.16) $$\begin{eqnarray}[V,P,{\it\rho},n_{g}]=[\bar{v},\bar{p},\bar{{\it\rho}},\bar{n}]+[v,p,{\it\rho}^{\prime },n],\end{eqnarray}$$

and we linearize the governing equations around a background or base state. As mentioned earlier, perturbations in conduit area $A$ associated with elastic deformation of the conduit walls by pressure changes are neglected in this study.

We begin by considering a fluid parcel initially at thermodynamic equilibrium at pressure $P=\bar{p}$ , for which the initial gas mass fraction is $\bar{n}=n_{eq}(\bar{p})$ and the initial density is, from definition (2.9), $\bar{{\it\rho}}={\it\rho}_{eq}(\bar{p})$ . The more general problem of perturbations about a non-equilibrium background state is left for future work.

Linearizing the equation of state (2.6) yields

(2.17) $$\begin{eqnarray}\displaystyle {\it\rho}^{\prime } & = & \displaystyle \frac{\partial {\it\rho}}{\partial P}(P-\bar{p})+\frac{\partial {\it\rho}}{\partial n_{g}}(n_{g}-\bar{n})\nonumber\\ \displaystyle & = & \displaystyle \bar{{\it\rho}}(p/K-an),\end{eqnarray}$$

using the definitions of the bulk modulus $K$ in (2.11) and $a$ in (2.12). These and all subsequent coefficients are understood to be evaluated in the base state. The equilibrium gas mass fraction is similarly linearized,

(2.18) $$\begin{eqnarray}\displaystyle n_{eq}(P) & = & \displaystyle \bar{n}+\frac{\partial n}{\partial P}(P-\bar{p})\nonumber\\ \displaystyle & = & \displaystyle \bar{n}-bp,\end{eqnarray}$$

using (2.8) to define $b$ . In the context of this linearization, and assumption of perturbations about a state of thermodynamic equilibrium, the evolution equation (2.15) is

(2.19) $$\begin{eqnarray}{\dot{m}}=-\frac{n+bp}{{\it\tau}}.\end{eqnarray}$$

Similarly, the equilibrium density (2.9) is linearized:

(2.20) $$\begin{eqnarray}\displaystyle {\it\rho}_{eq}(P) & = & \displaystyle \bar{{\it\rho}}+\frac{\text{d}{\it\rho}_{eq}}{\text{d}P}(P-\bar{p})\nonumber\\ \displaystyle & = & \displaystyle \bar{{\it\rho}}+\frac{\bar{{\it\rho}}}{K}(1+abK)p=\bar{{\it\rho}}\left(1+\frac{p}{K_{eq}}\right)\!,\end{eqnarray}$$

using (2.13) and the identity

(2.21) $$\begin{eqnarray}\frac{1}{K_{eq}}=\frac{1}{K}+ab.\end{eqnarray}$$

Taking the total time derivative of (2.17), and making use of the equilibrium variables and the linearized evolution equation in (2.18)–(2.21), we eliminate $n$ and obtain a single equation describing the evolution of density under conditions of changing pressure:

(2.22) $$\begin{eqnarray}\frac{\text{D}{\it\rho}}{\text{D}t}=\frac{\bar{{\it\rho}}}{K}\frac{\text{D}P}{\text{D}t}-\frac{{\it\rho}-{\it\rho}_{eq}(P)}{{\it\tau}},\end{eqnarray}$$

or in terms of density and pressure perturbations,

(2.23) $$\begin{eqnarray}\frac{\text{D}{\it\rho}^{\prime }}{\text{D}t}=\frac{\bar{{\it\rho}}}{K}\frac{\text{D}p}{\text{D}t}-\frac{{\it\rho}^{\prime }-\bar{{\it\rho}}p/K_{eq}}{{\it\tau}}.\end{eqnarray}$$

The significance of the two terms on the right-hand side of (2.22) or (2.23) can be understood by considering the response of the fluid parcel to a step increase in pressure (figure 2). That increase in pressure compresses both the gas and liquid phases, leading to an instantaneous increase in density described by the first term and quantified via the bulk modulus  $K$ . Mass exchange between the phases is negligible during this initial compression. However, the increase in pressure alters the equilibrium solubility and drives volatiles back into the melt over the time scale  ${\it\tau}$ . The second term captures this mass transfer and the associated reduction in mixture density towards the new, larger, equilibrium density.

Figure 2.

Response of gas mass fraction $n$ and density ${\it\rho}$ to a step increase in pressure from (2.4) and (2.22). Gas fraction relaxes to a new equilibrium value over time scale ${\it\tau}$ , while density changes include both the instantaneous response due to fluid compression and slower relaxation response from mass exchange between phases.

2.2 Energetics of a fluid parcel

Introduction of a finite adjustment time ${\it\tau}$ necessarily implies energy dissipation. To illustrate this, and to prepare ourselves for later discussion of the energy balance for waves in magmatic conduits, we first consider the energetics of a fluid parcel. In the context of our linearized model, the rate of work per unit reference volume of fluid done by the pressure change $p$ is

(2.24) $$\begin{eqnarray}-\frac{p}{1/\bar{{\it\rho}}}\frac{\text{D}(1/{\it\rho})}{\text{D}t}=\frac{p}{\bar{{\it\rho}}}\frac{\text{D}{\it\rho}^{\prime }}{\text{D}t},\end{eqnarray}$$

where $1/{\it\rho}$ is the specific volume. With no addition of heat, the first law of thermodynamics states that the work done on the fluid parcel is equal to the change in internal energy. Here we show that the internal energy change consists of a recoverable part, the free energy associated with compression of the fluid, and dissipation from non-equilibrium BGR. We anticipate that the free energy per unit reference volume will be $p^{2}/2K$ under non-equilibrium conditions (when mass exchange is negligible) and $p^{2}/2K_{eq}$ under equilibrium conditions.

Substituting the linearized evolution equation (2.23) into (2.24) and using the identity (2.21), we obtain

(2.25) $$\begin{eqnarray}\frac{p}{\bar{{\it\rho}}}\frac{\text{D}{\it\rho}^{\prime }}{\text{D}t}=\frac{\text{D}}{\text{D}t}\left[\frac{p^{2}}{2K}+\frac{({\it\rho}^{\prime }-\bar{{\it\rho}}p/K)^{2}}{2ab\bar{{\it\rho}}^{2}}\right]+\frac{({\it\rho}^{\prime }-\bar{{\it\rho}}p/K_{eq})^{2}}{ab\bar{{\it\rho}}^{2}{\it\tau}}.\end{eqnarray}$$

The quantity in brackets is the total free energy density associated with compression of the fluid. The first part, $p^{2}/2K$ , is the non-equilibrium free energy density, which is the only contribution in the absence of mass exchange (in which case ${\it\rho}^{\prime }=\bar{{\it\rho}}p/K$ ). The second part thus corresponds to an additional contribution to the free energy that arises as mass transfer lowers the fluid bulk modulus from $K$ towards $K_{eq}$ and additional work is done by pressure acting on the fluid parcel. The last term on the right-hand side of (2.25) is the rate of energy dissipation per unit reference volume. It is evident that dissipation occurs only when conditions depart from equilibrium. During equilibrium, ${\it\rho}^{\prime }=\bar{{\it\rho}}p/K_{eq}$ and dissipation vanishes. When BGR is limited by the ability of bubbles to expand or contract in a viscous melt, this energy loss is the usual viscous dissipation. Dissipation can also arise from diffusive processes (e.g. Landau & Lifshitz Reference Landau and Lifshitz1959, Ch. VI), such as when unequal chemical potentials drive transport of volatiles through the melt towards or away from the bubble walls.

We can also rewrite (2.25) by eliminating the density perturbation ${\it\rho}^{\prime }$ in favour of the gas mass fraction perturbation $n$ using (2.17):

(2.26) $$\begin{eqnarray}\frac{p}{\bar{{\it\rho}}}\frac{\text{D}{\it\rho}^{\prime }}{\text{D}t}=\frac{\text{D}}{\text{D}t}\left(\frac{p^{2}}{2K}+\frac{an^{2}}{2b}\right)+\frac{a(n+bp)^{2}}{b{\it\tau}}.\end{eqnarray}$$

The terms have the same interpretation, as $n=0$ in the absence of mass exchange and $n=-bp$ under equilibrium conditions. Furthermore, it is straightforward to verify that under equilibrium conditions the two free energy terms combine to $p^{2}/2K_{eq}$ , the equilibrium free energy density, as anticipated.

2.3 Perturbations about magmastatic base state

Having examined the response of an isolated fluid parcel, we now turn our attention to the response of the magma column. For all subsequent analysis we consider perturbations around a specific base state in which the magma is in both thermodynamic and mechanical equilibrium in the presence of a gravitational field. We refer to this as the magmastatic base state, for which

(2.27) $$\begin{eqnarray}\displaystyle & \bar{v}=0, & \displaystyle\end{eqnarray}$$

(2.28) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\bar{p}}{\text{d}z}=-\bar{{\it\rho}}g, & \displaystyle\end{eqnarray}$$

(2.29) $$\begin{eqnarray}\displaystyle & \bar{n}=n_{eq}(\bar{p}), & \displaystyle\end{eqnarray}$$

and hence from (2.29), with the definition (2.9),

(2.30) $$\begin{eqnarray}\bar{{\it\rho}}={\it\rho}_{eq}(\bar{p}).\end{eqnarray}$$

The magmastatic state implies a vertical density stratification that can be obtained by differentiating (2.30) with respect to $z$ and using (2.13), (2.14) and (2.28):

(2.31) $$\begin{eqnarray}\frac{\text{d}\bar{{\it\rho}}}{\text{d}z}=-\frac{\bar{{\it\rho}}g}{C_{eq}^{2}}.\end{eqnarray}$$

Similarly, differentiating (2.29) with respect to $z$ and using (2.8) and (2.28) yields

(2.32) $$\begin{eqnarray}\frac{\text{d}\bar{n}}{\text{d}z}=\bar{{\it\rho}}gb.\end{eqnarray}$$

Stratified fluids exhibit oscillations for which the restoring force is gravity, especially at length scales exceeding the scale height for stratification, $C_{eq}^{2}/g$ (Lighthill Reference Lighthill1978). To quantify these oscillations, we define the buoyancy or Brunt–Väisälä frequency

(2.33) $$\begin{eqnarray}N_{b}=\sqrt{-\frac{g}{\bar{{\it\rho}}}\frac{\text{d}\bar{{\it\rho}}}{\text{d}z}}=\frac{g}{C_{eq}}.\end{eqnarray}$$

This definition of the buoyancy frequency is valid when fluid compressibility is negligible. As we later show, this is an appropriate approximation in our problem at wavelengths that are much greater than the scale height.

Introducing the magmastatic base state into the linearized governing equations reveals that there are only three linearly independent variables. This can be seen by taking the linearized material derivative of (2.17):

(2.34) $$\begin{eqnarray}\frac{1}{\bar{{\it\rho}}}\left(\frac{\partial {\it\rho}^{\prime }}{\partial t}+v\frac{\text{d}\bar{{\it\rho}}}{\text{d}z}\right)=\frac{1}{K}\left(\frac{\partial p}{\partial t}+v\frac{\text{d}\bar{p}}{\text{d}z}\right)-a\left(\frac{\partial n}{\partial t}+v\frac{\text{d}\bar{n}}{\text{d}z}\right)\!.\end{eqnarray}$$

From (2.28), (2.31) and (2.32), and the identity (2.21),

(2.35) $$\begin{eqnarray}\frac{1}{\bar{{\it\rho}}}\frac{\text{d}\bar{{\it\rho}}}{\text{d}z}=\frac{1}{K}\frac{\text{d}\bar{p}}{\text{d}z}-a\frac{\text{d}\bar{n}}{\text{d}z},\end{eqnarray}$$

so all the advection terms (those involving $v$ ) in (2.34) cancel. We can then integrate (2.34) in time at fixed $z$ to find

(2.36) $$\begin{eqnarray}{\it\rho}^{\prime }=\bar{{\it\rho}}(p/K-an).\end{eqnarray}$$

This identity does not hold for more general base states in which a background flow is present.

Linearizing the momentum balance equation (2.5), mass balance equation (2.3), and BGR equation (2.4) yields

(2.37) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{{\it\rho}}\frac{\partial v}{\partial t}+\frac{\partial p}{\partial z}+{\it\rho}^{\prime }g=0, & \displaystyle\end{eqnarray}$$

(2.38) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{\bar{{\it\rho}}}\left(\frac{\partial {\it\rho}^{\prime }}{\partial t}+v\frac{\text{d}\bar{{\it\rho}}}{\text{d}z}\right)+\frac{1}{A}\frac{\partial (Av)}{\partial z}=0, & \displaystyle\end{eqnarray}$$

(2.39) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial n}{\partial t}+v\frac{\text{d}\bar{n}}{\text{d}z}+\frac{n+bp}{{\it\tau}}=0, & \displaystyle\end{eqnarray}$$

respectively. The following steps then simplify the governing equations: use (2.36) to eliminate ${\it\rho}^{\prime }$ in (2.37). Use (2.34) to replace the material derivative of ${\it\rho}^{\prime }$ in (2.38) with material derivatives of $p$ and $n$ ; then use (2.28) to eliminate $\text{d}\bar{p}/\text{d}z$ , and use (2.39) to eliminate the material derivative of $n$ . Finally use (2.32) to eliminate $\text{d}\bar{n}/\text{d}z$ in (2.39). What remains are three governing equations for $v,p$ , and $n$ :

(2.40) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{{\it\rho}}\frac{\partial v}{\partial t}+\frac{\partial p}{\partial z}+\frac{\bar{{\it\rho}}g}{K}p-\bar{{\it\rho}}gan=0, & \displaystyle\end{eqnarray}$$

(2.41) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{K}\frac{\partial p}{\partial t}+\frac{1}{A}\frac{\partial (Av)}{\partial z}-\frac{\bar{{\it\rho}}g}{K}v+\frac{a}{{\it\tau}}\left(n+bp\right)=0, & \displaystyle\end{eqnarray}$$

(2.42) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial n}{\partial t}+\bar{{\it\rho}}gbv+\frac{n+bp}{{\it\tau}}=0. & \displaystyle\end{eqnarray}$$

In the limit of thermodynamic equilibrium ( ${\it\tau}\rightarrow 0$ ), often taken in quasi-steady volcanic conduit flow studies (e.g. Wilson & Head Reference Wilson and Head1981; Melnik & Sparks Reference Melnik and Sparks1999; Mastin Reference Mastin2002; Koyaguchi Reference Koyaguchi2005; Gonnermann & Manga Reference Gonnermann and Manga2006), these equations may be simplified (by utilizing the identities $n=-bp$ and (2.21)) to

(2.43) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{{\it\rho}}\frac{\partial v}{\partial t}+\frac{\partial p}{\partial z}+\frac{\bar{{\it\rho}}g}{K_{eq}}p=0, & \displaystyle\end{eqnarray}$$

(2.44) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{K_{eq}}\frac{\partial p}{\partial t}+\frac{1}{A}\frac{\partial (Av)}{\partial z}-\frac{\bar{{\it\rho}}g}{K_{eq}}v=0. & \displaystyle\end{eqnarray}$$

2.4 Relative motion between gas and liquid

Having derived the linearized equations governing perturbations about the magmastatic state, we now discuss a basic assumption of this model, that of negligible relative motion between the gas and liquid phases. A full discussion and justification of this assumption is beyond the scope of this paper, but here we point out a few relevant issues. Relative motion is likely most important for the background state. Buoyancy drives upward motion of gas bubbles through the liquid melt, and out through the open vent and permeable conduit walls. Bubble dynamics plays a primary role in bulk magma motion, and may govern eruption style in low-viscosity settings such as we are studying (Houghton & Gonnermann Reference Houghton and Gonnermann2008). Gas loss will increase density, bulk modulus, and sound speed relative to the magmastatic state used in this study, especially in the shallow part of the conduit. Rise of gas bubbles will also cause departures from the hydrostatic balance condition (2.28). Furthermore, in some cases the gas volume fraction might increase to the point where interconnected pathways between bubbles develop to form a permeable bubble network (Kozono & Koyaguchi Reference Kozono and Koyaguchi2009), even in low-viscosity magmas (Rust & Cashman Reference Rust and Cashman2011).

Relative motion might also be important when considering perturbations about a base state. When that base state features upward gas flow, there will be advection terms in the governing equations. We suspect these terms will be negligible, given that the rise distance of even large bubbles or Darcian gas flow in a permeable network over VLP periods is small compared to vertical wavelengths of interest. Additionally, we can estimate the time scale over which perturbations in gas and liquid velocities equilibrate in response to interaction forces in the momentum balance equations. Assuming Stokes drag for the bubbles (providing an interaction force proportional to the difference in gas and liquid velocities), that adjustment time scale is of the order of $r_{B}^{2}{\it\rho}_{l}/{\it\eta}_{l}$ , where $r_{B}$ is the bubble radius, ${\it\rho}_{l}$ is liquid density, and ${\it\eta}_{l}$ is liquid viscosity (e.g. Landau & Lifshitz Reference Landau and Lifshitz1959). Even for large bubbles ( $r_{B}\sim 10^{-2}$  m) in very low viscosity melt ( ${\it\eta}_{l}/{\it\rho}_{l}\sim 10^{-2}$  m $^{2}$  s $^{-1}$ ), the adjustment time scale is only ${\sim}$ 1 s, still less than the ${\sim}$ 30 s periods of the VLP events of interest. For more typical parameters, the adjustment time scale is many orders of magnitude smaller than periods of interest.

In any case, a rigorous examination of the role of relative motion, particularly in the base state, is warranted.

2.5 Boundary conditions

We now consider a conduit of length $L$ , extending between $z=0$ at depth to $z=L$ (prior to perturbations) at the surface. One boundary condition is required at each end of the conduit. Motivated by our interest in open-vent systems with lava lakes (such as Kilauea and Erebus), the top surface is assumed to be at constant atmospheric pressure:

(2.45) $$\begin{eqnarray}P(L+h(t),t)=P_{a}.\end{eqnarray}$$

While the surface position varies in time, the amplitude of its motions, $h(t)$ , is small compared to the height of the column $L$ . Thus, we expand the pressure at the surface in a Taylor series to first order:

(2.46) $$\begin{eqnarray}P(L+h(t),t)\approx P(L,t)+h(t)\left.\frac{\partial P}{\partial z}\right|_{z=L}.\end{eqnarray}$$

It suffices to evaluate the pressure gradient in the magmastatic state using (2.28), yielding the linearized boundary condition

(2.47) $$\begin{eqnarray}P(L,t)-\bar{{\it\rho}}(L)gh(t)=P_{a}.\end{eqnarray}$$

The evolution of the surface position $h(t)$ follows from the kinematic condition

(2.48) $$\begin{eqnarray}\text{d}h/\text{d}t=v(L+h(t),t),\end{eqnarray}$$

which is similarly linearized to

(2.49) $$\begin{eqnarray}\text{d}h/\text{d}t=v(L,t).\end{eqnarray}$$

The base of the conduit is connected to a magma reservoir. In this initial study we idealize such connection in terms of either constant-pressure or zero-velocity (no-flow) boundary conditions at the bottom of the conduit: $p(0,t)=0$ or $v(0,t)=0$ , respectively. In § 5.4 we introduce a simple model for the magma reservoir, and establish conditions under which the coupled model reduces to one of the simpler boundary conditions used in this work. We discuss possible more sophisticated models of the coupled system in § 6. In the examples that follow, except when stated otherwise, we use the constant-pressure condition.

To summarize, the linearized boundary conditions are

(2.50) $$\begin{eqnarray}\displaystyle & p(0,t)=0, & \displaystyle\end{eqnarray}$$

(2.51) $$\begin{eqnarray}\displaystyle & p(L,t)-\bar{{\it\rho}}(L)gh(t)=0, & \displaystyle\end{eqnarray}$$

and we must add (2.49) to governing equations (2.40)–(2.42).

2.6 Energy balance

In this section, we derive the balance of mechanical energy $E$ for a finite length column of fluid, in anticipation of the application to magmatic conduits. The energy balance is obtained by multiplying the momentum balance (2.40) by $Av$ and integrating over the conduit. After utilizing the other governing equations, properties of the base state, and boundary conditions, we obtain

(2.52) $$\begin{eqnarray}E=\int _{0}^{L}\left(\frac{\bar{{\it\rho}}v^{2}}{2}+\frac{p^{2}}{2K}+\frac{an^{2}}{2b}\right)A\,\text{d}z+\frac{A(L)\bar{{\it\rho}}(L)gh^{2}}{2},\end{eqnarray}$$

which evolves in time according to

(2.53) $$\begin{eqnarray}\frac{\text{d}E}{\text{d}t}=-\int _{0}^{L}\frac{a(n+bp)^{2}}{b{\it\tau}}A\,\text{d}z.\end{eqnarray}$$

The total energy is the sum of kinetic energy and several potential energy terms that appeared in our discussion of the energetics of a fluid parcel (2.25) and (2.26). The last term in (2.52) is the gravitational potential energy associated with motions of the free surface. As before, we refer to $p^{2}/2K$ as the non-equilibrium free energy density and $p^{2}/2K+an^{2}/2b$ as the total free energy density. The latter is always less than or equal to the equilibrium free energy density $p^{2}/2K_{eq}$ . The right-hand side of (2.53) is the energy dissipation rate from non-equilibrium BGR.

In the case of equilibrium BGR, the energy simplifies to

(2.54) $$\begin{eqnarray}E=\int _{0}^{L}\left(\frac{\bar{{\it\rho}}v^{2}}{2}+\frac{p^{2}}{2K_{eq}}\right)A\,\text{d}z+\frac{A(L)\bar{{\it\rho}}(L)gh^{2}}{2},\end{eqnarray}$$

and is conserved, $\text{d}E/\text{d}t=0$ .

Equation (2.53) establishes that perturbations to the magmastatic base state uniformly decay in time. Energy loss is due to non-equilibrium effects (the last term on the right-hand side of (2.53)). Although we assume inviscid flow here, had we included viscous effects and wall drag there would be an additional dissipation term appearing in the energy balance.

2.7 Magmatic application

We must specify an appropriate mixture equation of state and a gas solubility law to close our problem. For the magmatic systems of interest, the gas phase will follow an ideal gas law to a good approximation, and a linearized model for a compressible liquid is appropriate. Density is given by the mixture equation (2.7), with

(2.55a,b ) $$\begin{eqnarray}{\it\rho}_{g}(P)=\frac{P}{RT}\quad \text{and}\quad {\it\rho}_{l}(P)={\it\rho}_{liq}\left(1+\frac{P-P_{0}}{K_{liq}}\right)\!,\end{eqnarray}$$

with $R=462$  J kg $^{-1}$  K the ideal gas constant and $T=1273$ K a temperature representative of basaltic magma, ${\it\rho}_{liq}=3000$ kg m $^{-3}$ a reference magma liquid density, $K_{liq}=10$  GPa the liquid bulk modulus, and $P_{0}$ is a reference pressure about which the liquid density is linearized. In what follows we deal with a linearized form of this equation of state as in (2.17).

We assume at first that water is the only volatile species present, in which case equilibrium solubility $n_{eq}(P)$ approximately follows from Henry’s law,

(2.56) $$\begin{eqnarray}n_{eq}(P)=n_{t}-S_{w}\sqrt{P},\end{eqnarray}$$

with $n_{t}$ the total water mass fraction. Equation (2.56) can be modified to include nucleation and supersaturation effects (Woods Reference Woods1995). $S_{w}\approx 4\times 10^{-6}$ Pa $^{-1/2}$ is an empirical constant, which is an approximation to the equilibrium solubility of water in silicate melts. In general there are multiple volatile species present in magma, but water is nearly always the dominant volatile by at least an order of magnitude (Wallace Reference Wallace2005). Volatile exsolution is dominantly pressure-dependent, with single-phase fluid (no gas bubbles) above a critical pressure $P_{\mathit{exsolve}}$ . For Henry’s law and water the only volatile species, $P_{\mathit{exsolve}}=(n_{t}/S_{w})^{2}$ .

Figure 3 illustrates the variation of background variables with depth for the case of total water content $n_{t}=0.01$ . Spatial gradients in gas content, and therefore density and pressure, are most pronounced just above the depth of exsolution, defined as the value of $z$ for which $\bar{p}(z)=P_{\mathit{exsolve}}$ . Differences between equilibrium and non-equilibrium sound speeds are also largest at this transition. The scale height for background stratification (the length scale at and above which effects of gravity are non-negligible) is $C_{eq}^{2}/g\sim 25{-}500$  m between the depth of gas exsolution and the surface.

Figure 3.

Depth dependence of background fluid properties in the magmastatic state, assuming total gas mass fraction $n_{t}=0.01$ . (a) Pressure $P$ , (b) density ${\it\rho}$ and gas mass fraction $n_{g}$ , (c) equilibrium and non-equilibrium sound speeds $C_{eq}$ and $C$ , and (d) depth distribution of the dimensionless parameters ${\it\gamma},{\it\alpha}$ , and ${\it\sigma}$ (3.5)–(3.7) for two representative wavenumbers and ${\it\tau}=1$ s. While ${\it\sigma}$ and ${\it\alpha}$ are not defined below the exsolution depth (at ${\approx}$ 1420 m), ${\it\gamma}$ takes a small but finite value below exsolution ( ${\it\gamma}=3\times 10^{-5}$ , $3\times 10^{-4}$ for $k=0.1$ , 0.01 m $^{-1}$ , respectively).

3.1 Zero-gravity limit

If ${\it\gamma}=0$ , gravitational restoring force in the momentum balance (3.2) is eliminated, as is stratification in the magmastatic state. Stratification enters through advection terms in (3.3) and (3.4). This limit is obtained by taking $g\rightarrow 0$ for finite wavenumber $k$ (for example by considering a horizontal instead of vertical conduit) such that ${\it\sigma}$ remains finite. Non-equilibrium BGR provides the only deviation from standard acoustics, where the solutions are well known (undamped waves propagating up and down the conduit at speed $C$ , equipartition between kinetic and non-equilibrium free energy, Lighthill Reference Lighthill1978). The characteristic equation reduces to

(3.10) $$\begin{eqnarray}1+{\it\sigma}s^{\ast }+(1+{\it\alpha})s^{\ast 2}+{\it\sigma}s^{\ast 3}=0.\end{eqnarray}$$

This cubic equation has three roots, one of which is always real and non-positive. This root quantifies the rate of mass exchange between the phases. The other two solutions correspond to waves, and there are two limits of interest. In the high-frequency limit there is not time for mass exchange over a wave period, so ${\it\sigma}\rightarrow \infty$ , (3.10) reduces to $s^{\ast }(1+s^{\ast 2})=0$ and the wave modes are $s^{\ast }=\mp \text{i}$ . These waves are undamped and propagate at the non-equilibrium sound speed $C$ without dispersion. In the low-frequency limit, ${\it\sigma}\rightarrow 0$ and thermodynamic equilibrium is maintained over a wave period. Equation (3.10) reduces to $1+(1+{\it\alpha})s^{\ast 2}=0$ and propagation occurs at the equilibrium sound speed $C_{eq}$ with $s^{\ast }=\mp \text{i}/\sqrt{1+{\it\alpha}}$ .

The differences between equilibrium and non-equilibrium fluid properties, as measured by ${\it\alpha}$ , become most significant as the gas exsolution pressure is approached (figure 3). Although wavelike solutions for most of the parameter space come in complex conjugate pairs, there are parameters for which all roots of (3.10) are purely real and negative. This overdamped regime of perturbations coincides with the transition between equilibrium and non-equilibrium wave propagation, and occurs where the discriminant of (3.10) is positive, or where

(3.11) $$\begin{eqnarray}-4{\it\sigma}^{4}+{\it\sigma}^{2}({\it\alpha}(20+{\it\alpha})-8)-4(1+{\it\alpha})^{3}\geqslant 0.\end{eqnarray}$$

The quadratic equation in ${\it\sigma}^{2}$ may only be satisfied for ${\it\alpha}\geqslant 8$ , which for realistic parameters is generally restricted to a narrow region above the exsolution depth. We note that similar overdamped behaviour is found in other two-phase models that consider relative motion between phases (e.g. Bercovici & Michaut Reference Bercovici and Michaut2010).

3.2 Long-wavelength limit

We now take ${\it\gamma}\rightarrow \infty$ (by taking $k\rightarrow 0$ with finite $g$ ), which corresponds to a long-wavelength limit in which the effects of gravity and stratification dominate. The pressure gradient term is negligible in the momentum balance and the velocity divergence term is negligible in the mass balance equation. Instead, terms involving $g$ become important. These describe the gravitational restoring force associated with density perturbation in the momentum balance, and the change in density arising from advection through a stratified background state. Compressibility is negligible. The purely real root approaches $s^{\ast }=-(1+{\it\alpha})/{\it\sigma}$ , which moves off towards $-\infty$ in the ${\it\sigma}\rightarrow 0$ thermodynamic equilibrium limit (corresponding to instantaneous adjustment to equilibrium).

For the two wavelike solutions, we anticipate (and can verify a posteriori) that $s^{\ast }$ is $O({\it\gamma})$ . Removing the purely real root and dropping all but the highest-order terms in ${\it\gamma}$ , the dispersion relation reduces to

(3.12) $$\begin{eqnarray}{\it\gamma}^{2}(1+{\it\alpha})+s^{\ast 2}=0,\end{eqnarray}$$

irrespective of the value of ${\it\sigma}$ . The solutions, $s^{\ast }=\mp \text{i}{\it\gamma}\sqrt{1+{\it\alpha}}$ , describe undamped oscillations at the Brunt–Väisälä frequency $N_{b}$ defined in (2.33).

3.3 General case

We now return to the general case. Figure 4 presents the wavelike solutions to the dispersion relation (3.9). In the ${\it\gamma}=0$ limit, ${\it\sigma}\gg 1$ is the non-equilibrium BGR regime, while ${\it\sigma}\ll 1$ is the equilibrium BGR regime (figure 4 a). The value ${\it\alpha}$ determines the extent to which equilibrium and non-equilibrium regimes are distinct, and as well whether there is an overdamped regime (dashed white contour). Figure 4(b) examines another slice through this parameter space, at ${\it\alpha}=10$ for non-zero ${\it\gamma}$ , and we see that for large ${\it\gamma}$ , phase velocity follows the gravity wave limit, independent of ${\it\sigma}$ .

Figure 4.

(a) Dimensionless phase velocity $\text{Im}[s^{\ast }]$ in the zero-gravity limit ( ${\it\gamma}=0$ ) described by (3.10), as a function of ${\it\alpha}$ and ${\it\sigma}$ . Within the dashed white contour, solutions are overdamped and non-propagating. The contour is defined by (3.11). The phase velocity is determined primarily by ${\it\alpha}$ , which quantifies differences between the equilibrium and non-equilibrium sound speeds $C_{eq}$ and $C$ . (b) Regime diagram for fixed ${\it\alpha}=10$ as a function of ${\it\gamma}$ and ${\it\sigma}$ , showing the parameter dependence of the three propagation regimes. The red curve is the trajectory through parameter space as wavenumber in decreased with all other properties fixed; solutions along this curve are plotted in figure 5.

In figure 5 we explore these three propagation regimes as wavenumber is varied. Changing $k$ alters both ${\it\sigma}$ and ${\it\gamma}$ , while leaving ${\it\alpha}$ and ${\it\sigma}/{\it\gamma}$ fixed. These three regimes are plotted explicitly for several choices of ${\it\alpha}$ in figure 5(a), with thin lines corresponding to the phase velocity limits discussed above.

Figure 5.

Effect of varying wavenumber while holding all other properties fixed. This changes both ${\it\sigma}$ and ${\it\gamma}$ with their ratio held fixed, as shown in figure 4 (and see text). (a) Phase velocity $s^{\ast }$ of eigenmodes as a function of ${\it\sigma}/(2{\rm\pi})$ , for several choices of ${\it\alpha}$ . Thin lines outline the three wave propagation regimes: non-equilibrium BGR, equilibrium BGR, and gravity wave limits. For ${\it\alpha}=100$ we see a region of overdamped waves (zero phase velocity). For all others there is a symmetric conjugate wave (not plotted). (b) Quality factor $Q$ (3.13) for the cases in (a). (c–e) Non-equilibrium (dashed) and total (solid) free energy density, relative to equilibrium free energy (§ 2.2) for the modes in (a). Ratio of non-equilibrium to equilibrium free energy density equals $K_{eq}/K=(1+{\it\alpha})^{-1}$ . Differences between total and equilibrium free energy indicate departure from thermodynamic equilibrium.

Even when there are no overdamped modes, the transitional regime around ${\it\sigma}/(2{\rm\pi})\sim 1$ , where the wave period is similar to ${\it\tau}$ , represents the most damped part of parameter space. As figure 5(a,b) shows with a family of curves, if we define a quality factor (damping over one oscillation cycle) for propagating modes as

(3.13) $$\begin{eqnarray}Q=\frac{\text{Im}(s)}{2\text{Re}(s)},\end{eqnarray}$$

the least damped (highest Q) modes correspond to the short- and long-wavelength limits. The region around ${\it\sigma}/(2{\rm\pi})\sim 1$ is the most highly damped, optimizing the dissipation from non-equilibrium phase exchange (2.53). Figure 5(a) plots only one of the complex conjugate roots, except for the region of overdamped propagation for large ${\it\alpha}$ (in which there are not conjugate roots). Curves in figure 5(a,b) can also be interpreted as the dispersion relation for coefficients frozen at increasing depths in the conduit (pressures) as ${\it\alpha}$ increases (figure 3 d). ${\it\gamma}$ increases by a factor of ${\sim}50$ from ${\it\alpha}=0.1$ to ${\it\alpha}=100$ at fixed wavenumber.

3.4 Energetics

For each of the solutions of the characteristic equation (3.9) at a given wavelength, we can evaluate the partition of total energy into kinetic and potential energy. The gravitational energy term does not appear, as the domain is periodic. Figure 5(c–e) plots contributions to the potential energy of the propagating modes from panel (a) ( ${\it\alpha}=0.1,1,10$ ) according to the partitioning defined by (2.52).

For wave periods much smaller than ${\it\tau}$ ( ${\it\sigma}\gg 1$ ), there is negligible BGR over a wave period and energy is equally partitioned between kinetic and non-equilibrium free energy, as is the case for simple acoustics. For ${\it\sigma}\ll 1$ , BGR keeps the system in thermodynamic equilibrium and the mixture is more compressible than in the non-equilibrium case; total energy is equally partitioned between kinetic and equilibrium free energy (2.54). Figure 5(c–e) illustrates this transition in fluid compressibility, and hence free energy, and highlights its dependence on the time scale of the wave motion relative to ${\it\tau}$ as quantified through ${\it\sigma}$ . It is also clear that ${\it\alpha}$ quantifies the difference in fluid compressibility; as ${\it\alpha}\rightarrow 0$ , the difference between non-equilibrium and equilibrium free energy becomes negligible.

4.1 Numerical method

We therefore develop a numerical method that can be used to both solve the eigenvalue problem and study excitation of these eigenmodes in the time domain. We return to dimensional governing equations, but now discretize in space using finite differences. Specifically, a field such as $v(z,t)$ is discretized on a finite domain $0\leqslant z\leqslant L$ at a set of $M+2$ possibly non-uniformly spaced grid points $z_{i}$ . The discontinuity in properties at the gas exsolution depth is handled using split nodes, that is, a pair of collocated points, one on either side of the material interface, with appropriate conditions enforcing mass and momentum balance across the interface. A field such as $v(z,t)$ is approximated at the grid points, with the grid values stored in a vector $\boldsymbol{v}(t)$ having components $v_{i}(t)\approx v(z_{i},t)$ .

We use summation-by-parts (SBP) finite differences (Kreiss & Scherer Reference Kreiss, Scherer and de Boor1974, Reference Kreiss and Scherer1977; Strand Reference Strand1994; Gustafsson, Kreiss & Oliger Reference Gustafsson, Kreiss and Oliger1995). Boundary conditions (and interface conditions across the material interface at the exsolution depth) are enforced weakly with the simultaneous approximation term (SAT) method (Carpenter, Gottlieb & Abarbanel Reference Carpenter, Gottlieb and Abarbanel1994). With the SBP–SAT framework (e.g. Svard & Nordstrom Reference Svard and Nordstrom2014) we can construct a discrete energy balance that closely mimics the continuous energy balance (2.53) and prove that numerical solutions are both stable and high-order accurate. The examples shown below use a globally fifth-order-accurate difference operator. Derivation of the discrete version of governing equations (2.40)–(2.42) and (2.49), boundary and interface conditions, and the energy balance (2.53) are presented in appendix A.

This method is first used to solve the eigenproblem that arises after Laplace transforming, with Laplace parameter $s$ , the semi-discrete governing equations:

(4.1) $$\begin{eqnarray}s\tilde{\boldsymbol{u}}=\unicode[STIX]{x1D63D}\tilde{\boldsymbol{u}},\end{eqnarray}$$

where $\boldsymbol{u}(t)=[v_{0}(t),\ldots ,v_{M+2}(t),p_{0}(t),\ldots ,p_{M+2}(t),n_{0}(t),\ldots ,n_{M+2}(t),h(t)]^{\text{T}}$ is a column vector of velocity, pressure perturbation, and gas fraction perturbation at the grid points along with the free surface displacement $h$ , and

(4.2) $$\begin{eqnarray}\tilde{\boldsymbol{u}}(s)=\int _{0}^{\infty }\boldsymbol{u}(t)\text{e}^{-st}\,\text{d}t\end{eqnarray}$$

defines the Laplace transform. $\unicode[STIX]{x1D63D}$ is a $3(M+2)+1\times 3(M+2)+1$ matrix containing spatially variable fluid properties evaluated at the grid points, the finite difference operator, and penalty terms that enforce interface and boundary conditions (see appendix A for details). There are $3(M+2)+1$ eigenvalues, which we store for later use in the diagonal matrix  ${\it\bf\Lambda}$ .

Diagonalization of $\unicode[STIX]{x1D63D}$ defines the eigenmodes of the system, with the corresponding eigenvalues characterizing the oscillation frequencies and decay rates of modes. In contrast to the frozen coefficient case, where eigenmodes took the form of Fourier modes, eigenmodes of the finite conduit exhibit a complex spatial structure that reflects the strong gradients in background fluid properties. In particular, a large impedance contrast at the exsolution interface serves to segment the conduit into two parts, with significantly different properties due to the presence of bubbles and increased compressibility in the upper portion of the conduit.

4.2 Eigenmode families

Results for a 3 km length conduit with $n_{t}=0.01$ , $p=0$ bottom boundary condition, and three choices of ${\it\tau}$ are shown in figures 6 and 7. The eigenfunction amplitudes are only defined up to an arbitrary constant. We choose this constant in figure 6 so that total free energy density has unit maximum amplitude, which facilitates comparison of mode families. For figure 7 the constant is chosen to make the total energy of each mode equal to 10 kJ. Note that $n$ is not defined, and hence not plotted, below the exsolution depth in figure 7.

Figure 6.

Families of eigenmodes for (a,e,h) ${\it\tau}=0$ s (equilibrium BGR), (b,f,i) ${\it\tau}=1$ s, (c,g,j) ${\it\tau}=1000$ s (non-equilibrium BGR). Ten longest-period eigenmodes in each family are shown, with the normalization and family naming scheme described in the text. (d) Background gas mass fraction, showing segmentation of conduit into bubbly magma above exsolution and liquid melt below. (a–c) Spatial distribution of total free energy density. Bright colours indicate regions where magma is compressed or expanded. (e–g) Mode period, dots coloured by decay rate. (h–j) Energetics. (a,e,h) shows contributions to potential energy from equilibrium free energy (which equals total free energy) and gravitational potential energy. In panels (b,f,i,c,g,j), potential energy is divided into total free energy and gravitational potential energy; also shown is non-equilibrium free energy (which differs from total free energy when BGR is appreciable over a mode period). White arrows in (b) refer to modes plotted in figure 7.

Figure 7.

Three eigenmodes, labelled with white arrows in figure 6(b), for $p=0$ bottom boundary condition and ${\it\tau}=1$  s, with normalization described in text. Modes are uniformly scaled to reflect typical excitation amplitudes presented here. (a,d,g) Kinetic and free energy density. (b,e,h) and (c,f,i) Eigenfunctions for pressure perturbation, velocity, and exsolved gas mass fraction perturbation. Particular labelled modes have been uniformly scaled for visualization purposes.

We find three families of modes, labelled in figure 6(a–c). Family I contains the longest-period modes of the conduit. Most of these modes, with the exception of the fundamental mode discussed below, are localized primarily above the exsolution depth. Family II contains shorter-period eigenmodes. When ${\it\tau}$ is comparable to the mode periods ( ${\sim}$ 1 s in this example), these modes exhibit dominant or even exclusive spatial localization below the exsolution depth. Family III exists only for non-zero ${\it\tau}$ , and contains modes that are non-oscillatory. These monotonically decaying modes are associated with mass exchange between the phases, similar to the frozen coefficient case. The modes shown in figure 6 are arranged in this order, from left to right, with the ten longest-period modes from each family shown. Note that the non-oscillatory modes in family III are absent for the ${\it\tau}=0$ (equilibrium BGR) case in figure 6(a).

The mode Ex 1 in figures 6(b,f,i) and 7(a) is the longest-period mode of the conduit for the $p=0$ bottom boundary condition. We refer to this mode as the fundamental mode. The $p$ , $v$ , and $n$ eigenfunctions have no nodes (zero crossings). We classify this mode as part of family I, though it is somewhat unique and spans the entire conduit. Furthermore, as discussed later, this mode does not exist for the $v=0$ bottom boundary condition. The $p=0$ bottom boundary condition permits flow in and out of the bottom of the conduit. Kinetic energy is largest below the exsolution depth, while free energy is largest immediately above this depth. Thus, liquid melt at depth oscillates up and down, compressing and expanding the bubbly magma just above the exsolution depth. Gravity provides an additional restoring force, although small compared to magma compressibility, for this mode.

The second-longest-period mode (Ex 2 in figures 6 b,f,i and 7 a) has somewhat different characteristics, and is more representative of other members of family I. The $p$ , $v$ , and $n$ eigenfunctions each have one node. Both kinetic and free energy are largest in the bubbly magma. The bubbly magma in the upper 1 km of the conduit serves as the mass, with restoring force coming from compressing and expanding the bubbly magma immediately above the exsolution depth and also just below the free surface. Motions are largely decoupled from the liquid melt region at depth. Gravity provides an additional restoring force for this mode of approximately equal importance relative to magma compressibility. Gravity becomes increasingly less important for higher-frequency modes in this family (figure 7 h) as the spatial wavelength of oscillations becomes small compared to the density scale height.

The mode Ex 3 (figures 6 b,f,i and 7 a) is localized below the exsolution depth. This mode, and others in family II, are somewhat akin to standing acoustic waves in the liquid melt region, relatively uninfluenced by gravity, with a $p=0$ boundary condition at the bottom and high degree of reflection off of the exsolution surface due to the large impedance contrast there. The mode is expressed as very short wavelength oscillations in the bubbly magma above the exsolution depth (due to the much lower sound speed there).

We next discuss how ${\it\tau}$ influences the oscillatory modes. For the fundamental mode (Ex 1) and other modes in family I (such as Ex 2), the differences between ${\it\tau}=0$ (equilibrium BGR, figure 6 a,e,h), ${\it\tau}=1$ s (figure 6 b,f,i), and ${\it\tau}=1000$ s (non-equilibrium BGR, figure 6 c,g,j) are relatively minor. These differences arise from the increased compressibility of the bubbly magma in the equilibrium limit. The differences between equilibrium and non-equilibrium properties are largest immediately above the exsolution depth (figure 3). When ${\it\tau}$ is comparable to the mode period (not shown), energy dissipation from BGR processes increases.

In contrast, ${\it\tau}$ exerts a much stronger influence on modes in family II (such as Ex 3). For values of ${\it\tau}$ close to the mode period (e.g. ${\it\tau}=1$ s, figure 6 b,f,i), there is little contribution to the overall kinetic or potential energy from the bubbly magma region, because this is where energy dissipation from non-equilibrium BGR occurs. Because dissipation vanishes as ${\it\tau}\rightarrow 0$ or $\infty$ , the amplitude of the eigenfunctions above the exsolution depth is larger in these limits (figure 6 a,e,h,c,g,j). The impedance contrast across the exsolution depth also changes with ${\it\tau}$ . This contrast is evidently more important for these shorter-period (and hence shorter-wavelength) modes than for the longer-period modes in family I, as illustrated by more substantial differences in the modes for the ${\it\tau}=0$ and 1000 s cases.

5.1 Connection to conduit eigenmodes and to observables

Before we present numerical solutions, we point out that time-dependent solutions to the inhomogeneous equations do not directly provide information about which of the conduit eigenmodes from (4.1) are excited. Neither do they reveal how excitation of waves within the conduit relates to potential observations. There is, however, a formal procedure for connecting these two calculations. We first Laplace transform the semi-discrete, inhomogeneous governing equations (A 2)–(A 5), yielding a linear system of the form

(5.3) $$\begin{eqnarray}s\tilde{\boldsymbol{u}}=\unicode[STIX]{x1D63D}\tilde{\boldsymbol{u}}+\boldsymbol{b}_{FI}\tilde{f}_{FI}(s)+\boldsymbol{b}_{FM}\tilde{f}_{FM}(s),\end{eqnarray}$$

where $\tilde{\boldsymbol{u}}$ and $\unicode[STIX]{x1D63D}$ are defined as before, and the two terms of the form $\boldsymbol{b}\tilde{f}(s)$ are the Laplace-transformed source vectors associated with surface impulses (subscript $FI$ ) and internal volumetric or mass sources (subscript $FM$ ). We have assumed a source process that permits us to factor the source into a vector $\boldsymbol{b}$ capturing the spatial dependence of the source and another function, $f(t)$ , carrying the source time dependence and magnitude. For this special case $\boldsymbol{b}$ is independent of Laplace parameter $s$ . Discussion of $\boldsymbol{b}$ and $f$ is presented in appendix A.

$\unicode[STIX]{x1D63D}$ can be decomposed as $\unicode[STIX]{x1D63D}=\unicode[STIX]{x1D652}{\it\bf\Lambda}\unicode[STIX]{x1D652}^{-1}$ , with $\unicode[STIX]{x1D652}$ a matrix containing the eigenvectors of $\unicode[STIX]{x1D63D}$ as columns and ${\it\bf\Lambda}=\text{diag}({\it\lambda}_{1},\ldots ,{\it\lambda}_{3(M+2)+1})$ a diagonal matrix of eigenvalues ${\it\lambda}_{i}$ (denoted as $s$ in the previous section on the eigenvalue problem). As mentioned in § 3, $\unicode[STIX]{x1D63D}$ is non-normal, so although $\unicode[STIX]{x1D63D}$ is diagonalizable, eigenfunctions are not necessarily orthogonal ( $\unicode[STIX]{x1D652}$ is not unitary). We have experimented with partially symmetrizing the eigenproblem by rescaling $\boldsymbol{u}(t)$ such that its $L_{2}$ norm is the total mechanical energy density, which greatly reduces the condition number of the resulting coefficient matrix and improves numerical accuracy. However, all calculations have shown that the longest-period modes are nearly orthogonal and we will show that we can accurately recover the time-domain solution with the method described below. This suggests that our problem is not far from normal, at least for the range of parameters we have studied thus far.

Suppose we are interested in a particular set of observations $\boldsymbol{y}(t)$ , which are a linear functional of the fields $\boldsymbol{u}(t)$ . Examples of observations include the ground displacements in an elastic half-space caused by pressure changes applied to the surface of the conduit, motions of the lava lake surface, and infrasound signals generated by oscillations of the lava lake surface. From linearity it follows that there exists an observation matrix linking the observables to the fields within the conduit, expressed in the Laplace domain as

(5.4) $$\begin{eqnarray}\tilde{\boldsymbol{y}}=\mathscr{C}\tilde{\boldsymbol{u}}.\end{eqnarray}$$

In general $\mathscr{C}=\mathscr{C}(s)$ (such that the observation at time $t$ depends on the fields at all times leading up to $t$ ), but in this study we restrict attention to observables for which $\mathscr{C}$ is independent of $s$ .

Using (5.3) the solution of (5.4) is then

(5.5) $$\begin{eqnarray}\tilde{\boldsymbol{y}}=\mathscr{H}_{FI}\tilde{f}_{FI}+\mathscr{H}_{FM}\tilde{f}_{FM},\end{eqnarray}$$

where $\mathscr{H}_{FI}(s)$ and $\mathscr{H}_{FM}(s)$ are transfer functions between the observables and surface impulse and internal volumetric/mass sources, respectively, which take the form

(5.6) $$\begin{eqnarray}\mathscr{H}(s)=\mathscr{C}\unicode[STIX]{x1D652}\text{diag}\left(\frac{1}{s-{\it\lambda}_{1}},\ldots ,\frac{1}{s-{\it\lambda}_{3(M+2)+1}}\right)\unicode[STIX]{x1D652}^{-1}\boldsymbol{b}.\end{eqnarray}$$

Figure 8.

Excitation of waves by a Gaussian pressure impulse on the otherwise traction-free upper boundary. Source is centred at 10 s, with amplitude 10 kPa and duration half-width 1 s. ${\it\tau}=1$ s and $n_{t}=0.01$ . (a) Pressure at free surface and (b) height of free surface as a function of time. (c) Excitation amplitude of free surface eigenmodes from the transfer function of $h$ (5.6) for the free surface height (blue line) and pressure at exsolution (red line). (d) Pressure perturbation and (e) velocity perturbation within the conduit. Dashed line is exsolution depth. (f) Quality factor $Q$ of conduit eigenmodes (figure 6) as a function of mode period. Panels (c) and (f) together express which modes are excited and which will contribute to some particular observable field. A mode can be excited but not expressed in some field.

The time-domain solution may be recovered by inverting the Laplace transform through use of the residue theorem, or alternatively by approximating the integral as a discrete inverse Fourier transform (e.g. Davies & Martin Reference Davies and Martin1979). We find that numerically the latter is the more accurate method. This method requires that the time series is chosen to be long enough that the least damped mode decays sufficiently to prevent wrap-around effects in the source excitation, and thus is not an efficient method in general. However, we can accurately reproduce the time evolution calculated by the time-domain method for a single or multiple observables.

5.2 Source excitation by top boundary pressure forcing

We first study the problem of a pressure impulse on the upper boundary of an initially static magma column. This could arise from perturbations to the surface of a lava lake by, for example, rockfall or shallow near-surface degassing bursts. Figure 8 illustrates results from a test case using $L=3$  km, ${\it\tau}=1$  s, and $n_{t}=0.01$ , with source pressure amplitude $F_{I}=10$  kPa centred in time at $t_{FI}=10$  s after the start of the simulation, with duration ${\it\sigma}_{FI}=1$  s. This pressure pulse initially pushes the free surface down (figure 8 b), generating a downward-propagating wave carrying a pressure increase and downward velocity perturbation. At the exsolution depth ( $z\approx 1420$ m) this wave is partially reflected. The transmitted wave continues into the single-phase region at depth, where it excites high-frequency oscillations. The reflected wave causes low-frequency oscillations in the upper, bubbly region. Damping in these simulations is only through energy lost by non-equilibrium BGR, as we have assumed inviscid magma. Inclusion of viscosity would result in smaller free surface height variations and more rapid decay of oscillations.

We then use the transfer function in (5.6) to obtain excitation amplitudes of individual eigenmodes. Figure 8(c) shows the magnitude of $\mathscr{H}_{FI}(s)$ as a function of excitation period ( $2{\rm\pi}/\text{Im}(s)$ with purely imaginary $s$ ), for two choices of $y$ : the free surface height $h$ and the pressure perturbation at the exsolution depth $p_{\mathit{exsolve}}$ . The latter is a crude proxy for excitation of seismic waves and ground displacements by pressure perturbations on the conduit walls. The former is inverted to the time domain to demonstrate that the Fourier transform inversion method reproduces finite difference calculations accurately and is plotted on top of the time-dependent solution in figure 8(b). Figure 8(f) then shows the period and quality factor of eigenmodes obtained during study of the homogeneous problem. Comparison to figure 8(c) shows that the top forcing excites long-period modes at much larger amplitude than short-period modes, as observed in the time-dependent simulations.

5.3 Source excitation by an impulsive internal mass source

The second type of source we study is an impulsive, localized volumetric source within the conduit. This could arise from influx or outflux of mass from a branching sill or dike, or any other discrete dilational event. We explore a case where dilation rate is $F_{M}=0.05$  1 s $^{-1}$ , with $t_{FM}=10$  s, ${\it\sigma}_{FM}=1$  s, and spatial half-width ${\it\sigma}_{Z}=10$  m, all other parameters as in the top forcing example of § 5.2. Figure 9 shows two examples, one where the volumetric source is below exsolution ( $z_{FM}=2$  km, figure 9 a) and one where the source is above exsolution ( $z_{FM}=1$  km, figure 9 b).

Figure 9.

(a,c,e) Excitation of waves by an internal mass source, a Gaussian source with spatial extent ${\it\sigma}_{Z}=10$ m, temporal duration ${\it\sigma}_{FM}=1$ s, injected at $z_{FM}=2$  km depth (below exsolution) at time $t_{FM}=10$ s. (a) is free surface height, (c) is pressure perturbation amplitude in kPa versus time. (e) is transfer function amplitude of pressure at the exsolution depth and surface height versus mode period in seconds. (b,d,f) Same as part (a,c,e), except the mass source is injected at $z_{FM}=1$  km depth (above exsolution).

Looking first at the amplitude of the pressure perturbations, we see that sources below exsolution generate much higher pressures than sources above exsolution. This can be understood in terms of the difference in fluid properties between bubbly and non-bubbly magma: the fluid bulk modulus is ${\sim}$ 4 orders of magnitude larger at 2 km depth than at 1 km depth. From the amplitude of pressure eigenmode excitation in figure 9(e,f) (amplitude of transfer function $\mathscr{H}_{FM}(s)$ for pressure at the exsolution depth) we also see that high-frequency modes are excited in the single-phase region when the source is below the exsolution depth, but not when the source is above it. By the time perturbations reach the free surface, however, in both cases the high-frequency modes have been damped out, and the free surface height still primarily records long-period modes from the bubbly upper conduit.

5.4 Influence of bottom boundary condition

All results so far assume a constant-pressure bottom boundary condition, so fluid can flow freely in and out of the conduit from some deeper magma reservoir or chamber. This is not the only possible boundary condition for real volcanic systems. Another possibility is a no-flow (zero-velocity) boundary condition, which would lead to variations in pressure during oscillatory motion of the magma column. We can assess the conditions under which either of these boundary conditions is realistic in the context of a simple model of the reservoir that is coupled to the bottom boundary of the conduit. We assume uniform pressure throughout the reservoir, which might be a chamber or sill within an elastic solid.

We combine a statement of mass balance for the reservoir, neglecting influx from below, with linearized expressions relating changes in fluid density and reservoir volume to pressure changes (Segall Reference Segall2010). This reservoir model is coupled to the conduit by setting reservoir pressure equal to the pressure at the bottom of the conduit, $P_{b}$ , and setting volume flux out of the reservoir into the conduit equal to $A(0)V_{b}$ , where $A(0)$ is the cross-sectional area at the bottom of the conduit and $V_{b}$ is the velocity there. The result is an ordinary differential equation linking $P_{b}$ and $V_{b}$ :

(5.7) $$\begin{eqnarray}\frac{\text{d}P_{b}}{\text{d}t}=-{\it\epsilon}_{r}V_{b},\end{eqnarray}$$

where ${\it\epsilon}_{r}=A(0)/{\it\beta}_{r}V_{r}$ contains the overall compressibility ${\it\beta}_{r}$ (the sum of chamber or sill compressibility and fluid compressibility) and the reservoir volume $V_{r}$ . We have assumed for simplicity that the conduit and chamber have the same bubble content.

Fourier transforming in time, with angular frequency ${\it\omega}$ , permits us to define the impedance of the reservoir system as the ratio of the Fourier-transformed pressure and velocity:

(5.8) $$\begin{eqnarray}\frac{\tilde{P}_{b}}{\tilde{V}_{b}}=-\frac{{\it\epsilon}_{r}}{\text{i}{\it\omega}}.\end{eqnarray}$$

We then compare this impedance to the impedance of a wave in the conduit, ${\it\rho}_{b}C_{b}$ , where ${\it\rho}_{b}$ and $C_{b}$ are the density and phase velocity of the wave at the bottom of the conduit. The impedance ratio,

(5.9) $$\begin{eqnarray}{\it\Gamma}=-\frac{{\it\epsilon}_{r}}{\text{i}{\it\omega}{\it\rho}_{b}C_{b}},\end{eqnarray}$$

determines how waves reflect from the reservoir (e.g. Lighthill Reference Lighthill1978). When ${\it\Gamma}\ll 1$ , then $P_{b}\approx \text{constant}$ ; when ${\it\Gamma}\gg 1$ , then $V_{b}\approx 0$ .

We consider two examples of reservoir systems to estimate the impedance ratio. First, assume the conduit, of radius ${\sim}$ 10 m, is connected to a magma-filled sill of radius ${\sim}$ 1 km and thickness ${\sim}$ 1 m. Compressibility of the system is principally from opening/closing the sill, such that ${\it\beta}_{r}$ is approximately equal to the overpressure ${\sim}$ 10 MPa in the sill (i.e. the difference between magma pressure and compressive normal stress acting to close the sill). We also take ${\it\rho}_{b}C_{b}\sim 10^{6}$  Pa s m $^{-1}$ . Thus, ${\it\epsilon}_{r}\sim 10^{3}$ Pa m $^{-1}$ . It follows that at periods less than ${\sim}$ 1000 s the constant-pressure condition is appropriate.

As a second example, we assume the reservoir is a spherical magma chamber of radius ${\sim}$ 100 m; we again take the conduit radius to be ${\sim}$ 10 m. Furthermore, ${\it\beta}_{r}$ is approximately equal to the reciprocal of the magma bulk modulus ( ${\sim}$ 10 GPa), with a slight correction due to the elastic compliance of the chamber. In this case, ${\it\epsilon}_{r}\sim 10^{6}$  Pa m $^{-1}$ , and here for modes with periods longer than ${\sim}$ 10 s, the no-flow boundary condition is appropriate.

The major difference between the constant-pressure and no-flow bottom boundary conditions lies in the longest-period eigenmode. For the constant-pressure boundary condition, this longest period has the majority of kinetic energy residing in the bubble-free region (figure 7). Thus, it represents oscillatory mass transport in and out of the conduit. For the no-flow boundary condition, this mode is absent, and the longest-period mode is one that resides in the bubbly region (figure 10 a, dashed curve).

Figure 10.

(a) Magnitude of transfer function between top pressure forcing and pressure at the exsolution depth, varying total conduit length $L$ , comparing constant-pressure (solid) and zero-velocity (dashed) bottom boundary conditions. $n_{t}=0.01$ and ${\it\tau}=1$ s. Note that modes localized above exsolution are not sensitive to $L$ or bottom boundary condition. (b) Same as panel (a) (and with $p=0$ bottom boundary condition), but varying total volatile content with $L=3$  km. We also include a case of joint H $_{2}$ O and CO $_{2}$ solubility (Papale, Moretti & Barbato Reference Papale, Moretti and Barbato2006) that represents typical volatile contents at Kilauea, HI. Inset shows depth variation of gas mass fraction in the base state.

5.5 Influence of conduit length, total volatile content, and solubility law

We next briefly explore how mode properties and excitation (the latter expressed via transfer functions) are influenced by the conduit length, total volatile content, and solubility law.

One might reasonably speculate that the longest eigenmode periods are simply set by the conduit length $L$ , as they are in organ pipes containing a homogeneous fluid. Figure 10 explores this in more detail for the $p=0$ bottom boundary condition, using the top boundary pressure excitation and transfer function for the pressure at the exsolution depth to probe the sensitivity of mode properties to  $L$ . As we increase $L$ , some of the modes vary systematically. The fundamental mode (Ex 1) shifts to longer periods as $L$ increases (figure 10 a), consistent with an increase in effective mass (due to the increase of the length of the liquid melt region), with the restoring force from magma compression immediately above the exsolution depth relatively unaltered. Similarly, the shorter-period modes in family II (such as Ex 3, with ${\sim}$ 1 s period) that are localized within the liquid melt region at depth shift to larger periods as $L$ increases. In contrast, modes with periods between ${\sim}$ 10 and 50 s that are oscillations within the bubbly region (modes in family I such as Ex 2) are largely insensitive to $L$ , since it is the length of the bubbly magma region (set by total volatile content $n_{t}$ ) that influences their period.

Finally, switching from a $p=0$ to $v=0$ bottom boundary condition eliminates flow in and out of the conduit, and the fundamental mode no longer exists. Modes such as Ex 2 are effectively unaltered (and they are insensitive to $L$ ), while those such as Ex 3 have longer periods (close to double those for $p=0$ ) that reflect the changed boundary conditions (figure 10 a).

We next explore, for fixed $L$ , the influence of the total volatile content $n_{t}$ , and hence the length of the bubbly magma region, on mode properties (figure 10 b). Modes of all periods are sensitive to the total volatile content $n_{t}$ , but in differing ways, as illustrated by comparison of the red and blue curves in figure 10(b), corresponding to $n_{t}=0.005$ and 0.01, respectively. The period of the fundamental open-ended pipe mode (Ex 1) decreases as $n_{t}$ increases, because the mass (of the liquid melt region at depth) decreases as the exsolution surface moves to greater depths. In contrast, the periods of the long-period modes such as Ex 2 increase as $n_{t}$ increases. These modes reside within the bubbly magma above the exsolution depth, and the extent of this region increases with  $n_{t}$ .

Of course, real magmas have a more complicated volatile content, with H $_{2}$ O as well as lesser amounts of CO $_{2}$ and SO $_{2}$ (as well as other volatiles) that have much lower solubilities. Although a detailed exploration of the consequences of more complex solubility laws is outside the scope of this paper, we include in figure 10(b) a calculation using the joint H $_{2}$ O and CO $_{2}$ solubility model of Papale et al. (Reference Papale, Moretti and Barbato2006) as the gas base state. Using volatile contents typical of shallowly stored Kilauea magmas (Gerlach & Graeber Reference Gerlach and Graeber1985) as a guide, we use a H $_{2}$ O mass fraction of $n_{t,w}=0.006$ and a CO $_{2}$ mass fraction of $n_{t,c}=4\times 10^{-4}$ . The inset demonstrates that the background gas profile exhibits a more gradual decline with depth than the case of pure water, due to the persistence of CO $_{2}$ bubbles at depths greater than that at which all H $_{2}$ O bubbles have dissolved. Nevertheless, this case still exhibits a large contrast in properties at the water exsolution depth. Such a result is not unexpected, since water is more than an order of magnitude more abundant that CO $_{2}$ in this example (a result that holds generally for magmas). Therefore the segmentation of conduit eigenmodes by water exsolution should be even more pronounced for magmas that contain greater total volatile contents, such as at subduction zones, where H $_{2}$ O may reach 6–8 % by mass and CO $_{2}$ up to 0.3 % (Wallace Reference Wallace2005).